Fork me on GitHub

source: svn/trunk/modules/Isolation.cc

Last change on this file was 1317, checked in by Pavel Demin, 11 years ago

fix loop in JetPileUpSubtractor and add rho eta ranges for Isolation

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