Fork me on GitHub

source: git/modules/Isolation.cc@ 0f8390e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0f8390e was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix GPLv3 header

  • Property mode set to 100644
File size: 5.9 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 * $Date$
28 * $Revision$
29 *
30 *
31 * \author P. Demin - UCL, Louvain-la-Neuve
32 *
33 */
34
35#include "modules/Isolation.h"
36
37#include "classes/DelphesClasses.h"
38#include "classes/DelphesFactory.h"
39#include "classes/DelphesFormula.h"
40
41#include "ExRootAnalysis/ExRootResult.h"
42#include "ExRootAnalysis/ExRootFilter.h"
43#include "ExRootAnalysis/ExRootClassifier.h"
44
45#include "TMath.h"
46#include "TString.h"
47#include "TFormula.h"
48#include "TRandom3.h"
49#include "TObjArray.h"
50#include "TDatabasePDG.h"
51#include "TLorentzVector.h"
52
53#include <algorithm>
54#include <stdexcept>
55#include <iostream>
56#include <sstream>
57
58using namespace std;
59
60//------------------------------------------------------------------------------
61
62class IsolationClassifier : public ExRootClassifier
63{
64public:
65
66 IsolationClassifier() {}
67
68 Int_t GetCategory(TObject *object);
69
70 Double_t fPTMin;
71};
72
73//------------------------------------------------------------------------------
74
75Int_t IsolationClassifier::GetCategory(TObject *object)
76{
77 Candidate *track = static_cast<Candidate*>(object);
78 const TLorentzVector &momentum = track->Momentum;
79
80 if(momentum.Pt() < fPTMin) return -1;
81
82 return 0;
83}
84
85//------------------------------------------------------------------------------
86
87Isolation::Isolation() :
88 fClassifier(0), fFilter(0),
89 fItIsolationInputArray(0), fItCandidateInputArray(0),
90 fItRhoInputArray(0)
91{
92 fClassifier = new IsolationClassifier;
93}
94
95//------------------------------------------------------------------------------
96
97Isolation::~Isolation()
98{
99}
100
101//------------------------------------------------------------------------------
102
103void Isolation::Init()
104{
105 const char *rhoInputArrayName;
106
107 fDeltaRMax = GetDouble("DeltaRMax", 0.5);
108
109 fPTRatioMax = GetDouble("PTRatioMax", 0.1);
110
111 fPTSumMax = GetDouble("PTSumMax", 5.0);
112
113 fUsePTSum = GetBool("UsePTSum", false);
114
115 fClassifier->fPTMin = GetDouble("PTMin", 0.5);
116
117 // import input array(s)
118
119 fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
120 fItIsolationInputArray = fIsolationInputArray->MakeIterator();
121
122 fFilter = new ExRootFilter(fIsolationInputArray);
123
124 fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
125 fItCandidateInputArray = fCandidateInputArray->MakeIterator();
126
127 rhoInputArrayName = GetString("RhoInputArray", "");
128 if(rhoInputArrayName[0] != '\0')
129 {
130 fRhoInputArray = ImportArray(rhoInputArrayName);
131 fItRhoInputArray = fRhoInputArray->MakeIterator();
132 }
133 else
134 {
135 fRhoInputArray = 0;
136 }
137
138 // create output array
139
140 fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
141}
142
143//------------------------------------------------------------------------------
144
145void Isolation::Finish()
146{
147 if(fItRhoInputArray) delete fItRhoInputArray;
148 if(fFilter) delete fFilter;
149 if(fItCandidateInputArray) delete fItCandidateInputArray;
150 if(fItIsolationInputArray) delete fItIsolationInputArray;
151}
152
153//------------------------------------------------------------------------------
154
155void Isolation::Process()
156{
157 Candidate *candidate, *isolation, *object;
158 TObjArray *isolationArray;
159 Double_t sum, ratio;
160 Int_t counter;
161 Double_t eta = 0.0;
162 Double_t rho = 0.0;
163
164 if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)
165 {
166 candidate = static_cast<Candidate*>(fRhoInputArray->At(0));
167 rho = candidate->Momentum.Pt();
168 }
169
170 // select isolation objects
171 fFilter->Reset();
172 isolationArray = fFilter->GetSubArray(fClassifier, 0);
173
174 if(isolationArray == 0) return;
175
176 TIter itIsolationArray(isolationArray);
177
178 // loop over all input jets
179 fItCandidateInputArray->Reset();
180 while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
181 {
182 const TLorentzVector &candidateMomentum = candidate->Momentum;
183 eta = TMath::Abs(candidateMomentum.Eta());
184
185 // loop over all input tracks
186 sum = 0.0;
187 counter = 0;
188 itIsolationArray.Reset();
189 while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
190 {
191 const TLorentzVector &isolationMomentum = isolation->Momentum;
192
193 if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
194 !candidate->Overlaps(isolation))
195 {
196 sum += isolationMomentum.Pt();
197 ++counter;
198 }
199 }
200
201 // find rho
202 rho = 0.0;
203 if(fRhoInputArray)
204 {
205 fItRhoInputArray->Reset();
206 while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
207 {
208 if(eta >= object->Edges[0] && eta < object->Edges[1])
209 {
210 rho = object->Momentum.Pt();
211 }
212 }
213 }
214
215 // correct sum for pile-up contamination
216 sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
217
218 ratio = sum/candidateMomentum.Pt();
219 if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
220
221 fOutputArray->Add(candidate);
222 }
223}
224
225//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.