Fork me on GitHub

source: git/modules/IdentificationMap.cc@ 835f0d5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 835f0d5 was 57e42c6, checked in by Michele <michele.selvaggi@…>, 10 years ago

minor change

  • Property mode set to 100644
File size: 5.4 KB
Line 
1
2/** \class IdentificationMap
3 *
4 * Converts particles with some PDG code into another particle, according to parametrized probability as fction of pt eta
5 given by user.
6 *
7 * $Date: 2014-08-07 14:57:44 +0100 (Thu, 07 Aug 2014) $
8 * $Revision: 905 $
9 *
10 *
11 * \author M. Selvaggi - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "modules/IdentificationMap.h"
16
17#include "classes/DelphesClasses.h"
18#include "classes/DelphesFactory.h"
19#include "classes/DelphesFormula.h"
20
21#include "ExRootAnalysis/ExRootResult.h"
22#include "ExRootAnalysis/ExRootFilter.h"
23#include "ExRootAnalysis/ExRootClassifier.h"
24
25#include "TMath.h"
26#include "TString.h"
27#include "TFormula.h"
28#include "TRandom3.h"
29#include "TObjArray.h"
30#include "TDatabasePDG.h"
31#include "TLorentzVector.h"
32
33#include <algorithm>
34#include <stdexcept>
35#include <iostream>
36#include <sstream>
37
38using namespace std;
39
40//------------------------------------------------------------------------------
41
42IdentificationMap::IdentificationMap() :
43 fItInputArray(0)
44{
45}
46
47//------------------------------------------------------------------------------
48
49IdentificationMap::~IdentificationMap()
50{
51}
52
53//------------------------------------------------------------------------------
54
55void IdentificationMap::Init()
56{
57 // read efficiency formula
58
59
60 TMisIDMap::iterator itEfficiencyMap;
61 ExRootConfParam param;
62 DelphesFormula *formula;
63 Int_t i, size, pdg;
64
65 // read efficiency formulas
66 param = GetParam("EfficiencyFormula");
67 size = param.GetSize();
68
69 fEfficiencyMap.clear();
70 for(i = 0; i < size/3; ++i)
71 {
72 formula = new DelphesFormula;
73 formula->Compile(param[i*3 + 2].GetString());
74 pdg = param[i*3].GetInt();
75 fEfficiencyMap.insert(make_pair(pdg,make_pair(param[i*3 + 1].GetInt(),formula)));
76
77 // cout<<param[i*3].GetInt()<<","<<param[i*3+1].GetInt()<<","<<param[i*3 + 2].GetString()<<endl;
78
79 }
80
81 // set default efficiency formula
82 itEfficiencyMap = fEfficiencyMap.find(0);
83 if(itEfficiencyMap == fEfficiencyMap.end())
84 {
85 formula = new DelphesFormula;
86 formula->Compile("1.0");
87
88 fEfficiencyMap.insert(make_pair(0,make_pair(0,formula)));
89 }
90
91 // import input array
92
93 fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles"));
94 fItInputArray = fInputArray->MakeIterator();
95
96 // create output array
97
98 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
99}
100
101//------------------------------------------------------------------------------
102
103void IdentificationMap::Finish()
104{
105 if(fItInputArray) delete fItInputArray;
106
107 TMisIDMap::iterator itEfficiencyMap;
108 DelphesFormula *formula;
109 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
110 {
111 formula = (itEfficiencyMap->second).second;
112 if(formula) delete formula;
113 }
114
115}
116
117//------------------------------------------------------------------------------
118
119void IdentificationMap::Process()
120{
121 Candidate *candidate;
122 Double_t pt, eta, phi;
123 TMisIDMap::iterator itEfficiencyMap;
124 pair <TMisIDMap::iterator, TMisIDMap::iterator> range;
125 DelphesFormula *formula;
126 Int_t pdgIn, pdgOut, charge;
127
128 Double_t P, Pi;
129
130 // cout<<"------------ New Event ------------"<<endl;
131
132 fItInputArray->Reset();
133 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
134 {
135 const TLorentzVector &candidatePosition = candidate->Position;
136 const TLorentzVector &candidateMomentum = candidate->Momentum;
137 eta = candidatePosition.Eta();
138 phi = candidatePosition.Phi();
139 pt = candidateMomentum.Pt();
140 pdgIn = candidate->PID;
141 charge = candidate->Charge;
142
143 // cout<<"------------ New Candidate ------------"<<endl;
144 // cout<<candidate->PID<<" "<<pt<<","<<eta<<","<<phi<<endl;
145
146 P = 1.0;
147
148 //first check that PID of this particle is specified in cfg, if not set look for PID=0
149
150 itEfficiencyMap = fEfficiencyMap.find(pdgIn);
151
152 range = fEfficiencyMap.equal_range(pdgIn);
153 if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgIn);
154 if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
155
156 //loop over submap for this pid
157 for (TMisIDMap::iterator it=range.first; it!=range.second; ++it)
158 {
159
160 formula = (it->second).second;
161 pdgOut = (it->second).first;
162
163 Pi = formula->Eval(pt, eta);
164
165 // cout<<"PdgIn: "<<pdgIn<<", PdgOut: "<<pdgOut<<", Pi = "<<Pi<<endl;
166
167 // check that sum of probabilities does not exceed 1.
168
169 P = (P - Pi)/P;
170
171 if( P < 0.0 )
172 {
173 // cout<<"Sum of probabilities exceeds 1 for particle ID "<<pdgIn<<", check formulas in IdentificationMap module"<<endl;
174 continue;
175 }
176 else
177 {
178
179 //randomly assign a PID to particle according to map
180
181 Double_t rndm = gRandom->Uniform();
182 // cout<<"Random number vs P // "<<rndm<<" "<<P<<endl;
183
184 if(rndm > P)
185 //if(gRandom->Uniform() > P)
186 {
187 // cout<<"particle "<<candidate->PID<<" --> ";
188
189 //change PID of particle
190 if(pdgOut != 0) candidate->PID = charge*pdgOut;
191 // cout<<"Final PID: "<<candidate->PID<<endl;
192
193 fOutputArray->Add(candidate);
194
195 }
196 }
197
198 }
199
200 }
201
202}
203
204//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.