Fork me on GitHub

source: git/modules/Isolation.cc@ f32ea92

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f32ea92 was fdfac34, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

fix sumCharged in Isolation

  • Property mode set to 100644
File size: 7.2 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 Isolation
21 *
22 * Sums transverse momenta of isolation objects (tracks, calorimeter towers, etc)
23 * within a DeltaR cone around a candidate and calculates fraction of this sum
24 * to the candidate's transverse momentum. outputs candidates that have
25 * the transverse momenta fraction within (PTRatioMin, PTRatioMax].
26 *
27 * \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
28 *
29 */
30
31#include "modules/Isolation.h"
32
33#include "classes/DelphesClasses.h"
34#include "classes/DelphesFactory.h"
35#include "classes/DelphesFormula.h"
36
37#include "ExRootAnalysis/ExRootResult.h"
38#include "ExRootAnalysis/ExRootFilter.h"
39#include "ExRootAnalysis/ExRootClassifier.h"
40
41#include "TMath.h"
42#include "TString.h"
43#include "TFormula.h"
44#include "TRandom3.h"
45#include "TObjArray.h"
46#include "TDatabasePDG.h"
47#include "TLorentzVector.h"
48
49#include <algorithm>
50#include <stdexcept>
51#include <iostream>
52#include <sstream>
53
54using namespace std;
55
56//------------------------------------------------------------------------------
57
58class IsolationClassifier : public ExRootClassifier
59{
60public:
61
62 IsolationClassifier() {}
63
64 Int_t GetCategory(TObject *object);
65
66 Double_t fPTMin;
67};
68
69//------------------------------------------------------------------------------
70
71Int_t IsolationClassifier::GetCategory(TObject *object)
72{
73 Candidate *track = static_cast<Candidate*>(object);
74 const TLorentzVector &momentum = track->Momentum;
75
76 if(momentum.Pt() < fPTMin) return -1;
77
78 return 0;
79}
80
81//------------------------------------------------------------------------------
82
83Isolation::Isolation() :
84 fClassifier(0), fFilter(0),
85 fItIsolationInputArray(0), fItCandidateInputArray(0),
86 fItRhoInputArray(0)
87{
88 fClassifier = new IsolationClassifier;
89}
90
91//------------------------------------------------------------------------------
92
93Isolation::~Isolation()
94{
95}
96
97//------------------------------------------------------------------------------
98
99void Isolation::Init()
100{
101 const char *rhoInputArrayName;
102
103 fDeltaRMax = GetDouble("DeltaRMax", 0.5);
104
105 fPTRatioMax = GetDouble("PTRatioMax", 0.1);
106
107 fPTSumMax = GetDouble("PTSumMax", 5.0);
108
109 fUsePTSum = GetBool("UsePTSum", false);
110
111 fClassifier->fPTMin = GetDouble("PTMin", 0.5);
112
113 // import input array(s)
114
115 fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
116 fItIsolationInputArray = fIsolationInputArray->MakeIterator();
117
118 fFilter = new ExRootFilter(fIsolationInputArray);
119
120 fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
121 fItCandidateInputArray = fCandidateInputArray->MakeIterator();
122
123 rhoInputArrayName = GetString("RhoInputArray", "");
124 if(rhoInputArrayName[0] != '\0')
125 {
126 fRhoInputArray = ImportArray(rhoInputArrayName);
127 fItRhoInputArray = fRhoInputArray->MakeIterator();
128 }
129 else
130 {
131 fRhoInputArray = 0;
132 }
133
134 // create output array
135
136 fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
137}
138
139//------------------------------------------------------------------------------
140
141void Isolation::Finish()
142{
143 if(fItRhoInputArray) delete fItRhoInputArray;
144 if(fFilter) delete fFilter;
145 if(fItCandidateInputArray) delete fItCandidateInputArray;
146 if(fItIsolationInputArray) delete fItIsolationInputArray;
147}
148
149//------------------------------------------------------------------------------
150
151void Isolation::Process()
152{
153 Candidate *candidate, *isolation, *object;
154 TObjArray *isolationArray;
155 Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU;
156 Double_t sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr, sum, ratio;
157 Int_t counter;
158 Double_t eta = 0.0;
159 Double_t rho = 0.0;
160
161 // select isolation objects
162 fFilter->Reset();
163 isolationArray = fFilter->GetSubArray(fClassifier, 0);
164
165 if(isolationArray == 0) return;
166
167 TIter itIsolationArray(isolationArray);
168
169 // loop over all input jets
170 fItCandidateInputArray->Reset();
171 while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
172 {
173 const TLorentzVector &candidateMomentum = candidate->Momentum;
174 eta = TMath::Abs(candidateMomentum.Eta());
175
176 // find rho
177 rho = 0.0;
178 if(fRhoInputArray)
179 {
180 fItRhoInputArray->Reset();
181 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
182 {
183 if(eta >= object->Edges[0] && eta < object->Edges[1])
184 {
185 rho = object->Momentum.Pt();
186 }
187 }
188 }
189
190 // loop over all input tracks
191
192 sumNeutral = 0.0;
193 sumCharged = 0.0;
194 sumChargedPU = 0.0;
195 sumAllParticles = 0.0;
196
197 counter = 0;
198 itIsolationArray.Reset();
199
200 while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
201 {
202 const TLorentzVector &isolationMomentum = isolation->Momentum;
203
204 if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
205 candidate->GetUniqueID() != isolation->GetUniqueID())
206 {
207 sumAllParticles += isolationMomentum.Pt();
208 if(isolation->Charge != 0)
209 {
210 if(isolation->IsRecoPU)
211 {
212 sumChargedPU += isolationMomentum.Pt();
213 }
214 else
215 {
216 sumCharged += isolationMomentum.Pt();
217 }
218 }
219 else
220 {
221 sumNeutral += isolationMomentum.Pt();
222 }
223 ++counter;
224 }
225 }
226
227 // find rho
228 rho = 0.0;
229 if(fRhoInputArray)
230 {
231 fItRhoInputArray->Reset();
232 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
233 {
234 if(eta >= object->Edges[0] && eta < object->Edges[1])
235 {
236 rho = object->Momentum.Pt();
237 }
238 }
239 }
240
241 // correct sum for pile-up contamination
242 sumDBeta = sumCharged + TMath::Max(sumNeutral - 0.5*sumChargedPU, 0.0);
243 sumRhoCorr = sumCharged + TMath::Max(sumNeutral - TMath::Max(rho, 0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(), 0.0);
244 ratioDBeta = sumDBeta/candidateMomentum.Pt();
245 ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
246
247 candidate->IsolationVar = ratioDBeta;
248 candidate->IsolationVarRhoCorr = ratioRhoCorr;
249 candidate->SumPtCharged = sumCharged;
250 candidate->SumPtNeutral = sumNeutral;
251 candidate->SumPtChargedPU = sumChargedPU;
252 candidate->SumPt = sumAllParticles;
253
254 sum = fRhoInputArray ? sumRhoCorr : sumDBeta;
255 if(fUsePTSum && sum > fPTSumMax) continue;
256
257 ratio = fRhoInputArray ? ratioRhoCorr : ratioDBeta;
258 if(!fUsePTSum && ratio > fPTRatioMax) continue;
259
260 fOutputArray->Add(candidate);
261 }
262}
263
264//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.