Split connected region into showers.
215{
216
217
219 double backgroundLevel = 0.0;
221 backgroundLevel =
static_cast<double>(bkgdcount) /
static_cast<double>(
m_fullBkgdCount);
222 }
223
224
226
227 B2DEBUG(170, "ECLCRSplitterModule::splitConnectedRegion: nLocalMaximums = " << nLocalMaximums);
228
229
230
231
232
233
234
235
236 if (nLocalMaximums == 1) {
237
238
240
241
242 aECLShower->addRelationTo(&aCR);
243
244
245 double weightSum = 0.0;
246
247
250
251 const int locmaxcellid = locmaxvector[0]->getCellId();
256 double highestEnergyTimeResolution = (
m_eclCalDigits[pos])->getTimeResolution();
257
258
260
261
266
267
268 std::vector<ECLCalDigit> digits;
269 std::vector<double> weights;
270 for (
auto& neighbourId : neighbourMap->
getNeighbours(highestEnergyID)) {
272 neighbourId);
274
277 weights.push_back(1.0);
278 weightSum += 1.0;
279
281 }
282
283
285 aECLShower->setTheta(showerposition.
Theta());
286 aECLShower->setPhi(showerposition.
Phi());
287 aECLShower->setR(showerposition.
Mag());
288
289
290 double showerEnergy = 0.0;
292
293
295 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
296 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
297
298
299
300 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
301 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
302 aECLShower->setNOptimalEnergy(energyEstimation);
303
304
305 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
306 listCrystalPairs.resize(digits.size());
307
308 std::vector < std::pair<double, double> > weighteddigits;
309 weighteddigits.resize(digits.size());
310 for (unsigned int i = 0; i < digits.size(); ++i) {
311 weighteddigits.at(i) = std::make_pair((digits.at(i)).getEnergy(), weights.at(i));
312 listCrystalPairs.at(i) = std::make_pair((digits.at(i)).getCellId(), weights.at(i) * (digits.at(i)).getEnergy());
313 }
314
315
316 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
317 return left.second > right.second;
318 });
319 std::vector< unsigned int> listCrystals;
320
321 for (unsigned int i = 0; i < digits.size(); ++i) {
322 if (i < nOptimal) {
323 listCrystals.push_back(listCrystalPairs[i].first);
324 }
325 }
326
327 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
328 aECLShower->setListOfCrystalsForEnergy(listCrystals);
329
331 B2DEBUG(175, "Shower Energy (1): " << showerEnergy);
332
333 } else {
334 showerEnergy = Belle2::ECL::computeEnergySum(digits, weights);
335 }
336
337 aECLShower->setEnergy(showerEnergy);
338 aECLShower->setEnergyRaw(showerEnergy);
339 aECLShower->setEnergyHighestCrystal(highestEnergy);
340 aECLShower->setTime(highestEnergyTime);
341 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
342 aECLShower->setNumberOfCrystals(weightSum);
343 aECLShower->setCentralCellId(highestEnergyID);
344
345 B2DEBUG(175,
"theta = " << showerposition.
Theta());
346 B2DEBUG(175,
"phi = " << showerposition.
Phi());
347 B2DEBUG(175,
"R = " << showerposition.
Mag());
348 B2DEBUG(175, "energy = " << showerEnergy);
349 B2DEBUG(175, "time = " << highestEnergyTime);
350 B2DEBUG(175, "time resolution = " << highestEnergyTimeResolution);
351 B2DEBUG(175, "neighbours = " << nNeighbours);
352 B2DEBUG(175, "backgroundLevel = " << backgroundLevel);
353
354
355 aECLShower->setShowerId(1);
357 aECLShower->setConnectedRegionId(aCR.
getCRId());
358
359
364 }
365
366 }
367 else {
368
369
370
371 std::vector<std::pair<ECLLocalMaximum, double>> lm_energy_vector;
373 const int cellid = aLocalMaximum.getCellId();
376 lm_energy_vector.push_back(std::pair<ECLLocalMaximum, double>(aLocalMaximum, digitenergy));
377 };
378
379
380 if (lm_energy_vector.size() >=
static_cast<size_t>(
m_maxSplits)) {
381 std::sort(lm_energy_vector.begin(), lm_energy_vector.end(), [](const std::pair<ECLLocalMaximum, double>& x,
382 const std::pair<ECLLocalMaximum, double>& y) {
383 return x.second > y.second;
384 });
385
387 }
388
389 std::vector<ECLCalDigit> digits;
390 std::vector<double> weights;
391 std::map<int, B2Vector3D> centroidList;
392 std::map<int, double> centroidEnergyList;
393 std::map<int, B2Vector3D> allPoints;
394 std::map<int, std::vector < double > > weightMap;
395 std::vector < ECLCalDigit > digitVector;
396
397
398 std::map<int, B2Vector3D> localMaximumsPoints;
399 std::map<int, B2Vector3D> centroidPoints;
400
401 for (auto& aLocalMaximum : lm_energy_vector) {
402
403 int cellid = aLocalMaximum.first.getCellId();
404
405
407 localMaximumsPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
408 centroidPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
409 }
410
411
412 bool iterateclusters = true;
413 do {
414 digits.clear();
415 weights.clear();
416 centroidList.clear();
417 centroidEnergyList.clear();
418 allPoints.clear();
419 weightMap.clear();
420 digitVector.clear();
421
422
424 const int cellid = aCalDigit.getCellId();
425
427 allPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
428 digits.push_back(aCalDigit);
429 }
430
431 for (const auto& digitpoint : allPoints) {
432 const int cellid = digitpoint.first;
435 }
436
437
438
439
440
441
442 int nIterations = 0;
443 double centroidShiftAverage = 0.0;
444
445
446 do {
447 B2DEBUG(175,
"Iteration: #" << nIterations <<
" (of max. " <<
m_maxIterations <<
")");
448
449 centroidShiftAverage = 0.0;
450 if (nIterations == 0) {
451 centroidList.clear();
452 centroidEnergyList.clear();
453 weightMap.clear();
454 }
455
456
457 for (const auto& locmaxpoint : localMaximumsPoints) {
458
459
460 const int locmaxcellid = locmaxpoint.first;
461
462
463 weights.clear();
464
465
466 if (nIterations == 0) {
469 centroidEnergyList[locmaxcellid] = 1.5 * locmaxenergy;
470 }
471
472 B2DEBUG(175, "local maximum cellid: " << locmaxcellid);
473
474
475
476 for (const auto& digitpoint : allPoints) {
477
478
479 const int digitcellid = digitpoint.first;
481
484
485 double weight = 0.0;
486 double energy = 0.0;
487 double distance = 0.0;
488 double distanceEnergySum = 0.0;
489
490
491 for (const auto& centroidpoint : centroidPoints) {
492
493
494 const int centroidcellid = centroidpoint.first;
495 B2Vector3D centroidpos = centroidpoint.second;
496
497 double thisdistance = 0.;
498
499
500 if (nIterations == 0 and digitcellid == centroidcellid) {
501 thisdistance = 0.0;
502 } else {
503 B2Vector3D vectorDistance = ((centroidpos) - (digitpos));
504 thisdistance = vectorDistance.
Mag();
505 }
506
507
510
511
512 if (locmaxcellid == centroidcellid) {
513 distance = thisdistance;
514 energy = thisenergy;
515 }
516
517
519 distanceEnergySum += (thisenergy * expfactor);
520
521 }
522
523
524 if (distanceEnergySum > 0.0) {
526 weight = energy * expfactor / distanceEnergySum;
527 } else {
528 weight = 0.0;
529 }
530
531
533 weight = 0.0;
534 }
535
536
537 if (weight > 1.0) {
538 B2WARNING("ECLCRSplitterModule::splitConnectedRegion: Floating point glitch, weight for this digit " << weight <<
539 ", resetting it to 1.0.");
540 weight = 1.0;
541 }
542
543
544 B2DEBUG(175, " cellid: " << digitcellid << ", energy: " << digitenergy << ", weight: " << weight << ", distance: " << distance);
545 weights.push_back(weight);
546
547 }
548
549
550 B2Vector3D oldCentroidPos = (centroidPoints.find(locmaxcellid))->second;
551
552
554
555
556 const double newEnergy = Belle2::ECL::computeEnergySum(digits, weights);
557
558
559 const B2Vector3D centroidShift = (oldCentroidPos - newCentroidPos);
560
561
562 centroidList[locmaxcellid] = newCentroidPos;
563 weightMap[locmaxcellid] = weights;
564
565 B2DEBUG(175, "--> inserting new energy: " << newEnergy << " for local maximum " << locmaxcellid);
566 centroidEnergyList[locmaxcellid] = newEnergy;
567 double showerenergy = (*centroidEnergyList.find(locmaxcellid)).second /
Belle2::Unit::MeV;
568 B2DEBUG(175, "--> new energy = " << showerenergy << " MeV");
569
570
571 centroidShiftAverage += centroidShift.
Mag();
572
573
574 B2DEBUG(175,
" old centroid: " << oldCentroidPos.
X() <<
" cm, " << oldCentroidPos.
Y() <<
" cm, " << oldCentroidPos.
Z() <<
575 "cm");
576 B2DEBUG(175,
" new centroid: " << newCentroidPos.
X() <<
" cm, " << newCentroidPos.
Y() <<
" cm, " << newCentroidPos.
Z() <<
577 "cm");
578 B2DEBUG(175,
" centroid shift: " << centroidShift.
Mag() <<
" cm");
579
580 }
581
582
583 centroidShiftAverage /= static_cast<double>(nLocalMaximums);
584
585 B2DEBUG(175,
"--> average centroid shift: " << centroidShiftAverage <<
" cm (tolerance is " <<
m_shiftTolerance <<
" cm)");
586
587
588 for (const auto& locmaxpoint : localMaximumsPoints) {
589 centroidPoints[locmaxpoint.first] = (centroidList.find(locmaxpoint.first))->second;
590 }
591
592 ++nIterations;
593
594 }
while (nIterations < m_maxIterations and centroidShiftAverage >
m_shiftTolerance);
595
596
597
598 std::vector<int> markfordeletion;
599 iterateclusters = false;
600 for (const auto& locmaxpoint : localMaximumsPoints) {
601
602
603 const int locmaxcellid = locmaxpoint.first;
605
606 B2DEBUG(175, "locmaxcellid: " << locmaxcellid);
608 B2DEBUG(175, "ok: ");
609
610
611 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
612
613 for (unsigned int i = 0; i < digitVector.size(); ++i) {
614
616 const double weight = myWeights[i];
619
620
621
622
623 if ((energy * weight > lmenergy and cellid != locmaxcellid) or (energy * weight <
m_threshold and cellid == locmaxcellid)) {
624 markfordeletion.push_back(locmaxcellid);
625 iterateclusters = true;
626 continue;
627 }
628 }
629 }
630
631
632 for (const auto lmid : markfordeletion) {
633 localMaximumsPoints.erase(lmid);
634 centroidPoints.erase(lmid);
635 }
636
637 } while (iterateclusters);
638
639
640 unsigned int iShower = 1;
641 for (const auto& locmaxpoint : localMaximumsPoints) {
642
643 const int locmaxcellid = locmaxpoint.first;
645
646
648
649
651
652
657
658
659 std::vector<short int> neighbourlist = neighbourMap->
getNeighbours(locmaxcellid);
660
661
662 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
663
664
665 std::vector<ECLCalDigit> newdigits;
666 std::vector<double> newweights;
667 double highestEnergy = 0.;
668 double highestEnergyTime = 0.;
669 double highestEnergyTimeResolution = 0.;
670 double weightSum = 0.0;
671
672 for (unsigned int i = 0; i < digitVector.size(); ++i) {
673
675 const double weight = myWeights[i];
676
679
680
682
683
684 if (weight > 0.0) {
685 if (std::find(neighbourlist.begin(), neighbourlist.end(), cellid) != neighbourlist.end()) {
686
688
689 newdigits.push_back(dig);
690 newweights.push_back(weight);
691
692 weightSum += weight;
693
695
696 if (energy * weight > highestEnergy) {
697 highestEnergy = energy * weight;
698 highestEnergyTime = dig.
getTime();
700 }
701 }
702 }
703 }
704
705
707
708 B2DEBUG(175,
"old theta: " << oldshowerposition->
Theta());
709 B2DEBUG(175,
"old phi: " << oldshowerposition->
Phi());
710 B2DEBUG(175,
"old R: " << oldshowerposition->
Mag());
711 B2DEBUG(175, "old energy: " << energyEstimation);
712 delete oldshowerposition;
713
714
715
716
718 aECLShower->setTheta(showerposition->
Theta());
719 aECLShower->setPhi(showerposition->
Phi());
720 aECLShower->setR(showerposition->
Mag());
721
722 B2DEBUG(175,
"new theta: " << showerposition->
Theta());
723 B2DEBUG(175,
"new phi: " << showerposition->
Phi());
724 B2DEBUG(175,
"new R: " << showerposition->
Mag());
725 delete showerposition;
726
727
728 double showerEnergy = 0.0;
730
731
732
734 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
735 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
736
737
738
739 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
740 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
741 aECLShower->setNOptimalEnergy(energyEstimation);
742
743
744 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
745 listCrystalPairs.resize(newdigits.size());
746
747 std::vector < std::pair<double, double> > weighteddigits;
748 weighteddigits.resize(newdigits.size());
749 for (unsigned int i = 0; i < newdigits.size(); ++i) {
750 weighteddigits.at(i) = std::make_pair((newdigits.at(i)).getEnergy(), newweights.at(i));
751 listCrystalPairs.at(i) = std::make_pair((newdigits.at(i)).getCellId(), newweights.at(i) * (newdigits.at(i)).getEnergy());
752 }
753
754
755 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
756 return left.second > right.second;
757 });
758
759 std::vector< unsigned int> listCrystals;
760
761 for (unsigned int i = 0; i < newdigits.size(); ++i) {
762 if (i < nOptimal) {
763 listCrystals.push_back(listCrystalPairs[i].first);
764 }
765 }
766
767 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
768 aECLShower->setListOfCrystalsForEnergy(listCrystals);
769
771 B2DEBUG(175, "Shower Energy (2): " << showerEnergy);
772
773 } else {
774 showerEnergy = Belle2::ECL::computeEnergySum(newdigits, newweights);
775 }
776
777 aECLShower->setEnergy(showerEnergy);
778 aECLShower->setEnergyRaw(showerEnergy);
779 aECLShower->setEnergyHighestCrystal(highestEnergy);
780 aECLShower->setTime(highestEnergyTime);
781 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
782 aECLShower->setNumberOfCrystals(weightSum);
783 aECLShower->setCentralCellId(locmaxcellid);
784
785 B2DEBUG(175, "new energy: " << showerEnergy);
786
787
788 aECLShower->setShowerId(iShower);
789 ++iShower;
791 aECLShower->setConnectedRegionId(aCR.
getCRId());
792
793
794 aECLShower->addRelationTo(&aCR);
795
796
798 }
799 }
800}
DataType Phi() const
The azimuth angle.
DataType Z() const
access variable Z (= .at(2) without boundary check)
DataType Theta() const
The polar angle.
DataType X() const
access variable X (= .at(0) without boundary check)
DataType Y() const
access variable Y (= .at(1) without boundary check)
DataType Mag() const
The magnitude (rho in spherical coordinate system).
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.
int getCRId() const
Get CR identifieer.
Class to store local maxima (LM)
@ c_nPhotons
CR is split into n photons (N1)
int m_maxIterations
Maximum number of iterations.
double m_minimumSharedEnergy
Minimum shared energy.
double m_threshold
Local maximum threshold after splitting.
double estimateEnergy(const int centerid)
Estimate energy using 3x3 around central crystal.
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...
double m_shiftTolerance
Tolerance level for centroid shifts.
int m_fullBkgdCount
Number of expected background digits at full background, FIXME: move to database.
const double c_molierRadius
Constant RM (Molier Radius) from exp(-a*dist/RM), http://pdg.lbl.gov/2009/AtomicNuclearProperties/HTM...
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 getEnergySum(std::vector< std::pair< double, double > > &weighteddigits, const unsigned int n)
Get energy sum for weighted entries.
double m_expConstant
Constant a from exp(-a*dist/RM), 1.5 to 2.5.
ROOT::Math::XYZVector GetCrystalPos(int cid)
The Position of crystal.
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 MeV
[megaelectronvolt]
B2Vector3< double > B2Vector3D
typedef for common usage with double