Belle II Software development
ECLLocalMaximumFinderModule.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/* Own header. */
10#include <ecl/modules/eclLocalMaximumFinder/ECLLocalMaximumFinderModule.h>
11
12/* ECL headers. */
13#include <ecl/dataobjects/ECLCalDigit.h>
14#include <ecl/dataobjects/ECLConnectedRegion.h>
15#include <ecl/dataobjects/ECLDigit.h>
16#include <ecl/dataobjects/ECLElementNumbers.h>
17#include <ecl/dataobjects/ECLHit.h>
18#include <ecl/dataobjects/ECLLocalMaximum.h>
19#include <ecl/geometry/ECLGeometryPar.h>
20#include <ecl/geometry/ECLNeighbours.h>
21
22/* ROOT headers. */
23#include <TFile.h>
24#include <TTree.h>
25
26/* Basf2 headers. */
27#include <framework/gearbox/Const.h>
28#include <framework/logging/Logger.h>
29#include <mdst/dataobjects/MCParticle.h>
30
31// NAMESPACE(S)
32using namespace Belle2;
33using namespace ECL;
34
35//-----------------------------------------------------------------
36// Register the Module
37//-----------------------------------------------------------------
38REG_MODULE(ECLLocalMaximumFinder);
39REG_MODULE(ECLLocalMaximumFinderPureCsI);
40
41//-----------------------------------------------------------------
42// Implementation
43//-----------------------------------------------------------------
45 m_mcParticles(mcParticleArrayName()),
46 m_eclHits(eclHitArrayName()),
47 m_eclDigits(eclDigitArrayName()),
48 m_eclCalDigits(eclCalDigitArrayName()),
49 m_eclConnectedRegions(eclConnectedRegionArrayName()),
50 m_eclLocalMaximums(eclLocalMaximumArrayName())
51{
52 // Set description.
53 setDescription("ECLLocalMaximumFinderModule");
54
55 // Parallel processing certification.
57
58 // Add module parameters.
59 addParam("energyCut", m_energyCut, "Seed energy cut [MeV], minimum is 5.0 MeV.", 10.0 * Belle2::Unit::MeV);
60 addParam("isTrainingMode", m_isTrainingMode,
61 "Run in training mode (i.e. fill file with MVA input variables and determine MC truth of LM.).", 0);
62 addParam("outfileName", m_outfileName, "Output file name for training file.", std::string("ECLLocalMaximumFinderOutput.root"));
63 addParam("method", m_method, "Method to determine the LM (cut, none, fastbdt).", std::string("none"));
64 addParam("truthFraction", m_truthFraction, "Minimum matched energy fraction truth/rec for the LM.", 0.51);
65 addParam("cutOffset", m_cutOffset, "Cut method specific: Offset. (BaBar: 2.5, high eff: 1.40)", 2.5);
66 addParam("cutSlope", m_cutSlope, "Cut method specific: Slope. (BaBar: 0.5, high eff: 3.0)", 0.5);
67 addParam("cutRatioCorrection", m_cutRatioCorrection, "Cut method specific: Ratio correction.", 0.0);
68
69}
70
72{
73 ;
74}
75
77{
78 B2DEBUG(200, "ECLLocalMaximumFinderModule::initialize()");
79
82 m_eclConnectedRegions.requireRelationTo(m_eclCalDigits);
83 m_eclLocalMaximums.registerInDataStore(eclLocalMaximumArrayName());
84 m_eclConnectedRegions.registerRelationTo(m_eclLocalMaximums);
85 m_eventLevelClusteringInfo.isRequired();
86
87 // Check user input.
89 B2WARNING("ECLLocalMaximumFinderModule::initialize: Energy threshold too small, resetting to " << c_minEnergyCut << " GeV");
91 }
92
93 if (m_truthFraction < 1e-6) {
94 B2WARNING("ECLLocalMaximumFinderModule::initialize: Matching fraction must be above 1e-6, input: " << m_truthFraction <<
95 " will be reset to 1e-6.");
96 m_truthFraction = 1e-6;
97 }
98
99 if (m_truthFraction > 1.0) {
100 B2WARNING("ECLLocalMaximumFinderModule::initialize: Matching fraction must be below 1.0, input: " << m_truthFraction <<
101 " will be reset to 1.0.");
102 m_truthFraction = 1.0;
103 }
104
105 // Geometry instance.
107
108 // Initialize neighbour map.
109 m_neighbourMap = new ECLNeighbours("N", 1);
110
111 // Reset all variables.
113
114 // Open output files and declare branches if in training mode. Each file will hold a flat ntuple of training data.
115 m_outfile = nullptr;
116 m_tree = nullptr;
117 if (m_isTrainingMode > 0) {
118 const int cBufferLength = 500;
119 char tmpBuffer[cBufferLength];
120 int tmpBufferLength = sprintf(tmpBuffer, "%s", m_outfileName.c_str());
121
122 if (tmpBufferLength >= cBufferLength) {
123 B2FATAL("ECLLocalMaximumFinderModule::initialize(): Output training filename length too long!");
124 }
125
126 m_outfile = new TFile(tmpBuffer, "RECREATE");
127 m_tree = new TTree("locmax", "locmax");
128 m_tree->Branch("energyRatioNeighbour0", &m_energyRatioNeighbour[0], "energyRatioNeighbour0/F");
129 m_tree->Branch("energyRatioNeighbour1", &m_energyRatioNeighbour[1], "energyRatioNeighbour1/F");
130 m_tree->Branch("energyRatioNeighbour2", &m_energyRatioNeighbour[2], "energyRatioNeighbour2/F");
131 m_tree->Branch("energyRatioNeighbour3", &m_energyRatioNeighbour[3], "energyRatioNeighbour3/F");
132 m_tree->Branch("energyRatioNeighbour4", &m_energyRatioNeighbour[4], "energyRatioNeighbour4/F");
133 m_tree->Branch("energyRatioNeighbour5", &m_energyRatioNeighbour[5], "energyRatioNeighbour5/F");
134 m_tree->Branch("energyRatioNeighbour6", &m_energyRatioNeighbour[6], "energyRatioNeighbour6/F");
135 m_tree->Branch("energyRatioNeighbour7", &m_energyRatioNeighbour[7], "energyRatioNeighbour7/F");
136 m_tree->Branch("energyRatioNeighbour8", &m_energyRatioNeighbour[8], "energyRatioNeighbour8/F");
137 m_tree->Branch("energyRatioNeighbour9", &m_energyRatioNeighbour[9], "energyRatioNeighbour9/F");
138 m_tree->Branch("energyRatioNeighbour10", &m_energyRatioNeighbour[10], "energyRatioNeighbour10/F");
139 m_tree->Branch("energyRatioNeighbour11", &m_energyRatioNeighbour[11], "energyRatioNeighbour11/F");
140 m_tree->Branch("energy", &m_energy, "energy/F");
141 m_tree->Branch("target", &m_target, "target/F");
142 m_tree->Branch("targetindex", &m_targetindex, "targetindex/F");
143 m_tree->Branch("targetpi0index", &m_targetpi0index, "targetpi0index/F");
144 m_tree->Branch("thetaId", &m_thetaId, "thetaId/F");
145 m_tree->Branch("phiId", &m_phiId, "phiId/F");
146 m_tree->Branch("cellId", &m_cellId, "cellId/F");
147 m_tree->Branch("maxNeighbourEnergy", &m_maxNeighbourEnergy, "maxNeighbourEnergy/F");
148 m_tree->Branch("nNeighbours10", &m_nNeighbours10, "nNeighbours10/F");
149 m_tree->Branch("time", &m_time, "time/F");
150 m_tree->Branch("timeResolution", &m_timeResolution, "timeResolution/F");
151 m_tree->Branch("timeFitFailed", &m_timeFitFailed, "timeFitFailed/F");
152 m_tree->Branch("CRId", &m_CRId, "CRId/F");
153 m_tree->Branch("LMId", &m_LMId, "LMId/F");
154 }
155
156 // initialize the vector that gives the relation between cellid and store array position
158
159}
160
162{
163 ;
164}
165
167{
168 B2DEBUG(200, "ECLLocalMaximumFinderModule::event()");
169
170 // Fill a vector that can be used to map cellid -> store array position
171 std::fill_n(m_StoreArrPosition.begin(), m_StoreArrPosition.size(), -1);
172 for (int i = 0; i < m_eclCalDigits.getEntries(); i++) {
173 m_StoreArrPosition[m_eclCalDigits[i]->getCellId()] = i;
174 }
175
176 // Vector with neighbour ids.
177 std::vector< double > vNeighourEnergies;
178 vNeighourEnergies.resize(c_nMaxNeighbours);
179
180 // Loop over connected regions.
181 for (const ECLConnectedRegion& aCR : m_eclConnectedRegions) {
182 const int crId = aCR.getCRId();
183 int iLM = 1;
184
185 // Loop over all entries in this CR.
186 for (const ECLCalDigit& aECLCalDigit : aCR.getRelationsTo<ECLCalDigit>()) {
187 // Check seed energy cut.
188 if (aECLCalDigit.getEnergy() >= m_energyCut) {
189
190 // Clean up for this candidate (MVA is trained per LM, regardless of CR)
191 std::fill_n(vNeighourEnergies.begin(), vNeighourEnergies.size(),
192 -999); // -999 means later: this digit is just not available in this neighbour definition.
195
196 // Check neighbours: Must be a local energy maximum.
197 bool isLocMax = 1;
198 int neighbourCount = 0;
199 for (auto& neighbourId : m_neighbourMap->getNeighbours(aECLCalDigit.getCellId())) {
200 if (neighbourId == aECLCalDigit.getCellId()) continue; // Skip the center cell to avoid possible floating point issues.
201
202 const int pos = m_StoreArrPosition[neighbourId]; // Get position in the store array for this digit.
203
204 double energyNeighbour = 0.0;
205 if (pos >= 0) {
206 energyNeighbour = m_eclCalDigits[pos]->getEnergy(); // Get the energy directly from the store array.
207 vNeighourEnergies[neighbourCount] = energyNeighbour;
208 } else {
209 // Digit does not belong to this CR
210 vNeighourEnergies[neighbourCount] = 0.0;
211 }
212 ++neighbourCount;
213
214 if (energyNeighbour > aECLCalDigit.getEnergy()) {
215 isLocMax = 0;
216 break;
217 }
218 }
219
220 // It is a local maximum. Get all variables needed for classification.
221 if (isLocMax) {
222
223 for (unsigned int npos = 0; npos < vNeighourEnergies.size(); ++npos) {
224 if (vNeighourEnergies[npos] >= 0) {
225 m_energyRatioNeighbour[npos] = static_cast < float >(vNeighourEnergies[npos] / aECLCalDigit.getEnergy());
226 if (vNeighourEnergies[npos] > m_maxNeighbourEnergy) m_maxNeighbourEnergy = vNeighourEnergies[npos];
227 if (vNeighourEnergies[npos] > 1.0 * Belle2::Unit::MeV) ++m_nNeighbours10;
228 } else m_energyRatioNeighbour[npos] = 0.0;
229 }
230 m_time = aECLCalDigit.getTime();
231 m_maxEnergyRatio = m_maxNeighbourEnergy / aECLCalDigit.getEnergy();
232
233 // Fill training monitoring variables and MC truth information.
234 if (m_isTrainingMode > 0) {
235
236 m_energy = static_cast < float >(aECLCalDigit.getEnergy());
237 m_cellId = static_cast < float >(aECLCalDigit.getCellId());
238 m_timeResolution = static_cast < float >(aECLCalDigit.getTimeResolution());
239 m_timeFitFailed = static_cast < float >(aECLCalDigit.isFailedFit());
240 m_CRId = static_cast < float >(crId);
241 m_LMId = static_cast < float >(iLM);
242
243 m_geom->Mapping(m_cellId - 1);
244 m_thetaId = static_cast < float >(m_geom->GetThetaID());
245 m_phiId = static_cast < float >(m_geom->GetPhiID());
246
247 // This requires MC matching before this stage!
248 int pi0index = -1;
249 int maxtype = 0;
250 int maxpos = 0;
251
252 auto relatedParticlePairs = aECLCalDigit.getRelationsWith<MCParticle>();
253
254 for (unsigned int irel = 0; irel < relatedParticlePairs.size(); irel++) {
255 const auto particle = relatedParticlePairs.object(irel);
256 const double weight = relatedParticlePairs.weight(irel);
257
258 int motherpdg = -1;
259 int motherindex = -1;
260 pi0index = -1;
261 getEnteringMother(*particle, motherpdg, motherindex, pi0index);
262 addToSignalEnergy(motherpdg, motherindex, pi0index, weight);
263 }
264
265 B2DEBUG(175, " -> digt energy: " << aECLCalDigit.getEnergy());
266 B2DEBUG(175, "photon: " << m_signalEnergy[0][0] << " " << m_signalEnergy[0][1]);
267 B2DEBUG(175, "pi0: " << m_signalEnergy[1][0] << " " << m_signalEnergy[1][1]);
268 B2DEBUG(175, "electron: " << m_signalEnergy[2][0] << " " << m_signalEnergy[2][1]);
269 B2DEBUG(175, "muon: " << m_signalEnergy[3][0] << " " << m_signalEnergy[3][1]);
270 B2DEBUG(175, "neutral hadron: " << m_signalEnergy[4][0] << " " << m_signalEnergy[4][1]);
271 B2DEBUG(175, "charged hadron: " << m_signalEnergy[5][0] << " " << m_signalEnergy[5][1]);
272 B2DEBUG(175, "other: " << m_signalEnergy[6][0]);
273
274 maxtype = 0;
275 maxpos = 0;
276 getMax(maxtype, maxpos);
277
278 if (maxtype >= 0) {
279 if (m_signalEnergy[maxtype][maxpos] >= m_truthFraction * aECLCalDigit.getEnergy()) {
280 m_target = maxtype;
281 m_targetindex = m_signalId[maxtype][maxpos];
282 m_targetpi0index = pi0index;
283 } else {
284 m_target = 7;
285 }
286 } else {
287 m_target = 7;
288 }
289
290 m_tree->Fill();
291
292 } // end training
293
294 if (m_method == "cut") {
295
296 B2DEBUG(200, "m_cutSlope: " << m_cutSlope << ", m_nNeighbours10: " << m_nNeighbours10 << ", m_cutOffset: " << m_cutOffset <<
297 ", m_maxNeighbourEnergy: " << m_maxNeighbourEnergy << ", m_cutRatioCorrection: " << m_cutRatioCorrection <<
298 ", aECLCalDigit.getEnergy(): " << aECLCalDigit.getEnergy());
299 B2DEBUG(200, "m_cutSlope * (m_nNeighbours10 - m_cutOffset): " << m_cutSlope * (m_nNeighbours10 - m_cutOffset));
300 B2DEBUG(200, "(m_maxNeighbourEnergy - m_cutRatioCorrection) / (aECLCalDigit.getEnergy() - m_cutRatioCorrection)" <<
301 (m_maxNeighbourEnergy - m_cutRatioCorrection) / (aECLCalDigit.getEnergy() - m_cutRatioCorrection) << "\n");
302
304 (aECLCalDigit.getEnergy() - m_cutRatioCorrection)) {
305 makeLocalMaximum(aCR, aECLCalDigit.getCellId(), iLM);
306 ++iLM;
307 }
308 } else if (m_method == "none") { // All energy local maximums will become local maximums.
309 makeLocalMaximum(aCR, aECLCalDigit.getCellId(), iLM);
310 ++iLM;
311 }
312
313 }
314
315 } // end check energy
316
317 } // end CalDigit loop
318
319 // Check if there is at least one local maximum in the CR. If not, make the highest energetic crystal one.
320 if (iLM == 1) {
321
322 int highestEnergyCellId = -1;
323 double highestEnergy = 0.0;
324
325 // Loop over all entries in this CR.
326 for (const ECLCalDigit& aECLCalDigit : aCR.getRelationsTo<ECLCalDigit>(eclCalDigitArrayName())) {
327 if (aECLCalDigit.getEnergy() > highestEnergy) {
328 highestEnergyCellId = aECLCalDigit.getCellId();
329 highestEnergy = aECLCalDigit.getEnergy();
330 }
331 } // end CalDigit loop
332
333 makeLocalMaximum(aCR, highestEnergyCellId, 0);
334
335 }
336
337 } // end CR loop
338
339 // Find the number of local maximums in each ECL region in mdst
340 uint16_t nLMPerRegion[3] = {};
341 for (const ECLLocalMaximum& aLM : m_eclLocalMaximums) {
342 const int iCellId = aLM.getCellId();
343 if (ECLElementNumbers::isForward(iCellId)) {nLMPerRegion[0]++;}
344 if (ECLElementNumbers::isBarrel(iCellId)) {nLMPerRegion[1]++;}
345 if (ECLElementNumbers::isBackward(iCellId)) {nLMPerRegion[2]++;}
346 }
347
348 // Store numbers in EventLevelClusteringInfo mdst object
350 m_eventLevelClusteringInfo->setNECLLocalMaximumsFWD(nLMPerRegion[0]);
351 m_eventLevelClusteringInfo->setNECLLocalMaximumsBarrel(nLMPerRegion[1]);
352 m_eventLevelClusteringInfo->setNECLLocalMaximumsBWD(nLMPerRegion[2]);
353
354}
355
356
358{
359 B2DEBUG(200, "ECLLocalMaximumFinderModule::endRun()");
360}
361
362
364{
365 B2DEBUG(200, "ECLLocalMaximumFinderModule::terminate()");
366
367 // If run in trainingmode, write the training file.
368 if (m_outfile and m_tree) {
369 m_outfile->cd();
370 m_tree->Write();
371 m_outfile->Write();
372 m_outfile->Close();
373 delete m_outfile;
374 }
375
376 if (m_neighbourMap) delete m_neighbourMap;
377
378}
379
380void ECLLocalMaximumFinderModule::makeLocalMaximum(const ECLConnectedRegion& aCR, const int cellId, const int lmId)
381{
382 // Set the local maximum dataobject.
383 const auto aLocalMaximum = m_eclLocalMaximums.appendNew();
384
385 B2DEBUG(175, "ECLLocalMaximumFinderModule::makeLocalMaximum(): local maximum cellid: " << cellId);
386
387 // Set the id of this local maximum.
388 aLocalMaximum->setLMId(lmId);
389
390 // Set the cell Id of the digit.
391 aLocalMaximum->setCellId(cellId);
392
393 // Add relations to ECLConnectedRegion.
394 aCR.addRelationTo(aLocalMaximum);
395
396}
397
399{
400 m_target = -1;
401 m_targetindex = -1;
402 m_targetpi0index = -1;
403
405
406 for (unsigned int i = 0; i < 10; ++i) {
407 for (unsigned int j = 0; j < 5; ++j) {
408 m_signalEnergy[i][j] = 0.;
409 m_signalId[i][j] = -1;
410 }
411 }
412
413}
414
416{
417
418 for (unsigned i = 0; i < c_nMaxNeighbours; ++i) {
420 }
421
422 m_energy = 0.;
423 m_cellId = -1;
424 m_thetaId = -1;
425 m_phiId = -1;
426 m_nNeighbours10 = 0;
428 m_time = -9999;
429 m_timeResolution = -9999.;
430 m_timeFitFailed = -9999.;
431 m_CRId = -1;
432 m_LMId = -1;
433 m_maxEnergyRatio = 0.0;
434}
435
436
437void ECLLocalMaximumFinderModule::getEnteringMother(const MCParticle& particle, int& pdg, int& arrayindex, int& pi0arrayindex)
438{
439
440 int index = particle.getArrayIndex();
441 int pi0index = -1;
442
443 while (!isEnteringECL(m_mcParticles[index]->getProductionVertex())) {
444 if (m_mcParticles[index]->getMother()) index = m_mcParticles[index]->getMother()->getArrayIndex();
445 else index = -1;
446 };
447
448 // For photon mother: are they from a pi0? This can be used to improved overlap/merged pi0 reconstruction.
449 if (m_mcParticles[index]->getPDG() == Const::photon.getPDGCode()) {
450 if (m_mcParticles[index]->getMother()->getPDG() == Const::pi0.getPDGCode()) {
451 pi0index = m_mcParticles[index]->getMother()->getArrayIndex();
452 }
453 }
454
455 // Dont include mother if its energy is too low or if its a photon from a neutron interaction.
456 if ((m_mcParticles[index]->getEnergy() < 5.0 * Belle2::Unit::MeV)
457 or (m_mcParticles[index]->getPDG() == Const::photon.getPDGCode()
458 and abs(m_mcParticles[index]->getMother()->getPDG()) == Const::proton.getPDGCode())) {
459 pdg = -1;
460 arrayindex = -1;
461 pi0arrayindex = -1;
462 } else {
463 pdg = m_mcParticles[index]->getPDG();
464 arrayindex = index;
465 pi0arrayindex = pi0index;
466 }
467}
468
470{
471
472 const double theta = vertex.Theta();
473
474 if (theta > 0.555015 and theta < 2.26369) { //barrel
475 double radius = vertex.Perp();
476 if (radius < 125 * Belle2::Unit::cm) return true;
477 } else if (theta <= 0.555015) { //fwd
478 if (vertex.Z() < 196.16 * Belle2::Unit::cm) return true;
479 } else if (theta >= 2.26369) { //bwd
480 if (vertex.Z() > -102.16 * Belle2::Unit::cm) return true;
481 }
482
483 return false;
484}
485
486int ECLLocalMaximumFinderModule::getIdPosition(const int type, const int id)
487{
488
489 for (int i = 0; i < 5; i++) {
490 if (m_signalId[type][i] == id) return i;
491 if (m_signalId[type][i] == -1) {
492 m_signalId[type][i] = id;
493 return i; //next free one
494 }
495 }
496
497 return -1;
498}
499
500// of all entries, get the maximum!
502{
503 double maxe = 0.;
504 int maxtype = -1;
505 int maxid = -1;
506
507 for (unsigned int i = 0; i < 10; ++i) {
508 for (unsigned int j = 0; j < 5; ++j) {
509 if (m_signalEnergy[i][j] > maxe) {
510 maxe = m_signalEnergy[i][j];
511 maxtype = i;
512 maxid = j;
513 }
514 }
515 }
516
517 type = maxtype;
518 id = maxid;
519
520}
521
522void ECLLocalMaximumFinderModule::addToSignalEnergy(int motherpdg, int motherindex, int pi0index, double weight)
523{
524
525 // for the LM training and CR/LM debugging
526 if (motherpdg == Const::photon.getPDGCode()) {
527 if (pi0index >= 0) { // photon from pi0
528 int idpos = getIdPosition(1, motherindex);
529 m_signalEnergy[1][idpos] += weight;
530 } else {
531 int idpos = getIdPosition(0, motherindex);
532 m_signalEnergy[0][idpos] += weight; // photon from another source
533 }
534 } else if (abs(motherpdg) == Const::electron.getPDGCode()) { // electron
535 int idpos = getIdPosition(2, motherindex);
536 m_signalEnergy[2][idpos] += weight;
537 } else if (abs(motherpdg) == Const::muon.getPDGCode()) { // muon
538 int idpos = getIdPosition(3, motherindex);
539 m_signalEnergy[3][idpos] += weight;
540 } else if (abs(motherpdg) == Const::Klong.getPDGCode() or abs(motherpdg) == Const::neutron.getPDGCode()) { // neutral hadron
541 int idpos = getIdPosition(4, motherindex);
542 m_signalEnergy[4][idpos] += weight;
543 } else if (abs(motherpdg) == Const::pion.getPDGCode() or abs(motherpdg) == Const::kaon.getPDGCode()
544 or abs(motherpdg) == Const::proton.getPDGCode()) { // charged hadron
545 int idpos = getIdPosition(5, motherindex);
546 m_signalEnergy[5][idpos] += weight;
547 } else { // everything else
548 int idpos = getIdPosition(6, motherindex);
549 m_signalEnergy[6][idpos] += weight;
550 }
551}
static const ParticleType neutron
neutron particle
Definition: Const.h:675
static const ParticleType pi0
neutral pion particle
Definition: Const.h:674
static const ChargedStable muon
muon particle
Definition: Const.h:660
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:678
static const ChargedStable proton
proton particle
Definition: Const.h:663
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:662
static const ParticleType photon
photon particle
Definition: Const.h:673
static const ChargedStable electron
electron particle
Definition: Const.h:659
Class to store calibrated ECLDigits: ECLCalDigits.
Definition: ECLCalDigit.h:23
Class to store connected regions (CRs)
TTree * m_tree
tree that contain information for MVA training
double m_cutRatioCorrection
correction for nominator and denominator of the ratio.
double m_signalEnergy[10][5]
total energy per MC matching type of this digit
void getEnteringMother(const MCParticle &particle, int &pdg, int &arrayindex, int &pi0arrayindex)
Get enterging mother of this particle.
StoreArray< ECLConnectedRegion > m_eclConnectedRegions
Store array: ECLConnectedRegion.
float m_energyRatioNeighbour[c_nMaxNeighbours]
energy ratio of neighbour 0..9 to center
double m_totalSignalEnergy
total energy of this digit
float m_phiId
local maximum center theta Id
virtual void initialize() override
Initialize.
StoreArray< ECLLocalMaximum > m_eclLocalMaximums
Store array: ECLLocalMaximum.
float m_cellId
local maximum center cell Id
virtual const char * eclLocalMaximumArrayName() const
Name to be used for default option: ECLLocalMaximums.
int getIdPosition(const int type, const int id)
Get Id position in the vector.
float m_energy
Variables to monitor the MVA training.
std::string m_method
Method to find the local maximum.
virtual void terminate() override
Terminate (close ROOT files here if you have opened any).
std::string m_outfileName
file name prefix of the training output file
ECL::ECLNeighbours * m_neighbourMap
Neighbour maps.
void makeLocalMaximum(const ECLConnectedRegion &aCR, const int cellId, const int lmId)
Make local maximum for a given connected region.
void resetClassifierVariables()
Reset Classifier Variables.
float m_thetaId
local maximum center theta Id
float m_maxNeighbourEnergy
highest energy of all neighbours
bool isEnteringECL(const B2Vector3D &vec)
Check if particle is produced outside of the ECL.
TFile * m_outfile
Output training files and trees.
void addToSignalEnergy(int motherpdg, int motherindex, int pi0index, double weight)
Add energy to vector.
virtual const char * eclConnectedRegionArrayName() const
Name to be used for default option: ECLConnectedRegions.
static const unsigned c_nMaxNeighbours
Variables (possibly) used for MVA classification.
StoreArray< MCParticle > m_mcParticles
Store array: MCParticle.
virtual const char * eclCalDigitArrayName() const
Name to be used for default or PureCsI option: ECLCalDigits.
float m_maxEnergyRatio
Highest energetic neighbour energy divided by LM energy.
int m_signalId[10][5]
total energy per MC matching type of this digit
const double c_minEnergyCut
Minimum LM energy.
std::vector< int > m_StoreArrPosition
vector (ECLElementNumbers::c_NCrystals + 1 entries) with cell id to store array positions
StoreObjPtr< EventLevelClusteringInfo > m_eventLevelClusteringInfo
EventLevelClusteringInfo.
void resetTrainingVariables()
Reset Debug Variables.
void getMax(int &maxtype, int &maxpos)
Get the highest energy deposition particle type.
StoreArray< ECLCalDigit > m_eclCalDigits
Store array: ECLCalDigit.
int m_isTrainingMode
training mode for MVA methods (i.e.
float m_nNeighbours10
Variables (possibly) used for cut classification.
Class to store local maxima (LM)
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
void Mapping(int cid)
Mapping theta, phi Id.
int GetThetaID()
Get Theta Id.
Class to get the neighbours for a given cell id.
Definition: ECLNeighbours.h:25
const std::vector< short int > & getNeighbours(short int cid) const
Return the neighbours for a given cell ID.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Base class for Modules.
Definition: Module.h:72
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
@ 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
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
static const double MeV
[megaelectronvolt]
Definition: Unit.h:114
static const double cm
Standard units with the value = 1.
Definition: Unit.h:47
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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
bool isForward(int cellId)
Check whether the crystal is in forward ECL.
bool isBarrel(int cellId)
Check whether the crystal is in barrel ECL.
const int c_NCrystals
Number of crystals.
bool isBackward(int cellId)
Check whether the crystal is in backward ECL.
Abstract base class for different kinds of events.