14 #include "ecl/utility/Position.h"
17 #include <ecl/geometry/ECLGeometryPar.h>
20 #include <framework/logging/Logger.h>
32 B2Vector3D computePositionLiLo(std::vector<ECLCalDigit>& digits, std::vector<double>& weights, std::vector<double>& parameters)
36 const double energySum = computeEnergySum(digits, weights);
37 if (energySum <= 0.0) {
38 B2WARNING(
"ECL computePositionLiLo() energy sum zero or negative: " << energySum <<
", return (0, 0, 0)");
40 return liloPointError;
44 const double offset = parameters[0] - parameters[1] * TMath::Exp(-parameters[2] * energySum);
53 double logWeightSum = 0.0;
55 int nInLogWeightSum = 0;
56 double firstPhi = -999.;
57 double firstTheta = -999.;
58 bool foundSecondPhi =
false;
59 bool foundSecondTheta =
false;
62 for (
const auto& digit : digits) {
63 const double energy = digit.getEnergy();
64 const int cellid = digit.getCellId();
65 const double weight = weights.at(nTotal);
66 const B2Vector3D position = geom->GetCrystalPos(cellid - 1);
67 const double theta = position.Theta();
68 const double phi = position.Phi();
72 const double linearWeight = energy * weight / energySum;
73 const double logWeight = offset + TMath::Log(linearWeight);
76 linearPositionVector += linearWeight * position;
79 if (logWeight > 0.0) {
82 logPositionVector += logWeight * position;
83 logWeightSum += logWeight;
87 if (nInLogWeightSum == 1) {
91 if (!foundSecondPhi and fabs(phi - firstPhi) > 1e-4) {
92 foundSecondPhi =
true;
95 and fabs(theta - firstTheta) > 1e-4) {
96 foundSecondTheta =
true;
103 if (nInLogWeightSum > 0) logPositionVector *= 1. / logWeightSum;
107 liloPoint.SetTheta(foundSecondTheta ? logPositionVector.Theta() : linearPositionVector.Theta());
108 liloPoint.SetPhi(foundSecondPhi ? logPositionVector.Phi() : linearPositionVector.Phi());
109 liloPoint.SetMag(nInLogWeightSum > 0 ? logPositionVector.Mag() : linearPositionVector.Mag());
115 double computeEnergySum(std::vector<ECLCalDigit>& digits, std::vector<double>& weights)
120 for (
const auto& digit : digits) {
121 sum += digit.getEnergy() * weights.at(n);