Fork me on GitHub

source: svn/trunk/modules/Isolation.cc@ 1364

Last change on this file since 1364 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
RevLine 
[694]1
[814]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
[1028]7 * the transverse momenta fraction within (PTRatioMin, PTRatioMax].
[814]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
[694]17#include "modules/Isolation.h"
18
19#include "classes/DelphesClasses.h"
[703]20#include "classes/DelphesFactory.h"
[766]21#include "classes/DelphesFormula.h"
[694]22
[703]23#include "ExRootAnalysis/ExRootResult.h"
[694]24#include "ExRootAnalysis/ExRootFilter.h"
25#include "ExRootAnalysis/ExRootClassifier.h"
26
27#include "TMath.h"
28#include "TString.h"
29#include "TFormula.h"
[703]30#include "TRandom3.h"
31#include "TObjArray.h"
32#include "TDatabasePDG.h"
[694]33#include "TLorentzVector.h"
34
[1028]35#include <algorithm>
[703]36#include <stdexcept>
[694]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);
[1028]51
[694]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),
[1317]71 fItIsolationInputArray(0), fItCandidateInputArray(0),
72 fItRhoInputArray(0)
[694]73{
74 fClassifier = new IsolationClassifier;
75}
76
77//------------------------------------------------------------------------------
78
79Isolation::~Isolation()
80{
81}
82
83//------------------------------------------------------------------------------
84
85void Isolation::Init()
86{
[1023]87 const char *rhoInputArrayName;
88
[878]89 fDeltaRMax = GetDouble("DeltaRMax", 0.5);
[694]90
[878]91 fPTRatioMax = GetDouble("PTRatioMax", 0.1);
[696]92
[1267]93 fPTSumMax = GetDouble("PTSumMax", 5.0);
94
95 fUsePTSum = GetBool("UsePTSum", false);
96
[694]97 fClassifier->fPTMin = GetDouble("PTMin", 0.5);
[1028]98
[694]99 // import input array(s)
100
101 fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
102 fItIsolationInputArray = fIsolationInputArray->MakeIterator();
103
104 fFilter = new ExRootFilter(fIsolationInputArray);
[1028]105
[694]106 fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
107 fItCandidateInputArray = fCandidateInputArray->MakeIterator();
[1028]108
[1023]109 rhoInputArrayName = GetString("RhoInputArray", "");
110 if(rhoInputArrayName[0] != '\0')
111 {
[1033]112 fRhoInputArray = ImportArray(rhoInputArrayName);
[1317]113 fItRhoInputArray = fRhoInputArray->MakeIterator();
[1023]114 }
115 else
116 {
117 fRhoInputArray = 0;
118 }
[1008]119
[694]120 // create output array
121
122 fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
123}
124
125//------------------------------------------------------------------------------
126
127void Isolation::Finish()
128{
[1317]129 if(fItRhoInputArray) delete fItRhoInputArray;
[694]130 if(fFilter) delete fFilter;
131 if(fItCandidateInputArray) delete fItCandidateInputArray;
132 if(fItIsolationInputArray) delete fItIsolationInputArray;
133}
134
135//------------------------------------------------------------------------------
136
137void Isolation::Process()
138{
[1317]139 Candidate *candidate, *isolation, *object;
[694]140 TObjArray *isolationArray;
[1267]141 Double_t sum, ratio;
[894]142 Int_t counter;
[1317]143 Double_t eta = 0.0;
[1028]144 Double_t rho = 0.0;
145
146 if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)
[1023]147 {
148 candidate = static_cast<Candidate*>(fRhoInputArray->At(0));
[1008]149 rho = candidate->Momentum.Pt();
150 }
[1028]151
[694]152 // select isolation objects
153 fFilter->Reset();
154 isolationArray = fFilter->GetSubArray(fClassifier, 0);
[1028]155
[694]156 if(isolationArray == 0) return;
157
158 TIter itIsolationArray(isolationArray);
[1028]159
[694]160 // loop over all input jets
161 fItCandidateInputArray->Reset();
162 while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
163 {
[753]164 const TLorentzVector &candidateMomentum = candidate->Momentum;
[1317]165 eta = TMath::Abs(candidateMomentum.Eta());
[694]166
167 // loop over all input tracks
[1267]168 sum = 0.0;
[894]169 counter = 0;
[694]170 itIsolationArray.Reset();
171 while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
172 {
173 const TLorentzVector &isolationMomentum = isolation->Momentum;
[894]174
175 if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
176 !candidate->Overlaps(isolation))
[694]177 {
[1267]178 sum += isolationMomentum.Pt();
[894]179 ++counter;
[694]180 }
181 }
[878]182
[1317]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
[1267]197 // correct sum for pile-up contamination
[1317]198 sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
[1028]199
[1267]200 ratio = sum/candidateMomentum.Pt();
201 if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
[753]202
203 fOutputArray->Add(candidate);
[694]204 }
205}
206
207//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.