Fork me on GitHub

source: svn/trunk/Delphes.cpp@ 258

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

resolution from delphes output
File size: 28.6 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#include "TFile.h"
18
19#include "ExRootTreeReader.h"
20#include "ExRootTreeWriter.h"
21#include "ExRootTreeBranch.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 "BFieldProp.h"
30#include "TriggerUtil.h"
31#include "VeryForward.h"
32#include "JetsUtil.h"
33#include "FrogUtil.h"
34
35#include <vector>
36#include <iostream>
37
38#include "Utilities/ExRootAnalysis/interface/ExRootProgressBar.h"
39
40using namespace std;
41
42//------------------------------------------------------------------------------
43void todo(string filename) {
44 ifstream infile(filename.c_str());
45 cout << "** TODO list ..." << endl;
46 while(infile.good()) {
47 string temp;
48 getline(infile,temp);
49 cout << "*" << temp << endl;
50 }
51 cout << "** done...\n";
52}
53
54//------------------------------------------------------------------------------
55
56int main(int argc, char *argv[])
57{
58 int appargc = 2;
59 char *appName= new char[20];
60 char *appOpt= new char[20];
61 sprintf(appName,"Delphes");
62 sprintf(appOpt,"-b");
63 char *appargv[] = {appName,appOpt};
64 TApplication app(appName, &appargc, appargv);
65 delete [] appName;
66 delete [] appOpt;
67
68
69 if(argc != 3 && argc != 4 && argc != 5) {
70 cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
71 cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
72 cout << " output_file - output file." << endl;
73 cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
74 cout << " trigger_card - Datacard containing the trigger algorithms (optional) "<<endl;
75 exit(1);
76 }
77
78
79
80
81cout << endl << endl;
82
83cout <<"*********************************************************************"<< endl;
84cout <<"*********************************************************************"<< endl;
85cout <<"** **"<< endl;
86cout <<"** Welcome to **"<< endl;
87cout <<"** **"<< endl;
88cout <<"** **"<< endl;
89cout <<"** .ddddddd- lL hH **"<< endl;
90cout <<"** -Dd` `dD: Ll hH` **"<< endl;
91cout <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;
92cout <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;
93cout <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;
94cout <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;
95cout <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;
96cout <<"** Pp **"<< endl;
97cout <<"** **"<< endl;
98cout <<"** Delphes, a framework for the fast simulation **"<< endl;
99cout <<"** of a generic collider experiment **"<< endl;
100cout <<"** **"<< endl;
101cout <<"** --- Version 1.3beta of Delphes --- **"<< endl;
102cout <<"** Last date of change: 29 January 2009 **"<< endl;
103cout <<"** **"<< endl;
104cout <<"** **"<< endl;
105cout <<"** This package uses: **"<< endl;
106cout <<"** ------------------ **"<< endl;
107cout <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
108cout <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;
109cout <<"** FROG: **"<< endl;
110cout <<"** **"<< endl;
111cout <<"**-----------------------------------------------------------------**"<< endl;
112cout <<"** **"<< endl;
113cout <<"** Main authors: **"<< endl;
114cout <<"** ------------- **"<< endl;
115cout <<"** **"<< endl;
116cout <<"** Séverine Ovyn Xavier Rouby **"<< endl;
117cout <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;
118cout <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;
119cout <<"** Universite Catholique de Louvain (UCL) **"<< endl;
120cout <<"** Louvain-la-Neuve, Belgium **"<< endl;
121cout <<"** **"<< endl;
122cout <<"**-----------------------------------------------------------------**"<< endl;
123cout <<"** **"<< endl;
124cout <<"** Former Delphes versions and documentation can be found on : **"<< endl;
125cout <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;
126cout <<"** **"<< endl;
127cout <<"** **"<< endl;
128cout <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl;
129cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
130cout <<"** give us your feedbacks about potential bugs **"<< endl;
131cout <<"** **"<< endl;
132cout <<"*********************************************************************"<< endl;
133cout <<"*********************************************************************"<< endl;
134
135// 1. ********** initialisation ***********
136
137 srand (time (NULL)); /* Initialisation du générateur */
138 TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch;
139 globalwatch.Start();
140
141
142 //read the output TROOT file
143 string inputFileList(argv[1]), outputfilename(argv[2]);
144 if(outputfilename.find(".root") > outputfilename.length()) {
145 cout <<"** ERROR: 'output_file' should be a .root file. Exiting... **"<< endl;
146 exit(1);
147 }
148 //create output log-file name
149 string forLog = outputfilename;
150 string LogName = forLog.erase(forLog.find(".root"));
151 LogName = LogName+"_run.log";
152
153 TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
154 outputFile->Close();
155
156 string line;
157 ifstream infile(inputFileList.c_str());
158 infile >> line; // the first line determines the type of input files
159
160 //read the datacard input file
161 string DetDatacard(""); //for detector smearing parameters
162 string TrigDatacard("data/trigger.dat"); //for trigger selection
163
164 string lineCard1,lineCard2;
165 bool detecCard=false,trigCard=false;
166 if(argv[3])
167 {
168 ifstream infile1(argv[3]);
169 infile1 >> lineCard1; // the first line determines the type of input files
170 if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==true)
171 cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl;
172 else if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[3]; detecCard=true;}
173 else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==true)
174 cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl;
175 else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[3]; trigCard=true;}
176 }
177 if(argv[4])
178 {
179 ifstream infile2(argv[4]);
180 infile2 >> lineCard2; // the first line determines the type of input files
181 if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==true)
182 cerr <<"** ERROR: A DETECTOR card has already been loaded **"<< endl;
183 else if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[4]; detecCard=true;}
184 else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==true)
185 cerr <<"** ERROR: A TRIGGER card has already been loaded **"<< endl;
186 else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[4]; trigCard=true;}
187 }
188
189 //Smearing information
190 RESOLution *DET = new RESOLution();
191 cout <<"** **"<< endl;
192 cout <<"** ####### Start reading DETECTOR parameters ####### **"<< endl;
193 cout << left << setw(40) <<"** Opening configuration card: "<<""
194 << left << setw(27) << DetDatacard <<""
195 << right << setw(2) <<"**"<<""<<endl;
196 DET->ReadDataCard(DetDatacard);
197 DET->Logfile(LogName);
198 cout << left << setw(42) <<"** Parameters summarised in file: "<<""
199 << left << setw(27) << LogName <<""
200 << right << setw(2) <<"**"<<""<<endl;
201 cout <<"** **"<< endl;
202
203 //Trigger information
204 cout <<"** ########### Start reading TRIGGER card ########## **"<< endl;
205 if(trigCard==false)
206 {
207 cout <<"** WARNING: Datadard not found, use default card **" << endl;
208 TrigDatacard="data/trigger.dat";
209 }
210 TriggerTable *TRIGT = new TriggerTable();
211 TRIGT->TriggerCardReader(TrigDatacard.c_str());
212 TRIGT->PrintTriggerTable(LogName);
213 if(DET->FLAG_trigger == 1)
214 {
215 cout << left << setw(40) <<"** Opening configuration card: "<<""
216 << left << setw(27) << TrigDatacard <<""
217 << right << setw(2) <<"**"<<""<<endl;
218 cout <<"** **"<< endl;
219 }
220
221 //Propagation of tracks in the B field
222 TrackPropagation *TRACP = new TrackPropagation(DET);
223 //TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
224
225 //Jet information
226 JetsUtil *JETRUN = new JetsUtil(DET);
227 //JetsUtil *JETRUN = new JetsUtil(DetDatacard);
228
229 //VFD information
230 VeryForward * VFD = new VeryForward(DET);
231 //VeryForward * VFD = new VeryForward(DetDatacard);
232
233 // data converters
234 DataConverter *converter=NULL;
235 cout <<"** **"<<endl;
236 cout <<"** ####### Start convertion to TRoot format ######## **"<< endl;
237
238 if(strstr(line.c_str(),".hep"))
239 {
240 cout <<"** StdHEP file format detected **"<<endl;
241 cout <<"** This can take several minutes **"<< endl;
242 converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
243 }
244 else if(strstr(line.c_str(),".lhe"))
245 {
246 cout <<"** LHEF file format detected **"<<endl;
247 cout <<"** This can take several minutes **"<< endl;
248 converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
249 }
250 else if(strstr(line.c_str(),".root"))
251 {
252 cout <<"** h2root file format detected **"<<endl;
253 cout <<"** This can take several minutes **"<< endl;
254 converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
255 }
256 else {
257 cerr << left << setw(4) <<"** "<<""
258 << left << setw(63) << line.c_str() <<""
259 << right << setw(2) <<"**"<<endl;
260 cerr <<"** ERROR: File format not identified -- Exiting... **"<< endl;
261 cout <<"** **"<< endl;
262 cout <<"*********************************************************************"<< endl;
263 return -1;};
264 cout <<"** Exiting conversion... **"<< endl;
265
266 TChain chain("GEN");
267 chain.Add(outputfilename.c_str());
268 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
269 const TClonesArray *branchGen = treeReader->UseBranch("Particle");
270
271 TIter itGen((TCollection*)branchGen);
272
273 //Output file : contents of the analysis object data
274 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
275 ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
276 ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
277 ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
278 ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
279 ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
280 ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
281 ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
282 ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
283 ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
284 ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
285 ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
286
287 TRootGenParticle *particle;
288 TRootETmis *elementEtmis;
289 TRootElectron *elementElec;
290 TRootMuon *elementMu;
291 TRootPhoton *elementPhoton;
292 TRootTracks *elementTracks;
293 TRootCalo *elementCalo;
294
295 TLorentzVector genMomentum(0,0,0,0); // four-momentum at the vertex
296 TLorentzVector genMomentumBfield(0,0,0,0); // four-momentum at the exit of the tracks
297 TLorentzVector momentumCaloSegmentation(0,0,0,0); // four-momentum in the calo, after applying the calo segmentation
298 LorentzVector jetMomentum;
299
300 vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
301 vector<fastjet::PseudoJet> sorted_jets;
302 vector<TLorentzVector> TrackCentral;
303 vector<PhysicsTower> towers;
304 vector<ParticleUtil> electron;
305 vector<ParticleUtil> muon;
306 vector<ParticleUtil> gamma;
307
308 TSimpleArray<TRootGenParticle> NFCentralQ;
309 float iPhi=0,iEta=0;
310
311
312
313// 2. ********** Loop over all events ***********
314 Long64_t entry, allEntries = treeReader->GetEntries();
315 cout <<"** **"<<endl;
316 cout <<"** ####### Start fast detector simulation ######## **"<< endl;
317 cout << left << setw(52) <<"** Total number of events to run: "<<""
318 << left << setw(15) << allEntries <<""
319 << right << setw(2) <<"**"<<endl;
320
321 ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
322
323 loopwatch.Start();
324
325 // loop on all events
326 for(entry = 0; entry < allEntries; ++entry)
327 {
328 Progress->Update(entry);
329 TLorentzVector PTmis(0,0,0,0);
330 treeReader->ReadEntry(entry);
331 treeWriter->Clear();
332 if((entry % 100) == 0 && entry > 0 ) {
333 // cout << left << setw(52) <<"** Processing element # "<<""
334 // << left << setw(15) << entry <<""
335 // << right << setw(2) <<"**"<<endl;
336 }
337
338 electron.clear();
339 muon.clear();
340 gamma.clear();
341 NFCentralQ.Clear();
342
343 TrackCentral.clear();
344 towers.clear();
345 input_particles.clear();
346
347 // 2.1 Loop over all particles in event
348 itGen.Reset();
349 while( (particle = (TRootGenParticle*) itGen.Next()) ) {
350 int pid = abs(particle->PID);
351
352
353 // 2.1.1********************* preparation for the b-tagging
354 //// This subarray is needed for the B-jet algorithm
355 // optimization for speed : put first PID condition, then ETA condition, then either pt or status
356 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 ?
357 fabs(particle->Eta) < DET->CEN_max_tracker &&
358 particle->Status != 1 &&
359 particle->PT > DET->PT_QUARKS_MIN ) {
360 NFCentralQ.Add(particle);
361 }
362
363 // 2.1.2 ********************* central detector: keep only visible particles
364 // keeps only final particles, visible by the central detector, including the fiducial volume
365 // the ordering of conditions have been optimised for speed : put first the STATUS condition
366 if( (particle->Status == 1) &&
367 ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
368 (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
369 )
370 {
371 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
372 genMomentumBfield = genMomentum;
373
374 // 2.1.2.1 ********************* central detector: magnetic field
375 // genMomentumBfield is then changed with respect to the magnetic field
376 if(DET->FLAG_bfield==1) TRACP->Propagation(particle,genMomentumBfield);
377 float eta=fabs(genMomentumBfield.Eta());
378
379
380 // 2.1.2.2 ********************* central detector: smearing (calorimeters, muon chambers)
381 switch(pid) {
382
383 case pE: // all electrons with eta < DET->MAX_CALO_FWD
384 DET->SmearElectron(genMomentumBfield);
385 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_elec){
386 electron.push_back(ParticleUtil(genMomentumBfield,particle->PID));
387 }
388 break; // case pE
389 case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
390 DET->SmearElectron(genMomentumBfield);
391 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_gamma) {
392 gamma.push_back(ParticleUtil(genMomentumBfield,particle->PID));
393 }
394 break; // case pGAMMA
395 case pMU: // all muons with eta < DET->MAX_MU
396 DET->SmearMu(genMomentumBfield);
397 if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_mu && genMomentumBfield.Pt() > DET->PTCUT_muon){
398 muon.push_back(ParticleUtil(genMomentumBfield,particle->PID));
399 }
400 break; // case pMU
401 case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
402 case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
403 DET->SmearHadron(genMomentumBfield, 0.7);
404 break; // case hadron
405 default: // all other final particles with eta < DET->MAX_CALO_FWD
406 DET->SmearHadron(genMomentumBfield, 1.0);
407 break;
408 } // 2.1.2.2 switch (pid)
409
410
411 // 2.1.2.3 ********************* central detector: calotowers
412 // all final particles but muons and neutrinos
413 // for calorimetric towers and missing PT
414 int charge=Charge(pid);
415 if(genMomentumBfield.E() !=0 && pid != pMU) {
416 // in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
417 if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
418 else { // particles reach calos
419 // applies the calo segmentation and returns iEta & iPhi
420 DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta);
421 if(iEta != -100 && iPhi != -100) {
422 momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
423 elementCalo = (TRootCalo*) branchCalo->NewEntry();
424 elementCalo->Set(momentumCaloSegmentation);
425 PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
426 towers.push_back(Tower);
427 } // if iEta != -100
428 } // else : when particles reach the calos
429 } // 2.1.2.3 calotowers
430
431
432 // 2.1.2.4 ********************* central detector: tracks
433 // all final charged particles
434 if(
435 (genMomentumBfield.E()!=0) &&
436 (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) &&
437 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&
438 // if bfield not simulated, pt should be high enough to be taken into account
439 ((rand()%100) < DET->TRACK_eff) &&
440 (charge!=0)
441 )
442 {
443 elementTracks = (TRootTracks*) branchTracks->NewEntry();
444 elementTracks->Set(genMomentum); // fills px,py,pz,pt,e,eta,phi at vertex
445 elementTracks->Etaout = genMomentumBfield.Eta();
446 elementTracks->Phiout = genMomentumBfield.Phi();
447 // TODO!!! apply a smearing on the position of the origin of the track
448 // elementTracks->SetPosition(particle->X,particle->Y,particle->Z);
449 // uses the output of the bfield computation : Xout=... Yout=... Zout...
450 // TODO!!! elementTrakcs->SetPositionOut(Xout,Yout,Zout);
451 TrackCentral.push_back(genMomentum); // tracks at vertex!
452 } // 2.1.2.4 tracks
453
454 } // 2.1.2 central detector
455
456 // 2.1.3 ********************* forward detectors: zdc
457 if(DET->FLAG_vfd==1) {
458 VFD->ZDC(treeWriter,branchZDC,particle);
459 VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
460 }
461
462 } // 2.1 while : loop on all particles of the event.
463
464
465
466 // 2.2 ********** Output preparation & complex objects ***********
467
468 // 2.2.1 ********************* sorting collections by decreasing pt
469 DET->SortedVector(electron);
470 for(unsigned int i=0; i < electron.size(); i++) {
471 elementElec = (TRootElectron*) branchElectron->NewEntry();
472 elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E());
473 elementElec->Charge = sign(electron[i].PID());
474 elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks
475 }
476 DET->SortedVector(muon);
477 for(unsigned int i=0; i < muon.size(); i++) {
478 elementMu = (TRootMuon*) branchMuon->NewEntry();
479 elementMu->Charge = sign(muon[i].PID());
480 elementMu->Set(muon[i].Px(),muon[i].Py(),muon[i].Pz(),muon[i].E());
481 elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
482 }
483 DET->SortedVector(gamma);
484 for(unsigned int i=0; i < gamma.size(); i++) {
485 elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
486 elementPhoton->Set(gamma[i].Px(),gamma[i].Py(),gamma[i].Pz(),gamma[i].E());
487 }
488
489 // 2.2.2 ************* computes the Missing Transverse Momentum
490 TLorentzVector Att(0.,0.,0.,0.);
491 for(unsigned int i=0; i < towers.size(); i++)
492 {
493 Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
494 if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
495 {
496 PTmis = PTmis + Att;
497 // create a fastjet::PseudoJet with these components and put it onto
498 // back of the input_particles vector
499 input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
500 }
501 }
502 elementEtmis = (TRootETmis*) branchETmis->NewEntry();
503 elementEtmis->ET = (PTmis).Pt();
504 elementEtmis->Phi = (-PTmis).Phi();
505 elementEtmis->Px = (-PTmis).Px();
506 elementEtmis->Py = (-PTmis).Py();
507
508 // 2.2.3 ************* B-tag, tau jets
509 sorted_jets=JETRUN->RunJets(input_particles);
510 JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
511 JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
512
513 treeWriter->Fill();
514 } // 2. Loop over all events ('for' loop)
515
516 cout <<"** **"<< endl;
517 cout <<"** Exiting detector simulation... **"<< endl;
518
519
520 treeWriter->Write();
521 delete treeWriter;
522 loopwatch.Stop();
523
524
525
526 // 3. ********** Trigger & Frog ***********
527 // 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard
528 triggerwatch.Start();
529 if(DET->FLAG_trigger == 1)
530 {
531 cout <<"** **"<<endl;
532 cout <<"** ########### Start Trigger selection ########### **"<< endl;
533
534 // input
535 TChain chainT("Analysis");
536 chainT.Add(outputfilename.c_str());
537 ExRootTreeReader *treeReaderT = new ExRootTreeReader(&chainT);
538
539 // output
540 TClonesArray *branchElecTrig = treeReaderT->UseBranch("Electron");
541 TClonesArray *branchMuonTrig = treeReaderT->UseBranch("Muon");
542 TClonesArray *branchJetTrig = treeReaderT->UseBranch("Jet");
543 TClonesArray *branchTauJetTrig = treeReaderT->UseBranch("TauJet");
544 TClonesArray *branchPhotonTrig = treeReaderT->UseBranch("Photon");
545 TClonesArray *branchETmisTrig = treeReaderT->UseBranch("ETmis");
546
547 ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger");
548 ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class());
549
550
551 Long64_t entryT, allEntriesT = treeReaderT->GetEntries();
552 // loop on all entries
553 for(entryT = 0; entryT < allEntriesT; ++entryT) {
554 treeWriterT->Clear();
555 treeReaderT->ReadEntry(entryT);
556 TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
557 treeWriterT->Fill();
558 } // loop on all entries
559 cout <<"** Exiting trigger simulation... **"<< endl;
560
561 treeWriterT->Write();
562 delete treeWriterT;
563 delete treeReaderT;
564 } // trigger
565 triggerwatch.Stop();
566
567
568 // 3.2 ************** FROG display
569 frogwatch.Start();
570 if(DET->FLAG_frog == 1) {
571 cout <<"** **"<<endl;
572 cout <<"** ################## Start FROG ################# **"<< endl;
573
574 FrogDisplay *FROG = new FrogDisplay(DET);
575 FROG->BuildEvents(outputfilename);
576 FROG->BuildGeom();
577 delete FROG;
578 }
579 frogwatch.Stop();
580
581
582
583
584// 4. ********** End & Exit ***********
585
586 globalwatch.Stop();
587 cout <<"** **"<< endl;
588 cout <<"** ################## Time report ################# **"<< endl;
589 cout << left << setw(32) <<"** Time report for "<<""
590 << left << setw(15) << allEntries <<""
591 << right << setw(22) <<"events **"<<endl;
592 cout <<"** **"<< endl;
593 cout << left << setw(10) <<"**"<<""
594 << left << setw(15) <<"Time (s):"<<""
595 << right << setw(15) <<"CPU"<<""
596 << right << setw(15) <<"Real"<<""
597 << right << setw(14) <<"**"<<endl;
598 cout << left << setw(10) <<"**"<<""
599 << left << setw(15) <<" + Global:"<<""
600 << right << setw(15) <<globalwatch.CpuTime()<<""
601 << right << setw(15) <<globalwatch.RealTime()<<""
602 << right << setw(14) <<"**"<<endl;
603 cout << left << setw(10) <<"**"<<""
604 << left << setw(15) <<" + Events:"<<""
605 << right << setw(15) <<loopwatch.CpuTime()<<""
606 << right << setw(15) <<loopwatch.RealTime()<<""
607 << right << setw(14) <<"**"<<endl;
608 if(DET->FLAG_trigger == 1)
609 {
610 cout << left << setw(10) <<"**"<<""
611 << left << setw(15) <<" + Trigger:"<<""
612 << right << setw(15) <<triggerwatch.CpuTime()<<""
613 << right << setw(15) <<triggerwatch.RealTime()<<""
614 << right << setw(14) <<"**"<<endl;
615 }
616 if(DET->FLAG_frog == 1)
617 {
618 cout << left << setw(10) <<"**"<<""
619 << left << setw(15) <<" + Frog:"<<""
620 << right << setw(15) <<frogwatch.CpuTime()<<""
621 << right << setw(15) <<frogwatch.RealTime()<<""
622 << right << setw(14) <<"**"<<endl;
623
624 }
625
626 cout <<"** **"<< endl;
627 cout <<"** Exiting Delphes ... **"<< endl;
628 cout <<"** **"<< endl;
629 cout <<"*********************************************************************"<< endl;
630 cout <<"*********************************************************************"<< endl;
631
632
633 delete treeReader;
634 delete DET;
635// delete TRIGT;
636 delete TRACP;
637 delete JETRUN;
638 delete VFD;
639 delete converter;
640
641// todo("TODO");
642}
Note: See TracBrowser for help on using the repository browser.