Fork me on GitHub

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

Last change on this file since 425 was 384, checked in by Xavier Rouby, 16 years ago

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