Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 256

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

include statements have been cleaned
File size: 14.9 KB
Line 
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"
16#include "TFile.h"
17
18#include "ExRootTreeReader.h"
19#include "ExRootTreeWriter.h"
20#include "ExRootTreeBranch.h"
21#include "TreeClasses.h"
22
23#include "DataConverter.h"
24#include "HEPEVTConverter.h"
25#include "LHEFConverter.h"
26#include "STDHEPConverter.h"
27
28#include "SmearUtil.h"
29#include "JetsUtil.h"
30#include "BFieldProp.h"
31
32//#include "PseudoJet.hh"
33//#include "ClusterSequence.hh"
34
35#include<vector>
36#include<iostream>
37
38using namespace std;
39
40//------------------------------------------------------------------------------
41
42// //********************************** PYTHIA INFORMATION*********************************
43
44TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
45 {
46 TIter it((TCollection*)GEN);
47 it.Reset();
48 TRootGenParticle *gen1;
49 TSimpleArray<TRootGenParticle> array,array2;
50
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)
82{
83 JETSm.SetPxPyPzE(0,0,0,0);
84 float deltaRtest=5000;
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 }
97}
98
99int main(int argc, char *argv[])
100{
101 int appargc = 2;
102 char *appName = "Smearing";
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 }
122
123 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
124 outputFile->Close();
125
126 string line;
127 ifstream infile(inputFileList.c_str());
128 infile >> line; // the first line determines the type of input files
129
130 DataConverter *converter=0;
131
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;};
157
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());
171
172 TRootGenParticle *particle;
173
174 RESOLELEC * elementElec;
175 RESOLMUON *elementMuon;
176 RESOLJET *elementJet;
177 TAUHAD *elementTaujet;
178 ETMIS *elementEtmis;
179
180int numTau=0;
181int numTauRec=0;
182
183 //read the datacard input file
184 string DetDatacard("data/DataCardDet.dat");
185 if(argc==4) DetDatacard =argv[3];
186 RESOLution *DET = new RESOLution();
187 DET->ReadDataCard(DetDatacard);
188
189 //Jet information
190 JetsUtil *JETRUN = new JetsUtil(DetDatacard);
191
192 //Propagation of tracks in the B field
193 //TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
194
195 TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
196 TLorentzVector recoMomentumCalo(0,0,0,0);
197 TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
198 LorentzVector jetMomentum;
199 vector<TLorentzVector> TrackCentral;
200
201 vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
202 vector<fastjet::PseudoJet> sorted_jetsGEN;
203
204 vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
205 vector<fastjet::PseudoJet> sorted_jetsReco;
206
207 vector<PhysicsTower> towers;
208
209 float iPhi=0,iEta=0;
210
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 {
216 TLorentzVector PTmisReco(0,0,0,0);
217 TLorentzVector PTmisGEN(0,0,0,0);
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;
226
227 input_particlesReco.clear();
228 input_particlesGEN.clear();
229 towers.clear();
230
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);
235
236 int pid = abs(particle->PID);
237 float eta = fabs(particle->Eta);
238
239 //input generator level particle for jet algorithm
240 if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
241 {
242 input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
243 }
244
245 //Calculate ETMIS from generated particles
246 if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
247
248 if( (particle->Status == 1) &&
249 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
250 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
251 )
252 {
253 recoMomentum = genMomentum;
254 //use of the magnetic field propagation
255 //TRACP->Propagation(particle,recoMomentum);
256 // cout<<"eta "<<eta<<endl;
257 eta=fabs(recoMomentum.Eta());
258 // cout<<"eta apres"<<eta<<endl;
259
260 switch(pid) {
261
262 case pE: // all electrons with eta < DET->MAX_CALO_FWD
263 DET->SmearElectron(recoMomentum);
264 if(recoMomentum.E() !=0 && eta < DET->CEN_max_tracker){
265 elementElec=(RESOLELEC*) branchelec->NewEntry();
266 elementElec->E = genMomentum.E();
267 elementElec->SmearedE = recoMomentum.E();}
268 break; // case pE
269 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
270 DET->SmearElectron(recoMomentum);
271 break; // case pGAMMA
272 case pMU: // all muons with eta < DET->MAX_MU
273 DET->SmearMu(recoMomentum);
274 if(recoMomentum.E() !=0){
275 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
276 elementMuon->OverPT = (1/genMomentum.Pt());
277 elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
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
281 DET->SmearHadron(recoMomentum, 0.7);
282 break; // case hadron
283 default: // all other final particles with eta < DET->MAX_CALO_FWD
284 DET->SmearHadron(recoMomentum, 1.0);
285 break;
286 } // switch (pid)
287
288 //information to reconstruct jets from reconstructed towers
289 int charge=Charge(pid);
290 if(recoMomentum.E() !=0 && pid != pMU) {
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
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 }
304
305 // all final charged particles
306 if(
307 (recoMomentum.E()!=0) &&
308 (fabs(recoMomentum.Eta()) < DET->CEN_max_tracker) &&
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
311 ((rand()%100) < DET->TRACK_eff) &&
312 (charge!=0)
313 )
314 {TrackCentral.push_back(recoMomentum);}
315
316 } // switch
317 } // while
318
319 //compute missing transverse energy from calo towers
320 TLorentzVector Att(0.,0.,0.,0.);
321 float ScalarEt=0;
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 }
334
335 elementEtmis= (ETMIS*) branchetmis->NewEntry();
336 elementEtmis->Et = (PTmisGEN).Pt();
337 elementEtmis->Ex = (-PTmisGEN).Px();
338 elementEtmis->SEt = ScalarEt;
339
340 elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
341 elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
342
343 //*****************************
344
345 sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
346 sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
347
348 TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
349
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)
356 {
357 elementJet= (RESOLJET*) branchjet->NewEntry();
358 elementJet->PT = JETgen.Et();
359 elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
360 }
361 }
362 numTau = numTau+TausHadr.GetEntries();
363 for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
364 TLorentzVector JETT(0,0,0,0);
365 JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
366 if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
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());
374 elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi());
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
378 }
379 }
380 }
381
382
383 } // for itJet : loop on all jets
384
385 treeWriter->Fill();
386 } // Loop over all events
387 treeWriter->Write();
388float frac = numTauRec/numTau;
389cout<<numTauRec<<endl;
390cout<<numTau<<endl;
391
392 cout << "** Exiting..." << endl;
393 cout<<frac<<endl;
394
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.