Fork me on GitHub

source: git/modules/PileUpJetID.cc@ 83eff6d

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 83eff6d was 54b6dfc, checked in by mselvaggi <mselvaggi@…>, 11 years ago

JetId v2

  • Property mode set to 100644
File size: 6.0 KB
Line 
1/** \class PileUpJetID
2 *
3 * CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
4 *
5 * \author S. Zenz, December 2013
6 *
7 *
8 */
9
10#include "modules/PileUpJetID.h"
11
12#include "classes/DelphesClasses.h"
13#include "classes/DelphesFactory.h"
14#include "classes/DelphesFormula.h"
15
16#include "ExRootAnalysis/ExRootResult.h"
17#include "ExRootAnalysis/ExRootFilter.h"
18#include "ExRootAnalysis/ExRootClassifier.h"
19
20#include "TMath.h"
21#include "TString.h"
22#include "TFormula.h"
23#include "TRandom3.h"
24#include "TObjArray.h"
25#include "TDatabasePDG.h"
26#include "TLorentzVector.h"
27
28#include <algorithm>
29#include <stdexcept>
30#include <iostream>
31#include <sstream>
32
33using namespace std;
34
35//------------------------------------------------------------------------------
36
37PileUpJetID::PileUpJetID() :
38 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
39{
40
41}
42
43//------------------------------------------------------------------------------
44
45PileUpJetID::~PileUpJetID()
46{
47
48}
49
50//------------------------------------------------------------------------------
51
52void PileUpJetID::Init()
53{
54 fJetPTMin = GetDouble("JetPTMin", 20.0);
55 fParameterR = GetDouble("ParameterR", 0.5);
56 fUseConstituents = GetInt("UseConstituents", 0);
57
58 fAverageEachTower = false; // for timing
59
60 // import input array(s)
61
62 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
63 fItJetInputArray = fJetInputArray->MakeIterator();
64
65 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
66 fItTrackInputArray = fTrackInputArray->MakeIterator();
67
68 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
69 fItNeutralInputArray = fNeutralInputArray->MakeIterator();
70
71 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
72 fItVertexInputArray = fVertexInputArray->MakeIterator();
73
74 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3;
75// create output array(s)
76
77 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
78
79}
80
81//------------------------------------------------------------------------------
82
83void PileUpJetID::Finish()
84{
85
86 if(fItJetInputArray) delete fItJetInputArray;
87 if(fItTrackInputArray) delete fItTrackInputArray;
88 if(fItNeutralInputArray) delete fItNeutralInputArray;
89 if(fItVertexInputArray) delete fItVertexInputArray;
90
91}
92
93//------------------------------------------------------------------------------
94
95void PileUpJetID::Process()
96{
97 Candidate *candidate, *constituent;
98 TLorentzVector momentum, area;
99 Double_t zvtx=0;
100
101 Candidate *trk;
102
103 // find z position of primary vertex
104
105 fItVertexInputArray->Reset();
106 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
107 {
108 if(!candidate->IsPU)
109 {
110 zvtx = candidate->Position.Z();
111 break;
112 }
113 }
114
115 // loop over all input candidates
116 fItJetInputArray->Reset();
117 while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
118 {
119 momentum = candidate->Momentum;
120 area = candidate->Area;
121
122 float sumpt = 0.;
123 float sumptch = 0.;
124 float sumptchpv = 0.;
125 float sumptchpu = 0.;
126 float sumdrsqptsq = 0.;
127 float sumptsq = 0.;
128 int nc = 0;
129 int nn = 0;
130 float pt_ann[5];
131
132 for (int i = 0 ; i < 5 ; i++) {
133 pt_ann[i] = 0.;
134 }
135
136 if (fUseConstituents) {
137 TIter itConstituents(candidate->GetCandidates());
138 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
139 float pt = constituent->Momentum.Pt();
140 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
141 sumpt += pt;
142 sumdrsqptsq += dr*dr*pt*pt;
143 sumptsq += pt*pt;
144 if (constituent->Charge == 0) {
145 // neutrals
146 nn++;
147 } else {
148 // charged
149 if (constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) {
150 sumptchpu += pt;
151 } else {
152 sumptchpv += pt;
153 }
154 sumptch += pt;
155 nc++;
156 }
157 for (int i = 0 ; i < 5 ; i++) {
158 if (dr > 0.1*i && dr < 0.1*(i+1)) {
159 pt_ann[i] += pt;
160 }
161 }
162 }
163 } else {
164 // Not using constituents, using dr
165 fItTrackInputArray->Reset();
166 while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
167 if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
168 float pt = trk->Momentum.Pt();
169 sumpt += pt;
170 sumptch += pt;
171 if (trk->IsPU && TMath::Abs(trk->Position.Z()-zvtx) > fZVertexResolution) {
172 sumptchpu += pt;
173 } else {
174 sumptchpv += pt;
175 }
176 float dr = candidate->Momentum.DeltaR(trk->Momentum);
177 sumdrsqptsq += dr*dr*pt*pt;
178 sumptsq += pt*pt;
179 nc++;
180 for (int i = 0 ; i < 5 ; i++) {
181 if (dr > 0.1*i && dr < 0.1*(i+1)) {
182 pt_ann[i] += pt;
183 }
184 }
185 }
186 }
187 fItNeutralInputArray->Reset();
188 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
189 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
190 float pt = constituent->Momentum.Pt();
191 sumpt += pt;
192 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
193 sumdrsqptsq += dr*dr*pt*pt;
194 sumptsq += pt*pt;
195 nn++;
196 for (int i = 0 ; i < 5 ; i++) {
197 if (dr > 0.1*i && dr < 0.1*(i+1)) {
198 pt_ann[i] += pt;
199 }
200 }
201 }
202 }
203 }
204
205 if (sumptch > 0.) {
206 candidate->Beta = sumptchpu/sumptch;
207 candidate->BetaStar = sumptchpv/sumptch;
208 } else {
209 candidate->Beta = -999.;
210 candidate->BetaStar = -999.;
211 }
212 if (sumptsq > 0.) {
213 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
214 } else {
215 candidate->MeanSqDeltaR = -999.;
216 }
217 candidate->NCharged = nc;
218 candidate->NNeutrals = nn;
219 if (sumpt > 0.) {
220 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
221 for (int i = 0 ; i < 5 ; i++) {
222 candidate->FracPt[i] = pt_ann[i]/sumpt;
223 }
224 } else {
225 candidate->PTD = -999.;
226 for (int i = 0 ; i < 5 ; i++) {
227 candidate->FracPt[i] = -999.;
228 }
229 }
230
231 fOutputArray->Add(candidate);
232 }
233}
234
235//------------------------------------------------------------------------------
236
Note: See TracBrowser for help on using the repository browser.