Context Navigation

Resolutions_ATLAS.cpp@ 942

Last change on this file since 942 was 514, checked in by Xavier Rouby, 15 years ago
cleaning
File size: 19.2 KB

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

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

Context Navigation

source: svn/trunk/Resolutions_ATLAS.cpp@ 942

Download in other formats: