ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/mschen/SusyAnalysis/code/utils.cc
Revision: 1.1
Committed: Mon Mar 28 09:23:54 2011 UTC (14 years, 1 month ago) by mschen
Content type: text/plain
Branch: MAIN
CVS Tags: V2010_data_analysis_effModelInPaper, V2010_data_analysis, HEAD
Log Message:
2010 same sign analysis codeing

File Contents

# User Rev Content
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