Fork me on GitHub

source: git/modules/PileUpJetID.cc@ 3c46e17

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3c46e17 was 84edab9, checked in by Michele Selvaggi <michele.selvaggi@…>, 9 years ago

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