9#include <analysis/modules/TrackIsoCalculator/TrackIsoCalculatorModule.h>
11#include <analysis/dbobjects/PIDDetectorWeights.h>
12#include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
13#include <analysis/VariableManager/Manager.h>
14#include <analysis/utility/DetectorSurface.h>
16#include <framework/database/DBObjPtr.h>
19#include <boost/algorithm/string.hpp>
30 R
"DOC(Calculate track isolation variables on the charged stable particles, or selected charged daughters, of the input ParticleList.)DOC");
36 "The name of the input charged stable particle list, e.g. ``mu+:all``, or a composite particle w/ charged stable daughters for which distances are to be calculated, e.g. ``D0 -> ^K- pi+``. Note that in the latter case we allow only one daughter to be selected in the decay string per module instance.");
39 "The name of the input ParticleList of reference tracks. Must be a charged stable particle as defined in Const::chargedStableSet.");
42 "The list of names of the detectors at whose (cylindrical) surface(s) we extrapolate each helix's polar and azimuthal angle. Allowed values: {CDC, TOP, ARICH, ECL, KLM}.",
43 std::vector<std::string>());
46 "If this option is set, the helix extrapolation for the target and reference particles will use the track fit result for the most probable mass hypothesis, namely, the one that gives the highest chi2Prob of the fit.",
50 "The name of the database payload object with the PID detector weights.",
51 std::string(
"PIDDetectorWeights"));
54 "If set to true, will not use the PID detector weights for the score definition.",
81 B2ERROR(
"More than one daughter is selected in the decay string " <<
m_decayString <<
".");
86 B2INFO(
"Calculating track-based isolation variables for the decay string: "
88 <<
" using the reference ParticleList: "
93 B2ERROR(
"Selected ParticleList in decay string:"
95 <<
" and/or reference ParticleList: "
97 <<
" is not that of a valid particle in Const::chargedStableSet!");
101 B2INFO(
"Will use track fit result for the most probable mass hypothesis in helix extrapolation.");
104 std::string detNamesConcat(
"");
109 boost::to_upper(detName);
111 detNamesConcat +=
"_" + detName;
117 auto iDetLayer = detName + std::to_string(iLayer);
121 distVarName +=
"__useHighestProbMassForExt";
151 std::unordered_map<int, const Particle*> targetParticles;
152 targetParticles.reserve(nMotherParticles);
154 for (
unsigned int iPart(0); iPart < nMotherParticles; ++iPart) {
164 auto iDetLayer = detName + std::to_string(iLayer);
173 targetParticles.insert({iDaughter->getMdstSource(), iDaughter});
179 targetParticles.insert({iParticle->getMdstSource(), iParticle});
186 const auto nParticlesTarget = targetParticles.size();
195 auto iDetLayer = detName + std::to_string(iLayer);
198 <<
"Detector surface: " << iDetLayer <<
"\n"
199 <<
"nMotherParticles: " << nMotherParticles <<
"\n"
200 <<
"nParticlesTarget: " << nParticlesTarget <<
"\n"
201 <<
"nParticlesReference: " << nParticlesReference);
203 double dummyDist(-1.0);
207 std::map<std::pair<int, int>,
double> particleMdstSourcePairsToDist;
210 for (
const auto& targetParticle : targetParticles) {
212 auto iMdstSource = targetParticle.first;
213 auto iParticle = targetParticle.second;
215 for (
unsigned int jPart(0); jPart < nParticlesReference; ++jPart) {
218 auto jMdstSource = jParticle->getMdstSource();
220 auto partMdstSourcePair = std::make_pair(iMdstSource, jMdstSource);
223 if (iMdstSource == jMdstSource) {
224 particleMdstSourcePairsToDist[partMdstSourcePair] = dummyDist;
234 if (particleMdstSourcePairsToDist.count({jMdstSource, iMdstSource})) {
235 particleMdstSourcePairsToDist[partMdstSourcePair] = particleMdstSourcePairsToDist[ {jMdstSource, iMdstSource}];
240 particleMdstSourcePairsToDist[partMdstSourcePair] = this->
getDistAtDetSurface(iParticle, jParticle, iDetLayer);
246 for (
const auto& targetParticle : targetParticles) {
248 auto iMdstSource = targetParticle.first;
249 auto iParticle = targetParticle.second;
252 std::vector<std::pair<double, int>> iDistancesAndRefMdstSources;
253 for (
const auto& [mdstIdxs, dist] : particleMdstSourcePairsToDist) {
254 if (mdstIdxs.first == iMdstSource) {
255 if (!std::isnan(dist) && dist >= 0) {
256 iDistancesAndRefMdstSources.push_back(std::make_pair(dist, mdstIdxs.second));
261 if (!iDistancesAndRefMdstSources.size()) {
262 B2DEBUG(12,
"The container of distances is empty. Perhaps the target and reference lists contain the same exact particles?");
266 const auto minDist = *std::min_element(iDistancesAndRefMdstSources.begin(), iDistancesAndRefMdstSources.end(),
267 [](
const auto & l,
const auto & r) {return l.first < r.first;});
269 auto jParticle =
m_pListReference->getParticleWithMdstSource(minDist.second);
272 <<
"Particle w/ mdstSource[" << iMdstSource <<
"] (PDG = "
273 << iParticle->getPDGCode() <<
"). Closest charged particle w/ mdstSource["
275 <<
"] (PDG = " << jParticle->getPDGCode()
276 <<
") at " << iDetLayer
277 <<
" surface is found at D = " << minDist.first
279 <<
"Storing extraInfo variables:\n"
296 for (
const auto& targetParticle : targetParticles) {
298 auto iMdstSource = targetParticle.first;
299 auto iParticle = targetParticle.second;
302 double isoScore(0.0);
303 double isoScoreAsWeightedAvg(0.0);
304 double sumWeights(0.0);
306 B2DEBUG(11,
"Particle w/ mdstSource[" << iMdstSource <<
"] (PDG = " << iParticle->getPDGCode() <<
").");
315 isoScore += weightedSumNonIsoLayers;
316 isoScoreAsWeightedAvg += weightedSumInvDists;
317 sumWeights += detWeight;
324 isoScore = (isoScore - minScore) / (maxScore - minScore);
326 isoScoreAsWeightedAvg /= sumWeights;
330 isoScoreAsWeightedAvg = 1. - ((std::min(isoScoreAsWeightedAvg, 1e2) - minScore) / (1e2 - minScore));
333 <<
"Isolation score: " << isoScore
335 <<
"Isolation score (as weighted avg): " << isoScoreAsWeightedAvg
337 <<
"Storing extraInfo variable:\n"
360 auto p = iParticle->
getP();
364 auto detWeight = (*
m_DBWeights.get())->getWeight(hypo, det, p, theta);
368 detWeight = (detWeight < 0 || std::isnan(detWeight)) ? detWeight : 0.0;
376 const float detWeight)
const
384 auto iDetLayer = detName + std::to_string(iLayer);
397 B2DEBUG(12,
"\nNo close-enough neighbours to this particle in the " << detName <<
" were found.");
401 B2FATAL(
"\nTotal layers in " << detName <<
": N=" << nLayers
403 <<
"Layers w/ a close-enough particle: n=" << n
405 <<
"n > N ?? Abort...");
408 return 1. - (-detWeight * (n / nLayers));
414 const float detWeight)
const
419 double sumInvDists(0.);
422 auto iDetLayer = detName + std::to_string(iLayer);
427 sumInvDists += threshDist / iParticle->
getExtraInfo(distVar);
432 return detWeight * sumInvDists;
439 const std::string& detLayerName)
const
452 std::string nameExtTheta =
"helixExtTheta(" + std::to_string(rho) +
"," + std::to_string(zfwd) +
"," + std::to_string(zbwd) +
")";
454 nameExtTheta.replace(nameExtTheta.size() - 1, 1,
", 1)");
459 std::string nameExtPhi =
"helixExtPhi(" + std::to_string(rho) +
"," + std::to_string(zfwd) +
"," + std::to_string(zbwd) +
")";
461 nameExtPhi.replace(nameExtPhi.size() - 1, 1,
", 1)");
466 const auto iExtInBarrel = (iExtTheta >= th_fwd_brl && iExtTheta < th_bwd_brl);
467 const auto jExtInBarrel = (jExtTheta >= th_fwd_brl && jExtTheta < th_bwd_brl);
469 const auto iExtInFWDEndcap = (iExtTheta >= th_fwd && iExtTheta < th_fwd_brl);
470 const auto jExtInFWDEndcap = (jExtTheta >= th_fwd && jExtTheta < th_fwd_brl);
472 const auto iExtInBWDEndcap = (iExtTheta >= th_bwd_brl && iExtTheta < th_bwd);
473 const auto jExtInBWDEndcap = (jExtTheta >= th_bwd_brl && jExtTheta < th_bwd);
475 if (boost::contains(detLayerName,
"CDC") || boost::contains(detLayerName,
"TOP")) {
478 if (!iExtInBarrel || !jExtInBarrel) {
479 return std::numeric_limits<double>::quiet_NaN();
486 auto diffPhi = jExtPhi - iExtPhi;
487 if (std::abs(diffPhi) > M_PI) {
488 diffPhi = (diffPhi > M_PI) ? diffPhi - 2 * M_PI : 2 * M_PI + diffPhi;
491 return 2.0 * rho * sin(std::abs(diffPhi) / 2.0);
495 if (boost::contains(detLayerName,
"ARICH")) {
497 if (!iExtInFWDEndcap || !jExtInFWDEndcap) {
498 return std::numeric_limits<double>::quiet_NaN();
501 if (boost::contains(detLayerName,
"ECL") || boost::contains(detLayerName,
"KLM")) {
506 !(iExtInBarrel && jExtInBarrel) &&
507 !(iExtInFWDEndcap && jExtInFWDEndcap) &&
508 !(iExtInBWDEndcap && jExtInBWDEndcap)
510 return std::numeric_limits<double>::quiet_NaN();
524 const auto iExtR = (iExtTheta >= th_fwd_brl && iExtTheta < th_bwd_brl) ? rho / sin(iExtTheta) : ((iExtTheta >= th_fwd
525 && iExtTheta < th_fwd_brl) ? zfwd / cos(iExtTheta) : zbwd / cos(iExtTheta));
526 const auto jExtR = (jExtTheta >= th_fwd_brl && jExtTheta < th_bwd_brl) ? rho / sin(jExtTheta) : ((jExtTheta >= th_fwd
527 && jExtTheta < th_fwd_brl) ? zfwd / cos(jExtTheta) : zbwd / cos(jExtTheta));
529 return sqrt((iExtR * iExtR) + (jExtR * jExtR) - 2 * iExtR * jExtR * (sin(iExtTheta) * sin(jExtTheta) * cos(iExtPhi - jExtPhi)
530 + cos(iExtTheta) * cos(jExtTheta)));
534 return std::numeric_limits<double>::quiet_NaN();
541 bool checkPList =
false;
555 bool checkPListReference =
false;
560 return (checkPList and checkPListReference);
Provides a type-safe way to pass members of the chargedStableSet set.
static const ParticleSet chargedStableSet
set of charged stable particles
static const ParticleType invalidParticle
Invalid particle, used internally.
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
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Class to store reconstructed particles.
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
int getPDGCode(void) const
Returns PDG code.
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
double getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
double getExtraInfo(const std::string &name) const
Return given value if set.
StoreObjPtr< EventMetaData > m_event_metadata
The event information.
std::unique_ptr< DBObjPtr< PIDDetectorWeights > > m_DBWeights
Interface to get the database payload with the PID detector weights.
~TrackIsoCalculatorModule() override
Destructor, use this to clean up anything you created in the constructor.
void initialize() override
Use this to initialize resources or memory your module needs.
double getDistAtDetSurface(const Particle *iParticle, const Particle *jParticle, const std::string &detLayerName) const
Calculate the distance between the points where the two input extrapolated track helices cross the gi...
void event() override
Called once for each event.
std::unordered_map< std::string, std::string > m_detLayerToDistVariable
Map that associates to each detector layer (e.g., 'CDC6') the name of the variable representing the d...
std::string m_decayString
The name of the input charged stable particle list, or composite particle w/ charged stable daughters...
StoreObjPtr< ParticleList > m_pListTarget
The input ParticleList object for which distances are to be calculated.
std::string m_isoScoreVariable
The name of the variable representing the track isolation score.
double getWeightedSumInvDists(const Particle *iParticle, const std::string &detName, const float detWeight) const
Get the sum of the inverse (scaled) minimum distances over the given detector layers,...
StoreArray< Particle > m_particles
StoreArray of Particles.
std::string m_payloadName
The name of the database payload object with the MVA weights.
std::vector< std::string > m_detNames
The list of names of the detectors at whose inner (cylindrical) surface we extrapolate each track's p...
Const::EDetector getDetEnum(const std::string &detName) const
Get the enum type for this detector name.
std::string m_pListReferenceName
The name of the input ParticleList of reference tracks.
std::unordered_map< std::string, std::string > m_detLayerToRefPartIdxVariable
Map that associates to each detector layer (e.g, 'CDC6') the name of the variable representing the md...
double getDistThreshold(Const::EDetector det, int layer) const
Get the threshold value per detector layer for the distance to closest ext.
bool onlySelectedStdChargedInDecay()
Check whether input particle list and reference list are of a valid charged stable particle.
double getDetectorWeight(const Particle *iParticle, const std::string &detName) const
Get the PID weight, , for this particle and detector reading it from the payload, if selected.
bool m_excludePIDDetWeights
Exclude the PID detector weights for the isolation score definition.
DecayDescriptor m_decaydescriptor
< Decay descriptor of decays to look for.
TrackIsoCalculatorModule()
Constructor: Sets the description, the properties and the parameters of the module.
double getWeightedSumNonIsoLayers(const Particle *iParticle, const std::string &detName, const float detWeight) const
Get the sum of layers with a close-by track, divided by the total number of layers,...
unsigned short m_nSelectedDaughters
The number of selected daughters in the decay string.
StoreObjPtr< ParticleList > m_pListReference
The input ParticleList object of reference tracks.
std::string m_isoScoreVariableAsWeightedAvg
The name of the variable representing the track isolation score.
bool m_useHighestProbMassForExt
If this option is set, the helix extrapolation for the target and reference particles will use the tr...
static Manager & Instance()
get singleton instance.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
Abstract base class for different kinds of events.
static const std::unordered_map< std::string, DetSurfCylBoundaries > detLayerToSurfBoundaries
Map that associates to each detector layer its valid cylindrical surface's boundaries.
static const std::unordered_map< std::string, std::vector< int > > detToLayers
Map that associates to each detector its list of valid layers.