source: HiSusy/trunk/Delphes/Delphes-3.0.9/modules/Isolation.cc @ 5

Last change on this file since 5 was 5, checked in by zerwas, 11 years ago

update to Delphes-3.0.9

File size: 4.6 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-03-06 18:11:35 +0100 (Wed, 06 Mar 2013) $
10 *  $Revision: 1033 $
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{
73  fClassifier = new IsolationClassifier;
74}
75
76//------------------------------------------------------------------------------
77
78Isolation::~Isolation()
79{
80}
81
82//------------------------------------------------------------------------------
83
84void Isolation::Init()
85{
86  const char *rhoInputArrayName;
87
88  fDeltaRMax = GetDouble("DeltaRMax", 0.5);
89
90  fPTRatioMax = GetDouble("PTRatioMax", 0.1);
91
92  fClassifier->fPTMin = GetDouble("PTMin", 0.5);
93
94  // import input array(s)
95
96  fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));
97  fItIsolationInputArray = fIsolationInputArray->MakeIterator();
98
99  fFilter = new ExRootFilter(fIsolationInputArray);
100
101  fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));
102  fItCandidateInputArray = fCandidateInputArray->MakeIterator();
103
104  rhoInputArrayName = GetString("RhoInputArray", "");
105  if(rhoInputArrayName[0] != '\0')
106  {
107    fRhoInputArray = ImportArray(rhoInputArrayName);
108  }
109  else
110  {
111    fRhoInputArray = 0;
112  }
113
114  // create output array
115
116  fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
117}
118
119//------------------------------------------------------------------------------
120
121void Isolation::Finish()
122{
123  if(fFilter) delete fFilter;
124  if(fItCandidateInputArray) delete fItCandidateInputArray;
125  if(fItIsolationInputArray) delete fItIsolationInputArray;
126}
127
128//------------------------------------------------------------------------------
129
130void Isolation::Process()
131{
132  Candidate *candidate, *isolation;
133  TObjArray *isolationArray;
134  Double_t sumPT, ratio;
135  Int_t counter;
136  Double_t rho = 0.0;
137
138  if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)
139  {
140    candidate = static_cast<Candidate*>(fRhoInputArray->At(0));
141    rho = candidate->Momentum.Pt();
142  }
143
144  // select isolation objects
145  fFilter->Reset();
146  isolationArray = fFilter->GetSubArray(fClassifier, 0);
147
148  if(isolationArray == 0) return;
149
150  TIter itIsolationArray(isolationArray);
151
152  // loop over all input jets
153  fItCandidateInputArray->Reset();
154  while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
155  {
156    const TLorentzVector &candidateMomentum = candidate->Momentum;
157
158    // loop over all input tracks
159    sumPT = 0.0;
160    counter = 0;
161    itIsolationArray.Reset();
162    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
163    {
164      const TLorentzVector &isolationMomentum = isolation->Momentum;
165
166      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
167         !candidate->Overlaps(isolation))
168      {
169        sumPT += isolationMomentum.Pt();
170        ++counter;
171      }
172    }
173
174    // correct sumPT for pile-up contamination
175    sumPT = sumPT - rho*fDeltaRMax*fDeltaRMax*TMath::Pi(); 
176
177    ratio = sumPT/candidateMomentum.Pt();
178    if(ratio > fPTRatioMax) continue;
179
180    fOutputArray->Add(candidate);
181  }
182}
183
184//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.