/*********************************************************************** ** ** ** /----------------------------------------------\ ** ** | 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. ** ** ** ***********************************************************************/ #include "JetsUtil.h" using namespace std; JetsUtil::JetsUtil() : plugins(NULL) { DET = new RESOLution(); init(); } JetsUtil::JetsUtil(const string & DetDatacard) : plugins(NULL) { DET = new RESOLution(); DET->ReadDataCard(DetDatacard); init(); } JetsUtil::JetsUtil(const RESOLution * DetDatacard) : plugins(NULL) { DET = new RESOLution(*DetDatacard); init(); } JetsUtil::JetsUtil(const JetsUtil& jet){ DET = new RESOLution(*(jet.DET)); } JetsUtil& JetsUtil::operator=(const JetsUtil& jet) { if(this==&jet) return *this; DET = new RESOLution(*(jet.DET)); return *this; } void JetsUtil::init() { if(plugins) delete plugins; switch(DET->JET_jetalgo) { default: case 1: { // set up a CDF midpoint jet definition #ifdef ENABLE_PLUGIN_CDFCONES plugins = new fastjet::CDFJetCluPlugin(DET->JET_seed,DET->JET_coneradius,DET->JET_C_adjacencycut,DET->JET_C_maxiterations,DET->JET_C_iratch,DET->JET_overlap); jet_def = fastjet::JetDefinition(plugins); #else plugins = NULL; #endif } break; case 2: { // set up a CDF midpoint jet definition #ifdef ENABLE_PLUGIN_CDFCONES plugins = new fastjet::CDFMidPointPlugin (DET->JET_seed,DET->JET_coneradius,DET->JET_M_coneareafraction,DET->JET_M_maxpairsize,DET->JET_M_maxiterations,DET->JET_overlap); jet_def = fastjet::JetDefinition(plugins); #else plugins = NULL; #endif } break; case 3: { // set up a siscone jet definition #ifdef ENABLE_PLUGIN_SISCONE plugins = new fastjet::SISConePlugin (DET->JET_coneradius,DET->JET_overlap,DET->JET_S_npass, DET->JET_S_protojet_ptmin); jet_def = fastjet::JetDefinition(plugins); #else plugins = NULL; #endif } break; case 4: { jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->JET_coneradius); } break; case 5: { jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->JET_coneradius); } break; case 6: { jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->JET_coneradius); } break; } } vector JetsUtil::RunJets(const vector& input_particles, const vector & TrackCentral, vector &NTrackJet, vector &EHADEEM,D_CaloTowerList list_of_active_towers) { inclusive_jets.clear(); sorted_jets.clear(); // run the jet clustering with the above jet definition if(input_particles.size()!=0) { fastjet::ClusterSequence clust_seq(input_particles, jet_def); // extract the inclusive jets with pt > 5 GeV double ptmin = 5.0; inclusive_jets = clust_seq.inclusive_jets(ptmin); // sort jets into increasing pt sorted_jets = sorted_by_pt(inclusive_jets); //Bin tracks to make the link float iEtaTrack[TrackCentral.size()]; float iPhiTrack[TrackCentral.size()]; for(unsigned int t = 0; t < TrackCentral.size(); t++) { if(!DET->JET_Eflow) DET->BinEtaPhi(TrackCentral[t].PhiOuter,TrackCentral[t].EtaOuter,iPhiTrack[t],iEtaTrack[t]); else { iPhiTrack[t] = TrackCentral[t].PhiOuter; iEtaTrack[t] = TrackCentral[t].EtaOuter; } } int numTrackJet; for (unsigned int i = 0; i < sorted_jets.size(); i++) { //check number of tracks in jets numTrackJet=0; vector constituents = clust_seq.constituents(sorted_jets[i]); for(unsigned int it = 0; it < TrackCentral.size(); it++) { for (unsigned int j = 0; j < constituents.size(); j++) { if(DeltaR(iPhiTrack[it], iEtaTrack[it], constituents[j].phi(), constituents[j].eta())<0.001)numTrackJet++; } } NTrackJet.push_back(numTrackJet); //now, get EHad over EEm float EmVal=0,HadVal=0; for (unsigned int j = 0; j < constituents.size(); j++) { //D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi())); // bug fix! bad association of phi [-pi;pi] <-> [0 ; 2pi] D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi_std())); EmVal += calConsti.getEem(); HadVal += calConsti.getEhad(); } EHADEEM.push_back(HadVal/EmVal); } } return sorted_jets; } vector JetsUtil::RunJetsResol(const vector& input_particles) { inclusive_jets.clear(); sorted_jets.clear(); // run the jet clustering with the above jet definition if(input_particles.size()!=0) { fastjet::ClusterSequence clust_seq(input_particles, jet_def); // extract the inclusive jets with pt > 5 GeV double ptmin = 5.0; inclusive_jets = clust_seq.inclusive_jets(ptmin); // sort jets into increasing pt sorted_jets = sorted_by_pt(inclusive_jets); } return sorted_jets; } void JetsUtil::RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector & sorted_jets,const TSimpleArray& NFCentralQ, const vector &NTrackJet, const vector &EHADEEM) { TRootJet *elementJet; TLorentzVector JET; for (unsigned int i = 0; i < sorted_jets.size(); i++) { JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E()); if(JET.Pt() > DET->PTCUT_jet) { elementJet = (TRootJet*) branchJet->NewEntry(); elementJet->Set(JET); elementJet->NTracks = NTrackJet[i]; elementJet->EHoverEE = EHADEEM[i]; // b-jets bool btag=false; if((fabs(JET.Eta()) < DET->CEN_max_tracker && DET->Btaggedjet(JET, NFCentralQ)))btag=true; elementJet->Btag = btag; } } // for itJet : loop on all jets } void JetsUtil::RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector & sorted_jets,const vector & towers, const vector & TrackCentral, const vector &NTrackJet, const vector &EHADEEM) { TRootTauJet *elementTauJet; TLorentzVector JET; float charge=0; for (unsigned int i = 0; i < sorted_jets.size(); i++) { JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E()); // Tau jet identification : 1! or 3! track and electromagnetic collimation if(fabs(JET.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) { double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi()); int NumTrackTau = DET->NumTracks(charge,TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi()); if( ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) && (( NumTrackTau == 1 ) || ( NumTrackTau == 3 )) && ( JET.Pt() > DET->PTCUT_taujet) ) { elementTauJet = (TRootTauJet*) branchTauJet->NewEntry(); elementTauJet->Set(JET); elementTauJet->NTracks=NTrackJet[i]; elementTauJet->Charge = charge; elementTauJet->EHoverEE = EHADEEM[i]; } // if tau jet } // if JET.eta < tracker - tau_cone : Tau jet identification } // for itJet : loop on all jets }