Fork me on GitHub

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

Last change on this file since 84 was 84, checked in by uid677, 16 years ago

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