Split connected region into showers.
217{
218
219
221 double backgroundLevel = 0.0;
223 backgroundLevel =
static_cast<double>(bkgdcount) /
static_cast<double>(
m_fullBkgdCount);
224 }
225
226
228
229 B2DEBUG(170, "ECLCRSplitterModule::splitConnectedRegion: nLocalMaximums = " << nLocalMaximums);
230
231
232
233
234
235
236
237
238 if (nLocalMaximums == 1) {
239
240
242
243
244 aECLShower->addRelationTo(&aCR);
245
246
247 double weightSum = 0.0;
248
249
252
253 const int locmaxcellid = locmaxvector[0]->getCellId();
258 double highestEnergyTimeResolution = (
m_eclCalDigits[pos])->getTimeResolution();
259
260
262
263
268
269
270 std::vector<ECLCalDigit> digits;
271 std::vector<double> weights;
272 for (
auto& neighbourId : neighbourMap->
getNeighbours(highestEnergyID)) {
274 neighbourId);
276
279 weights.push_back(1.0);
280 weightSum += 1.0;
281
283 }
284
285
287 aECLShower->setTheta(showerposition.
Theta());
288 aECLShower->setPhi(showerposition.
Phi());
289 aECLShower->setR(showerposition.
Mag());
290
291
292 double showerEnergy = 0.0;
294
295
297 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
298 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
299
300
301
302 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
303 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
304 aECLShower->setNOptimalEnergy(energyEstimation);
305
306
307 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
308 listCrystalPairs.resize(digits.size());
309
310 std::vector < std::pair<double, double> > weighteddigits;
311 weighteddigits.resize(digits.size());
312 for (unsigned int i = 0; i < digits.size(); ++i) {
313 weighteddigits.at(i) = std::make_pair((digits.at(i)).getEnergy(), weights.at(i));
314 listCrystalPairs.at(i) = std::make_pair((digits.at(i)).getCellId(), weights.at(i) * (digits.at(i)).getEnergy());
315 }
316
317
318 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
319 return left.second > right.second;
320 });
321 std::vector< unsigned int> listCrystals;
322
323 for (unsigned int i = 0; i < digits.size(); ++i) {
324 if (i < nOptimal) {
325 listCrystals.push_back(listCrystalPairs[i].first);
326 }
327 }
328
329 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
330 aECLShower->setListOfCrystalsForEnergy(listCrystals);
331
333 B2DEBUG(175, "Shower Energy (1): " << showerEnergy);
334
335 } else {
336 showerEnergy = Belle2::ECL::computeEnergySum(digits, weights);
337 }
338
339 aECLShower->setEnergy(showerEnergy);
340 aECLShower->setEnergyRaw(showerEnergy);
341 aECLShower->setEnergyHighestCrystal(highestEnergy);
342 aECLShower->setTime(highestEnergyTime);
343 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
344 aECLShower->setNumberOfCrystals(weightSum);
345 aECLShower->setCentralCellId(highestEnergyID);
346
347 B2DEBUG(175,
"theta = " << showerposition.
Theta());
348 B2DEBUG(175,
"phi = " << showerposition.
Phi());
349 B2DEBUG(175,
"R = " << showerposition.
Mag());
350 B2DEBUG(175, "energy = " << showerEnergy);
351 B2DEBUG(175, "time = " << highestEnergyTime);
352 B2DEBUG(175, "time resolution = " << highestEnergyTimeResolution);
353 B2DEBUG(175, "neighbours = " << nNeighbours);
354 B2DEBUG(175, "backgroundLevel = " << backgroundLevel);
355
356
357 aECLShower->setShowerId(1);
359 aECLShower->setConnectedRegionId(aCR.
getCRId());
360
361
366 }
367
368 }
369 else {
370
371
372
373 std::vector<std::pair<ECLLocalMaximum, double>> lm_energy_vector;
375 const int cellid = aLocalMaximum.getCellId();
378 lm_energy_vector.push_back(std::pair<ECLLocalMaximum, double>(aLocalMaximum, digitenergy));
379 };
380
381
382 if (lm_energy_vector.size() >=
static_cast<size_t>(
m_maxSplits)) {
383 std::sort(lm_energy_vector.begin(), lm_energy_vector.end(), [](const std::pair<ECLLocalMaximum, double>& x,
384 const std::pair<ECLLocalMaximum, double>& y) {
385 return x.second > y.second;
386 });
387
389 }
390
391 std::vector<ECLCalDigit> digits;
392 std::vector<double> weights;
393 std::map<int, B2Vector3D> centroidList;
394 std::map<int, double> centroidEnergyList;
395 std::map<int, B2Vector3D> allPoints;
396 std::map<int, std::vector < double > > weightMap;
397 std::vector < ECLCalDigit > digitVector;
398
399
400 std::map<int, B2Vector3D> localMaximumsPoints;
401 std::map<int, B2Vector3D> centroidPoints;
402
403 for (auto& aLocalMaximum : lm_energy_vector) {
404
405 int cellid = aLocalMaximum.first.getCellId();
406
407
409 localMaximumsPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
410 centroidPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
411 }
412
413
414 bool iterateclusters = true;
415 do {
416 digits.clear();
417 weights.clear();
418 centroidList.clear();
419 centroidEnergyList.clear();
420 allPoints.clear();
421 weightMap.clear();
422 digitVector.clear();
423
424
426 const int cellid = aCalDigit.getCellId();
427
429 allPoints.insert(std::map<int, B2Vector3D>::value_type(cellid, vectorPosition));
430 digits.push_back(aCalDigit);
431 }
432
433 for (const auto& digitpoint : allPoints) {
434 const int cellid = digitpoint.first;
437 }
438
439
440
441
442
443
444 int nIterations = 0;
445 double centroidShiftAverage = 0.0;
446
447
448 do {
449 B2DEBUG(175,
"Iteration: #" << nIterations <<
" (of max. " <<
m_maxIterations <<
")");
450
451 centroidShiftAverage = 0.0;
452 if (nIterations == 0) {
453 centroidList.clear();
454 centroidEnergyList.clear();
455 weightMap.clear();
456 }
457
458
459 for (const auto& locmaxpoint : localMaximumsPoints) {
460
461
462 const int locmaxcellid = locmaxpoint.first;
463
464
465 weights.clear();
466
467
468 if (nIterations == 0) {
471 centroidEnergyList[locmaxcellid] = 1.5 * locmaxenergy;
472 }
473
474 B2DEBUG(175, "local maximum cellid: " << locmaxcellid);
475
476
477
478 for (const auto& digitpoint : allPoints) {
479
480
481 const int digitcellid = digitpoint.first;
483
486
487 double weight = 0.0;
488 double energy = 0.0;
489 double distance = 0.0;
490 double distanceEnergySum = 0.0;
491
492
493 for (const auto& centroidpoint : centroidPoints) {
494
495
496 const int centroidcellid = centroidpoint.first;
497 B2Vector3D centroidpos = centroidpoint.second;
498
499 double thisdistance = 0.;
500
501
502 if (nIterations == 0 and digitcellid == centroidcellid) {
503 thisdistance = 0.0;
504 } else {
505 B2Vector3D vectorDistance = ((centroidpos) - (digitpos));
506 thisdistance = vectorDistance.
Mag();
507 }
508
509
512
513
514 if (locmaxcellid == centroidcellid) {
515 distance = thisdistance;
516 energy = thisenergy;
517 }
518
519
521 distanceEnergySum += (thisenergy * expfactor);
522
523 }
524
525
526 if (distanceEnergySum > 0.0) {
528 weight = energy * expfactor / distanceEnergySum;
529 } else {
530 weight = 0.0;
531 }
532
533
535 weight = 0.0;
536 }
537
538
539 if (weight > 1.0) {
540 B2WARNING("ECLCRSplitterModule::splitConnectedRegion: Floating point glitch, weight for this digit " << weight <<
541 ", resetting it to 1.0.");
542 weight = 1.0;
543 }
544
545
546 B2DEBUG(175, " cellid: " << digitcellid << ", energy: " << digitenergy << ", weight: " << weight << ", distance: " << distance);
547 weights.push_back(weight);
548
549 }
550
551
552 B2Vector3D oldCentroidPos = (centroidPoints.find(locmaxcellid))->second;
553
554
556
557
558 const double newEnergy = Belle2::ECL::computeEnergySum(digits, weights);
559
560
561 const B2Vector3D centroidShift = (oldCentroidPos - newCentroidPos);
562
563
564 centroidList[locmaxcellid] = newCentroidPos;
565 weightMap[locmaxcellid] = weights;
566
567 B2DEBUG(175, "--> inserting new energy: " << newEnergy << " for local maximum " << locmaxcellid);
568 centroidEnergyList[locmaxcellid] = newEnergy;
569 double showerenergy = (*centroidEnergyList.find(locmaxcellid)).second /
Belle2::Unit::MeV;
570 B2DEBUG(175, "--> new energy = " << showerenergy << " MeV");
571
572
573 centroidShiftAverage += centroidShift.
Mag();
574
575
576 B2DEBUG(175,
" old centroid: " << oldCentroidPos.
X() <<
" cm, " << oldCentroidPos.
Y() <<
" cm, " << oldCentroidPos.
Z() <<
577 "cm");
578 B2DEBUG(175,
" new centroid: " << newCentroidPos.
X() <<
" cm, " << newCentroidPos.
Y() <<
" cm, " << newCentroidPos.
Z() <<
579 "cm");
580 B2DEBUG(175,
" centroid shift: " << centroidShift.
Mag() <<
" cm");
581
582 }
583
584
585 centroidShiftAverage /= static_cast<double>(nLocalMaximums);
586
587 B2DEBUG(175,
"--> average centroid shift: " << centroidShiftAverage <<
" cm (tolerance is " <<
m_shiftTolerance <<
" cm)");
588
589
590 for (const auto& locmaxpoint : localMaximumsPoints) {
591 centroidPoints[locmaxpoint.first] = (centroidList.find(locmaxpoint.first))->second;
592 }
593
594 ++nIterations;
595
596 }
while (nIterations < m_maxIterations and centroidShiftAverage >
m_shiftTolerance);
597
598
599
600 std::vector<int> markfordeletion;
601 iterateclusters = false;
602 for (const auto& locmaxpoint : localMaximumsPoints) {
603
604
605 const int locmaxcellid = locmaxpoint.first;
607
608 B2DEBUG(175, "locmaxcellid: " << locmaxcellid);
610 B2DEBUG(175, "ok: ");
611
612
613 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
614
615 for (unsigned int i = 0; i < digitVector.size(); ++i) {
616
618 const double weight = myWeights[i];
621
622
623
624
625 if ((energy * weight > lmenergy and cellid != locmaxcellid) or (energy * weight <
m_threshold and cellid == locmaxcellid)) {
626 markfordeletion.push_back(locmaxcellid);
627 iterateclusters = true;
628 continue;
629 }
630 }
631 }
632
633
634 for (const auto lmid : markfordeletion) {
635 localMaximumsPoints.erase(lmid);
636 centroidPoints.erase(lmid);
637 }
638
639 } while (iterateclusters);
640
641
642 unsigned int iShower = 1;
643 for (const auto& locmaxpoint : localMaximumsPoints) {
644
645 const int locmaxcellid = locmaxpoint.first;
647
648
650
651
653
654
659
660
661 std::vector<short int> neighbourlist = neighbourMap->
getNeighbours(locmaxcellid);
662
663
664 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
665
666
667 std::vector<ECLCalDigit> newdigits;
668 std::vector<double> newweights;
669 double highestEnergy = 0.;
670 double highestEnergyTime = 0.;
671 double highestEnergyTimeResolution = 0.;
672 double weightSum = 0.0;
673
674 for (unsigned int i = 0; i < digitVector.size(); ++i) {
675
677 const double weight = myWeights[i];
678
681
682
684
685
686 if (weight > 0.0) {
687 if (std::find(neighbourlist.begin(), neighbourlist.end(), cellid) != neighbourlist.end()) {
688
690
691 newdigits.push_back(dig);
692 newweights.push_back(weight);
693
694 weightSum += weight;
695
697
698 if (energy * weight > highestEnergy) {
699 highestEnergy = energy * weight;
700 highestEnergyTime = dig.
getTime();
702 }
703 }
704 }
705 }
706
707
709
710 B2DEBUG(175,
"old theta: " << oldshowerposition->
Theta());
711 B2DEBUG(175,
"old phi: " << oldshowerposition->
Phi());
712 B2DEBUG(175,
"old R: " << oldshowerposition->
Mag());
713 B2DEBUG(175, "old energy: " << energyEstimation);
714 delete oldshowerposition;
715
716
717
718
720 aECLShower->setTheta(showerposition->
Theta());
721 aECLShower->setPhi(showerposition->
Phi());
722 aECLShower->setR(showerposition->
Mag());
723
724 B2DEBUG(175,
"new theta: " << showerposition->
Theta());
725 B2DEBUG(175,
"new phi: " << showerposition->
Phi());
726 B2DEBUG(175,
"new R: " << showerposition->
Mag());
727 delete showerposition;
728
729
730 double showerEnergy = 0.0;
732
733
734
736 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
737 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
738
739
740
741 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
742 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
743 aECLShower->setNOptimalEnergy(energyEstimation);
744
745
746 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
747 listCrystalPairs.resize(newdigits.size());
748
749 std::vector < std::pair<double, double> > weighteddigits;
750 weighteddigits.resize(newdigits.size());
751 for (unsigned int i = 0; i < newdigits.size(); ++i) {
752 weighteddigits.at(i) = std::make_pair((newdigits.at(i)).getEnergy(), newweights.at(i));
753 listCrystalPairs.at(i) = std::make_pair((newdigits.at(i)).getCellId(), newweights.at(i) * (newdigits.at(i)).getEnergy());
754 }
755
756
757 std::sort(listCrystalPairs.begin(), listCrystalPairs.end(), [](const auto & left, const auto & right) {
758 return left.second > right.second;
759 });
760
761 std::vector< unsigned int> listCrystals;
762
763 for (unsigned int i = 0; i < newdigits.size(); ++i) {
764 if (i < nOptimal) {
765 listCrystals.push_back(listCrystalPairs[i].first);
766 }
767 }
768
769 aECLShower->setNumberOfCrystalsForEnergy(static_cast<double>(listCrystals.size()));
770 aECLShower->setListOfCrystalsForEnergy(listCrystals);
771
773 B2DEBUG(175, "Shower Energy (2): " << showerEnergy);
774
775 } else {
776 showerEnergy = Belle2::ECL::computeEnergySum(newdigits, newweights);
777 }
778
779 aECLShower->setEnergy(showerEnergy);
780 aECLShower->setEnergyRaw(showerEnergy);
781 aECLShower->setEnergyHighestCrystal(highestEnergy);
782 aECLShower->setTime(highestEnergyTime);
783 aECLShower->setDeltaTime99(highestEnergyTimeResolution);
784 aECLShower->setNumberOfCrystals(weightSum);
785 aECLShower->setCentralCellId(locmaxcellid);
786
787 B2DEBUG(175, "new energy: " << showerEnergy);
788
789
790 aECLShower->setShowerId(iShower);
791 ++iShower;
793 aECLShower->setConnectedRegionId(aCR.
getCRId());
794
795
796 aECLShower->addRelationTo(&aCR);
797
798
800 }
801 }
802}
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 identifier.
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