/*********************************************************************** ** ** ** /----------------------------------------------\ ** ** | Delphes, a framework for the fast simulation | ** ** | of a generic collider experiment | ** ** \------------- arXiv:0903.2225v1 ------------/ ** ** ** ** ** ** This package uses: ** ** ------------------ ** ** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 ** ** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] ** ** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] ** ** FROG: [hep-ex/0901.2718v1] ** ** HepMC: Comput. Phys. Commun.134 (2001) 41 ** ** ** ** ------------------------------------------------------------------ ** ** ** ** Main authors: ** ** ------------- ** ** ** ** Severine Ovyn Xavier Rouby ** ** severine.ovyn@uclouvain.be xavier.rouby@cern ** ** ** ** Center for Particle Physics and Phenomenology (CP3) ** ** Universite catholique de Louvain (UCL) ** ** Louvain-la-Neuve, Belgium ** ** ** ** Copyright (C) 2008-2009, ** ** All rights reserved. ** ** ** ***********************************************************************/ /// \file Resolution.cpp /// \brief Resolution for CMS #include "TChain.h" #include "TApplication.h" #include "TFile.h" #include "ExRootTreeReader.h" #include "ExRootTreeWriter.h" #include "ExRootTreeBranch.h" #include "TreeClasses.h" #include "DataConverter.h" #include "HEPEVTConverter.h" #include "LHEFConverter.h" #include "STDHEPConverter.h" #include "SmearUtil.h" #include "JetsUtil.h" #include "BFieldProp.h" //#include "PseudoJet.hh" //#include "ClusterSequence.hh" #include #include using namespace std; //------------------------------------------------------------------------------ // //********************************** PYTHIA INFORMATION********************************* TSimpleArray TauHadr(const TClonesArray *GEN) { TIter it((TCollection*)GEN); it.Reset(); TRootGenParticle *gen1; TSimpleArray array,array2; while((gen1 = (TRootGenParticle*) it.Next())) { array.Add(gen1); } it.Reset(); bool tauhad; while((gen1 = (TRootGenParticle*) it.Next())) { tauhad=false; if(abs(gen1->PID)==15) { int d1=gen1->D1; int d2=gen1->D2; if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0)) { tauhad=true; for(int d=d1; d < d2+1; d++) { if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false; } } } if(tauhad)array2.Add(gen1); } return array2; } double EnergySmallCone(const vector &towers, const float eta, const float phi,float energy_scone,float JET_seed) { double Energie=0; for(unsigned int i=0; i < towers.size(); i++) { if(towers[i].Pt() < JET_seed) continue; if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) { Energie += towers[i].E(); } } return Energie; } void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet) { JETSm.SetPxPyPzE(0,0,0,0); float deltaRtest=5000; TIter itJet((TCollection*)branchJet); TRootJet *jet; itJet.Reset(); while( (jet = (TRootJet*) itJet.Next()) ) { TLorentzVector Att; Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E); if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest) { deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()); if(deltaRtest < 0.25) { JETSm = Att; } } } } void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec) { ELECSm.SetPxPyPzE(0,0,0,0); float deltaRtest=5000; TIter itElec((TCollection*)branchElec); TRootElectron *elec; itElec.Reset(); while( (elec = (TRootElectron*) itElec.Next()) ) { TLorentzVector Att; Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E); if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest) { deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()); if(deltaRtest < 0.025) { ELECSm = Att; } } } } void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon) { MUONSm.SetPxPyPzE(0,0,0,0); float deltaRtest=5000; TIter itMuon((TCollection*)branchMuon); TRootMuon *muon; itMuon.Reset(); while( (muon = (TRootMuon*) itMuon.Next()) ) { TLorentzVector Att; Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E); if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest) { deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()); if(deltaRtest < 0.025) { MUONSm = Att; } } } } unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) { unsigned int numtrack=0; TIter itTrack((TCollection*)branchTracks); TRootTracks *track; itTrack.Reset(); while( (track = (TRootTracks*) itTrack.Next()) ) { if((track->PT < pt_track )|| (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone) )continue; numtrack++; } return numtrack; } int main(int argc, char *argv[]) { int appargc = 2; char *appName = "Resolution"; char *appargv[] = {appName, "-b"}; TApplication app(appName, &appargc, appargv); if(argc != 3) { cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl; cout << " input_list - list of files in root format," << endl; cout << " output_file - output file." << endl; exit(1); } srand (time (NULL)); /* Initialisation du générateur */ //read the input TROOT file string inputfilename(argv[1]), outputfilename(argv[2]); if(outputfilename.find(".root") > outputfilename.length() ) { cout << "output_file should be a .root file!\n"; return -1; } TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after outputFile->Close(); TChain chainGEN("GEN"); chainGEN.Add(inputfilename.c_str()); ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN); TChain chain("Analysis"); chain.Add(inputfilename.c_str()); ExRootTreeReader *treeReader = new ExRootTreeReader(&chain); const TClonesArray *branchJet = treeReader->UseBranch("Jet"); const TClonesArray *branchElec = treeReader->UseBranch("Electron"); const TClonesArray *branchMuon = treeReader->UseBranch("Muon"); const TClonesArray *branchTracks = treeReader->UseBranch("Tracks"); const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower"); const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle"); TIter itGen((TCollection*)branchGen); //write the output root file ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis"); ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class()); ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class()); ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class()); ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class()); ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class()); TRootGenParticle *particle; RESOLELEC * elementElec; RESOLMUON *elementMuon; RESOLJET *elementJet; TAUHAD *elementTaujet; ETMIS *elementEtmis; int numTau=0; int numTauRec=0; RESOLution *DET = new RESOLution(); DET->ReadDataCard("data/DetectorCard_CMS.dat"); //Jet information JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat"); TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information LorentzVector jetMomentum; vector input_particlesGEN;//for FastJet algorithm vector sorted_jetsGEN; vector NTrackJet; vector towers; // Loop over all events Long64_t entry, allEntries = treeReader->GetEntries(); cout << "** Chain contains " << allEntries << " events" << endl; for(entry = 0; entry < allEntries; ++entry) { TLorentzVector PTmisReco(0,0,0,0); TLorentzVector PTmisGEN(0,0,0,0); treeReader->ReadEntry(entry); treeReaderGEN->ReadEntry(entry); treeWriter->Clear(); if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl; TSimpleArray bGen; itGen.Reset(); TSimpleArray NFCentralQ; input_particlesGEN.clear(); towers.clear(); // Loop over all particles in event while( (particle = (TRootGenParticle*) itGen.Next()) ) { genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E); int pid = abs(particle->PID); float eta = fabs(particle->Eta); //input generator level particle for jet algorithm if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd) { input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); } //Calculate ETMIS from generated particles if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum; if( (particle->Status == 1) && ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) && (fabs(particle->Eta) < DET->CEN_max_calo_fwd) ) { eta=fabs(genMomentum.Eta()); switch(pid) { case pE: // all electrons with eta < DET->MAX_CALO_FWD PairingElec(recoMomentum,genMomentum,branchElec); if(recoMomentum.E()!=0){ elementElec=(RESOLELEC*) branchelec->NewEntry(); elementElec->E = genMomentum.E(); elementElec->SmearedE = recoMomentum.E();} break; // case pE case pMU: // all muons with eta < DET->MAX_MU PairingMuon(recoMomentum,genMomentum,branchMuon); if(recoMomentum.E() !=0){ elementMuon = (RESOLMUON*) branchmuon->NewEntry(); elementMuon->OverPT = (1/genMomentum.Pt()); elementMuon->OverSmearedPT = (1/recoMomentum.Pt());} break; // case pMU default: break; } // switch (pid) } } // while //compute missing transverse energy from calo towers TIter itCalo((TCollection*)branchTowers); TRootCalo *calo; itCalo.Reset(); TLorentzVector Att(0.,0.,0.,0.); float ScalarEt=0; while( (calo = (TRootCalo*) itCalo.Next()) ) { if(calo->E !=0){ Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E); towers.push_back(Att); if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd) { ScalarEt = ScalarEt + calo->getET(); PTmisReco = PTmisReco + Att; } } } elementEtmis= (ETMIS*) branchetmis->NewEntry(); elementEtmis->Et = (PTmisGEN).Pt(); elementEtmis->Ex = (-PTmisGEN).Px(); elementEtmis->SEt = ScalarEt; elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt(); elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px(); //***************************** sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN); TSimpleArray TausHadr = TauHadr(branchGen); TLorentzVector JETreco(0,0,0,0); for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) { TLorentzVector JETgen(0,0,0,0); JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E()); PairingJet(JETreco,JETgen,branchJet); if(JETreco.Pt()>1) { elementJet= (RESOLJET*) branchjet->NewEntry(); elementJet->PT = JETgen.Et(); elementJet->SmearedPT = JETreco.Et()/JETgen.Et(); } } numTau = numTau+TausHadr.GetEntries(); TIter itJet((TCollection*)branchJet); TRootJet *jet; itJet.Reset(); while( (jet = (TRootJet*) itJet.Next()) ) { TLorentzVector JETT(0,0,0,0); JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E); if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) { for(Int_t i=0; iPhi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1) { elementTaujet= (TAUHAD*) branchtaujet->NewEntry(); elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()); elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone); if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 && (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++; } } } } // for itJet : loop on all jets //cout<<"i"<Fill(); } // Loop over all events treeWriter->Write(); float frac = numTauRec/numTau; cout<