Fork me on GitHub

source: git/modules/JetFakeParticle.cc@ 08d82fe

Last change on this file since 08d82fe was 527f67a, checked in by Michele Selvaggi <michele.selvaggi@…>, 4 years ago

fixed formula in JetFakeParticle

  • Property mode set to 100644
File size: 5.5 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/** \class JetFakeParticle
20 *
21 * Converts jet into particle with some PID,
22 * according to parametrized probability.
23 *
24 * \author M. Selvaggi - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "modules/JetFakeParticle.h"
29
30#include "classes/DelphesClasses.h"
31#include "classes/DelphesFactory.h"
32#include "classes/DelphesFormula.h"
33
34#include "ExRootAnalysis/ExRootClassifier.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootResult.h"
37
38#include "TDatabasePDG.h"
39#include "TFormula.h"
40#include "TLorentzVector.h"
41#include "TMath.h"
42#include "TObjArray.h"
43#include "TRandom3.h"
44#include "TString.h"
45
46#include <algorithm>
47#include <iostream>
48#include <sstream>
49#include <stdexcept>
50
51using namespace std;
52
53//------------------------------------------------------------------------------
54
55JetFakeParticle::JetFakeParticle() :
56 fItInputArray(0)
57{
58}
59
60//------------------------------------------------------------------------------
61
62JetFakeParticle::~JetFakeParticle()
63{
64}
65
66//------------------------------------------------------------------------------
67
68void JetFakeParticle::Init()
69{
70 TFakeMap::iterator itEfficiencyMap;
71 ExRootConfParam param;
72 DelphesFormula *formula;
73 Int_t i, size, pdgCode;
74
75 // read efficiency formulas
76 param = GetParam("EfficiencyFormula");
77 size = param.GetSize();
78
79 fEfficiencyMap.clear();
80
81 for(i = 0; i < size / 2; ++i)
82 {
83 formula = new DelphesFormula;
84 formula->Compile(param[i * 2 + 1].GetString());
85 pdgCode = param[i * 2].GetInt();
86
87 if(TMath::Abs(pdgCode) != 11 && TMath::Abs(pdgCode) != 13 && TMath::Abs(pdgCode) != 22)
88 {
89 throw runtime_error("Jets can only fake into electrons, muons or photons. Other particles are not authorized.");
90 }
91
92 fEfficiencyMap[param[i * 2].GetInt()] = formula;
93 }
94
95 // set default efficiency formula
96 itEfficiencyMap = fEfficiencyMap.find(0);
97 if(itEfficiencyMap == fEfficiencyMap.end())
98 {
99 formula = new DelphesFormula;
100 formula->Compile("0.0");
101
102 fEfficiencyMap[0] = formula;
103 }
104
105 // import input array
106
107 fInputArray = ImportArray(GetString("InputArray", "FastJetFinder/jets"));
108 fItInputArray = fInputArray->MakeIterator();
109
110 // create output array
111
112 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "fakeElectrons"));
113 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "fakeMuons"));
114 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "fakePhotons"));
115 fJetOutputArray = ExportArray(GetString("JetOutputArray", "jets"));
116}
117
118//------------------------------------------------------------------------------
119
120void JetFakeParticle::Finish()
121{
122 if(fItInputArray) delete fItInputArray;
123
124 TFakeMap::iterator itEfficiencyMap;
125 DelphesFormula *formula;
126 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
127 {
128 formula = itEfficiencyMap->second;
129 if(formula) delete formula;
130 }
131}
132
133//------------------------------------------------------------------------------
134
135void JetFakeParticle::Process()
136{
137 Candidate *candidate, *fake = 0;
138 Double_t pt, eta, phi, e;
139 TFakeMap::iterator itEfficiencyMap;
140 DelphesFormula *formula;
141 Int_t pdgCodeOut;
142
143 Double_t p, r, rs, total;
144
145 fItInputArray->Reset();
146 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
147 {
148 const TLorentzVector &candidateMomentum = candidate->Momentum;
149 eta = candidateMomentum.Eta();
150 phi = candidateMomentum.Phi();
151 pt = candidateMomentum.Pt();
152 e = candidateMomentum.E();
153
154 r = gRandom->Uniform();
155 total = 0.0;
156 fake = 0;
157
158 // loop over map for this jet
159 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
160 {
161 formula = itEfficiencyMap->second;
162 pdgCodeOut = itEfficiencyMap->first;
163
164 p = formula->Eval(pt, eta, phi, e);
165
166 if(total <= r && r < total + p)
167 {
168 fake = static_cast<Candidate *>(candidate->Clone());
169
170 // convert jet
171
172 if(TMath::Abs(pdgCodeOut) == 11 || TMath::Abs(pdgCodeOut) == 13)
173 {
174 if(candidate->Charge != 0)
175 {
176 fake->Charge = candidate->Charge / TMath::Abs(candidate->Charge);
177 }
178 else
179 {
180 rs = gRandom->Uniform();
181 fake->Charge = (rs < 0.5) ? -1 : 1;
182 }
183 }
184
185 if(TMath::Abs(pdgCodeOut) == 22) fake->PID = 22;
186
187 if(TMath::Abs(pdgCodeOut) == 11) fElectronOutputArray->Add(fake);
188 if(TMath::Abs(pdgCodeOut) == 13) fMuonOutputArray->Add(fake);
189 if(TMath::Abs(pdgCodeOut) == 22) fPhotonOutputArray->Add(fake);
190
191 break;
192 }
193
194 total += p;
195 }
196
197 if(!fake) fJetOutputArray->Add(candidate);
198 }
199}
200
201//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.