Fork me on GitHub

Context Navigation

source: svn/trunk/Delphes.cpp@ 457

Last change on this file since 457 was 457, checked in by Xavier Rouby, 15 years ago
new class for inputing directly rootfiles from Delphes (with GEN tree)
File size: 35.3 KB

Line
1	/***********************************************************************
2
3	/----------------------------------------------\
4	\| Delphes, a framework for the fast simulation \|
5	\| of a generic collider experiment \|
6	\------------- arXiv:0903.2225v1 ------------/
7
8
9	This package uses:
10	------------------
11	ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86
12	FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]
13	Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]
14	FROG: [hep-ex/0901.2718v1]
15	HepMC: Comput. Phys. Commun.134 (2001) 41
16
17	------------------------------------------------------------------
18
19	Main authors:
20	-------------
21
22	Severine Ovyn Xavier Rouby
23	severine.ovyn@uclouvain.be xavier.rouby@cern
24
25	Center for Particle Physics and Phenomenology (CP3)
26	Universite catholique de Louvain (UCL)
27	Louvain-la-Neuve, Belgium
28
29	Copyright (C) 2008-2009,
30	All rights reserved.
31
32	***********************************************************************/
33
34	/// \file Delphes.cpp
35	/// \brief Executable for Delphes
36
37	#include "TChain.h"
38	#include "TApplication.h"
39	#include "TStopwatch.h"
40	#include "TFile.h"
41
42	#include "ExRootTreeReader.h"
43	#include "ExRootTreeWriter.h"
44	#include "ExRootTreeBranch.h"
45	#include "ExRootProgressBar.h"
46
47	#include "DataConverter.h"
48	#include "LHEFConverter.h"
49	#include "HepMCConverter.h"
50	#include "HEPEVTConverter.h"
51	#include "STDHEPConverter.h"
52	#include "LHCOConverter.h"
53	#include "DelphesRootConverter.h"
54
55	#include "SmearUtil.h"
56	#include "CaloUtil.h"
57	#include "BFieldProp.h"
58	#include "TriggerUtil.h"
59	#include "VeryForward.h"
60	#include "JetsUtil.h"
61	#include "FrogUtil.h"
62
63	#include <vector>
64	#include <iostream>
65
66	using namespace std;
67
68	//------------------------------------------------------------------------------
69	void todo(string filename) {
70	ifstream infile(filename.c_str());
71	cout << "** TODO list ..." << endl;
72	while(infile.good()) {
73	string temp;
74	getline(infile,temp);
75	cout << "*" << temp << endl;
76	}
77	cout << "** done...\n";
78	}
79
80	//------------------------------------------------------------------------------
81
82	int main(int argc, char *argv[])
83	{
84
85	int appargc = 2;
86	char *appName= new char[20];
87	char *appOpt= new char[20];
88	sprintf(appName,"Delphes");
89	sprintf(appOpt,"-b");
90	char *appargv[] = {appName,appOpt};
91	TApplication app(appName, &appargc, appargv);
92	delete [] appName;
93	delete [] appOpt;
94
95	if(argc != 3 && argc != 4 && argc != 5) {
96	cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
97	cout << " input_list - list of files in Ntpl, StdHep, HepMC or LHEF format," << endl;
98	cout << " output_file - output file." << endl;
99	cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
100	cout << " trigger_card - Datacard containing the trigger algorithms (optional) "<<endl;
101	exit(1);
102	}
103
104	print_header();
105
106	// 1. ******** initialisation *********
107
108	srand (time (NULL)); /* Initialisation du gÃ©nÃ©rateur */
109	TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch, lhcowatch;
110	globalwatch.Start();
111
112
113	//read the output TROOT file
114	string inputFileList(argv[1]), outputfilename(argv[2]);
115	if(outputfilename.find(".root") > outputfilename.length()) {
116	cout <<" ERROR: 'output_file' should be a .root file. Exiting... "<< endl;
117	exit(1);
118	}
119	//create output log-file name
120	string forLog = outputfilename;
121	string LogName = forLog.erase(forLog.find(".root"));
122	LogName = LogName+"_run.log";
123
124	TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
125	outputFile->Close();
126
127	string line;
128	ifstream infile(inputFileList.c_str());
129	if(!infile.good()) {
130	cout << " ERROR: Input list (" << left << setw(13) << inputFileList << ") not found. Exiting... "<< endl;
131	cout <<"*********************************************************************"<< endl;
132	exit(1);
133	}
134	infile >> line; // the first line determines the type of input files
135
136	//read the datacard input file
137	string DetDatacard("data/DetectorCard.dat"); //for detector smearing parameters
138	string TrigDatacard("data/TriggerCard.dat"); //for trigger selection
139
140	string lineCard1,lineCard2;
141	bool detecCard=false,trigCard=false;
142	if(argv[3])
143	{
144	ifstream infile1(argv[3]);
145	infile1 >> lineCard1; // the first line determines the type of input files
146	if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==true)
147	cerr <<" ERROR: A DETECTOR card has already been loaded "<< endl;
148	else if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[3]; detecCard=true;}
149	else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==true)
150	cerr <<" ERROR: A TRIGGER card has already been loaded "<< endl;
151	else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[3]; trigCard=true;}
152	}
153	if(argv[4])
154	{
155	ifstream infile2(argv[4]);
156	infile2 >> lineCard2; // the first line determines the type of input files
157	if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==true)
158	cerr <<" ERROR: A DETECTOR card has already been loaded "<< endl;
159	else if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[4]; detecCard=true;}
160	else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==true)
161	cerr <<" ERROR: A TRIGGER card has already been loaded "<< endl;
162	else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[4]; trigCard=true;}
163	}
164
165	//Smearing information
166	RESOLution *DET = new RESOLution();
167
168	cout <<" "<< endl;
169	cout <<" ####### Start reading DETECTOR parameters ####### "<< endl;
170	cout << left << setw(40) <<"** Opening configuration card: "<<""
171	<< left << setw(27) << DetDatacard <<""
172	<< right << setw(2) <<"**"<<""<<endl;
173	DET->ReadDataCard(DetDatacard);
174	cout << left << setw(40) <<"** Parameters summarised in: "<<""
175	<< left << setw(27) << LogName <<""
176	<< right << setw(2) <<"**"<<""<<endl;
177	cout <<" "<< endl;
178	DET->ReadParticleDataGroupTable();
179	// DET->PDGtable.print();
180
181	//Trigger information
182	cout <<" ########### Start reading TRIGGER card ########## "<< endl;
183	if(trigCard==false)
184	{
185	cout <<" WARNING: Datacard not found, use default card " << endl;
186	TrigDatacard="data/TriggerCard.dat";
187	}
188	TriggerTable *TRIGT = new TriggerTable();
189	TRIGT->TriggerCardReader(TrigDatacard.c_str());
190	TRIGT->PrintTriggerTable(LogName);
191	if(DET->FLAG_trigger == 1)
192	{
193	cout << left << setw(40) <<"** Opening configuration card: "<<""
194	<< left << setw(27) << TrigDatacard <<""
195	<< right << setw(2) <<"**"<<""<<endl;
196	cout <<" "<< endl;
197	}
198
199	// Logfile
200	DET->setNames(inputFileList,DetDatacard,TrigDatacard);
201	DET->Logfile(LogName);
202
203	//Propagation of tracks in the B field
204	TrackPropagation *TRACP = new TrackPropagation(DET);
205
206	//Jet information
207	JetsUtil *JETRUN = new JetsUtil(DET);
208
209	//VFD information
210	VeryForward * VFD = new VeryForward(DET);
211
212	// data converters
213	cout <<" "<<endl;
214	cout <<" ####### Start conversion to TRoot format ######## "<< endl;
215
216	if(line.rfind(".hepmc") < line.length())
217	{
218	cout <<" HepMC ASCII file format detected "<<endl;
219	cout <<" This can take several minutes "<< endl;
220	HepMCConverter converter(inputFileList,outputfilename,DET->PDGtable,DET->NEvents);
221	}
222	else if(line.rfind(".hep") < line.length())
223	{
224	cout <<" StdHEP file format detected "<<endl;
225	cout <<" This can take several minutes "<< endl;
226	STDHEPConverter converter(inputFileList,outputfilename,DET->PDGtable,DET->NEvents);
227	}
228	else if(line.rfind(".lhe") < line.length())
229	{
230	cout <<" LHEF file format detected "<<endl;
231	cout <<" This can take several minutes "<< endl;
232	LHEFConverter converter(inputFileList,outputfilename,DET->PDGtable,DET->NEvents);
233	}
234	else if(line.rfind(".root") < line.length())
235	// can be either a root file from h2root (i.e. with "h101" tree)
236	// or a root file from Delphes (i.e. with "GEN" tree)
237	{
238	TFile f(line.c_str());
239	if (f.FindKey("GEN")) {
240	cout <<" Delphes ROOT file format detected "<<endl;
241	cout <<" This can take several minutes "<< endl;
242	DelphesRootConverter converter(inputFileList,outputfilename,DET->NEvents);
243	}
244	else
245	if (f.FindKey("h101")) {
246	cout <<" h2root file format detected "<<endl;
247	cout <<" This can take several minutes "<< endl;
248	HEPEVTConverter converter(inputFileList,outputfilename,DET->PDGtable,DET->NEvents);
249	}
250	else {
251	cerr << left << setw(4) <<"** "<<""
252	<< left << setw(63) << line.c_str() <<""
253	<< right << setw(2) <<"**"<<endl;
254	cerr <<" ERROR: File format not identified -- Exiting... "<< endl;
255	cout <<" "<< endl;
256	cout <<"*********************************************************************"<< endl;
257	return -1;
258	} // not found any interesting input tree
259	f.Close();
260	} // .root file
261	else {
262	cerr << left << setw(4) <<"** "<<""
263	<< left << setw(63) << line.c_str() <<""
264	<< right << setw(2) <<"**"<<endl;
265	cerr <<" ERROR: File format not identified -- Exiting... "<< endl;
266	cout <<" "<< endl;
267	cout <<"*********************************************************************"<< endl;
268	return -1;};
269	cout <<" Exiting conversion... "<< endl;
270
271	TChain chain("GEN");
272	chain.Add(outputfilename.c_str());
273	ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
274	const TClonesArray *branchGen = treeReader->UseBranch("Particle");
275
276	TIter itGen((TCollection*)branchGen);
277
278	//Output file : contents of the analysis object data
279	ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
280	ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
281	ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
282	ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
283	ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
284	ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
285	ExRootTreeBranch *branchTrack = treeWriter->NewBranch("Tracks", TRootTracks::Class());
286	ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
287	ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
288	ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
289	ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
290	//ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootForwardTaggerHits::Class());
291	ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
292
293	TRootETmis *elementEtmis;
294	TRootElectron *elementElec;
295	TRootMuon *elementMu;
296	TRootPhoton *elementPhoton;
297	TRootTracks * elementTrack;
298	TRootCalo *elementCalo;
299
300	TLorentzVector genMomentum(0,0,0,0); // four-momentum at the vertex
301	TLorentzVector genMomentumBfield(0,0,0,0); // four-momentum at the exit of the tracks
302	TLorentzVector momentumCaloSegmentation(0,0,0,0); // four-momentum in the calo, after applying the calo segmentation
303	LorentzVector jetMomentum;
304
305	vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
306	vector<fastjet::PseudoJet> sorted_jets;
307	vector<TRootTracks> TrackCentral;
308	vector<PhysicsTower> towers;
309	vector<D_Particle> electron;
310	vector<D_Particle> muon;
311	vector<D_Particle> gamma;
312
313	vector<int> NTrackJet;
314
315	TSimpleArray<TRootC::GenParticle> NFCentralQ;
316
317	D_CaloList list_of_calorimeters;
318	D_CaloElement CentralCalo("centralcalo",
319	-DET->CEN_max_calo_cen, DET->CEN_max_calo_cen,
320	DET->ELG_Ccen, DET->ELG_Ncen, DET->ELG_Scen,
321	DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal);
322	D_CaloElement ForwardCalo("forwardcalo",
323	DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd,
324	DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd,
325	DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf );
326	D_CaloElement BackwardCalo("backwardcalo",
327	-DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen,
328	DET->ELG_Cfwd, DET->ELG_Nfwd, DET->ELG_Sfwd,
329	DET->HAD_Chf, DET->HAD_Nhf, DET->HAD_Shf );
330	//D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0);
331	list_of_calorimeters.addElement(CentralCalo);
332	list_of_calorimeters.addElement(ForwardCalo);
333	list_of_calorimeters.addElement(BackwardCalo);
334	//list_of_calorimeters.addElement(CastorCalo);
335	list_of_calorimeters.sortElements();
336
337
338	// 2. ******** Loop over all events *********
339	Long64_t entry, allEntries = treeReader->GetEntries();
340	cout <<" "<<endl;
341	cout <<" ####### Start fast detector simulation ######## "<< endl;
342	cout << left << setw(52) <<"** Total number of events to run: "<<""
343	<< left << setw(15) << allEntries <<""
344	<< right << setw(2) <<"**"<<endl;
345
346	ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
347
348	loopwatch.Start();
349
350	// loop on all events
351	for(entry = 0; entry < allEntries; ++entry)
352	{
353	Progress->Update(entry);
354	TLorentzVector PTmis(0,0,0,0);
355	treeReader->ReadEntry(entry);
356	treeWriter->Clear();
357
358	electron.clear();
359	muon.clear();
360	gamma.clear();
361	NFCentralQ.Clear();
362
363	TrackCentral.clear();
364	towers.clear();
365	input_particles.clear();
366	NTrackJet.clear();
367
368	// 'list_of_active_towers' contains the exact list of calorimetric towers which have some deposits inside (E>0).
369	// The towers of this list will be smeared according to the calo resolution, afterwards
370	D_CaloTowerList list_of_active_towers;
371
372	// 'list_of_towers_with_photon' and 'list_of_centowers_with_neutrals' are list of towers, whose energy is not computed.
373	// They are only used to store the eta/phi of some towers, in order to search later inside 'list_of_active_towers'.
374	// 'list_of_towers_with_photon' contains the towers hit by photons only
375	// 'list_of_centowers_with_neutrals' is used to the jet-E-flow calculation: contains the towers with eta < CEN_max_tracker,
376	// i.e. towers behind the tracker.
377	D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates
378
379	D_CaloTowerList list_of_centowers_with_neutrals; // list of towers with neutral particles : for jet E-flow
380	float etamax_calocoverage_behindtracker = DET->CEN_max_tracker; // finds the extension in eta of the furthest
381	for (unsigned int i=1; i< DET->TOWER_number+1; i++) { // cell (at least) partially behind the tracker
382	if(DET->TOWER_eta_edges[i] > DET->CEN_max_tracker) break;
383	etamax_calocoverage_behindtracker = DET->TOWER_eta_edges[i];
384	}
385	// 2.1a Loop over all particles in event, to fill the towers
386	itGen.Reset();
387	TRootC::GenParticle *particleG;
388	while( (particleG = (TRootC::GenParticle*) itGen.Next()) )
389	{
390	TRootGenParticle *particle = new TRootGenParticle(particleG);
391	PdgParticle pdg_part = DET->PDGtable[particle->PID];
392	particle->Charge = pdg_part.charge();
393	particle->M = pdg_part.mass();
394	//particle->Charge=ChargeVal(particle->PID);
395	particle->setFractions(); // init
396	int pid = abs(particle->PID);
397
398	// 2.1a.1********************* preparation for the b-tagging
399	//// This subarray is needed for the B-jet algorithm
400	// optimization for speed : put first PID condition, then ETA condition, then either pt or status
401	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 ?
402	fabs(particle->Eta) < DET->CEN_max_tracker &&
403	particle->Status != 1 &&
404	particle->PT > DET->PT_QUARKS_MIN )
405	{
406	NFCentralQ.Add(particleG);
407	}
408
409	// 2.1a.2********************* visible particles only
410	if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) )
411	{
412	// 2.1a.2.1 Central solenoidal magnetic field
413	TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo
414	// 2.1a.2.2 Filling the calorimetric towers -- includes also forward detectors ?
415	// first checks if the charged particles reach the calo!
416	if( DET->FLAG_bfield \|\|
417	particle->Charge==0 \|\|
418	(!DET->FLAG_bfield && particle->Charge!=0 && particle->PT > DET->TRACK_ptmin))
419	if(
420	(particle->EtaCalo > list_of_calorimeters.getEtamin() ) &&
421	(particle->EtaCalo < list_of_calorimeters.getEtamax() )
422	)
423	{
424	float iEta=UNDEFINED, iPhi=UNDEFINED;
425	DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
426	if (iEta != UNDEFINED && iPhi != UNDEFINED)
427	{
428	D_CaloTower tower(iEta,iPhi); // new tower
429	tower.Set_Eem_Ehad_E_ET(particle->Eparticle->getFem() , particle->Eparticle->getFhad() );
430	list_of_active_towers.addTower(tower);
431	// this list may contain several times the same calotower, as several particles
432	// may leave some energy in the same calotower
433	// After the loop on particles, identical cells in the list should be merged
434	} // iEta and iPhi must be defined
435	}
436
437	// 2.1a.2.3 charged particles in tracker: energy flow
438	// if bfield not simulated, pt should be high enough to be taken into account
439	// it is supposed here that DET->MAX_calo > DET->CEN_max_tracker > DET->CEN_max_mu > 0
440	if( particle->Charge !=0 &&
441	fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available
442	( DET->FLAG_bfield \|\|
443	(!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin)
444	)
445	)
446	{
447	// 2.1a.2.3.1 Filling the particle properties + smearing
448	// Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction
449	// This is the EnergyFlow hypothesis
450	particle->SetEtaPhi(particle->Eta,particle->Phi);
451	float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks
452
453	// 2.1a.2.3.2 Muons
454	if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu)
455	{
456	TLorentzVector p;
457	float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT );
458	if (sPT > 0 && sPT > DET->PTCUT_muon)
459	{
460	p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta));
461	muon.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
462	}
463	sET = (sPT >0)? sPT : 0;
464	}
465	// 2.1a.2.3.3 Electrons
466	else if (pid == pE)
467	{
468	// Finds in which calorimeter the particle has gone, to know its resolution
469
470	D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
471	if(currentCalo.getName() == dummyCalo.getName())
472	{
473	cout << " Warning: the calo coverage behind the tracker is not complete! " << endl;
474	}
475
476	// final smeared EM energy // electromagnetic fraction F_em =1 for electrons;
477	float sE = currentCalo.getElectromagneticResolution().Smear(particle->E);
478	if (sE>0)
479	{
480	sET = sE/cosh(particle->Eta);
481	// NB: ET is found via the calorimetry and not via the track curvature
482
483	TLorentzVector p;
484	p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE);
485	if (sET > DET->PTCUT_elec)
486	electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
487	//if(DET->JET_Eflow) input_particles.push_back(fastjet::PseudoJet(p.Px(),p.Py(),p.Pz(),p.E()));
488	}
489	else { sET=0;} // if negative smeared energy -- needed for the tracks
490	}
491	// 2.1a.2.3.4 Other charged particles : smear them for the tracks!
492	else
493	{ //other particles
494	D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
495	float sEem = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem());
496	float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad());
497	float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 );
498	sET = sE/cosh(particle->EtaCalo);
499	}
500
501	// 2.1a.2.3.5 Tracks
502	if( (rand()%100) < DET->TRACK_eff && sET!=0)
503	{
504	elementTrack = (TRootTracks*) branchTrack->NewEntry();
505	elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
506	TrackCentral.push_back(*elementTrack); // tracks at vertex!
507	if(DET->JET_Eflow)
508	input_particles.push_back(fastjet::PseudoJet(particle->Px,particle->Py,particle->Pz,particle->E));
509	// TODO!!! apply a smearing on the position of the origin of the track
510	// TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout);
511	}
512	} // 2.1a.2.3 : if tracker/energy-flow
513	// 2.1a.2.4 Photons
514	// stays in the tracker -> track available -> gamma ID
515	else if( (pid == pGAMMA) && fabs(particle->EtaCalo)< DET->CEN_max_tracker )
516	{
517	float iEta=UNDEFINED, iPhi=UNDEFINED;
518	DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
519	D_CaloTower tower(iEta,iPhi);
520	// stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search
521	list_of_towers_with_photon.addTower(tower);
522	}
523	// 2.1a.2.5 Neutrals within tracker -- for jet energy flow
524	else if( particle->Charge ==0 && fabs(particle->EtaCalo)< etamax_calocoverage_behindtracker)
525	{
526	float iEta=UNDEFINED, iPhi=UNDEFINED;
527	DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
528	D_CaloTower tower(iEta,iPhi);
529	list_of_centowers_with_neutrals.addTower(tower);
530	}
531	// 2.1a.2.6 : very forward detectors
532	else
533	{
534	if (DET->FLAG_RP==1)
535	{
536	// for the moment, only protons are transported
537	// BUT !!! could be a beam of other particles! (heavy ions?)
538	// BUT ALSO !!! if very forward muons, or others!
539	VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
540	}
541	// 2.1a.2.6: Zero degree calorimeter
542	if(DET->FLAG_vfd==1)
543	{
544	VFD->ZDC(treeWriter,branchZDC,particle);
545	}
546	}
547
548	} // 2.1a.2 : if visible particle
549	delete particle;
550	}// loop on all particles 2.1a
551
552	// 2.1b loop on all (activated) towers
553	// at this stage, list_of_active_towers may contain several times the same tower
554	// first step is to merge identical towers, by matching their (iEta,iPhi)
555
556	list_of_active_towers.sortElements();
557	list_of_active_towers.mergeDuplicates();
558
559	// Calotower smearing
560	list_of_active_towers.smearTowers(list_of_calorimeters);
561
562	for(unsigned int i=0; i<list_of_active_towers.size(); i++)
563	{
564	float iEta = list_of_active_towers[i].getEta();
565	float iPhi = list_of_active_towers[i].getPhi();
566	float e = list_of_active_towers[i].getE();
567	if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0)
568	{
569	elementCalo = (TRootCalo*) branchCalo->NewEntry();
570	elementCalo->set(list_of_active_towers[i]);
571	// not beautiful : should be improved!
572	TLorentzVector p;
573	p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e );
574	PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E()));
575	towers.push_back(Tower);
576	}
577	} // loop on towers
578
579	// 2.1c photon ID
580	// list_of_towers_with_photon is the list of towers with photon candidates
581	// already smeared !
582	// sorts the vector and smears duplicates
583	list_of_towers_with_photon.mergeDuplicates();
584	for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) {
585	float eta = list_of_towers_with_photon[i].getEta();
586	float phi = list_of_towers_with_photon[i].getPhi();
587	D_CaloTower cal(list_of_active_towers.getElement(eta,phi)); //// <---------- buh???????
588	if(cal.getEta() != UNDEFINED && cal.getPhi() != UNDEFINED && cal.getE() > 0)
589	{
590	TLorentzVector p;
591	p.SetPtEtaPhiE(cal.getET(), eta,phi,cal.getE() );
592	if (cal.getET() > DET->PTCUT_gamma) { gamma.push_back(D_Particle(p,pGAMMA,p.Eta(),p.Phi())); }
593	}
594	} // for -- list of photons
595
596	// 2.1d jet-E-flow -- taking into account the neutrals within tracker
597	if(DET->JET_Eflow) {
598	list_of_centowers_with_neutrals.mergeDuplicates();
599	for(unsigned int i=0; i<list_of_centowers_with_neutrals.size(); i++) {
600	float eta = list_of_centowers_with_neutrals[i].getEta();
601	float phi = list_of_centowers_with_neutrals[i].getPhi();
602	D_CaloTower cal(list_of_active_towers.getElement(eta,phi));
603	if(cal.getEta() != UNDEFINED && cal.getPhi() != UNDEFINED && cal.getE() > 0)
604	{
605	TLorentzVector p;
606	p.SetPtEtaPhiE(cal.getET(), eta,phi,cal.getE() );
607	//cout << "**************list: " << p.Px() << " " << p.Py() << " " << p.Pz() << " " << p.E() << endl;
608	input_particles.push_back(fastjet::PseudoJet(p.Px(),p.Py(),p.Pz(),p.E()));
609	}
610	} // for - list_of_centowers
611	} // JET_Eflow
612
613	// 2.2 ******** Output preparation & complex objects *********
614	// 2.2.1 ********************* sorting collections by decreasing pt
615	DET->SortedVector(electron);
616	float iPhiEl=0,iEtaEl=0,ptisoEl=0;
617	for(unsigned int i=0; i < electron.size(); i++)
618	{
619	elementElec = (TRootElectron*) branchElectron->NewEntry();
620	elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E());
621	elementElec->EtaCalo = electron[i].EtaCalo();
622	elementElec->PhiCalo = electron[i].PhiCalo();
623	elementElec->Charge = sign(electron[i].PID());
624	elementElec->IsolFlag = DET->Isolation(electron[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoEl);
625	elementElec->IsolPt = ptisoEl;
626	DET->BinEtaPhi(elementElec->PhiCalo,elementElec->EtaCalo,iPhiEl,iEtaEl);
627	D_CaloTower calElec(list_of_active_towers.getElement(iEtaEl,iPhiEl));
628	elementElec->EHoverEE = calElec.getEhad()/calElec.getEem();
629	}
630
631	DET->SortedVector(muon);
632	float iPhiMu=0,iEtaMu=0,ptisoMu=0;
633	for(unsigned int i=0; i < muon.size(); i++)
634	{
635	elementMu = (TRootMuon*) branchMuon->NewEntry();
636	elementMu->Charge = sign(muon[i].PID());
637	elementMu->Set(muon[i].Px(),muon[i].Py(),muon[i].Pz(),muon[i].E());
638	elementMu->EtaCalo = muon[i].EtaCalo();
639	elementMu->PhiCalo = muon[i].PhiCalo();
640	elementMu->IsolFlag = DET->Isolation(muon[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoMu);
641	elementMu->IsolPt = ptisoMu;
642	DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
643	D_CaloTower calMuon(list_of_active_towers.getElement(iEtaMu,iPhiMu));
644	if( calMuon.getEem() !=0 ) elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
645	else elementMu->EHoverEE = UNDEFINED;
646	elementMu->EtRatio = DET->CaloIsolation(muon[i], list_of_active_towers,iPhiMu,iEtaMu);
647	}
648
649	DET->SortedVector(gamma);
650	for(unsigned int i=0; i < gamma.size(); i++)
651	{
652	elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
653	elementPhoton->Set(gamma[i].Px(),gamma[i].Py(),gamma[i].Pz(),gamma[i].E());
654	D_CaloTower calGamma(list_of_active_towers.getElement(gamma[i].EtaCalo(),gamma[i].PhiCalo()));
655	elementPhoton->EHoverEE = calGamma.getEhad()/calGamma.getEem();
656	}
657
658	// 2.2.2 ************* computes the Missing Transverse Momentum
659	TLorentzVector Att(0.,0.,0.,0.);
660	for(unsigned int i=0; i < towers.size(); i++)
661	{
662	Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
663	if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
664	{
665	PTmis = PTmis + Att;
666	// create a fastjet::PseudoJet with these components and put it onto
667	// back of the input_particles vector
668	if(!DET->JET_Eflow)
669	input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
670	else { if(fabs(Att.Eta()) > DET->CEN_max_tracker)
671	input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
672	}
673	}
674	}
675	elementEtmis = (TRootETmis*) branchETmis->NewEntry();
676	elementEtmis->ET = (PTmis).Pt();
677	elementEtmis->Phi = (-PTmis).Phi();
678	elementEtmis->Px = (-PTmis).Px();
679	elementEtmis->Py = (-PTmis).Py();
680
681	// 2.2.3 ************* jets, B-tag, tau jets
682	vector<int> NTrackJet; //for number of tracks
683	vector<float> EHADEEM; //for energyHad over energyEm
684	sorted_jets=JETRUN->RunJets(input_particles, TrackCentral,NTrackJet,EHADEEM,list_of_active_towers);
685	JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ,NTrackJet,EHADEEM);
686	JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral,NTrackJet,EHADEEM);
687
688	treeWriter->Fill();
689	} // 2. Loop over all events ('for' loop)
690
691	cout <<" Exiting detector simulation... "<< endl;
692
693
694	treeWriter->Write();
695	delete treeWriter;
696	loopwatch.Stop();
697
698
699
700	// 3. ******** Trigger & Frog *********
701	// 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard
702	triggerwatch.Start();
703	if(DET->FLAG_trigger == 1)
704	{
705	cout <<" "<<endl;
706	cout <<" ########### Start Trigger selection ########### "<< endl;
707
708	// input
709	TChain chainT("Analysis");
710	chainT.Add(outputfilename.c_str());
711	ExRootTreeReader *treeReaderT = new ExRootTreeReader(&chainT);
712
713	// output
714	TClonesArray *branchElecTrig = treeReaderT->UseBranch("Electron");
715	TClonesArray *branchMuonTrig = treeReaderT->UseBranch("Muon");
716	TClonesArray *branchJetTrig = treeReaderT->UseBranch("Jet");
717	TClonesArray *branchTauJetTrig = treeReaderT->UseBranch("TauJet");
718	TClonesArray *branchPhotonTrig = treeReaderT->UseBranch("Photon");
719	TClonesArray *branchETmisTrig = treeReaderT->UseBranch("ETmis");
720
721	ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger");
722	ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class());
723
724
725	Long64_t entryT, allEntriesT = treeReaderT->GetEntries();
726	// loop on all entries
727	for(entryT = 0; entryT < allEntriesT; ++entryT) {
728	treeWriterT->Clear();
729	treeReaderT->ReadEntry(entryT);
730	TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
731	treeWriterT->Fill();
732	} // loop on all entries
733	cout <<" Exiting trigger simulation... "<< endl;
734
735	treeWriterT->Write();
736	delete treeWriterT;
737	delete treeReaderT;
738	} // trigger
739	triggerwatch.Stop();
740
741
742	// 3.2 ************** FROG display
743	frogwatch.Start();
744	if(DET->FLAG_frog == 1) {
745	cout <<" "<<endl;
746	cout <<" ################## Start FROG ################# "<< endl;
747
748	FrogDisplay *FROG = new FrogDisplay(DET);
749	FROG->BuildEvents(outputfilename);
750	FROG->BuildGeom();
751	delete FROG;
752	cout <<" Exiting FROG preparation... "<< endl;
753	}
754	frogwatch.Stop();
755
756	// 3.3 *************** LHCO output
757	lhcowatch.Start();
758	if(DET->FLAG_lhco == 1){
759	cout <<" "<<endl;
760	cout <<" ############ Start LHCO conversion ############ "<< endl;
761
762	//LHCOConverter *LHCO = new LHCOConverter(outputfilename,LogNameLHCO);
763	LHCOConverter *LHCO = new LHCOConverter(outputfilename,"");
764	LHCO->CopyRunLogFile();
765	LHCO->ConvertExRootAnalysisToLHCO();
766	delete LHCO;
767	cout <<" Exiting LHCO conversion... "<< endl;
768	}
769	lhcowatch.Stop();
770
771
772
773	// 4. ******** End & Exit *********
774
775	globalwatch.Stop();
776	time_report(globalwatch,loopwatch,triggerwatch,frogwatch,lhcowatch,DET->FLAG_frog,DET->FLAG_trigger,DET->FLAG_lhco,LogName,allEntries);
777
778	cout <<" "<< endl;
779	cout <<" Exiting Delphes ... "<< endl;
780	cout <<" "<< endl;
781	cout <<"*********************************************************************"<< endl;
782	cout <<"*********************************************************************"<< endl;
783
784	delete treeReader;
785	delete DET;
786	delete TRIGT;
787	delete TRACP;
788	delete JETRUN;
789	delete VFD;
790
791	// todo("TODO");
792	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: