ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/ZHOfflineAnalysis.C
Revision: 1.14
Committed: Thu Jun 13 15:46:29 2013 UTC (11 years, 11 months ago) by dasu
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.13: +0 -0 lines
State: FILE REMOVED
Error occurred while calculating annotation data.
Log Message:
Pruned to the minimum needed for future

File Contents

# Content
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 #include <vector>
5 #include <utility>
6 #include <map>
7 #include <string>
8 #include <math.h>
9 using namespace std;
10 #include "TH1F.h"
11 #include "TTree.h"
12 #include "TFile.h"
13 #include "TSystem.h"
14 #include "UltraFastSim.h"
15 #include "LinkDef.h"
16 #include "TMath.h" //M_PI is in TMath
17
18 double deltaR(double eta1, double eta2, double phi1, double phi2)
19 {
20 double dphi = fabs(phi1-phi2);
21 if(dphi > M_PI)
22 {
23 dphi = (2*M_PI - dphi);
24 }
25 double deta = fabs(eta1-eta2);
26 return sqrt(dphi*dphi + deta*deta);
27 }
28 double deltaphi(double phi1, double phi2)
29 {
30 double dphi = fabs(phi1-phi2);
31 if(dphi > M_PI)
32 {
33 dphi = (2*M_PI - dphi);
34 }
35 return dphi;
36 }
37 double CorrectPhi(double phi, double x, double y)
38 {
39 if(phi>0 && x<0 && y<0)phi+=M_PI;
40 if(phi<0 && x>0 && y<0)phi=(2*M_PI-fabs(phi));
41 if(phi<0 && x<0 && y>0)phi=(M_PI-fabs(phi));
42 if(y<0)phi-=(2*M_PI);//without this you have 0<phi<2pi, but with this
43 return phi;//you get -pi<phi<pi
44 }
45
46 int main(){
47
48 vector<TParticle> muons;
49 vector<TLorentzVector> bjets[6];
50
51
52 int AnalyzedEvents=0;
53 int LeptonPairPreSelection=0;
54 int LeptonPairSelection=0;
55 int ZMassWindow=0;
56 int BJetPairPreSelection[6]={0};
57 int BJetPairSelection[6]={0};
58 int ZptSelection[6]={0};
59 int HptSelection[6]={0};
60 int DphiSelection[6]={0};
61 int HMassWindow[6]={0};
62
63 double Weight=1.;
64
65 //================================================
66 //event loop
67
68 ifstream ConfigFile;
69 ConfigFile.open("Config.txt");
70 if(!ConfigFile){cout<<"Unable to read the config file";exit(1);}
71
72 string OutputRootFile;
73 double CrossSection;
74 double Lumi;
75
76 ConfigFile>>OutputRootFile;
77 ConfigFile>>CrossSection;
78 ConfigFile>>Lumi;
79
80 //file should be created before tree to avoid memory resident trees
81 TFile fout(OutputRootFile.c_str(),"recreate");
82
83 TH1F* hptmuons = new TH1F("Muons Pt","Muons Pt",300,0.,300.);
84 TH1F* hetamuons = new TH1F("Muons Eta","Muons Eta",100,-5.,5.);
85 TH1F* hphimuons = new TH1F("Muons Phi","Muons Phi",100,-7.,7.);
86 TH1F* hzinvmass= new TH1F("Z invmass","Z invmass",300,0.,300.);
87 TH1F* hphiz= new TH1F("Z Phi","Z Phi",100,-7.,7.);
88
89 TH1F* hbjet_mult[6];
90 hbjet_mult[0] = new TH1F("B-jets Mult. Loose StdGeom","B-jets Mult. Loose StdGeom",10,0.,10.);
91 hbjet_mult[1] = new TH1F("B-jets Mult. Loose Phase 1","B-jets Mult. Loose Phase 1",10,0.,10.);
92 hbjet_mult[2] = new TH1F("B-jets Mult. Medium StdGeom","B-jets Mult. Medium StdGeom",10,0.,10.);
93 hbjet_mult[3] = new TH1F("B-jets Mult. Medium Phase 1","B-jets Mult. Medium Phase 1",10,0.,10.);
94 hbjet_mult[4] = new TH1F("B-jets Mult. Tight StdGeom","B-jets Mult. Tight StdGeom",10,0.,10.);
95 hbjet_mult[5] = new TH1F("B-jets Mult. Tight Phase 1","B-jets Mult. Tight Phase 1",10,0.,10.);
96
97 TH1F* hetbjets[6];
98 hetbjets[0] = new TH1F("B-jets Et Loose StdGeom","B-jets Et Loose StdGeom",300,0.,300.);
99 hetbjets[1] = new TH1F("B-jets Et Loose Phase 1","B-jets Et Loose Phase 1",300,0.,300.);
100 hetbjets[2] = new TH1F("B-jets Et Medium StdGeom","B-jets Et Medium StdGeom",300,0.,300.);
101 hetbjets[3] = new TH1F("B-jets Et Medium Phase 1","B-jets Et Medium Phase 1",300,0.,300.);
102 hetbjets[4] = new TH1F("B-jets Et Tight StdGeom","B-jets Et Tight StdGeom",300,0.,300.);
103 hetbjets[5] = new TH1F("B-jets Et Tight Phase 1","B-jets Et Tight Phase 1",300,0.,300.);
104
105 TH1F* hetabjets[6];
106 hetabjets[0] = new TH1F("B-jets Eta Loose StdGeom","B-jets Eta Loose StdGeom",100,-5.,5.);
107 hetabjets[1] = new TH1F("B-jets Eta Loose Phase 1","B-jets Eta Loose Phase 1",100,-5.,5.);
108 hetabjets[2] = new TH1F("B-jets Eta Medium StdGeom","B-jets Eta Medium StdGeom",100,-5.,5.);
109 hetabjets[3] = new TH1F("B-jets Eta Medium Phase 1","B-jets Eta Medium Phase 1",100,-5.,5.);
110 hetabjets[4] = new TH1F("B-jets Eta Tight StdGeom","B-jets Eta Tight StdGeom",100,-5.,5.);
111 hetabjets[5] = new TH1F("B-jets Eta Tight Phase 1","B-jets Eta Tight Phase 1",100,-5.,5.);
112
113 TH1F* hphibjets[6];
114 hphibjets[0] = new TH1F("B-jets Phi Loose StdGeom","B-jets Phi Loose StdGeom",100,-7.,7.);
115 hphibjets[1] = new TH1F("B-jets Phi Loose Phase 1","B-jets Phi Loose Phase 1",100,-7.,7.);
116 hphibjets[2] = new TH1F("B-jets Phi Medium StdGeom","B-jets Phi Medium StdGeom",100,-7.,7.);
117 hphibjets[3] = new TH1F("B-jets Phi Medium Phase 1","B-jets Phi Medium Phase 1",100,-7.,7.);
118 hphibjets[4] = new TH1F("B-jets Phi Tight StdGeom","B-jets Phi Tight StdGeom",100,-7.,7.);
119 hphibjets[5] = new TH1F("B-jets Phi Tight Phase 1","B-jets Phi Tight Phase 1",100,-7.,7.);
120
121 TH1F* hptz[6];
122 hptz[0] = new TH1F("Z Pt Loose StdGeom","Z Pt Loose StdGeom",300,0.,300.);
123 hptz[1] = new TH1F("Z Pt Loose Phase 1","Z Pt Loose Phase 1",300,0.,300.);
124 hptz[2] = new TH1F("Z Pt Medium StdGeom","Z Pt Medium StdGeom",300,0.,300.);
125 hptz[3] = new TH1F("Z Pt Medium Phase 1","Z Pt Medium Phase 1",300,0.,300.);
126 hptz[4] = new TH1F("Z Pt Tight StdGeom","Z Pt Tight StdGeom",300,0.,300.);
127 hptz[5] = new TH1F("Z Pt Tight Phase 1","Z Pt Tight Phase 1",300,0.,300.);
128
129 TH1F* hpth[6];
130 hpth[0] = new TH1F("H Pt Loose StdGeom","H Pt Loose StdGeom",300,0.,300.);
131 hpth[1] = new TH1F("H Pt Loose Phase 1","H Pt Loose Phase 1",300,0.,300.);
132 hpth[2] = new TH1F("H Pt Medium StdGeom","H Pt Medium StdGeom",300,0.,300.);
133 hpth[3] = new TH1F("H Pt Medium Phase 1","H Pt Medium Phase 1",300,0.,300.);
134 hpth[4] = new TH1F("H Pt Tight StdGeom","H Pt Tight StdGeom",300,0.,300.);
135 hpth[5] = new TH1F("H Pt Tight Phase 1","H Pt Tight Phase 1",300,0.,300.);
136
137 TH1F* hhinvmass[6];
138 hhinvmass[0] = new TH1F("H invmass preselection Loose StdGeom","H invmass preselection Loose StdGeom",300,0.,300.);
139 hhinvmass[1] = new TH1F("H invmass preselection Loose Phase 1","H invmass preselection Loose Phase 1",300,0.,300.);
140 hhinvmass[2] = new TH1F("H invmass preselection Medium StdGeom","H invmass preselection Medium StdGeom",300,0.,300.);
141 hhinvmass[3] = new TH1F("H invmass preselection Medium Phase 1","H invmass preselection Medium Phase 1",300,0.,300.);
142 hhinvmass[4] = new TH1F("H invmass preselection Tight StdGeom","H invmass preselection Tight StdGeom",300,0.,300.);
143 hhinvmass[5] = new TH1F("H invmass preselection Tight Phase 1","H invmass preselection Tight Phase 1",300,0.,300.);
144
145 TH1F* hhinvmass2[6];
146 hhinvmass2[0] = new TH1F("H invmass Loose StdGeom","H invmass Loose StdGeom",300,0.,300.);
147 hhinvmass2[1] = new TH1F("H invmass Loose Phase 1","H invmass Loose Phase 1",300,0.,300.);
148 hhinvmass2[2] = new TH1F("H invmass Medium StdGeom","H invmass Medium StdGeom",300,0.,300.);
149 hhinvmass2[3] = new TH1F("H invmass Medium Phase 1","H invmass Medium Phase 1",300,0.,300.);
150 hhinvmass2[4] = new TH1F("H invmass Tight StdGeom","H invmass Tight StdGeom",300,0.,300.);
151 hhinvmass2[5] = new TH1F("H invmass Tight Phase 1","H invmass Tight Phase 1",300,0.,300.);
152
153 TH1F* hphih[6];
154 hphih[0] = new TH1F("H Phi Loose StdGeom","H Phi Loose StdGeom",100,-7.,7.);
155 hphih[1] = new TH1F("H Phi Loose Phase 1","H Phi Loose Phase 1",100,-7.,7.);
156 hphih[2] = new TH1F("H Phi Medium StdGeom","H Phi Medium StdGeom",100,-7.,7.);
157 hphih[3] = new TH1F("H Phi Medium Phase 1","H Phi Medium Phase 1",100,-7.,7.);
158 hphih[4] = new TH1F("H Phi Tight StdGeom","H Phi Tight StdGeom",100,-7.,7.);
159 hphih[5] = new TH1F("H Phi Tight Phase 1","H Phi Tight Phase 1",100,-7.,7.);
160
161 TH1F* hdphiZH[6];
162 hdphiZH[0] = new TH1F("Dphi(Z,H) Loose StdGeom","Dphi(Z,H) Loose StdGeom",100,-7.,7.);
163 hdphiZH[1] = new TH1F("Dphi(Z,H) Loose Phase 1","Dphi(Z,H) Loose Phase 1",100,-7.,7.);
164 hdphiZH[2] = new TH1F("Dphi(Z,H) Medium StdGeom","Dphi(Z,H) Medium StdGeom",100,-7.,7.);
165 hdphiZH[3] = new TH1F("Dphi(Z,H) Medium Phase 1","Dphi(Z,H) Medium Phase 1",100,-7.,7.);
166 hdphiZH[4] = new TH1F("Dphi(Z,H) Tight StdGeom","Dphi(Z,H) Tight StdGeom",100,-7.,7.);
167 hdphiZH[5] = new TH1F("Dphi(Z,H) Tight Phase 1","Dphi(Z,H) Tight Phase 1",100,-7.,7.);
168
169
170 //Just to know the total statistics of the root files
171 //This is needed for the event normalization
172 string filename;
173 vector<string> fns;//file name container
174 int SumEvents=0;
175 while(ConfigFile>>filename && filename!="EOF"){
176 fns.push_back(filename);
177 cout<<"checking "<<filename<<endl;
178 TFile *file=TFile::Open(filename.c_str());
179 TTree *rootTree = (TTree*)file->Get("UltraFastSim");
180 int nev = int(rootTree->GetEntries());
181 cout<<"It has "<<nev<<" events."<<endl;
182 SumEvents+=nev;
183 delete rootTree;
184 delete file;
185 }
186
187 cout<<"Will run on total "<<SumEvents<<" events ... "<<endl;
188 Weight=CrossSection*Lumi/SumEvents;
189 cout<<"Event weight is "<<Weight<<endl;
190
191 // ========== >>>> Loop over files starts here
192 for(vector<string>::iterator iterstr=fns.begin();iterstr!=fns.end();iterstr++){
193
194 cout<<"Reading "<<*iterstr<<endl;
195 TFile *file=TFile::Open((*iterstr).c_str());
196 UltraFastSim *ufs=new UltraFastSim();
197 TTree *rootTree = (TTree*)file->Get("UltraFastSim");
198 int nev = int(rootTree->GetEntries());
199 cout<<"number of entries is : "<<nev<<endl;
200 TBranch *branch = rootTree->GetBranch("UltraFastSim");
201 branch->SetAddress(&ufs);
202
203
204 for(int i = 0; i < nev; i++){
205
206 AnalyzedEvents++;
207 if(!(AnalyzedEvents%1000))cout<<"event number "<<AnalyzedEvents<<endl;
208
209 rootTree->GetEvent(i);
210 muons=ufs->muonList();
211 bjets[0]=ufs->bJetListLooseStdGeom();
212 bjets[1]=ufs->bJetListLoose();
213 bjets[2]=ufs->bJetListMediumStdGeom();
214 bjets[3]=ufs->bJetListMedium();
215 bjets[4]=ufs->bJetListTightStdGeom();
216 bjets[5]=ufs->bJetListTight();
217
218 for(vector<TParticle>::iterator itm=muons.begin();itm!=muons.end();itm++){
219 hptmuons->Fill(itm->Pt(),Weight);
220 hetamuons->Fill(itm->Eta(),Weight);
221 hphimuons->Fill(itm->Phi(),Weight);
222 }
223 for(int i=0;i<6;i++){
224 for(vector<TLorentzVector>::iterator itj=bjets[i].begin();itj!=bjets[i].end();itj++){
225 hetbjets[i]->Fill(itj->Pt(),Weight);
226 hetabjets[i]->Fill(itj->Eta(),Weight);
227 hphibjets[i]->Fill(itj->Phi(),Weight);
228 }
229 }
230
231 double Zinvmass=0;
232 double Zpt=0;
233 double Zphi=0;
234
235 double Hinvmass[6]={0};
236 double Hpt[6]={0};
237 double Hpx[6]={0};
238 double Hpy[6]={0};
239 double Hphi[6]={0};
240 double DphiZH[6]={0};
241
242 if(muons.size()>1){
243
244 LeptonPairPreSelection++;
245
246 if(muons[0].Pt()>20. && muons[1].Pt()>20.){
247
248 LeptonPairSelection++;
249
250 for(int i=0;i<6;i++)hbjet_mult[i]->Fill(bjets[i].size());
251
252 Zinvmass=(pow((muons[0].Energy()+muons[1].Energy()),2)-
253 pow((muons[0].Px()+muons[1].Px()),2)-
254 pow((muons[0].Py()+muons[1].Py()),2)-
255 pow((muons[0].Pz()+muons[1].Pz()),2));
256 if(Zinvmass>0)Zinvmass=sqrt(Zinvmass);
257 else Zinvmass=0;
258 Zpt=sqrt(pow((muons[0].Px()+muons[1].Px()),2)+pow((muons[0].Py()+muons[1].Py()),2));
259
260 double Zpx=0;
261 double Zpy=0;
262 Zpx=muons[0].Px()+muons[1].Px();
263 Zpy=muons[0].Py()+muons[1].Py();
264 Zphi=atan(Zpy/Zpx);
265 Zphi=CorrectPhi(Zphi,Zpx,Zpy);
266
267 hzinvmass->Fill(Zinvmass,Weight);
268 hphiz->Fill(Zphi,Weight);
269
270 if(Zinvmass>70. && Zinvmass<110.){
271
272 ZMassWindow++;
273
274 //=======>>>>> loop over b-jet algo, and tk. geometry
275 for(int i=0;i<6;i++){
276
277 if(bjets[i].size()>1){
278
279 BJetPairPreSelection[i]++;
280
281 if(bjets[i][0].Pt()>30. && bjets[i][1].Pt()>30.){
282
283 BJetPairSelection[i]++;
284
285 Hinvmass[i]=(pow((bjets[i][0].Energy()+bjets[i][1].Energy()),2)-
286 pow((bjets[i][0].Px()+bjets[i][1].Px()),2)-
287 pow((bjets[i][0].Py()+bjets[i][1].Py()),2)-
288 pow((bjets[i][0].Pz()+bjets[i][1].Pz()),2));
289 if(Hinvmass[i]>0)Hinvmass[i]=sqrt(Hinvmass[i]);
290 else Hinvmass[i]=0;
291
292 Hpt[i]=sqrt(pow((bjets[i][0].Px()+bjets[i][1].Px()),2)
293 +pow((bjets[i][0].Py()+bjets[i][1].Py()),2));
294
295 Hpx[i]=bjets[i][0].Px()+bjets[i][1].Px();
296 Hpy[i]=bjets[i][0].Py()+bjets[i][1].Py();
297 Hphi[i]=atan(Hpy[i]/Hpx[i]);
298 Hphi[i]=CorrectPhi(Hphi[i],Hpx[i],Hpy[i]);
299
300 hhinvmass[i]->Fill(Hinvmass[i],Weight);
301 hphih[i]->Fill(Hphi[i],Weight);
302
303 DphiZH[i]=fabs(Hphi[i]-Zphi);
304 if(DphiZH[i]>M_PI)DphiZH[i] = (2.*M_PI-DphiZH[i]);
305
306 hptz[i]->Fill(Zpt,Weight);
307 if(Zpt>150.){
308
309 ZptSelection[i]++;
310
311 hpth[i]->Fill(Hpt[i],Weight);
312 if(Hpt[i]>150.){
313
314 HptSelection[i]++;
315
316 hdphiZH[i]->Fill(DphiZH[i],Weight);
317 if(DphiZH[i]>2.5){
318
319 DphiSelection[i]++;
320
321 hhinvmass2[i]->Fill(Hinvmass[i],Weight);
322 if(Hinvmass[i]<140. && Hinvmass[i]>100.){
323
324 HMassWindow[i]++;
325
326 }// H mass window
327 }//dphi cut
328 }//h pt cut
329 }//z pt cut
330 }//bjet et > 30 GeV
331 }//bjet pair
332 }//loop over b-tagging algo and tk geometries
333 }//Z mass window
334 }//muon pt>20 GeV
335 }//muon.size==2
336 }//for loop on events
337 delete rootTree;
338 delete ufs;
339 delete file;
340 }//for loop on files
341
342 string algoname[6];
343 algoname[0]="Loose b-tagging, Tk. Std. Geom.";
344 algoname[1]="Loose b-tagging, Tk. Phase 1";
345 algoname[2]="Medium b-tagging, Tk. Std. Geom.";
346 algoname[3]="Medium b-tagging, Tk. Phase 1";
347 algoname[4]="Tight b-tagging, Tk. Std. Geom.";
348 algoname[5]="Tight b-tagging, Tk. Phase 1";
349
350 for(int i=0;i<6;i++){
351 cout<<" ====================================== " <<endl;
352 cout<<" Algo : "<<algoname[i]<<endl;
353 cout<<" Analyzed Events : "<<AnalyzedEvents<<endl;
354 cout<<" ====================================== " <<endl;
355 cout<<" Lepton pair preselection : "<<LeptonPairPreSelection<<endl;
356 cout<<" Lepton pair selection : "<<LeptonPairSelection<<endl;
357 cout<<" Z Invariant Mass Window : "<<ZMassWindow<<endl;
358 cout<<" B-jet pair preselection : "<<BJetPairPreSelection[i]<<endl;
359 cout<<" B-jet pair selection : "<<BJetPairSelection[i]<<endl;
360 cout<<" Z pt : "<<ZptSelection[i]<<endl;
361 cout<<" H pt : "<<HptSelection[i]<<endl;
362 cout<<" Dphi(H,Z) : "<<DphiSelection[i]<<endl;
363 cout<<" H Invariant Mass Window : "<<HMassWindow[i]<<endl;
364 cout<<" Total efficiency (%) : "<<100*double(HMassWindow[i])/AnalyzedEvents<<endl;
365 cout<<" Number of Events : "<<double(HMassWindow[i])/AnalyzedEvents*CrossSection*Lumi<<endl;
366 cout<<" ====================================== " <<endl;
367 }
368 fout.cd();
369 fout.Write();
370 fout.Close();
371
372 return 0;
373 }