Split connected region into showers.
220{
221
222
224 if (maxDistance < 0.0) {
225 maxDistance = 9999999.;
226 }
227
228
229
231 double backgroundLevel = 0.0;
233 backgroundLevel =
static_cast<double>(bkgdcount) /
static_cast<double>(
m_fullBkgdCount);
234 }
235
236
238
239 B2DEBUG(170, "ECLCRSplitterModule::splitConnectedRegion: nLocalMaximums = " << nLocalMaximums);
240
241
242
243
244
245
246
247
248 if (nLocalMaximums == 1) {
249
250
252
253
254 aECLShower->addRelationTo(&aCR);
255
256
257 double weightSum = 0.0;
258
259
261 aECLShower->addRelationTo(locmaxvector[0]);
262
263 const int locmaxcellid = locmaxvector[0]->getCellId();
268 double highestEnergyTimeResolution = (
m_eclCalDigits[pos])->getTimeResolution();
269
270
272
273
274 ECLNeighbours* neighbourMap;
278
279
280 std::vector<ECLCalDigit> digits;
281 std::vector<double> weights;
282 for (
auto& neighbourId : neighbourMap->
getNeighbours(highestEnergyID)) {
284 neighbourId);
286
289 weights.push_back(1.0);
290 weightSum += 1.0;
291
293 }
294
295
297 aECLShower->setTheta(showerposition.
Theta());
298 aECLShower->setPhi(showerposition.
Phi());
299 aECLShower->setR(showerposition.
Mag());
300
301
302 double showerEnergy = 0.0;
304
305
307 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
308 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
309
310
311
312 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
313 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
314 aECLShower->setNOptimalEnergy(energyEstimation);
315
316
317 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
318 listCrystalPairs.resize(digits.size());
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
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;
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
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
367 aECLShower->setShowerId(1);
369 aECLShower->setConnectedRegionId(aCR.
getCRId());
370
371
376 }
377
378 }
379 else {
380
381
382
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));
389 };
390
391
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
399 }
400
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;
408
409
410 std::map<int, B2Vector3D> localMaximumsPoints;
411 std::map<int, B2Vector3D> centroidPoints;
412
413 for (auto& aLocalMaximum : lm_energy_vector) {
414
415 int cellid = aLocalMaximum.first.getCellId();
416
417
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
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
436 const int cellid = aCalDigit.getCellId();
437
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;
447 }
448
449
450
451
452
453
454 int nIterations = 0;
455 double centroidShiftAverage = 0.0;
456
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
469 for (const auto& locmaxpoint : localMaximumsPoints) {
470
471
472 const int locmaxcellid = locmaxpoint.first;
473
474
475 weights.clear();
476
477
478 if (nIterations == 0) {
481 centroidEnergyList[locmaxcellid] = 1.5 * locmaxenergy;
482 }
483
484 B2DEBUG(175, "local maximum cellid: " << locmaxcellid);
485
486
487
488 for (const auto& digitpoint : allPoints) {
489
490
491 const int digitcellid = digitpoint.first;
493
496
497 double weight = 0.0;
498 double energy = 0.0;
499 double distance = 0.0;
500 double distanceEnergySum = 0.0;
501
502
503 for (const auto& centroidpoint : centroidPoints) {
504
505
506 const int centroidcellid = centroidpoint.first;
507 B2Vector3D centroidpos = centroidpoint.second;
508
509 double thisdistance = 0.;
510
511
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
522
523
524
525 if ((locmaxcellid == centroidcellid) and (thisdistance < maxDistance)) {
526 distance = thisdistance;
527 energy = thisenergy;
528 }
529
530
531 if (thisdistance < maxDistance) {
533 distanceEnergySum += (thisenergy * expfactor);
534 }
535
536 }
537
538
539 if (distanceEnergySum > 0.0) {
541 weight = energy * expfactor / distanceEnergySum;
542 } else {
543 weight = 0.0;
544 }
545
546
548 weight = 0.0;
549 }
550
551
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
559 B2DEBUG(175, " cellid: " << digitcellid << ", energy: " << digitenergy << ", weight: " << weight << ", distance: " << distance);
560 weights.push_back(weight);
561
562 }
563
564
565 B2Vector3D oldCentroidPos = (centroidPoints.find(locmaxcellid))->second;
566
567
569
570
571 const double newEnergy = Belle2::ECL::computeEnergySum(digits, weights);
572
573
574 const B2Vector3D centroidShift = (oldCentroidPos - newCentroidPos);
575
576
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
586 centroidShiftAverage += centroidShift.
Mag();
587
588
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 }
596
597
598 centroidShiftAverage /= static_cast<double>(nLocalMaximums);
599
600 B2DEBUG(175,
"--> average centroid shift: " << centroidShiftAverage <<
" cm (tolerance is " <<
m_shiftTolerance <<
" cm)");
601
602
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
611
612
613 std::vector<int> markfordeletion;
614 iterateclusters = false;
615 for (const auto& locmaxpoint : localMaximumsPoints) {
616
617
618 const int locmaxcellid = locmaxpoint.first;
620
621 B2DEBUG(175, "locmaxcellid: " << locmaxcellid);
623 B2DEBUG(175, "ok: ");
624
625
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];
634
635
636
637
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 }
644 }
645 }
646
647
648 for (const auto lmid : markfordeletion) {
649 localMaximumsPoints.erase(lmid);
650 centroidPoints.erase(lmid);
651 }
652
653 } while (iterateclusters);
654
655
656 unsigned int iShower = 1;
657 for (const auto& locmaxpoint : localMaximumsPoints) {
658
659 const int locmaxcellid = locmaxpoint.first;
661
662
664
665
667
668
669 ECLNeighbours* neighbourMap;
673
674
675 std::vector<short int> neighbourlist = neighbourMap->
getNeighbours(locmaxcellid);
676
677
678 std::vector < double > myWeights = (*weightMap.find(locmaxcellid)).second;
679
680
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
695
696
698
699
700 if (weight > 0.0) {
701 if (std::find(neighbourlist.begin(), neighbourlist.end(), cellid) != neighbourlist.end()) {
702
704
705 newdigits.push_back(dig);
706 newweights.push_back(weight);
707
708 weightSum += weight;
709
711
712 if (energy * weight > highestEnergy) {
713 highestEnergy = energy * weight;
714 highestEnergyTime = dig.
getTime();
716 }
717 }
718 }
719 }
720
721
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
731
732
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
744 double showerEnergy = 0.0;
746
747
748
750 const unsigned int nOptimal = static_cast<unsigned int>(nOptimalVec[0]);
751 aECLShower->setNominalNumberOfCrystalsForEnergy(static_cast<double>(nOptimal));
752
753
754
755 aECLShower->setNOptimalGroupIndex(nOptimalVec[1]);
756 aECLShower->setNOptimalEnergyBin(nOptimalVec[2]);
757 aECLShower->setNOptimalEnergy(energyEstimation);
758
759
760 std::vector< std::pair<unsigned int, double>> listCrystalPairs;
761 listCrystalPairs.resize(newdigits.size());
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
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;
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
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
804 aECLShower->setShowerId(iShower);
805 ++iShower;
807 aECLShower->setConnectedRegionId(aCR.
getCRId());
808
809
810 aECLShower->addRelationTo(&aCR);
811
812
814 }
815 }
816}
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.
@ c_nPhotons
CR is split into n photons (N1)
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...
const double c_molierRadius
Constant RM (Molier Radius) from exp(-a*dist/RM), http://pdg.lbl.gov/2009/AtomicNuclearProperties/HTM...
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.
const std::vector< short int > & getNeighbours(short int cid) const
Return the neighbours for a given cell ID.
RelationVector< T > getRelationsWith(const std::string &name="", const std::string &namedRelation="") const
Get the relations between this object and another store array.
B2Vector3< double > B2Vector3D
typedef for common usage with double
1.13.2