9#include <analysis/modules/ChargedParticleIdentificator/ChargedPidMVAModule.h>
12#include <analysis/dataobjects/Particle.h>
13#include <analysis/dataobjects/ParticleList.h>
14#include <analysis/dbobjects/ChargedPidMVAWeights.h>
15#include <analysis/DecayDescriptor/DecayDescriptor.h>
16#include <analysis/VariableManager/Utility.h>
17#include <analysis/variables/ECLVariables.h>
20#include <framework/logging/LogConfig.h>
21#include <framework/logging/LogSystem.h>
22#include <mva/interface/Interface.h>
30 setDescription(
"This module evaluates the response of an MVA trained for binary charged particle identification between two hypotheses, S and B. For a given input set of (S,B) mass hypotheses, it takes the Particle objects in the appropriate charged stable particle's ParticleLists, calculates the MVA score using the appropriate xml weight file, and adds it as ExtraInfo to the Particle objects.");
32 setPropertyFlags(c_ParallelProcessingCertified);
34 addParam(
"sigHypoPDGCode",
36 "The input signal mass hypothesis' pdgId.",
38 addParam(
"bkgHypoPDGCode",
40 "The input background mass hypothesis' pdgId.",
42 addParam(
"particleLists",
44 "The input list of DecayStrings, where each selected (^) daughter should correspond to a standard charged ParticleList, e.g. ['Lambda0:sig -> ^p+ ^pi-', 'J/psi:sig -> ^mu+ ^mu-']. One can also directly pass a list of standard charged ParticleLists, e.g. ['e+:my_electrons', 'pi+:my_pions']. Note that charge-conjugated ParticleLists will automatically be included.",
45 std::vector<std::string>());
46 addParam(
"payloadName",
48 "The name of the database payload object with the MVA weights.",
49 std::string(
"ChargedPidMVAWeights"));
50 addParam(
"chargeIndependent",
52 "Specify whether to use a charge-independent training of the MVA.",
54 addParam(
"useECLOnlyTraining",
56 "Specify whether to use an ECL-only training of the MVA.",
72 " of the signal mass hypothesis is not that of a valid particle in Const::chargedStableSet! Aborting...");
76 " of the background mass hypothesis is not that of a valid particle in Const::chargedStableSet! Aborting...");
104 std::map<int, std::string> debugStr = {
113 DecayDescriptor decayDescriptor;
114 decayDescriptor.init(decayString);
115 auto pListName = decayDescriptor.getMother()->getFullName();
117 unsigned short m_nSelectedDaughters = decayDescriptor.getSelectionNames().size();
118 StoreObjPtr<ParticleList> pList(pListName);
121 B2FATAL(
"ParticleList: " << pListName <<
" could not be found. Aborting...");
124 auto pListSize = pList->getListSize();
126 B2DEBUG(11,
"ParticleList: " << pList->getParticleListName() <<
" - N = " << pListSize <<
" particles.");
128 const auto nTargetParticles = (m_nSelectedDaughters == 0) ? pListSize : pListSize * m_nSelectedDaughters;
131 std::vector<int> pdgs;
132 if (m_nSelectedDaughters == 0) {
133 pdgs.push_back(pList->getPDGCode());
135 pdgs = decayDescriptor.getSelectionPDGCodes();
137 for (
auto pdg : pdgs) {
140 B2FATAL(
"PDG: " << pdg <<
" of ParticleList: " << pListName <<
141 " is not that of a valid particle in Const::chargedStableSet! Aborting...");
144 std::vector<const Particle*> targetParticles;
145 if (m_nSelectedDaughters > 0) {
146 for (
unsigned int iPart(0); iPart < pListSize; ++iPart) {
147 auto* iParticle = pList->getParticle(iPart);
148 auto daughters = decayDescriptor.getSelectionParticles(iParticle);
149 for (
auto* iDaughter : daughters) {
150 targetParticles.push_back(iDaughter);
155 for (
unsigned int ipart(0); ipart < nTargetParticles; ++ipart) {
157 const Particle* particle = (m_nSelectedDaughters == 0) ? pList->getParticle(ipart) : targetParticles[ipart];
161 const ECLCluster* eclCluster = particle->getECLCluster();
163 B2DEBUG(11,
"\nParticle [" << ipart <<
"] has invalid Track-ECLCluster relation, skip MVA application...");
172 auto p = particle->getP();
175 if (std::isnan(theta) or std::isnan(p) or std::isnan(charge)) {
176 B2DEBUG(11,
"\nParticle [" << ipart <<
"] has invalid input variable, skip MVA application..." <<
177 " polar angle: " << theta <<
", p: " << p <<
", charge: " << charge);
181 int idx_theta, idx_p, idx_charge;
184 auto hasMatch = std::isnormal(Variable::eclClusterTrackMatched(particle));
186 debugStr[11] +=
"\n";
187 debugStr[11] += (
"Particle [" + std::to_string(ipart) +
"]\n");
188 debugStr[11] += (
"Has ECL cluster match? " + std::to_string(hasMatch) +
"\n");
189 debugStr[11] += (
"polar angle: " + thVarName +
" = " + std::to_string(theta) +
" [rad]\n");
190 debugStr[11] += (
"p = " + std::to_string(p) +
" [GeV/c]\n");
192 debugStr[11] += (
"charge = " + std::to_string(charge) +
"\n");
194 debugStr[11] += (
"Is brems corrected ? " + std::to_string(particle->hasExtraInfo(
"bremsCorrected")) +
"\n");
195 debugStr[11] += (
"Weightfile idx = " + std::to_string(index) +
" - (polar angle, p, charge) = (" + std::to_string(
196 idx_theta) +
", " + std::to_string(idx_p) +
", " +
197 std::to_string(idx_charge) +
")\n");
199 debugStr[11] += (
"Category cut: " +
m_cuts.at(index)->decompile() +
"\n");
202 B2DEBUG(11, debugStr[11]);
203 debugStr[11].clear();
207 if (!
m_cuts.at(index)->check(particle)) {
208 B2DEBUG(11,
"\nParticle [" << ipart <<
"] didn't pass MVA category cut, skip MVA application...");
215 debugStr[11] +=
"\n";
216 debugStr[11] +=
"MVA variables:\n";
219 for (
unsigned int ivar(0); ivar < nvars; ++ivar) {
223 double var = std::numeric_limits<double>::quiet_NaN();
224 auto var_result = varobj->function(particle);
225 if (std::holds_alternative<double>(var_result)) {
226 var = std::get<double>(var_result);
227 }
else if (std::holds_alternative<int>(var_result)) {
228 var = std::get<int>(var_result);
229 }
else if (std::holds_alternative<bool>(var_result)) {
230 var = std::get<bool>(var_result);
232 B2ERROR(
"Variable '" << varobj->name <<
"' has wrong data type! It must be one of double, integer, or bool.");
237 var = (std::isnan(var)) ? -999.0 : var;
240 debugStr[11] += (
"\tvar[" + std::to_string(ivar) +
"] : " + varobj->name +
" = " + std::to_string(var) +
"\n");
246 B2DEBUG(11, debugStr[11]);
247 debugStr[11].clear();
252 debugStr[12] +=
"\n";
253 debugStr[12] +=
"MVA spectators:\n";
256 for (
unsigned int ispec(0); ispec < nspecs; ++ispec) {
260 double spec = std::numeric_limits<double>::quiet_NaN();
261 auto spec_result = specobj->function(particle);
262 if (std::holds_alternative<double>(spec_result)) {
263 spec = std::get<double>(spec_result);
264 }
else if (std::holds_alternative<int>(spec_result)) {
265 spec = std::get<int>(spec_result);
266 }
else if (std::holds_alternative<bool>(spec_result)) {
267 spec = std::get<bool>(spec_result);
269 B2ERROR(
"Variable '" << specobj->name <<
"' has wrong data type! It must be one of double, integer, or bool.");
272 debugStr[12] += (
"\tspec[" + std::to_string(ispec) +
"] : " + specobj->name +
" = " + std::to_string(spec) +
"\n");
274 m_datasets.at(index)->m_spectators[ispec] = spec;
278 B2DEBUG(12, debugStr[12]);
279 debugStr[12].clear();
285 debugStr[11] +=
"\n";
286 debugStr[12] +=
"\n";
287 debugStr[11] +=
"MVA response:\n";
291 debugStr[11] += (
"\tscore = " + std::to_string(score));
297 B2DEBUG(11, debugStr[11]);
298 B2DEBUG(12, debugStr[12]);
299 debugStr[11].clear();
300 debugStr[12].clear();
315 std::map<std::string, std::string> aliasesLegacy;
317 aliasesLegacy.insert(std::make_pair(
"__event__",
"evtNum"));
324 aliasesLegacy.insert(std::make_pair(
"missingLogL_" + detName,
"pidMissingProbabilityExpert(" + detName +
")"));
328 std::string alias =
"deltaLogL_" + std::get<0>(info) +
"_" + std::get<1>(info) +
"_" + detName;
329 std::string var =
"pidDeltaLogLikelihoodValueExpert(" + std::to_string(pdgId) +
", " + std::to_string(std::get<2>
330 (info)) +
"," + detName +
")";
332 aliasesLegacy.insert(std::make_pair(alias, var));
334 if (it.getIndex() == 0) {
335 alias =
"deltaLogL_" + std::get<0>(info) +
"_" + std::get<1>(info) +
"_ALL";
336 var =
"pidDeltaLogLikelihoodValueExpert(" + std::to_string(pdgId) +
", " + std::to_string(std::get<2>(info)) +
", ALL)";
337 aliasesLegacy.insert(std::make_pair(alias, var));
344 B2INFO(
"Setting hard-coded aliases for the ChargedPidMVA algorithm.");
346 std::string debugStr(
"\n");
347 for (
const auto& [alias, variable] : aliasesLegacy) {
348 debugStr += (alias +
" --> " + variable +
"\n");
350 B2ERROR(
"Something went wrong with setting alias: " << alias <<
" for variable: " << variable);
353 B2DEBUG(10, debugStr);
363 if (!aliases->empty()) {
365 B2INFO(
"Setting aliases for the ChargedPidMVA algorithm read from the payload.");
367 std::string debugStr(
"\n");
368 for (
const auto& [alias, variable] : *aliases) {
369 if (alias != variable) {
370 debugStr += (alias +
" --> " + variable +
"\n");
372 B2ERROR(
"Something went wrong with setting alias: " << alias <<
" for variable: " << variable);
376 B2DEBUG(10, debugStr);
391 B2INFO(
"Run: " <<
m_event_metadata->getRun() <<
". Load supported MVA interfaces for binary charged particle identification...");
400 B2INFO(
"\tLoading weightfiles from the payload class for SIGNAL particle hypothesis: " <<
m_sig_pdg);
403 auto nfiles = serialized_weightfiles->size();
405 B2INFO(
"\tConstruct the MVA experts and datasets from N = " << nfiles <<
" weightfiles...");
415 for (
unsigned int idx(0); idx < nfiles; idx++) {
417 B2DEBUG(12,
"\t\tweightfile[" << idx <<
"]");
420 std::stringstream ss(serialized_weightfiles->at(idx));
423 MVA::GeneralOptions general_options;
424 weightfile.getOptions(general_options);
428 m_variables[idx] = manager.getVariables(general_options.m_variables);
429 m_spectators[idx] = manager.getVariables(general_options.m_spectators);
431 B2DEBUG(12,
"\t\tRetrieved N = " << general_options.m_variables.size()
432 <<
" variables, N = " << general_options.m_spectators.size()
436 m_experts[idx] = supported_interfaces[general_options.m_method]->getExpert();
439 B2DEBUG(12,
"\t\tweightfile loaded successfully into expert[" << idx <<
"]!");
442 std::vector<float> v(general_options.m_variables.size(), 0.0);
443 std::vector<float> s(general_options.m_spectators.size(), 0.0);
444 m_datasets[idx] = std::make_unique<MVA::SingleDataset>(general_options, v, 1.0, s);
446 B2DEBUG(12,
"\t\tdataset[" << idx <<
"] created successfully!");
450 const auto cutstr = (!cuts->empty()) ? cuts->at(idx) :
"";
453 B2DEBUG(12,
"\t\tcut[" << idx <<
"] created successfully!");
StoreObjPtr< EventMetaData > m_event_metadata
The event information.
int m_bkg_pdg
The input background mass hypothesis' pdgId.
std::vector< std::string > m_decayStrings
The input list of DecayStrings, where each selected (^) daughter should correspond to a standard char...
bool m_ecl_only
Flag to specify if we use an ECL-only based training.
std::string m_score_varname
The lookup name of the MVA score variable, given the input S, B mass hypotheses for the algorithm.
std::unique_ptr< DBObjPtr< ChargedPidMVAWeights > > m_weightfiles_representation
Interface to get the database payload with the MVA weight files.
StoreArray< Particle > m_particles
StoreArray of Particles.
DatasetsList m_datasets
List of MVA::SingleDataset objects.
std::map< int, std::tuple< std::string, std::string, int > > m_stdChargedInfo
Map with standard charged particles' info.
bool m_charge_independent
Flag to specify if we use a charge-independent training.
virtual void event() override
Called once for each event.
void registerAliases()
Set variable aliases needed by the MVA.
virtual void initialize() override
Use this to initialize resources or memory your module needs.
int m_sig_pdg
The input signal mass hypothesis' pdgId.
virtual void beginRun() override
Called once before a new run begins.
VariablesLists m_variables
List of lists of feature variables.
void registerAliasesLegacy()
Set variable aliases needed by the MVA.
ChargedPidMVAModule()
Constructor, for setting module description and parameters.
CutsList m_cuts
List of Cut objects.
VariablesLists m_spectators
List of lists of spectator variables.
virtual ~ChargedPidMVAModule()
Destructor, use this to clean up anything you created in the constructor.
ExpertsList m_experts
List of MVA::Expert objects.
std::string m_payload_name
The name of the database payload object with the MVA weights.
Iterator end() const
Ending iterator.
EDetector
Enum for identifying the detector components (detector and subdetector).
static std::string parseDetectors(EDetector det)
Converts Const::EDetector object to string.
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
@ c_Debug
Debug: for code development.
static LogSystem & Instance()
Static method to get a reference to the LogSystem instance.
static void initSupportedInterfaces()
Static function which initliazes all supported interfaces, has to be called once before getSupportedI...
static std::map< std::string, AbstractInterface * > getSupportedInterfaces()
Returns interfaces supported by the MVA Interface.
static Weightfile loadFromStream(std::istream &stream)
Static function which deserializes a Weightfile from a stream.
static Manager & Instance()
get singleton instance.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Abstract base class for different kinds of events.