Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 250

Last change on this file since 250 was 227, checked in by Xavier Rouby, 16 years ago

include statements have been cleaned
File size: 14.9 KB
RevLine 
[9]1/*
2 ---- FastSim ----
3 A Fast Simulator for general purpose LHC detector
4 S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5
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
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
[227]32//#include "PseudoJet.hh"
33//#include "ClusterSequence.hh"
[88]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
65 if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
66 {
67 tauhad=true;
68 for(int d=d1; d < d2+1; d++)
69 {
70 if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
71 }
72 }
73 }
74 if(tauhad)array2.Add(gen1);
75 }
76 return array2;
77 }
78
79
80
81void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, vector<fastjet::PseudoJet> sorted_jetsS)
[9]82{
83 JETSm.SetPxPyPzE(0,0,0,0);
84 float deltaRtest=5000;
[19]85 for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
86 TLorentzVector Att;
87 Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
88 if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
89 {
90 deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
91 if(deltaRtest < 0.25)
92 {
93 JETSm = Att;
94 }
95 }
96 }
[9]97}
98
99int main(int argc, char *argv[])
100{
101 int appargc = 2;
[26]102 char *appName = "Smearing";
[9]103 char *appargv[] = {appName, "-b"};
104 TApplication app(appName, &appargc, appargv);
105
106 if(argc != 4 && argc != 3) {
107 cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
108 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
109 cout << " output_file - output file." << endl;
110 cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
111 exit(1);
112 }
113
114 srand (time (NULL)); /* Initialisation du générateur */
115
116 //read the input TROOT file
117 string inputFileList(argv[1]), outputfilename(argv[2]);
118 if(outputfilename.find(".root") > outputfilename.length() ) {
119 cout << "output_file should be a .root file!\n";
120 return -1;
121 }
[88]122
[9]123 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
124 outputFile->Close();
[26]125
[9]126 string line;
127 ifstream infile(inputFileList.c_str());
128 infile >> line; // the first line determines the type of input files
[26]129
[9]130 DataConverter *converter=0;
[26]131
[9]132 if(strstr(line.c_str(),".hep"))
133 {
134 cout<<"*************************************************************************"<<endl;
135 cout<<"************ StdHEP file format detected **************"<<endl;
136 cout<<"************ Starting convertion to TRoot format **************"<<endl;
137 cout<<"*************************************************************************"<<endl;
138 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
139 }
140 else if(strstr(line.c_str(),".lhe"))
141 {
142 cout<<"*************************************************************************"<<endl;
143 cout<<"************ LHEF file format detected **************"<<endl;
144 cout<<"************ Starting convertion to TRoot format **************"<<endl;
145 cout<<"*************************************************************************"<<endl;
146 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
147 }
148 else if(strstr(line.c_str(),".root"))
149 {
150 cout<<"*************************************************************************"<<endl;
151 cout<<"************ h2root file format detected **************"<<endl;
152 cout<<"************ Starting convertion to TRoot format **************"<<endl;
153 cout<<"*************************************************************************"<<endl;
154 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
155 }
156 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
[26]157
[9]158 TChain chain("GEN");
159 chain.Add(outputfilename.c_str());
160 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
161 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
162 TIter itGen((TCollection*)branchGen);
163
164 //write the output root file
165 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
166 ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
167 ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
168 ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
169 ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
170 ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
[26]171
[9]172 TRootGenParticle *particle;
[26]173
[39]174 RESOLELEC * elementElec;
[9]175 RESOLMUON *elementMuon;
176 RESOLJET *elementJet;
177 TAUHAD *elementTaujet;
178 ETMIS *elementEtmis;
[26]179
[152]180int numTau=0;
181int numTauRec=0;
[26]182
[9]183 //read the datacard input file
[89]184 string DetDatacard("data/DataCardDet.dat");
[9]185 if(argc==4) DetDatacard =argv[3];
186 RESOLution *DET = new RESOLution();
187 DET->ReadDataCard(DetDatacard);
[88]188
189 //Jet information
[100]190 JetsUtil *JETRUN = new JetsUtil(DetDatacard);
[88]191
192 //Propagation of tracks in the B field
[202]193 //TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
[88]194
195 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
[152]196 TLorentzVector recoMomentumCalo(0,0,0,0);
[88]197 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
[9]198 LorentzVector jetMomentum;
[19]199 vector<TLorentzVector> TrackCentral;
[26]200
[88]201 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
202 vector<fastjet::PseudoJet> sorted_jetsGEN;
[26]203
[88]204 vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
205 vector<fastjet::PseudoJet> sorted_jetsReco;
[39]206
[23]207 vector<PhysicsTower> towers;
[152]208
209 float iPhi=0,iEta=0;
210
[9]211 // Loop over all events
212 Long64_t entry, allEntries = treeReader->GetEntries();
213 cout << "** Chain contains " << allEntries << " events" << endl;
214 for(entry = 0; entry < allEntries; ++entry)
215 {
[88]216 TLorentzVector PTmisReco(0,0,0,0);
217 TLorentzVector PTmisGEN(0,0,0,0);
[9]218 treeReader->ReadEntry(entry);
219 treeWriter->Clear();
220 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
221
222 TSimpleArray<TRootGenParticle> bGen;
223 itGen.Reset();
224 TrackCentral.clear();
225 TSimpleArray<TRootGenParticle> NFCentralQ;
[89]226
[88]227 input_particlesReco.clear();
228 input_particlesGEN.clear();
[23]229 towers.clear();
[26]230
[9]231 // Loop over all particles in event
232 while( (particle = (TRootGenParticle*) itGen.Next()) )
233 {
234 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
[26]235
[9]236 int pid = abs(particle->PID);
237 float eta = fabs(particle->Eta);
[89]238
[88]239 //input generator level particle for jet algorithm
[94]240 if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
[26]241 {
[88]242 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
[26]243 }
[89]244
[88]245 //Calculate ETMIS from generated particles
246 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
[89]247
[88]248 if( (particle->Status == 1) &&
249 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
[94]250 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
[88]251 )
252 {
253 recoMomentum = genMomentum;
254 //use of the magnetic field propagation
[152]255 //TRACP->Propagation(particle,recoMomentum);
[89]256 // cout<<"eta "<<eta<<endl;
257 eta=fabs(recoMomentum.Eta());
258 // cout<<"eta apres"<<eta<<endl;
259
[26]260 switch(pid) {
[89]261
[26]262 case pE: // all electrons with eta < DET->MAX_CALO_FWD
[88]263 DET->SmearElectron(recoMomentum);
[94]264 if(recoMomentum.E() !=0 && eta < DET->CEN_max_tracker){
[89]265 elementElec=(RESOLELEC*) branchelec->NewEntry();
266 elementElec->E = genMomentum.E();
267 elementElec->SmearedE = recoMomentum.E();}
[26]268 break; // case pE
269 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
[88]270 DET->SmearElectron(recoMomentum);
[26]271 break; // case pGAMMA
272 case pMU: // all muons with eta < DET->MAX_MU
[88]273 DET->SmearMu(recoMomentum);
[89]274 if(recoMomentum.E() !=0){
[39]275 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
[88]276 elementMuon->OverPT = (1/genMomentum.Pt());
[89]277 elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
[26]278 break; // case pMU
279 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
280 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
[88]281 DET->SmearHadron(recoMomentum, 0.7);
[26]282 break; // case hadron
283 default: // all other final particles with eta < DET->MAX_CALO_FWD
[88]284 DET->SmearHadron(recoMomentum, 1.0);
[26]285 break;
286 } // switch (pid)
[89]287
288 //information to reconstruct jets from reconstructed towers
[152]289 int charge=Charge(pid);
290 if(recoMomentum.E() !=0 && pid != pMU) {
[182]291 // in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
292 if ( !DET->FLAG_bfield && charge!=0 && genMomentum.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
293 else { // particles reach calos
[152]294 DET->BinEtaPhi(recoMomentum.Phi(), recoMomentum.Eta(), iPhi, iEta);
295 if(iEta != -100 && iPhi != -100)
296 {
297 recoMomentumCalo.SetPtEtaPhiE(recoMomentum.Pt(),iEta,iPhi,recoMomentum.E());
298
299 PhysicsTower Tower(LorentzVector(recoMomentumCalo.Px(),recoMomentumCalo.Py(),recoMomentumCalo.Pz(), recoMomentumCalo.E()));
300 towers.push_back(Tower);
301 }
302 }
303 }
[182]304
[26]305 // all final charged particles
[89]306 if(
[152]307 (recoMomentum.E()!=0) &&
308 (fabs(recoMomentum.Eta()) < DET->CEN_max_tracker) &&
[182]309 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentum.Pt() > DET->TRACK_ptmin )) &&
310 // if bfield not simulated, pt should be high enough to be taken into account
[94]311 ((rand()%100) < DET->TRACK_eff) &&
[89]312 (charge!=0)
313 )
[152]314 {TrackCentral.push_back(recoMomentum);}
[89]315
[9]316 } // switch
317 } // while
[89]318
319 //compute missing transverse energy from calo towers
320 TLorentzVector Att(0.,0.,0.,0.);
321 float ScalarEt=0;
[152]322 for(unsigned int i=0; i < towers.size(); i++)
323 {
324 Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
325 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
326 {
327 ScalarEt = ScalarEt + Att.Et();
328 PTmisReco = PTmisReco + Att;
329 // create a fastjet::PseudoJet with these components and put it onto
330 // back of the input_particles vector
331 input_particlesReco.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
332 }
333 }
[89]334
335 elementEtmis= (ETMIS*) branchetmis->NewEntry();
336 elementEtmis->Et = (PTmisGEN).Pt();
337 elementEtmis->Ex = (-PTmisGEN).Px();
338 elementEtmis->SEt = ScalarEt;
[88]339
[89]340 elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
341 elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
[9]342
343 //*****************************
[26]344
[88]345 sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
346 sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
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());
354 PairingJet(JETreco,JETgen,sorted_jetsReco);
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();
[88]363 for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
[39]364 TLorentzVector JETT(0,0,0,0);
[88]365 JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
[94]366 if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
[39]367 {
368 for(Int_t i=0; i<TausHadr.GetEntries();i++)
369 {
370 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
371 {
372 elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
373 elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
[94]374 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi());
[152]375 if( (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E()) > 0.95
376 && (DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi()))==1)numTauRec++;
377
[39]378 }
379 }
380 }
381
[26]382
383 } // for itJet : loop on all jets
384
[9]385 treeWriter->Fill();
386 } // Loop over all events
387 treeWriter->Write();
[152]388float frac = numTauRec/numTau;
389cout<<numTauRec<<endl;
390cout<<numTau<<endl;
391
[9]392 cout << "** Exiting..." << endl;
[152]393 cout<<frac<<endl;
394
[9]395
396 delete treeWriter;
397 delete treeReader;
398 delete DET;
399 if(converter) delete converter;
400}
401
Note: See TracBrowser for help on using the repository browser.