Belle II Software  light-2212-foldex
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/Vector3D.h>
27 #include <Math/Vector4D.h>
28 #include <Math/VectorUtil.h>
29 
30 #include <vector>
31 
32 using namespace std;
33 using namespace Belle2;
34 
35 //-----------------------------------------------------------------
36 // Register module
37 //-----------------------------------------------------------------
38 
39 REG_MODULE(BelleBremRecovery);
40 
41 //-----------------------------------------------------------------
42 // Implementation
43 //-----------------------------------------------------------------
44 
45 BelleBremRecoveryModule::BelleBremRecoveryModule() :
46  Module(), m_pdgCode(0)
47 
48 {
49  // set module description (e.g. insert text)
50  setDescription(R"DOC(
51  Takes the charged particle from the given charged particle list (``inputListName``) and
52  copies it to the output list (``outputListName``). The 4-vector of the nearest (all) photon(s)
53  from ``gammaListName`` (considered as radiative) is added to the charged particle, if it is
54  found 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("writeOut", m_writeOut,
67  "Set to true if you want to write out the output list to a root file", false);
68 
69 }
70 
72 {
73  // check the validity of output ParticleList name
75  if (!valid)
76  B2ERROR("[BelleBremRecoveryModule] Invalid output ParticleList name: " << m_outputListName);
77 
78  // output particle
80  m_pdgCode = mother->getPDGCode();
82 
83  // get existing particle lists
85  B2ERROR("[BelleBremRecoveryModule] Input and output particle list names are the same: " << m_inputListName);
86  } else if (!m_decaydescriptor.init(m_inputListName)) {
87  B2ERROR("[BelleBremRecoveryModule] Invalid input particle list name: " << m_inputListName);
88  } else {
90  }
91 
93  B2ERROR("[BelleBremRecoveryModule] Invalid gamma particle list name: " << m_gammaListName);
94  } else {
95  m_gammaList.isRequired(m_gammaListName);
96  }
97 
98  // make output list
100  m_outputparticleList.registerInDataStore(m_outputListName, flags);
101  m_outputAntiparticleList.registerInDataStore(m_outputAntiListName, flags);
102 
103  m_particles.registerRelationTo(m_pidlikelihoods);
104 }
105 
106 
108 {
109  RelationArray particlesToMCParticles(m_particles, m_mcParticles);
110 
111  // new output particle list
112  m_outputparticleList.create();
114 
115  m_outputAntiparticleList.create();
117  m_outputAntiparticleList->bindAntiParticleList(*(m_outputparticleList));
118 
119  // loop over charged particles, correct them and add them to the output list
120  const unsigned int nLep = m_inputparticleList->getListSize();
121  const unsigned int nGam = m_gammaList->getListSize();
122  for (unsigned i = 0; i < nLep; i++) {
123  const Particle* lepton = m_inputparticleList->getParticle(i);
124  ROOT::Math::PxPyPzEVector new4Vec = lepton->get4Vector();
125  std::vector<Particle*> selectedGammas;
126  double bremsGammaEnergySum = 0.0;
127  // look for all possible (radiative) gamma
128  for (unsigned j = 0; j < nGam; j++) {
129  Particle* gamma = m_gammaList->getParticle(j);
130  // get angle (in lab system)
131  ROOT::Math::XYZVector pi = lepton->getMomentum();
132  ROOT::Math::XYZVector pj = gamma->getMomentum();
133  double angle = ROOT::Math::VectorUtil::Angle(pi, pj);
134  //Instead of first-come-first-serve, serve all charged particle equally.
135  // https://indico.belle2.org/event/946/contributions/4007/attachments/1946/2967/SCunliffe190827.pdf
136  if (m_angleThres > angle) {
137  gamma->writeExtraInfo("theta_e_gamma", angle);
138  selectedGammas.push_back(gamma);
139  }
140  }
141  //sorting the bremphotons in ascending order of the angle with the charged particle
142  std::sort(selectedGammas.begin(), selectedGammas.end(), [](const Particle * photon1, const Particle * photon2) {
143  return photon1->getExtraInfo("theta_e_gamma") < photon2->getExtraInfo("theta_e_gamma");
144  });
145  //Preparing 4-momentum vector of charged particle by adding bremphoton momenta
146  for (auto const& fsrgamma : selectedGammas) {
147  new4Vec += fsrgamma->get4Vector();
148  bremsGammaEnergySum += Variable::eclClusterE(fsrgamma);
149  if (!m_isMultiPho) break;
150  }
151  Particle correctedLepton(new4Vec, lepton->getPDGCode(), Particle::EFlavorType::c_Flavored, Particle::c_Track,
152  lepton->getTrack()->getArrayIndex());
153  correctedLepton.appendDaughter(lepton, false);
154 
155  //Calculating matrix element of corrected charged particle
156  const TMatrixFSym& lepErrorMatrix = lepton->getMomentumVertexErrorMatrix();
157  TMatrixFSym corLepMatrix = lepErrorMatrix;
158  for (auto const& fsrgamma : selectedGammas) {
159  const TMatrixFSym& fsrErrorMatrix = fsrgamma->getMomentumVertexErrorMatrix();
160  for (int irow = 0; irow <= 3 ; irow++)
161  for (int icol = irow; icol <= 3; icol++)
162  corLepMatrix(irow, icol) += fsrErrorMatrix(irow, icol);
163  correctedLepton.appendDaughter(fsrgamma, false);
164  B2DEBUG(19, "[BelleBremRecoveryModule] Found a radiative gamma and added its 4-vector to the charge particle");
165  if (!m_isMultiPho) break;
166  }
167  correctedLepton.setMomentumVertexErrorMatrix(corLepMatrix);
168  correctedLepton.setVertex(lepton->getVertex());
169  correctedLepton.setPValue(lepton->getPValue());
170  correctedLepton.addExtraInfo("bremsCorrected", float(selectedGammas.size() > 0));
171  correctedLepton.addExtraInfo("bremsCorrectedPhotonEnergy", bremsGammaEnergySum);
172  // add the mc relation
173  Particle* newLepton = m_particles.appendNew(correctedLepton);
174  const MCParticle* mcLepton = lepton->getRelated<MCParticle>();
175  const PIDLikelihood* pid = lepton->getPIDLikelihood();
176  if (pid)
177  newLepton->addRelationTo(pid);
178  if (mcLepton != nullptr) newLepton->addRelationTo(mcLepton);
179  m_outputparticleList->addParticle(newLepton);
180  }
181 }
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: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: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:875
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:708
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1335
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:310
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:640
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:604
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:901
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:441
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:532
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1363
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:425
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:541
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:353
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:452
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