Belle II Software development
KLMMuonIDDNNExpertModule.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
10#include <reconstruction/modules/KLMMuonIDDNNExpert/KLMMuonIDDNNExpertModule.h>
11
12#include <klm/bklm/geometry/GeometryPar.h>
13#include <klm/bklm/geometry/Module.h>
14#include <klm/eklm/geometry/TransformDataGlobalAligned.h>
15#include <klm/dataobjects/KLMHit2d.h>
16#include <klm/dataobjects/KLMMuidLikelihood.h>
17
18#include <reconstruction/dataobjects/KLMMuonIDDNNInputVariable.h>
19
20#include <CLHEP/Units/SystemOfUnits.h>
21
22#include <framework/logging/Logger.h>
23#include <framework/gearbox/Const.h>
24
25#include <tracking/dataobjects/ExtHit.h>
26
27#include <mdst/dataobjects/PIDLikelihood.h>
28
29#include <mdst/dataobjects/Track.h>
30
31#include <mva/interface/Interface.h>
32#include <mva/dataobjects/DatabaseRepresentationOfWeightfile.h>
33#include <mva/interface/Expert.h>
34#include <mva/interface/Weightfile.h>
35
36#include <iostream>
37#include <fstream>
38#include <sstream>
39
40using namespace Belle2;
41
42REG_MODULE(KLMMuonIDDNNExpert);
43
45{
46 // Set module properties
47 setDescription(R"DOC(Get information from KLMMuIDLikelihood)DOC");
49}
50
54
56{
57 m_tracks.isRequired();
58
59 m_inputVariable.registerInDataStore();
60 m_tracks.registerRelationTo(m_inputVariable);
61
62 // setup KLM geometry
65
66 m_EndcapScintWidth = eklmGeometry.getStripGeometry()->getWidth() / CLHEP::cm; // in G4e units (cm)
67
68 for (int layer = 1; layer <= m_maxBKLMLayers; ++layer) {
69 const bklm::Module* module =
70 bklmGeometry->findModule(BKLMElementNumbers::c_ForwardSection, 1, layer);
71 m_BarrelPhiStripWidth[layer - 1] = module->getPhiStripWidth(); // in G4e units (cm)
72 m_BarrelZStripWidth[layer - 1] = module->getZStripWidth(); // in G4e units (cm)
73 }
74
76}
77
79{
80 m_expert.reset();
81 m_dataset.reset();
82}
83
85{
86 if (m_weightfile_representation.isValid()) {
87 if (m_weightfile_representation.hasChanged()) {
88 std::stringstream ss(m_weightfile_representation->m_data);
89 auto weightfile = MVA::Weightfile::loadFromStream(ss);
90 initializeMVA(weightfile);
91 }
92 } else {
93 B2FATAL("Payload " << m_identifier << " is not valid!");
94 }
95}
96
100
102{
103 auto supported_interfaces = MVA::AbstractInterface::getSupportedInterfaces();
104 MVA::GeneralOptions general_options;
105 weightfile.getOptions(general_options);
106 m_expert = supported_interfaces[general_options.m_method]->getExpert();
107 m_expert->load(weightfile);
108 std::vector<float> dummy;
109 int nInputVariables = general_options.m_variables.size();
110 if (nInputVariables != 5 + 4 * (m_maxBKLMLayers + m_maxEKLMLayers)) {
111 B2FATAL("Number of input variables mismatch. Required " << 5 + 4 * (m_maxBKLMLayers + m_maxEKLMLayers) << " but " << nInputVariables
112 << " given. ");
113 }
114 dummy.resize(nInputVariables, 0);
115 m_dataset = std::unique_ptr<MVA::SingleDataset>(new MVA::SingleDataset(general_options, std::move(dummy), 0));
116}
117
118
120{
121 for (Track& track : m_tracks) {
122
123 KLMMuidLikelihood* klmll = track.getRelatedTo<KLMMuidLikelihood>();
124
125 if (!klmll) continue;
126
127 // initialize hit pattern arrays
128 m_hitpattern_steplength.fill(-1.);
129 m_hitpattern_width.fill(-1.);
130 m_hitpattern_chi2.fill(-1.);
131 m_hitpattern_hasext.fill(0);
132
133 bool hasExtInKLM = false;
134 for (const ExtHit& exthit : track.getRelationsTo<ExtHit>()) {
135
136 if (exthit.getDetectorID() != Const::EDetector::BKLM
137 and exthit.getDetectorID() != Const::EDetector::EKLM) continue;
138
139 int layer;
140 bool inBKLM = (exthit.getDetectorID() == Const::EDetector::BKLM);
141 int copyid = exthit.getCopyID();
142
143 int section, sector, plane, strip;
144
145 if (inBKLM) {
146 BKLMElementNumbers::moduleNumberToElementNumbers(copyid, &section, &sector, &layer);
147 if (layer > m_maxBKLMLayers) continue;
148 m_hitpattern_hasext.at(layer - 1) = 1;
149 } else {
150 EKLMElementNumbers::Instance().stripNumberToElementNumbers(copyid, &section, &layer, &sector, &plane, &strip);
151 if (layer > m_maxEKLMLayers) continue;
152 m_hitpattern_hasext.at(m_maxBKLMLayers + layer - 1) = 1;
153 }
154
155 hasExtInKLM = true;
156
157 }
158
159 RelationVector<KLMHit2d> KLMHit2drelation = track.getRelationsTo<KLMHit2d>();
160
161 // only apply NN muonID to tracks with at least one KLMHit2d or one ExtHit in KLM.
162 if (not(hasExtInKLM || KLMHit2drelation.size())) continue;
163
164 std::map<int, int> Hit2dMap; // arrange KLMHit2d in the order of layer
165 for (long unsigned int ii = 0; ii < KLMHit2drelation.size(); ii++) {
166 KLMHit2d* klmhit = KLMHit2drelation[ii];
167 bool hit_inBKLM = (klmhit->getSubdetector() == KLMElementNumbers::c_BKLM);
168 unsigned long int hit_layer = klmhit->getLayer();
169
170 int index = hit_layer - 1 + m_maxBKLMLayers * (1 - hit_inBKLM); // BKLM hits are in front of EKLM hits
171 if (index > (m_maxBKLMLayers + m_maxEKLMLayers)) continue;
172 Hit2dMap.insert(std::pair<int, int> {index, ii});
173 }
174
175 int nklmhits = 0;
176 ROOT::Math::XYZVector previousPosition(0., 0., 0.);
177 for (auto itermap = Hit2dMap.begin(); itermap != Hit2dMap.end(); itermap ++) {
178
179 nklmhits += 1;
180
181 KLMHit2d* klmhit = KLMHit2drelation[itermap->second];
182
183 float KFchi2 = KLMHit2drelation.weight(itermap->second);
184 float width = getHitWidth(klmhit);
185
186 ROOT::Math::XYZVector hitPosition = klmhit->getPosition();
187 float steplength = 0.;
188 if (nklmhits > 1) {
189 steplength = (hitPosition - previousPosition).R();
190 }
191 previousPosition = hitPosition;
192
193 // hit pattern creation.
194 int hitpatternindex = itermap->first;
195 m_hitpattern_chi2.at(hitpatternindex) = KFchi2;
196 m_hitpattern_steplength.at(hitpatternindex) = steplength;
197 m_hitpattern_width.at(hitpatternindex) = width;
198 } // loop of Hit2dMap
199
200 Hit2dMap.clear();
201
202 double muprob_nn = getNNmuProbability(&track, klmll);
203 PIDLikelihood* pid = track.getRelated<PIDLikelihood>();
204 pid->addPreOfficialLikelihood("klmMuonIDDNN", muprob_nn);
205 } // loop of tracks
206}
207
209{
210 m_dataset->m_input[0] = klmll->getChiSquared();
211 m_dataset->m_input[1] = klmll->getDegreesOfFreedom();
212 m_dataset->m_input[2] = klmll->getExtLayer() - klmll->getHitLayer();
213 m_dataset->m_input[3] = klmll->getExtLayer();
214 m_dataset->m_input[4] = track->getTrackFitResultWithClosestMass(Const::muon)->getTransverseMomentum();
215
216 for (int layer = 0; layer < (m_maxBKLMLayers + m_maxEKLMLayers); layer ++) {
217 m_dataset->m_input[5 + 4 * layer + 0] = m_hitpattern_width.at(layer); // width
218 m_dataset->m_input[5 + 4 * layer + 1] = m_hitpattern_steplength.at(layer); // steplength
219 m_dataset->m_input[5 + 4 * layer + 2] = m_hitpattern_chi2.at(layer); // chi2
220 m_dataset->m_input[5 + 4 * layer + 3] = m_hitpattern_hasext.at(layer); // hasext
221 }
222
223 KLMMuonIDDNNInputVariable* inputVariable = m_inputVariable.appendNew();
224 inputVariable->setKLMMuonIDDNNInputVariable(m_dataset->m_input);
225 track->addRelationTo(inputVariable);
226
227 float muprob_nn = m_expert->apply(*m_dataset)[0];
228 return muprob_nn;
229}
230
232{
233 float stripwidth1 = 0; // strip width of phi or X direction
234 float stripwidth2 = 0; // strip width of Z or Y direction
235 float stripdiff1 = 0; // max minus min strip number in phi or X direction
236 float stripdiff2 = 0; // max minus min strip number in Z or Y direction
237 if (klmhit->getSubdetector() == KLMElementNumbers::c_BKLM) {
238 stripwidth1 = m_BarrelPhiStripWidth[klmhit->getLayer() - 1];
239 stripwidth2 = m_BarrelZStripWidth[klmhit->getLayer() - 1];
240 stripdiff1 = (klmhit->getPhiStripMax() - klmhit->getPhiStripMin() + 1) * 0.5;
241 stripdiff2 = (klmhit->getZStripMax() - klmhit->getZStripMin() + 1) * 0.5;
242 } else {
243 stripwidth1 = m_EndcapScintWidth;
244 stripwidth2 = m_EndcapScintWidth;
245 stripdiff1 = (klmhit->getXStripMax() - klmhit->getXStripMin() + 1) * 0.5;
246 stripdiff2 = (klmhit->getYStripMax() - klmhit->getYStripMin() + 1) * 0.5;
247 }
248 float width1 = stripwidth1 * stripdiff1;
249 float width2 = stripwidth2 * stripdiff2;
250 return std::sqrt(width1 * width1 + width2 * width2);
251}
252
253
double R
typedef autogenerated by FFTW
static void moduleNumberToElementNumbers(KLMModuleNumber module, int *section, int *sector, int *layer)
Get element numbers by module number.
static const ChargedStable muon
muon particle
Definition Const.h:660
static const EKLMElementNumbers & Instance()
Instantiation.
void stripNumberToElementNumbers(int stripGlobal, int *section, int *layer, int *sector, int *plane, int *strip) const
Get element numbers by strip global number.
double getWidth() const
Get width.
const StripGeometry * getStripGeometry() const
Get strip geometry data.
EKLM geometry data.
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
Store one Ext hit as a ROOT object.
Definition ExtHit.h:31
KLM 2d hit.
Definition KLMHit2d.h:33
Class to store the likelihoods from KLM with additional information related to the extrapolation.
int getDegreesOfFreedom() const
Get the number of degrees of freedom (= 2 times the number of KLM hits) for the chi-squared computati...
int getHitLayer() const
Get the outermost KLM layer actually crossed by the track.
int getExtLayer() const
Get the outermost KLM layer crossed in the extrapolation.
double getChiSquared() const
Get the chi-squared of the extrapolation.
std::array< float, m_maxBKLMLayers+m_maxEKLMLayers > m_hitpattern_width
Container of hit widths of one track.
std::unique_ptr< MVA::SingleDataset > m_dataset
Pointer to the current dataset.
float getHitWidth(const KLMHit2d *klmhit)
Get Hit width (cluster size) of a KLMHit2d.
void event() override
This method is called for each event.
float m_EndcapScintWidth
EKLM scintillator strip width (cm).
void endRun() override
This method is called if the current run ends.
static constexpr int m_maxBKLMLayers
Total BKLM layers.
std::array< float, m_maxBKLMLayers+m_maxEKLMLayers > m_hitpattern_chi2
Container of hit chi2 of one track.
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.
void beginRun() override
Called when entering a new run.
StoreArray< Track > m_tracks
Required array for Tracks.
std::array< float, m_maxBKLMLayers+m_maxEKLMLayers > m_hitpattern_steplength
Container of hit steplength of one track.
float getNNmuProbability(const Track *track, const KLMMuidLikelihood *klmll)
Get the NN-based muon probability.
std::array< float, m_maxBKLMLayers > m_BarrelPhiStripWidth
BKLM phi-measuring strip width (cm) by layer.
DBObjPtr< DatabaseRepresentationOfWeightfile > m_weightfile_representation
Database pointer to the database representation of the weightfile.
void initializeMVA(MVA::Weightfile &weightfile)
Initialize mva expert, dataset and features.
std::array< float, m_maxBKLMLayers > m_BarrelZStripWidth
BKLM Z-measuring strip width (cm) by layer.
std::array< bool, m_maxBKLMLayers+m_maxEKLMLayers > m_hitpattern_hasext
Container of extrapolation situation at each KLM layer of one track.
const std::string m_identifier
Database identifier or file used to load the weights.
StoreArray< KLMMuonIDDNNInputVariable > m_inputVariable
Input variables of DNN.
static constexpr int m_maxEKLMLayers
Total EKLM layers.
KLM MuonID DNN input variables datastore object to store the input variables for retraining KLMMuonID...
void setKLMMuonIDDNNInputVariable(const std::vector< float > &inputdata)
set the DNN input variables.
static void initSupportedInterfaces()
Static function which initliazes all supported interfaces, has to be called once before getSupportedI...
Definition Interface.cc:46
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.
void getOptions(Options &options) const
Fills an Option object from the xml tree.
Definition Weightfile.cc:67
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
Module()
Constructor.
Definition Module.cc:30
@ 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 collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
float weight(int index) const
Get weight with index.
Class that bundles various TrackFitResults.
Definition Track.h:25
Provides BKLM geometry parameters for simulation, reconstruction etc (from Gearbox or DataBase)
Definition GeometryPar.h:37
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
static GeometryPar * instance(void)
Static method to get a reference to the singleton GeometryPar instance.
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition Module.h:76
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.