Fork me on GitHub

source: svn/trunk/Resolutions.cpp@ 206

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

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