Belle II Software  light-2205-abys
KlongMomentumCalculatorExpertModule.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/KlongDecayReconstructor/KlongMomentumCalculatorExpertModule.h>
11 
12 // framework aux
13 #include <framework/logging/Logger.h>
14 
15 // dataobjects
16 #include <analysis/dataobjects/Particle.h>
17 
18 // decay descriptor
19 #include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
20 
21 // utilities
22 #include <analysis/DecayDescriptor/ParticleListName.h>
23 #include <analysis/utility/ParticleCopy.h>
24 
25 #include <Math/Vector4D.h>
26 
27 #include <memory>
28 
29 using namespace std;
30 using namespace Belle2;
31 
32 //-----------------------------------------------------------------
33 // Register module
34 //-----------------------------------------------------------------
35 
36 REG_MODULE(KlongMomentumCalculatorExpert);
37 
38 //-----------------------------------------------------------------
39 // Implementation
40 //-----------------------------------------------------------------
41 
42 KlongMomentumCalculatorExpertModule::KlongMomentumCalculatorExpertModule() :
43  Module(), m_pdgCode(0)
44 
45 {
46  // set module description (e.g. insert text)
47  setDescription("This module is used to employ kinematic constraints to determine the momentum of Klongs for two body B decays containing a K_L0 and something else. The module creates a list of K_L0 candidates whose K_L0 momentum is reconstructed by combining the reconstructed direction (from either the ECL or KLM) of the K_L0 and kinematic constraints of the intial state.");
49 
50  // Add parameters
51  addParam("decayString", m_decayString,
52  "Input DecayDescriptor string.");
53  addParam("cut", m_cutParameter, "Selection criteria to be applied", std::string(""));
54  addParam("maximumNumberOfCandidates", m_maximumNumberOfCandidates,
55  "Don't reconstruct channel if more candidates than given are produced.", -1);
56  addParam("decayMode", m_decayModeID, "User-specified decay mode identifier (saved in 'decayModeID' extra-info for each Particle)",
57  0);
58  addParam("writeOut", m_writeOut,
59  "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
60  addParam("recoList", m_recoList,
61  "Suffix attached to the output K_L list, if not defined it is set to '_reco' \n", std::string("_reco"));
62 
63 }
64 
66 {
68 
69  // clear everything, initialize private members
70  m_pdgCode = 0;
71  m_listName = "";
73  m_generator = nullptr;
74 
75  // obtain the input and output particle lists from the decay string
76  bool valid = m_decaydescriptor.init(m_decayString);
77  if (!valid)
78  B2ERROR("Invalid input DecayString: " << m_decayString);
79 
80  // Mother particle
82 
83  m_pdgCode = mother->getPDGCode();
84  m_listName = mother->getFullName();
85 
88 
90 
91  // Daughters
92  bool k_check = false;
93  int nProducts = m_decaydescriptor.getNDaughters();
94  for (int i = 0; i < nProducts; ++i) {
95  const DecayDescriptorParticle* daughter =
97  StoreObjPtr<ParticleList>().isRequired(daughter->getFullName());
98  if (daughter->getPDGCode() == Const::Klong.getPDGCode()) {
99  m_klistName = daughter->getFullName() + m_klistName;
100  k_check = true;
101  }
102  }
103 
104  if (!k_check)
105  B2FATAL("This module is meant to reconstruct decays with a K_L0 in the final state. There is no K_L0 in this decay!");
106 
107  m_generator = std::make_unique<ParticleGenerator>(m_decayString, m_cutParameter);
108 
110  m_koutputList.registerInDataStore(m_klistName, flags);
111 
113 
114 }
115 
117 {
118  m_koutputList.create();
119  m_koutputList->initialize(Const::Klong.getPDGCode(), m_klistName);
120 
121  m_generator->init();
122 
123  int numberOfCandidates = 0;
124  while (m_generator->loadNext()) {
125 
126  Particle particle = m_generator->getCurrentParticle();
127 
128  Particle* kparticle = nullptr;
129 
130  bool is_physical = true;
131 
132  const std::vector<Particle*> daughters = particle.getDaughters();
133 
134  if (daughters.size() < 2)
135  B2FATAL("Reconstructing particle as a daughter of a decay with less then 2 daughters!");
136 
137  if (daughters.size() > 3)
138  B2FATAL("Higher multiplicity (>2) missing momentum decays not implemented yet!");
139 
140  ROOT::Math::PxPyPzEVector pDaughters;
141  for (auto daughter : daughters) {
142  if (daughter->getPDGCode() != Const::Klong.getPDGCode()) {
143  pDaughters += daughter->get4Vector();
144  }
145  }
146 
147  ROOT::Math::PxPyPzEVector klDaughters;
148  for (auto daughter : daughters) {
149  if (daughter->getPDGCode() == Const::Klong.getPDGCode()) {
150  kparticle = ParticleCopy::copyParticle(daughter);
151  klDaughters += daughter->get4Vector();
152  }
153  }
154 
155  double k_mag1 = 0.;
156  double k_mag2 = 0.;
157  double m_b = particle.getPDGMass();
158  double m_k = Const::Klong.getMass();
159  double m_j = 0;
160 
161  int idx = 0;
162 
163  for (auto daughter : daughters) {
164  if (daughter->getPDGCode() != Const::Klong.getPDGCode()) {
165  m_j = daughter->getPDGMass();
166  idx = daughter->getArrayIndex() + idx * 100;
167  }
168  }
169 
170  if (daughters.size() == 3) {
171  m_j = pDaughters.M();
172  }
173 
174  double s_p = (klDaughters.Vect().Unit()).Dot(pDaughters.Vect());
175  double m_sum = (m_b * m_b) - (m_j * m_j) - (m_k * m_k);
176  double e_j = pDaughters.E();
177 
178  double s_p2 = s_p * s_p;
179  double m_sum2 = m_sum * m_sum;
180  double s_pm = s_p * m_sum;
181  double e_j2 = e_j * e_j;
182  double m_k2 = m_k * m_k;
183 
184  k_mag1 = (s_pm + std::sqrt((s_p2) * (m_sum2) - 4 * ((e_j2) - (s_p2)) * ((e_j2) * (m_k2) - (m_sum2) / 4))) / (2 * (e_j2 - s_p2));
185  k_mag2 = (s_pm - std::sqrt((s_p2) * (m_sum2) - 4 * ((e_j2) - (s_p2)) * ((e_j2) * (m_k2) - (m_sum2) / 4))) / (2 * (e_j2 - s_p2));
186 
187  ROOT::Math::PxPyPzEVector missDaughters;
188 
189  if (k_mag1 > 0)
190  missDaughters = k_mag1 * klDaughters / klDaughters.P();
191  else
192  missDaughters = k_mag2 * klDaughters / klDaughters.P();
193  missDaughters.SetE(std::sqrt(m_k * m_k + missDaughters.P2()));
194 
195  for (auto daughter : daughters) {
196  if (daughter->getPDGCode() == Const::Klong.getPDGCode()) {
197  if (!isnan(missDaughters.P())) {
198  kparticle->set4Vector(missDaughters);
199  } else
200  is_physical = false;
201  }
202  }
203 
204  if (is_physical) {
205  ROOT::Math::PxPyPzEVector mom = pDaughters + missDaughters;
206  mom.SetE(std::sqrt(m_b * m_b + mom.P2()));
207  particle.set4Vector(mom);
208  if (isnan(mom.P()))
209  is_physical = false;
210  }
211 
212  if (!m_cut->check(&particle))
213  continue;
214 
215  if (!is_physical)
216  continue;
217 
218  numberOfCandidates++;
219 
220  if (m_maximumNumberOfCandidates > 0 and numberOfCandidates > m_maximumNumberOfCandidates) {
221  m_koutputList->clear();
222  break;
223  }
224 
225  m_koutputList->addParticle(kparticle);
226  kparticle->addExtraInfo("permID", idx);
227 
228  } //while
229 
230 } //event
int getPDGCode() const
PDG code.
Definition: Const.h:354
double getMass() const
Particle mass.
Definition: UnitConst.cc:323
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:558
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.
std::string getFullName() const
returns the full name of the particle full_name = name:label
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
int getNDaughters() const
return number of direct daughters.
const DecayDescriptor * getDaughter(int i) const
return i-th daughter (0 based index).
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Creates an instance of a cut and returns a unique_ptr to it, if you need a copy-able object instead y...
Definition: GeneralCut.h:84
bool m_isSelfConjugatedParticle
flag that indicates whether an anti-particle mother does not exist and should not be reconstructed as...
StoreObjPtr< ParticleList > m_koutputList
Klong output particle list.
int m_maximumNumberOfCandidates
drop all candidates if more candidates than this parameter are produced
std::string m_antiListName
output anti-particle list name
virtual void initialize() override
Initialize the Module.
std::string m_decayString
Input DecayString specifying the decay being reconstructed.
std::string m_recoList
Suffix attached to the output K_L list, if not defined it is set to '_reco'
std::unique_ptr< ParticleGenerator > m_generator
Generates the combinations.
std::unique_ptr< Variable::Cut > m_cut
cut object which performs the cuts
DecayDescriptor m_decaydescriptor
Decay descriptor of the decay being reconstructed.
int m_pdgCode
PDG code of the combined mother particle.
std::string m_klistName
output K_L0 particle list name
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 store reconstructed particles.
Definition: Particle.h:74
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:660
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:636
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 set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector.
Definition: Particle.h:276
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
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
Particle * copyParticle(const Particle *original)
Function takes argument Particle and creates a copy of it and copies of all its (grand-)^n-daughters.
Definition: ParticleCopy.cc:18
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