Fork me on GitHub

Context Navigation

source: svn/trunk/Delphes.cpp@ 217

Last change on this file since 217 was 212, checked in by severine ovyn, 16 years ago
my improvement of output
File size: 26.8 KB

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

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

Download in other formats: