source: HiSusy/trunk/Delphes-3.0.0/modules/BTagging.cc @ 1

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

first import of structure, PYTHIA8 and DELPHES

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: 2012-12-06 01:10:00 +0100 (Thu, 06 Dec 2012) $
9 *  $Revision: 870 $
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.