Fork me on GitHub

source: git/modules/PileUpJetID.cc@ a3c2d52

Last change on this file since a3c2d52 was 6cdc544, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix formatting

  • Property mode set to 100644
File size: 7.3 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class PileUpJetID
21 *
22 * CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
23 *
24 * \author S. Zenz, December 2013
25 *
26 */
27
28#include "modules/PileUpJetID.h"
29
30#include "classes/DelphesClasses.h"
31#include "classes/DelphesFactory.h"
32#include "classes/DelphesFormula.h"
33
34#include "ExRootAnalysis/ExRootResult.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootClassifier.h"
37
38#include "TMath.h"
39#include "TString.h"
40#include "TFormula.h"
41#include "TRandom3.h"
42#include "TObjArray.h"
43#include "TDatabasePDG.h"
44#include "TLorentzVector.h"
45
46#include <algorithm>
47#include <stdexcept>
48#include <iostream>
49#include <sstream>
50
51using namespace std;
52
53//------------------------------------------------------------------------------
54
55PileUpJetID::PileUpJetID() :
56 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
57{
58
59}
60
61//------------------------------------------------------------------------------
62
63PileUpJetID::~PileUpJetID()
64{
65
66}
67
68//------------------------------------------------------------------------------
69
70void PileUpJetID::Init()
71{
72 fJetPTMin = GetDouble("JetPTMin", 20.0);
73 fParameterR = GetDouble("ParameterR", 0.5);
74 fUseConstituents = GetInt("UseConstituents", 0);
75
76 fAverageEachTower = false; // for timing
77
78 // import input array(s)
79
80 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
81 fItJetInputArray = fJetInputArray->MakeIterator();
82
83 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
84 fItTrackInputArray = fTrackInputArray->MakeIterator();
85
86 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
87 fItNeutralInputArray = fNeutralInputArray->MakeIterator();
88
89 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
90 fItVertexInputArray = fVertexInputArray->MakeIterator();
91
92 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3;
93
94 // create output array(s)
95
96 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
97}
98
99//------------------------------------------------------------------------------
100
101void PileUpJetID::Finish()
102{
103 if(fItJetInputArray) delete fItJetInputArray;
104 if(fItTrackInputArray) delete fItTrackInputArray;
105 if(fItNeutralInputArray) delete fItNeutralInputArray;
106 if(fItVertexInputArray) delete fItVertexInputArray;
107}
108
109//------------------------------------------------------------------------------
110
111void PileUpJetID::Process()
112{
113 Candidate *candidate, *constituent;
114 TLorentzVector momentum, area;
115 Int_t i, nc, nn;
116 Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
117 Double_t dr, pt, pt_ann[5];
118 Double_t zvtx = 0.0;
119
120 Candidate *track;
121
122 // find z position of primary vertex
123
124 fItVertexInputArray->Reset();
125 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
126 {
127 if(!candidate->IsPU)
128 {
129 zvtx = candidate->Position.Z();
130 break;
131 }
132 }
133
134 // loop over all input candidates
135 fItJetInputArray->Reset();
136 while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
137 {
138 momentum = candidate->Momentum;
139 area = candidate->Area;
140
141 sumpt = 0.0;
142 sumptch = 0.0;
143 sumptchpv = 0.0;
144 sumptchpu = 0.0;
145 sumdrsqptsq = 0.0;
146 sumptsq = 0.0;
147 nc = 0;
148 nn = 0;
149
150 for(i = 0; i < 5; ++i)
151 {
152 pt_ann[i] = 0.0;
153 }
154
155 if(fUseConstituents)
156 {
157 TIter itConstituents(candidate->GetCandidates());
158 while((constituent = static_cast<Candidate*>(itConstituents.Next())))
159 {
160 pt = constituent->Momentum.Pt();
161 dr = candidate->Momentum.DeltaR(constituent->Momentum);
162 sumpt += pt;
163 sumdrsqptsq += dr*dr*pt*pt;
164 sumptsq += pt*pt;
165 if(constituent->Charge == 0)
166 {
167 // neutrals
168 ++nn;
169 }
170 else
171 {
172 // charged
173 if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
174 {
175 sumptchpu += pt;
176 }
177 else
178 {
179 sumptchpv += pt;
180 }
181 sumptch += pt;
182 ++nc;
183 }
184 for(i = 0; i < 5; ++i)
185 {
186 if(dr > 0.1*i && dr < 0.1*(i + 1))
187 {
188 pt_ann[i] += pt;
189 }
190 }
191 }
192 }
193 else
194 {
195 // Not using constituents, using dr
196 fItTrackInputArray->Reset();
197 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
198 {
199 if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
200 {
201 pt = track->Momentum.Pt();
202 sumpt += pt;
203 sumptch += pt;
204 if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
205 {
206 sumptchpu += pt;
207 }
208 else
209 {
210 sumptchpv += pt;
211 }
212 dr = candidate->Momentum.DeltaR(track->Momentum);
213 sumdrsqptsq += dr*dr*pt*pt;
214 sumptsq += pt*pt;
215 nc++;
216 for(i = 0; i < 5; ++i)
217 {
218 if(dr > 0.1*i && dr < 0.1*(i + 1))
219 {
220 pt_ann[i] += pt;
221 }
222 }
223 }
224 }
225
226 fItNeutralInputArray->Reset();
227 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
228 {
229 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
230 {
231 pt = constituent->Momentum.Pt();
232 sumpt += pt;
233 dr = candidate->Momentum.DeltaR(constituent->Momentum);
234 sumdrsqptsq += dr*dr*pt*pt;
235 sumptsq += pt*pt;
236 nn++;
237 for(i = 0; i < 5; ++i)
238 {
239 if(dr > 0.1*i && dr < 0.1*(i + 1))
240 {
241 pt_ann[i] += pt;
242 }
243 }
244 }
245 }
246 }
247
248 if(sumptch > 0.0)
249 {
250 candidate->Beta = sumptchpu/sumptch;
251 candidate->BetaStar = sumptchpv/sumptch;
252 }
253 else
254 {
255 candidate->Beta = -999.0;
256 candidate->BetaStar = -999.0;
257 }
258 if(sumptsq > 0.0)
259 {
260 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
261 }
262 else
263 {
264 candidate->MeanSqDeltaR = -999.0;
265 }
266 candidate->NCharged = nc;
267 candidate->NNeutrals = nn;
268 if(sumpt > 0.0)
269 {
270 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
271 for(i = 0; i < 5; ++i)
272 {
273 candidate->FracPt[i] = pt_ann[i]/sumpt;
274 }
275 }
276 else
277 {
278 candidate->PTD = -999.0;
279 for(i = 0; i < 5; ++i)
280 {
281 candidate->FracPt[i] = -999.0;
282 }
283 }
284
285 fOutputArray->Add(candidate);
286 }
287}
288
289//------------------------------------------------------------------------------
290
Note: See TracBrowser for help on using the repository browser.