Fork me on GitHub

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

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

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