Fork me on GitHub

source: git/modules/PileUpJetID.cc@ 9040259

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 9040259 was 2862770, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago

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