Fork me on GitHub

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

Last change on this file since 657 was 443, checked in by Xavier Rouby, 15 years ago

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