18 #include <ecl/modules/eclLocalMaximumFinder/ECLLocalMaximumFinderModule.h>
25 #include <framework/datastore/StoreArray.h>
26 #include <framework/logging/Logger.h>
29 #include <mdst/dataobjects/MCParticle.h>
32 #include <ecl/dataobjects/ECLHit.h>
33 #include <ecl/dataobjects/ECLDigit.h>
34 #include <ecl/dataobjects/ECLCalDigit.h>
35 #include <ecl/dataobjects/ECLLocalMaximum.h>
36 #include <ecl/dataobjects/ECLConnectedRegion.h>
37 #include <ecl/geometry/ECLNeighbours.h>
38 #include <ecl/geometry/ECLGeometryPar.h>
54 m_mcParticles(mcParticleArrayName()),
55 m_eclHits(eclHitArrayName()),
56 m_eclDigits(eclDigitArrayName()),
57 m_eclCalDigits(eclCalDigitArrayName()),
58 m_eclConnectedRegions(eclConnectedRegionArrayName()),
59 m_eclLocalMaximums(eclLocalMaximumArrayName())
62 setDescription(
"ECLLocalMaximumFinderModule");
65 setPropertyFlags(c_ParallelProcessingCertified);
68 addParam(
"energyCut", m_energyCut,
"Seed energy cut [MeV], minimum is 5.0 MeV.", 10.0 *
Belle2::Unit::MeV);
69 addParam(
"isTrainingMode", m_isTrainingMode,
70 "Run in training mode (i.e. fill file with MVA input variables and determine MC truth of LM.).", 0);
71 addParam(
"outfileName", m_outfileName,
"Output file name for training file.", std::string(
"ECLLocalMaximumFinderOutput.root"));
72 addParam(
"method", m_method,
"Method to determine the LM (cut, none, fastbdt).", std::string(
"none"));
73 addParam(
"truthFraction", m_truthFraction,
"Minimum matched energy fraction truth/rec for the LM.", 0.51);
74 addParam(
"cutOffset", m_cutOffset,
"Cut method specific: Offset. (BaBar: 2.5, high eff: 1.40)", 2.5);
75 addParam(
"cutSlope", m_cutSlope,
"Cut method specific: Slope. (BaBar: 0.5, high eff: 3.0)", 0.5);
76 addParam(
"cutRatioCorrection", m_cutRatioCorrection,
"Cut method specific: Ratio correction.", 0.0);
87 B2DEBUG(200,
"ECLLocalMaximumFinderModule::initialize()");
89 m_eclLocalMaximums.registerInDataStore(eclLocalMaximumArrayName());
90 m_eclCalDigits.registerInDataStore(eclCalDigitArrayName());
91 m_eclConnectedRegions.registerInDataStore(eclConnectedRegionArrayName());
92 m_eclConnectedRegions.registerRelationTo(m_eclLocalMaximums);
93 m_eclCalDigits.registerRelationTo(m_eclLocalMaximums);
96 if (m_energyCut < c_minEnergyCut) {
97 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Energy threshold too small, resetting to " << c_minEnergyCut <<
" GeV");
98 m_energyCut = c_minEnergyCut;
101 if (m_truthFraction < 1e-6) {
102 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Matching fraction must be above 1e-6, input: " << m_truthFraction <<
103 " will be reset to 1e-6.");
104 m_truthFraction = 1e-6;
107 if (m_truthFraction > 1.0) {
108 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Matching fraction must be below 1.0, input: " << m_truthFraction <<
109 " will be reset to 1.0.");
110 m_truthFraction = 1.0;
120 resetClassifierVariables();
125 if (m_isTrainingMode > 0) {
126 const int cBufferLength = 500;
127 char tmpBuffer[cBufferLength];
128 int tmpBufferLength = sprintf(tmpBuffer,
"%s", m_outfileName.c_str());
130 if (tmpBufferLength >= cBufferLength) {
131 B2FATAL(
"ECLLocalMaximumFinderModule::initialize(): Output training filename length too long!");
134 m_outfile =
new TFile(tmpBuffer,
"RECREATE");
135 m_tree =
new TTree(
"locmax",
"locmax");
136 m_tree->Branch(
"energyRatioNeighbour0", &m_energyRatioNeighbour[0],
"energyRatioNeighbour0/F");
137 m_tree->Branch(
"energyRatioNeighbour1", &m_energyRatioNeighbour[1],
"energyRatioNeighbour1/F");
138 m_tree->Branch(
"energyRatioNeighbour2", &m_energyRatioNeighbour[2],
"energyRatioNeighbour2/F");
139 m_tree->Branch(
"energyRatioNeighbour3", &m_energyRatioNeighbour[3],
"energyRatioNeighbour3/F");
140 m_tree->Branch(
"energyRatioNeighbour4", &m_energyRatioNeighbour[4],
"energyRatioNeighbour4/F");
141 m_tree->Branch(
"energyRatioNeighbour5", &m_energyRatioNeighbour[5],
"energyRatioNeighbour5/F");
142 m_tree->Branch(
"energyRatioNeighbour6", &m_energyRatioNeighbour[6],
"energyRatioNeighbour6/F");
143 m_tree->Branch(
"energyRatioNeighbour7", &m_energyRatioNeighbour[7],
"energyRatioNeighbour7/F");
144 m_tree->Branch(
"energyRatioNeighbour8", &m_energyRatioNeighbour[8],
"energyRatioNeighbour8/F");
145 m_tree->Branch(
"energyRatioNeighbour9", &m_energyRatioNeighbour[9],
"energyRatioNeighbour9/F");
146 m_tree->Branch(
"energyRatioNeighbour10", &m_energyRatioNeighbour[10],
"energyRatioNeighbour10/F");
147 m_tree->Branch(
"energyRatioNeighbour11", &m_energyRatioNeighbour[11],
"energyRatioNeighbour11/F");
148 m_tree->Branch(
"energy", &m_energy,
"energy/F");
149 m_tree->Branch(
"target", &m_target,
"target/F");
150 m_tree->Branch(
"targetindex", &m_targetindex,
"targetindex/F");
151 m_tree->Branch(
"targetpi0index", &m_targetpi0index,
"targetpi0index/F");
152 m_tree->Branch(
"thetaId", &m_thetaId,
"thetaId/F");
153 m_tree->Branch(
"phiId", &m_phiId,
"phiId/F");
154 m_tree->Branch(
"cellId", &m_cellId,
"cellId/F");
155 m_tree->Branch(
"maxNeighbourEnergy", &m_maxNeighbourEnergy,
"maxNeighbourEnergy/F");
156 m_tree->Branch(
"nNeighbours10", &m_nNeighbours10,
"nNeighbours10/F");
157 m_tree->Branch(
"time", &m_time,
"time/F");
158 m_tree->Branch(
"timeResolution", &m_timeResolution,
"timeResolution/F");
159 m_tree->Branch(
"timeFitFailed", &m_timeFitFailed,
"timeFitFailed/F");
160 m_tree->Branch(
"CRId", &m_CRId,
"CRId/F");
161 m_tree->Branch(
"LMId", &m_LMId,
"LMId/F");
165 m_StoreArrPosition.resize(8736 + 1);
176 B2DEBUG(200,
"ECLLocalMaximumFinderModule::event()");
179 std::fill_n(m_StoreArrPosition.begin(), m_StoreArrPosition.size(), -1);
180 for (
int i = 0; i < m_eclCalDigits.getEntries(); i++) {
181 m_StoreArrPosition[m_eclCalDigits[i]->getCellId()] = i;
185 std::vector< double > vNeighourEnergies;
186 vNeighourEnergies.resize(c_nMaxNeighbours);
190 const int crId = aCR.getCRId();
197 if (aECLCalDigit.getEnergy() >= m_energyCut) {
200 std::fill_n(vNeighourEnergies.begin(), vNeighourEnergies.size(),
202 resetTrainingVariables();
203 resetClassifierVariables();
207 int neighbourCount = 0;
208 for (
auto& neighbourId : m_neighbourMap->getNeighbours(aECLCalDigit.getCellId())) {
209 if (neighbourId == aECLCalDigit.getCellId())
continue;
211 const int pos = m_StoreArrPosition[neighbourId];
213 double energyNeighbour = 0.0;
215 energyNeighbour = m_eclCalDigits[pos]->getEnergy();
216 vNeighourEnergies[neighbourCount] = energyNeighbour;
219 vNeighourEnergies[neighbourCount] = 0.0;
223 if (energyNeighbour > aECLCalDigit.getEnergy()) {
232 for (
unsigned int npos = 0; npos < vNeighourEnergies.size(); ++npos) {
233 if (vNeighourEnergies[npos] >= 0) {
234 m_energyRatioNeighbour[npos] =
static_cast < float >(vNeighourEnergies[npos] / aECLCalDigit.getEnergy());
235 if (vNeighourEnergies[npos] > m_maxNeighbourEnergy) m_maxNeighbourEnergy = vNeighourEnergies[npos];
237 }
else m_energyRatioNeighbour[npos] = 0.0;
239 m_time = aECLCalDigit.getTime();
240 m_maxEnergyRatio = m_maxNeighbourEnergy / aECLCalDigit.getEnergy();
243 if (m_isTrainingMode > 0) {
245 m_energy =
static_cast < float >(aECLCalDigit.getEnergy());
246 m_cellId =
static_cast < float >(aECLCalDigit.getCellId());
247 m_timeResolution =
static_cast < float >(aECLCalDigit.getTimeResolution());
248 m_timeFitFailed =
static_cast < float >(aECLCalDigit.isFailedFit());
249 m_CRId =
static_cast < float >(crId);
250 m_LMId =
static_cast < float >(iLM);
252 m_geom->Mapping(m_cellId - 1);
253 m_thetaId =
static_cast < float >(m_geom->GetThetaID());
254 m_phiId =
static_cast < float >(m_geom->GetPhiID());
257 if (m_isTrainingMode > 0) {
262 auto relatedParticlePairs = aECLCalDigit.getRelationsWith<
MCParticle>();
264 for (
unsigned int irel = 0; irel < relatedParticlePairs.size(); irel++) {
265 const auto particle = relatedParticlePairs.object(irel);
266 const double weight = relatedParticlePairs.weight(irel);
269 int motherindex = -1;
271 getEnteringMother(*particle, motherpdg, motherindex, pi0index);
272 addToSignalEnergy(motherpdg, motherindex, pi0index, weight);
275 B2DEBUG(175,
" -> digt energy: " << aECLCalDigit.getEnergy());
276 B2DEBUG(175,
"photon: " << m_signalEnergy[0][0] <<
" " << m_signalEnergy[0][1]);
277 B2DEBUG(175,
"pi0: " << m_signalEnergy[1][0] <<
" " << m_signalEnergy[1][1]);
278 B2DEBUG(175,
"electron: " << m_signalEnergy[2][0] <<
" " << m_signalEnergy[2][1]);
279 B2DEBUG(175,
"muon: " << m_signalEnergy[3][0] <<
" " << m_signalEnergy[3][1]);
280 B2DEBUG(175,
"neutral hadron: " << m_signalEnergy[4][0] <<
" " << m_signalEnergy[4][1]);
281 B2DEBUG(175,
"charged hadron: " << m_signalEnergy[5][0] <<
" " << m_signalEnergy[5][1]);
282 B2DEBUG(175,
"other: " << m_signalEnergy[6][0]);
286 getMax(maxtype, maxpos);
289 if (m_signalEnergy[maxtype][maxpos] >= m_truthFraction * aECLCalDigit.getEnergy()) {
291 m_targetindex = m_signalId[maxtype][maxpos];
292 m_targetpi0index = pi0index;
304 if (m_method ==
"cut") {
306 B2DEBUG(200,
"m_cutSlope: " << m_cutSlope <<
", m_nNeighbours10: " << m_nNeighbours10 <<
", m_cutOffset: " << m_cutOffset <<
307 ", m_maxNeighbourEnergy: " << m_maxNeighbourEnergy <<
", m_cutRatioCorrection: " << m_cutRatioCorrection <<
308 ", aECLCalDigit.getEnergy(): " << aECLCalDigit.getEnergy());
309 B2DEBUG(200,
"m_cutSlope * (m_nNeighbours10 - m_cutOffset): " << m_cutSlope * (m_nNeighbours10 - m_cutOffset));
310 B2DEBUG(200,
"(m_maxNeighbourEnergy - m_cutRatioCorrection) / (aECLCalDigit.getEnergy() - m_cutRatioCorrection)" <<
311 (m_maxNeighbourEnergy - m_cutRatioCorrection) / (aECLCalDigit.getEnergy() - m_cutRatioCorrection) <<
"\n");
313 if (m_cutSlope * (m_nNeighbours10 - m_cutOffset) >= (m_maxNeighbourEnergy - m_cutRatioCorrection) /
314 (aECLCalDigit.getEnergy() - m_cutRatioCorrection)) {
315 makeLocalMaximum(aCR, aECLCalDigit.getCellId(), iLM);
318 }
else if (m_method ==
"none") {
319 makeLocalMaximum(aCR, aECLCalDigit.getCellId(), iLM);
332 int highestEnergyCellId = -1;
333 double highestEnergy = 0.0;
337 if (aECLCalDigit.getEnergy() > highestEnergy) {
338 highestEnergyCellId = aECLCalDigit.getCellId();
339 highestEnergy = aECLCalDigit.getEnergy();
343 makeLocalMaximum(aCR, highestEnergyCellId, 0);
354 B2DEBUG(200,
"ECLLocalMaximumFinderModule::endRun()");
360 B2DEBUG(200,
"ECLLocalMaximumFinderModule::terminate()");
363 if (m_outfile and m_tree) {
371 if (m_neighbourMap)
delete m_neighbourMap;
378 const auto aLocalMaximum = m_eclLocalMaximums.appendNew();
380 B2DEBUG(175,
"ECLLocalMaximumFinderModule::makeLocalMaximum(): local maximum cellid: " << cellId);
383 aLocalMaximum->setLMId(lmId);
386 aLocalMaximum->setCellId(cellId);
397 m_targetpi0index = -1;
399 m_totalSignalEnergy = 0.;
401 for (
unsigned int i = 0; i < 10; ++i) {
402 for (
unsigned int j = 0; j < 5; ++j) {
403 m_signalEnergy[i][j] = 0.;
404 m_signalId[i][j] = -1;
413 for (
unsigned i = 0; i < c_nMaxNeighbours; ++i) {
414 m_energyRatioNeighbour[i] = 0.;
422 m_maxNeighbourEnergy = 0.0;
424 m_timeResolution = -9999.;
425 m_timeFitFailed = -9999.;
428 m_maxEnergyRatio = 0.0;
435 int index = particle.getArrayIndex();
438 while (!isEnteringECL(m_mcParticles[index]->getProductionVertex())) {
439 if (m_mcParticles[index]->getMother()) index = m_mcParticles[index]->getMother()->getArrayIndex();
444 if (m_mcParticles[index]->getPDG() == 22) {
445 if (m_mcParticles[index]->getMother()->getPDG() == 111) {
446 pi0index = m_mcParticles[index]->getMother()->getArrayIndex();
452 or (m_mcParticles[index]->getPDG() == 22 and abs(m_mcParticles[index]->getMother()->getPDG()) == 2212)) {
457 pdg = m_mcParticles[index]->getPDG();
459 pi0arrayindex = pi0index;
466 const double theta = vertex.Theta();
468 if (theta > 0.555015 and theta < 2.26369) {
469 double radius = vertex.Perp();
471 }
else if (theta <= 0.555015) {
473 }
else if (theta >= 2.26369) {
483 for (
int i = 0; i < 5; i++) {
484 if (m_signalId[type][i] ==
id)
return i;
485 if (m_signalId[type][i] == -1) {
486 m_signalId[type][i] = id;
501 for (
unsigned int i = 0; i < 10; ++i) {
502 for (
unsigned int j = 0; j < 5; ++j) {
503 if (m_signalEnergy[i][j] > maxe) {
504 maxe = m_signalEnergy[i][j];
520 if (motherpdg == 22) {
522 int idpos = getIdPosition(1, motherindex);
523 m_signalEnergy[1][idpos] += weight;
525 int idpos = getIdPosition(0, motherindex);
526 m_signalEnergy[0][idpos] += weight;
528 }
else if (abs(motherpdg) == 11) {
529 int idpos = getIdPosition(2, motherindex);
530 m_signalEnergy[2][idpos] += weight;
531 }
else if (abs(motherpdg) == 13) {
532 int idpos = getIdPosition(3, motherindex);
533 m_signalEnergy[3][idpos] += weight;
534 }
else if (abs(motherpdg) == 130 or abs(motherpdg) == 2112) {
535 int idpos = getIdPosition(4, motherindex);
536 m_signalEnergy[4][idpos] += weight;
537 }
else if (abs(motherpdg) == 211 or abs(motherpdg) == 321 or abs(motherpdg) == 2212) {
538 int idpos = getIdPosition(5, motherindex);
539 m_signalEnergy[5][idpos] += weight;
541 int idpos = getIdPosition(6, motherindex);
542 m_signalEnergy[6][idpos] += weight;