Belle II Software  light-2205-abys
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 include
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/Vector4D.h>
27 
28 #include <vector>
29 
30 using namespace std;
31 using namespace Belle2;
32 
33 //-----------------------------------------------------------------
34 // Register module
35 //-----------------------------------------------------------------
36 
37 REG_MODULE(BelleBremRecovery);
38 
39 //-----------------------------------------------------------------
40 // Implementation
41 //-----------------------------------------------------------------
42 
43 BelleBremRecoveryModule::BelleBremRecoveryModule() :
44  Module(), m_pdgCode(0)
45 
46 {
47  // set module description (e.g. insert text)
48  setDescription(R"DOC(
49  Takes the charged particle from the given charged particle list (``inputListName``) and
50  copies it to the output list (``outputListName``). The 4-vector of the nearest (all) photon(s)
51  from ``gammaListName`` (considered as radiative) is added to the charged particle, if it is
52  found inside the cone around the charged particle with the given maximum angle (``angleThreshold``).
53  )DOC");
55 
56  // Add parameters
57  addParam("inputListName", m_inputListName,
58  "The initial charged particle list containing the charged particles to correct, should already exist.");
59  addParam("outputListName", m_outputListName, "The output charged particle list containing the corrected charged particles.");
60  addParam("gammaListName", m_gammaListName, "The gammas list containing possibly radiative gammas, should already exist.");
61  addParam("angleThreshold", m_angleThres,
62  "The maximum angle in radians between the charged particle and the (radiative) gamma to be accepted.", 0.05);
63  addParam("multiplePhotons", m_isMultiPho, "If only the nearest photon to add then make it False otherwise true", true);
64  addParam("writeOut", m_writeOut,
65  "Set to true if you want to write out the output list to a root file", false);
66 
67 }
68 
70 {
71  // check the validity of output ParticleList name
73  if (!valid)
74  B2ERROR("[BelleBremRecoveryModule] Invalid output ParticleList name: " << m_outputListName);
75 
76  // output particle
78  m_pdgCode = mother->getPDGCode();
80 
81  // get existing particle lists
83  B2ERROR("[BelleBremRecoveryModule] Input and output particle list names are the same: " << m_inputListName);
84  } else if (!m_decaydescriptor.init(m_inputListName)) {
85  B2ERROR("[BelleBremRecoveryModule] Invalid input particle list name: " << m_inputListName);
86  } else {
88  }
89 
91  B2ERROR("[BelleBremRecoveryModule] Invalid gamma particle list name: " << m_gammaListName);
92  } else {
93  m_gammaList.isRequired(m_gammaListName);
94  }
95 
96  // make output list
98  m_outputparticleList.registerInDataStore(m_outputListName, flags);
99  m_outputAntiparticleList.registerInDataStore(m_outputAntiListName, flags);
100 
101  m_particles.registerRelationTo(m_pidlikelihoods);
102 }
103 
104 
106 {
107  RelationArray particlesToMCParticles(m_particles, m_mcParticles);
108 
109  // new output particle list
110  m_outputparticleList.create();
112 
113  m_outputAntiparticleList.create();
115  m_outputAntiparticleList->bindAntiParticleList(*(m_outputparticleList));
116 
117  // loop over charged particles, correct them and add them to the output list
118  const unsigned int nLep = m_inputparticleList->getListSize();
119  const unsigned int nGam = m_gammaList->getListSize();
120  for (unsigned i = 0; i < nLep; i++) {
121  const Particle* lepton = m_inputparticleList->getParticle(i);
122  ROOT::Math::PxPyPzEVector new4Vec = lepton->get4Vector();
123  std::vector<Particle*> selectedGammas;
124  double bremsGammaEnergySum = 0.0;
125  // look for all possible (radiative) gamma
126  for (unsigned j = 0; j < nGam; j++) {
127  Particle* gamma = m_gammaList->getParticle(j);
128  // get angle (in lab system)
129  B2Vector3D pi = lepton->getMomentum();
130  B2Vector3D pj = gamma->getMomentum();
131  double angle = (pi.Angle(pj));
132  //Instead of first-come-first-serve, serve all charged particle equally.
133  // https://indico.belle2.org/event/946/contributions/4007/attachments/1946/2967/SCunliffe190827.pdf
134  if (m_angleThres > angle) {
135  gamma->writeExtraInfo("theta_e_gamma", angle);
136  selectedGammas.push_back(gamma);
137  }
138  }
139  //sorting the bremphotons in ascending order of the angle with the charged particle
140  std::sort(selectedGammas.begin(), selectedGammas.end(), [](const Particle * photon1, const Particle * photon2) {
141  return photon1->getExtraInfo("theta_e_gamma") < photon2->getExtraInfo("theta_e_gamma");
142  });
143  //Preparing 4-momentum vector of charged particle by adding bremphoton momenta
144  for (auto const& fsrgamma : selectedGammas) {
145  new4Vec += fsrgamma->get4Vector();
146  bremsGammaEnergySum += Variable::eclClusterE(fsrgamma);
147  if (!m_isMultiPho) break;
148  }
149  Particle correctedLepton(new4Vec, lepton->getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
150  lepton->getTrack()->getArrayIndex());
151  correctedLepton.appendDaughter(lepton, false);
152 
153  //Calculating matrix element of corrected charged particle
154  const TMatrixFSym& lepErrorMatrix = lepton->getMomentumVertexErrorMatrix();
155  TMatrixFSym corLepMatrix = lepErrorMatrix;
156  for (auto const& fsrgamma : selectedGammas) {
157  const TMatrixFSym& fsrErrorMatrix = fsrgamma->getMomentumVertexErrorMatrix();
158  for (int irow = 0; irow <= 3 ; irow++)
159  for (int icol = irow; icol <= 3; icol++)
160  corLepMatrix(irow, icol) += fsrErrorMatrix(irow, icol);
161  correctedLepton.appendDaughter(fsrgamma, false);
162  B2DEBUG(19, "[BelleBremRecoveryModule] Found a radiative gamma and added its 4-vector to the charge particle");
163  if (!m_isMultiPho) break;
164  }
165  correctedLepton.setMomentumVertexErrorMatrix(corLepMatrix);
166  correctedLepton.setVertex(lepton->getVertex());
167  correctedLepton.setPValue(lepton->getPValue());
168  correctedLepton.addExtraInfo("bremsCorrected", float(selectedGammas.size() > 0));
169  correctedLepton.addExtraInfo("bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
170  // add the mc relation
171  Particle* newLepton = m_particles.appendNew(correctedLepton);
172  const MCParticle* mcLepton = lepton->getRelated<MCParticle>();
173  const PIDLikelihood* pid = lepton->getPIDLikelihood();
174  if (pid)
175  newLepton->addRelationTo(pid);
176  if (mcLepton != nullptr) newLepton->addRelationTo(mcLepton);
177  m_outputparticleList->addParticle(newLepton);
178  }
179 }
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.
StoreArray< MCParticle > m_mcParticles
StoreArray of MCParticle objects.
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:34
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:26
Class to store reconstructed particles.
Definition: Particle.h:74
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:837
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1297
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:288
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:596
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:560
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:863
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:403
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:488
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1325
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:425
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:497
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:331
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:452
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:700
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
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.
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:23