Belle II Software  release-05-02-19
BremsFinderModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Alejandro Mora *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <analysis/modules/BremsCorrection/BremsFinderModule.h>
13 
14 // framework aux
15 #include <framework/gearbox/Const.h>
16 #include <framework/logging/Logger.h>
17 #include <framework/datastore/RelationArray.h>
18 #include <framework/datastore/RelationVector.h>
19 #include <framework/datastore/StoreArray.h>
20 
21 // dataobjects
22 #include <analysis/dataobjects/Particle.h>
23 #include <mdst/dataobjects/MCParticle.h>
24 #include <mdst/dataobjects/PIDLikelihood.h>
25 #include <mdst/dataobjects/Track.h>
26 #include <mdst/dataobjects/ECLCluster.h>
27 
28 // utilities
29 #include <analysis/DecayDescriptor/ParticleListName.h>
30 #include <analysis/DecayDescriptor/DecayDescriptor.h>
31 
32 // variables
33 #include <analysis/variables/ECLVariables.h>
34 
35 #include <cmath>
36 #include <algorithm>
37 #include <TMatrixFSym.h>
38 
39 using namespace std;
40 
41 namespace Belle2 {
46 //-----------------------------------------------------------------
47 // Register module
48 //-----------------------------------------------------------------
49 
50  REG_MODULE(BremsFinder)
51 
52 //-----------------------------------------------------------------
53 // Implementation
54 //-----------------------------------------------------------------
55 
57  Module(), m_pdgCode(0)
58  {
59  // set module description (e.g. insert text)
60  setDescription(R"DOC(
61  This module copies each particle in the ``inputList`` to the ``outputList`` and uses
62  the results of the **eclTrackBremFinder** module to look for possible bremsstrahlung photons; if these
63  photons exists, it adds their four momentum to the particle in the ``outputList``.
64  It also adds the original particle and these photons as daughters of the new, corrected particle.
65  Track and PID information of the original particle are copied onto the new one to facilitate their access
66  in the analysis scripts.
67 
68  The **eclTrackBremFinder** module uses the lepton track PXD and SVD hits and extrapolates them to the ECL;
69  then looks for ECL clusters with energies between 0.2 and 1 times the track energy and without associated
70  tracks, and checks if the distance between each of these clusters
71  and the extrapolated hit is smaller than 0.5 mm. If it is, a *Bremsstrahlung* weighted relation
72  between said cluster and the track is established. The weight is determined as
73 
74  .. math::
75 
76  \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)
77 
78  where :math:`\phi_i` and :math:`\theta_i` are the azimutal and polar angles of the ECL cluster and the
79  extrapolated hit, and :math:`\Delta x` represents the uncertainty of the value :math:`x`. The details of
80  the calculation of these quantities are `here`_. By default, only relations with a weight
81  smaller than 3.0 are stored. The user can further determine the maximum value of this weight required in order
82  to perform the bremsstrahlung correction.
83 
84  This module looks for photons in the ``gammaList`` whose clusters have a *Bremsstrahlung* relation with the track
85  of one of the particles in the ``inputList``, and adds their 4-momentum to the particle's one. It also stores the value
86  of each relation weight as ``extraInfo`` of the corrected particle, under the name ``"bremsWeightWithPhotonN"``, where
87  N is the index of the photon as daughter of the corrected particle; thus ``"bremsWeightWithPhoton0"`` gives the weight
88  of the Bremsstrahlung relation between the new, corrected particle, and the first photon daughter.
89 
90  Warning:
91  Even in the event of no bremsstrahlung photons found, a new particle is still created, and the original one is still
92  added as its daughter.
93 
94  See also:
95  `eclTrackBremFinder module`_
96 
97  .. _eclTrackBremFinder module: https://stash.desy.de/projects/B2/repos/software/browse/ecl/modules/eclTrackBremFinder
98  .. _here: https://stash.desy.de/projects/B2/repos/software/browse/ecl/modules/eclTrackBremFinder/src/BremFindingMatchCompute.cc)DOC");
99  setPropertyFlags(c_ParallelProcessingCertified);
100 
101  // Add parameters
102  addParam("outputList", m_outputListName, "The output particle list name.");
103  addParam("inputList", m_inputListName,
104  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");
105  addParam("gammaList", m_gammaListName,
106  R"DOC(The photon list containing the preselected bremsstrahlung candidates. *It should already exist* and *the particles in the list must be photons*)DOC");
107  addParam("maximumAcceptance", m_maximumAcceptance,
108  "The maximum value of the relation weight between a bremsstrahlung cluster and a particle track",
109  m_maximumAcceptance);
110  addParam("multiplePhotons", m_addMultiplePhotons, "If true, use all possible photons to correct the particle's 4-momentum",
111  m_addMultiplePhotons);
112  addParam("usePhotonOnlyOnce", m_usePhotonOnlyOnce,
113  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",
114  m_usePhotonOnlyOnce);
115  addParam("writeOut", m_writeOut,
116  R"DOC(If true, the output ``ParticleList`` will be saved by `RootOutput`. If false, it will be ignored when writing the file.)DOC",
117  m_writeOut);
118  }
119 
120  void BremsFinderModule::initialize()
121  {
122  //Get input lepton particle list
123  m_inputList.isRequired(m_inputListName);
124 
125  DecayDescriptor decDes;
126  decDes.init(m_inputListName);
127 
128  m_pdgCode = decDes.getMother()->getPDGCode();
129  //Check for the particles in the lepton list to have a track associated
130  if (Const::chargedStableSet.find(abs(m_pdgCode)) == Const::invalidParticle)
131  B2ERROR("[BremsFinderModule] Invalid particle list. the particles in " << m_inputListName << " should have an associated track.");
132 
133  //Get input photon particle list
134  m_gammaList.isRequired(m_gammaListName);
135 
136  decDes.init(m_gammaListName);
137  int temp_pdg = decDes.getMother()->getPDGCode();
138 
139  //Check that this is a gamma list
140  if (temp_pdg != Const::photon.getPDGCode())
141  B2ERROR("[BremsFinderModule] Invalid particle list. the particles in " << m_gammaListName << " should be photons!");
142 
143  decDes.init(m_outputListName);
144  temp_pdg = decDes.getMother()->getPDGCode();
145  // check the validity of output ParticleList name
146  if (m_inputListName == m_outputListName)
147  B2ERROR("[BremsFinderModule] Input and output particle list names are the same: " << m_inputListName);
148  else if (temp_pdg != m_pdgCode) {
149  B2ERROR("[BremsFinderModule] The input and output particle list correspond to different particles: " << m_inputListName << " --> "
150  << m_outputListName);
151  }
152 
153  // output particle
154  m_outputAntiListName = ParticleListName::antiParticleListName(m_outputListName);
155 
156  // make output lists
157  DataStore::EStoreFlags flags = m_writeOut ? DataStore::c_WriteOut : DataStore::c_DontWriteOut;
158  m_outputList.registerInDataStore(m_outputListName, flags);
159  m_outputAntiList.registerInDataStore(m_outputAntiListName, flags);
160 
161  StoreArray<Particle> particles;
162  StoreArray<PIDLikelihood> pidlikelihoods;
163  particles.registerRelationTo(pidlikelihoods);
164  }
165 
166  void BremsFinderModule::event()
167  {
168  StoreArray<Particle> particles;
169  StoreArray<MCParticle> mcParticles;
170 
171  RelationArray particlesToMCParticles(particles, mcParticles);
172 
173  // new output particle list
174  m_outputList.create();
175  m_outputList->initialize(m_pdgCode, m_outputListName);
176 
177  m_outputAntiList.create();
178  m_outputAntiList->initialize(-1 * m_pdgCode, m_outputAntiListName);
179  m_outputAntiList->bindAntiParticleList(*(m_outputList));
180 
181  // Number of photons (calculate it here only once)
182  const unsigned int nGamma = m_gammaList->getListSize();
183 
184  // Number of leptons (calculate it here only once)
185  const unsigned int nLep = m_inputList->getListSize();
186 
187  const std::string relationName = "Bremsstrahlung";
188 
189  //In the case of only one track per photon
190  if (m_usePhotonOnlyOnce) {
191  for (unsigned n = 0; n < nGamma; n++) {
192  Particle* gamma = m_gammaList->getParticle(n);
193  //Skip this if it the best match has already been asigned (pathological case: happens only if you use the same gamma list
194  //to correct more than once. Performance studies, in which the same list is used with different options, are an example
195  if (gamma->hasExtraInfo("bestMatchIndex")) continue;
196 
197  auto cluster = gamma->getECLCluster();
198  //Get the tracks related to each photon...
199  RelationVector<Track> relatedTracks = cluster->getRelationsTo<Track>("", relationName);
200  double bestWeight = m_maximumAcceptance;
201  unsigned bestMatchIndex = 0;
202  unsigned trkIndex = 0;
203  //Loop over the related tracks...
204  for (auto trk = relatedTracks.begin(); trk != relatedTracks.end(); trk++, trkIndex++) {
205  //... and over the input particles' tracks...
206  for (unsigned i = 0; i < nLep; i++) {
207  const Particle* lepton = m_inputList->getParticle(i);
208  auto leptonTrack = lepton->getTrack();
209  //... check that the particle track corresponds to the related track....
210  if (leptonTrack->getArrayIndex() == trk->getArrayIndex()) {
211  double weight = relatedTracks.weight(trkIndex);
212  if (weight < bestWeight) {
213  bestWeight = weight;
214  //... and only select the best match among the tracks in the input list
215  bestMatchIndex = trk->getArrayIndex();
216  }
217  break; //If the particle corresponding to the related track is found, break the loop over the particles and go for the next related track
218  }
219  }
220  }
221  //... finally, add the best match index as an extra info for the photon
222  gamma->addExtraInfo("bestMatchIndex", bestMatchIndex);
223  }
224  }
225 
226  // loop over charged particles, correct them and add them to the output list
227 
228  for (unsigned i = 0; i < nLep; i++) {
229  const Particle* lepton = m_inputList->getParticle(i);
230 
231  //Get the track of this lepton...
232  auto track = lepton->getTrack();
233 
234  //... and get the bremsstrahlung clusters related to this track
235  RelationVector<ECLCluster> bremClusters = track->getRelationsFrom<ECLCluster>("", relationName);
236 
237  std::vector<std::pair <double, Particle*> > selectedGammas;
238 
239  unsigned j = 0;
240  for (auto bremCluster = bremClusters.begin(); bremCluster != bremClusters.end(); bremCluster++, j++) {
241  double weight = bremClusters.weight(j);
242 
243  if (weight > m_maximumAcceptance) continue;
244 
245  for (unsigned k = 0; k < nGamma; k++) {
246 
247  Particle* gamma = m_gammaList->getParticle(k);
248  auto cluster = gamma->getECLCluster();
249 
250  if (bremCluster->getClusterId() == cluster->getClusterId()) {
251  if (m_usePhotonOnlyOnce) { //If only one track per photon should be used...
252  if (track->getArrayIndex() == static_cast<int>
253  (gamma->getExtraInfo("bestMatchIndex"))) //... check if this track is the best match ...
254  selectedGammas.push_back(std::make_pair(weight, gamma)); //... and if it is, add it to the selected gammas
255  } else {
256  selectedGammas.push_back(std::make_pair(weight, gamma));
257  }
258  }
259 
260  } // Closes for loop on gammas
261 
262  } // Closes for loop on brem clusters
263 
264  //The 4-momentum of the new lepton in the output particle list
265  TLorentzVector new4Vec = lepton->get4Vector();
266 
267  //Sort weight-particle pairs by weight. Smaller weights go first
268  std::sort(selectedGammas.begin(), selectedGammas.end());
269 
270  //Add to this 4-momentum those of the selected photon(s)
271  for (auto const& bremsPair : selectedGammas) {
272  Particle* g = bremsPair.second;
273  new4Vec += g->get4Vector();
274  if (! m_addMultiplePhotons) break; //stop after adding the first photon
275  }
276 
277  //Create the new particle with the 4-momentum calculated before
278  Particle correctedLepton(new4Vec, lepton->getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
279  track->getArrayIndex());
280 
281  //And add the original lepton as its daughter
282  correctedLepton.appendDaughter(lepton, false);
283 
284  const TMatrixFSym& lepErrorMatrix = lepton->getMomentumVertexErrorMatrix();
285  TMatrixFSym corLepMatrix(lepErrorMatrix);
286 
287  double bremsGammaEnergySum = 0.0;
288  //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
289  int photonIndex = 0;
290  for (auto const& bremsPair : selectedGammas) {
291  //Add the weights as extra info of the mother
292  Particle* bremsGamma = bremsPair.second;
293  std::string extraInfoName = "bremsWeightWithPhoton" + std::to_string(photonIndex);
294  correctedLepton.addExtraInfo(extraInfoName, bremsPair.first);
295  photonIndex++;
296  bremsGammaEnergySum += Variable::eclClusterE(bremsGamma);
297 
298  const TMatrixFSym& gammaErrorMatrix = bremsGamma->getMomentumVertexErrorMatrix();
299  for (int irow = 0; irow <= 3; irow++) {
300  for (int icol = irow; icol <= 3; icol++) corLepMatrix(irow, icol) += gammaErrorMatrix(irow, icol);
301  }
302  correctedLepton.appendDaughter(bremsGamma, false);
303  B2DEBUG(10, "[BremsFinderModule] Found a bremsstrahlung gamma and added its 4-vector to the charged particle");
304  if (! m_addMultiplePhotons) break; //stop after adding the first photon
305  }
306 
307  correctedLepton.setMomentumVertexErrorMatrix(corLepMatrix);
308 
309  // add the info from original lepton to the new lepton
310  correctedLepton.setVertex(lepton->getVertex());
311  correctedLepton.setPValue(lepton->getPValue());
312  correctedLepton.addExtraInfo("bremsCorrected", float(selectedGammas.size() > 0));
313  correctedLepton.addExtraInfo("bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
314 
315  // add the mc relation
316  Particle* newLepton = particles.appendNew(correctedLepton);
317  const MCParticle* mcLepton = lepton->getRelated<MCParticle>();
318  const PIDLikelihood* pid = lepton->getPIDLikelihood();
319 
320  if (pid) newLepton->addRelationTo(pid);
321 
322  if (mcLepton != nullptr) newLepton->addRelationTo(mcLepton);
323 
324  m_outputList->addParticle(newLepton);
325 
326  } //Closes for loop on leptons
327 
328  } //Close event()
329 
330 
332 } // end Belle2 namespace
333 
Belle2::Particle::getPDGCode
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:380
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::Particle::getVertex
TVector3 getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:529
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::Particle::addExtraInfo
void addExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1224
Belle2::DataStore::EStoreFlags
EStoreFlags
Flags describing behaviours of objects etc.
Definition: DataStore.h:71
Belle2::DecayDescriptorParticle::getPDGCode
int getPDGCode() const
Return PDG code.
Definition: DecayDescriptorParticle.h:89
Belle2::RelationsInterface::addRelationTo
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).
Definition: RelationsObject.h:144
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::DecayDescriptor::init
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
Definition: DecayDescriptor.cc:47
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::Particle::getPIDLikelihood
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:796
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::RelationVector::begin
iterator begin()
Return iterator to first entry.
Definition: RelationVector.h:151
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Particle::getMomentumVertexErrorMatrix
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:386
Belle2::Particle::appendDaughter
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:633
Belle2::Particle::setMomentumVertexErrorMatrix
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:373
Belle2::Particle::setVertex
void setVertex(const TVector3 &vertex)
Sets position (decay vertex)
Definition: Particle.h:291
Belle2::BremsFinderModule
Bremsstrahlung finder correction module For each lepton in the input particle list,...
Definition: BremsFinderModule.h:53
Belle2::Particle::getTrack
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:770
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::RelationVector::weight
float weight(int index) const
Get weight with index.
Definition: RelationVector.h:120
Belle2::DecayDescriptor
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
Definition: DecayDescriptor.h:43
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::RelationVector::end
iterator end()
Return iterator to last entry +1.
Definition: RelationVector.h:153
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray< Particle >
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::DecayDescriptor::getMother
const DecayDescriptorParticle * getMother() const
return mother.
Definition: DecayDescriptor.h:136
Belle2::Particle::setPValue
void setPValue(float pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:308
Belle2::Particle::getPValue
float getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:565