Belle II Software development
ECLSplitterN1Module.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/eclSplitterN1/ECLSplitterN1Module.h>
11
12/* ECL header. */
13#include <ecl/dataobjects/ECLCalDigit.h>
14#include <ecl/dataobjects/ECLConnectedRegion.h>
15#include <ecl/dataobjects/ECLElementNumbers.h>
16#include <ecl/dataobjects/ECLLocalMaximum.h>
17#include <ecl/dataobjects/ECLShower.h>
18#include <ecl/dbobjects/ECLnOptimal.h>
19#include <ecl/geometry/ECLGeometryPar.h>
20#include <ecl/geometry/ECLNeighbours.h>
21#include <ecl/utility/Position.h>
22
23/* Basf2 headers. */
24#include <framework/geometry/B2Vector3.h>
25#include <framework/logging/Logger.h>
26#include <mdst/dataobjects/EventLevelClusteringInfo.h>
27
28/* C++ headers. */
29#include <algorithm>
30#include <string>
31#include <utility>
32
33// NAMESPACES
34using namespace Belle2;
35using namespace ECL;
36
37//-----------------------------------------------------------------
38// Register the Module(s)
39//-----------------------------------------------------------------
40REG_MODULE(ECLSplitterN1);
41REG_MODULE(ECLSplitterN1PureCsI);
42
43//-----------------------------------------------------------------
44// Implementation
45//-----------------------------------------------------------------
46
53{
54 // Set description.
55 setDescription("ECLSplitterN1Module: Baseline reconstruction splitter code for the n photon hypothesis.");
56 addParam("fullBkgdCount", m_fullBkgdCount,
57 "Number of background digits at full background (as provided by EventLevelClusteringInfo).",
58 182);
59
60 // Set module parameters.
61
62 // Splitter.
63 addParam("threshold", m_threshold, "Threshold energy after splitting.", 7.5 * Belle2::Unit::MeV);
64 addParam("expConstant", m_expConstant, "Constant a from exp(-a*dist/RM), typical: 1.5 to 3.5?", 2.5); // to be optimized!
65 addParam("maxIterations", m_maxIterations, "Maximum number of iterations for centroid shifts.", 100);
66 addParam("shiftTolerance", m_shiftTolerance, "Tolerance level for centroid shifts.", 1.0 * Belle2::Unit::mm);
67 addParam("minimumSharedEnergy", m_minimumSharedEnergy, "Minimum shared energy.", 25.0 * Belle2::Unit::keV);
68 addParam("maxSplits", m_maxSplits, "Maximum number of splits within one connected region.", 10);
69 addParam("cutDigitEnergyForEnergy", m_cutDigitEnergyForEnergy,
70 "Minimum digit energy to be included in the shower energy calculation. (NOT USED)", 0.5 * Belle2::Unit::MeV);
71 addParam("cutDigitTimeResidualForEnergy", m_cutDigitTimeResidualForEnergy,
72 "Maximum time residual to be included in the shower energy calculation. (NOT USED)", 5.0);
73 addParam("removeShiftedLMs", m_removeShiftedLMs, "Remove LMs if they shifted too much during energy sharing.", 0);
74 addParam("sharingDistanceMolierMultiplier", m_sharingDistanceMolierMultiplier,
75 "Maximum distance d to use when sharing energy d = m_sharingDistanceMolierMultiplier*c_molierRadius.", 50.);
76
77 // Neighbour definitions
78 addParam("useOptimalNumberOfDigitsForEnergy", m_useOptimalNumberOfDigitsForEnergy,
79 "Optimize the number of digits for energy calculations.", 1);
80
81 // Position.
82 addParam("positionMethod", m_positionMethod, "Position determination method.", std::string("lilo"));
83 addParam("liloParameterA", m_liloParameterA, "Position determination linear-log. parameter A.", 4.0);
84 addParam("liloParameterB", m_liloParameterB, "Position determination linear-log. parameter B.", 0.0);
85 addParam("liloParameterC", m_liloParameterC, "Position determination linear-log. parameter C.", 0.0);
86
87 // Set parallel processing flag.
89}
90
92{
93 // do not delete objects here, do it in terminate()!
94}
95
97{
98
99 // Geometry instance.
101
102 // Check and format user input.
103 m_liloParameters.resize(3);
107
108 // ECL dataobjects.
112 m_eclShowers.registerInDataStore(eclShowerArrayName());
113
114 // mDST dataobjects.
116
117 // Register relations (we probably dont need all, but keep them for now for debugging).
118 m_eclShowers.registerRelationTo(m_eclConnectedRegions);
119 m_eclShowers.registerRelationTo(m_eclCalDigits);
120 m_eclShowers.registerRelationTo(m_eclLocalMaximums);
121 m_eclLocalMaximums.registerRelationTo(m_eclCalDigits);
122 m_eclConnectedRegions.requireRelationTo(m_eclLocalMaximums);
123 m_eclConnectedRegions.requireRelationTo(m_eclCalDigits);
124
125 // Initialize neighbour maps (we will optimize the endcaps later, there is more than just a certain energy containment to be considered)
126 m_NeighbourMap9 = new ECLNeighbours("N", 1); // N: 3x3 = 9
127 m_NeighbourMap21 = new ECLNeighbours("NC", 2); // NC: 5x5 excluding corners = 21
128
129 // initialize the vector that gives the relation between cellid and store array position
132
133}
134
136{
137 //..Read in nOptimal crystal payload from database
138 if (m_eclNOptimal.hasChanged()) {
139
140 //..Vectors of energy boundaries for each region
141 std::vector<float> eBoundariesFwd = m_eclNOptimal->getUpperBoundariesFwd();
142 std::vector<float> eBoundariesBrl = m_eclNOptimal->getUpperBoundariesBrl();
143 std::vector<float> eBoundariesBwd = m_eclNOptimal->getUpperBoundariesBwd();
144
145 //..Adjust the size of the vector of boundaries to match the number in the payload
146 m_nEnergyBins = eBoundariesBrl.size();
147 B2INFO("ECLSplitterN1 beginRun: number of nOptimal payload energies = " << m_nEnergyBins);
148 m_eBoundaries.resize(m_nLeakReg, std::vector<float>(m_nEnergyBins, 0.));
149
150 //..Copy values to m_eBoundaries
151 for (int ie = 0; ie < m_nEnergyBins; ie++) {
152 m_eBoundaries[0][ie] = eBoundariesFwd[ie];
153 m_eBoundaries[1][ie] = eBoundariesBrl[ie];
154 m_eBoundaries[2][ie] = eBoundariesBwd[ie];
155 B2INFO(" upper boundaries for energy point " << ie << " " << m_eBoundaries[0][ie] << " " << m_eBoundaries[1][ie] << " " <<
156 m_eBoundaries[2][ie]);
157 }
158
159 //..Group number of each crystal
160 m_groupNumber = m_eclNOptimal->getGroupNumber();
161
162 //..2D histogram of nOptimal for each group and energy point
163 m_nOptimal2D = m_eclNOptimal->getNOptimal();
164 }
165}
166
168{
169 B2DEBUG(175, "ECLCRSplitterModule::event()");
170
171 // Fill a vector that can be used to map cellid -> store array position for eclCalDigits.
172 memset(&m_StoreArrPosition[0], -1, m_StoreArrPosition.size() * sizeof m_StoreArrPosition[0]);
173 for (int i = 0; i < m_eclCalDigits.getEntries(); i++) {
174 m_StoreArrPosition[m_eclCalDigits[i]->getCellId()] = i;
175 }
176
177 // Fill a vector that can be used to map cellid -> store array position for eclLocalMaximums.
178 memset(&m_StoreArrPositionLM[0], -1, m_StoreArrPositionLM.size() * sizeof m_StoreArrPositionLM[0]);
179 for (int i = 0; i < m_eclLocalMaximums.getEntries(); i++) {
180 m_StoreArrPositionLM[m_eclLocalMaximums[i]->getCellId()] = i;
181 }
182
183 // Loop over all connected regions
184 for (auto& aCR : m_eclConnectedRegions) {
185 // list theat will hold all cellids in this connected region
186 m_cellIdInCR.clear();
187
188 const unsigned int entries = (aCR.getRelationsWith<ECLCalDigit>(eclCalDigitArrayName())).size();
189
190 m_cellIdInCR.resize(entries);
191
192 // Fill all calDigits ids in this CR into a vector to make them 'find'-able.
193 int i = 0;
194 for (const auto& caldigit : aCR.getRelationsWith<ECLCalDigit>(eclCalDigitArrayName())) {
195 m_cellIdInCR[i] = caldigit.getCellId();
196 ++i;
197 }
198
199 // Split and reconstruct the showers in this connected regions.
201
202 } // end auto& aCR
203
204}
205
206
208{
209// if (m_tg2OptimalNumberOfDigitsForEnergy) delete m_tg2OptimalNumberOfDigitsForEnergy;
210}
211
212
218
220{
221
222 //calculate the maximum distance for later weighting
224 if (maxDistance < 0.0) {
225 maxDistance = 9999999.;
226 }
227
228
229 // Get the event background level
230 const int bkgdcount = m_eventLevelClusteringInfo->getNECLCalDigitsOutOfTime();
231 double backgroundLevel = 0.0; // from out of time digit counting
232 if (m_fullBkgdCount > 0) {
233 backgroundLevel = static_cast<double>(bkgdcount) / static_cast<double>(m_fullBkgdCount);
234 }
235
236 // Get the number of LMs in this CR
237 const int nLocalMaximums = aCR.getRelationsWith<ECLLocalMaximum>(eclLocalMaximumArrayName()).size();
238
239 B2DEBUG(170, "ECLCRSplitterModule::splitConnectedRegion: nLocalMaximums = " << nLocalMaximums);
240
241 // Three cases:
242 // 1) There is no local maximum (can only happen if the CR seed is not a LM itself).
243 // 2) There is exactly one local maximum.
244 // 3) There are more than one, typically two or three, local maxima and we have to share energy between them.
245 // If there are more than m_maxSplits local maxima, the m_maxSplits highest energy local maxima will be used.
246
247 // ---------------------------------------------------------------------
248 if (nLocalMaximums == 1) {
249
250 // Create a shower.
251 const auto aECLShower = m_eclShowers.appendNew();
252
253 // Add relation to the CR.
254 aECLShower->addRelationTo(&aCR);
255
256 // Find the highest energetic crystal in this CR or use the LM.
257 double weightSum = 0.0;
258
259 // Add relation to the LM.
261 aECLShower->addRelationTo(locmaxvector[0]);
262
263 const int locmaxcellid = locmaxvector[0]->getCellId();
264 const int pos = m_StoreArrPosition[locmaxcellid];
265 double highestEnergyID = (m_eclCalDigits[pos])->getCellId();
266 double highestEnergy = (m_eclCalDigits[pos])->getEnergy();
267 double highestEnergyTime = (m_eclCalDigits[pos])->getTime();
268 double highestEnergyTimeResolution = (m_eclCalDigits[pos])->getTimeResolution();
269
270 // Get a first estimation of the energy using 3x3 neighbours.
271 const double energyEstimation = estimateEnergy(highestEnergyID);
272
273 // Check if 21 would be better in the present background conditions:
274 ECLNeighbours* neighbourMap; // FIXME pointer needed?
275 int nNeighbours = getNeighbourMap(energyEstimation, backgroundLevel);
276 if (nNeighbours == 9 and !m_useOptimalNumberOfDigitsForEnergy) neighbourMap = m_NeighbourMap9;
277 else neighbourMap = m_NeighbourMap21;
278
279 // Add neighbours and weights for the shower.
280 std::vector<ECLCalDigit> digits;
281 std::vector<double> weights;
282 for (auto& neighbourId : neighbourMap->getNeighbours(highestEnergyID)) {
283 const auto it = std::find(m_cellIdInCR.begin(), m_cellIdInCR.end(),
284 neighbourId); // check if the neighbour is in the list for this CR
285 if (it == m_cellIdInCR.end()) continue; // not in this CR
286
287 const int neighbourpos = m_StoreArrPosition[neighbourId];
288 digits.push_back(*m_eclCalDigits[neighbourpos]); // list of digits for position reconstruction
289 weights.push_back(1.0); // list of weights (all 1 in this case for now)
290 weightSum += 1.0;
291
292 aECLShower->addRelationTo(m_eclCalDigits[neighbourpos], 1.0); // add digits to this shower, weight = 1
293 }
294
295 // Get position.
296 const B2Vector3D& showerposition = Belle2::ECL::computePositionLiLo(digits, weights, m_liloParameters);
297 aECLShower->setTheta(showerposition.Theta());
298 aECLShower->setPhi(showerposition.Phi());
299 aECLShower->setR(showerposition.Mag());
300
301 // Get Energy, if requested, set some weights to zero for energy calculation.
302 double showerEnergy = 0.0;
304
305 // Get the optimal number of neighbours for this crystal and energy
306 std::vector<int> nOptimalVec = getOptimalNumberOfDigits(highestEnergyID, energyEstimation);
307 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
308 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
309
310 // Store the indices used; will be needed later for energy corrections.
311 // Also store energy, used along with cellID to find leakage corrections.
312 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
313 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
314 aECLShower->setNOptimalEnergy(energyEstimation);
315
316 // Get the list of crystals used for the energy calculation
317 std::vector< std::pair<unsigned int, double>> listCrystalPairs; // cell id and weighted reconstructed energy
318 listCrystalPairs.resize(digits.size()); //resize to number of all crystals in cluster
319
320 std::vector < std::pair<double, double> > weighteddigits;
321 weighteddigits.resize(digits.size());
322 for (unsigned int i = 0; i < digits.size(); ++i) {
323 weighteddigits.at(i) = std::make_pair((digits.at(i)).getEnergy(), weights.at(i));
324 listCrystalPairs.at(i) = std::make_pair((digits.at(i)).getCellId(), weights.at(i) * (digits.at(i)).getEnergy());
325 }
326
327 // sort the listCrystals and keep the n highest in descending order
328 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
329 return left.second > right.second;
330 });
331 std::vector< unsigned int> listCrystals; //cell id
332
333 for (unsigned int i = 0; i < digits.size(); ++i) {
334 if (i < nOptimal) {
335 listCrystals.push_back(listCrystalPairs[i].first);
336 }
337 }
338
339 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
340 aECLShower->setListOfCrystalsForEnergy(listCrystals);
341
342 showerEnergy = getEnergySum(weighteddigits, nOptimal);
343 B2DEBUG(175, "Shower Energy (1): " << showerEnergy);
344
345 } else {
346 showerEnergy = Belle2::ECL::computeEnergySum(digits, weights);
347 }
348
349 aECLShower->setEnergy(showerEnergy);
350 aECLShower->setEnergyRaw(showerEnergy);
351 aECLShower->setEnergyHighestCrystal(highestEnergy);
352 aECLShower->setTime(highestEnergyTime);
353 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
354 aECLShower->setNumberOfCrystals(weightSum);
355 aECLShower->setCentralCellId(highestEnergyID);
356
357 B2DEBUG(175, "theta = " << showerposition.Theta());
358 B2DEBUG(175, "phi = " << showerposition.Phi());
359 B2DEBUG(175, "R = " << showerposition.Mag());
360 B2DEBUG(175, "energy = " << showerEnergy);
361 B2DEBUG(175, "time = " << highestEnergyTime);
362 B2DEBUG(175, "time resolution = " << highestEnergyTimeResolution);
363 B2DEBUG(175, "neighbours = " << nNeighbours);
364 B2DEBUG(175, "backgroundLevel = " << backgroundLevel);
365
366 // Fill shower Ids
367 aECLShower->setShowerId(1); // always one (only this single shower in the CR)
368 aECLShower->setHypothesisId(Belle2::ECLShower::c_nPhotons);
369 aECLShower->setConnectedRegionId(aCR.getCRId());
370
371 // Add relations of all CalDigits of the CR to the local maximum (here: all weights = 1).
372 const int posLM = m_StoreArrPositionLM[locmaxcellid];
373 for (const auto& aDigit : aCR.getRelationsWith<ECLCalDigit>()) {
374 const int posDigit = m_StoreArrPosition[aDigit.getCellId()];
375 m_eclLocalMaximums[posLM]->addRelationTo(m_eclCalDigits[posDigit], 1.0);
376 }
377
378 } // end case with one LM
379 else { // More than one LM, energy must be split. This algorithm is inspired by BaBar code.
380
381 // check if we have too many local maximums. if yes: limit to user set maximum
382 // create a vector with all local maximums and its crystal energies
383 std::vector<std::pair<ECLLocalMaximum, double>> lm_energy_vector;
384 for (auto& aLocalMaximum : aCR.getRelationsWith<ECLLocalMaximum>(eclLocalMaximumArrayName())) {
385 const int cellid = aLocalMaximum.getCellId();
386 const int pos = m_StoreArrPosition[cellid];
387 const double digitenergy = m_eclCalDigits[pos]->getEnergy();
388 lm_energy_vector.push_back(std::pair<ECLLocalMaximum, double>(aLocalMaximum, digitenergy));
389 };
390
391 // sort this vector in descending order and keep only up to m_maxSplits entries
392 if (lm_energy_vector.size() >= static_cast<size_t>(m_maxSplits)) {
393 std::sort(lm_energy_vector.begin(), lm_energy_vector.end(), [](const std::pair<ECLLocalMaximum, double>& x,
394 const std::pair<ECLLocalMaximum, double>& y) {
395 return x.second > y.second;
396 });
397
398 lm_energy_vector.resize(m_maxSplits);
399 }
400
401 std::vector<ECLCalDigit> digits;
402 std::vector<double> weights;
403 std::map<int, B2Vector3D> centroidList; // key = cellid, value = centroid position
404 std::map<int, double> centroidEnergyList; // key = cellid, value = centroid position
405 std::map<int, B2Vector3D> allPoints; // key = cellid, value = digit position
406 std::map<int, std::vector < double > > weightMap; // key = locmaxid, value = vector of weights
407 std::vector < ECLCalDigit > digitVector; // the order of weights in weightMap must be the same
408
409 // Fill the maxima positions in a map
410 std::map<int, B2Vector3D> localMaximumsPoints; // key = locmaxid, value = maximum position
411 std::map<int, B2Vector3D> centroidPoints; // key = locmaxid (as index), value = centroid position
412
413 for (auto& aLocalMaximum : lm_energy_vector) {
414
415 int cellid = aLocalMaximum.first.getCellId();
416
417 // Get the position of this crystal and fill it in two maps.
418 B2Vector3D vectorPosition = m_geom->GetCrystalPos(cellid - 1);
419 localMaximumsPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
420 centroidPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
421 }
422
423 // The following will be done iteratively. Empty clusters after splitting will be removed, and the procedure will be repeated.
424 bool iterateclusters = true;
425 do {
426 digits.clear();
427 weights.clear();
428 centroidList.clear();
429 centroidEnergyList.clear();
430 allPoints.clear();
431 weightMap.clear();
432 digitVector.clear();
433
434 // Fill all digits from this CR in a map
435 for (auto& aCalDigit : aCR.getRelationsWith<ECLCalDigit>(eclCalDigitArrayName())) {
436 const int cellid = aCalDigit.getCellId();
437 // get the position of this crystal and fill them in a map
438 B2Vector3D vectorPosition = m_geom->GetCrystalPos(cellid - 1);
439 allPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
440 digits.push_back(aCalDigit);
441 }
442
443 for (const auto& digitpoint : allPoints) {
444 const int cellid = digitpoint.first;
445 const int pos = m_StoreArrPosition[cellid];
446 digitVector.push_back(*m_eclCalDigits[pos]);
447 }
448
449
450 // -----------------------------------------------------------------------------------------
451 // The 'heart' of the splitter
452 // Start with each local maximum and set it to the first centroid position. Then iterate over all ECLCalDigits
453 // in this CR and calculate the weighted distances to the local maximum
454 int nIterations = 0;
455 double centroidShiftAverage = 0.0;
456 // double lastcentroidShiftAverage = 0.0;
457
458 do {
459 B2DEBUG(175, "Iteration: #" << nIterations << " (of max. " << m_maxIterations << ")");
460
461 centroidShiftAverage = 0.0;
462 if (nIterations == 0) {
463 centroidList.clear();
464 centroidEnergyList.clear();
465 weightMap.clear();
466 }
467
468 // Loop over all local maximums points, each one will become a shower!
469 for (const auto& locmaxpoint : localMaximumsPoints) {
470
471 // cell id of this local maximum
472 const int locmaxcellid = locmaxpoint.first;
473
474 // clear weights vector
475 weights.clear();
476
477 // if this is the first iteration the shower energy is not know, take the local maximum energy * 1.5.
478 if (nIterations == 0) {
479 const int pos = m_StoreArrPosition[locmaxcellid];
480 const double locmaxenergy = m_eclCalDigits[pos]->getEnergy();
481 centroidEnergyList[locmaxcellid] = 1.5 * locmaxenergy;
482 }
483
484 B2DEBUG(175, "local maximum cellid: " << locmaxcellid);
485
486 //-------------------------------------------------------------------
487 // Loop over all digits. They get a weight using the distance to the respective centroid.
488 for (const auto& digitpoint : allPoints) {
489
490 // cellid and position of this digit
491 const int digitcellid = digitpoint.first;
492 B2Vector3D digitpos = digitpoint.second;
493
494 const int pos = m_StoreArrPosition[digitcellid];
495 const double digitenergy = m_eclCalDigits[pos]->getEnergy();
496
497 double weight = 0.0;
498 double energy = 0.0;
499 double distance = 0.0;
500 double distanceEnergySum = 0.0;
501
502 // Loop over all centroids to get the normalization for the weight
503 for (const auto& centroidpoint : centroidPoints) {
504
505 // cell id and position of this centroid
506 const int centroidcellid = centroidpoint.first;
507 B2Vector3D centroidpos = centroidpoint.second;
508
509 double thisdistance = 0.;
510
511 // in the first iteration, this distance is really zero, avoid floating point problems
512 if (nIterations == 0 and digitcellid == centroidcellid) {
513 thisdistance = 0.0;
514 } else {
515 B2Vector3D vectorDistance = ((centroidpos) - (digitpos));
516 thisdistance = vectorDistance.Mag();
517 }
518
519 // energy of the centroid aka locmax
520 const int thispos = m_StoreArrPosition[centroidcellid];
521 const double thisenergy = m_eclCalDigits[thispos]->getEnergy();
522
523 // get distance of the current local maximum (its centroid position) to this caldigit
524 // Not the most efficienct way to get this information, but not worth the thinking yet:
525 if ((locmaxcellid == centroidcellid) and (thisdistance < maxDistance)) {
526 distance = thisdistance;
527 energy = thisenergy;
528 }
529
530 // Get the product of distance and energy
531 if (thisdistance < maxDistance) {
532 const double expfactor = exp(-m_expConstant * thisdistance / c_molierRadius);
533 distanceEnergySum += (thisenergy * expfactor);
534 }
535
536 } // end centroidPoints
537
538 // Calculate the weight for this digit for this local maximum.
539 if (distanceEnergySum > 0.0) {
540 const double expfactor = exp(-m_expConstant * distance / c_molierRadius);
541 weight = energy * expfactor / distanceEnergySum;
542 } else {
543 weight = 0.0;
544 }
545
546 // Check if the weighted energy is above threshold
547 if ((digitenergy * weight) < m_minimumSharedEnergy) {
548 weight = 0.0;
549 }
550
551 // Check for rounding problems larger than unity
552 if (weight > 1.0) {
553 B2WARNING("ECLCRSplitterModule::splitConnectedRegion: Floating point glitch, weight for this digit " << weight <<
554 ", resetting it to 1.0.");
555 weight = 1.0;
556 }
557
558 // Fill the weight for this digits and this local maximum.
559 B2DEBUG(175, " cellid: " << digitcellid << ", energy: " << digitenergy << ", weight: " << weight << ", distance: " << distance);
560 weights.push_back(weight);
561
562 } // end allPoints
563
564 // Get the old centroid position.
565 B2Vector3D oldCentroidPos = (centroidPoints.find(locmaxcellid))->second;
566
567 // Calculate the new centroid position.
568 B2Vector3D newCentroidPos = Belle2::ECL::computePositionLiLo(digits, weights, m_liloParameters);
569
570 // Calculate new energy
571 const double newEnergy = Belle2::ECL::computeEnergySum(digits, weights);
572
573 // Calculate the shift of the centroid position for this local maximum.
574 const B2Vector3D centroidShift = (oldCentroidPos - newCentroidPos);
575
576 // Save the new centroid position (but dont update yet!), also save the weights and energy.
577 centroidList[locmaxcellid] = newCentroidPos;
578 weightMap[locmaxcellid] = weights;
579
580 B2DEBUG(175, "--> inserting new energy: " << newEnergy << " for local maximum " << locmaxcellid);
581 centroidEnergyList[locmaxcellid] = newEnergy;
582 double showerenergy = (*centroidEnergyList.find(locmaxcellid)).second / Belle2::Unit::MeV;
583 B2DEBUG(175, "--> new energy = " << showerenergy << " MeV");
584
585 // Add this to the average centroid shift.
586 centroidShiftAverage += centroidShift.Mag();
587
588 // Debugging output
589 B2DEBUG(175, " old centroid: " << oldCentroidPos.X() << " cm, " << oldCentroidPos.Y() << " cm, " << oldCentroidPos.Z() <<
590 "cm");
591 B2DEBUG(175, " new centroid: " << newCentroidPos.X() << " cm, " << newCentroidPos.Y() << " cm, " << newCentroidPos.Z() <<
592 "cm");
593 B2DEBUG(175, " centroid shift: " << centroidShift.Mag() << " cm");
594
595 } // end localMaximumsPoints
596
597 // Get the average centroid shift.
598 centroidShiftAverage /= static_cast<double>(nLocalMaximums);
599 // lastcentroidShiftAverage = centroidShiftAverage;
600 B2DEBUG(175, "--> average centroid shift: " << centroidShiftAverage << " cm (tolerance is " << m_shiftTolerance << " cm)");
601
602 // Update centroid positions for the next round
603 for (const auto& locmaxpoint : localMaximumsPoints) {
604 centroidPoints[locmaxpoint.first] = (centroidList.find(locmaxpoint.first))->second;
605 }
606
607 ++nIterations;
608
609 } while (nIterations < m_maxIterations and centroidShiftAverage > m_shiftTolerance);
610 // DONE!
611
612 // check that local maxima are still local maxima
613 std::vector<int> markfordeletion;
614 iterateclusters = false;
615 for (const auto& locmaxpoint : localMaximumsPoints) {
616
617 // Get locmax cellid
618 const int locmaxcellid = locmaxpoint.first;
619 const int pos = m_StoreArrPosition[locmaxcellid];
620
621 B2DEBUG(175, "locmaxcellid: " << locmaxcellid);
622 const double LMEnergy = m_eclCalDigits[pos]->getEnergy();
623 B2DEBUG(175, "ok: ");
624
625 // Get the weight vector.
626 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
627
628 for (unsigned int i = 0; i < digitVector.size(); ++i) {
629
630 const ECLCalDigit dig = digitVector[i];
631 const double weightInShower = myWeights[i];
632 const int cellid = dig.getCellId();
633 const double energy = dig.getEnergy();
634
635 // two ways to fail:
636 // 1) another cell in this shower has more energy: energy*weight > LMEnergy and cellid != locmaxcellid
637 // 2) local maximum has cell has less than threshold energy left: energy*weight < m_threshold and cellid == locmaxcellid
638 if ((energy * weightInShower > LMEnergy and cellid != locmaxcellid and m_removeShiftedLMs > 0) or
639 (energy * weightInShower < m_threshold and cellid == locmaxcellid)) {
640 markfordeletion.push_back(locmaxcellid);
641 iterateclusters = true;
642 continue;
643 } // end check deletion for this digit
644 }// end digit loop
645 }
646
647 // delete LMs
648 for (const auto lmid : markfordeletion) {
649 localMaximumsPoints.erase(lmid);
650 centroidPoints.erase(lmid);
651 }
652
653 } while (iterateclusters);
654
655 // Create the ECLShower objects, one per LocalMaximumPoints
656 unsigned int iShower = 1;
657 for (const auto& locmaxpoint : localMaximumsPoints) {
658
659 const int locmaxcellid = locmaxpoint.first;
660 const int posLM = m_StoreArrPositionLM[locmaxcellid];
661
662 // Create a shower
663 const auto aECLShower = m_eclShowers.appendNew();
664
665 // Use the same method for the estimate (3x3).
666 const double energyEstimation = estimateEnergy(locmaxcellid);
667
668 // Get the neighbour list.
669 ECLNeighbours* neighbourMap; // FIXME need pointer?
670 int nNeighbours = getNeighbourMap(energyEstimation, backgroundLevel);
671 if (nNeighbours == 9 and !m_useOptimalNumberOfDigitsForEnergy) neighbourMap = m_NeighbourMap9;
672 else neighbourMap = m_NeighbourMap21;
673
674 // Get the neighbour list.
675 std::vector<short int> neighbourlist = neighbourMap->getNeighbours(locmaxcellid);
676
677 // Get the weight vector.
678 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
679
680 // Loop over all digits.
681 std::vector<ECLCalDigit> newdigits;
682 std::vector<double> newweights;
683 double highestEnergy = 0.;
684 double highestEnergyTime = 0.;
685 double highestEnergyTimeResolution = 0.;
686 double weightSum = 0.0;
687
688 for (unsigned int i = 0; i < digitVector.size(); ++i) {
689
690 const ECLCalDigit dig = digitVector[i];
691 const double weight = myWeights[i];
692
693 const int cellid = dig.getCellId();
694 const int pos = m_StoreArrPosition[cellid];
695
696 // Add weighted relations of all CalDigits to the local maximum.
697 m_eclLocalMaximums[posLM]->addRelationTo(m_eclCalDigits[pos], weight);
698
699 // Positive weight and in allowed neighbour list?
700 if (weight > 0.0) {
701 if (std::find(neighbourlist.begin(), neighbourlist.end(), cellid) != neighbourlist.end()) {
702
703 aECLShower->addRelationTo(m_eclCalDigits[pos], weight);
704
705 newdigits.push_back(dig);
706 newweights.push_back(weight);
707
708 weightSum += weight;
709
710 const double energy = dig.getEnergy();
711
712 if (energy * weight > highestEnergy) {
713 highestEnergy = energy * weight;
714 highestEnergyTime = dig.getTime();
715 highestEnergyTimeResolution = dig.getTimeResolution();
716 }
717 }
718 }
719 }
720
721 // Old position:
722 B2Vector3D* oldshowerposition = new B2Vector3D((centroidList.find(locmaxcellid))->second);
723
724 B2DEBUG(175, "old theta: " << oldshowerposition->Theta());
725 B2DEBUG(175, "old phi: " << oldshowerposition->Phi());
726 B2DEBUG(175, "old R: " << oldshowerposition->Mag());
727 B2DEBUG(175, "old energy: " << energyEstimation);
728 delete oldshowerposition;
729
730 // New position (with reduced number of neighbours)
731 // There are some cases where high backgrounds fake local maxima and the split centroid position is far
732 // away from the original LM cell... this will throw a (non fatal) error, and create a cluster with zero energy now).
733 B2Vector3D* showerposition = new B2Vector3D(Belle2::ECL::computePositionLiLo(newdigits, newweights, m_liloParameters));
734 aECLShower->setTheta(showerposition->Theta());
735 aECLShower->setPhi(showerposition->Phi());
736 aECLShower->setR(showerposition->Mag());
737
738 B2DEBUG(175, "new theta: " << showerposition->Theta());
739 B2DEBUG(175, "new phi: " << showerposition->Phi());
740 B2DEBUG(175, "new R: " << showerposition->Mag());
741 delete showerposition;
742
743 // Get Energy, if requested, set weights to zero for energy calculation.
744 double showerEnergy = 0.0;
746
747
748 // Get the optimal number of neighbours for this crystal and energy
749 std::vector<int> nOptimalVec = getOptimalNumberOfDigits(locmaxcellid, energyEstimation);
750 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
751 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
752
753 // Store the indices used; will be needed later for energy corrections.
754 // Also store energy, used along with cellID to find leakage corrections.
755 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
756 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
757 aECLShower->setNOptimalEnergy(energyEstimation);
758
759 // Get the list of crystals used for the energy calculation
760 std::vector< std::pair<unsigned int, double>> listCrystalPairs; // cell id and weighted reconstructed energy
761 listCrystalPairs.resize(newdigits.size()); //resize to number of all crystals in cluster
762
763 std::vector < std::pair<double, double> > weighteddigits;
764 weighteddigits.resize(newdigits.size());
765 for (unsigned int i = 0; i < newdigits.size(); ++i) {
766 weighteddigits.at(i) = std::make_pair((newdigits.at(i)).getEnergy(), newweights.at(i));
767 listCrystalPairs.at(i) = std::make_pair((newdigits.at(i)).getCellId(), newweights.at(i) * (newdigits.at(i)).getEnergy());
768 }
769
770 // sort the listCrystals and keep the n highest in descending order
771 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
772 return left.second > right.second;
773 });
774
775 std::vector< unsigned int> listCrystals; //cell id
776
777 for (unsigned int i = 0; i < newdigits.size(); ++i) {
778 if (i < nOptimal) {
779 listCrystals.push_back(listCrystalPairs[i].first);
780 }
781 }
782
783 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
784 aECLShower->setListOfCrystalsForEnergy(listCrystals);
785
786 showerEnergy = getEnergySum(weighteddigits, nOptimal);
787 B2DEBUG(175, "Shower Energy (2): " << showerEnergy);
788
789 } else {
790 showerEnergy = Belle2::ECL::computeEnergySum(newdigits, newweights);
791 }
792
793 aECLShower->setEnergy(showerEnergy);
794 aECLShower->setEnergyRaw(showerEnergy);
795 aECLShower->setEnergyHighestCrystal(highestEnergy);
796 aECLShower->setTime(highestEnergyTime);
797 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
798 aECLShower->setNumberOfCrystals(weightSum);
799 aECLShower->setCentralCellId(locmaxcellid);
800
801 B2DEBUG(175, "new energy: " << showerEnergy);
802
803 // Get unique ID
804 aECLShower->setShowerId(iShower);
805 ++iShower;
806 aECLShower->setHypothesisId(Belle2::ECLShower::c_nPhotons);
807 aECLShower->setConnectedRegionId(aCR.getCRId());
808
809 // Add relation to the CR.
810 aECLShower->addRelationTo(&aCR);
811
812 // Add relation to the LM.
813 aECLShower->addRelationTo(m_eclLocalMaximums[posLM]);
814 }
815 }
816}
817
818int ECLSplitterN1Module::getNeighbourMap(const double energy, const double background)
819{
820 if (background <= 0.1) return 21;
821 else {
822 if (energy > 0.06 + 0.4 * background) return 21; // based on preliminary study TF, valid in barrel only (TF).
823 else return 9;
824 }
825}
826
827std::vector<int> ECLSplitterN1Module::getOptimalNumberOfDigits(const int cellid, const double energy)
828{
829
830 //..nOptimal depends on the energy bin, which in turn depends on ECL region
831 int iRegion = 1; // barrel
832 if (ECLElementNumbers::isForward(cellid)) {iRegion = 0;}
833 if (ECLElementNumbers::isBackward(cellid)) {iRegion = 2;}
834
835 //..Find the energy bin
836 int iEnergy = 0;
837 while (energy > m_eBoundaries[iRegion][iEnergy] and iEnergy < m_nEnergyBins - 1) {iEnergy++;}
838
839 //..Group number just depends on the cellID
840 int iGroup = m_groupNumber[cellid - 1];
841
842 //..Optimal number of crystals from energy and group numbers
843 int nOptimalNeighbours = (int)(0.5 + m_nOptimal2D.GetBinContent(iGroup + 1, iEnergy + 1));
844
845 //..Store these in a vector to return
846 std::vector<int> nOptimalVector;
847 nOptimalVector.push_back(nOptimalNeighbours);
848 nOptimalVector.push_back(iGroup);
849 nOptimalVector.push_back(iEnergy);
850
851 B2DEBUG(175, "ECLSplitterN1Module::getOptimalNumberOfDigits: cellID: " << cellid << " energy: " << energy << " iG: " << iGroup <<
852 " iE: " << iEnergy << " nOpt: " << nOptimalNeighbours);
853
854 return nOptimalVector;
855
856}
857
858double ECLSplitterN1Module::getEnergySum(std::vector < std::pair<double, double> >& weighteddigits, const unsigned int n)
859{
860
861 double energysum = 0.;
862
863 std::sort(weighteddigits.begin(), weighteddigits.end(), [](const auto & left, const auto & right) {
864 return left.first * left.second > right.first * right.second;
865 });
866
867
868 unsigned int min = n;
869 if (weighteddigits.size() < n) min = weighteddigits.size();
870
871 for (unsigned int i = 0; i < min; ++i) {
872 B2DEBUG(175, "getEnergySum: " << weighteddigits.at(i).first << " " << weighteddigits.at(i).second);
873 energysum += (weighteddigits.at(i).first * weighteddigits.at(i).second);
874 }
875 B2DEBUG(175, "getEnergySum: energysum=" << energysum);
876
877 return energysum;
878}
879
880
881double ECLSplitterN1Module::estimateEnergy(const int centerid)
882{
883
884 double energyEstimation = 0.0;
885
886 for (auto& neighbourId : m_NeighbourMap9->getNeighbours(centerid)) {
887
888 // Check if this neighbour is in this CR
889 const auto it = std::find(m_cellIdInCR.begin(), m_cellIdInCR.end(),
890 neighbourId); // check if the neighbour is in the list for this CR
891 if (it == m_cellIdInCR.end()) continue; // not in this CR
892
893 const int pos = m_StoreArrPosition[neighbourId];
894 const double energyNeighbour = m_eclCalDigits[pos]->getEnergy();
895
896 energyEstimation += energyNeighbour;
897 }
898
899 return energyEstimation;
900}
DataType Phi() const
The azimuth angle.
Definition B2Vector3.h:151
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition B2Vector3.h:435
DataType Theta() const
The polar angle.
Definition B2Vector3.h:153
DataType X() const
access variable X (= .at(0) without boundary check)
Definition B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition B2Vector3.h:159
Class to store calibrated ECLDigits: ECLCalDigits.
Definition ECLCalDigit.h:23
int getCellId() const
Get Cell ID.
double getEnergy() const
Get Calibrated Energy.
double getTimeResolution() const
Get Calibrated Time Resolution.
double getTime() const
Get Calibrated Time.
Class to store connected regions (CRs)
int getCRId() const
Get CR identifier.
Class to store local maxima (LM)
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLShower.h:42
virtual const char * eclShowerArrayName() const
Default name ECLShowers.
ECL::ECLNeighbours * m_NeighbourMap9
Neighbour maps.
virtual const char * eventLevelClusteringInfoName() const
Name to be used for default option: EventLevelClusteringInfo.
double m_liloParameterB
lin-log parameter B
int m_maxIterations
Maximum number of iterations.
StoreArray< ECLShower > m_eclShowers
Store array: ECLShower.
ECL::ECLNeighbours * m_NeighbourMap21
5x5 neighbours excluding corners = 21
double m_sharingDistanceMolierMultiplier
Maximum distance d to use when sharing energy d = m_sharingDistanceMolierMultiplier*c_molierRadius.
StoreArray< ECLConnectedRegion > m_eclConnectedRegions
Store array: ECLConnectedRegion.
int m_removeShiftedLMs
Remove LMs if the splitter shifted the centroid position too much.
double m_minimumSharedEnergy
Minimum shared energy.
const int m_nLeakReg
3 ECL regions: 0 = forward, 1 = barrel, 2 = backward
std::string m_positionMethod
Position calculation: lilo or linear.
void splitConnectedRegion(ECLConnectedRegion &aCR)
Split connected region into showers.
virtual void initialize() override
Initialize.
double m_cutDigitTimeResidualForEnergy
Maximum time residual to be included in the shower energy calculation.
StoreArray< ECLLocalMaximum > m_eclLocalMaximums
Store array: ECLLocalMaximum.
virtual void event() override
Event.
virtual const char * eclLocalMaximumArrayName() const
Default name ECLLocalMaximums.
double m_threshold
Local maximum threshold after splitting.
double estimateEnergy(const int centerid)
Estimate energy using 3x3 around central crystal.
virtual void endRun() override
End run.
virtual void terminate() override
Terminate.
std::vector< int > getOptimalNumberOfDigits(const int cellid, const double energy)
Get optimal number of digits (out of 21) based on first energy estimation and background level per ev...
ECL::ECLGeometryPar * m_geom
Geometry.
std::vector< int > m_cellIdInCR
list with all cellid of this connected region
double m_shiftTolerance
Tolerance level for centroid shifts.
DBObjPtr< ECLnOptimal > m_eclNOptimal
nOptimal payload
virtual void beginRun() override
Begin run.
int m_fullBkgdCount
Number of expected background digits at full background, FIXME: move to database.
int m_nEnergyBins
number of energies bins in nOptimal payload
const double c_molierRadius
Constant RM (Molier Radius) from exp(-a*dist/RM), http://pdg.lbl.gov/2009/AtomicNuclearProperties/HTM...
virtual const char * eclConnectedRegionArrayName() const
Default name ECLConnectedRegions.
int m_maxSplits
Maximum number of splits.
int m_useOptimalNumberOfDigitsForEnergy
Optimize the number of neighbours for energy calculations.
int getNeighbourMap(const double energy, const double background)
Get number of neighbours based on first energy estimation and background level per event.
double m_cutDigitEnergyForEnergy
Minimum digit energy to be included in the shower energy calculation.
virtual const char * eclCalDigitArrayName() const
Default name ECLCalDigits.
std::vector< double > m_liloParameters
lin-log parameters A, B, and C
std::vector< int > m_StoreArrPosition
vector (ECLElementNumbers::c_NCrystals + 1 entries) with cell id to store array positions
std::vector< int > m_StoreArrPositionLM
vector (ECLElementNumbers::c_NCrystals + 1 entries) with cell id to store array positions for LM
double m_liloParameterA
lin-log parameter A
StoreObjPtr< EventLevelClusteringInfo > m_eventLevelClusteringInfo
Store object pointer: EventLevelClusteringInfo.
std::vector< std::vector< float > > m_eBoundaries
energy boundaries each region
TH2F m_nOptimal2D
2D hist of nOptimal for Ebin vs groupID
std::vector< int > m_groupNumber
group number for each crystal
double getEnergySum(std::vector< std::pair< double, double > > &weighteddigits, const unsigned int n)
Get energy sum for weighted entries.
StoreArray< ECLCalDigit > m_eclCalDigits
Store array: ECLCalDigit.
double m_expConstant
Constant a from exp(-a*dist/RM), 1.5 to 2.5.
double m_liloParameterC
lin-log parameter C
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Class to get the neighbours for a given cell id.
const std::vector< short int > & getNeighbours(short int cid) const
Return the neighbours for a given cell ID.
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 for type safe access to objects that are referred to in relations.
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< T > getRelationsWith(const std::string &name="", const std::string &namedRelation="") const
Get the relations between this object and another store array.
static const double mm
[millimeters]
Definition Unit.h:70
static const double keV
[kiloelectronvolt]
Definition Unit.h:113
static const double MeV
[megaelectronvolt]
Definition Unit.h:114
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
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
bool isForward(int cellId)
Check whether the crystal is in forward 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.