Fork me on GitHub

source: git/modules/IdentificationMap.cc@ 50d7259

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 50d7259 was a5bbe8a, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix probability calculation in IdentificationMap

  • Property mode set to 100644
File size: 5.0 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 IdentificationMap
21 *
22 * Converts particles with some PDG code into another particle,
23 * according to parametrized probability.
24 *
25 * \author M. Selvaggi - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/IdentificationMap.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
56IdentificationMap::IdentificationMap() :
57 fItInputArray(0)
58{
59}
60
61//------------------------------------------------------------------------------
62
63IdentificationMap::~IdentificationMap()
64{
65}
66
67//------------------------------------------------------------------------------
68
69void IdentificationMap::Init()
70{
71 TMisIDMap::iterator itEfficiencyMap;
72 ExRootConfParam param;
73 DelphesFormula *formula;
74 Int_t i, size, pdg;
75
76 // read efficiency formulas
77 param = GetParam("EfficiencyFormula");
78 size = param.GetSize();
79
80 fEfficiencyMap.clear();
81 for(i = 0; i < size/3; ++i)
82 {
83 formula = new DelphesFormula;
84 formula->Compile(param[i*3 + 2].GetString());
85 pdg = param[i*3].GetInt();
86 fEfficiencyMap.insert(make_pair(pdg, make_pair(param[i*3 + 1].GetInt(), formula)));
87 }
88
89 // set default efficiency formula
90 itEfficiencyMap = fEfficiencyMap.find(0);
91 if(itEfficiencyMap == fEfficiencyMap.end())
92 {
93 formula = new DelphesFormula;
94 formula->Compile("1.0");
95
96 fEfficiencyMap.insert(make_pair(0, make_pair(0, formula)));
97 }
98
99 // import input array
100
101 fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles"));
102 fItInputArray = fInputArray->MakeIterator();
103
104 // create output array
105
106 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
107}
108
109//------------------------------------------------------------------------------
110
111void IdentificationMap::Finish()
112{
113 if(fItInputArray) delete fItInputArray;
114
115 TMisIDMap::iterator itEfficiencyMap;
116 DelphesFormula *formula;
117 for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
118 {
119 formula = (itEfficiencyMap->second).second;
120 if(formula) delete formula;
121 }
122}
123
124//------------------------------------------------------------------------------
125
126void IdentificationMap::Process()
127{
128 Candidate *candidate;
129 Double_t pt, eta, phi;
130 TMisIDMap::iterator itEfficiencyMap;
131 pair <TMisIDMap::iterator, TMisIDMap::iterator> range;
132 DelphesFormula *formula;
133 Int_t pdgCodeIn, pdgCodeOut, charge;
134
135 Double_t p, r, total;
136
137 fItInputArray->Reset();
138 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
139 {
140 const TLorentzVector &candidatePosition = candidate->Position;
141 const TLorentzVector &candidateMomentum = candidate->Momentum;
142 eta = candidatePosition.Eta();
143 phi = candidatePosition.Phi();
144 pt = candidateMomentum.Pt();
145 pdgCodeIn = candidate->PID;
146 charge = candidate->Charge;
147
148 // first check that PID of this particle is specified in the map
149 // otherwise, look for PID = 0
150
151 itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn);
152
153 range = fEfficiencyMap.equal_range(pdgCodeIn);
154 if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn);
155 if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
156
157 r = gRandom->Uniform();
158 total = 0.0;
159
160 // loop over sub-map for this PID
161 for(TMisIDMap::iterator it = range.first; it != range.second; ++it)
162 {
163 formula = (it->second).second;
164 pdgCodeOut = (it->second).first;
165
166 p = formula->Eval(pt, eta);
167
168 if(total <= r && r < total + p)
169 {
170 // change PID of particle
171 if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut;
172 fOutputArray->Add(candidate);
173 break;
174 }
175
176 total += p;
177 }
178 }
179}
180
181//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.