Fork me on GitHub

source: svn/trunk/Examples/src/Analysis_Ex.cc@ 245

Last change on this file since 245 was 227, checked in by Xavier Rouby, 16 years ago

include statements have been cleaned
File size: 13.8 KB
Line 
1#include "Examples/interface/Analysis_Ex.h"
2#include <iostream>
3#include <sstream>
4#include <fstream>
5#include <iomanip>
6
7using namespace std;
8
9//******************************Debut de l'analyse****************************************
10//*****************************************************************************************
11
12Analysis_Ex::Analysis_Ex(string CardWithCuts,string LogName)
13{
14 string temp_string;
15 istringstream curstring;
16
17 ifstream fichier_a_lire(CardWithCuts.c_str());
18 if(!fichier_a_lire.good()) {
19 cout << "DataCardname " << CardWithCuts << " not found, use default values" << endl;
20 return;
21 }
22
23 while (getline(fichier_a_lire,temp_string)) {
24 curstring.clear(); // needed when using several times istringstream::str(string)
25 curstring.str(temp_string);
26 string varname;
27 float value;
28
29 if(strstr(temp_string.c_str(),"#")) { }//remove comments
30 else if(strstr(temp_string.c_str(),"PT_ELEC")){curstring >> varname >> value; PT_ELEC = value;}
31 else if(strstr(temp_string.c_str(),"PT_MUON")){curstring >> varname >> value; PT_MUON = value;}
32 else if(strstr(temp_string.c_str(),"INV_MASS_LL")){curstring >> varname >> value; INV_MASS_LL = value;}
33 }
34
35 ofstream f_out(LogName.c_str(),ofstream::app);
36
37 f_out<<"*******************************************************************"<<endl;
38 f_out << left << setw(30) <<"Cut values used in the analysis: "<<""
39 << right << setw(37) <<"-------------------------------------"<<"\n";
40 f_out << left <<setw(50) << "Invariant mass of the leptons: "<<""
41 << right <<setw(17) << INV_MASS_LL <<"\n";
42 f_out<<"*******************************************************************"<<endl;
43
44}
45
46Analysis_Ex::~Analysis_Ex()
47{
48}
49
50void Analysis_Ex::Run(ExRootTreeReader *treeReaderGen, ExRootTreeReader *treeReaderRec, ExRootTreeReader *treeReaderTrig, ExRootTreeWriter *treeWriter)
51{
52 total=0;//initialisation of total number of events
53 cut_trig=0;
54 cut_1=0;//initialisation of counter for cut 1
55 cut_2=0;
56 //access the branches************************
57 //to get the generator level information
58 const TClonesArray *GEN = treeReaderGen->UseBranch("Particle");
59
60 //to get the reconstructed level information
61 const TClonesArray *JET = treeReaderRec->UseBranch("Jet");
62 const TClonesArray *TAUJET = treeReaderRec->UseBranch("TauJet");
63 const TClonesArray *PHOTO = treeReaderRec->UseBranch("Photon");
64 const TClonesArray *ELEC = treeReaderRec->UseBranch("Electron");
65 const TClonesArray *MUON = treeReaderRec->UseBranch("Muon");
66 const TClonesArray *TRACKS = treeReaderRec->UseBranch("Tracks");
67 const TClonesArray *CALO = treeReaderRec->UseBranch("CaloTower");
68
69 //to get the VFD reconstructed level information
70 const TClonesArray *ZDC = treeReaderRec->UseBranch("ZDChits");
71 const TClonesArray *RP220 = treeReaderRec->UseBranch("RP220hits");
72 const TClonesArray *FP420 = treeReaderRec->UseBranch("FP420hits");
73
74 //to get the trigger information
75 const TClonesArray *TRIGGER = treeReaderTrig->UseBranch("TrigResult");
76
77 //Define the branches that will be filled during the analysis
78 ExRootTreeBranch *INVMASS = treeWriter->NewBranch("INVMass", TRootInvm::Class());
79 TRootInvm *inv_mass;
80 //*******************************************
81
82 //run on the events
83 Long64_t entry, allEntries = treeReaderRec->GetEntries();
84 cout << "** Chain contains " << allEntries << " events" << endl;
85 total=allEntries;
86
87 //general information
88 float E,Px,Py,Pz;
89 float PT,Eta,Phi;
90
91 //lepton information
92 bool IsolFlag;
93
94 //bjet information
95 bool Btag;
96
97 //Particle level information
98 int PID, Status, M1,M2,D1,D2;
99 float Charge, T, X, Y, Z, M;
100
101 //VFD information
102 float S,q2,Tx,Ty;
103 int side;
104
105 for(entry = 0; entry < allEntries; ++entry)
106 {
107 treeReaderGen->ReadEntry(entry);//access information of generated information
108 treeReaderRec->ReadEntry(entry);//access information of reconstructed information
109 treeReaderTrig->ReadEntry(entry);//access information of Trigger information
110
111 //*****************************************************
112 //Example how to run on the generator level information
113 //*****************************************************
114 TIter itGen((TCollection*)GEN);
115 TRootGenParticle *gen;
116 itGen.Reset();
117 while( (gen = (TRootGenParticle*) itGen.Next()) )
118 {
119 PID = gen->PID; // particle HEP ID number
120 Status = gen->Status; // particle status
121 M1 = gen->M1; // particle 1st mother
122 M2 = gen->M2; // particle 2nd mother
123 D1 = gen->D1; // particle 1st daughter
124 D2 = gen->D2; // particle 2nd daughter
125 Charge = gen->Charge; // electrical charge
126
127 T = gen->T; // particle vertex position (t component)
128 X = gen->X; // particle vertex position (x component)
129 Y = gen->Y; // particle vertex position (y component)
130 Z = gen->Z; // particle vertex position (z component)
131 M = gen->M; // particle mass
132 }
133
134
135 //***********************************************
136 //Example how to run on the reconstructed objects
137 //***********************************************
138
139 //access the Electron branch
140 TIter itElec((TCollection*)ELEC);
141 TRootElectron *elec;
142 itElec.Reset();
143 while( (elec = (TRootElectron*) itElec.Next()) )
144 {
145 E = elec->E; // particle energy in GeV
146 Px = elec->Px; // particle momentum vector (x component) in GeV
147 Py = elec->Py; // particle momentum vector (y component) in GeV
148 Pz = elec->Pz; // particle momentum vector (z component) in GeV
149
150 PT = elec->PT; // particle transverse momentum in GeV
151 Eta = elec->Eta; // particle pseudorapidity
152 Phi = elec->Phi; // particle azimuthal angle in rad
153 IsolFlag = elec->IsolFlag; // is the particule isolated?
154 }
155 //Running on the muon branch is identical:
156 TIter itMuon((TCollection*)MUON);
157 TRootMuon *muon;
158 itMuon.Reset();
159 while( (muon = (TRootMuon*) itMuon.Next()) ){}
160
161 //access the Photon branch
162 TIter itGam((TCollection*)PHOTO);
163 TRootPhoton *gam;
164 itGam.Reset();
165 while( (gam = (TRootPhoton*) itGam.Next()) )
166 {
167 E = gam->E; // particle energy in GeV
168 Px = gam->Px; // particle momentum vector (x component) in GeV
169 Py = gam->Py; // particle momentum vector (y component) in GeV
170 Pz = gam->Pz; // particle momentum vector (z component) in GeV
171
172 PT = gam->PT; // particle transverse momentum in GeV
173 Eta = gam->Eta; // particle pseudorapidity
174 Phi = gam->Phi; // particle azimuthal angle in rad
175 }
176
177 //access the jet branch
178 TIter itJet((TCollection*)JET);
179 TRootJet *jet;
180 itJet.Reset();
181 while( (jet = (TRootJet*) itJet.Next()) )
182 {
183 E = jet->E; // particle energy in GeV
184 Px = jet->Px; // particle momentum vector (x component) in GeV
185 Py = jet->Py; // particle momentum vector (y component) in GeV
186 Pz = jet->Pz; // particle momentum vector (z component) in GeV
187
188 PT = jet->PT; // particle transverse momentum in GeV
189 Eta = jet->Eta; // particle pseudorapidity
190 Phi = jet->Phi; // particle azimuthal angle in rad
191 Btag = jet->Btag; // is the jet BTagged
192 }
193 //Running on the tau-jet branch is identical:
194 TIter itTaujet((TCollection*)TAUJET);
195 TRootTauJet *taujet;
196 itTaujet.Reset();
197 while( (taujet = (TRootTauJet*) itTaujet.Next()) ){}
198
199 //access the track branch
200 TIter itTrack((TCollection*)TRACKS);
201 TRootTracks *tracks;
202 itTrack.Reset();
203 while( (tracks = (TRootTracks*) itTrack.Next()) )
204 {
205 E = tracks->E; // particle energy in GeV
206 Px = tracks->Px; // particle momentum vector (x component) in GeV
207 Py = tracks->Py; // particle momentum vector (y component) in GeV
208 Pz = tracks->Pz; // particle momentum vector (z component) in GeV
209
210 PT = tracks->PT; // particle transverse momentum in GeV
211 Eta = tracks->Eta; // particle pseudorapidity
212 Phi = tracks->Phi; // particle azimuthal angle in rad
213 }
214
215 //Running on the calo branch is identical:
216 TIter itCalo((TCollection*)CALO);
217 TRootCalo *calo;
218 itCalo.Reset();
219 while( (calo = (TRootCalo*) itCalo.Next()) ){}
220
221 //***************************************************
222 //Example how to run on the VFD reconstructed objects
223 //***************************************************
224
225 //access the ZDC branch
226 TIter itZdc((TCollection*)ZDC);
227 TRootZdcHits *zdc;
228 itZdc.Reset();
229 while( (zdc = (TRootZdcHits*) itZdc.Next()) )
230 {
231 E = zdc->E; // particle energy in GeV
232 Px = zdc->Px; // particle momentum vector (x component) in GeV
233 Py = zdc->Py; // particle momentum vector (y component) in GeV
234 Pz = zdc->Pz; // particle momentum vector (z component) in GeV
235
236 PT = zdc->PT; // particle transverse momentum in GeV
237 Eta = zdc->Eta; // particle pseudorapidity
238 Phi = zdc->Phi; // particle azimuthal angle in rad
239
240 T = zdc->T; // time of flight [s]
241 side = zdc->side; // -1 or +1
242 }
243
244 //access the RP220 branch
245 TIter itRp220((TCollection*)RP220);
246 TRootRomanPotHits *rp220;
247 itRp220.Reset();
248 //TRootRomanPotHits.MakeClass("test_xav");
249 while( (rp220 = (TRootRomanPotHits*) itRp220.Next()) )
250 {
251 T = rp220->T; // time of flight to the detector [s]
252 S = rp220->S; // distance to the IP [m]
253 E = rp220->E; // reconstructed energy [GeV]
254 q2 = rp220->q2; // reconstructed squared momentum transfer [GeV^2]
255
256 X = rp220->X; // horizontal distance to the beam [um]
257 Y = rp220->Y; // vertical distance to the beam [um]
258
259 Tx = rp220->Tx; // angle of the momentum in the horizontal (x,z) plane [urad]
260 Ty = rp220->Ty; // angle of the momentum in the verical (y,z) plane [urad]
261
262 side = rp220->side; // -1 or 1
263 }
264 //running on FP420 branch is identical
265 TIter itFp420((TCollection*)FP420);
266 TRootRomanPotHits *fp420;
267 itFp420.Reset();
268 while( (fp420 = (TRootRomanPotHits*) itFp420.Next()) ){}
269
270 //*********************************************
271 //Example how to run on the trigger information
272 //*********************************************
273
274 TRootTrigger *trig;
275 int NumTrigBit = TRIGGER->GetEntries();
276 //get the global response of the trigger
277
278 bool GlobalResponse=false;
279 if(NumTrigBit!=0)GlobalResponse=true;
280 cout<<"GlobalResponse "<<GlobalResponse<<endl;
281 for(int i=0; i < NumTrigBit-1; i++){
282 trig = (TRootTrigger*)TRIGGER->At(i);
283 cout<<"The event has been accepted by the trigger number: "<<trig->Accepted<<endl;
284 }
285
286 //********************************
287 //Example of a very small analysis
288 //********************************
289
290 TLorentzVector Lept[2];
291
292 if(NumTrigBit==0)continue; //event not accepted by the trigger
293 cut_trig++;//event accepted
294
295 TSimpleArray<TRootElectron> el=SubArrayEl(ELEC,PT_ELEC);//the central isolated electrons, pt > PT_ELEC GeV
296 TSimpleArray<TRootMuon> mu=SubArrayMu(MUON,PT_MUON);//the central isolated electrons, pt > PT_MUON GeV
297
298 Int_t numElec=el.GetEntries();
299
300 if(el.GetEntries()+mu.GetEntries()!=2)continue;//Exactly 2 isolated leptons are needed
301 cut_1++;//event accepted
302 for(Int_t i=0;i < numElec; i++)Lept[i].SetPxPyPzE(el[i]->Px,el[i]->Py,el[i]->Pz,el[i]->E);
303 for(Int_t k = numElec; k < (numElec+mu.GetEntries()); k++)Lept[k].SetPxPyPzE(mu[k-numElec]->Px,mu[k-numElec]->Py,mu[k-numElec]->Pz,mu[k-numElec]->E);
304
305 //Example how to white a branch in the output file
306 inv_mass=(TRootInvm*) INVMASS->NewEntry();
307 inv_mass->M=(Lept[0]+Lept[1]).M();
308
309 if((Lept[0]+Lept[1]).M() > INV_MASS_LL )continue;// the invariant mass should be < INV_MASS_LL
310 cut_2++;//event accepted
311
312 treeWriter->Fill();
313 }
314 treeWriter->Write();
315
316}
317
318void Analysis_Ex::WriteOutput(string LogName)
319{
320 ofstream f_out(LogName.c_str(),ofstream::app);
321
322 f_out<<"*******************************************************************"<<endl;
323 f_out << left << setw(20) << "Numer of Events "<<""
324 << right << setw(15) << total <<"\n";
325 f_out << left << setw(17) << " Accepted by the trigger "<<""
326 << right << setw(20) << cut_trig <<"\n";
327 f_out << left << setw(17) <<" 2 leptons "<< ""
328 << right << setw(20) << cut_1 << ""
329 << right << setw(15) << cut_1/total << "\n";
330 f_out << left << setw(17) <<" Invariant mass "<< ""
331 << right << setw(20) << cut_2 << ""
332 << right << setw(15) << cut_2/total << "\n";
333 f_out<<"*******************************************************************"<<endl;
334 f_out<<" "<<endl;
335
336}
337
338TSimpleArray<TRootElectron> Analysis_Ex::SubArrayEl(const TClonesArray *ELEC,float pt)
339{
340 TIter itElec((TCollection*)ELEC);
341 TRootElectron *elec;
342 itElec.Reset();
343 TSimpleArray<TRootElectron> array;
344 while( (elec = (TRootElectron*) itElec.Next()) )
345 {
346 if(elec->PT<pt)continue;
347 array.Add(elec);
348 }
349 return array;
350}
351
352TSimpleArray<TRootMuon> Analysis_Ex::SubArrayMu(const TClonesArray *MUON,float pt)
353{
354 TIter itMuon((TCollection*)MUON);
355 TRootMuon *muon;
356 itMuon.Reset();
357 TSimpleArray<TRootMuon> array;
358 while( (muon = (TRootMuon*) itMuon.Next()) )
359 {
360 if(muon->PT<pt)continue;
361 array.Add(muon);
362 }
363 return array;
364}
365
Note: See TracBrowser for help on using the repository browser.