Fork me on GitHub

source: git/modules/BTagging.cc@ 00bcbe6

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 00bcbe6 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 5.2 KB
Line 
1
2/** \class BTagging
3 *
4 * Determines origin of jet,
5 * applies b-tagging efficiency (miss identification rate) formulas
6 * and sets b-tagging flags
7 *
8 * $Date$
9 * $Revision$
10 *
11 *
12 * \author P. Demin - UCL, Louvain-la-Neuve
13 *
14 */
15
16#include "modules/BTagging.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 BTaggingPartonClassifier : public ExRootClassifier
44{
45public:
46
47 BTaggingPartonClassifier() {}
48
49 Int_t GetCategory(TObject *object);
50
51 Double_t fEtaMax, fPTMin;
52};
53
54//------------------------------------------------------------------------------
55
56Int_t BTaggingPartonClassifier::GetCategory(TObject *object)
57{
58 Candidate *parton = static_cast<Candidate*>(object);
59 const TLorentzVector &momentum = parton->Momentum;
60 Int_t pdgCode;
61
62 if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
63
64 pdgCode = TMath::Abs(parton->PID);
65 if(pdgCode != 21 && pdgCode > 5) return -1;
66
67 return 0;
68}
69
70//------------------------------------------------------------------------------
71
72BTagging::BTagging() :
73 fClassifier(0), fFilter(0),
74 fItPartonInputArray(0), fItJetInputArray(0)
75{
76 fClassifier = new BTaggingPartonClassifier;
77}
78
79//------------------------------------------------------------------------------
80
81BTagging::~BTagging()
82{
83 if(fClassifier) delete fClassifier;
84}
85
86//------------------------------------------------------------------------------
87
88void BTagging::Init()
89{
90 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
91 ExRootConfParam param;
92 DelphesFormula *formula;
93 Int_t i, size;
94
95 fDeltaR = GetDouble("DeltaR", 0.5);
96
97 fClassifier->fPTMin = GetDouble("PartonPTMin", 1.0);
98 fClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
99
100 // read efficiency formulas
101 param = GetParam("EfficiencyFormula");
102 size = param.GetSize();
103
104 fEfficiencyMap.clear();
105 for(i = 0; i < size/2; ++i)
106 {
107 formula = new DelphesFormula;
108 formula->Compile(param[i*2 + 1].GetString());
109
110 fEfficiencyMap[param[i*2].GetInt()] = formula;
111 }
112
113 // set default efficiency formula
114 itEfficiencyMap = fEfficiencyMap.find(0);
115 if(itEfficiencyMap == fEfficiencyMap.end())
116 {
117 formula = new DelphesFormula;
118 formula->Compile("0.0");
119
120 fEfficiencyMap[0] = formula;
121 }
122
123 // import input array(s)
124
125 fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
126 fItPartonInputArray = fPartonInputArray->MakeIterator();
127
128 fFilter = new ExRootFilter(fPartonInputArray);
129
130 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
131 fItJetInputArray = fJetInputArray->MakeIterator();
132}
133
134//------------------------------------------------------------------------------
135
136void BTagging::Finish()
137{
138 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
139 DelphesFormula *formula;
140
141 if(fFilter) delete fFilter;
142 if(fItJetInputArray) delete fItJetInputArray;
143 if(fItPartonInputArray) delete fItPartonInputArray;
144
145 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
146 {
147 formula = itEfficiencyMap->second;
148 if(formula) delete formula;
149 }
150}
151
152//------------------------------------------------------------------------------
153
154void BTagging::Process()
155{
156 Candidate *jet, *parton;
157 Double_t pt, eta, phi;
158 TObjArray *partonArray;
159 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
160 DelphesFormula *formula;
161 Int_t pdgCode, pdgCodeMax;
162
163 // select quark and gluons
164 fFilter->Reset();
165 partonArray = fFilter->GetSubArray(fClassifier, 0);
166
167 if(partonArray == 0) return;
168
169 TIter itPartonArray(partonArray);
170
171 // loop over all input jets
172 fItJetInputArray->Reset();
173 while((jet = static_cast<Candidate*>(fItJetInputArray->Next())))
174 {
175 const TLorentzVector &jetMomentum = jet->Momentum;
176 pdgCodeMax = -1;
177 eta = jetMomentum.Eta();
178 phi = jetMomentum.Phi();
179 pt = jetMomentum.Pt();
180
181 // loop over all input partons
182 itPartonArray.Reset();
183 while((parton = static_cast<Candidate*>(itPartonArray.Next())))
184 {
185 pdgCode = TMath::Abs(parton->PID);
186 if(pdgCode == 21) pdgCode = 0;
187 if(jetMomentum.DeltaR(parton->Momentum) <= fDeltaR)
188 {
189 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
190 }
191 }
192 if(pdgCodeMax == 0) pdgCodeMax = 21;
193 if(pdgCodeMax == -1) pdgCodeMax = 0;
194
195 // find an efficency formula
196 itEfficiencyMap = fEfficiencyMap.find(pdgCodeMax);
197 if(itEfficiencyMap == fEfficiencyMap.end())
198 {
199 itEfficiencyMap = fEfficiencyMap.find(0);
200 }
201 formula = itEfficiencyMap->second;
202
203 // apply an efficency formula
204 jet->BTag = gRandom->Uniform() <= formula->Eval(pt, eta);
205 }
206}
207
208//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.