Fork me on GitHub

source: git/modules/Isolation.cc@ 7834e65f

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 7834e65f was 65dc4e7, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

clean iso code

  • Property mode set to 100644
File size: 7.5 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 fUseRhoCorrection = GetBool("UseRhoCorrection", true);
112
113 fDeltaRMin = GetDouble("DeltaRMin", 0.01);
114 fUseMiniCone = GetBool("UseMiniCone", false);
115
116 fClassifier->fPTMin = GetDouble("PTMin", 0.5);
117
118 // import input array(s)
119
120 fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
121 fItIsolationInputArray = fIsolationInputArray->MakeIterator();
122
123 fFilter = new ExRootFilter(fIsolationInputArray);
124
125 fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
126 fItCandidateInputArray = fCandidateInputArray->MakeIterator();
127
128 rhoInputArrayName = GetString("RhoInputArray", "");
129 if(rhoInputArrayName[0] != '\0')
130 {
131 fRhoInputArray = ImportArray(rhoInputArrayName);
132 fItRhoInputArray = fRhoInputArray->MakeIterator();
133 }
134 else
135 {
136 fRhoInputArray = 0;
137 }
138
139 // create output array
140
141 fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
142}
143
144//------------------------------------------------------------------------------
145
146void Isolation::Finish()
147{
148 if(fItRhoInputArray) delete fItRhoInputArray;
149 if(fFilter) delete fFilter;
150 if(fItCandidateInputArray) delete fItCandidateInputArray;
151 if(fItIsolationInputArray) delete fItIsolationInputArray;
152}
153
154//------------------------------------------------------------------------------
155
156void Isolation::Process()
157{
158 Candidate *candidate, *isolation, *object;
159 TObjArray *isolationArray;
160 Double_t sumChargedNoPU, sumChargedPU, sumNeutral, sumAllParticles;
161 Double_t sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr, sum, ratio;
162 Bool_t pass = kFALSE;
163 Double_t eta = 0.0;
164 Double_t rho = 0.0;
165
166 // select isolation objects
167 fFilter->Reset();
168 isolationArray = fFilter->GetSubArray(fClassifier, 0);
169 TIter itIsolationArray(isolationArray);
170
171 // loop over all input jets
172 fItCandidateInputArray->Reset();
173 while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
174 {
175 const TLorentzVector &candidateMomentum = candidate->Momentum;
176 eta = TMath::Abs(candidateMomentum.Eta());
177
178 // find rho
179 rho = 0.0;
180 if(fRhoInputArray)
181 {
182 fItRhoInputArray->Reset();
183 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
184 {
185 if(eta >= object->Edges[0] && eta < object->Edges[1])
186 {
187 rho = object->Momentum.Pt();
188 }
189 }
190 }
191
192 // loop over all input tracks
193
194 sumNeutral = 0.0;
195 sumChargedNoPU = 0.0;
196 sumChargedPU = 0.0;
197 sumAllParticles = 0.0;
198
199 itIsolationArray.Reset();
200 while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
201 {
202 const TLorentzVector &isolationMomentum = isolation->Momentum;
203
204 if(fUseMiniCone)
205 {
206 pass = candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
207 candidateMomentum.DeltaR(isolationMomentum) > fDeltaRMin;
208 }
209 else
210 {
211 pass = candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
212 candidate->GetUniqueID() != isolation->GetUniqueID();
213 }
214
215 if(pass)
216 {
217
218 sumAllParticles += isolationMomentum.Pt();
219 if(isolation->Charge != 0)
220 {
221 if(isolation->IsRecoPU)
222 {
223 sumChargedPU += isolationMomentum.Pt();
224 }
225 else
226 {
227 sumChargedNoPU += isolationMomentum.Pt();
228 }
229 }
230 else
231 {
232 sumNeutral += isolationMomentum.Pt();
233 }
234 }
235
236 }
237
238 // find rho
239 rho = 0.0;
240 if(fRhoInputArray)
241 {
242 fItRhoInputArray->Reset();
243 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
244 {
245 if(eta >= object->Edges[0] && eta < object->Edges[1])
246 {
247 rho = object->Momentum.Pt();
248 }
249 }
250 }
251
252
253
254 // correct sum for pile-up contamination
255 sumDBeta = sumChargedNoPU + TMath::Max(sumNeutral - 0.5*sumChargedPU, 0.0);
256 sumRhoCorr = sumChargedNoPU + TMath::Max(sumNeutral - TMath::Max(rho, 0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(), 0.0);
257 ratioDBeta = sumDBeta/candidateMomentum.Pt();
258 ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
259
260 candidate->IsolationVar = ratioDBeta;
261 candidate->IsolationVarRhoCorr = ratioRhoCorr;
262 candidate->SumPtCharged = sumChargedNoPU;
263 candidate->SumPtNeutral = sumNeutral;
264 candidate->SumPtChargedPU = sumChargedPU;
265 candidate->SumPt = sumAllParticles;
266
267 sum = fUseRhoCorrection ? sumRhoCorr : sumDBeta;
268 if(fUsePTSum && sum > fPTSumMax) continue;
269
270 ratio = fUseRhoCorrection ? ratioRhoCorr : ratioDBeta;
271 if(!fUsePTSum && ratio > fPTRatioMax) continue;
272
273 fOutputArray->Add(candidate);
274 }
275}
276
277//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.