Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 973

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

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