Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 326

Last change on this file since 326 was 310, checked in by severine ovyn, 16 years ago

add Ehad/Eem

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