Belle II Software  release-05-01-25
ECLClusterPSD.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * This module computes shower variables using pulse shape *
6  * information from offline two component fits. Using pulse *
7  * shape discrimination, these shower variables can be used *
8  * for particle id. *
9  * *
10  * Author: The Belle II Collaboration *
11  * Contributors: Savino Longo (longos@uvic.ca) *
12  * *
13  * This software is provided "as is" without any warranty. *
14  **************************************************************************/
15 //This module
16 #include <ecl/modules/eclClusterPSD/ECLClusterPSD.h>
17 
18 //BOOST
19 #include <boost/algorithm/string/predicate.hpp>
20 
21 // ECL
22 #include <ecl/dataobjects/ECLShower.h>
23 #include <ecl/dataobjects/ECLCalDigit.h>
24 #include <ecl/geometry/ECLGeometryPar.h>
25 
26 //MVA
27 #include <mva/interface/Expert.h>
28 #include <mva/interface/Weightfile.h>
29 #include <mva/interface/Interface.h>
30 #include <mva/dataobjects/DatabaseRepresentationOfWeightfile.h>
31 
32 #include<math.h>
33 
34 // FRAMEWORK
35 #include <framework/logging/Logger.h>
36 #include <framework/geometry/B2Vector3.h>
37 
38 using namespace Belle2;
39 
40 //-----------------------------------------------------------------
41 // Register the Modules
42 //-----------------------------------------------------------------
43 REG_MODULE(ECLClusterPSD)
44 //-----------------------------------------------------------------
45 // Implementation
46 //-----------------------------------------------------------------
47 
48 // constructor
50 {
51  // Set module properties
52  setDescription("Module uses offline two component fit results to compute pulse shape discrimation variables for particle identification.");
53  setPropertyFlags(c_ParallelProcessingCertified);
54  addParam("CrystalHadronEnergyThreshold", m_CrystalHadronEnergyThreshold,
55  "Hadron component energy threshold to identify as hadron digit.(GeV)", 0.003);
56  addParam("CrystalHadronIntensityThreshold", m_CrystalHadronIntensityThreshold,
57  "Hadron component intensity threshold to identify as hadron digit.", 0.005);
58  addParam("MVAidentifier", m_MVAidentifier, "MVA database identifier.", std::string{"eclClusterPSD_MVA"});
59 }
60 
61 // destructor
63 {
64 }
65 
66 // initialize MVA weightFile
67 void ECLClusterPSDModule::initializeMVAweightFile(const std::string& identifier,
68  std::unique_ptr<DBObjPtr<DatabaseRepresentationOfWeightfile>>& weightFileRepresentation)
69 {
70  if (not(boost::ends_with(identifier, ".root") or boost::ends_with(identifier, ".xml"))) {
71  weightFileRepresentation = std::unique_ptr<DBObjPtr<DatabaseRepresentationOfWeightfile>>(new
73  }
75 }
76 
77 // initialize
79 {
80  // ECL dataobjects
81  m_eclShowers.registerInDataStore(eclShowerArrayName());
82  m_eclCalDigits.registerInDataStore(eclCalDigitArrayName());
84 }
85 
86 // initialize MVA
87 void ECLClusterPSDModule::initializeMVA(const std::string& identifier,
88  std::unique_ptr<DBObjPtr<DatabaseRepresentationOfWeightfile>>& weightFileRepresentation, std::unique_ptr<MVA::Expert>& expert)
89 {
90  MVA::Weightfile weightfile;
91  //Load MVA weight file
92  if (weightFileRepresentation) {
93 
94  if (weightFileRepresentation->hasChanged()) {
95  std::stringstream ss((*weightFileRepresentation)->m_data);
96  weightfile = MVA::Weightfile::loadFromStream(ss);
97  } else
98  return;
99  } else {
100  weightfile = MVA::Weightfile::loadFromFile(identifier);
101  }
102 
103 
104  auto supported_interfaces = MVA::AbstractInterface::getSupportedInterfaces();
105  MVA::GeneralOptions general_options;
106  weightfile.getOptions(general_options);
107 
108  //Check number of variables in weight file
109  if (m_numMVAvariables != general_options.m_variables.size())
110  B2FATAL("Expecting " << m_numMVAvariables << " variables, found " << general_options.m_variables.size());
111 
112  expert = supported_interfaces[general_options.m_method]->getExpert();
113  expert->load(weightfile);
114 
115  //create new dataset
116  if (weightFileRepresentation == m_weightfile_representation) {
117  std::vector<float> dummy(general_options.m_variables.size(), 0);
118  m_dataset = std::unique_ptr<MVA::SingleDataset>(new MVA::SingleDataset(general_options, dummy, 0));
119  }
120 }
121 
122 // begin run
124 {
126 }
127 
128 // evaluates mva
130 {
131 
132  //geometry for cell id position
134 
135  auto relatedDigits = cluster->getRelationsTo<ECLCalDigit>();
136 
137  //EnergyToSort vector is used for sorting digits by offline two component energy
138  std::vector<std::tuple<double, unsigned int>> EnergyToSort;
139 
140  for (unsigned int iRel = 0; iRel < relatedDigits.size(); iRel++) {
141 
142  const auto caldigit = relatedDigits.object(iRel);
143 
144  //exclude digits without waveforms
145  const double digitChi2 = caldigit->getTwoComponentChi2();
146  if (digitChi2 < 0) continue;
147 
148  ECLDsp::TwoComponentFitType digitFitType1 = caldigit->getTwoComponentFitType();
149 
150  //exclude digits digits with poor chi2
151  if (digitFitType1 == ECLDsp::poorChi2) continue;
152 
153  //exclude digits with diode-crossing fits
154  if (digitFitType1 == ECLDsp::photonDiodeCrossing) continue;
155 
156  EnergyToSort.emplace_back(caldigit->getTwoComponentTotalEnergy(), iRel);
157 
158  }
159 
160  //sorting by energy
161  std::sort(EnergyToSort.begin(), EnergyToSort.end(), std::greater<>());
162 
163  //get cluster position information
164  const double showerR = cluster->getR();
165  const double showerTheta = cluster->getTheta();
166  const double showerPhi = cluster->getPhi();
167 
168  B2Vector3D showerPosition;
169  showerPosition.SetMagThetaPhi(showerR, showerTheta, showerPhi);
170 
171  size_t input_index{0};
172  auto& input = m_dataset->m_input;
173 
174  for (unsigned int digit = 0; digit < maxdigits; ++digit) {
175 
176  if (digit >= EnergyToSort.size()) break;
177 
178  const auto [digitEnergy, next] = EnergyToSort[digit];
179 
180  const auto caldigit = relatedDigits.object(next);
181  const double digitHadronEnergy = caldigit->getTwoComponentHadronEnergy();
182  const double digitOnlineEnergy = caldigit->getEnergy();
183  const double digitWeight = relatedDigits.weight(next);
184  ECLDsp::TwoComponentFitType digitFitType1 = caldigit->getTwoComponentFitType();
185  const int digitFitType = digitFitType1;
186  const int cellId = caldigit->getCellId();
187  B2Vector3D calDigitPosition = geometry->GetCrystalPos(cellId - 1);
188  TVector3 tempP = showerPosition - calDigitPosition;
189  const double Rval = tempP.Mag();
190  const double theVal = tempP.CosTheta();
191  const double phiVal = cos(tempP.Phi());
192 
193  input[input_index++] = theVal;
194  input[input_index++] = phiVal;
195  input[input_index++] = Rval;
196  input[input_index++] = digitOnlineEnergy;
197  input[input_index++] = digitEnergy;
198  input[input_index++] = (digitHadronEnergy / digitEnergy);
199  input[input_index++] = digitFitType;
200  input[input_index++] = digitWeight;
201 
202  }
203 
204  //fill remainder with defaults
205  while (input_index < input.size()) {
206  if (((input_index - 6) % 8) != 0) {
207  input[input_index++] = 0.0;
208  } else {
209  input[input_index++] = -1.0; //Fit Type
210  }
211  }
212 
213  //compute mva from input variables
214  const double MVAout = m_expert->apply(*m_dataset)[0];
215 
216  return MVAout;
217 }
218 
219 
221 {
222 
223  for (auto& shower : m_eclShowers) {
224 
225 
226  auto relatedDigits = shower.getRelationsTo<ECLCalDigit>();
227 
228  double cluster2CTotalEnergy = 0;
229  double cluster2CHadronEnergy = 0;
230  double numberofHadronDigits = 0;
231  double nWaveforminCluster = 0;
232 
233  for (unsigned int iRel = 0; iRel < relatedDigits.size(); iRel++) {
234 
235  const auto weight = relatedDigits.weight(iRel);
236 
237  const auto caldigit = relatedDigits.object(iRel);
238  const double digit2CChi2 = caldigit->getTwoComponentChi2();
239 
240  if (digit2CChi2 < 0) continue; //only digits with waveforms
241 
242  ECLDsp::TwoComponentFitType digitFitType1 = caldigit->getTwoComponentFitType();
243 
244  //exclude digits digits with poor chi2
245  if (digitFitType1 == ECLDsp::poorChi2) continue;
246 
247  //exclude digits with diode-crossing fits
248  if (digitFitType1 == ECLDsp::photonDiodeCrossing) continue;
249 
250  const double digit2CTotalEnergy = caldigit->getTwoComponentTotalEnergy();
251  const double digit2CHadronComponentEnergy = caldigit->getTwoComponentHadronEnergy();
252 
253  cluster2CTotalEnergy += digit2CTotalEnergy;
254  cluster2CHadronEnergy += digit2CHadronComponentEnergy;
255 
256  if (digit2CTotalEnergy < 0.6) {
257  if (digit2CHadronComponentEnergy > m_CrystalHadronEnergyThreshold) numberofHadronDigits += weight;
258  } else {
259  const double digitHadronComponentIntensity = digit2CHadronComponentEnergy / digit2CTotalEnergy;
260  if (digitHadronComponentIntensity > m_CrystalHadronIntensityThreshold) numberofHadronDigits += weight;
261  }
262 
263  nWaveforminCluster += weight;
264 
265  }
266 
267  if (nWaveforminCluster > 0) {
268  if (cluster2CTotalEnergy != 0) shower.setShowerHadronIntensity(cluster2CHadronEnergy / cluster2CTotalEnergy);
269  //evaluates mva classifier only if waveforms are available in the cluster
270  const double mvaout = evaluateMVA(&shower);
271  shower.setPulseShapeDiscriminationMVA(mvaout);
272 
273  shower.setNumberOfHadronDigits(numberofHadronDigits);
275 
276  } else {
277  shower.setShowerHadronIntensity(0);
278  shower.setPulseShapeDiscriminationMVA(0.5);
279  shower.setNumberOfHadronDigits(0);
280  }
281  }
282 }
283 
284 // end run
286 {
287 }
288 
289 // terminate
291 {
292 }
Belle2::ECLClusterPSDModule::m_MVAidentifier
std::string m_MVAidentifier
MVA - weight-file.
Definition: ECLClusterPSD.h:89
Belle2::ECLClusterPSDModule::evaluateMVA
double evaluateMVA(const ECLShower *cluster)
Evaluates mva.
Definition: ECLClusterPSD.cc:129
Belle2::ECLCalDigit::getTwoComponentChi2
double getTwoComponentChi2() const
Get two componnent chi2.
Definition: ECLCalDigit.h:154
Belle2::ECLClusterPSDModule::~ECLClusterPSDModule
~ECLClusterPSDModule()
Destructor.
Definition: ECLClusterPSD.cc:62
Belle2::MVA::AbstractInterface::getSupportedInterfaces
static std::map< std::string, AbstractInterface * > getSupportedInterfaces()
Returns interfaces supported by the MVA Interface.
Definition: Interface.h:55
Belle2::ECLCalDigit
Class to store calibrated ECLDigits: ECLCalDigits.
Definition: ECLCalDigit.h:38
Belle2::ECLDsp::TwoComponentFitType
TwoComponentFitType
Offline two component fit type.
Definition: ECLDsp.h:42
Belle2::MVA::Weightfile::getOptions
void getOptions(Options &options) const
Fills an Option object from the xml tree.
Definition: Weightfile.cc:76
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECLClusterPSDModule::m_numMVAvariables
const unsigned int m_numMVAvariables
number of variables expected in the MVA weightfile
Definition: ECLClusterPSD.h:88
Belle2::ECLClusterPSDModule::initializeMVAweightFile
void initializeMVAweightFile(const std::string &identifier, std::unique_ptr< DBObjPtr< DatabaseRepresentationOfWeightfile >> &weightFileRepresentation)
initialize MVA weight file from DB
Definition: ECLClusterPSD.cc:67
Belle2::MVA::Weightfile
The Weightfile class serializes all information about a training into an xml tree.
Definition: Weightfile.h:40
Belle2::ECLClusterPSDModule::m_eclShowers
StoreArray< ECLShower > m_eclShowers
ECLShower's.
Definition: ECLClusterPSD.h:85
Belle2::ECLClusterPSDModule::initialize
virtual void initialize() override
Initialize variables.
Definition: ECLClusterPSD.cc:78
Belle2::ECLShower::c_hasPulseShapeDiscrimination
@ c_hasPulseShapeDiscrimination
bit 2: Shower has pulse shape discrimination variables.
Definition: ECLShower.h:73
Belle2::ECLClusterPSDModule::beginRun
virtual void beginRun() override
begin run.
Definition: ECLClusterPSD.cc:123
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::B2Vector3< double >
Belle2::ECLClusterPSDModule::m_eclCalDigits
StoreArray< ECLCalDigit > m_eclCalDigits
ECLCalDigit's.
Definition: ECLClusterPSD.h:84
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition: ECLGeometryPar.cc:151
Belle2::ECLClusterPSDModule::m_dataset
std::unique_ptr< MVA::SingleDataset > m_dataset
Pointer to the current dataset.
Definition: ECLClusterPSD.h:93
Belle2::ECLDsp::poorChi2
@ poorChi2
All offline fit attempts were greater than chi2 threshold.
Definition: ECLDsp.h:43
Belle2::ECLClusterPSDModule
This module computes shower variables using pulse shape information from offline two component fits.
Definition: ECLClusterPSD.h:45
Belle2::ECLClusterPSDModule::terminate
virtual void terminate() override
terminate.
Definition: ECLClusterPSD.cc:290
Belle2::ECLClusterPSDModule::event
virtual void event() override
event per event.
Definition: ECLClusterPSD.cc:220
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLClusterPSDModule::m_weightfile_representation
std::unique_ptr< DBObjPtr< DatabaseRepresentationOfWeightfile > > m_weightfile_representation
Database pointer to the Database representation of the MVA weightfile.
Definition: ECLClusterPSD.h:91
Belle2::MVA::Weightfile::loadFromFile
static Weightfile loadFromFile(const std::string &filename)
Static function which loads a Weightfile from a file.
Definition: Weightfile.cc:215
Belle2::ECLDsp::photonDiodeCrossing
@ photonDiodeCrossing
photon + diode template fit
Definition: ECLDsp.h:46
Belle2::ECLClusterPSDModule::maxdigits
const unsigned int maxdigits
Max number of digits mva can include.
Definition: ECLClusterPSD.h:87
Belle2::MVA::GeneralOptions
General options which are shared by all MVA trainings.
Definition: Options.h:64
Belle2::MVA::AbstractInterface::initSupportedInterfaces
static void initSupportedInterfaces()
Static function which initliazes all supported interfaces, has to be called once before getSupportedI...
Definition: Interface.cc:55
Belle2::MVA::SingleDataset
Wraps the data of a single event into a Dataset.
Definition: Dataset.h:136
Belle2::B2Vector3::SetMagThetaPhi
void SetMagThetaPhi(DataType mag, DataType theta, DataType phi)
setter with mag, theta, phi
Definition: B2Vector3.h:258
Belle2::ECLClusterPSDModule::eclCalDigitArrayName
virtual const char * eclCalDigitArrayName() const
ECLCalDigits array name.
Definition: ECLClusterPSD.h:74
Belle2::ECLClusterPSDModule::eclShowerArrayName
virtual const char * eclShowerArrayName() const
ECLShowers array name.
Definition: ECLClusterPSD.h:78
Belle2::MVA::Weightfile::loadFromStream
static Weightfile loadFromStream(std::istream &stream)
Static function which deserializes a Weightfile from a stream.
Definition: Weightfile.cc:260
Belle2::ECLClusterPSDModule::m_expert
std::unique_ptr< MVA::Expert > m_expert
Pointer to the current MVA Expert.
Definition: ECLClusterPSD.h:92
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition: ECLGeometryPar.h:45
Belle2::ECLClusterPSDModule::m_CrystalHadronEnergyThreshold
double m_CrystalHadronEnergyThreshold
hadron component energy threshold to classify as hadron.
Definition: ECLClusterPSD.h:82
Belle2::ECLClusterPSDModule::endRun
virtual void endRun() override
end run.
Definition: ECLClusterPSD.cc:285
Belle2::ECLClusterPSDModule::m_CrystalHadronIntensityThreshold
double m_CrystalHadronIntensityThreshold
hadron component intensity threshold to classify as hadron.
Definition: ECLClusterPSD.h:83
Belle2::B2Vector3::Mag
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:158
Belle2::ECLClusterPSDModule::initializeMVA
void initializeMVA(const std::string &identifier, std::unique_ptr< DBObjPtr< DatabaseRepresentationOfWeightfile >> &weightFileRepresentation, std::unique_ptr< MVA::Expert > &expert)
Load MVA weight file and set pointer of expert.
Definition: ECLClusterPSD.cc:87
Belle2::ECLShower
Class to store ECL Showers.
Definition: ECLShower.h:42