Belle II Software light-2406-ragdoll
NeutralHadron4MomentumCalculatorModule.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#include <analysis/modules/NeutralHadron4MomentumCalculator/NeutralHadron4MomentumCalculatorModule.h>
10
11#include <framework/logging/Logger.h>
12#include <analysis/dataobjects/Particle.h>
13#include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
14#include <analysis/utility/ParticleCopy.h>
15#include <vector>
16#include <cmath>
17
18using namespace Belle2;
19
20REG_MODULE(NeutralHadron4MomentumCalculator);
21
23{
24 // Set module properties
26 R"DOC(Calculates 4-momentum of a neutral hadron in a given decay chain e.g. B0 -> J/Psi K_L0, or anti-B0 -> p+ K- anti-n0.)DOC");
27
28 // Parameter definitions
29 addParam("decayString", m_decayString, "Decay string for which one wants to perform the calculation", std::string(""));
30 addParam("allowGamma", m_allowGamma, "Whether allow the selected particle to be gamma", false);
31 addParam("allowAnyParticleSource", m_allowAnyParticleSource, "Whether allow the selected particle to be from any ParticleSource",
32 false);
33
34}
35
37{
38 B2DEBUG(17, "Neutralhadron4MomentumCalculator: Use particle list: " << m_decayString);
40 if (!valid)
41 B2ERROR("NeutralHadron4MomentumCalculatorModule::initialize Invalid Decay Descriptor " << m_decayString);
42
44 if (hierarchy.size() != 1)
45 B2ERROR("NeutralHadron4MomentumCalculatorModule::initialize Only one particle can be selected in " << m_decayString);
46 if (hierarchy[0].size() != 2)
47 B2ERROR("NeutralHadron4MomentumCalculatorModule::initialize The selected particle must be a direct daughter " << m_decayString);
48
49 std::string neutralHadronName = hierarchy[0][1].second;
50 if (neutralHadronName != "n0" and neutralHadronName != "K_L0") {
51 if (m_allowGamma and neutralHadronName == "gamma")
52 B2WARNING("NeutralHadron4MomentumCalculatorModule::initialize The selected particle is gamma but you allowed so; be aware.");
54 B2WARNING("NeutralHadron4momentumCalculatorModule::initialize The selected particle can be from anything; the magnitude of "
55 "the momentum will be overridden by that calculated from the mother mass constraint; be aware.");
56 else
57 B2ERROR("NeutralHadron4MomentumCalculatorModule::initialize The selected particle must be a long-lived neutral hadron "
58 "i.e. (anti-)n0 or K_L0, unless you set allowGamma to true and selected a photon (gamma), or you set "
59 "allowAnyParticleSource to true."
60 "Input particle: " << m_decayString);
61 }
62
63 m_iNeutral = hierarchy[0][1].first;
64
65 std::string motherFullName = m_decayDescriptor.getMother()->getFullName();
66 m_plist.isRequired(motherFullName);
67}
68
70{
71 const unsigned int n = m_plist->getListSize();
72 std::vector<unsigned int> toRemove;
73 for (unsigned i = 0; i < n; i++) {
74 Particle* particle = m_plist->getParticle(i);
75 std::vector<Particle*> daughters = particle->getDaughters();
76 ROOT::Math::PxPyPzEVector others4Momentum;
77 for (int j = 0; j < m_decayDescriptor.getNDaughters(); j++) {
78 if (j != m_iNeutral) {
79 others4Momentum += daughters[j]->get4Vector();
80 }
81 }
82 const Particle* originalNeutral = daughters[m_iNeutral];
83 Particle* neutral = ParticleCopy::copyParticle(originalNeutral);
84 particle->replaceDaughter(originalNeutral, neutral);
85 ROOT::Math::XYZVector neutralDirection = neutral->getMomentum().Unit();
86 double a = others4Momentum.Vect().Dot(neutralDirection);
87 double b = (std::pow(particle->getPDGMass(), 2) - std::pow(neutral->getMass(), 2) - others4Momentum.mag2()) / 2.;
88 double c = others4Momentum.E();
89 double d = std::pow(neutral->getMass(), 2);
90 double D = (a * a - c * c) * d + b * b;
91 if (D >= 0) {
92 double neutralP = (-1. * a * b - c * std::sqrt(D)) / (a * a - c * c);
93 double neutralPx = neutralP * neutralDirection.X();
94 double neutralPy = neutralP * neutralDirection.Y();
95 double neutralPz = neutralP * neutralDirection.Z();
96 double neutralE = std::sqrt(neutralP * neutralP + d);
97 const ROOT::Math::PxPyPzEVector newNeutral4Momentum(neutralPx, neutralPy, neutralPz, neutralE);
98
99 double motherPx = newNeutral4Momentum.Px() + others4Momentum.Px();
100 double motherPy = newNeutral4Momentum.Py() + others4Momentum.Py();
101 double motherPz = newNeutral4Momentum.Pz() + others4Momentum.Pz();
102 double motherE = newNeutral4Momentum.E() + others4Momentum.E();
103 const ROOT::Math::PxPyPzEVector newMother4Momentum(motherPx, motherPy, motherPz, motherE);
104
105 neutral->set4Vector(newNeutral4Momentum);
106 particle->set4Vector(newMother4Momentum);
107 } else {
108 toRemove.push_back(particle->getArrayIndex());
109 }
110 }
111 m_plist->removeParticles(toRemove);
112}
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.
int getNDaughters() const
return number of direct daughters.
std::vector< std::vector< std::pair< int, std::string > > > getHierarchyOfSelected()
Function to get hierarchy of selected particles and their names (for python use)
const DecayDescriptorParticle * getMother() const
return mother.
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
bool m_allowAnyParticleSource
Whether allow the selected particle to be from any ParticleSource.
void event() override
Action to perform for each event.
std::string m_decayString
Decay string for which one wants to perform the calculation.
int m_iNeutral
Index of the neutral hadron in the decay string.
NeutralHadron4MomentumCalculatorModule()
Constructor: Sets the description, the properties and the parameters of the module.
DecayDescriptor m_decayDescriptor
Decay Descriptor to be initialized with m_decayString.
bool m_allowGamma
Whether allow the selected particle to be gamma.
StoreObjPtr< ParticleList > m_plist
ParticleList that one wants to perform the calculation.
Class to store reconstructed particles.
Definition: Particle.h:75
bool replaceDaughter(const Particle *oldDaughter, Particle *newDaughter)
Replace index of given daughter with new daughter, return true if a replacement is made.
Definition: Particle.cc:705
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:637
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:604
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector.
Definition: Particle.h:271
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
double getMass() const
Returns invariant mass (= nominal for FS particles)
Definition: Particle.h:507
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
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
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24