Fork me on GitHub

source: git/modules/Isolation.cc@ 4f728cd

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

add parameter UseRhoCorrection to module Isolation

  • 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 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 fClassifier->fPTMin = GetDouble("PTMin", 0.5);
114
115 // import input array(s)
116
117 fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
118 fItIsolationInputArray = fIsolationInputArray->MakeIterator();
119
120 fFilter = new ExRootFilter(fIsolationInputArray);
121
122 fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
123 fItCandidateInputArray = fCandidateInputArray->MakeIterator();
124
125 rhoInputArrayName = GetString("RhoInputArray", "");
126 if(rhoInputArrayName[0] != '\0')
127 {
128 fRhoInputArray = ImportArray(rhoInputArrayName);
129 fItRhoInputArray = fRhoInputArray->MakeIterator();
130 }
131 else
132 {
133 fRhoInputArray = 0;
134 }
135
136 // create output array
137
138 fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
139}
140
141//------------------------------------------------------------------------------
142
143void Isolation::Finish()
144{
145 if(fItRhoInputArray) delete fItRhoInputArray;
146 if(fFilter) delete fFilter;
147 if(fItCandidateInputArray) delete fItCandidateInputArray;
148 if(fItIsolationInputArray) delete fItIsolationInputArray;
149}
150
151//------------------------------------------------------------------------------
152
153void Isolation::Process()
154{
155 Candidate *candidate, *isolation, *object;
156 TObjArray *isolationArray;
157 Double_t sumChargedNoPU, sumChargedPU, sumNeutral, sumAllParticles;
158 Double_t sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr, sum, ratio;
159 Int_t counter;
160 Double_t eta = 0.0;
161 Double_t rho = 0.0;
162
163 // select isolation objects
164 fFilter->Reset();
165 isolationArray = fFilter->GetSubArray(fClassifier, 0);
166
167 if(isolationArray == 0) return;
168
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 counter = 0;
200 itIsolationArray.Reset();
201
202 while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
203 {
204 const TLorentzVector &isolationMomentum = isolation->Momentum;
205
206 if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
207 candidate->GetUniqueID() != isolation->GetUniqueID())
208 {
209 sumAllParticles += isolationMomentum.Pt();
210 if(isolation->Charge != 0)
211 {
212 if(isolation->IsRecoPU)
213 {
214 sumChargedPU += isolationMomentum.Pt();
215 }
216 else
217 {
218 sumChargedNoPU += isolationMomentum.Pt();
219 }
220 }
221 else
222 {
223 sumNeutral += isolationMomentum.Pt();
224 }
225 ++counter;
226 }
227 }
228
229 // find rho
230 rho = 0.0;
231 if(fRhoInputArray)
232 {
233 fItRhoInputArray->Reset();
234 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
235 {
236 if(eta >= object->Edges[0] && eta < object->Edges[1])
237 {
238 rho = object->Momentum.Pt();
239 }
240 }
241 }
242
243 // correct sum for pile-up contamination
244 sumDBeta = sumChargedNoPU + TMath::Max(sumNeutral - 0.5*sumChargedPU, 0.0);
245 sumRhoCorr = sumChargedNoPU + TMath::Max(sumNeutral - TMath::Max(rho, 0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(), 0.0);
246 ratioDBeta = sumDBeta/candidateMomentum.Pt();
247 ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
248
249 candidate->IsolationVar = ratioDBeta;
250 candidate->IsolationVarRhoCorr = ratioRhoCorr;
251 candidate->SumPtCharged = sumChargedNoPU;
252 candidate->SumPtNeutral = sumNeutral;
253 candidate->SumPtChargedPU = sumChargedPU;
254 candidate->SumPt = sumAllParticles;
255
256 sum = fUseRhoCorrection ? sumRhoCorr : sumDBeta;
257 if(fUsePTSum && sum > fPTSumMax) continue;
258
259 ratio = fUseRhoCorrection ? ratioRhoCorr : ratioDBeta;
260 if(!fUsePTSum && ratio > fPTRatioMax) continue;
261
262 fOutputArray->Add(candidate);
263 }
264}
265
266//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.