ATLAS Offline Software
MultiComponentStateModeCalculator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
15 //
16 #include "CxxUtils/phihelper.h"
17 #include <cmath>
18 
19 #include<boost/container/static_vector.hpp>
20 
21 
22 namespace {
23 constexpr double invsqrt2PI =
24  M_2_SQRTPI / (2. * M_SQRT2); // 1. / sqrt(2. * M_PI);
25 
26 // Simple representation of 1D component
27 struct Component
28 {
29  Component() = default;
30  ~Component() = default;
31  Component(const Component&) = default;
32  Component& operator=(const Component&) = default;
33  Component(Component&&) = default;
34  Component& operator=(Component&&) = default;
35  // Constructor with arguments
36  Component(double aWeight, double aMean, double aSigma)
37  : weight(aWeight)
38  , mean(aMean)
39  , sigma(aSigma)
40  {}
41  double weight = 0;
42  double mean = 0;
43  double sigma = 0;
44 };
45 
46 using VecOfComponents =
47  boost::container::static_vector<Component,
49 struct pdfAndDeriv
50 {
51  double value = 0.;
52  double deriv1 = 0.;
53  double deriv2 = .0;
54 };
55 
58 double
59 gaus(double x, double mean, double sigma)
60 {
61  const double invertsigma = 1. / sigma;
62  const double z = (x - mean) * invertsigma;
63  return (invsqrt2PI * invertsigma) * exp(-0.5 * z * z);
64 }
65 
67 double
68 pdf(double x, int i, const std::array<VecOfComponents, 5>& mixture)
69 {
70  double pdf(0.);
71  auto component = mixture[i].begin();
72  for (; component != mixture[i].end(); ++component) {
73  pdf += component->weight * gaus(x, component->mean, component->sigma);
74  }
75  return pdf;
76 }
77 
80 pdfAndDeriv
81 fullPdf(double x, int i, const std::array<VecOfComponents, 5>& mixture)
82 {
83  pdfAndDeriv pdf{};
84  auto component = mixture[i].begin();
85  for (; component != mixture[i].end(); ++component) {
86  double componentgaus = gaus(x, component->mean, component->sigma);
87  pdf.value += component->weight * componentgaus;
88  const double invertSigma = 1. / component->sigma;
89  const double z = (x - component->mean) * invertSigma;
90  pdf.deriv1 += -1. * component->weight * z * componentgaus * invertSigma;
91  pdf.deriv2 += component->weight * invertSigma * invertSigma * (z * z - 1.) *
92  componentgaus;
93  }
94  return pdf;
95 }
96 
99 double
100 width(int i, const std::array<VecOfComponents, 5>& mixture)
101 {
102  double pdf(0.);
103  auto component = mixture[i].begin();
104  for (; component != mixture[i].end(); ++component) {
105  pdf += component->weight * component->sigma;
106  }
107  return pdf;
108 }
109 
110 void
111 fillMixture(const Trk::MultiComponentState& multiComponentState,
112  std::array<VecOfComponents, 5>& mixture)
113 {
114  constexpr Trk::ParamDefs parameter[5] = {
116  };
117 
118  // Loop over all the components in the multi-component state
119  for (const Trk::ComponentParameters& component : multiComponentState) {
120 
121  // And then for each component over each 5 parameters
122  for (size_t i = 0; i < 5; ++i) {
123  const Trk::TrackParameters* componentParameters = component.params.get();
124  const AmgSymMatrix(5)* measuredCov = componentParameters->covariance();
125  if (!measuredCov) {
126  return;
127  }
128  // Enums for Perigee
129  // d0=0, z0=1, phi0=2, theta=3, qOverP=4,
130  double weight = component.weight;
131  double mean = componentParameters->parameters()[parameter[i]];
132  // FIXME ATLASRECTS-598 this std::abs() should not be necessary... for
133  // some reason cov(qOverP,qOverP) can be negative
134  double sigma = sqrt(std::abs((*measuredCov)(parameter[i], parameter[i])));
135  // Ensure that we don't have any problems with the cyclical nature of
136  // phi Use first state as reference point
137  if (i == 2) { // phi
138  double deltaPhi =
139  multiComponentState.begin()->params->parameters()[2] - mean;
140  if (deltaPhi > M_PI) {
141  mean += 2 * M_PI;
142  } else if (deltaPhi < -M_PI) {
143  mean -= 2 * M_PI;
144  }
145  }
146  mixture[i].emplace_back(weight, mean, sigma);
147  }
148  }
149 }
150 
151 double
152 findMode(double xStart,
153  int i,
154  const std::array<VecOfComponents, 5>& mixture)
155 {
156 
157  bool converged = false;
158  double tolerance(1.);
159  // start position for mode
160  double currentMode(xStart);
161  double nextMode(currentMode);
162  // pdf at current mode and next mode
163  pdfAndDeriv currentPdf = fullPdf(currentMode, i, mixture);
164 
165  // Allow up to 20 iterations for convergence
166  for (int iteration = 0; iteration < 20; ++iteration) {
167  // calculate next mode point
168  if (currentPdf.deriv2 != 0.0) {
169  nextMode = currentMode - currentPdf.deriv1 / currentPdf.deriv2;
170  } else {
171  return xStart;
172  }
173 
174  // Calculate the mixture pdf at next mode point
175  pdfAndDeriv nextPdf = fullPdf(nextMode, i, mixture);
176  // check if we have converged
177  if ((nextPdf.value + currentPdf.value) != 0.0) {
178  tolerance = std::abs(nextPdf.value - currentPdf.value) /
179  (nextPdf.value + currentPdf.value);
180  } else {
181  return xStart;
182  }
183  if (tolerance < 1.e-8) {
184  converged = true;
185  break;
186  }
187  // if we have not yet converged
188  // next becomes current and we retry
189  currentPdf = nextPdf;
190  currentMode = nextMode;
191  }
192 
193  if (!converged) {
194  return xStart;
195  }
196 
197  return currentMode;
198 }
199 
200 double
201 findRoot(double& result,
202  double xlo,
203  double xhi,
204  double value,
205  double i,
206  const std::array<VecOfComponents, 5>& mixture)
207 {
208  // Do the root finding using the Brent-Decker method. Returns a boolean
209  // status and loads 'result' with our best guess at the root if true.
210  double a(xlo);
211  double b(xhi);
212  double fa = pdf(a, i, mixture) - value;
213  double fb = pdf(b, i, mixture) - value;
214 
215  if (fb * fa > 0) {
216  return false;
217  }
218  bool ac_equal(false);
219  double fc = fb;
220  double c(0);
221  double d(0);
222  double e(0);
223  int MaxIterations = 20;
224  double tolerance = 1.e-6;
225 
226  for (int iter = 0; iter <= MaxIterations; iter++) {
227 
228  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
229  // Rename a,b,c and adjust bounding interval d
230  ac_equal = true;
231  c = a;
232  fc = fa;
233  d = b - a;
234  e = b - a;
235  }
236 
237  if (std::abs(fc) < std::abs(fb)) {
238  ac_equal = true;
239  a = b;
240  b = c;
241  c = a;
242  fa = fb;
243  fb = fc;
244  fc = fa;
245  }
246 
247  double tol = 0.5 * tolerance * std::abs(b);
248  double m = 0.5 * (c - b);
249 
250  if (fb == 0 || std::abs(m) <= tol) {
251  result = b;
252  return true;
253  }
254 
255  if (std::abs(e) < tol || std::abs(fa) <= std::abs(fb)) {
256  // Bounds decreasing too slowly: use bisection
257  d = m;
258  e = m;
259  } else {
260  // Attempt inverse cubic interpolation
261  double p = 0;
262  double q = 0;
263  double r = 0;
264  double s = fb / fa;
265 
266  if (ac_equal) {
267  p = 2 * m * s;
268  q = 1 - s;
269  } else {
270  q = fa / fc;
271  r = fb / fc;
272  p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
273  q = (q - 1) * (r - 1) * (s - 1);
274  }
275  // Check whether we are in bounds
276  if (p > 0) {
277  q = -q;
278  } else {
279  p = -p;
280  }
281 
282  double min1 = 3 * m * q - std::abs(tol * q);
283  double min2 = std::abs(e * q);
284  if (2 * p < (min1 < min2 ? min1 : min2)) {
285  // Accept the interpolation
286  e = d;
287  d = p / q;
288  } else {
289  // Interpolation failed: use bisection.
290  d = m;
291  e = m;
292  }
293  }
294  // Move last best guess to a
295  a = b;
296  fa = fb;
297  // Evaluate new trial root
298  if (std::abs(d) > tol) {
299  b += d;
300  } else {
301  b += (m > 0 ? +tol : -tol);
302  }
303  fb = pdf(b, i, mixture) - value;
304  }
305  // Return our best guess if we run out of iterations
306  result = b;
307 
308  return false;
309 }
310 
314 std::array<double, 10>
315 evaluateMode(const std::array<VecOfComponents, 5>& mixture)
316 {
317  std::array<double, 10> modes{};
318  /* loop over the 5 direction , d0,z0,phi,theta,qOverP*/
319 
320  for (int i = 0; i < 5; i++) {
321 
322  double largerPdfComponent = 0.0;
323  double largerMeanComponent = 0.0;
324  /*
325  * Loop over the mixture in the ith direction and find the component
326  * whose mean give the larger value for the Gaussian Mixture pdf.
327  * This should be a good enough starting point for the mode
328  * finding in this direction
329  */
330  for (const Component& component : mixture[i]) {
331  double pdfValue = pdf(component.mean, i, mixture);
332  if (pdfValue > largerPdfComponent) {
333  largerPdfComponent = pdfValue;
334  largerMeanComponent = component.mean;
335  }
336  }
337  modes[i] = findMode(largerMeanComponent, i, mixture);
338  // Calculate the FWHM and return this back so that it can be used to correct
339  // the covariance matrix
340  if (largerMeanComponent != modes[i]) {
341  // mode calculation was successful now calulate FWHM
342  double currentWidth = width(i, mixture);
343  modes[i + 5] = -1; // Failure is flagged with a value less than 0;
344 
345  double pdfVal = pdf(modes[i], i, mixture);
346  double highX(0);
347  double lowX(0);
348 
349  double upperbound = modes[i] + 1.5 * currentWidth;
350  while (true) {
351  if (pdf(upperbound, i, mixture) > pdfVal * 0.5) {
352  upperbound += currentWidth;
353  } else {
354  break;
355  }
356  }
357 
358  bool highXFound =
359  findRoot(highX, modes[i], upperbound, pdfVal * 0.5, i, mixture);
360 
361  double lowerbound = modes[i] - 1.5 * currentWidth;
362  while (true) {
363  if (pdf(lowerbound, i, mixture) > pdfVal * 0.5) {
364  lowerbound -= currentWidth;
365  } else {
366  break;
367  }
368  }
369  bool lowXFound =
370  findRoot(lowX, lowerbound, modes[i], pdfVal * 0.5, i, mixture);
371  if (highXFound && lowXFound) {
372  double FWHM = highX - lowX;
373  modes[i + 5] = FWHM / 2.35482; // 2 * sqrt( 2* log(2))
374  }
375  // Ensure that phi is between -pi and pi
376  if (i == 2) {
377  modes[i] = CxxUtils::wrapToPi(modes[i]);
378  }
379  }
380  }
381  return modes;
382 }
383 } // end of anonymous namespace
384 
385 std::array<double, 10>
387  const Trk::MultiComponentState& multiComponentState)
388 {
389  // Check to see if all components have covariance
390  if (!MultiComponentStateHelpers::allHaveCovariance(multiComponentState)) {
391  return {};
392  }
393  std::array<VecOfComponents, 5> mixture;
394 
395  fillMixture(multiComponentState, mixture);
396  return evaluateMode(mixture);
397 }
398 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::MultiComponentStateHelpers::allHaveCovariance
bool allHaveCovariance(const MultiComponentState &in)
Check to see if all components have covariance Matrix.
Definition: ComponentParameters.cxx:66
beamspotman.r
def r
Definition: beamspotman.py:676
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
tolerance
constexpr double tolerance
Definition: runMdtGeoComparison.cxx:104
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
get_generator_info.result
result
Definition: get_generator_info.py:21
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
hist_file_dump.d
d
Definition: hist_file_dump.py:137
xAOD::deltaPhi
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
Definition: L2StandAloneMuon_v1.cxx:160
Trk::ParamDefs
ParamDefs
Definition: ParamDefs.h:32
CxxUtils::wrapToPi
T wrapToPi(T phi)
Wrap angle in radians to [-pi, pi].
Definition: phihelper.h:24
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::z0
@ z0
Definition: ParamDefs.h:64
athena.value
value
Definition: athena.py:124
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
Trk::MultiComponentStateModeCalculator::calculateMode
std::array< double, 10 > calculateMode(const MultiComponentState &)
Method to calculate mode with MultiComponentState state as input.
Definition: MultiComponentStateModeCalculator.cxx:386
GSFConstants::maxNumberofStateComponents
constexpr int8_t maxNumberofStateComponents
Note the Gaussian sum approach as describe e.g in " Optimal Filtering" Anderson and Moore "Track Fitt...
Definition: GsfConstants.h:47
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
WritePulseShapeToCool.xhi
xhi
Definition: WritePulseShapeToCool.py:152
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
Trk::theta
@ theta
Definition: ParamDefs.h:66
Trk::MultiComponentState
std::vector< ComponentParameters > MultiComponentState
Definition: ComponentParameters.h:27
MultiComponentStateModeCalculator.h
Trk::ParametersBase
Definition: ParametersBase.h:55
WritePulseShapeToCool.xlo
xlo
Definition: WritePulseShapeToCool.py:133
tolerance
Definition: suep_shower.h:17
Trk::d0
@ d0
Definition: ParamDefs.h:63
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
phihelper.h
Helper for azimuthal angle calculations.
Trk::ComponentParameters
Definition: ComponentParameters.h:22
a
TList * a
Definition: liststreamerinfos.cxx:10
Base_Fragment.width
width
Definition: Sherpa_i/share/common/Base_Fragment.py:59
lwtDev::Component
Component
Definition: NNLayerConfig.h:25
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
PowhegPythia8EvtGen_jetjet.pdf
pdf
Definition: PowhegPythia8EvtGen_jetjet.py:4
extractSporadic.q
list q
Definition: extractSporadic.py:98
GsfConstants.h
Trk::phi
@ phi
Definition: ParamDefs.h:75
python.compressB64.c
def c
Definition: compressB64.py:93