Fork me on GitHub

source: svn/trunk/modules/TauTagging.cc

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

check tau's daughter index

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 6.5 KB
Line 
1
2/** \class TauTagging
3 *
4 * Determines origin of jet,
5 * applies b-tagging efficiency (miss identification rate) formulas
6 * and sets b-tagging flags
7 *
8 * $Date: 2013-07-11 22:22:27 +0000 (Thu, 11 Jul 2013) $
9 * $Revision: 1217 $
10 *
11 *
12 * \author P. Demin - UCL, Louvain-la-Neuve
13 *
14 */
15
16#include "modules/TauTagging.h"
17
18#include "classes/DelphesClasses.h"
19#include "classes/DelphesFactory.h"
20#include "classes/DelphesFormula.h"
21
22#include "ExRootAnalysis/ExRootResult.h"
23#include "ExRootAnalysis/ExRootFilter.h"
24#include "ExRootAnalysis/ExRootClassifier.h"
25
26#include "TMath.h"
27#include "TString.h"
28#include "TFormula.h"
29#include "TRandom3.h"
30#include "TObjArray.h"
31#include "TDatabasePDG.h"
32#include "TLorentzVector.h"
33
34#include <algorithm>
35#include <stdexcept>
36#include <iostream>
37#include <sstream>
38
39using namespace std;
40
41//------------------------------------------------------------------------------
42
43class TauTaggingPartonClassifier : public ExRootClassifier
44{
45public:
46
47 TauTaggingPartonClassifier(const TObjArray *array);
48
49 Int_t GetCategory(TObject *object);
50
51 Double_t fEtaMax, fPTMin;
52
53 const TObjArray *fParticleInputArray;
54};
55
56//------------------------------------------------------------------------------
57TauTaggingPartonClassifier::TauTaggingPartonClassifier(const TObjArray *array) :
58 fParticleInputArray(array)
59{
60}
61
62//------------------------------------------------------------------------------
63
64Int_t TauTaggingPartonClassifier::GetCategory(TObject *object)
65{
66 Candidate *tau = static_cast<Candidate *>(object);
67 Candidate *daughter = 0;
68 const TLorentzVector &momentum = tau->Momentum;
69 Int_t pdgCode, i;
70
71 pdgCode = TMath::Abs(tau->PID);
72 if(pdgCode != 15) return -1;
73
74 if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
75
76 if(tau->D1 < 0) return -1;
77
78 if(tau->D1 >= fParticleInputArray->GetEntriesFast() ||
79 tau->D2 >= fParticleInputArray->GetEntriesFast())
80 {
81 throw runtime_error("tau's daughter index is greater than the ParticleInputArray size");
82 }
83
84 for(i = tau->D1; i <= tau->D2; ++i)
85 {
86 daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
87 pdgCode = TMath::Abs(daughter->PID);
88 if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15 || pdgCode == 24) return -1;
89 }
90
91 return 0;
92}
93
94//------------------------------------------------------------------------------
95
96TauTagging::TauTagging() :
97 fClassifier(0), fFilter(0),
98 fItPartonInputArray(0), fItJetInputArray(0)
99{
100}
101
102//------------------------------------------------------------------------------
103
104TauTagging::~TauTagging()
105{
106}
107
108//------------------------------------------------------------------------------
109
110void TauTagging::Init()
111{
112 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
113 ExRootConfParam param;
114 DelphesFormula *formula;
115 Int_t i, size;
116
117 fDeltaR = GetDouble("DeltaR", 0.5);
118
119 // read efficiency formulas
120 param = GetParam("EfficiencyFormula");
121 size = param.GetSize();
122
123 fEfficiencyMap.clear();
124 for(i = 0; i < size/2; ++i)
125 {
126 formula = new DelphesFormula;
127 formula->Compile(param[i*2 + 1].GetString());
128
129 fEfficiencyMap[param[i*2].GetInt()] = formula;
130 }
131
132 // set default efficiency formula
133 itEfficiencyMap = fEfficiencyMap.find(0);
134 if(itEfficiencyMap == fEfficiencyMap.end())
135 {
136 formula = new DelphesFormula;
137 formula->Compile("0.0");
138
139 fEfficiencyMap[0] = formula;
140 }
141
142 // import input array(s)
143
144 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
145
146 fClassifier = new TauTaggingPartonClassifier(fParticleInputArray);
147 fClassifier->fPTMin = GetDouble("TauPTMin", 1.0);
148 fClassifier->fEtaMax = GetDouble("TauEtaMax", 2.5);
149
150 fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
151 fItPartonInputArray = fPartonInputArray->MakeIterator();
152
153 fFilter = new ExRootFilter(fPartonInputArray);
154
155 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
156 fItJetInputArray = fJetInputArray->MakeIterator();
157}
158
159//------------------------------------------------------------------------------
160
161void TauTagging::Finish()
162{
163 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
164 DelphesFormula *formula;
165
166 if(fFilter) delete fFilter;
167 if(fClassifier) delete fClassifier;
168 if(fItJetInputArray) delete fItJetInputArray;
169 if(fItPartonInputArray) delete fItPartonInputArray;
170
171 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
172 {
173 formula = itEfficiencyMap->second;
174 if(formula) delete formula;
175 }
176}
177
178//------------------------------------------------------------------------------
179
180void TauTagging::Process()
181{
182 Candidate *jet, *tau, *daughter;
183 TLorentzVector tauMomentum;
184 Double_t pt, eta, phi;
185 TObjArray *tauArray;
186 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
187 DelphesFormula *formula;
188 Int_t pdgCode, charge, i;
189
190 // select taus
191 fFilter->Reset();
192 tauArray = fFilter->GetSubArray(fClassifier, 0);
193
194 if(tauArray == 0) return;
195
196 TIter itTauArray(tauArray);
197
198 // loop over all input jets
199 fItJetInputArray->Reset();
200 while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
201 {
202 const TLorentzVector &jetMomentum = jet->Momentum;
203 pdgCode = 0;
204 charge = gRandom->Uniform() > 0.5 ? 1 : -1;
205 eta = jetMomentum.Eta();
206 phi = jetMomentum.Phi();
207 pt = jetMomentum.Pt();
208
209 // loop over all input taus
210 itTauArray.Reset();
211 while((tau = static_cast<Candidate *>(itTauArray.Next())))
212 {
213 if(tau->D1 < 0) continue;
214
215 if(tau->D1 >= fParticleInputArray->GetEntriesFast() ||
216 tau->D2 >= fParticleInputArray->GetEntriesFast())
217 {
218 throw runtime_error("tau's daughter index is greater than the ParticleInputArray size");
219 }
220
221 tauMomentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
222
223 for(i = tau->D1; i <= tau->D2; ++i)
224 {
225 daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
226 if(TMath::Abs(daughter->PID) == 16) continue;
227 tauMomentum += daughter->Momentum;
228 }
229
230 if(jetMomentum.DeltaR(tauMomentum) <= fDeltaR)
231 {
232 pdgCode = 15;
233 charge = tau->Charge;
234 }
235 }
236 // find an efficency formula
237 itEfficiencyMap = fEfficiencyMap.find(pdgCode);
238 if(itEfficiencyMap == fEfficiencyMap.end())
239 {
240 itEfficiencyMap = fEfficiencyMap.find(0);
241 }
242 formula = itEfficiencyMap->second;
243
244 // apply an efficency formula
245 jet->TauTag = gRandom->Uniform() <= formula->Eval(pt, eta);
246 // set tau charge
247 jet->Charge = charge;
248 }
249}
250
251//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.