Fork me on GitHub

Context Navigation

source: svn/trunk/Resolutions.cpp@ 445

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

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

Download in other formats: