Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 449

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

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