Fork me on GitHub

Context Navigation

source: svn/trunk/Delphes.cpp@ 254

Last change on this file since 254 was 251, checked in by severine ovyn, 16 years ago
add progress bar
File size: 28.6 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"
[228]	17	#include "TFile.h"
[2]	18
[228]	19	#include "ExRootTreeReader.h"
	20	#include "ExRootTreeWriter.h"
	21	#include "ExRootTreeBranch.h"
[2]	22
[228]	23	#include "DataConverter.h"
	24	#include "HEPEVTConverter.h"
	25	#include "LHEFConverter.h"
	26	#include "STDHEPConverter.h"
[2]	27
[228]	28	#include "SmearUtil.h"
	29	#include "BFieldProp.h"
	30	#include "TriggerUtil.h"
	31	#include "VeryForward.h"
	32	#include "JetsUtil.h"
	33	#include "FrogUtil.h"
[2]	34
[55]	35	#include <vector>
	36	#include <iostream>
[11]	37
[251]	38	#include "Utilities/ExRootAnalysis/interface/ExRootProgressBar.h"
	39
[2]	40	using namespace std;
	41
	42	//------------------------------------------------------------------------------
	43	void todo(string filename) {
	44	ifstream infile(filename.c_str());
	45	cout << "** TODO list ..." << endl;
	46	while(infile.good()) {
[94]	47	string temp;
	48	getline(infile,temp);
	49	cout << "*" << temp << endl;
[2]	50	}
	51	cout << "** done...\n";
	52	}
	53
	54	//------------------------------------------------------------------------------
	55
	56	int main(int argc, char *argv[])
	57	{
	58	int appargc = 2;
[228]	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};
[2]	64	TApplication app(appName, &appargc, appargv);
[228]	65	delete [] appName;
	66	delete [] appOpt;
	67
	68
[249]	69	if(argc != 3 && argc != 4 && argc != 5) {
[94]	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);
[2]	76	}
[228]	77
	78
	79
[212]	80
	81	cout << endl << endl;
	82
	83	cout <<"*********************************************************************"<< endl;
	84	cout <<"*********************************************************************"<< endl;
	85	cout <<" "<< endl;
	86	cout <<" Welcome to "<< endl;
	87	cout <<" "<< endl;
	88	cout <<" "<< endl;
	89	cout <<" .ddddddd- lL hH "<< endl;
	90	cout <<" -Dd` `dD: Ll hH` "<< endl;
	91	cout <<" dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss "<< endl;
	92	cout <<" -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs "<< endl;
	93	cout <<" dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. "<< endl;
	94	cout <<" .Dd :dd eE. LlL PpppPP Hh Hh eE sSS "<< endl;
	95	cout <<" dddddd:. eee+: lL. pp. hh. hh eee+ sssssS "<< endl;
	96	cout <<" Pp "<< endl;
	97	cout <<" "<< endl;
	98	cout <<" Delphes, a framework for the fast simulation "<< endl;
	99	cout <<" of a generic collider experiment "<< endl;
	100	cout <<" "<< endl;
	101	cout <<" --- Version 1.3beta of Delphes --- "<< endl;
	102	cout <<" Last date of change: 29 January 2009 "<< endl;
	103	cout <<" "<< endl;
	104	cout <<" "<< endl;
	105	cout <<" This package uses: "<< endl;
	106	cout <<" ------------------ "<< endl;
	107	cout <<" FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] "<< endl;
	108	cout <<" Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] "<< endl;
	109	cout <<" FROG: "<< endl;
	110	cout <<" "<< endl;
	111	cout <<"-----------------------------------------------------------------"<< endl;
	112	cout <<" "<< endl;
	113	cout <<" Main authors: "<< endl;
	114	cout <<" ------------- "<< endl;
	115	cout <<" "<< endl;
	116	cout <<" SÃ©verine Ovyn Xavier Rouby "<< endl;
	117	cout <<" severine.ovyn@uclouvain.be xavier.rouby@cern "<< endl;
	118	cout <<" Center for Particle Physics and Phenomenology (CP3) "<< endl;
	119	cout <<" Universite Catholique de Louvain (UCL) "<< endl;
	120	cout <<" Louvain-la-Neuve, Belgium "<< endl;
	121	cout <<" "<< endl;
	122	cout <<"-----------------------------------------------------------------"<< endl;
	123	cout <<" "<< endl;
	124	cout <<" Former Delphes versions and documentation can be found on : "<< endl;
	125	cout <<" http://www.fynu.ucl.ac.be/delphes.html "<< endl;
	126	cout <<" "<< endl;
	127	cout <<" "<< endl;
	128	cout <<" Disclaimer: this program is a beta version of Delphes and "<< endl;
	129	cout <<" therefore comes without guarantees. Beware of errors and please "<< endl;
	130	cout <<" give us your feedbacks about potential bugs "<< endl;
	131	cout <<" "<< endl;
	132	cout <<"*********************************************************************"<< endl;
	133	cout <<"*********************************************************************"<< endl;
	134
[178]	135	// 1. ******** initialisation *********
	136
[2]	137	srand (time (NULL)); /* Initialisation du gÃ©nÃ©rateur */
[191]	138	TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch;
	139	globalwatch.Start();
	140
[2]	141
[249]	142	//read the output TROOT file
[2]	143	string inputFileList(argv[1]), outputfilename(argv[2]);
[249]	144	if(outputfilename.find(".root") > outputfilename.length()) {
	145	cout <<" ERROR: 'output_file' should be a .root file. Exiting... "<< endl;
[94]	146	exit(1);
[2]	147	}
[44]	148	//create output log-file name
[45]	149	string forLog = outputfilename;
	150	string LogName = forLog.erase(forLog.find(".root"));
[44]	151	LogName = LogName+"_run.log";
[94]	152
[2]	153	TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
	154	outputFile->Close();
[94]	155
[2]	156	string line;
	157	ifstream infile(inputFileList.c_str());
	158	infile >> line; // the first line determines the type of input files
[94]	159
[44]	160	//read the datacard input file
[249]	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	}
[94]	188
[55]	189	//Smearing information
[44]	190	RESOLution *DET = new RESOLution();
[212]	191	cout <<" "<< endl;
	192	cout <<" ####### Start reading DETECTOR parameters ####### "<< endl;
	193	cout << left << setw(40) <<"** Opening configuration card: "<<""
	194	<< left << setw(20) << DetDatacard <<""
	195	<< right << setw(9) <<" **"<<""<<endl;
[44]	196	DET->ReadDataCard(DetDatacard);
	197	DET->Logfile(LogName);
[212]	198	cout << left << setw(42) <<"** Parameters summarised in file: "<<""
	199	<< left << setw(20) << LogName <<""
	200	<< right << setw(7) <<"**"<<""<<endl;
	201	cout <<" "<< endl;
[228]	202
[55]	203	//Trigger information
[249]	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	}
[72]	210	TriggerTable *TRIGT = new TriggerTable();
[80]	211	TRIGT->TriggerCardReader(TrigDatacard.c_str());
[72]	212	TRIGT->PrintTriggerTable(LogName);
[212]	213	if(DET->FLAG_trigger == 1)
	214	{
	215	cout << left << setw(40) <<"** Opening configuration card: "<<""
	216	<< left << setw(20) << TrigDatacard <<""
	217	<< right << setw(9) <<" **"<<""<<endl;
	218	cout <<" "<< endl;
	219	}
[94]	220
[55]	221	//Propagation of tracks in the B field
[228]	222	TrackPropagation *TRACP = new TrackPropagation(DET);
	223	//TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
[94]	224
[55]	225	//Jet information
[228]	226	JetsUtil *JETRUN = new JetsUtil(DET);
	227	//JetsUtil *JETRUN = new JetsUtil(DetDatacard);
[94]	228
[55]	229	//VFD information
[228]	230	VeryForward * VFD = new VeryForward(DET);
	231	//VeryForward * VFD = new VeryForward(DetDatacard);
[178]	232
	233	// data converters
[228]	234	DataConverter *converter=NULL;
[212]	235	cout <<" "<<endl;
	236	cout <<" ####### Start convertion to TRoot format ######## "<< endl;
	237
[2]	238	if(strstr(line.c_str(),".hep"))
	239	{
[212]	240	cout <<" StdHEP file format detected "<<endl;
	241	cout <<" This can take several minutes "<< endl;
[2]	242	converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
	243	}
	244	else if(strstr(line.c_str(),".lhe"))
	245	{
[212]	246	cout <<" LHEF file format detected "<<endl;
	247	cout <<" This can take several minutes "<< endl;
[2]	248	converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
	249	}
	250	else if(strstr(line.c_str(),".root"))
	251	{
[212]	252	cout <<" h2root file format detected "<<endl;
	253	cout <<" This can take several minutes "<< endl;
[2]	254	converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
	255	}
[212]	256	else {
[246]	257	cerr << left << setw(4) <<"** "<<""
[212]	258	<< left << setw(63) << line.c_str() <<""
	259	<< right << setw(2) <<"**"<<endl;
[246]	260	cerr <<" ERROR: File format not identified -- Exiting... "<< endl;
[212]	261	cout <<" "<< endl;
	262	cout <<"*********************************************************************"<< endl;
	263	return -1;};
	264	cout <<" Exiting conversion... "<< endl;
	265
[2]	266	TChain chain("GEN");
	267	chain.Add(outputfilename.c_str());
	268	ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
	269	const TClonesArray *branchGen = treeReader->UseBranch("Particle");
[249]	270
[2]	271	TIter itGen((TCollection*)branchGen);
	272
[178]	273	//Output file : contents of the analysis object data
[2]	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());
[30]	286
[2]	287	TRootGenParticle *particle;
	288	TRootETmis *elementEtmis;
	289	TRootElectron *elementElec;
	290	TRootMuon *elementMu;
	291	TRootPhoton *elementPhoton;
	292	TRootTracks *elementTracks;
	293	TRootCalo *elementCalo;
	294
[184]	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
[2]	298	LorentzVector jetMomentum;
[94]	299
[55]	300	vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
	301	vector<fastjet::PseudoJet> sorted_jets;
[2]	302	vector<TLorentzVector> TrackCentral;
	303	vector<PhysicsTower> towers;
[73]	304	vector<ParticleUtil> electron;
	305	vector<ParticleUtil> muon;
[74]	306	vector<ParticleUtil> gamma;
[98]	307
[30]	308	TSimpleArray<TRootGenParticle> NFCentralQ;
[100]	309	float iPhi=0,iEta=0;
	310
[178]	311
	312
	313	// 2. ******** Loop over all events *********
[2]	314	Long64_t entry, allEntries = treeReader->GetEntries();
[212]	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
[251]	321	ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
	322
[191]	323	loopwatch.Start();
[178]	324
	325	// loop on all events
	326	for(entry = 0; entry < allEntries; ++entry)
[2]	327	{
[251]	328	Progress->Update(entry);
[2]	329	TLorentzVector PTmis(0,0,0,0);
	330	treeReader->ReadEntry(entry);
	331	treeWriter->Clear();
[228]	332	if((entry % 100) == 0 && entry > 0 ) {
[212]	333	cout << left << setw(52) <<"** Processing element # "<<""
	334	<< left << setw(15) << entry <<""
	335	<< right << setw(2) <<"**"<<endl;
	336	}
[2]	337
[30]	338	electron.clear();
	339	muon.clear();
[74]	340	gamma.clear();
[30]	341	NFCentralQ.Clear();
	342
[2]	343	TrackCentral.clear();
	344	towers.clear();
[11]	345	input_particles.clear();
[30]	346
[178]	347	// 2.1 Loop over all particles in event
[74]	348	itGen.Reset();
[178]	349	while( (particle = (TRootGenParticle*) itGen.Next()) ) {
[55]	350	int pid = abs(particle->PID);
[178]	351
	352
	353	// 2.1.1********************* preparation for the b-tagging
[94]	354	//// This subarray is needed for the B-jet algorithm
[2]	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 ?
[94]	357	fabs(particle->Eta) < DET->CEN_max_tracker &&
[2]	358	particle->Status != 1 &&
	359	particle->PT > DET->PT_QUARKS_MIN ) {
	360	NFCentralQ.Add(particle);
	361	}
[94]	362
[178]	363	// 2.1.2 ********************* central detector: keep only visible particles
[2]	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) &&
[59]	367	((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
[94]	368	(fabs(particle->Eta) < DET->CEN_max_calo_fwd)
	369	)
[59]	370	{
[94]	371	genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
[184]	372	genMomentumBfield = genMomentum;
[178]	373
	374	// 2.1.2.1 ********************* central detector: magnetic field
[184]	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());
[178]	378
	379
	380	// 2.1.2.2 ********************* central detector: smearing (calorimeters, muon chambers)
[94]	381	switch(pid) {
[74]	382
[94]	383	case pE: // all electrons with eta < DET->MAX_CALO_FWD
[184]	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));
[94]	387	}
	388	break; // case pE
	389	case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
[184]	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));
[94]	393	}
	394	break; // case pGAMMA
	395	case pMU: // all muons with eta < DET->MAX_MU
[184]	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));
[94]	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
[184]	403	DET->SmearHadron(genMomentumBfield, 0.7);
[94]	404	break; // case hadron
	405	default: // all other final particles with eta < DET->MAX_CALO_FWD
[184]	406	DET->SmearHadron(genMomentumBfield, 1.0);
[94]	407	break;
[178]	408	} // 2.1.2.2 switch (pid)
	409
	410
	411	// 2.1.2.3 ********************* central detector: calotowers
[94]	412	// all final particles but muons and neutrinos
[100]	413	// for calorimetric towers and missing PT
[94]	414	int charge=Charge(pid);
[184]	415	if(genMomentumBfield.E() !=0 && pid != pMU) {
[178]	416	// in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
[184]	417	if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
[178]	418	else { // particles reach calos
	419	// applies the calo segmentation and returns iEta & iPhi
[184]	420	DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta);
[178]	421	if(iEta != -100 && iPhi != -100) {
[184]	422	momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
[100]	423	elementCalo = (TRootCalo*) branchCalo->NewEntry();
[184]	424	elementCalo->Set(momentumCaloSegmentation);
	425	PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
[100]	426	towers.push_back(Tower);
[178]	427	} // if iEta != -100
	428	} // else : when particles reach the calos
	429	} // 2.1.2.3 calotowers
[94]	430
[178]	431
	432	// 2.1.2.4 ********************* central detector: tracks
[94]	433	// all final charged particles
	434	if(
[184]	435	(genMomentumBfield.E()!=0) &&
	436	(fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) &&
	437	(DET->FLAG_bfield \|\| ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&
[178]	438	// if bfield not simulated, pt should be high enough to be taken into account
[94]	439	((rand()%100) < DET->TRACK_eff) &&
	440	(charge!=0)
	441	)
	442	{
	443	elementTracks = (TRootTracks*) branchTracks->NewEntry();
[184]	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!
[178]	452	} // 2.1.2.4 tracks
[94]	453
[178]	454	} // 2.1.2 central detector
[30]	455
[178]	456	// 2.1.3 ********************* forward detectors: zdc
	457	if(DET->FLAG_vfd==1) {
[100]	458	VFD->ZDC(treeWriter,branchZDC,particle);
	459	VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
[178]	460	}
[11]	461
[178]	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
[74]	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());
[184]	474	elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks
[30]	475	}
[74]	476	DET->SortedVector(muon);
[30]	477	for(unsigned int i=0; i < muon.size(); i++) {
[74]	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);
[30]	482	}
[74]	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	}
[30]	488
[178]	489	// 2.2.2 ************* computes the Missing Transverse Momentum
[71]	490	TLorentzVector Att(0.,0.,0.,0.);
	491	for(unsigned int i=0; i < towers.size(); i++)
	492	{
[107]	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	}
[71]	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();
[74]	507
[178]	508	// 2.2.3 ************* B-tag, tau jets
[55]	509	sorted_jets=JETRUN->RunJets(input_particles);
	510	JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
	511	JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
[74]	512
[72]	513	treeWriter->Fill();
[178]	514	} // 2. Loop over all events ('for' loop)
[212]	515
	516	cout <<" "<< endl;
	517	cout <<" Exiting detector simulation... "<< endl;
	518
	519
[2]	520	treeWriter->Write();
[77]	521	delete treeWriter;
[191]	522	loopwatch.Stop();
[178]	523
	524
	525
[212]	526	// 3. ******** Trigger & Frog *********
[178]	527	// 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard
[191]	528	triggerwatch.Start();
[94]	529	if(DET->FLAG_trigger == 1)
[72]	530	{
[212]	531	cout <<" "<<endl;
	532	cout <<" ########### Start Trigger selection ########### "<< endl;
	533
[178]	534	// input
[72]	535	TChain chainT("Analysis");
	536	chainT.Add(outputfilename.c_str());
	537	ExRootTreeReader *treeReaderT = new ExRootTreeReader(&chainT);
[74]	538
[178]	539	// output
[72]	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");
[74]	546
[72]	547	ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger");
	548	ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class());
[178]	549
	550
[72]	551	Long64_t entryT, allEntriesT = treeReaderT->GetEntries();
[178]	552	// loop on all entries
	553	for(entryT = 0; entryT < allEntriesT; ++entryT) {
[72]	554	treeWriterT->Clear();
	555	treeReaderT->ReadEntry(entryT);
	556	TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
	557	treeWriterT->Fill();
[178]	558	} // loop on all entries
[212]	559	cout <<" Exiting trigger simulation... "<< endl;
[94]	560
[72]	561	treeWriterT->Write();
	562	delete treeWriterT;
[228]	563	delete treeReaderT;
[178]	564	} // trigger
[191]	565	triggerwatch.Stop();
[178]	566
	567
	568	// 3.2 ************** FROG display
[191]	569	frogwatch.Start();
[228]	570	if(DET->FLAG_frog == 1) {
[212]	571	cout <<" "<<endl;
	572	cout <<" ################## Start FROG ################# "<< endl;
	573
[228]	574	FrogDisplay *FROG = new FrogDisplay(DET);
	575	FROG->BuildEvents(outputfilename);
	576	FROG->BuildGeom();
	577	delete FROG;
	578	}
	579	frogwatch.Stop();
[94]	580
[178]	581
	582
	583
	584	// 4. ******** End & Exit *********
[212]	585
[191]	586	globalwatch.Stop();
[212]	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;
[191]	608	if(DET->FLAG_trigger == 1)
[212]	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	}
[191]	616	if(DET->FLAG_frog == 1)
[212]	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;
[191]	623
[212]	624	}
[2]	625
[212]	626	cout <<" "<< endl;
	627	cout <<" Exiting Delphes ... "<< endl;
	628	cout <<" "<< endl;
	629	cout <<"*********************************************************************"<< endl;
	630	cout <<"*********************************************************************"<< endl;
	631
	632
[2]	633	delete treeReader;
	634	delete DET;
[251]	635	// delete TRIGT;
[74]	636	delete TRACP;
	637	delete JETRUN;
	638	delete VFD;
[228]	639	delete converter;
[94]	640
[191]	641	// todo("TODO");
[2]	642	}

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

Download in other formats: