Fork me on GitHub

source: svn/trunk/Resolutions_ATLAS.cpp@ 513

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

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