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 |
|
---|
39 | using namespace std;
|
---|
40 |
|
---|
41 | //------------------------------------------------------------------------------
|
---|
42 |
|
---|
43 | class BTaggingPartonClassifier : public ExRootClassifier
|
---|
44 | {
|
---|
45 | public:
|
---|
46 |
|
---|
47 | BTaggingPartonClassifier() {}
|
---|
48 |
|
---|
49 | Int_t GetCategory(TObject *object);
|
---|
50 |
|
---|
51 | Double_t fEtaMax, fPTMin;
|
---|
52 | };
|
---|
53 |
|
---|
54 | //------------------------------------------------------------------------------
|
---|
55 |
|
---|
56 | Int_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 |
|
---|
72 | BTagging::BTagging() :
|
---|
73 | fClassifier(0), fFilter(0),
|
---|
74 | fItPartonInputArray(0), fItJetInputArray(0)
|
---|
75 | {
|
---|
76 | fClassifier = new BTaggingPartonClassifier;
|
---|
77 | }
|
---|
78 |
|
---|
79 | //------------------------------------------------------------------------------
|
---|
80 |
|
---|
81 | BTagging::~BTagging()
|
---|
82 | {
|
---|
83 | if(fClassifier) delete fClassifier;
|
---|
84 | }
|
---|
85 |
|
---|
86 | //------------------------------------------------------------------------------
|
---|
87 |
|
---|
88 | void 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 |
|
---|
136 | void 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 |
|
---|
154 | void 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 | //------------------------------------------------------------------------------
|
---|