Fork me on GitHub

source: svn/trunk/Resolutions_ATLAS.cpp@ 460

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

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