Belle II Software light-2406-ragdoll
NeutralHadronMatcherModule.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/NeutralHadronMatcher/NeutralHadronMatcherModule.h>
10#include <framework/dataobjects/EventExtraInfo.h>
11#include <framework/core/Environment.h>
12#include <framework/datastore/StoreArray.h>
13#include <framework/datastore/StoreObjPtr.h>
14#include <TRandom.h>
15#include <analysis/variables/MCTruthVariables.h>
16
17using namespace Belle2;
18
19//-----------------------------------------------------------------
20// Register the Module
21//-----------------------------------------------------------------
22REG_MODULE(NeutralHadronMatcher);
23
24NeutralHadronMatcherModule::NeutralHadronMatcherModule() : Module()
25{
26 setDescription("Perform geometrical match between MC neutral hadron (given by mcPDG) and ECL clusters from the particleLists");
27 addParam("particleLists", m_ParticleLists, "Input particle list");
28 addParam("efficiencyCorrection", m_effcorr, "data/mc efficiency ratio", 0.83);
29 addParam("distanceCut", m_distance, "Matching distance", 15.0);
30 addParam("mcPDGcode", m_mcPDG, "MC PDG code of the neutral hadron", Const::Klong.getPDGCode());
31 addParam("ignoreClustersWithPDGcodes", m_PDGignore, "Do not attempt to match clusters that are already matched with specific codes",
32 std::vector<int> {Const::photon.getPDGCode()});
33}
34
36{
37 if (m_mcPDG == Const::Klong.getPDGCode()) {
38 m_infoName = "mcdistanceKL";
39 m_matchedId = "mdstIndexTruthKL";
40 } else if (abs(m_mcPDG) == Const::neutron.getPDGCode()) {
41 m_infoName = "mcdistanceNeutron";
42 m_matchedId = "mdstIndexTruthNeutron";
43 } else {
44 B2FATAL("Unsupported mcPDG value: " << m_mcPDG);
45 }
46
47 for (auto& iList : m_ParticleLists)
49
50 if (! StoreArray<MCParticle>().isValid())
51 B2WARNING("No MCParticles array found. Do you call the NeutralHadronMatcher with real data?");
52
53}
54
56{
57
58 StoreArray<MCParticle> mcparticles;
59 if (!mcparticles.isValid())
60 return;
61 if (mcparticles.getEntries() < 1)
62 B2FATAL("Missing mcparticles list for MC");
63
64 // Initialize extra info
65 for (auto& iList : m_ParticleLists) {
66 StoreObjPtr<ParticleList> particleList(iList);
67
68 //check particle List has particles
69 const size_t nPart = particleList->getListSize();
70 for (size_t iPart = 0; iPart < nPart; iPart++) {
71 auto particle = particleList->getParticle(iPart);
72 particle->addExtraInfo(m_infoName, 1.e10);
73 particle->addExtraInfo(m_matchedId, -1);
74 }
75 }
76
77 // collect only primary mcparticles that match m_mcPDG
78 std::vector< std::pair<MCParticle*, bool> > primaryMCParticles;
79 for (int i = 0; i < mcparticles.getEntries(); i++) {
80 auto mcPart = mcparticles[i];
81 if (mcPart->isPrimaryParticle() && abs(mcPart->getPDG()) == m_mcPDG) {
82 bool AddInefficiency = (gRandom->Uniform() > m_effcorr);
83 primaryMCParticles.push_back(std::pair(mcPart, AddInefficiency));
84 }
85 }
86
87 const size_t nMCPart = primaryMCParticles.size();
88
89 // if no primary MCParticles matching m_mcPDG exists, nothing to do.
90 if (nMCPart == 0)
91 return;
92
93 for (auto& iList : m_ParticleLists) {
94 StoreObjPtr<ParticleList> particleList(iList);
95
96 // loop over particles first
97 const size_t nPart = particleList->getListSize();
98 for (size_t iPart = 0; iPart < nPart; iPart++) {
99
100 auto particle = particleList->getParticle(iPart);
101 auto bestMatch = Variable::particleClusterBestMCPDGCode(particle);
102
103 if ((bestMatch) && (std::find(m_PDGignore.begin(), m_PDGignore.end(), int(bestMatch)) != m_PDGignore.end())) {
104 continue;
105 }
106
107 const ECLCluster* eclcluster = particle->getECLCluster();
108 if (!eclcluster) continue;
109
110 // create a vector of distances from MCParticles
111 std::vector< std::pair <double, bool> > distances;
112 distances.resize(nMCPart);
113
114 // loop over mcparticles
115 for (size_t iMCPart = 0; iMCPart < nMCPart; iMCPart++) {
116
117 auto mcPart = primaryMCParticles[iMCPart].first;
118
119 auto vtx = mcPart->getProductionVertex();
120 auto momentum = mcPart->getMomentum();
121
122 // got potential KL candidate
123 double R = eclcluster->getR();
124 double phi = eclcluster->getPhi();
125 double theta = eclcluster->getTheta();
126 double zcl = R / tan(theta);
127 double xcl = R * cos(phi);
128 double ycl = R * sin(phi);
129
130 double rECL = R; //cm
131 double z = vtx.Z() + rECL / tan(momentum.Theta());
132 double x = vtx.X() + rECL * cos(momentum.Phi());
133 double y = vtx.Y() + rECL * sin(momentum.Phi());
134
135 double dist = sqrt((x - xcl) * (x - xcl) + (y - ycl) * (y - ycl) + (z - zcl) * (z - zcl));
136 distances[iMCPart] = std::pair(dist, primaryMCParticles[iMCPart].second);
137 }
138
139 auto it_distMin = std::min_element(distances.begin(), distances.end(),
140 [](const auto & a, const auto & b) { return a.first < b.first;});
141 double distMin = (*it_distMin).first;
142 bool AddInefficiency = (*it_distMin).second;
143 if ((distMin < m_distance) && AddInefficiency)
144 distMin = -distMin;
145
146 particle->setExtraInfo(m_infoName, distMin);
147
148 if ((distMin < m_distance) && (distMin > 0))
149 particle->setExtraInfo(m_matchedId, primaryMCParticles[std::distance(distances.begin(), it_distMin)].first->getArrayIndex());
150
151 }
152 }
153
154}
155
156
int getPDGCode() const
PDG code.
Definition: Const.h:473
static const ParticleType neutron
neutron particle
Definition: Const.h:675
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:678
static const ParticleType photon
photon particle
Definition: Const.h:673
ECL cluster data.
Definition: ECLCluster.h:27
double getPhi() const
Return Corrected Phi of Shower (radian).
Definition: ECLCluster.h:304
double getR() const
Return R.
Definition: ECLCluster.h:310
double getTheta() const
Return Corrected Theta of Shower (radian).
Definition: ECLCluster.h:307
Base class for Modules.
Definition: Module.h:72
std::vector< std::string > m_ParticleLists
input particle lists name
virtual void initialize() override final
Overridden initialize method.
double m_effcorr
input efficiency correction
std::string m_infoName
extra info variable for distance
virtual void event() override final
Overridden event method.
std::string m_matchedId
extra info variable for matched MC
std::vector< int > m_PDGignore
ignore clusters that are matched with the following PDG codes
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
bool isValid() const
Check wether the array was registered.
Definition: StoreArray.h:288
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
#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