Belle II Software development
GRLNeuroModule.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 <framework/datastore/StoreArray.h>
11#include <framework/datastore/StoreObjPtr.h>
12#include <framework/dataobjects/EventMetaData.h>
13#include <framework/core/ModuleParam.templateDetails.h>
14#include <framework/gearbox/Unit.h>
15#include <framework/utilities/FileSystem.h>
16#include <analysis/utility/PCmsLabTransform.h>
17#include <trg/ecl/TrgEclMapping.h>
18#include <trg/ecl/dataobjects/TRGECLCluster.h>
19#include <trg/grl/dataobjects/GRLMLPData.h>
20#include <trg/grl/dataobjects/TRGGRLInfo.h>
21#include "trg/grl/dataobjects/TRGGRLUnpackerStore.h"
22#include "trg/grl/modules/trggrlneuralnet/GRLNeuroModule.h"
23
24#include <fstream>
25#include <cmath>
26#include <TFile.h>
27
28using namespace Belle2;
29using namespace std;
30
31//-----------------------------------------------------------------
32// Register the Module
33//-----------------------------------------------------------------
35
36//-----------------------------------------------------------------
37// Implementation
38//-----------------------------------------------------------------
39//
41{
43 "The NeuroTrigger module of the GRL.\n"
44 "Takes CDC track and ECL cluster to prepare input data\n"
45 "for the training of a neural network.\n"
46 "Networks are trained after the event loop and saved."
47 );
49 // parameters for saving / loading
50 addParam("TrgGrlInformation", m_TrgGrlInformationName,
51 "Name of the StoreArray holding the information of tracks and clusters from cdc ecl klm.",
52 string("TRGGRLObjects"));
53 addParam("HistFileName", m_HistFileName,
54 "Name of the root file saving the output histogram.",
55 string("hist.root"));
56 addParam("SaveHist", m_saveHist,
57 "Save the output histogram to root file.",
58 false);
59 addParam("nMLP", m_parameters.nMLP,
60 "Number of expert MLPs.", 1u);
61 addParam("n_cdc_sector", m_parameters.n_cdc_sector,
62 "Number of expert CDC MLPs.", 0u);
63 addParam("n_ecl_sector", m_parameters.n_ecl_sector,
64 "Number of expert ECL MLPs.", 1u);
65 addParam("i_cdc_sector", m_parameters.i_cdc_sector,
66 "#cdc track of expert MLPs.", {0});
67 addParam("i_ecl_sector", m_parameters.i_ecl_sector,
68 "#ecl cluster of expert MLPs.", {24});
69 addParam("nHidden", m_parameters.nHidden,
70 "Number of nodes in each hidden layer for all networks "
71 "or factor to multiply with number of inputs (1 list or nMLP lists). "
72 "The number of layers is derived from the shape.",
73 {{64, 64}});
74 addParam("multiplyHidden", m_parameters.multiplyHidden,
75 "If true, multiply nHidden with number of input nodes.", false);
76 addParam("outputScale", m_parameters.outputScale,
77 "Output scale for all networks (1 value list or nMLP value lists). "
78 "Output[i] of the MLP is scaled from [-1, 1] "
79 "to [outputScale[2*i], outputScale[2*i+1]]. "
80 "(units: z[cm] / theta[degree])",
81 {{ -1., 1.}});
82 addParam("weightFiles", m_weightFileNames,
83 "Name of the file where the weights of MLPs are saved. "
84 "the default file is $BELLE2_LOCAL_DIR/data/trg/grl/weights.dat",
85 {FileSystem::findFile("/data/trg/grl/weights.dat", true)});
86 addParam("biasFiles", m_biasFileNames,
87 "Name of the file where the biases of MLPs are saved. "
88 "the default file is $BELLE2_LOCAL_DIR/data/trg/grl/bias.dat",
89 {FileSystem::findFile("/data/trg/grl/bias.dat", true)});
90 addParam("TRGECLClusters", m_TrgECLClusterName,
91 "Name of the StoreArray holding the information of trigger ecl clusters ",
92 string("TRGECLClusters"));
93 addParam("MVACut", m_nn_thres,
94 "Cut value applied to the MLP output",
95 {-1});
96 addParam("useDB", m_useDB,
97 "Flag to use database to set config", true);
98}
99
100
101void
103{
104
105 TrgEclMapping* trgecl_obj = new TrgEclMapping();
106 for (int tc = 1; tc <= 576; tc++) {
107 TCThetaID.push_back(trgecl_obj->getTCThetaIdFromTCId(tc));
108 }
109
110 //-----------------------------------------------------------------------------------------
111 //..ECL look up tables
112 PCmsLabTransform boostrotate;
113 for (int tc = 1; tc <= 576; tc++) {
114
115 //..Four vector of a 1 GeV lab photon at this TC
116 ROOT::Math::XYZVector CellPosition = trgecl_obj->getTCPosition(tc);
117 ROOT::Math::PxPyPzEVector CellLab;
118 CellLab.SetPx(CellPosition.Unit().X());
119 CellLab.SetPy(CellPosition.Unit().Y());
120 CellLab.SetPz(CellPosition.Unit().Z());
121 CellLab.SetE(1.);
122
123 TCThetaLab.push_back(CellLab.Theta()*TMath::RadToDeg());
124 TCPhiLab.push_back(CellLab.Phi()*TMath::RadToDeg() + 180.0);
125
126 }
127
128 delete trgecl_obj;
129
130}
131
132void
134{
135 if (m_useDB) {
136 if (not m_db_trggrlconfig.isValid()) {
137 StoreObjPtr<EventMetaData> evtMetaData;
138 B2FATAL("No database for TRG GRL config. exp " << evtMetaData->getExperiment() << " run "
139 << evtMetaData->getRun());
140 } else {
141 m_parameters.nMLP = m_db_trggrlconfig->get_ecltaunn_nMLP();
142 m_parameters.multiplyHidden = m_db_trggrlconfig->get_ecltaunn_multiplyHidden();
143 m_parameters.nHidden = m_db_trggrlconfig->get_ecltaunn_nHidden();
144 m_parameters.n_cdc_sector = m_db_trggrlconfig->get_ecltaunn_n_cdc_sector();
145 m_parameters.n_ecl_sector = m_db_trggrlconfig->get_ecltaunn_n_ecl_sector();
146 m_parameters.i_cdc_sector = m_db_trggrlconfig->get_ecltaunn_i_cdc_sector();
147 m_parameters.i_ecl_sector = m_db_trggrlconfig->get_ecltaunn_i_ecl_sector();
148 m_nn_thres[0] = m_db_trggrlconfig->get_ecltaunn_threshold();
149
150 m_GRLNeuro.initialize(m_parameters);
151
152 for (unsigned int isector = 0; isector < m_parameters.nMLP; isector++) {
153 if (!m_GRLNeuro.load(isector, m_db_trggrlconfig->get_ecltaunn_weight()[isector], m_db_trggrlconfig->get_ecltaunn_bias()[isector])) {
154 B2ERROR("weight of GRL ecltaunn could not be loaded correctly.");
155 }
156 }
157 }
158 }
159
160 if (m_saveHist) {
161 for (unsigned int isector = 0; isector < m_parameters.nMLP; isector++) {
162 h_target.push_back(new TH1D(("h_target_" + to_string(isector)).c_str(),
163 ("h_target_" + to_string(isector)).c_str(), 100, 0.0, 1.0));
164 }
165 }
166}
167
168void
170{
171
172 //ECL input
173 //..Use only clusters within 100 ns of event timing (from ECL).
174 //StoreArray<TRGECLTrg> trgArray;
178 std::vector<std::tuple<float, float, float, float>> eclClusters;
180 int necl = eclTrgClusterArray.getEntries();
181 for (int ic = 0; ic < necl; ic++) {
182 int TC = eclTrgClusterArray[ic]->getMaxTCId();
183 float energy = eclTrgClusterArray[ic]->getEnergyDep() * 1000.0;
184 float theta = TCThetaLab[TC - 1];
185 float phi = TCPhiLab[TC - 1];
186 float time = eclTrgClusterArray[ic]->getTimeAve();
187 eclClusters.emplace_back(energy, theta, phi, time);
188 }
189
190 std::sort(eclClusters.begin(), eclClusters.end(),
191 [](const std::tuple<float, float, float, float>& a,
192 const std::tuple<float, float, float, float>& b) {
193 return std::get<0>(a) > std::get<0>(b);
194 });
195
196 // MLPinput
197 std::vector<float> MLPinput;
198 MLPinput.clear();
199 MLPinput.assign(24, 0);
200
201 for (size_t i = 0; i < std::min(eclClusters.size(), size_t(6)); i++) {
202 MLPinput[i] = std::get<0>(eclClusters[i]); // E
203 MLPinput[6 + i] = std::get<1>(eclClusters[i]); // Theta
204 MLPinput[12 + i] = std::get<2>(eclClusters[i]); // Phi
205 MLPinput[18 + i] = std::get<3>(eclClusters[i]); // time
206 }
207
208
209
210 // //////////////////////////////////////////////////////////////////////////////////////////////////////////
211
212 //Energy, theta and phi data are the quantized.
213 //Energy: LSB = ADC count (5.5MeV)
214 //theta : 0 ~ 180∘, LSB = 1.6025∘
215 //phi : 0 ~ 360∘, LSB = 1.6025∘
216 float LSB_ADC = 1 / 5.5;
217 float LSB_angle = 1 / 1.6025;
219 std::for_each(MLPinput.begin() + 0, MLPinput.begin() + 6, [LSB_ADC](float & x) { x = std::ceil(x * LSB_ADC); });
220 std::for_each(MLPinput.begin() + 6, MLPinput.begin() + 12, [LSB_angle](float & x) { x = std::ceil(x * LSB_angle);});
221 std::for_each(MLPinput.begin() + 12, MLPinput.begin() + 18, [LSB_angle](float & x) { x = std::ceil(x * LSB_angle); });
222
223//new add
224// float LSB_time = 1.0f;
225// std::for_each(MLPinput.begin() + 18, MLPinput.begin() + 24, [](float &x) {
226// x = std::ceil(x);
227// x = std::min(x, 255.0f);
228// x = std::max(x, 0.0f);
229// });
230
231
232 float target = m_GRLNeuro.runMLP(0, MLPinput);
233 if (m_saveHist) {
234 h_target[0]->Fill(target);
235 }
236 if (target > m_nn_thres[0]) {
237 trgInfo->setTauNN(true);
238 } else trgInfo->setTauNN(false);
239
240}
241
242void
244{
245 if (m_saveHist) {
246 TFile* file = new TFile(TString(m_HistFileName), "recreate");
247 for (unsigned int isector = 0; isector < m_parameters.nMLP; isector++) {
248 h_target[isector]->Write();
249 }
250 file->Write();
251 file->Close();
252 }
253
254}
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
std::vector< double > TCThetaLab
Polar angle of a given TRGcluster.
DBObjPtr< TRGGRLConfig > m_db_trggrlconfig
dbobject to store grl config
std::vector< double > TCPhiLab
Azimuthal angle of a given TRGcluster.
std::string m_HistFileName
Name of root file to save the histogram.
virtual void initialize() override
Initialize the module.
std::string m_TrgECLClusterName
Name of the StoreArray containing the ECL clusters.
std::vector< int > TCThetaID
TCID of a given TRGcluster.
GRLNeuro m_GRLNeuro
Instance of the NeuroTrigger.
virtual void event() override
Called once for each event.
GRLNeuro::Parameters m_parameters
Parameters for the NeuroTrigger.
GRLNeuroModule()
Constructor, for setting module description and parameters.
virtual void terminate() override
This method is called at the end of the event processing.
virtual void beginRun() override
Register run-dependent DataStore arrays.
std::string m_TrgGrlInformationName
name of TRG GRL information object
bool m_saveHist
save the output histogram
std::vector< float > m_nn_thres
cut on MVA to separate signal
std::vector< TH1D * > h_target
Histograms to save the NN classifiers.
bool m_useDB
flag to use database to load config
Class to represent the GRL Neuro.
Definition GRLNeuro.h:35
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 hold Lorentz transformations from/to CMS and boost vector.
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
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
A class of TC Mapping.
int getTCThetaIdFromTCId(int)
get [TC Theta ID] from [TC ID]
ROOT::Math::XYZVector getTCPosition(int)
TC position (cm)
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:559
#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.
STL namespace.