Fork me on GitHub

source: svn/trunk/src/JetsUtil.cc@ 559

Last change on this file since 559 was 547, checked in by severine ovyn, 15 years ago

bug fix(jet Ntracks), 3-prong taus, vertex coordinates for tracks
File size: 9.5 KB
RevLine 
[260]1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
[443]6** \------------- arXiv:0903.2225v1 ------------/ **
[260]7** **
8** **
9** This package uses: **
10** ------------------ **
[443]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] **
[260]14** FROG: [hep-ex/0901.2718v1] **
[443]15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
[260]16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
[443]22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
[260]24** **
[443]25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
[260]29** Copyright (C) 2008-2009, **
[443]30** All rights reserved. **
[260]31** **
32***********************************************************************/
[215]33
34#include "JetsUtil.h"
35
36using namespace std;
37
38JetsUtil::JetsUtil() : plugins(NULL) {
39 DET = new RESOLution();
40 init();
41}
42
43
44JetsUtil::JetsUtil(const string & DetDatacard) : plugins(NULL) {
45 DET = new RESOLution();
46 DET->ReadDataCard(DetDatacard);
47 init();
48}
49
50
51JetsUtil::JetsUtil(const RESOLution * DetDatacard) : plugins(NULL) {
52 DET = new RESOLution(*DetDatacard);
53 init();
54}
55
56
57JetsUtil::JetsUtil(const JetsUtil& jet){
58 DET = new RESOLution(*(jet.DET));
59}
60
61JetsUtil& JetsUtil::operator=(const JetsUtil& jet) {
62 if(this==&jet) return *this;
63 DET = new RESOLution(*(jet.DET));
64 return *this;
65}
66
67void JetsUtil::init() {
68 if(plugins) delete plugins;
69
70 switch(DET->JET_jetalgo) {
71 default:
72 case 1: {
73 // set up a CDF midpoint jet definition
74#ifdef ENABLE_PLUGIN_CDFCONES
75 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);
76 jet_def = fastjet::JetDefinition(plugins);
77#else
78 plugins = NULL;
79#endif
80 }
81 break;
82
83 case 2: {
84 // set up a CDF midpoint jet definition
85#ifdef ENABLE_PLUGIN_CDFCONES
86 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);
87 jet_def = fastjet::JetDefinition(plugins);
88#else
89 plugins = NULL;
90#endif
91 }
92 break;
93
94 case 3: {
95 // set up a siscone jet definition
96#ifdef ENABLE_PLUGIN_SISCONE
97 plugins = new fastjet::SISConePlugin (DET->JET_coneradius,DET->JET_overlap,DET->JET_S_npass, DET->JET_S_protojet_ptmin);
98 jet_def = fastjet::JetDefinition(plugins);
99#else
100 plugins = NULL;
101#endif
102 }
103 break;
104
105 case 4: {
106
107 jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->JET_coneradius);
108 }
109 break;
110
111 case 5: {
112 jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->JET_coneradius);
113 }
114 break;
115
116 case 6: {
117 jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->JET_coneradius);
118 }
119 break;
120 }
121
122}
123
[310]124vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet>& input_particles, const vector<TRootTracks> & TrackCentral, vector<int> &NTrackJet, vector<float> &EHADEEM,D_CaloTowerList list_of_active_towers)
[215]125{
126 inclusive_jets.clear();
127 sorted_jets.clear();
128 // run the jet clustering with the above jet definition
129 if(input_particles.size()!=0)
130 {
131 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
132 // extract the inclusive jets with pt > 5 GeV
133 double ptmin = 5.0;
134 inclusive_jets = clust_seq.inclusive_jets(ptmin);
[307]135
[215]136 // sort jets into increasing pt
137 sorted_jets = sorted_by_pt(inclusive_jets);
[310]138 //Bin tracks to make the link
[307]139 float iEtaTrack[TrackCentral.size()];
140 float iPhiTrack[TrackCentral.size()];
141 for(unsigned int t = 0; t < TrackCentral.size(); t++)
[384]142 {
143 if(!DET->JET_Eflow)
144 DET->BinEtaPhi(TrackCentral[t].PhiOuter,TrackCentral[t].EtaOuter,iPhiTrack[t],iEtaTrack[t]);
145 else {
146 iPhiTrack[t] = TrackCentral[t].PhiOuter; iEtaTrack[t] = TrackCentral[t].EtaOuter;
147 }
148 }
[307]149 int numTrackJet;
150 for (unsigned int i = 0; i < sorted_jets.size(); i++)
151 {
[310]152 //check number of tracks in jets
[307]153 numTrackJet=0;
154 vector<fastjet::PseudoJet> constituents = clust_seq.constituents(sorted_jets[i]);
155 for(unsigned int it = 0; it < TrackCentral.size(); it++)
156 {
157 for (unsigned int j = 0; j < constituents.size(); j++)
158 {
159 if(DeltaR(iPhiTrack[it], iEtaTrack[it], constituents[j].phi(), constituents[j].eta())<0.001)numTrackJet++;
160 }
161 }
162 NTrackJet.push_back(numTrackJet);
[310]163 //now, get EHad over EEm
164 float EmVal=0,HadVal=0;
165 for (unsigned int j = 0; j < constituents.size(); j++)
166 {
[547]167 //D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi()));
168 // bug fix! bad association of phi [-pi;pi] <-> [0 ; 2pi]
169 D_CaloTower calConsti(list_of_active_towers.getElement(constituents[j].eta(),constituents[j].phi_std()));
[310]170 EmVal += calConsti.getEem();
171 HadVal += calConsti.getEhad();
172 }
173 EHADEEM.push_back(HadVal/EmVal);
[307]174 }
[215]175 }
176
177 return sorted_jets;
178}
179
[310]180vector<fastjet::PseudoJet> JetsUtil::RunJetsResol(const vector<fastjet::PseudoJet>& input_particles)
181{
182 inclusive_jets.clear();
183 sorted_jets.clear();
184 // run the jet clustering with the above jet definition
185 if(input_particles.size()!=0)
186 {
187 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
188 // extract the inclusive jets with pt > 5 GeV
189 double ptmin = 5.0;
190 inclusive_jets = clust_seq.inclusive_jets(ptmin);
[307]191
[310]192 // sort jets into increasing pt
193 sorted_jets = sorted_by_pt(inclusive_jets);
194 }
195 return sorted_jets;
196}
[307]197
[310]198
[350]199void JetsUtil::RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector<fastjet::PseudoJet> & sorted_jets,const TSimpleArray<TRootC::GenParticle>& NFCentralQ, const vector<int> &NTrackJet, const vector<float> &EHADEEM)
[215]200{
201 TRootJet *elementJet;
202 TLorentzVector JET;
203 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
204 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
205 if(JET.Pt() > DET->PTCUT_jet)
206 {
207 elementJet = (TRootJet*) branchJet->NewEntry();
208 elementJet->Set(JET);
[307]209 elementJet->NTracks = NTrackJet[i];
[310]210 elementJet->EHoverEE = EHADEEM[i];
[307]211
[215]212 // b-jets
213 bool btag=false;
214 if((fabs(JET.Eta()) < DET->CEN_max_tracker && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
215 elementJet->Btag = btag;
216 }
217 } // for itJet : loop on all jets
218}
219
[310]220void JetsUtil::RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector<fastjet::PseudoJet> & sorted_jets,const vector<PhysicsTower> & towers, const vector<TRootTracks> & TrackCentral, const vector<int> &NTrackJet, const vector<float> &EHADEEM)
[215]221{
222 TRootTauJet *elementTauJet;
223 TLorentzVector JET;
[286]224 float charge=0;
[315]225
[215]226 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
227 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
[547]228 // Tau jet identification : 1! or 3! track and electromagnetic collimation
[215]229 if(fabs(JET.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) {
230 double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
[547]231 int NumTrackTau = DET->NumTracks(charge,TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi());
[215]232 if(
233 ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) &&
[547]234 (( NumTrackTau == 1 ) || ( NumTrackTau == 3 )) &&
[215]235 ( JET.Pt() > DET->PTCUT_taujet)
236 ) {
237 elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
238 elementTauJet->Set(JET);
[307]239 elementTauJet->NTracks=NTrackJet[i];
[286]240 elementTauJet->Charge = charge;
[310]241 elementTauJet->EHoverEE = EHADEEM[i];
[215]242 } // if tau jet
243 } // if JET.eta < tracker - tau_cone : Tau jet identification
244 } // for itJet : loop on all jets
245}
Note: See TracBrowser for help on using the repository browser.