224 if (maxDistance < 0.0) {
225 maxDistance = 9999999.;
231 double backgroundLevel = 0.0;
233 backgroundLevel =
static_cast<double>(bkgdcount) /
static_cast<double>(
m_fullBkgdCount);
239 B2DEBUG(170,
"ECLCRSplitterModule::splitConnectedRegion: nLocalMaximums = " << nLocalMaximums);
248 if (nLocalMaximums == 1) {
254 aECLShower->addRelationTo(&aCR);
257 double weightSum = 0.0;
263 const int locmaxcellid = locmaxvector[0]->getCellId();
268 double highestEnergyTimeResolution = (
m_eclCalDigits[pos])->getTimeResolution();
280 std::vector<ECLCalDigit> digits;
281 std::vector<double> weights;
282 for (
auto& neighbourId : neighbourMap->
getNeighbours(highestEnergyID)) {
289 weights.push_back(1.0);
297 aECLShower->setTheta(showerposition.
Theta());
298 aECLShower->setPhi(showerposition.
Phi());
299 aECLShower->setR(showerposition.
Mag());
302 double showerEnergy = 0.0;
307 const unsigned int nOptimal =
static_cast<unsigned int>(nOptimalVec[0]);
308 aECLShower->setNominalNumberOfCrystalsForEnergy(
static_cast<double>(nOptimal));
312 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
313 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
314 aECLShower->setNOptimalEnergy(energyEstimation);
317 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
318 listCrystalPairs.resize(digits.size());
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());
328 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](
const auto & left,
const auto & right) {
329 return left.second > right.second;
331 std::vector< unsigned int> listCrystals;
333 for (
unsigned int i = 0; i < digits.size(); ++i) {
335 listCrystals.push_back(listCrystalPairs[i].first);
339 aECLShower->setNumberOfCrystalsForEnergy(
static_cast<double>(listCrystals.size()));
340 aECLShower->setListOfCrystalsForEnergy(listCrystals);
343 B2DEBUG(175,
"Shower Energy (1): " << showerEnergy);
346 showerEnergy = Belle2::ECL::computeEnergySum(digits, weights);
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);
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);
367 aECLShower->setShowerId(1);
369 aECLShower->setConnectedRegionId(aCR.
getCRId());
383 std::vector<std::pair<ECLLocalMaximum, double>> lm_energy_vector;
385 const int cellid = aLocalMaximum.getCellId();
388 lm_energy_vector.push_back(std::pair<ECLLocalMaximum, double>(aLocalMaximum, digitenergy));
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;
401 std::vector<ECLCalDigit> digits;
402 std::vector<double> weights;
403 std::map<int, B2Vector3D> centroidList;
404 std::map<int, double> centroidEnergyList;
405 std::map<int, B2Vector3D> allPoints;
406 std::map<int, std::vector < double > > weightMap;
407 std::vector < ECLCalDigit > digitVector;
410 std::map<int, B2Vector3D> localMaximumsPoints;
411 std::map<int, B2Vector3D> centroidPoints;
413 for (
auto& aLocalMaximum : lm_energy_vector) {
415 int cellid = aLocalMaximum.first.getCellId();
419 localMaximumsPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
420 centroidPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
424 bool iterateclusters =
true;
428 centroidList.clear();
429 centroidEnergyList.clear();
436 const int cellid = aCalDigit.getCellId();
439 allPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
440 digits.push_back(aCalDigit);
443 for (
const auto& digitpoint : allPoints) {
444 const int cellid = digitpoint.first;
455 double centroidShiftAverage = 0.0;
459 B2DEBUG(175,
"Iteration: #" << nIterations <<
" (of max. " <<
m_maxIterations <<
")");
461 centroidShiftAverage = 0.0;
462 if (nIterations == 0) {
463 centroidList.clear();
464 centroidEnergyList.clear();
469 for (
const auto& locmaxpoint : localMaximumsPoints) {
472 const int locmaxcellid = locmaxpoint.first;
478 if (nIterations == 0) {
481 centroidEnergyList[locmaxcellid] = 1.5 * locmaxenergy;
484 B2DEBUG(175,
"local maximum cellid: " << locmaxcellid);
488 for (
const auto& digitpoint : allPoints) {
491 const int digitcellid = digitpoint.first;
499 double distance = 0.0;
500 double distanceEnergySum = 0.0;
503 for (
const auto& centroidpoint : centroidPoints) {
506 const int centroidcellid = centroidpoint.first;
507 B2Vector3D centroidpos = centroidpoint.second;
509 double thisdistance = 0.;
512 if (nIterations == 0 and digitcellid == centroidcellid) {
515 B2Vector3D vectorDistance = ((centroidpos) - (digitpos));
516 thisdistance = vectorDistance.
Mag();
525 if ((locmaxcellid == centroidcellid) and (thisdistance < maxDistance)) {
526 distance = thisdistance;
531 if (thisdistance < maxDistance) {
533 distanceEnergySum += (thisenergy * expfactor);
539 if (distanceEnergySum > 0.0) {
541 weight = energy * expfactor / distanceEnergySum;
553 B2WARNING(
"ECLCRSplitterModule::splitConnectedRegion: Floating point glitch, weight for this digit " << weight <<
554 ", resetting it to 1.0.");
559 B2DEBUG(175,
" cellid: " << digitcellid <<
", energy: " << digitenergy <<
", weight: " << weight <<
", distance: " << distance);
560 weights.push_back(weight);
565 B2Vector3D oldCentroidPos = (centroidPoints.find(locmaxcellid))->second;
571 const double newEnergy = Belle2::ECL::computeEnergySum(digits, weights);
574 const B2Vector3D centroidShift = (oldCentroidPos - newCentroidPos);
577 centroidList[locmaxcellid] = newCentroidPos;
578 weightMap[locmaxcellid] = weights;
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");
586 centroidShiftAverage += centroidShift.
Mag();
589 B2DEBUG(175,
" old centroid: " << oldCentroidPos.
X() <<
" cm, " << oldCentroidPos.
Y() <<
" cm, " << oldCentroidPos.
Z() <<
591 B2DEBUG(175,
" new centroid: " << newCentroidPos.
X() <<
" cm, " << newCentroidPos.
Y() <<
" cm, " << newCentroidPos.
Z() <<
593 B2DEBUG(175,
" centroid shift: " << centroidShift.
Mag() <<
" cm");
598 centroidShiftAverage /=
static_cast<double>(nLocalMaximums);
600 B2DEBUG(175,
"--> average centroid shift: " << centroidShiftAverage <<
" cm (tolerance is " <<
m_shiftTolerance <<
" cm)");
603 for (
const auto& locmaxpoint : localMaximumsPoints) {
604 centroidPoints[locmaxpoint.first] = (centroidList.find(locmaxpoint.first))->second;
609 }
while (nIterations < m_maxIterations and centroidShiftAverage >
m_shiftTolerance);
613 std::vector<int> markfordeletion;
614 iterateclusters =
false;
615 for (
const auto& locmaxpoint : localMaximumsPoints) {
618 const int locmaxcellid = locmaxpoint.first;
621 B2DEBUG(175,
"locmaxcellid: " << locmaxcellid);
623 B2DEBUG(175,
"ok: ");
626 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
628 for (
unsigned int i = 0; i < digitVector.size(); ++i) {
631 const double weightInShower = myWeights[i];
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;
648 for (
const auto lmid : markfordeletion) {
649 localMaximumsPoints.erase(lmid);
650 centroidPoints.erase(lmid);
653 }
while (iterateclusters);
656 unsigned int iShower = 1;
657 for (
const auto& locmaxpoint : localMaximumsPoints) {
659 const int locmaxcellid = locmaxpoint.first;
675 std::vector<short int> neighbourlist = neighbourMap->
getNeighbours(locmaxcellid);
678 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
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;
688 for (
unsigned int i = 0; i < digitVector.size(); ++i) {
691 const double weight = myWeights[i];
701 if (std::find(neighbourlist.begin(), neighbourlist.end(), cellid) != neighbourlist.end()) {
705 newdigits.push_back(dig);
706 newweights.push_back(weight);
712 if (energy * weight > highestEnergy) {
713 highestEnergy = energy * weight;
714 highestEnergyTime = dig.
getTime();
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;
734 aECLShower->setTheta(showerposition->
Theta());
735 aECLShower->setPhi(showerposition->
Phi());
736 aECLShower->setR(showerposition->
Mag());
738 B2DEBUG(175,
"new theta: " << showerposition->
Theta());
739 B2DEBUG(175,
"new phi: " << showerposition->
Phi());
740 B2DEBUG(175,
"new R: " << showerposition->
Mag());
741 delete showerposition;
744 double showerEnergy = 0.0;
750 const unsigned int nOptimal =
static_cast<unsigned int>(nOptimalVec[0]);
751 aECLShower->setNominalNumberOfCrystalsForEnergy(
static_cast<double>(nOptimal));
755 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
756 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
757 aECLShower->setNOptimalEnergy(energyEstimation);
760 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
761 listCrystalPairs.resize(newdigits.size());
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());
771 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](
const auto & left,
const auto & right) {
772 return left.second > right.second;
775 std::vector< unsigned int> listCrystals;
777 for (
unsigned int i = 0; i < newdigits.size(); ++i) {
779 listCrystals.push_back(listCrystalPairs[i].first);
783 aECLShower->setNumberOfCrystalsForEnergy(
static_cast<double>(listCrystals.size()));
784 aECLShower->setListOfCrystalsForEnergy(listCrystals);
787 B2DEBUG(175,
"Shower Energy (2): " << showerEnergy);
790 showerEnergy = Belle2::ECL::computeEnergySum(newdigits, newweights);
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);
801 B2DEBUG(175,
"new energy: " << showerEnergy);
804 aECLShower->setShowerId(iShower);
807 aECLShower->setConnectedRegionId(aCR.
getCRId());
810 aECLShower->addRelationTo(&aCR);