Fork me on GitHub

Context Navigation

source: svn/trunk/Resolutions.cpp@ 489

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

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

Download in other formats: