Fork me on GitHub

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

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3a73e6d was b443089, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix EOL characters in GPLv3 header

  • Property mode set to 100644
File size: 6.8 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/** \class PileUpJetID
20 *
21 * CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
22 *
23 * \author S. Zenz, December 2013
24 *
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// create output array(s)
94
95 fOutputArray = ExportArray(GetString("OutputArray", "jets"));
96
97}
98
99//------------------------------------------------------------------------------
100
101void PileUpJetID::Finish()
102{
103
104 if(fItJetInputArray) delete fItJetInputArray;
105 if(fItTrackInputArray) delete fItTrackInputArray;
106 if(fItNeutralInputArray) delete fItNeutralInputArray;
107 if(fItVertexInputArray) delete fItVertexInputArray;
108
109}
110
111//------------------------------------------------------------------------------
112
113void PileUpJetID::Process()
114{
115 Candidate *candidate, *constituent;
116 TLorentzVector momentum, area;
117 Double_t zvtx=0;
118
119 Candidate *trk;
120
121 // find z position of primary vertex
122
123 fItVertexInputArray->Reset();
124 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
125 {
126 if(!candidate->IsPU)
127 {
128 zvtx = candidate->Position.Z();
129 break;
130 }
131 }
132
133 // loop over all input candidates
134 fItJetInputArray->Reset();
135 while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
136 {
137 momentum = candidate->Momentum;
138 area = candidate->Area;
139
140 float sumpt = 0.;
141 float sumptch = 0.;
142 float sumptchpv = 0.;
143 float sumptchpu = 0.;
144 float sumdrsqptsq = 0.;
145 float sumptsq = 0.;
146 int nc = 0;
147 int nn = 0;
148 float pt_ann[5];
149
150 for (int i = 0 ; i < 5 ; i++) {
151 pt_ann[i] = 0.;
152 }
153
154 if (fUseConstituents) {
155 TIter itConstituents(candidate->GetCandidates());
156 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
157 float pt = constituent->Momentum.Pt();
158 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
159 sumpt += pt;
160 sumdrsqptsq += dr*dr*pt*pt;
161 sumptsq += pt*pt;
162 if (constituent->Charge == 0) {
163 // neutrals
164 nn++;
165 } else {
166 // charged
167 if (constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) {
168 sumptchpu += pt;
169 } else {
170 sumptchpv += pt;
171 }
172 sumptch += pt;
173 nc++;
174 }
175 for (int i = 0 ; i < 5 ; i++) {
176 if (dr > 0.1*i && dr < 0.1*(i+1)) {
177 pt_ann[i] += pt;
178 }
179 }
180 }
181 } else {
182 // Not using constituents, using dr
183 fItTrackInputArray->Reset();
184 while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
185 if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
186 float pt = trk->Momentum.Pt();
187 sumpt += pt;
188 sumptch += pt;
189 if (trk->IsPU && TMath::Abs(trk->Position.Z()-zvtx) > fZVertexResolution) {
190 sumptchpu += pt;
191 } else {
192 sumptchpv += pt;
193 }
194 float dr = candidate->Momentum.DeltaR(trk->Momentum);
195 sumdrsqptsq += dr*dr*pt*pt;
196 sumptsq += pt*pt;
197 nc++;
198 for (int i = 0 ; i < 5 ; i++) {
199 if (dr > 0.1*i && dr < 0.1*(i+1)) {
200 pt_ann[i] += pt;
201 }
202 }
203 }
204 }
205 fItNeutralInputArray->Reset();
206 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
207 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
208 float pt = constituent->Momentum.Pt();
209 sumpt += pt;
210 float dr = candidate->Momentum.DeltaR(constituent->Momentum);
211 sumdrsqptsq += dr*dr*pt*pt;
212 sumptsq += pt*pt;
213 nn++;
214 for (int i = 0 ; i < 5 ; i++) {
215 if (dr > 0.1*i && dr < 0.1*(i+1)) {
216 pt_ann[i] += pt;
217 }
218 }
219 }
220 }
221 }
222
223 if (sumptch > 0.) {
224 candidate->Beta = sumptchpu/sumptch;
225 candidate->BetaStar = sumptchpv/sumptch;
226 } else {
227 candidate->Beta = -999.;
228 candidate->BetaStar = -999.;
229 }
230 if (sumptsq > 0.) {
231 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
232 } else {
233 candidate->MeanSqDeltaR = -999.;
234 }
235 candidate->NCharged = nc;
236 candidate->NNeutrals = nn;
237 if (sumpt > 0.) {
238 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
239 for (int i = 0 ; i < 5 ; i++) {
240 candidate->FracPt[i] = pt_ann[i]/sumpt;
241 }
242 } else {
243 candidate->PTD = -999.;
244 for (int i = 0 ; i < 5 ; i++) {
245 candidate->FracPt[i] = -999.;
246 }
247 }
248
249 fOutputArray->Add(candidate);
250 }
251}
252
253//------------------------------------------------------------------------------
254
Note: See TracBrowser for help on using the repository browser.