Belle II Software light-2406-ragdoll
BelleBremRecoveryModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// Own header.
10#include <analysis/modules/BremsCorrection/BelleBremRecoveryModule.h>
11// framework aux
12#include <framework/logging/Logger.h>
13#include <framework/datastore/RelationArray.h>
14
15// dataobjects
16#include <mdst/dataobjects/Track.h>
17
18// utilities
19#include <analysis/DecayDescriptor/ParticleListName.h>
20
21// variables
22#include <analysis/variables/ECLVariables.h>
23
24#include <algorithm>
25#include <TMatrixFSym.h>
26#include <Math/Vector3D.h>
27#include <Math/Vector4D.h>
28#include <Math/VectorUtil.h>
29
30#include <vector>
31
32using namespace std;
33using namespace Belle2;
34
35//-----------------------------------------------------------------
36// Register module
37//-----------------------------------------------------------------
38
39REG_MODULE(BelleBremRecovery);
40
41//-----------------------------------------------------------------
42// Implementation
43//-----------------------------------------------------------------
44
46 Module(), m_pdgCode(0)
47
48{
49 // set module description (e.g. insert text)
50 setDescription(R"DOC(
51Takes the charged particle from the given charged particle list (``inputListName``) and
52copies it to the output list (``outputListName``). The 4-vector of the nearest (all) photon(s)
53from ``gammaListName`` (considered as radiative) is added to the charged particle, if it is
54found inside the cone around the charged particle with the given maximum angle (``angleThreshold``).
55 )DOC");
57
58 // Add parameters
59 addParam("inputListName", m_inputListName,
60 "The initial charged particle list containing the charged particles to correct, should already exist.");
61 addParam("outputListName", m_outputListName, "The output charged particle list containing the corrected charged particles.");
62 addParam("gammaListName", m_gammaListName, "The gammas list containing possibly radiative gammas, should already exist.");
63 addParam("angleThreshold", m_angleThres,
64 "The maximum angle in radians between the charged particle and the (radiative) gamma to be accepted.", 0.05);
65 addParam("multiplePhotons", m_isMultiPho, "If only the nearest photon to add then make it False otherwise true", true);
66 addParam("usePhotonOnlyOnce", m_usePhotonOnlyOnce,
67 "If false, a photon is used for correction of the closest charged particle in the inputList. "
68 "If true, a photon is allowed to be used for correction multiple times. "
69 "WARNING: One cannot use a photon twice to reconstruct a composite particle. Thus, for example, if e+ and e- are corrected with a gamma, they cannot form a J/psi -> e+ e- candidate.",
70 false);
71
72 addParam("writeOut", m_writeOut,
73 "Set to true if you want to write out the output list to a root file", false);
74
75}
76
78{
79 // check the validity of output ParticleList name
81 if (!valid)
82 B2ERROR("[BelleBremRecoveryModule] Invalid output ParticleList name: " << m_outputListName);
83
84 // output particle
86 m_pdgCode = mother->getPDGCode();
88
89 // get existing particle lists
91 B2ERROR("[BelleBremRecoveryModule] Input and output particle list names are the same: " << m_inputListName);
93 B2ERROR("[BelleBremRecoveryModule] Invalid input particle list name: " << m_inputListName);
94 } else {
96 }
97
99 B2ERROR("[BelleBremRecoveryModule] Invalid gamma particle list name: " << m_gammaListName);
100 } else {
101 m_gammaList.isRequired(m_gammaListName);
102 }
103
104 // make output list
106 m_outputparticleList.registerInDataStore(m_outputListName, flags);
107 m_outputAntiparticleList.registerInDataStore(m_outputAntiListName, flags);
108
110}
111
112
114{
115 // new output particle list
116 m_outputparticleList.create();
118
121 m_outputAntiparticleList->bindAntiParticleList(*(m_outputparticleList));
122
123 const unsigned int nLep = m_inputparticleList->getListSize();
124 const unsigned int nGam = m_gammaList->getListSize();
125
126 if (nLep == 0) return;
127
128 const std::string extraInfoSuffix = "BelleBremRecovery_" + m_inputListName + "_" + m_gammaListName;
129
130 std::vector<std::vector<Particle*>> selectedGammas_nLep; // only used if m_usePhotonOnlyOnce is true
131 if (m_usePhotonOnlyOnce) selectedGammas_nLep.resize(nLep);
132
133 // loop over charged particles
134 for (unsigned iLep = 0; iLep < nLep; iLep++) {
135 std::vector<Particle*> selectedGammas;
136
137 const Particle* lepton = m_inputparticleList->getParticle(iLep);
138 const ROOT::Math::XYZVector pLep = lepton->getMomentum();
139
140 // look for all possible (radiative) gamma
141 for (unsigned iGam = 0; iGam < nGam; iGam++) {
142 Particle* gamma = m_gammaList->getParticle(iGam); // to write extraInfo
143 const ROOT::Math::XYZVector pGam = gamma->getMomentum();
144
145 const double angle = ROOT::Math::VectorUtil::Angle(pLep, pGam);
146
147 if (angle < m_angleThres) {
148 if (!m_usePhotonOnlyOnce) { // if multiple use is allowed, keep all selected gammas
149 gamma->writeExtraInfo("angle_" + extraInfoSuffix, angle);
150 gamma->writeExtraInfo("indexLep_" + extraInfoSuffix, iLep);
151 selectedGammas.push_back(gamma);
152 } else {
153 if (gamma->hasExtraInfo("angle_" + extraInfoSuffix)) { // check if the gamma has already selected
154 if (angle < gamma->getExtraInfo("angle_" + extraInfoSuffix)) {
155 // this lepton is closer to the gamma than previous lepton
156 int indexLep_prev = gamma->getExtraInfo("indexLep_" + extraInfoSuffix);
157 selectedGammas_nLep[indexLep_prev].erase(std::remove(selectedGammas_nLep[indexLep_prev].begin(),
158 selectedGammas_nLep[indexLep_prev].end(),
159 gamma),
160 selectedGammas_nLep[indexLep_prev].end());
161 gamma->writeExtraInfo("angle_" + extraInfoSuffix, angle);
162 gamma->writeExtraInfo("indexLep_" + extraInfoSuffix, iLep);
163 selectedGammas.push_back(gamma);
164 } // else, the gamma will not be used for correction
165 } else {
166 gamma->writeExtraInfo("angle_" + extraInfoSuffix, angle);
167 gamma->writeExtraInfo("indexLep_" + extraInfoSuffix, iLep);
168 selectedGammas.push_back(gamma);
169 }
170 }
171 }
172 }
173
174 // sorting the bremphotons in ascending order of the angle with the charged particle
175 std::sort(selectedGammas.begin(), selectedGammas.end(),
176 [extraInfoSuffix](const Particle * photon1, const Particle * photon2) {
177 return photon1->getExtraInfo("angle_" + extraInfoSuffix) < photon2->getExtraInfo("angle_" + extraInfoSuffix);
178 });
179
180 if (m_usePhotonOnlyOnce) // store the selectedGammas to check the closest lepton to gammas
181 selectedGammas_nLep[iLep] = selectedGammas;
182 else // perform correction here.
183 correctLepton(lepton, selectedGammas);
184
185 }
186
188 for (unsigned iLep = 0; iLep < nLep; iLep++) {
189 const Particle* lepton = m_inputparticleList->getParticle(iLep);
190 correctLepton(lepton, selectedGammas_nLep[iLep]);
191 }
192 }
193
194}
195
196void BelleBremRecoveryModule::correctLepton(const Particle* lepton, std::vector<Particle*> selectedGammas)
197{
198 // Updates 4-vector and errorMatrix
199 ROOT::Math::PxPyPzEVector new4Vec = lepton->get4Vector();
200 const TMatrixFSym& lepErrorMatrix = lepton->getMomentumVertexErrorMatrix();
201 TMatrixFSym corLepMatrix = lepErrorMatrix;
202
203 // Sum of energy of correction gammas for the extraInfo
204 double bremsGammaEnergySum = 0.0;
205
206 // Create a correctedLepton. 4-vector will be updated
207 Particle correctedLepton(new4Vec, lepton->getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
208 lepton->getTrack()->getArrayIndex());
209
210 correctedLepton.setVertex(lepton->getVertex());
211 correctedLepton.setPValue(lepton->getPValue());
212
213 correctedLepton.appendDaughter(lepton, false);
214
215 for (auto gamma : selectedGammas) {
216
217 new4Vec += gamma->get4Vector();
218 bremsGammaEnergySum += Variable::eclClusterE(gamma);
219
220 const TMatrixFSym& fsrErrorMatrix = gamma->getMomentumVertexErrorMatrix();
221 for (int irow = 0; irow <= 3 ; irow++)
222 for (int icol = irow; icol <= 3; icol++)
223 corLepMatrix(irow, icol) += fsrErrorMatrix(irow, icol);
224
225 correctedLepton.appendDaughter(gamma, false);
226
227 B2DEBUG(19, "[BelleBremRecoveryModule] Found a radiative gamma and added its 4-vector to the charge particle");
228
229 if (!m_isMultiPho) break;
230 }
231
232 correctedLepton.set4Vector(new4Vec);
233 correctedLepton.updateMass(m_pdgCode);
234 correctedLepton.setMomentumVertexErrorMatrix(corLepMatrix);
235 correctedLepton.addExtraInfo("bremsCorrected", float(bremsGammaEnergySum > 0));
236 correctedLepton.addExtraInfo("bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
237
238 // add the mc relation
239 Particle* newLepton = m_particles.appendNew(correctedLepton);
240 const MCParticle* mcLepton = lepton->getRelated<MCParticle>();
241 const PIDLikelihood* pid = lepton->getPIDLikelihood();
242 if (pid)
243 newLepton->addRelationTo(pid);
244 if (mcLepton)
245 newLepton->addRelationTo(mcLepton);
246
247 m_outputparticleList->addParticle(newLepton);
248
249}
void correctLepton(const Particle *lepton, std::vector< Particle * > selectedGammas)
Correct lepton kinematics using the selectedGammas.
StoreObjPtr< ParticleList > m_outputparticleList
StoreObjptr for output particlelist.
virtual void initialize() override
Initialize the Module.
std::string m_gammaListName
input ParticleList names
virtual void event() override
Event processor.
bool m_isMultiPho
multiple or one bremphoton addition option
StoreArray< Particle > m_particles
StoreArray of Particle objects.
double m_angleThres
input max angle to be accepted (in radian)
StoreArray< PIDLikelihood > m_pidlikelihoods
StoreArray of PIDLikelihood objects.
StoreObjPtr< ParticleList > m_gammaList
StoreObjptr for gamma list.
StoreObjPtr< ParticleList > m_inputparticleList
StoreObjptr for input charged particlelist.
StoreObjPtr< ParticleList > m_outputAntiparticleList
StoreObjptr for output antiparticlelist.
DecayDescriptor m_decaydescriptorGamma
Decay descriptor of the decay being reconstructed.
DecayDescriptor m_decaydescriptor
Decay descriptor of the charged particle decay.
bool m_writeOut
toggle output particle list btw.
bool m_usePhotonOnlyOnce
Each brems photon can be used to correct only one particle (the one with the smallest relation weight...
int m_pdgCode
PDG code of the combined mother particle.
std::string m_outputAntiListName
output anti-particle list name
std::string m_inputListName
input ParticleList names
std::string m_outputListName
output ParticleList name
EStoreFlags
Flags describing behaviours of objects etc.
Definition: DataStore.h:69
@ c_WriteOut
Object/array should be saved by output modules.
Definition: DataStore.h:70
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:71
Represents a particle in the DecayDescriptor.
int getPDGCode() const
Return PDG code.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:80
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:29
Class to store reconstructed particles.
Definition: Particle.h:75
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:845
void appendDaughter(const Particle *daughter, const bool updateType=true, const int daughterProperty=c_Ordinary)
Appends index of daughter to daughters index array.
Definition: Particle.cc:676
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1308
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:295
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:667
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:631
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1266
const PIDLikelihood * getPIDLikelihood() const
Returns the pointer to the PIDLikelihood object that is related to the Track, which was used to creat...
Definition: Particle.cc:871
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1336
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector.
Definition: Particle.h:271
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:393
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:366
void updateMass(const int pdgCode)
Updates particle mass with the mass of the particle corresponding to the given PDG.
Definition: Particle.cc:597
@ c_Flavored
Is either particle or antiparticle.
Definition: Particle.h:96
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:420
double getExtraInfo(const std::string &name) const
Return given value if set.
Definition: Particle.cc:1289
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
std::string antiParticleListName(const std::string &listName)
Returns name of anti-particle-list corresponding to listName.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.