Fork me on GitHub

source: svn/trunk/src/JetUtils.cc@ 145

Last change on this file since 145 was 100, checked in by severine ovyn, 16 years ago

Remove datacard bug + CaloTowers OK

File size: 4.2 KB
Line 
1/*
2 * ---- Delphes ----
3 * A Fast Simulator for general purpose LHC detector
4 * S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5 *
6 * Center for Particle Physics and Phenomenology (CP3)
7 * Universite Catholique de Louvain (UCL)
8 * Louvain-la-Neuve, Belgium
9 * */
10
11// \brief Trigger class, and some generic definitions
12
13#include "interface/JetUtils.h"
14
15using namespace std;
16
17JetsUtil::JetsUtil(const string DetDatacard) {
18
19 DET = new RESOLution();
20 DET->ReadDataCard(DetDatacard);
21
22 switch(DET->JET_jetalgo) {
23 default:
24 case 1: {
25 // set up a CDF midpoint jet definition
26#ifdef ENABLE_PLUGIN_CDFCONES
27 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);
28 jet_def = fastjet::JetDefinition(plugins);
29#else
30 plugins = NULL;
31#endif
32 }
33 break;
34
35 case 2: {
36 // set up a CDF midpoint jet definition
37#ifdef ENABLE_PLUGIN_CDFCONES
38 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);
39 jet_def = fastjet::JetDefinition(plugins);
40#else
41 plugins = NULL;
42#endif
43 }
44 break;
45
46 case 3: {
47 // set up a siscone jet definition
48#ifdef ENABLE_PLUGIN_SISCONE
49 plugins = new fastjet::SISConePlugin (DET->JET_coneradius,DET->JET_overlap,DET->JET_S_npass, DET->JET_S_protojet_ptmin);
50 jet_def = fastjet::JetDefinition(plugins);
51#else
52 plugins = NULL;
53#endif
54 }
55 break;
56
57 case 4: {
58
59 jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->JET_coneradius);
60 }
61 break;
62
63 case 5: {
64 jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->JET_coneradius);
65 }
66 break;
67
68 case 6: {
69 jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->JET_coneradius);
70 }
71 break;
72 }
73
74}
75
76vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet>& input_particles)
77{
78 inclusive_jets.clear();
79 sorted_jets.clear();
80 // run the jet clustering with the above jet definition
81 if(input_particles.size()!=0)
82 {
83 fastjet::ClusterSequence clust_seq(input_particles, jet_def);
84 // extract the inclusive jets with pt > 5 GeV
85 double ptmin = 5.0;
86 inclusive_jets = clust_seq.inclusive_jets(ptmin);
87 // sort jets into increasing pt
88 sorted_jets = sorted_by_pt(inclusive_jets);
89 }
90
91 return sorted_jets;
92}
93
94void JetsUtil::RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector<fastjet::PseudoJet> & sorted_jets,const TSimpleArray<TRootGenParticle>& NFCentralQ)
95{
96 TRootJet *elementJet;
97 TLorentzVector JET;
98 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
99 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
100 if(JET.Pt() > DET->PTCUT_jet)
101 {
102 elementJet = (TRootJet*) branchJet->NewEntry();
103 elementJet->Set(JET);
104 // b-jets
105 bool btag=false;
106 if((fabs(JET.Eta()) < DET->CEN_max_tracker && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
107 elementJet->Btag = btag;
108 }
109 } // for itJet : loop on all jets
110
111}
112
113void JetsUtil::RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector<fastjet::PseudoJet> & sorted_jets,const vector<PhysicsTower> & towers, const vector<TLorentzVector> & TrackCentral)
114{
115 TRootTauJet *elementTauJet;
116 TLorentzVector JET;
117 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
118 JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
119 // Tau jet identification : 1! track and electromagnetic collimation
120 if(fabs(JET.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) {
121 double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
122 if(
123 ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) &&
124 ( DET->NumTracks(TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi()) == 1 ) &&
125 ( JET.Pt() > DET->PTCUT_taujet)
126 ) {
127 elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
128 elementTauJet->Set(JET);
129 } // if tau jet
130 } // if JET.eta < tracker - tau_cone : Tau jet identification
131
132 } // for itJet : loop on all jets
133
134
135}
136
Note: See TracBrowser for help on using the repository browser.