Fork me on GitHub

source: svn/trunk/Delphes.cpp@ 203

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

Performance of timing included

File size: 18.8 KB
Line 
1/*
2 ---- Delphes ----
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 Delphes.cpp
12/// \brief executable for the Delphes
13
14#include "TChain.h"
15#include "TApplication.h"
16#include "TStopwatch.h"
17
18#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
19#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
20#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
21
22#include "interface/DataConverter.h"
23#include "interface/HEPEVTConverter.h"
24#include "interface/LHEFConverter.h"
25#include "interface/STDHEPConverter.h"
26
27#include "interface/SmearUtil.h"
28#include "interface/BFieldProp.h"
29#include "interface/TriggerUtil.h"
30#include "interface/VeryForward.h"
31#include "interface/JetUtils.h"
32#include "interface/FrogUtil.h"
33
34#include <vector>
35#include <iostream>
36
37using namespace std;
38
39//------------------------------------------------------------------------------
40void todo(string filename) {
41 ifstream infile(filename.c_str());
42 cout << "** TODO list ..." << endl;
43 while(infile.good()) {
44 string temp;
45 getline(infile,temp);
46 cout << "*" << temp << endl;
47 }
48 cout << "** done...\n";
49}
50
51//------------------------------------------------------------------------------
52
53int main(int argc, char *argv[])
54{
55 int appargc = 2;
56 char *appName = "Delphes";
57 char *appargv[] = {appName, "-b"};
58 TApplication app(appName, &appargc, appargv);
59
60 if(argc != 4 && argc != 3 && argc != 5) {
61 cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
62 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
63 cout << " output_file - output file." << endl;
64 cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
65 cout << " trigger_card - Datacard containing the trigger algorithms (optional) "<<endl;
66 exit(1);
67 }
68
69// 1. ********** initialisation ***********
70
71 srand (time (NULL)); /* Initialisation du générateur */
72 TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch;
73 globalwatch.Start();
74
75
76 //read the input TROOT file
77 string inputFileList(argv[1]), outputfilename(argv[2]);
78 if(outputfilename.find(".root") > outputfilename.length() ) {
79 cout << "output_file should be a .root file!\n";
80 exit(1);
81 }
82 //create output log-file name
83 string forLog = outputfilename;
84 string LogName = forLog.erase(forLog.find(".root"));
85 LogName = LogName+"_run.log";
86
87 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
88 outputFile->Close();
89
90 string line;
91 ifstream infile(inputFileList.c_str());
92 infile >> line; // the first line determines the type of input files
93
94 //read the datacard input file
95 string DetDatacard("");
96 if(argc>=4) DetDatacard =argv[3];
97
98 //Smearing information
99 RESOLution *DET = new RESOLution();
100 DET->ReadDataCard(DetDatacard);
101 DET->Logfile(LogName);
102
103 //read the trigger input file
104 string TrigDatacard("data/trigger.dat");
105 if(argc==5) TrigDatacard =argv[4];
106
107 //Trigger information
108 TriggerTable *TRIGT = new TriggerTable();
109 TRIGT->TriggerCardReader(TrigDatacard.c_str());
110 TRIGT->PrintTriggerTable(LogName);
111
112 //Propagation of tracks in the B field
113 TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
114
115 //Jet information
116 JetsUtil *JETRUN = new JetsUtil(DetDatacard);
117
118 //VFD information
119 VeryForward * VFD = new VeryForward(DetDatacard);
120
121 // data converters
122 DataConverter *converter=0;
123
124 if(strstr(line.c_str(),".hep"))
125 {
126 cout<<"#**********************************************************************"<<endl;
127 cout<<"#********** StdHEP file format detected *************"<<endl;
128 cout<<"#*********** Starting convertion to TRoot format **************"<<endl;
129 cout<<"#**********************************************************************"<<endl;
130 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
131 }
132 else if(strstr(line.c_str(),".lhe"))
133 {
134 cout<<"#**********************************************************************"<<endl;
135 cout<<"#*********** LHEF file format detected ************"<<endl;
136 cout<<"#*********** Starting convertion to TRoot format ************"<<endl;
137 cout<<"#**********************************************************************"<<endl;
138 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
139 }
140 else if(strstr(line.c_str(),".root"))
141 {
142 cout<<"#**********************************************************************"<<endl;
143 cout<<"#********** h2root file format detected *************"<<endl;
144 cout<<"#********** Starting convertion to TRoot format *************"<<endl;
145 cout<<"#**********************************************************************"<<endl;
146 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
147 }
148 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
149
150 TChain chain("GEN");
151 chain.Add(outputfilename.c_str());
152 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
153 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
154 TIter itGen((TCollection*)branchGen);
155
156 //Output file : contents of the analysis object data
157 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
158 ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
159 ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
160 ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
161 ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
162 ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
163 ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
164 ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
165 ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
166 ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
167 ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
168 ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
169
170 TRootGenParticle *particle;
171 TRootETmis *elementEtmis;
172 TRootElectron *elementElec;
173 TRootMuon *elementMu;
174 TRootPhoton *elementPhoton;
175 TRootTracks *elementTracks;
176 TRootCalo *elementCalo;
177
178 TLorentzVector genMomentum(0,0,0,0); // four-momentum at the vertex
179 TLorentzVector genMomentumBfield(0,0,0,0); // four-momentum at the exit of the tracks
180 TLorentzVector momentumCaloSegmentation(0,0,0,0); // four-momentum in the calo, after applying the calo segmentation
181 LorentzVector jetMomentum;
182
183 vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
184 vector<fastjet::PseudoJet> sorted_jets;
185 vector<TLorentzVector> TrackCentral;
186 vector<PhysicsTower> towers;
187 vector<ParticleUtil> electron;
188 vector<ParticleUtil> muon;
189 vector<ParticleUtil> gamma;
190
191 TSimpleArray<TRootGenParticle> NFCentralQ;
192 float iPhi=0,iEta=0;
193
194
195
196// 2. ********** Loop over all events ***********
197 Long64_t entry, allEntries = treeReader->GetEntries();
198 cout << "** The input list contains " << allEntries << " events" << endl;
199 loopwatch.Start();
200
201 // loop on all events
202 for(entry = 0; entry < allEntries; ++entry)
203 {
204 TLorentzVector PTmis(0,0,0,0);
205 treeReader->ReadEntry(entry);
206 treeWriter->Clear();
207 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
208
209 electron.clear();
210 muon.clear();
211 gamma.clear();
212 NFCentralQ.Clear();
213
214 TrackCentral.clear();
215 towers.clear();
216 input_particles.clear();
217
218 // 2.1 Loop over all particles in event
219 itGen.Reset();
220 while( (particle = (TRootGenParticle*) itGen.Next()) ) {
221 int pid = abs(particle->PID);
222
223
224 // 2.1.1********************* preparation for the b-tagging
225 //// This subarray is needed for the B-jet algorithm
226 // optimization for speed : put first PID condition, then ETA condition, then either pt or status
227 if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
228 fabs(particle->Eta) < DET->CEN_max_tracker &&
229 particle->Status != 1 &&
230 particle->PT > DET->PT_QUARKS_MIN ) {
231 NFCentralQ.Add(particle);
232 }
233
234 // 2.1.2 ********************* central detector: keep only visible particles
235 // keeps only final particles, visible by the central detector, including the fiducial volume
236 // the ordering of conditions have been optimised for speed : put first the STATUS condition
237 if( (particle->Status == 1) &&
238 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
239 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
240 )
241 {
242 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
243 genMomentumBfield = genMomentum;
244
245 // 2.1.2.1 ********************* central detector: magnetic field
246 // genMomentumBfield is then changed with respect to the magnetic field
247 if(DET->FLAG_bfield==1) TRACP->Propagation(particle,genMomentumBfield);
248 float eta=fabs(genMomentumBfield.Eta());
249
250
251 // 2.1.2.2 ********************* central detector: smearing (calorimeters, muon chambers)
252 switch(pid) {
253
254 case pE: // all electrons with eta < DET->MAX_CALO_FWD
255 DET->SmearElectron(genMomentumBfield);
256 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_elec){
257 electron.push_back(ParticleUtil(genMomentumBfield,particle->PID));
258 }
259 break; // case pE
260 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
261 DET->SmearElectron(genMomentumBfield);
262 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_gamma) {
263 gamma.push_back(ParticleUtil(genMomentumBfield,particle->PID));
264 }
265 break; // case pGAMMA
266 case pMU: // all muons with eta < DET->MAX_MU
267 DET->SmearMu(genMomentumBfield);
268 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_mu && genMomentumBfield.Pt() > DET->PTCUT_muon){
269 muon.push_back(ParticleUtil(genMomentumBfield,particle->PID));
270 }
271 break; // case pMU
272 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
273 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
274 DET->SmearHadron(genMomentumBfield, 0.7);
275 break; // case hadron
276 default: // all other final particles with eta < DET->MAX_CALO_FWD
277 DET->SmearHadron(genMomentumBfield, 1.0);
278 break;
279 } // 2.1.2.2 switch (pid)
280
281
282 // 2.1.2.3 ********************* central detector: calotowers
283 // all final particles but muons and neutrinos
284 // for calorimetric towers and missing PT
285 int charge=Charge(pid);
286 if(genMomentumBfield.E() !=0 && pid != pMU) {
287 // in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
288 if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
289 else { // particles reach calos
290 // applies the calo segmentation and returns iEta & iPhi
291 DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta);
292 if(iEta != -100 && iPhi != -100) {
293 momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
294 elementCalo = (TRootCalo*) branchCalo->NewEntry();
295 elementCalo->Set(momentumCaloSegmentation);
296 PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
297 towers.push_back(Tower);
298 } // if iEta != -100
299 } // else : when particles reach the calos
300 } // 2.1.2.3 calotowers
301
302
303 // 2.1.2.4 ********************* central detector: tracks
304 // all final charged particles
305 if(
306 (genMomentumBfield.E()!=0) &&
307 (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) &&
308 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&
309 // if bfield not simulated, pt should be high enough to be taken into account
310 ((rand()%100) < DET->TRACK_eff) &&
311 (charge!=0)
312 )
313 {
314 elementTracks = (TRootTracks*) branchTracks->NewEntry();
315 elementTracks->Set(genMomentum); // fills px,py,pz,pt,e,eta,phi at vertex
316 elementTracks->Etaout = genMomentumBfield.Eta();
317 elementTracks->Phiout = genMomentumBfield.Phi();
318 // TODO!!! apply a smearing on the position of the origin of the track
319 // elementTracks->SetPosition(particle->X,particle->Y,particle->Z);
320 // uses the output of the bfield computation : Xout=... Yout=... Zout...
321 // TODO!!! elementTrakcs->SetPositionOut(Xout,Yout,Zout);
322 TrackCentral.push_back(genMomentum); // tracks at vertex!
323 } // 2.1.2.4 tracks
324
325 } // 2.1.2 central detector
326
327 // 2.1.3 ********************* forward detectors: zdc
328 if(DET->FLAG_vfd==1) {
329 VFD->ZDC(treeWriter,branchZDC,particle);
330 VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
331 }
332
333 } // 2.1 while : loop on all particles of the event.
334
335
336
337 // 2.2 ********** Output preparation & complex objects ***********
338
339 // 2.2.1 ********************* sorting collections by decreasing pt
340 DET->SortedVector(electron);
341 for(unsigned int i=0; i < electron.size(); i++) {
342 elementElec = (TRootElectron*) branchElectron->NewEntry();
343 elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E());
344 elementElec->Charge = sign(electron[i].PID());
345 elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks
346 }
347 DET->SortedVector(muon);
348 for(unsigned int i=0; i < muon.size(); i++) {
349 elementMu = (TRootMuon*) branchMuon->NewEntry();
350 elementMu->Charge = sign(muon[i].PID());
351 elementMu->Set(muon[i].Px(),muon[i].Py(),muon[i].Pz(),muon[i].E());
352 elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
353 }
354 DET->SortedVector(gamma);
355 for(unsigned int i=0; i < gamma.size(); i++) {
356 elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
357 elementPhoton->Set(gamma[i].Px(),gamma[i].Py(),gamma[i].Pz(),gamma[i].E());
358 }
359
360 // 2.2.2 ************* computes the Missing Transverse Momentum
361 TLorentzVector Att(0.,0.,0.,0.);
362 for(unsigned int i=0; i < towers.size(); i++)
363 {
364 Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
365 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
366 {
367 PTmis = PTmis + Att;
368 // create a fastjet::PseudoJet with these components and put it onto
369 // back of the input_particles vector
370 input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
371 }
372 }
373 elementEtmis = (TRootETmis*) branchETmis->NewEntry();
374 elementEtmis->ET = (PTmis).Pt();
375 elementEtmis->Phi = (-PTmis).Phi();
376 elementEtmis->Px = (-PTmis).Px();
377 elementEtmis->Py = (-PTmis).Py();
378
379 // 2.2.3 ************* B-tag, tau jets
380 sorted_jets=JETRUN->RunJets(input_particles);
381 JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
382 JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
383
384 treeWriter->Fill();
385 } // 2. Loop over all events ('for' loop)
386
387 treeWriter->Write();
388 delete treeWriter;
389 loopwatch.Stop();
390
391
392
393// 3. ********** Trigger & Frog ***********
394 // 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard
395 triggerwatch.Start();
396 if(DET->FLAG_trigger == 1)
397 {
398 // input
399 TChain chainT("Analysis");
400 chainT.Add(outputfilename.c_str());
401 ExRootTreeReader *treeReaderT = new ExRootTreeReader(&chainT);
402
403 // output
404 TClonesArray *branchElecTrig = treeReaderT->UseBranch("Electron");
405 TClonesArray *branchMuonTrig = treeReaderT->UseBranch("Muon");
406 TClonesArray *branchJetTrig = treeReaderT->UseBranch("Jet");
407 TClonesArray *branchTauJetTrig = treeReaderT->UseBranch("TauJet");
408 TClonesArray *branchPhotonTrig = treeReaderT->UseBranch("Photon");
409 TClonesArray *branchETmisTrig = treeReaderT->UseBranch("ETmis");
410
411 ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger");
412 ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class());
413
414
415 Long64_t entryT, allEntriesT = treeReaderT->GetEntries();
416 cout << "** Trigger: the 'Analysis' tree contains " << allEntriesT << " events" << endl;
417 // loop on all entries
418 for(entryT = 0; entryT < allEntriesT; ++entryT) {
419 treeWriterT->Clear();
420 treeReaderT->ReadEntry(entryT);
421 TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
422 treeWriterT->Fill();
423 } // loop on all entries
424
425 treeWriterT->Write();
426 delete treeWriterT;
427 } // trigger
428 triggerwatch.Stop();
429
430
431 // 3.2 ************** FROG display
432 frogwatch.Start();
433 if(DET->FLAG_frog == 1)
434 {
435 FrogDisplay *FROG = new FrogDisplay();
436 FROG->BuidEvents(outputfilename,DET->NEvents_Frog);
437 FROG->BuildGeom(DetDatacard);
438 }
439 frogwatch.Stop();
440
441
442
443
444// 4. ********** End & Exit ***********
445 cout << "** Exiting..." << endl;
446 globalwatch.Stop();
447 cout << "** Time report for " << allEntries << " events.\n";
448 cout << " + Time (s): \tCPU \t real"<< endl;
449 cout << " + Global: \t" << globalwatch.CpuTime() << " \t " << globalwatch.RealTime() << endl;
450 cout << " + Events: \t" << loopwatch.CpuTime() << " \t " << loopwatch.RealTime() << endl;
451 if(DET->FLAG_trigger == 1)
452 cout << " + Trigger: \t" << triggerwatch.CpuTime() << " \t " << triggerwatch.RealTime() << endl;
453 if(DET->FLAG_frog == 1)
454 cout << " + Frog: \t" << frogwatch.CpuTime() << " \t " << frogwatch.RealTime() << endl;
455
456
457 delete treeReader;
458 delete DET;
459 delete TRIGT;
460 delete TRACP;
461 delete JETRUN;
462 delete VFD;
463 if(converter) delete converter;
464
465// todo("TODO");
466}
Note: See TracBrowser for help on using the repository browser.