Belle II Software light-2406-ragdoll
BremsFinderModule.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/BremsFinderModule.h>
11
12// framework aux
13#include <framework/gearbox/Const.h>
14#include <framework/logging/Logger.h>
15#include <framework/datastore/RelationArray.h>
16#include <framework/datastore/RelationVector.h>
17
18// dataobjects
19#include <mdst/dataobjects/Track.h>
20#include <mdst/dataobjects/ECLCluster.h>
21
22// utilities
23#include <analysis/DecayDescriptor/ParticleListName.h>
24#include <analysis/DecayDescriptor/DecayDescriptor.h>
25
26// variables
27#include <analysis/variables/ECLVariables.h>
28
29#include <cmath>
30#include <algorithm>
31#include <TMatrixFSym.h>
32
33using namespace std;
34using namespace Belle2;
35
36//-----------------------------------------------------------------
37// Register module
38//-----------------------------------------------------------------
39
40REG_MODULE(BremsFinder);
41
42//-----------------------------------------------------------------
43// Implementation
44//-----------------------------------------------------------------
45
47 Module(), m_pdgCode(0)
48{
49 // set module description (e.g. insert text)
50 setDescription(R"DOC(
51This module copies each particle in the ``inputList`` to the ``outputList`` and uses
52the results of the **eclTrackBremFinder** module to look for possible bremsstrahlung photons; if these
53photons exist, it adds their four momentum to the particle in the ``outputList``.
54It also adds the original particle and these photons as daughters of the new, corrected particle.
55Track and PID information of the original particle are copied onto the new one to facilitate their access
56in the analysis scripts.
57
58The **eclTrackBremFinder** module uses the lepton track PXD and SVD hits and extrapolates them to the ECL;
59then looks for ECL clusters with energies between 0.02 and 1 times the track energy and without associated
60tracks, and checks if the normalized distance between each of these clusters
61and the extrapolated hit is smaller than 0.05. If it is, a *Bremsstrahlung* weighted relation
62between said cluster and the track is established. The weight is determined as
63
64.. math::
65
66 \text{max}\left(\frac{\left|\phi_{\text{cluster}}-\phi_{\text{hit}}\right|}{\Delta\phi_{\text{cluster}}+\Delta\phi_{\text{hit}}}, \, \frac{\left|\theta_{\text{cluster}}-\theta_{\text{hit}}\right|}{\Delta\theta_{\text{cluster}}+\Delta\theta_{\text{hit}}}\right)
67
68where :math:`\phi_i` and :math:`\theta_i` are the azimuthal and polar angles of the ECL cluster and the
69extrapolated hit, and :math:`\Delta x` represents the uncertainty of the value :math:`x`. The details of the
70calculation of these quantities are `here`_. By default, only relations with a weight smaller than 3.0 are stored.
71The user can further reduce the maximally allowed value of this weight to remove unwanted photons from the
72bremsstrahlung correction.
73
74This module looks for photons in the ``gammaList`` whose clusters have a *Bremsstrahlung* relation with the track
75of one of the particles in the ``inputList``, and adds their 4-momentum to the particle's one. It also stores the value
76of each relation weight as ``extraInfo`` of the corrected particle, under the name ``"bremsWeightWithPhotonN"``, where
77N is the index of the photon as daughter of the corrected particle; thus ``"bremsWeightWithPhoton0"`` gives the weight
78of the Bremsstrahlung relation between the new, corrected particle, and the first photon daughter.
79
80Warning:
81 Even in the event of no bremsstrahlung photons found, a new particle is still created, and the original one is still
82 added as its daughter.
83
84Warning:
85 Studies have shown that the requirements that the energy of the photon must be between 0.2 and 1 times the track energy and
86 that only track-photon relations with a weight below three are considered, are too tight. Until these are relaxed and a new
87 processing is done (MC15 and proc 13) it might be better to use the alternative `BelleBremRecovery` module.
88
89See also:
90 `eclTrackBremFinder module`_
91
92.. _eclTrackBremFinder module: https://gitlab.desy.de/belle2/software/basf2/-/tree/main/ecl/modules/eclTrackBremFinder
93.. _here: https://gitlab.desy.de/belle2/software/basf2/-/tree/main/ecl/modules/eclTrackBremFinder/src/BremFindingMatchCompute.cc
94)DOC");
96
97 // Add parameters
98 addParam("outputList", m_outputListName, "The output particle list name.");
99 addParam("inputList", m_inputListName,
100 R"DOC(The initial particle list name containing the particles to correct. *It should already exist* and *the particles must have an associated track.*)DOC");
101 addParam("gammaList", m_gammaListName,
102 R"DOC(The photon list containing the preselected bremsstrahlung candidates. *It should already exist* and *the particles in the list must be photons*)DOC");
103 addParam("maximumAcceptance", m_maximumAcceptance,
104 "The maximum value of the relation weight between a bremsstrahlung cluster and a particle track",
106 addParam("multiplePhotons", m_addMultiplePhotons, "If true, use all possible photons to correct the particle's 4-momentum",
108 addParam("usePhotonOnlyOnce", m_usePhotonOnlyOnce,
109 R"DOC(If true, each brems candidate is used to correct maximum 1 particle (the one with the lowest relation weight among all in the ``inputList``).)DOC",
111 addParam("writeOut", m_writeOut,
112 R"DOC(If true, the output ``ParticleList`` will be saved by `RootOutput`. If false, it will be ignored when writing the file.)DOC",
113 m_writeOut);
114}
115
117{
118 //Get input lepton particle list
119 m_inputList.isRequired(m_inputListName);
120
121 DecayDescriptor decDes;
122 decDes.init(m_inputListName);
123
124 m_pdgCode = decDes.getMother()->getPDGCode();
125 //Check for the particles in the lepton list to have a track associated
127 B2ERROR("[BremsFinderModule] Invalid particle list. the particles in " << m_inputListName << " should have an associated track.");
128
129 //Get input photon particle list
130 m_gammaList.isRequired(m_gammaListName);
131
132 decDes.init(m_gammaListName);
133 int temp_pdg = decDes.getMother()->getPDGCode();
134
135 //Check that this is a gamma list
136 if (temp_pdg != Const::photon.getPDGCode())
137 B2ERROR("[BremsFinderModule] Invalid particle list. the particles in " << m_gammaListName << " should be photons!");
138
139 decDes.init(m_outputListName);
140 temp_pdg = decDes.getMother()->getPDGCode();
141 // check the validity of output ParticleList name
143 B2ERROR("[BremsFinderModule] Input and output particle list names are the same: " << m_inputListName);
144 else if (temp_pdg != m_pdgCode) {
145 B2ERROR("[BremsFinderModule] The input and output particle list correspond to different particles: " << m_inputListName << " --> "
147 }
148
149 // output particle
151
152 // make output lists
154 m_outputList.registerInDataStore(m_outputListName, flags);
155 m_outputAntiList.registerInDataStore(m_outputAntiListName, flags);
156
158}
159
161{
162 RelationArray particlesToMCParticles(m_particles, m_mcParticles);
163
164 // new output particle list
165 m_outputList.create();
167
168 m_outputAntiList.create();
170 m_outputAntiList->bindAntiParticleList(*(m_outputList));
171
172 // Number of photons (calculate it here only once)
173 const unsigned int nGamma = m_gammaList->getListSize();
174
175 // Number of leptons (calculate it here only once)
176 const unsigned int nLep = m_inputList->getListSize();
177
178 const std::string relationName = "Bremsstrahlung";
179
180 //In the case of only one track per photon
182 for (unsigned n = 0; n < nGamma; n++) {
183 Particle* gamma = m_gammaList->getParticle(n);
184 //Skip this if it the best match has already been assigned (pathological case: happens only if you use the same gamma list
185 //to correct more than once. Performance studies, in which the same list is used with different options, are an example
186 if (gamma->hasExtraInfo("bestMatchIndex")) continue;
187
188 auto cluster = gamma->getECLCluster();
189 //Get the tracks related to each photon...
190 RelationVector<Track> relatedTracks = cluster->getRelationsTo<Track>("", relationName);
191 double bestWeight = m_maximumAcceptance;
192 unsigned bestMatchIndex = 0;
193 unsigned trkIndex = 0;
194 //Loop over the related tracks...
195 for (auto trk = relatedTracks.begin(); trk != relatedTracks.end(); trk++, trkIndex++) {
196 //... and over the input particles' tracks...
197 for (unsigned i = 0; i < nLep; i++) {
198 const Particle* lepton = m_inputList->getParticle(i);
199 auto leptonTrack = lepton->getTrack();
200 //... check that the particle track corresponds to the related track....
201 if (leptonTrack->getArrayIndex() == trk->getArrayIndex()) {
202 double weight = relatedTracks.weight(trkIndex);
203 if (weight < bestWeight) {
204 bestWeight = weight;
205 //... and only select the best match among the tracks in the input list
206 bestMatchIndex = trk->getArrayIndex();
207 }
208 break; //If the particle corresponding to the related track is found, break the loop over the particles and go for the next related track
209 }
210 }
211 }
212 //... finally, add the best match index as an extra info for the photon
213 gamma->addExtraInfo("bestMatchIndex", bestMatchIndex);
214 }
215 }
216
217 // loop over charged particles, correct them and add them to the output list
218
219 for (unsigned i = 0; i < nLep; i++) {
220 const Particle* lepton = m_inputList->getParticle(i);
221
222 //Get the track of this lepton...
223 auto track = lepton->getTrack();
224
225 //... and get the bremsstrahlung clusters related to this track
226 RelationVector<ECLCluster> bremClusters = track->getRelationsFrom<ECLCluster>("", relationName);
227
228 std::vector<std::pair <double, Particle*> > selectedGammas;
229
230 unsigned j = 0;
231 for (auto bremCluster = bremClusters.begin(); bremCluster != bremClusters.end(); bremCluster++, j++) {
232 double weight = bremClusters.weight(j);
233
234 if (weight > m_maximumAcceptance) continue;
235
236 for (unsigned k = 0; k < nGamma; k++) {
237
238 Particle* gamma = m_gammaList->getParticle(k);
239 auto cluster = gamma->getECLCluster();
240
241 if (bremCluster->getUniqueClusterId() == cluster->getUniqueClusterId()) {
242 if (m_usePhotonOnlyOnce) { //If only one track per photon should be used...
243 if (track->getArrayIndex() == static_cast<int>
244 (gamma->getExtraInfo("bestMatchIndex"))) //... check if this track is the best match ...
245 selectedGammas.push_back(std::make_pair(weight, gamma)); //... and if it is, add it to the selected gammas
246 } else {
247 selectedGammas.push_back(std::make_pair(weight, gamma));
248 }
249 }
250
251 } // Closes for loop on gammas
252
253 } // Closes for loop on brem clusters
254
255 //The 4-momentum of the new lepton in the output particle list
256 ROOT::Math::PxPyPzEVector new4Vec = lepton->get4Vector();
257
258 //Sort weight-particle pairs by weight. Smaller weights go first
259 std::sort(selectedGammas.begin(), selectedGammas.end());
260
261 //Add to this 4-momentum those of the selected photon(s)
262 for (auto const& bremsPair : selectedGammas) {
263 Particle* g = bremsPair.second;
264 new4Vec += g->get4Vector();
265 if (! m_addMultiplePhotons) break; //stop after adding the first photon
266 }
267
268 //Create the new particle with the 4-momentum calculated before
269 Particle correctedLepton(new4Vec, lepton->getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
270 track->getArrayIndex());
271
272 //And add the original lepton as its daughter
273 correctedLepton.appendDaughter(lepton, false);
274
275 const TMatrixFSym& lepErrorMatrix = lepton->getMomentumVertexErrorMatrix();
276 TMatrixFSym corLepMatrix(lepErrorMatrix);
277
278 double bremsGammaEnergySum = 0.0;
279 //Now, if there are any, add the brems photons as daughters as well. As before, we distinguish between the multiple and only one brems photon cases
280 int photonIndex = 0;
281 for (auto const& bremsPair : selectedGammas) {
282 //Add the weights as extra info of the mother
283 Particle* bremsGamma = bremsPair.second;
284 std::string extraInfoName = "bremsWeightWithPhoton" + std::to_string(photonIndex);
285 correctedLepton.addExtraInfo(extraInfoName, bremsPair.first);
286 photonIndex++;
287 bremsGammaEnergySum += Variable::eclClusterE(bremsGamma);
288
289 const TMatrixFSym& gammaErrorMatrix = bremsGamma->getMomentumVertexErrorMatrix();
290 for (int irow = 0; irow <= 3; irow++) {
291 for (int icol = irow; icol <= 3; icol++) corLepMatrix(irow, icol) += gammaErrorMatrix(irow, icol);
292 }
293 correctedLepton.appendDaughter(bremsGamma, false);
294 B2DEBUG(10, "[BremsFinderModule] Found a bremsstrahlung gamma and added its 4-vector to the charged particle");
295 if (! m_addMultiplePhotons) break; //stop after adding the first photon
296 }
297
298 correctedLepton.setMomentumVertexErrorMatrix(corLepMatrix);
299
300 // add the info from original lepton to the new lepton
301 correctedLepton.setVertex(lepton->getVertex());
302 correctedLepton.setPValue(lepton->getPValue());
303 correctedLepton.addExtraInfo("bremsCorrected", float(selectedGammas.size() > 0));
304 correctedLepton.addExtraInfo("bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
305
306 // add the mc relation
307 Particle* newLepton = m_particles.appendNew(correctedLepton);
308 const MCParticle* mcLepton = lepton->getRelated<MCParticle>();
309 const PIDLikelihood* pid = lepton->getPIDLikelihood();
310
311 if (pid) newLepton->addRelationTo(pid);
312
313 if (mcLepton != nullptr) newLepton->addRelationTo(mcLepton);
314
315 m_outputList->addParticle(newLepton);
316
317 } //Closes for loop on leptons
318
319} //Close event()
virtual void initialize() override
Use this to initialize resources or memory your module needs.
double m_maximumAcceptance
photons whose clusters have relation weights higher than this will not be used for bremsstrahlung cor...
std::string m_gammaListName
input gamma list name
virtual void event() override
Called once for each event.
bool m_addMultiplePhotons
In case there is more than one brems photon, use only the best one (based on the weight of the relati...
StoreObjPtr< ParticleList > m_outputAntiList
StoreObjptr for output antiparticlelist.
StoreArray< Particle > m_particles
StoreArray of Particle objects.
StoreObjPtr< ParticleList > m_outputList
StoreObjptr for output particlelist.
StoreArray< PIDLikelihood > m_pidlikelihoods
StoreArray of PIDLikelihood objects.
StoreObjPtr< ParticleList > m_inputList
StoreObjptr for input charged particle list.
StoreObjPtr< ParticleList > m_gammaList
StoreObjptr for gamma list.
bool m_writeOut
Write the output particle list in the final file?
StoreArray< MCParticle > m_mcParticles
StoreArray of MCParticle objects.
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 particle to be corrected.
std::string m_outputAntiListName
output anti-particle list name
std::string m_inputListName
input particle list name
BremsFinderModule()
Constructor, for setting module description and parameters.
std::string m_outputListName
output particle list name
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
static const ParticleType invalidParticle
Invalid particle, used internally.
Definition: Const.h:681
static const ParticleType photon
photon particle
Definition: Const.h:673
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
int getPDGCode() const
Return PDG code.
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
ECL cluster data.
Definition: ECLCluster.h:27
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
const ECLCluster * getECLCluster() const
Returns the pointer to the ECLCluster object that was used to create this Particle (if ParticleType =...
Definition: Particle.cc:891
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 setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:393
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:366
@ 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
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
Class for type safe access to objects that are referred to in relations.
iterator end()
Return iterator to last entry +1.
iterator begin()
Return iterator to first entry.
float weight(int index) const
Get weight with index.
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).
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
Class that bundles various TrackFitResults.
Definition: Track.h:25
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.