Belle II Software light-2406-ragdoll
LowEnergyPi0VetoExpertModule.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 header. */
10#include <analysis/modules/LowEnergyPi0VetoExpert/LowEnergyPi0VetoExpertModule.h>
11
12/* Analysis headers. */
13#include <analysis/variables/ECLVariables.h>
14#include <analysis/variables/HelicityVariables.h>
15
16/* Basf2 headers. */
17#include <mva/interface/Interface.h>
18
19/* Boost headers. */
20#include <boost/algorithm/string/predicate.hpp>
21
22/* ROOT headers. */
23#include <Math/Vector3D.h>
24#include <Math/Vector4D.h>
25#include <Math/VectorUtil.h>
26
27using namespace Belle2;
28
29REG_MODULE(LowEnergyPi0VetoExpert);
30
32{
33 setDescription("Low-energy pi0 veto.");
35 "VetoPi0Daughters", m_VetoPi0Daughters,
36 "Veto for pi0 daughters (maximum over all pairs excluding this pi0).",
37 false);
38 addParam("GammaListName", m_GammaListName, "Gamma particle list name.",
39 std::string("gamma"));
40 addParam("Pi0ListName", m_Pi0ListName, "Pi0 particle list name.",
41 std::string("pi0"));
42 addParam("Belle1", m_Belle1, "Belle 1 data analysis.", false);
43 addParam("identifier", m_identifier,
44 "Database identifier or file used to load the weights.",
47}
48
50{
51}
52
54{
55 m_ListGamma.isRequired(m_GammaListName);
57 m_ListPi0.isRequired(m_Pi0ListName);
58 if (not(boost::ends_with(m_identifier, ".root") or boost::ends_with(m_identifier, ".xml"))) {
59 m_weightfile_representation = std::unique_ptr<DBObjPtr<DatabaseRepresentationOfWeightfile>>(new
61 }
63}
64
66{
67 m_expert.reset();
68 m_dataset.reset();
69}
70
72{
74 if (m_weightfile_representation->hasChanged()) {
75 std::stringstream ss((*m_weightfile_representation)->m_data);
76 auto weightfile = MVA::Weightfile::loadFromStream(ss);
77 init_mva(weightfile);
78 }
79 } else {
81 init_mva(weightfile);
82 }
83}
84
86{
87}
88
90{
91 auto supported_interfaces = MVA::AbstractInterface::getSupportedInterfaces();
92 MVA::GeneralOptions general_options;
93 weightfile.getOptions(general_options);
94 weightfile.addSignalFraction(0.5);
95 m_expert = supported_interfaces[general_options.m_method]->getExpert();
96 m_expert->load(weightfile);
97 std::vector<float> dummy;
98 /* The number of input variables depends on the experiment. */
99 int nInputVariables;
100 if (m_Belle1)
101 nInputVariables = 7;
102 else
103 nInputVariables = 11;
104 dummy.resize(nInputVariables, 0);
105 m_dataset = std::unique_ptr<MVA::SingleDataset>(new MVA::SingleDataset(general_options, std::move(dummy), 0));
106}
107
109 const Particle* pi0Gamma)
110{
111 float maxVeto = 0;
112 int n = m_ListGamma->getListSize();
113 for (int i = 0; i < n; ++i) {
114 Particle* gamma2 = m_ListGamma->getParticle(i);
115 if (gamma1 == gamma2)
116 continue;
117 if (pi0Gamma != nullptr) {
118 if (pi0Gamma == gamma2)
119 continue;
120 }
121 double pi0Mass = (gamma1->get4Vector() + gamma2->get4Vector()).M();
122 if (pi0Mass < 0.07 || pi0Mass > 0.20)
123 continue;
124 const Particle* gammaLowEnergy, *gammaHighEnergy;
125 if (gamma1->getEnergy() > gamma2->getEnergy()) {
126 gammaLowEnergy = gamma2;
127 gammaHighEnergy = gamma1;
128 } else {
129 gammaLowEnergy = gamma1;
130 gammaHighEnergy = gamma2;
131 }
132 double gammaLowEnergyEnergy, gammaHighEnergyEnergy;
133 double gammaLowEnergyE9E21, gammaHighEnergyE9E21;
134 double gammaLowEnergyClusterTheta, gammaHighEnergyClusterTheta;
135 double gammaLowEnergyZernikeMVA, gammaHighEnergyZernikeMVA;
136 double gammaLowEnergyIsolation, gammaHighEnergyIsolation;
137 double cosHelicityAngleMomentum;
138 gammaLowEnergyEnergy = gammaLowEnergy->getEnergy();
139 gammaHighEnergyEnergy = gammaHighEnergy->getEnergy();
140 ROOT::Math::PxPyPzEVector gammaHighEnergyMomentum(
141 gammaHighEnergy->getPx(), gammaHighEnergy->getPy(),
142 gammaHighEnergy->getPz(), gammaHighEnergyEnergy);
143 ROOT::Math::PxPyPzEVector gammaLowEnergyMomentum(
144 gammaLowEnergy->getPx(), gammaLowEnergy->getPy(),
145 gammaLowEnergy->getPz(), gammaLowEnergyEnergy);
146 ROOT::Math::PxPyPzEVector momentum = gammaHighEnergyMomentum +
147 gammaLowEnergyMomentum;
148 ROOT::Math::XYZVector boost = momentum.BoostToCM();
149 gammaHighEnergyMomentum =
150 ROOT::Math::VectorUtil::boost(gammaHighEnergyMomentum, boost);
151 cosHelicityAngleMomentum =
152 fabs(ROOT::Math::VectorUtil::CosTheta(momentum.Vect(),
153 gammaHighEnergyMomentum.Vect()));
154 gammaLowEnergyE9E21 = Variable::eclClusterE9E21(gammaLowEnergy);
155 gammaHighEnergyE9E21 = Variable::eclClusterE9E21(gammaHighEnergy);
156 gammaLowEnergyClusterTheta = Variable::eclClusterTheta(gammaLowEnergy);
157 gammaHighEnergyClusterTheta = Variable::eclClusterTheta(gammaHighEnergy);
158 if (!m_Belle1) {
159 gammaLowEnergyZernikeMVA =
160 Variable::eclClusterZernikeMVA(gammaLowEnergy);
161 gammaHighEnergyZernikeMVA =
162 Variable::eclClusterZernikeMVA(gammaHighEnergy);
163 gammaLowEnergyIsolation = Variable::eclClusterIsolation(gammaLowEnergy);
164 gammaHighEnergyIsolation =
165 Variable::eclClusterIsolation(gammaHighEnergy);
166 }
167 m_dataset->m_input[0] = gammaLowEnergyEnergy;
168 m_dataset->m_input[1] = pi0Mass;
169 m_dataset->m_input[2] = cosHelicityAngleMomentum;
170 m_dataset->m_input[3] = gammaLowEnergyE9E21;
171 m_dataset->m_input[4] = gammaHighEnergyE9E21;
172 m_dataset->m_input[5] = gammaLowEnergyClusterTheta;
173 m_dataset->m_input[6] = gammaHighEnergyClusterTheta;
174 if (!m_Belle1) {
175 m_dataset->m_input[7] = gammaLowEnergyZernikeMVA;
176 m_dataset->m_input[8] = gammaHighEnergyZernikeMVA;
177 m_dataset->m_input[9] = gammaLowEnergyIsolation;
178 m_dataset->m_input[10] = gammaHighEnergyIsolation;
179 }
180 float veto = m_expert->apply(*m_dataset)[0];
181 if (veto > maxVeto)
182 maxVeto = veto;
183 }
184 return maxVeto;
185}
186
188{
189 if (m_VetoPi0Daughters) {
190 int n = m_ListPi0->getListSize();
191 for (int i = 0; i < n; ++i) {
192 Particle* pi0 = m_ListPi0->getParticle(i);
193 const Particle* gamma1 = pi0->getDaughter(0);
194 const Particle* gamma2 = pi0->getDaughter(1);
195 const Particle* gammaLowEnergy, *gammaHighEnergy;
196 if (gamma1->getEnergy() > gamma2->getEnergy()) {
197 gammaLowEnergy = gamma2;
198 gammaHighEnergy = gamma1;
199 } else {
200 gammaLowEnergy = gamma1;
201 gammaHighEnergy = gamma2;
202 }
203 float maxVeto = getMaximumVeto(gammaLowEnergy, gammaHighEnergy);
204 pi0->addExtraInfo("lowEnergyPi0VetoGammaLowEnergy", maxVeto);
205 maxVeto = getMaximumVeto(gammaHighEnergy, gammaLowEnergy);
206 pi0->addExtraInfo("lowEnergyPi0VetoGammaHighEnergy", maxVeto);
207 }
208 } else {
209 int n = m_ListGamma->getListSize();
210 for (int i = 0; i < n; ++i) {
211 Particle* gamma = m_ListGamma->getParticle(i);
212 float maxVeto = getMaximumVeto(gamma, nullptr);
213 gamma->addExtraInfo("lowEnergyPi0Veto", maxVeto);
214 }
215 }
216}
Class for accessing objects in the database.
Definition: DBObjPtr.h:21
std::unique_ptr< MVA::SingleDataset > m_dataset
Pointer to the current dataset.
StoreObjPtr< ParticleList > m_ListGamma
Gamma candidates.
void event() override
This method is called for each event.
void endRun() override
This method is called if the current run ends.
void terminate() override
This method is called at the end of the event processing.
std::unique_ptr< MVA::Expert > m_expert
Pointer to the current MVA expert.
std::unique_ptr< DBObjPtr< DatabaseRepresentationOfWeightfile > > m_weightfile_representation
Database pointer to the database representation of the weightfile.
void beginRun() override
Called when entering a new run.
bool m_VetoPi0Daughters
Calculate veto for pi0 daughter photons (maximum over all pairs excluding this pi0).
StoreObjPtr< ParticleList > m_ListPi0
Pi0 candidates.
void init_mva(MVA::Weightfile &weightfile)
Initialize mva expert, dataset and features Called every time the weightfile in the database changes ...
std::string m_GammaListName
Gamma particle list name.
std::string m_Pi0ListName
Pi0 particle list name.
float getMaximumVeto(const Particle *gamma1, const Particle *pi0Gamma)
Get maximum veto value over all gamma pairs including the photon gamma1.
std::string m_identifier
Database identifier or file used to load the weights.
static void initSupportedInterfaces()
Static function which initliazes all supported interfaces, has to be called once before getSupportedI...
Definition: Interface.cc:45
static std::map< std::string, AbstractInterface * > getSupportedInterfaces()
Returns interfaces supported by the MVA Interface.
Definition: Interface.h:53
General options which are shared by all MVA trainings.
Definition: Options.h:62
Wraps the data of a single event into a Dataset.
Definition: Dataset.h:135
The Weightfile class serializes all information about a training into an xml tree.
Definition: Weightfile.h:38
static Weightfile loadFromStream(std::istream &stream)
Static function which deserializes a Weightfile from a stream.
Definition: Weightfile.cc:251
void getOptions(Options &options) const
Fills an Option object from the xml tree.
Definition: Weightfile.cc:67
static Weightfile loadFromFile(const std::string &filename)
Static function which loads a Weightfile from a file.
Definition: Weightfile.cc:206
void addSignalFraction(float signal_fraction)
Saves the signal fraction in the xml tree.
Definition: Weightfile.cc:95
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:75
double getPx() const
Returns x component of momentum.
Definition: Particle.h:587
double getPz() const
Returns z component of momentum.
Definition: Particle.h:605
double getEnergy() const
Returns total energy.
Definition: Particle.h:535
double getPy() const
Returns y component of momentum.
Definition: Particle.h:596
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1336
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
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
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24