Fork me on GitHub

Context Navigation

source: svn/trunk/Delphes.cpp@ 250

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

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

Download in other formats: