Context Navigation

Resolutions.cpp@ 686

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

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

Context Navigation

source: svn/trunk/Resolutions.cpp@ 686

Download in other formats: