Fork me on GitHub

source: git/modules/BTaggingCMS.cc@ 1d9974e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 1d9974e was 9040259, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

add BTaggingCMS and rename JetFlavourAssociation to JetFlavorAssociation
  • Property mode set to 100644
File size: 13.2 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class BTaggingCMS
21 *
22 * Applies b-tagging efficiency (miss identification rate) formulas
23 * and sets b-tagging flags
24 *
25 * \author P. Demin - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/BTaggingCMS.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56BTaggingCMS::BTaggingCMS() :
57 fItJetInputArray(0)
58{
59}
60
61//------------------------------------------------------------------------------
62
63BTaggingCMS::~BTaggingCMS()
64{
65}
66
67//------------------------------------------------------------------------------
68
69void BTaggingCMS::Init()
70{
71 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
72 ExRootConfParam param;
73 DelphesFormula *formula;
74 Int_t i, size;
75
76 // read efficiency formulas
77 param = GetParam("EfficiencyFormulaLoose");
78 size = param.GetSize();
79
80 fEfficiencyMapLoose.clear();
81 for(i = 0; i < size/2; ++i)
82 {
83 formula = new DelphesFormula;
84 formula->Compile(param[i*2 + 1].GetString());
85 fEfficiencyMapLoose[param[i*2].GetInt()] = formula;
86 }
87
88 // set default efficiency formula
89 itEfficiencyMap = fEfficiencyMapLoose.find(0);
90 if(itEfficiencyMap == fEfficiencyMapLoose.end())
91 {
92 formula = new DelphesFormula;
93 formula->Compile("0.0");
94 fEfficiencyMapLoose[0] = formula;
95 }
96
97 // read efficiency formulas
98 param = GetParam("EfficiencyFormulaMedium");
99 size = param.GetSize();
100
101 fEfficiencyMapMedium.clear();
102 for(i = 0; i < size/2; ++i)
103 {
104 formula = new DelphesFormula;
105 formula->Compile(param[i*2 + 1].GetString());
106 fEfficiencyMapMedium[param[i*2].GetInt()] = formula;
107 }
108
109 // set default efficiency formula
110 itEfficiencyMap = fEfficiencyMapMedium.find(0);
111 if(itEfficiencyMap == fEfficiencyMapMedium.end())
112 {
113 formula = new DelphesFormula;
114 formula->Compile("0.0");
115 fEfficiencyMapMedium[0] = formula;
116 }
117
118 // read efficiency formulas
119 param = GetParam("EfficiencyFormulaTight");
120 size = param.GetSize();
121
122 fEfficiencyMapTight.clear();
123 for(i = 0; i < size/2; ++i)
124 {
125 formula = new DelphesFormula;
126 formula->Compile(param[i*2 + 1].GetString());
127 fEfficiencyMapTight[param[i*2].GetInt()] = formula;
128 }
129
130 // set default efficiency formula
131 itEfficiencyMap = fEfficiencyMapTight.find(0);
132 if(itEfficiencyMap == fEfficiencyMapTight.end())
133 {
134 formula = new DelphesFormula;
135 formula->Compile("0.0");
136 fEfficiencyMapTight[0] = formula;
137 }
138
139 // import input array(s)
140 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
141 fItJetInputArray = fJetInputArray->MakeIterator();
142}
143
144//------------------------------------------------------------------------------
145
146void BTaggingCMS::Finish()
147{
148 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
149 DelphesFormula *formula;
150
151 if(fItJetInputArray) delete fItJetInputArray;
152
153 for(itEfficiencyMap = fEfficiencyMapLoose.begin(); itEfficiencyMap != fEfficiencyMapLoose.end(); ++itEfficiencyMap)
154 {
155 formula = itEfficiencyMap->second;
156 if(formula) delete formula;
157 }
158
159 for(itEfficiencyMap = fEfficiencyMapMedium.begin(); itEfficiencyMap != fEfficiencyMapMedium.end(); ++itEfficiencyMap)
160 {
161 formula = itEfficiencyMap->second;
162 if(formula) delete formula;
163 }
164
165 for(itEfficiencyMap = fEfficiencyMapTight.begin(); itEfficiencyMap != fEfficiencyMapTight.end(); ++itEfficiencyMap)
166 {
167 formula = itEfficiencyMap->second;
168 if(formula) delete formula;
169 }
170}
171
172//------------------------------------------------------------------------------
173
174void BTaggingCMS::Process()
175{
176 Candidate *jet;
177 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
178 DelphesFormula *formula;
179 float randomNumber;
180
181 // loop over all input jets
182 fItJetInputArray->Reset();
183
184 while((jet = static_cast<Candidate*>(fItJetInputArray->Next())))
185 {
186 const TLorentzVector &jetMomentum = jet->Momentum;
187
188 randomNumber = gRandom->Uniform();
189
190 // --------------------------
191 // Loose b-tagging
192 // --------------------------
193
194 // find heaviest flavor and b-tag
195 if(fEfficiencyMapLoose.size() != 0)
196 {
197 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorHeaviest);
198 if(itEfficiencyMap == fEfficiencyMapLoose.end())
199 {
200 itEfficiencyMap = fEfficiencyMapLoose.find(0);
201 }
202 formula = itEfficiencyMap->second;
203
204 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 1;
205
206 // find hisghestpt flavor and b-tag
207 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorHighestPt);
208 if(itEfficiencyMap == fEfficiencyMapLoose.end())
209 {
210 itEfficiencyMap = fEfficiencyMapLoose.find(0);
211 }
212 formula = itEfficiencyMap->second;
213
214 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 1;
215
216 // find nearest2 flavor and b-tag
217 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorNearest2);
218 if(itEfficiencyMap == fEfficiencyMapLoose.end())
219 {
220 itEfficiencyMap = fEfficiencyMapLoose.find(0);
221 }
222 formula = itEfficiencyMap->second;
223
224 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 1;
225
226 // find nearest3 flavor and b-tag
227 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorNearest3);
228 if(itEfficiencyMap == fEfficiencyMapLoose.end())
229 {
230 itEfficiencyMap = fEfficiencyMapLoose.find(0);
231 }
232 formula = itEfficiencyMap->second;
233
234 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 1;
235
236 // find nearest3 flavor and b-tag
237 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorAlgo);
238 if(itEfficiencyMap == fEfficiencyMapLoose.end())
239 {
240 itEfficiencyMap = fEfficiencyMapLoose.find(0);
241 }
242 formula = itEfficiencyMap->second;
243
244 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 1;
245
246 // find nearest3 flavor and b-tag
247 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorPhysics);
248 if(itEfficiencyMap == fEfficiencyMapLoose.end())
249 {
250 itEfficiencyMap = fEfficiencyMapLoose.find(0);
251 }
252 formula = itEfficiencyMap->second;
253
254 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 1;
255
256 // default
257 itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorDefault);
258 if(itEfficiencyMap == fEfficiencyMapLoose.end())
259 {
260 itEfficiencyMap = fEfficiencyMapLoose.find(0);
261 }
262 formula = itEfficiencyMap->second;
263
264 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 1;
265 }
266
267 // --------------------------
268 // Medium b-tagging
269 // --------------------------
270
271 // find heaviest flavor and b-tag
272 if(fEfficiencyMapMedium.size() != 0)
273 {
274 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorHeaviest);
275 if(itEfficiencyMap == fEfficiencyMapMedium.end())
276 {
277 itEfficiencyMap = fEfficiencyMapMedium.find(0);
278 }
279 formula = itEfficiencyMap->second;
280
281 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 2;
282
283 // find hisghestpt flavor and b-tag
284 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorHighestPt);
285 if(itEfficiencyMap == fEfficiencyMapMedium.end())
286 {
287 itEfficiencyMap = fEfficiencyMapMedium.find(0);
288 }
289 formula = itEfficiencyMap->second;
290
291 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 2;
292
293 // find nearest2 flavor and b-tag
294 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorNearest2);
295 if(itEfficiencyMap == fEfficiencyMapMedium.end())
296 {
297 itEfficiencyMap = fEfficiencyMapMedium.find(0);
298 }
299 formula = itEfficiencyMap->second;
300
301 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 2;
302
303 // find nearest3 flavor and b-tag
304 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorNearest3);
305 if(itEfficiencyMap == fEfficiencyMapMedium.end())
306 {
307 itEfficiencyMap = fEfficiencyMapMedium.find(0);
308 }
309 formula = itEfficiencyMap->second;
310
311 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 2;
312
313 // find nearest3 flavor and b-tag
314 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorAlgo);
315 if(itEfficiencyMap == fEfficiencyMapMedium.end())
316 {
317 itEfficiencyMap = fEfficiencyMapMedium.find(0);
318 }
319 formula = itEfficiencyMap->second;
320
321 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 2;
322
323 // find nearest3 flavor and b-tag
324 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorPhysics);
325 if(itEfficiencyMap == fEfficiencyMapMedium.end())
326 {
327 itEfficiencyMap = fEfficiencyMapMedium.find(0);
328 }
329 formula = itEfficiencyMap->second;
330
331 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 2;
332
333 // default
334 itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorDefault);
335 if(itEfficiencyMap == fEfficiencyMapMedium.end())
336 {
337 itEfficiencyMap = fEfficiencyMapMedium.find(0);
338 }
339 formula = itEfficiencyMap->second;
340
341 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 2;
342 }
343
344 // --------------------------
345 // Tight b-tagging
346 // --------------------------
347
348 // find heaviest flavor and b-tag
349 if(fEfficiencyMapTight.size() != 0)
350 {
351 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorHeaviest);
352 if(itEfficiencyMap == fEfficiencyMapTight.end())
353 {
354 itEfficiencyMap = fEfficiencyMapTight.find(0);
355 }
356 formula = itEfficiencyMap->second;
357
358 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 3;
359
360 // find hisghestpt flavor and b-tag
361 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorHighestPt);
362 if(itEfficiencyMap == fEfficiencyMapTight.end())
363 {
364 itEfficiencyMap = fEfficiencyMapTight.find(0);
365 }
366 formula = itEfficiencyMap->second;
367
368 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 3;
369
370 // find nearest2 flavor and b-tag
371 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorNearest2);
372 if(itEfficiencyMap == fEfficiencyMapTight.end())
373 {
374 itEfficiencyMap = fEfficiencyMapTight.find(0);
375 }
376 formula = itEfficiencyMap->second;
377
378 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 3;
379
380 // find nearest3 flavor and b-tag
381 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorNearest3);
382 if(itEfficiencyMap == fEfficiencyMapTight.end())
383 {
384 itEfficiencyMap = fEfficiencyMapTight.find(0);
385 }
386 formula = itEfficiencyMap->second;
387
388 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 3;
389
390 // find nearest3 flavor and b-tag
391 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorAlgo);
392 if(itEfficiencyMap == fEfficiencyMapTight.end())
393 {
394 itEfficiencyMap = fEfficiencyMapTight.find(0);
395 }
396 formula = itEfficiencyMap->second;
397
398 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 3;
399
400 // find nearest3 flavor and b-tag
401 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorPhysics);
402 if(itEfficiencyMap == fEfficiencyMapTight.end())
403 {
404 itEfficiencyMap = fEfficiencyMapTight.find(0);
405 }
406 formula = itEfficiencyMap->second;
407
408 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 3;
409
410 // default
411 itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorDefault);
412 if(itEfficiencyMap == fEfficiencyMapTight.end())
413 {
414 itEfficiencyMap = fEfficiencyMapTight.find(0);
415 }
416 formula = itEfficiencyMap->second;
417
418 if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 3;
419 }
420 }
421}
Note: See TracBrowser for help on using the repository browser.