Fork me on GitHub

source: git/modules/Isolation.cc@ c24ad79

Last change on this file since c24ad79 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

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