Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 468

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

cleaning, comments, progressbar, cout
File size: 18.2 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/// \file Resolution.cpp
35/// \brief Resolution for CMS
36
37#include "TChain.h"
38#include "TApplication.h"
39#include "TFile.h"
40#include "TStopwatch.h"
41
42#include "ExRootTreeReader.h"
43#include "ExRootTreeWriter.h"
44#include "ExRootTreeBranch.h"
45#include "ExRootProgressBar.h"
46#include "TreeClasses.h"
47
48#include "SmearUtil.h"
49#include "JetsUtil.h"
50#include "BFieldProp.h"
51
52#include<vector>
53#include<iostream>
54
55using namespace std;
56
57//------------------------------------------------------------------------------
58
59// //********************************** PYTHIA INFORMATION*********************************
60
61TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
62 {
63 TIter it((TCollection*)GEN);
64 it.Reset();
65 TRootGenParticle *gen1;
66 TSimpleArray<TRootGenParticle> array,array2;
67
68 while((gen1 = (TRootGenParticle*) it.Next()))
69 {
70 array.Add(gen1);
71 }
72 it.Reset();
73 bool tauhad;
74 while((gen1 = (TRootGenParticle*) it.Next()))
75 {
76 tauhad=false;
77 if(abs(gen1->PID)==15)
78 {
79 int d1=gen1->D1;
80 int d2=gen1->D2;
81 if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
82 {
83 tauhad=true;
84 for(int d=d1; d < d2+1; d++)
85 {
86 if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
87 }
88 }
89 }
90 if(tauhad)array2.Add(gen1);
91 }
92 return array2;
93 }
94
95double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
96 double Energie=0;
97 for(unsigned int i=0; i < towers.size(); i++) {
98 if(towers[i].Pt() < JET_seed) continue;
99 if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
100 Energie += towers[i].E();
101 }
102 }
103 return Energie;
104}
105
106
107void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
108{
109 JETSm.SetPxPyPzE(0,0,0,0);
110 float deltaRtest=5000;
111 TIter itJet((TCollection*)branchJet);
112 TRootJet *jet;
113 itJet.Reset();
114 while( (jet = (TRootJet*) itJet.Next()) )
115 {
116 TLorentzVector Att;
117 Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E);
118 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
119 {
120 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
121 if(deltaRtest < 0.25)
122 {
123 JETSm = Att;
124 }
125 }
126 }
127}
128
129void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
130{
131 ELECSm.SetPxPyPzE(0,0,0,0);
132 float deltaRtest=5000;
133 TIter itElec((TCollection*)branchElec);
134 TRootElectron *elec;
135 itElec.Reset();
136 while( (elec = (TRootElectron*) itElec.Next()) )
137 {
138 TLorentzVector Att;
139 Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E);
140 if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
141 {
142 deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
143 if(deltaRtest < 0.025)
144 {
145 ELECSm = Att;
146 }
147 }
148 }
149}
150
151void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
152{
153 MUONSm.SetPxPyPzE(0,0,0,0);
154 float deltaRtest=5000;
155 TIter itMuon((TCollection*)branchMuon);
156 TRootMuon *muon;
157 itMuon.Reset();
158 while( (muon = (TRootMuon*) itMuon.Next()) )
159 {
160 TLorentzVector Att;
161 Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
162 if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
163 {
164 deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
165 if(deltaRtest < 0.025)
166 {
167 MUONSm = Att;
168 }
169 }
170 }
171}
172
173unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
174 unsigned int numtrack=0;
175 TIter itTrack((TCollection*)branchTracks);
176 TRootTracks *track;
177 itTrack.Reset();
178 while( (track = (TRootTracks*) itTrack.Next()) )
179 {
180 if((track->PT < pt_track )||
181 (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
182 )continue;
183 numtrack++;
184 }
185 return numtrack;
186}
187
188
189
190int main(int argc, char *argv[])
191{
192 int appargc = 2;
193 char *appName= new char[20];
194 char *appOpt= new char[20];
195 sprintf(appName,"Resolution_ATLAS");
196 sprintf(appOpt,"-b");
197 char *appargv[] = {appName,appOpt};
198 TApplication app(appName, &appargc, appargv);
199 delete [] appName;
200 delete [] appOpt;
201
202 if(argc != 3) {
203 cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
204 cout << " input_file - input file in Delphes-root format," << endl;
205 cout << " output_file - output file." << endl;
206 exit(1);
207 }
208
209
210 // 1. ********** initialisation ***********
211
212 srand (time (NULL));
213 TStopwatch globalwatch;
214 globalwatch.Start();
215
216 string inputfilename(argv[1]), outputfilename(argv[2]);
217
218 // 1.1 checks the input file
219 if(inputfilename.find(".root") > inputfilename.length() ) {
220 cout << "input_file should be a .root file from Delphes!\n";
221 return -1;
222 }
223 TFile f(inputfilename.c_str());
224 if (!f.FindKey("GEN") || !f.FindKey("Analysis") ) {
225 cout << "input_file should be a .root file from Delphes!\n";
226 cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
227 return -1;
228 }
229 f.Close();
230
231 // 1.2 checks the output file
232 if(outputfilename.find(".root") > outputfilename.length() ) {
233 cout << "output_file should be a .root file!\n";
234 return -1;
235 }
236 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
237 outputFile->Close();
238
239 // 1.3 Reads the trees in input file
240 TChain chainGEN("GEN");
241 chainGEN.Add(inputfilename.c_str());
242 ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
243 TChain chain("Analysis");
244 chain.Add(inputfilename.c_str());
245 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
246 const TClonesArray *branchJet = treeReader->UseBranch("Jet");
247 const TClonesArray *branchElec = treeReader->UseBranch("Electron");
248 const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
249 const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
250 const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
251 const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
252 TIter itGen((TCollection*)branchGen);
253
254
255 // 1.4 Prepares the output root file
256 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
257 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
258 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
259 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
260 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
261 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
262
263 // 1.5 other initialisations
264 TRootGenParticle *particle;
265 RESOLELEC * elementElec;
266 RESOLMUON *elementMuon;
267 RESOLJET *elementJet;
268 TAUHAD *elementTaujet;
269 ETMIS *elementEtmis;
270
271 int numTau=0, numTauRec=0;
272
273 RESOLution *DET = new RESOLution();
274 /*
275 string detectorcard = "data/DetectorCard_CMS.dat";
276 const float dR_jetpairing = 0.25;
277 const float jet_pt_cut = 1;
278 */
279 string detectorcard = "data/DetectorCard_CMS.dat";
280 DET->ReadDataCard(detectorcard);
281
282 JetsUtil *JETRUN = new JetsUtil(detectorcard);
283
284 TLorentzVector genMomentum(0,0,0,0); //generator level information
285 TLorentzVector recoMomentum(0,0,0,0);//reconstructed level information
286 LorentzVector jetMomentum;
287
288 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
289 vector<fastjet::PseudoJet> sorted_jetsGEN;
290 vector<int> NTrackJet;
291 vector<TLorentzVector> towers;
292
293 // 2. Loop over all events
294 Long64_t entry, allEntries = treeReader->GetEntries();
295 cout << endl << endl;
296 cout <<"*********************************************************************"<< endl;
297 cout <<"** **"<< endl;
298 cout <<"** **"<< endl;
299 cout <<"** **"<< endl;
300 cout <<"** ####### Start resolution processing ######## **"<< endl;
301 cout << left << setw(52) <<"** Total number of events to run: "<<""
302 << left << setw(15) << allEntries <<""
303 << right << setw(2) <<"**"<<endl;
304
305 ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
306
307 for(entry = 0; entry < allEntries; ++entry)
308 {
309 Progress->Update(entry);
310 TLorentzVector PTmisReco(0,0,0,0);
311 TLorentzVector PTmisGEN(0,0,0,0);
312 treeReader->ReadEntry(entry);
313 treeReaderGEN->ReadEntry(entry);
314 treeWriter->Clear();
315
316 TSimpleArray<TRootGenParticle> bGen;
317 itGen.Reset();
318 TSimpleArray<TRootGenParticle> NFCentralQ;
319
320 input_particlesGEN.clear();
321 towers.clear();
322
323 // Loop over all particles in event
324 while( (particle = (TRootGenParticle*) itGen.Next()) )
325 {
326 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
327
328 int pid = abs(particle->PID);
329 float eta = fabs(particle->Eta);
330
331 //input generator level particle for jet algorithm
332 if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
333 {
334 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
335 }
336
337 //Calculate ETMIS from generated particles
338 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
339
340 //Electrons and muons
341 if( (particle->Status == 1) &&
342 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
343 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
344 )
345 {
346 eta=fabs(genMomentum.Eta());
347
348 switch(pid) {
349
350 case pE: // all electrons with eta < DET->MAX_CALO_FWD
351 PairingElec(recoMomentum,genMomentum,branchElec);
352 if(recoMomentum.E()!=0){
353 elementElec=(RESOLELEC*) branchelec->NewEntry();
354 elementElec->E = genMomentum.E();
355 elementElec->SmearedE = recoMomentum.E();}
356 break; // case pE
357 case pMU: // all muons with eta < DET->MAX_MU
358 PairingMuon(recoMomentum,genMomentum,branchMuon);
359 if(recoMomentum.E()!=0){
360 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
361 elementMuon->OverPT = 1./genMomentum.Pt();
362 elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
363 break; // case pMU
364 default:
365 break;
366 } // switch (pid)
367 }
368
369 } // while
370
371 //compute missing transverse energy from calo towers
372 TIter itCalo((TCollection*)branchTowers);
373 TRootCalo *calo;
374 itCalo.Reset();
375 TLorentzVector Att(0.,0.,0.,0.);
376 float ScalarEt=0;
377 while( (calo = (TRootCalo*) itCalo.Next()) )
378 {
379 if(calo->E !=0){
380 Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E);
381 towers.push_back(Att);
382 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
383 {
384 ScalarEt = ScalarEt + calo->getET();
385 PTmisReco = PTmisReco + Att;
386 }
387 }
388 }
389 elementEtmis= (ETMIS*) branchetmis->NewEntry();
390 elementEtmis->Et = (PTmisGEN).Pt();
391 elementEtmis->Ex = (-PTmisGEN).Px();
392 elementEtmis->SEt = ScalarEt;
393 elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
394 elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
395
396 //*****************************
397 sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN);
398
399 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
400 TLorentzVector JETreco(0,0,0,0);
401 for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
402 TLorentzVector JETgen(0,0,0,0);
403 JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
404 PairingJet(JETreco,JETgen,branchJet);
405 if(JETreco.Pt()>1)
406 {
407 elementJet= (RESOLJET*) branchjet->NewEntry();
408 elementJet->PT = JETgen.Et();
409 elementJet->SmearedPT = JETreco.Et()/JETgen.Et(); //// la difference pourrait être ici
410 }
411 }
412 numTau += TausHadr.GetEntries();
413
414 TIter itJet((TCollection*)branchJet);
415 TRootJet *jet;
416 itJet.Reset();
417 while( (jet = (TRootJet*) itJet.Next()) )
418 {
419 TLorentzVector JETT(0,0,0,0);
420 JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
421 if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
422 {
423 for(Int_t i=0; i<TausHadr.GetEntries();i++)
424 {
425 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
426 {
427 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
428 elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
429 elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
430 if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95
431 && (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
432 }
433 }
434 }
435
436
437 } // for itJet : loop on all jets
438
439 treeWriter->Fill();
440 } // Loop over all events
441 treeWriter->Write();
442 globalwatch.Stop();
443
444
445 // Screen output
446 cout <<"** **"<< endl;
447 cout <<"** **"<< endl;
448 cout <<"** ################## Time report ################# **"<< endl;
449 cout << left << setw(32) <<"** Time report for "<<""
450 << left << setw(15) << allEntries <<""
451 << right << setw(22) <<"events **"<<endl;
452 cout <<"** **"<< endl;
453 cout << left << setw(10) <<"**"<<""
454 << left << setw(15) <<"Time (s):"<<""
455 << right << setw(15) <<"CPU"<<""
456 << right << setw(15) <<"Real"<<""
457 << right << setw(14) <<"**"<<endl;
458 cout << left << setw(10) <<"**"<<""
459 << left << setw(15) <<" + Global:"<<""
460 << right << setw(15) <<globalwatch.CpuTime()<<""
461 << right << setw(15) <<globalwatch.RealTime()<<""
462 << right << setw(14) <<"**"<<endl;
463
464
465 string tempstring = detectorcard + " has been used.";
466 cout <<"** ################## Configuration ############### **"<< endl;
467 cout << left << setw(16) << "** " << ""
468 << left << setw(51) << tempstring << "**" << endl;
469
470 tempstring = outputfilename + " has been created.";
471 cout << left << setw(16) << "** " << ""
472 << left << setw(51) << tempstring << "**" << endl;
473
474 cout << left << setw(16) << "** " << ""
475 << left << setw(51) << get_time_date() << "**" << endl;
476
477
478 cout <<"** **"<< endl;
479 cout <<"** Exiting ... **"<< endl;
480 cout <<"** **"<< endl;
481 cout <<"** **"<< endl;
482 cout <<"*********************************************************************"<< endl;
483
484
485
486 delete treeWriter;
487 delete treeReader;
488 delete DET;
489}
490
Note: See TracBrowser for help on using the repository browser.