Fork me on GitHub

Context Navigation

source: svn/trunk/Resolutions.cpp@ 467

Last change on this file since 467 was 467, checked in by Xavier Rouby, 15 years ago
cleaning, comments, progressbar, cout
File size: 18.2 KB

Rev	Line
[443]	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	***********************************************************************/
[9]	33
[445]	34	/// \file Resolution.cpp
[447]	35	/// \brief Resolution for CMS
[264]	36
[9]	37	#include "TChain.h"
	38	#include "TApplication.h"
[227]	39	#include "TFile.h"
[467]	40	#include "TStopwatch.h"
[9]	41
[227]	42	#include "ExRootTreeReader.h"
	43	#include "ExRootTreeWriter.h"
	44	#include "ExRootTreeBranch.h"
[467]	45	#include "ExRootProgressBar.h"
[227]	46	#include "TreeClasses.h"
[9]	47
[227]	48	#include "SmearUtil.h"
	49	#include "JetsUtil.h"
	50	#include "BFieldProp.h"
[19]	51
	52	#include<vector>
	53	#include<iostream>
	54
[9]	55	using namespace std;
	56
	57	//------------------------------------------------------------------------------
	58
[39]	59	// //******************************** PYTHIA INFORMATION*******************************
[26]	60
[39]	61	TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
	62	{
	63	TIter it((TCollection*)GEN);
	64	it.Reset();
	65	TRootGenParticle *gen1;
	66	TSimpleArray<TRootGenParticle> array,array2;
[26]	67
[39]	68	while((gen1 = (TRootGenParticle*) it.Next()))
	69	{
	70	array.Add(gen1);
	71	}
	72	it.Reset();
	73	bool tauhad;
	74	while((gen1 = (TRootGenParticle*) it.Next()))
	75	{
	76	tauhad=false;
	77	if(abs(gen1->PID)==15)
	78	{
	79	int d1=gen1->D1;
	80	int d2=gen1->D2;
	81	if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
	82	{
	83	tauhad=true;
	84	for(int d=d1; d < d2+1; d++)
[264]	85	{
[39]	86	if(abs(array[d]->PID)== pE \|\| abs(array[d]->PID)== pMU)tauhad=false;
	87	}
	88	}
	89	}
	90	if(tauhad)array2.Add(gen1);
	91	}
	92	return array2;
	93	}
	94
[258]	95	double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
	96	double Energie=0;
	97	for(unsigned int i=0; i < towers.size(); i++) {
	98	if(towers[i].Pt() < JET_seed) continue;
	99	if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
	100	Energie += towers[i].E();
	101	}
	102	}
	103	return Energie;
	104	}
[39]	105
	106
[447]	107	void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
[9]	108	{
	109	JETSm.SetPxPyPzE(0,0,0,0);
	110	float deltaRtest=5000;
[258]	111	TIter itJet((TCollection*)branchJet);
	112	TRootJet *jet;
	113	itJet.Reset();
	114	while( (jet = (TRootJet*) itJet.Next()) )
	115	{
[19]	116	TLorentzVector Att;
[264]	117	Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E);
[19]	118	if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
	119	{
	120	deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
[447]	121	if(deltaRtest < 0.25)
[19]	122	{
	123	JETSm = Att;
	124	}
	125	}
	126	}
[9]	127	}
	128
[258]	129	void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
	130	{
	131	ELECSm.SetPxPyPzE(0,0,0,0);
	132	float deltaRtest=5000;
	133	TIter itElec((TCollection*)branchElec);
	134	TRootElectron *elec;
	135	itElec.Reset();
	136	while( (elec = (TRootElectron*) itElec.Next()) )
	137	{
	138	TLorentzVector Att;
[264]	139	Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E);
[258]	140	if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
	141	{
	142	deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
	143	if(deltaRtest < 0.025)
	144	{
	145	ELECSm = Att;
	146	}
	147	}
	148	}
	149	}
	150
	151	void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
	152	{
	153	MUONSm.SetPxPyPzE(0,0,0,0);
	154	float deltaRtest=5000;
	155	TIter itMuon((TCollection*)branchMuon);
	156	TRootMuon *muon;
	157	itMuon.Reset();
	158	while( (muon = (TRootMuon*) itMuon.Next()) )
	159	{
	160	TLorentzVector Att;
	161	Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
	162	if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
	163	{
	164	deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
	165	if(deltaRtest < 0.025)
	166	{
	167	MUONSm = Att;
	168	}
	169	}
	170	}
	171	}
	172
	173	unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
	174	unsigned int numtrack=0;
	175	TIter itTrack((TCollection*)branchTracks);
	176	TRootTracks *track;
	177	itTrack.Reset();
	178	while( (track = (TRootTracks*) itTrack.Next()) )
	179	{
	180	if((track->PT < pt_track )\|\|
	181	(DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
	182	)continue;
	183	numtrack++;
	184	}
	185	return numtrack;
	186	}
	187
	188
	189
[9]	190	int main(int argc, char *argv[])
	191	{
	192	int appargc = 2;
[467]	193	char *appName= new char[20];
	194	char *appOpt= new char[20];
	195	sprintf(appName,"Resolution_ATLAS");
	196	sprintf(appOpt,"-b");
	197	char *appargv[] = {appName,appOpt};
[9]	198	TApplication app(appName, &appargc, appargv);
[467]	199	delete [] appName;
	200	delete [] appOpt;
	201
[258]	202	if(argc != 3) {
	203	cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
[467]	204	cout << " input_file - input file in Delphes-root format," << endl;
[9]	205	cout << " output_file - output file." << endl;
	206	exit(1);
	207	}
	208
[467]	209
	210	// 1. ******** initialisation *********
	211
	212	srand (time (NULL));
	213	TStopwatch globalwatch;
	214	globalwatch.Start();
[9]	215
[258]	216	string inputfilename(argv[1]), outputfilename(argv[2]);
	217
[467]	218	// 1.1 checks the input file
	219	if(inputfilename.find(".root") > inputfilename.length() ) {
	220	cout << "input_file should be a .root file from Delphes!\n";
	221	return -1;
	222	}
	223	TFile f(inputfilename.c_str());
	224	if (!f.FindKey("GEN") \|\| !f.FindKey("Analysis") ) {
	225	cout << "input_file should be a .root file from Delphes!\n";
	226	cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
	227	return -1;
	228	}
	229	f.Close();
	230
	231	// 1.2 checks the output file
[9]	232	if(outputfilename.find(".root") > outputfilename.length() ) {
	233	cout << "output_file should be a .root file!\n";
	234	return -1;
	235	}
[447]	236	TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
[9]	237	outputFile->Close();
[467]	238
	239	// 1.3 Reads the trees in input file
[258]	240	TChain chainGEN("GEN");
	241	chainGEN.Add(inputfilename.c_str());
	242	ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
	243	TChain chain("Analysis");
	244	chain.Add(inputfilename.c_str());
[9]	245	ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
[258]	246	const TClonesArray *branchJet = treeReader->UseBranch("Jet");
	247	const TClonesArray *branchElec = treeReader->UseBranch("Electron");
	248	const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
	249	const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
	250	const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
	251	const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
[9]	252	TIter itGen((TCollection*)branchGen);
	253
[467]	254
	255	// 1.4 Prepares the output root file
[9]	256	ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
	257	ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
	258	ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
	259	ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
	260	ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
	261	ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
[467]	262
	263	// 1.5 other initialisations
[9]	264	TRootGenParticle *particle;
[39]	265	RESOLELEC * elementElec;
[9]	266	RESOLMUON *elementMuon;
	267	RESOLJET *elementJet;
	268	TAUHAD *elementTaujet;
	269	ETMIS *elementEtmis;
[26]	270
[467]	271	int numTau=0, numTauRec=0;
[26]	272
[258]	273	RESOLution *DET = new RESOLution();
[467]	274	/*
	275	string detectorcard = "data/DetectorCard_CMS.dat";
	276	const float dR_jetpairing = 0.25;
	277	const float jet_pt_cut = 1;
	278	*/
[447]	279	DET->ReadDataCard("data/DetectorCard_CMS.dat");
[88]	280
[258]	281
[447]	282	JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
[88]	283
[467]	284	TLorentzVector genMomentum(0,0,0,0); //generator level information
	285	TLorentzVector recoMomentum(0,0,0,0);//reconstructed level information
[9]	286	LorentzVector jetMomentum;
[26]	287
[88]	288	vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
	289	vector<fastjet::PseudoJet> sorted_jetsGEN;
[307]	290	vector<int> NTrackJet;
[258]	291	vector<TLorentzVector> towers;
[152]	292
[467]	293	// 2. Loop over all events
[9]	294	Long64_t entry, allEntries = treeReader->GetEntries();
[467]	295	cout << endl << endl;
	296	cout <<"*********************************************************************"<< endl;
	297	cout <<" "<< endl;
	298	cout <<" "<< endl;
	299	cout <<" "<< endl;
	300	cout <<" ####### Start resolution processing ######## "<< endl;
	301	cout << left << setw(52) <<"** Total number of events to run: "<<""
	302	<< left << setw(15) << allEntries <<""
	303	<< right << setw(2) <<"**"<<endl;
	304
	305	ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
	306
[9]	307	for(entry = 0; entry < allEntries; ++entry)
	308	{
[467]	309	Progress->Update(entry);
[88]	310	TLorentzVector PTmisReco(0,0,0,0);
	311	TLorentzVector PTmisGEN(0,0,0,0);
[9]	312	treeReader->ReadEntry(entry);
[258]	313	treeReaderGEN->ReadEntry(entry);
[9]	314	treeWriter->Clear();
	315
	316	TSimpleArray<TRootGenParticle> bGen;
	317	itGen.Reset();
	318	TSimpleArray<TRootGenParticle> NFCentralQ;
[89]	319
[88]	320	input_particlesGEN.clear();
[23]	321	towers.clear();
[26]	322
[9]	323	// Loop over all particles in event
	324	while( (particle = (TRootGenParticle*) itGen.Next()) )
	325	{
	326	genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
[26]	327
[9]	328	int pid = abs(particle->PID);
	329	float eta = fabs(particle->Eta);
[89]	330
[88]	331	//input generator level particle for jet algorithm
[94]	332	if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
[26]	333	{
[88]	334	input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
[26]	335	}
[89]	336
[88]	337	//Calculate ETMIS from generated particles
	338	if((pid == pNU1) \|\| (pid == pNU2) \|\| (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
[467]	339
	340	//Electrons and muons
[258]	341	if( (particle->Status == 1) &&
	342	((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
	343	(fabs(particle->Eta) < DET->CEN_max_calo_fwd)
	344	)
[88]	345	{
[258]	346	eta=fabs(genMomentum.Eta());
[152]	347
[258]	348	switch(pid) {
	349
	350	case pE: // all electrons with eta < DET->MAX_CALO_FWD
	351	PairingElec(recoMomentum,genMomentum,branchElec);
[264]	352	if(recoMomentum.E()!=0){
[258]	353	elementElec=(RESOLELEC*) branchelec->NewEntry();
	354	elementElec->E = genMomentum.E();
	355	elementElec->SmearedE = recoMomentum.E();}
	356	break; // case pE
	357	case pMU: // all muons with eta < DET->MAX_MU
	358	PairingMuon(recoMomentum,genMomentum,branchMuon);
[467]	359	if(recoMomentum.E()!=0){
[258]	360	elementMuon = (RESOLMUON*) branchmuon->NewEntry();
[467]	361	elementMuon->OverPT = 1./genMomentum.Pt();
	362	elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
[258]	363	break; // case pMU
	364	default:
	365	break;
	366	} // switch (pid)
	367	}
	368
[9]	369	} // while
[264]	370
[89]	371	//compute missing transverse energy from calo towers
[258]	372	TIter itCalo((TCollection*)branchTowers);
	373	TRootCalo *calo;
	374	itCalo.Reset();
[89]	375	TLorentzVector Att(0.,0.,0.,0.);
	376	float ScalarEt=0;
[258]	377	while( (calo = (TRootCalo*) itCalo.Next()) )
	378	{
[264]	379	if(calo->E !=0){
	380	Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E);
	381	towers.push_back(Att);
	382	if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
	383	{
	384	ScalarEt = ScalarEt + calo->getET();
	385	PTmisReco = PTmisReco + Att;
	386	}
	387	}
	388	}
[89]	389	elementEtmis= (ETMIS*) branchetmis->NewEntry();
	390	elementEtmis->Et = (PTmisGEN).Pt();
	391	elementEtmis->Ex = (-PTmisGEN).Px();
	392	elementEtmis->SEt = ScalarEt;
	393	elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
	394	elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
[9]	395
	396	//*****************************
[310]	397	sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN);
[88]	398
[39]	399	TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
[88]	400	TLorentzVector JETreco(0,0,0,0);
	401	for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
	402	TLorentzVector JETgen(0,0,0,0);
	403	JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
[258]	404	PairingJet(JETreco,JETgen,branchJet);
[447]	405	if(JETreco.Pt()>1)
[19]	406	{
	407	elementJet= (RESOLJET*) branchjet->NewEntry();
[88]	408	elementJet->PT = JETgen.Et();
[467]	409	elementJet->SmearedPT = JETreco.Et()/JETgen.Et(); //// la difference pourrait Ãªtre ici
[19]	410	}
[39]	411	}
[467]	412	numTau += TausHadr.GetEntries();
[258]	413
	414	TIter itJet((TCollection*)branchJet);
	415	TRootJet *jet;
	416	itJet.Reset();
	417	while( (jet = (TRootJet*) itJet.Next()) )
	418	{
	419	TLorentzVector JETT(0,0,0,0);
	420	JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
	421	if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
[39]	422	{
	423	for(Int_t i=0; i<TausHadr.GetEntries();i++)
	424	{
	425	if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
	426	{
	427	elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
[467]	428	elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
[258]	429	elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
	430	if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95
	431	&& (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
[39]	432	}
	433	}
	434	}
	435
[26]	436
	437	} // for itJet : loop on all jets
	438
[9]	439	treeWriter->Fill();
	440	} // Loop over all events
	441	treeWriter->Write();
[467]	442	globalwatch.Stop();
[152]	443
[447]	444
[467]	445	// Screen output
	446	cout <<" "<< endl;
	447	cout <<" "<< endl;
	448	cout <<" ################## Time report ################# "<< endl;
	449	cout << left << setw(32) <<"** Time report for "<<""
	450	<< left << setw(15) << allEntries <<""
	451	<< right << setw(22) <<"events **"<<endl;
	452	cout <<" "<< endl;
	453	cout << left << setw(10) <<"**"<<""
	454	<< left << setw(15) <<"Time (s):"<<""
	455	<< right << setw(15) <<"CPU"<<""
	456	<< right << setw(15) <<"Real"<<""
	457	<< right << setw(14) <<"**"<<endl;
	458	cout << left << setw(10) <<"**"<<""
	459	<< left << setw(15) <<" + Global:"<<""
	460	<< right << setw(15) <<globalwatch.CpuTime()<<""
	461	<< right << setw(15) <<globalwatch.RealTime()<<""
	462	<< right << setw(14) <<"**"<<endl;
	463
	464
	465	string tempstring = detectorcard + " has been used.";
	466	cout <<" ################## Configuration ############### "<< endl;
	467	cout << left << setw(16) << "** " << ""
	468	<< left << setw(51) << tempstring << "**" << endl;
	469
	470	tempstring = outputfilename + " has been created.";
	471	cout << left << setw(16) << "** " << ""
	472	<< left << setw(51) << tempstring << "**" << endl;
	473
	474	cout << left << setw(16) << "** " << ""
	475	<< left << setw(51) << get_time_date() << "**" << endl;
	476
	477
	478	cout <<" "<< endl;
	479	cout <<" Exiting ... "<< endl;
	480	cout <<" "<< endl;
	481	cout <<" "<< endl;
	482	cout <<"*********************************************************************"<< endl;
	483
	484
[9]	485
	486	delete treeWriter;
	487	delete treeReader;
	488	delete DET;
	489	}
	490

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

Download in other formats: