Fork me on GitHub

source: git/modules/IdentificationMap.cc@ 0983d48

Last change on this file since 0983d48 was a3261d7, checked in by GitHub <noreply@…>, 4 years ago

Clone candidate before altering its PID property

Currently, if a particle is accepted for mapping onto a new particle ID (PID), the original copy of the candidate (which might be owned by a module upstream of the IdentificationMap instance) is altered. This may have unintended consequences. For example, what if the intent is not to alter the original particle, but rather to create a new list of particles with alternate PID settings? If one alters the original candidate's PID setting, then after-the-fact it's not possible to compare how the PID algorithm performed (e.g. you lose the ability to do efficiency and mis-identification studies to cross-check the implementation of the IdentificaionMap in the first place).

This change clones the candidate before then altering its PID variable and storing it in the output list.

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