Belle II Software  release-05-01-25
Position.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Calculation of shower positions using lists of digits and weights. *
6  * *
7  * Author: The Belle II Collaboration *
8  * Contributors: Torben Ferber (ferber@physics.ubc.ca) *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
13 // THIS
14 #include "ecl/utility/Position.h"
15 
16 // GEOMETRY
17 #include <ecl/geometry/ECLGeometryPar.h>
18 
19 // FRAMEWORK
20 #include <framework/logging/Logger.h>
21 
22 // OTHER
23 #include "TMath.h"
24 
25 namespace Belle2 {
30  namespace ECL {
31  // computePositionLiLo
32  B2Vector3D computePositionLiLo(std::vector<ECLCalDigit>& digits, std::vector<double>& weights, std::vector<double>& parameters)
33  {
34 
35  // Total weighted sum.
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)");
39  B2Vector3D liloPointError(0, 0, 0);
40  return liloPointError;
41  }
42 
43  // Offset from user input constants: offset = A-B * exp(-C * energysum).
44  const double offset = parameters[0] - parameters[1] * TMath::Exp(-parameters[2] * energySum);
45 
46  // Geometry information for crystal positions.
47  ECLGeometryPar* geom = ECLGeometryPar::Instance();
48 
49  // Log-weighted position calculation
50  B2Vector3D liloPoint(1, 1, 1); // anything is fine as long its not zero length
51  B2Vector3D linearPositionVector(0, 0, 0);
52  B2Vector3D logPositionVector(0, 0, 0);
53  double logWeightSum = 0.0;
54  int nTotal = 0;
55  int nInLogWeightSum = 0;
56  double firstPhi = -999.;
57  double firstTheta = -999.;
58  bool foundSecondPhi = false;
59  bool foundSecondTheta = false;
60 
61  // Loop over all digits in the vector.
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); // v = crystal center - (0, 0, 0)
67  const double theta = position.Theta();
68  const double phi = position.Phi();
69  ++nTotal;
70 
71  // Weights.
72  const double linearWeight = energy * weight / energySum; // fraction of this digit energy to the total shower energy
73  const double logWeight = offset + TMath::Log(linearWeight);
74 
75  // Update the linear position vector.
76  linearPositionVector += linearWeight * position;
77 
78  // Only digits with a positive weight are included.
79  if (logWeight > 0.0) {
80 
81  // Update the log position vector.
82  logPositionVector += logWeight * position;
83  logWeightSum += logWeight;
84  ++nInLogWeightSum;
85 
86  // We have to avoid that only crystals with the same theta or phi values are used (this will bias the position towards the crystal center
87  if (nInLogWeightSum == 1) {
88  firstTheta = theta;
89  firstPhi = phi;
90  } else {
91  if (!foundSecondPhi and fabs(phi - firstPhi) > 1e-4) { // we cant use phi id in the endcaps
92  foundSecondPhi = true;
93  }
94  if (!foundSecondTheta
95  and fabs(theta - firstTheta) > 1e-4) { // we could probably use thetaid since this works in the endcaps as well
96  foundSecondTheta = true;
97  }
98  }
99  } // end if logWeight
100  } // end digit
101 
102  // Check if at least one digit has a positive weight in the logweightsum.
103  if (nInLogWeightSum > 0) logPositionVector *= 1. / logWeightSum;
104 
105  // The direction is filled, now get the distance to the center. It can happen, that it is outside of the crystals
106  // if the shower is in the barrel/endcap overlap regions?!
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());
110 
111  return liloPoint;
112  }
113 
114  // helper: computeEnergySum
115  double computeEnergySum(std::vector<ECLCalDigit>& digits, std::vector<double>& weights)
116  {
117  int n = 0;
118  double sum = 0.0;
119 
120  for (const auto& digit : digits) {
121  sum += digit.getEnergy() * weights.at(n);
122  ++n;
123  }
124 
125  return sum;
126  }
127 
128 // // helper: findClosestCrystal
129 // int findClosestCrystal(std::vector<ECLCalDigit>& digits, TVector3& direction)
130 // {
131 // ECLGeometryPar* geom = ECLGeometryPar::Instance();
132 //
133 // int best = -1;
134 // double min = 999.;
135 // for (const auto& digit : digits) {
136 // const int crystalid = digit.getCellId() - 1;
137 // const double alpha = direction.Angle(geom->GetCrystalVec(crystalid));
138 // const double R = (geom->GetCrystalPos(crystalid)).Mag();
139 // const double distance = 2.0 * R * TMath::Sin(alpha / 2.0); // chord distance
140 //
141 // if (distance < min) {
142 // min = distance;
143 // best = crystalid;
144 // }
145 // }
146 // return best;
147 // }
148 
149  // ...
150 
151  }
153 }
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition: ECLGeometryPar.cc:151
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:507
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19