Fork me on GitHub

source: git/modules/PileUpJetID.cc@ a5af1df

Last change on this file since a5af1df was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files
  • Property mode set to 100644
File size: 10.3 KB
Line 
1
2/** \class PileUpJetID
3 *
4 * CMS PileUp Jet ID Variables
5 *
6 * \author S. Zenz
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/ExRootClassifier.h"
17#include "ExRootAnalysis/ExRootFilter.h"
18#include "ExRootAnalysis/ExRootResult.h"
19
20#include "TFormula.h"
21#include "TMath.h"
22#include "TObjArray.h"
23#include "TRandom3.h"
24#include "TString.h"
25//#include "TDatabasePDG.h"
26#include "TLorentzVector.h"
27
28#include <algorithm>
29#include <iostream>
30#include <sstream>
31#include <stdexcept>
32
33using namespace std;
34
35//------------------------------------------------------------------------------
36
37PileUpJetID::PileUpJetID() :
38 fItJetInputArray(0), fTrackInputArray(0), fNeutralInputArray(0)
39{
40}
41
42//------------------------------------------------------------------------------
43
44PileUpJetID::~PileUpJetID()
45{
46}
47
48//------------------------------------------------------------------------------
49
50void PileUpJetID::Init()
51{
52
53 fJetPTMin = GetDouble("JetPTMin", 20.0);
54 fParameterR = GetDouble("ParameterR", 0.5);
55 fUseConstituents = GetInt("UseConstituents", 0);
56
57 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel", 0.1);
58 fBetaMinBarrel = GetDouble("BetaMinBarrel", 0.1);
59 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap", 0.1);
60 fBetaMinEndcap = GetDouble("BetaMinEndcap", 0.1);
61 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward", 0.1);
62 fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
63 fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
64
65 fAverageEachTower = false; // for timing
66
67 // import input array(s)
68
69 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
70 fItJetInputArray = fJetInputArray->MakeIterator();
71
72 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
73 fItTrackInputArray = fTrackInputArray->MakeIterator();
74
75 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
76 fItNeutralInputArray = fNeutralInputArray->MakeIterator();
77
78 // create output array(s)
79
80 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
81
82 fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets", "eflowtowers"));
83}
84
85//------------------------------------------------------------------------------
86
87void PileUpJetID::Finish()
88{
89 // cout << "In finish" << endl;
90
91 if(fItJetInputArray) delete fItJetInputArray;
92 if(fItTrackInputArray) delete fItTrackInputArray;
93 if(fItNeutralInputArray) delete fItNeutralInputArray;
94}
95
96//------------------------------------------------------------------------------
97
98void PileUpJetID::Process()
99{
100 Candidate *candidate, *constituent;
101 TLorentzVector momentum, area;
102
103 Candidate *trk;
104
105 // loop over all input candidates
106 fItJetInputArray->Reset();
107 while((candidate = static_cast<Candidate *>(fItJetInputArray->Next())))
108 {
109 momentum = candidate->Momentum;
110 area = candidate->Area;
111
112 float sumT0 = 0.;
113 float sumT1 = 0.;
114 float sumT10 = 0.;
115 float sumT20 = 0.;
116 float sumT30 = 0.;
117 float sumT40 = 0.;
118 float sumWeightsForT = 0.;
119 candidate->NTimeHits = 0;
120
121 float sumpt = 0.;
122 float sumptch = 0.;
123 float sumptchpv = 0.;
124 float sumptchpu = 0.;
125 float sumdrsqptsq = 0.;
126 float sumptsq = 0.;
127 int nc = 0;
128 int nn = 0;
129 float pt_ann[5];
130
131 for(int i = 0; i < 5; i++)
132 {
133 pt_ann[i] = 0.;
134 }
135
136 if(fUseConstituents)
137 {
138 TIter itConstituents(candidate->GetCandidates());
139 while((constituent = static_cast<Candidate *>(itConstituents.Next())))
140 {
141 float pt = constituent->Momentum.Pt();
142 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
143 // cout << " There exists a constituent with dr=" << dr << endl;
144 sumpt += pt;
145 sumdrsqptsq += dr * dr * pt * pt;
146 sumptsq += pt * pt;
147 if(constituent->Charge == 0)
148 {
149 nn++;
150 }
151 else
152 {
153 if(constituent->IsRecoPU)
154 {
155 sumptchpu += pt;
156 }
157 else
158 {
159 sumptchpv += pt;
160 }
161 sumptch += pt;
162 nc++;
163 }
164 for(int i = 0; i < 5; i++)
165 {
166 if(dr > 0.1 * i && dr < 0.1 * (i + 1))
167 {
168 pt_ann[i] += pt;
169 }
170 }
171 float tow_sumT = 0;
172 float tow_sumW = 0;
173 for(int i = 0; i < constituent->ECalEnergyTimePairs.size(); i++)
174 {
175 float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first);
176 if(fAverageEachTower)
177 {
178 tow_sumT += w * constituent->ECalEnergyTimePairs[i].second;
179 tow_sumW += w;
180 }
181 else
182 {
183 sumT0 += w * constituent->ECalEnergyTimePairs[i].second;
184 sumT1 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.001);
185 sumT10 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.010);
186 sumT20 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.020);
187 sumT30 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.030);
188 sumT40 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.040);
189 sumWeightsForT += w;
190 candidate->NTimeHits++;
191 }
192 }
193 if(fAverageEachTower && tow_sumW > 0.)
194 {
195 sumT0 += tow_sumT;
196 sumT1 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.001);
197 sumT10 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0010);
198 sumT20 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0020);
199 sumT30 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0030);
200 sumT40 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0040);
201 sumWeightsForT += tow_sumW;
202 candidate->NTimeHits++;
203 }
204 }
205 }
206 else
207 {
208 // Not using constituents, using dr
209 fItTrackInputArray->Reset();
210 while((trk = static_cast<Candidate *>(fItTrackInputArray->Next())))
211 {
212 if(trk->Momentum.DeltaR(candidate->Momentum) < fParameterR)
213 {
214 float pt = trk->Momentum.Pt();
215 sumpt += pt;
216 sumptch += pt;
217 if(trk->IsRecoPU)
218 {
219 sumptchpu += pt;
220 }
221 else
222 {
223 sumptchpv += pt;
224 }
225 float dr = candidate->Momentum.DeltaR(trk->Momentum);
226 sumdrsqptsq += dr * dr * pt * pt;
227 sumptsq += pt * pt;
228 nc++;
229 for(int i = 0; i < 5; i++)
230 {
231 if(dr > 0.1 * i && dr < 0.1 * (i + 1))
232 {
233 pt_ann[i] += pt;
234 }
235 }
236 }
237 }
238 fItNeutralInputArray->Reset();
239 while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next())))
240 {
241 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
242 {
243 float pt = constituent->Momentum.Pt();
244 sumpt += pt;
245 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
246 sumdrsqptsq += dr * dr * pt * pt;
247 sumptsq += pt * pt;
248 nn++;
249 for(int i = 0; i < 5; i++)
250 {
251 if(dr > 0.1 * i && dr < 0.1 * (i + 1))
252 {
253 pt_ann[i] += pt;
254 }
255 }
256 }
257 }
258 }
259
260 if(sumptch > 0.)
261 {
262 candidate->Beta = sumptchpv / sumptch;
263 candidate->BetaStar = sumptchpu / sumptch;
264 }
265 else
266 {
267 candidate->Beta = -999.;
268 candidate->BetaStar = -999.;
269 }
270 if(sumptsq > 0.)
271 {
272 candidate->MeanSqDeltaR = sumdrsqptsq / sumptsq;
273 }
274 else
275 {
276 candidate->MeanSqDeltaR = -999.;
277 }
278 candidate->NCharged = nc;
279 candidate->NNeutrals = nn;
280 if(sumpt > 0.)
281 {
282 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
283 for(int i = 0; i < 5; i++)
284 {
285 candidate->FracPt[i] = pt_ann[i] / sumpt;
286 }
287 }
288 else
289 {
290 candidate->PTD = -999.;
291 for(int i = 0; i < 5; i++)
292 {
293 candidate->FracPt[i] = -999.;
294 }
295 }
296
297 fOutputArray->Add(candidate);
298
299 // New stuff
300 /*
301 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
302 fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
303 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
304 fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
305 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
306 */
307
308 bool passId = false;
309 if(candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1)
310 {
311 if(fabs(candidate->Momentum.Eta()) < 1.5)
312 {
313 passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
314 }
315 else if(fabs(candidate->Momentum.Eta()) < 4.0)
316 {
317 passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
318 }
319 else
320 {
321 passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
322 }
323 }
324
325 // cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
326 // << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
327
328 if(passId)
329 {
330 if(fUseConstituents)
331 {
332 TIter itConstituents(candidate->GetCandidates());
333 while((constituent = static_cast<Candidate *>(itConstituents.Next())))
334 {
335 if(constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin)
336 {
337 fNeutralsInPassingJets->Add(constituent);
338 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
339 }
340 }
341 }
342 else
343 { // use DeltaR
344 fItNeutralInputArray->Reset();
345 while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next())))
346 {
347 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin)
348 {
349 fNeutralsInPassingJets->Add(constituent);
350 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
351 }
352 }
353 }
354 }
355 }
356}
357
358//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.