1 |
mschen |
1.1 |
#include "utils.h"
|
2 |
|
|
#include <math.h>
|
3 |
|
|
#include <stdio.h>
|
4 |
|
|
#include <iostream>
|
5 |
|
|
|
6 |
|
|
double sortMaxToMin(const double* a, const int size, const int j)
|
7 |
|
|
{
|
8 |
|
|
vector<double> b;
|
9 |
|
|
for (int i=0; i<size; i++)
|
10 |
|
|
{
|
11 |
|
|
b.push_back(a[i]);
|
12 |
|
|
}
|
13 |
|
|
sort(b.begin(), b.end());
|
14 |
|
|
reverse(b.begin(), b.end());
|
15 |
|
|
return b[j-1];
|
16 |
|
|
}
|
17 |
|
|
|
18 |
|
|
|
19 |
|
|
|
20 |
|
|
double correctZZweight(double weight,const double m4l)
|
21 |
|
|
{
|
22 |
|
|
const double NLObox = 0.2;
|
23 |
|
|
|
24 |
|
|
const double x = m4l;
|
25 |
|
|
double corr=0;
|
26 |
|
|
const double par[9] =
|
27 |
|
|
{
|
28 |
|
|
-4.9468200E-01,
|
29 |
|
|
3.0085800E-01,
|
30 |
|
|
-1.7096700E-02,
|
31 |
|
|
4.5142500E-04,
|
32 |
|
|
-6.4673600E-06,
|
33 |
|
|
5.3404600E-08,
|
34 |
|
|
-2.5399165E-10,
|
35 |
|
|
6.4625000E-13,
|
36 |
|
|
-6.8131800E-16
|
37 |
|
|
};
|
38 |
|
|
if (x<=200 && x>30)
|
39 |
|
|
{
|
40 |
|
|
corr = 0;
|
41 |
|
|
for (int p=0; p<9; p++)
|
42 |
|
|
corr = corr + par[p]*pow(x,double(p));
|
43 |
|
|
}
|
44 |
|
|
if (x>200)
|
45 |
|
|
corr = 1.595*(1-exp((-0.007888)*x));
|
46 |
|
|
|
47 |
|
|
weight = weight*(corr/1.35 + NLObox);
|
48 |
|
|
|
49 |
|
|
return weight;
|
50 |
|
|
}
|
51 |
|
|
|
52 |
|
|
|
53 |
|
|
|
54 |
|
|
|
55 |
|
|
double dR(double eta1, double phi1, double eta2, double phi2)
|
56 |
|
|
{
|
57 |
|
|
return sqrt((eta1-eta2)*(eta1-eta2) + (phi1-phi2)*(phi1-phi2));
|
58 |
|
|
}
|
59 |
|
|
|
60 |
|
|
double eta(double pz,double p)
|
61 |
|
|
{
|
62 |
|
|
if (p == 0)
|
63 |
|
|
{
|
64 |
|
|
cout << "ERROR: eta, p==0\n";
|
65 |
|
|
return 0;
|
66 |
|
|
}
|
67 |
|
|
double theta = acos(pz/p);
|
68 |
|
|
double eta = -log(tan(theta/2.0));
|
69 |
|
|
|
70 |
|
|
//if (pz<0) eta = -eta;
|
71 |
|
|
|
72 |
|
|
return eta;
|
73 |
|
|
}
|
74 |
|
|
|
75 |
|
|
double pt(double px,double py)
|
76 |
|
|
{
|
77 |
|
|
return sqrt(px*px + py*py);
|
78 |
|
|
}
|
79 |
|
|
|
80 |
|
|
double p(double px,double py,double pz)
|
81 |
|
|
{
|
82 |
|
|
return sqrt(px*px + py*py + pz*pz);
|
83 |
|
|
}
|
84 |
|
|
|
85 |
|
|
double phi(double px,double py)
|
86 |
|
|
{
|
87 |
|
|
double phi = 0;
|
88 |
|
|
|
89 |
|
|
if (px==0 && py>0)
|
90 |
|
|
phi = pi/2;
|
91 |
|
|
if (px==0 && py<0)
|
92 |
|
|
phi = -pi/2;
|
93 |
|
|
if (px>0 && py==0)
|
94 |
|
|
phi = 0;
|
95 |
|
|
if (px<0 && py==0)
|
96 |
|
|
phi = pi;
|
97 |
|
|
|
98 |
|
|
double px_ = fabs(px), py_ = fabs(py);
|
99 |
|
|
if (px>0 && py>0)
|
100 |
|
|
phi = atan(py_/px_);
|
101 |
|
|
if (px>0 && py<0)
|
102 |
|
|
phi = -atan(py_/px_);
|
103 |
|
|
if (px<0 && py>0)
|
104 |
|
|
phi = pi - atan(py_/px_);
|
105 |
|
|
if (px<0 && py<0)
|
106 |
|
|
phi = -pi + atan(py_/px_);
|
107 |
|
|
|
108 |
|
|
|
109 |
|
|
if (phi<-pi || phi>pi)
|
110 |
|
|
{
|
111 |
|
|
cout << "ERROR: phi out of range\n";
|
112 |
|
|
return 0;
|
113 |
|
|
}
|
114 |
|
|
return phi;
|
115 |
|
|
}
|
116 |
|
|
int FindTrueLepMatchedToThisRECLep(Electron* recMu,
|
117 |
|
|
vector<mcParticle*> genMus, double deltar, bool requiresamecharge)
|
118 |
|
|
{
|
119 |
|
|
int trueMuonIndex = -1; // -1 : unmatched
|
120 |
|
|
double dr = 999;
|
121 |
|
|
for(int i=0; i<(int)genMus.size(); i++)
|
122 |
|
|
{
|
123 |
|
|
if(abs(genMus[i]->id())!=11) continue;
|
124 |
|
|
if(requiresamecharge){
|
125 |
|
|
if(genMus[i]->q()!=recMu->q()) continue;
|
126 |
|
|
}
|
127 |
|
|
double tmpdr = GetDeltaR(recMu->eta(),genMus[i]->eta(),recMu->phi(),genMus[i]->phi());
|
128 |
|
|
if (tmpdr < dr ) { dr = tmpdr; trueMuonIndex = i ; }
|
129 |
|
|
}
|
130 |
|
|
if( dr < deltar ) return trueMuonIndex;
|
131 |
|
|
return -1;
|
132 |
|
|
}
|
133 |
|
|
int FindTrueLepMatchedToThisRECLep(Muon* recMu,
|
134 |
|
|
vector<mcParticle*> genMus, double deltar, bool requiresamecharge)
|
135 |
|
|
{
|
136 |
|
|
int trueMuonIndex = -1; // -1 : unmatched
|
137 |
|
|
double dr = 999;
|
138 |
|
|
for(int i=0; i<(int)genMus.size(); i++)
|
139 |
|
|
{
|
140 |
|
|
// if(abs(genMus[i]->id())!=13) continue;
|
141 |
|
|
if(requiresamecharge){
|
142 |
|
|
if(genMus[i]->q()!=recMu->q()) continue;
|
143 |
|
|
}
|
144 |
|
|
double tmpdr = GetDeltaR(recMu->eta(),genMus[i]->eta(),recMu->phi(),genMus[i]->phi());
|
145 |
|
|
if (tmpdr < dr ) { dr = tmpdr; trueMuonIndex = i ; }
|
146 |
|
|
}
|
147 |
|
|
if( dr < deltar ) return trueMuonIndex;
|
148 |
|
|
return -1;
|
149 |
|
|
}
|
150 |
|
|
int FindRecLepMatchedToThisGenLep(mcParticle* genMu,
|
151 |
|
|
vector<Electron*> recMus, int used, double deltar, bool requiresamecharge)
|
152 |
|
|
{
|
153 |
|
|
int recMuonIndex = -1; // -1 : unmatched
|
154 |
|
|
double dr = 999;
|
155 |
|
|
for(int i=0; i<(int)recMus.size(); i++)
|
156 |
|
|
{
|
157 |
|
|
if(i==used) continue;
|
158 |
|
|
if(abs(recMus[i]->id())!=11) continue;
|
159 |
|
|
if(requiresamecharge){
|
160 |
|
|
if(recMus[i]->q()!=genMu->q()) continue;
|
161 |
|
|
}
|
162 |
|
|
double tmpdr = GetDeltaR(genMu->eta(),recMus[i]->eta(),genMu->phi(),recMus[i]->phi());
|
163 |
|
|
if (tmpdr < dr ) { dr = tmpdr; recMuonIndex = i ; }
|
164 |
|
|
}
|
165 |
|
|
if( dr < deltar ) return recMuonIndex;
|
166 |
|
|
return -1;
|
167 |
|
|
}
|
168 |
|
|
int FindRecLepMatchedToThisGenLep(mcParticle* genMu,
|
169 |
|
|
vector<Muon*> recMus, int used, double deltar, bool requiresamecharge)
|
170 |
|
|
{
|
171 |
|
|
int recMuonIndex = -1; // -1 : unmatched
|
172 |
|
|
double dr = 999;
|
173 |
|
|
for(int i=0; i<(int)recMus.size(); i++)
|
174 |
|
|
{
|
175 |
|
|
if(i==used) continue;
|
176 |
|
|
if(abs(recMus[i]->id())!=13) continue;
|
177 |
|
|
if(requiresamecharge){
|
178 |
|
|
if(recMus[i]->q()!=genMu->q()) continue;
|
179 |
|
|
}
|
180 |
|
|
double tmpdr = GetDeltaR(genMu->eta(),recMus[i]->eta(),genMu->phi(),recMus[i]->phi());
|
181 |
|
|
if (tmpdr < dr ) { dr = tmpdr; recMuonIndex = i ; }
|
182 |
|
|
}
|
183 |
|
|
if( dr < deltar ) return recMuonIndex;
|
184 |
|
|
return -1;
|
185 |
|
|
}
|
186 |
|
|
int FindRecLepMatchedToThisGenLep(mcParticle* genMu,
|
187 |
|
|
vector<Electron*> recMus, double deltar, bool requiresamecharge)
|
188 |
|
|
{
|
189 |
|
|
int recMuonIndex = -1; // -1 : unmatched
|
190 |
|
|
double dr = 999;
|
191 |
|
|
for(int i=0; i<(int)recMus.size(); i++)
|
192 |
|
|
{
|
193 |
|
|
if(abs(recMus[i]->id())!=11) continue;
|
194 |
|
|
if(requiresamecharge){
|
195 |
|
|
if(recMus[i]->q()!=genMu->q()) continue;
|
196 |
|
|
}
|
197 |
|
|
double tmpdr = GetDeltaR(genMu->eta(),recMus[i]->eta(),genMu->phi(),recMus[i]->phi());
|
198 |
|
|
if (tmpdr < dr ) { dr = tmpdr; recMuonIndex = i ; }
|
199 |
|
|
}
|
200 |
|
|
if( dr < deltar ) return recMuonIndex;
|
201 |
|
|
return -1;
|
202 |
|
|
}
|
203 |
|
|
int FindRecLepMatchedToThisGenLep(mcParticle* genMu,
|
204 |
|
|
vector<Muon*> recMus, double deltar, bool requiresamecharge)
|
205 |
|
|
{
|
206 |
|
|
int recMuonIndex = -1; // -1 : unmatched
|
207 |
|
|
double dr = 999;
|
208 |
|
|
for(int i=0; i<(int)recMus.size(); i++)
|
209 |
|
|
{
|
210 |
|
|
if(abs(recMus[i]->id())!=13) continue;
|
211 |
|
|
if(requiresamecharge){
|
212 |
|
|
if(recMus[i]->q()!=genMu->q()) continue;
|
213 |
|
|
}
|
214 |
|
|
double tmpdr = GetDeltaR(genMu->eta(),recMus[i]->eta(),genMu->phi(),recMus[i]->phi());
|
215 |
|
|
if (tmpdr < dr ) { dr = tmpdr; recMuonIndex = i ; }
|
216 |
|
|
}
|
217 |
|
|
if( dr < deltar ) return recMuonIndex;
|
218 |
|
|
return -1;
|
219 |
|
|
}
|
220 |
|
|
void PrintRecMuonInfo(vector<Muon*> recMu){
|
221 |
|
|
printf(" muonSize = %3d\n",recMu.size());
|
222 |
|
|
printf(" pt eta phi d0 |calEm calHm RelIso Chi2N NHits Charge |\n");
|
223 |
|
|
for(Int_t i=0; i<(int)recMu.size(); i++){
|
224 |
|
|
printf(" %10.3f %10.3f %6.3f %6.3f | %5.1f %8.2f %10.1f %8.3f %6.1f %3d |\n",
|
225 |
|
|
recMu[i]->pt(),
|
226 |
|
|
recMu[i]->eta(),
|
227 |
|
|
recMu[i]->phi(),
|
228 |
|
|
recMu[i]->d0(),
|
229 |
|
|
recMu[i]->calEm(),
|
230 |
|
|
recMu[i]->calHm(),
|
231 |
|
|
recMu[i]->isolSumNorm(),
|
232 |
|
|
recMu[i]->Chi2N(),
|
233 |
|
|
recMu[i]->nHit(),
|
234 |
|
|
recMu[i]->q()
|
235 |
|
|
);
|
236 |
|
|
|
237 |
|
|
}
|
238 |
|
|
}
|
239 |
|
|
void PrintRecMuonInfo(vector<Muon*> recMu,ofstream& fout){
|
240 |
|
|
|
241 |
|
|
fout<<"rec muon Size "<<setw(8)<<recMu.size()<<endl;
|
242 |
|
|
fout<<" pt eta phi d0 |calEm calHm RelIso Chi2N NHits Charge "<<endl;
|
243 |
|
|
for(Int_t i=0; i<(int)recMu.size(); i++){
|
244 |
|
|
fout<<setw(6)<< recMu[i]->pt() << " "
|
245 |
|
|
<<setw(6)<< recMu[i]->eta()<< " "
|
246 |
|
|
<<setw(6)<< recMu[i]->phi()<< " "
|
247 |
|
|
<<setw(6)<< recMu[i]->d0()<< " | "
|
248 |
|
|
<<setw(6)<< recMu[i]->calEm()<< " "
|
249 |
|
|
<<setw(6)<< recMu[i]->calHm()<< " "
|
250 |
|
|
<<setw(6)<< recMu[i]->isolSumNorm()<< " "
|
251 |
|
|
<<setw(6)<< recMu[i]->Chi2N()<< " "
|
252 |
|
|
<<setw(6)<< recMu[i]->nHit()<<" "
|
253 |
|
|
<<setw(6)<< recMu[i]->q()<<endl;
|
254 |
|
|
}
|
255 |
|
|
}
|
256 |
|
|
float DeltaPhi(float v1, float v2)
|
257 |
|
|
{ // Computes the correctly normalized phi difference
|
258 |
|
|
// v1, v2 = phi of object 1 and 2
|
259 |
|
|
float diff = fabs(v2 - v1);
|
260 |
|
|
float corr = 2*acos(-1.) - diff;
|
261 |
|
|
if (diff < acos(-1.)){ return diff;} else { return corr;}
|
262 |
|
|
}
|
263 |
|
|
float GetDeltaR(float eta1, float eta2, float phi1, float phi2)
|
264 |
|
|
{ // Computes the DeltaR of two objects from their eta and phi values
|
265 |
|
|
return sqrt( (eta1-eta2)*(eta1-eta2)
|
266 |
|
|
+ DeltaPhi(phi1, phi2)*DeltaPhi(phi1, phi2) );
|
267 |
|
|
}
|
268 |
|
|
void PrintRecElectronInfo(vector<Electron*> recMu,ofstream& fout){
|
269 |
|
|
|
270 |
|
|
fout<<"rec electron Size "<<setw(8)<<recMu.size()<<endl;
|
271 |
|
|
fout<<" pt eta phi d0 | RelIso Chi2N NHits Charge "<<endl;
|
272 |
|
|
for(Int_t i=0; i<(int)recMu.size(); i++){
|
273 |
|
|
fout<<setw(6)<< recMu[i]->pt() << " "
|
274 |
|
|
<<setw(6)<< recMu[i]->eta()<< " "
|
275 |
|
|
<<setw(6)<< recMu[i]->phi()<< " "
|
276 |
|
|
<<setw(6)<< recMu[i]->d0()<< " | "
|
277 |
|
|
<<setw(6)<< recMu[i]->isolSumNorm()<< " "
|
278 |
|
|
<<setw(6)<< recMu[i]->Chi2N()<< " "
|
279 |
|
|
<<setw(6)<< recMu[i]->nHit()<<" "
|
280 |
|
|
<<setw(6)<< recMu[i]->q()<<endl;
|
281 |
|
|
}
|
282 |
|
|
}
|
283 |
|
|
int FindTrueElectronMatchedToThisRECElectron(Electron* recMu,
|
284 |
|
|
vector<mcParticle*> genMus)
|
285 |
|
|
{
|
286 |
|
|
int trueElectronIndex = -1; // -1 : unmatched
|
287 |
|
|
double dr = 999;
|
288 |
|
|
for(int i=0; i<(int)genMus.size(); i++)
|
289 |
|
|
{
|
290 |
|
|
if(abs(genMus[i]->id())!=11) continue;
|
291 |
|
|
double tmpdr = GetDeltaR(recMu->eta(),genMus[i]->eta(),recMu->phi(),genMus[i]->phi());
|
292 |
|
|
if (tmpdr < dr ) { dr = tmpdr; trueElectronIndex = i ; }
|
293 |
|
|
}
|
294 |
|
|
// FIXME here is a hard-coded cut value for matching
|
295 |
|
|
if( dr < 0.3 ) return trueElectronIndex;
|
296 |
|
|
return -1;
|
297 |
|
|
}
|
298 |
|
|
void PrintBriefGenChain(vector<mcParticle*> genMu, ofstream& fout)
|
299 |
|
|
{
|
300 |
|
|
fout<<"gen muon Size "<<setw(8)<<genMu.size()<<endl;
|
301 |
|
|
fout<<"Charge parentID | pt eta phi"<<endl;
|
302 |
|
|
for(Int_t i=0; i<(int)genMu.size(); i++){
|
303 |
|
|
fout<<setw(6)<< genMu[i]->q()<<" "
|
304 |
|
|
<<setw(6)<< genMu[i]->getParentPDG()<<" | "
|
305 |
|
|
<<setw(6)<< genMu[i]->pt() << " "
|
306 |
|
|
<<setw(6)<< genMu[i]->eta()<< " "
|
307 |
|
|
<<setw(6)<< genMu[i]->phi()<< " "<<endl;
|
308 |
|
|
}
|
309 |
|
|
}
|
310 |
|
|
char * asctime_mine(const struct tm *timeptr) {
|
311 |
|
|
static char mon_name[12][4] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
|
312 |
|
|
static char result[19];
|
313 |
|
|
// Result: "06-Dec-09 09:45:37"
|
314 |
|
|
sprintf(result, "%.2d-%.3s-%.2d %.2d:%.2d:%.2d",
|
315 |
|
|
timeptr->tm_mday,
|
316 |
|
|
mon_name[timeptr->tm_mon],
|
317 |
|
|
(timeptr->tm_year)%100,
|
318 |
|
|
timeptr->tm_hour, timeptr->tm_min, timeptr->tm_sec);
|
319 |
|
|
return result;
|
320 |
|
|
}
|
321 |
|
|
|
322 |
|
|
bool filterRun(int run, int lumi){
|
323 |
|
|
if(
|
324 |
|
|
(run!=123596) &&
|
325 |
|
|
(run!=123615) &&
|
326 |
|
|
(run!=123732) &&
|
327 |
|
|
(run!=123815) &&
|
328 |
|
|
(run!=123818) &&
|
329 |
|
|
(run!=123906) &&
|
330 |
|
|
(run!=123908) &&
|
331 |
|
|
(run!=123987) &&
|
332 |
|
|
(run!=124006) &&
|
333 |
|
|
(run!=124009) &&
|
334 |
|
|
(run!=124020) &&
|
335 |
|
|
(run!=124022) &&
|
336 |
|
|
(run!=124023) &&
|
337 |
|
|
(run!=124024) &&
|
338 |
|
|
(run!=124025) &&
|
339 |
|
|
(run!=124027) &&
|
340 |
|
|
(run!=124030) &&
|
341 |
|
|
(run!=124120) &&
|
342 |
|
|
(run!=124230)
|
343 |
|
|
)
|
344 |
|
|
return true;
|
345 |
|
|
|
346 |
|
|
bool accepted = false;
|
347 |
|
|
|
348 |
|
|
if (run==123151) {
|
349 |
|
|
if(lumi>=3 && lumi<=19)
|
350 |
|
|
accepted=true;
|
351 |
|
|
} else if (run==123596) {
|
352 |
|
|
if (/*(lumi>=4 && lumi<=26) ||*/ //pixel timing scan
|
353 |
|
|
(lumi>=69 && lumi<=144) )
|
354 |
|
|
accepted=true;
|
355 |
|
|
} else if (run==123615) {
|
356 |
|
|
if(lumi>=71)
|
357 |
|
|
accepted=true;
|
358 |
|
|
} else if (run==123732) {
|
359 |
|
|
if(lumi>=62 && lumi<=112) //though phys bit starting 57
|
360 |
|
|
accepted=true;
|
361 |
|
|
} else if (run==123815) {
|
362 |
|
|
if(lumi>=8 && lumi<=16)
|
363 |
|
|
accepted=true;
|
364 |
|
|
} else if (run==123818) {
|
365 |
|
|
if(lumi>=2 && lumi<=18) //RunRegistry says lumi scan starting at lumi 19 until 42
|
366 |
|
|
accepted=true;
|
367 |
|
|
} else if (run==123906) {
|
368 |
|
|
if(lumi>=17 && lumi<=28)
|
369 |
|
|
accepted=true;
|
370 |
|
|
} else if (run==123908) {
|
371 |
|
|
if(lumi>=2 && lumi<=13)
|
372 |
|
|
accepted=true;
|
373 |
|
|
/* else if (run==123987) //3T run
|
374 |
|
|
if(lumi>=1 && lumi<=21)
|
375 |
|
|
accepted=true;*/
|
376 |
|
|
} else if (run==124006) {
|
377 |
|
|
if(lumi>=1 && lumi<=6) //though Phys bit set from lumi 6
|
378 |
|
|
accepted=true;
|
379 |
|
|
} else if (run==124009) { //lumi scan through 29-63
|
380 |
|
|
if(lumi>=1 && lumi<=68)
|
381 |
|
|
accepted=true;
|
382 |
|
|
} else if (run==124020) {
|
383 |
|
|
if(lumi>=12 && lumi<=94)
|
384 |
|
|
accepted=true;
|
385 |
|
|
} else if (run==124022) {
|
386 |
|
|
if(lumi>=65 && lumi<=160) //lumi scan through 161-183
|
387 |
|
|
accepted=true;
|
388 |
|
|
} else if (run==124023) {
|
389 |
|
|
if(lumi>=41 && lumi<=96)
|
390 |
|
|
accepted=true;
|
391 |
|
|
} else if (run==124024) {
|
392 |
|
|
if(lumi>=2 && lumi<=83)
|
393 |
|
|
accepted=true;
|
394 |
|
|
} else if (run==124025) {
|
395 |
|
|
if(lumi>=3 && lumi<=13)
|
396 |
|
|
accepted=true;
|
397 |
|
|
} else if (run==124027) {
|
398 |
|
|
if(lumi>=23 && lumi<=39)
|
399 |
|
|
accepted=true;
|
400 |
|
|
} else if (run==124030) {
|
401 |
|
|
if(lumi>=1 && lumi<=32)
|
402 |
|
|
accepted=true;
|
403 |
|
|
} else if (run==124120) { //2.36 TeV run
|
404 |
|
|
if(lumi>=1)
|
405 |
|
|
accepted=true;
|
406 |
|
|
} else if (run==124230) {
|
407 |
|
|
if(lumi>=26 && lumi<=68) //paused at 47, resumed 48
|
408 |
|
|
accepted=true;
|
409 |
|
|
}
|
410 |
|
|
}
|
411 |
|
|
////////////////////////////////////////////////////////////////////////
|
412 |
|
|
//calculates theta from ONIA
|
413 |
|
|
//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/HeavyFlavorAnalysis/Onia2MuMu/src/Onia2MuMu.cc?hideattic=0&revision=1.37&view=markup
|
414 |
|
|
////////////////////////////////////////////////////////////////////////
|
415 |
|
|
double GetTheta( TLorentzVector a, TLorentzVector b) {
|
416 |
|
|
TLorentzVector c = a+b;
|
417 |
|
|
TVector3 bv=c.BoostVector();
|
418 |
|
|
a.Boost(-bv);
|
419 |
|
|
b.Boost(-bv);
|
420 |
|
|
double theta = c.Vect().Angle(a.Vect());
|
421 |
|
|
return theta;
|
422 |
|
|
}
|
423 |
|
|
double TransverseMass(double metx, double mety, TLorentzVector lepton ){
|
424 |
|
|
|
425 |
|
|
// Candidates have a mt() function which computes the tranverse mass from E & pz.
|
426 |
|
|
// As MET does not have pz information... WMuNuCandidates have an alternative function to compute the mt quantity
|
427 |
|
|
// used in the WMuNu Inclusive analysis just from px, py
|
428 |
|
|
//
|
429 |
|
|
|
430 |
|
|
TLorentzVector neutrino(metx, mety, 0., sqrt(metx*metx+mety*mety));
|
431 |
|
|
|
432 |
|
|
double eT=0;
|
433 |
|
|
eT=lepton.Pt()+neutrino.Pt();
|
434 |
|
|
|
435 |
|
|
TLorentzVector w = lepton + neutrino;
|
436 |
|
|
|
437 |
|
|
double wpx=w.Px(); double wpy=w.Py();
|
438 |
|
|
double mt = eT*eT - wpx*wpx - wpy*wpy;
|
439 |
|
|
mt = (mt>0) ? sqrt(mt) : 0;
|
440 |
|
|
return mt;
|
441 |
|
|
}
|
442 |
|
|
|