LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
TrackMomentumCalculator.cxx
Go to the documentation of this file.
1 // \file TrackMomentumCalculator.cxx
2 //
3 // \author sowjanyag@phys.ksu.edu
4 
6 #include "cetlib/pow.h"
8 
9 #include <array>
10 #include <cassert>
11 #include <cmath>
12 #include <iostream>
13 #include <limits>
14 #include <string>
15 #include <tuple>
16 
17 #include "Math/Functor.h"
18 #include "Math/GenVector/PositionVector3D.h"
19 #include "Minuit2/Minuit2Minimizer.h"
20 #include "Rtypes.h"
21 #include "TAxis.h"
22 #include "TGraphErrors.h"
23 #include "TMath.h"
24 #include "TMatrixDSymEigen.h"
25 #include "TMatrixDSymfwd.h"
26 #include "TMatrixDfwd.h"
27 #include "TMatrixT.h"
28 #include "TMatrixTSym.h"
29 #include "TPolyLine3D.h"
30 #include "TSpline.h"
31 #include "TVectorDfwd.h"
32 #include "TVectorT.h"
35 
36 using std::cout;
37 using std::endl;
38 
39 namespace {
40 
41  //TODO: Replace for some global variable present in the art file
42  constexpr double LAr_density{1.396};
43  constexpr double rad_length{14.0};
44  constexpr double m_muon{0.1057}; // muon mass
45  constexpr auto range_gramper_cm()
46  {
47  std::array<double, 73> Range_grampercm{
48  {0.9833, 1.36, 1.786, 2.507, 3.321, 4.859, 6.598, 8.512, 10.58, 12.78,
49  15.1, 17.52, 20.04, 25.31, 30.84, 36.59, 42.5, 54.73, 67.32, 86.66,
50  106.3, 139.4, 172.5, 205.6, 238.5, 271.1, 303.5, 335.7, 367.7, 431.0,
51  493.4, 555.2, 616.3, 736.8, 855.2, 1030.0, 1202.0, 1482.0, 1758.0, 2029.0,
52  2297.0, 2562.0, 2825.0, 3085.0, 3343.0, 3854.0, 4359.0, 4859.0, 5354.0, 6333.0,
53  7298.0, 8726.0, 10130.0, 12430.0, 14690.0, 16920.0, 19100.0, 21260.0, 23380.0, 25480.0,
54  27550.0, 31610.0, 35580.0, 39460.0, 43260.0, 50620.0, 57680.0, 67780.0, 77340.0, 92220.0,
55  1.06e+05, 1.188e+05, 1.307e+05}};
56  for (double& value : Range_grampercm) {
57  value /= LAr_density; // convert to cm
58  }
59  return Range_grampercm;
60  }
61 
62  constexpr auto Range_grampercm = range_gramper_cm();
63  constexpr std::array<double, 73> KE_MeV{
64  {10.0, 12.0, 14.0, 17.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0,
65  50.0, 55.0, 60.0, 70.0, 80.0, 90.0, 100.0, 120.0, 140.0, 170.0,
66  200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 700.0,
67  800.0, 900.0, 1000.0, 1200.0, 1400.0, 1700.0, 2000.0, 2500.0, 3000.0, 3500.0,
68  4000.0, 4500.0, 5000.0, 5500.0, 6000.0, 7000.0, 8000.0, 9000.0, 10000.0, 12000.0,
69  14000.0, 17000.0, 20000.0, 25000.0, 30000.0, 35000.0, 40000.0, 45000.0, 50000.0, 55000.0,
70  60000.0, 70000.0, 80000.0, 90000.0, 1e+05, 1.2e+05, 1.4e+05, 1.7e+05, 2e+05, 2.5e+05,
71  3e+05, 3.5e+05, 4e+05}};
72  TGraph const KEvsR{73, Range_grampercm.data(), KE_MeV.data()};
73  TSpline3 const KEvsR_spline3{"KEvsRS", &KEvsR};
74 
75  // Stopping power data from pdg, to be used with MCS
76  std::vector<double> dedx_GeV_per_cm()
77  {
78  // original table (in MeV cm2/g)
79  std::vector<double> dEdx{
80  {5.687, 4.979, 4.461, 3.902, 3.502, 3.042, 2.731, 2.508, 2.34, 2.21, 2.107, 2.023, 1.954,
81  1.848, 1.771, 1.713, 1.67, 1.609, 1.57, 1.536, 1.519, 1.508, 1.51, 1.517, 1.526, 1.537,
82  1.548, 1.559, 1.57, 1.591, 1.61, 1.628, 1.645, 1.675, 1.7, 1.733, 1.761, 1.799, 1.829,
83  1.855, 1.877, 1.897, 1.914, 1.93, 1.944, 1.969, 1.991, 2.01, 2.028, 2.058, 2.084, 2.119,
84  2.149, 2.193, 2.232, 2.269, 2.304, 2.337, 2.369, 2.4, 2.43, 2.49, 2.548, 2.606, 2.663,
85  2.776, 2.888, 3.055, 3.224, 3.498, 3.774, 4.052, 4.332}};
86  for (double& value : dEdx) {
87  value *= LAr_density * 1e-3; // convert to GeV/cm
88  }
89  return dEdx;
90  }
91 
92  const int ndedx = 73;
93  const std::vector<double> dEdx_GeV_per_cm = dedx_GeV_per_cm();
94  const std::vector<double> E_GeV{
95  {0.115695, 0.117697, 0.119693, 0.122694, 0.125695, 0.13069, 0.135694, 0.14069, 0.145714,
96  0.150689, 0.155682, 0.160666, 0.165693, 0.17566, 0.185714, 0.1957, 0.205644, 0.225683,
97  0.245698, 0.275669, 0.305658, 0.355669, 0.405711, 0.45563, 0.505671, 0.555646, 0.605694,
98  0.655676, 0.705661, 0.805664, 0.905689, 1.00557, 1.10606, 1.30529, 1.50571, 1.8061,
99  2.10565, 2.60614, 3.1058, 3.60555, 4.10536, 4.60521, 5.10609, 5.606, 6.10591,
100  7.10579, 8.10569, 9.10561, 10.1106, 12.1105, 14.1104, 17.1103, 20.1103, 25.1102,
101  30.1102, 35.1102, 40.1101, 45.1101, 50.1101, 55.1101, 60.1101, 70.1101, 80.1101,
102  90.1101, 100.1, 120.1, 140.1, 170.1, 200.1, 250.1, 300.1, 350.1,
103  400.1}};
104  TGraph const dEdx_vs_E{ndedx, E_GeV.data(), dEdx_GeV_per_cm.data()};
105  TSpline3 const dEdx_vs_E_spline3{"dEdx_vs_E", &dEdx_vs_E};
106 
107  TVector3 const basex{1, 0, 0};
108  TVector3 const basey{0, 1, 0};
109  TVector3 const basez{0, 0, 1};
110  constexpr double kcal{0.0024}; // Approximation of dE/dx for mip muon in LAr
111 
112  constexpr double MomentumDependentConstant(const double p)
113  {
114  // values measured with MC: https://indico.fnal.gov/event/62283/contributions/299163/attachments/181289/248600/20240910_SimReco_Update_CM.pdf
115  double a = 0.079;
116  double c = 10.435;
117  return (a / (p * p)) + c;
118  }
119  double ComputeExpetecteRMS(const double p, const double red_length)
120  {
121  double beta = std::sqrt(1 - ((m_muon * m_muon) / (p * p + m_muon * m_muon)));
122  return (MomentumDependentConstant(p) / (p * beta)) *
123  (1.0 + 0.038 * TMath::Log(red_length / cet::square(beta))) * std::sqrt(red_length);
124  }
125  class FcnWrapper {
126  public:
127  explicit FcnWrapper(std::vector<double>&& xmeas,
128  std::vector<double>&& ymeas,
129  std::vector<double>&& eymeas,
130  double correction)
131  : xmeas_{std::move(xmeas)}
132  , ymeas_{std::move(ymeas)}
133  , eymeas_{std::move(eymeas)}
134  , correction_{correction}
135  {}
136 
137  double my_mcs_chi2(double const* x) const
138  {
139  double result = 0.0;
140 
141  double p = x[0];
142  double theta0 = x[1];
143 
144  auto const n = xmeas_.size();
145  assert(n == ymeas_.size());
146  assert(n == eymeas_.size());
147 
148  for (std::size_t i = 0; i < n; ++i) {
149  double const xx = xmeas_[i];
150  double const yy = ymeas_[i];
151  double const ey = eymeas_[i];
152 
153  if (std::abs(ey) < std::numeric_limits<double>::epsilon()) {
154  std::cout << " Zero denominator in my_mcs_chi2 ! " << std::endl;
155  return -1;
156  }
157 
158  double const l0 = xx / rad_length;
159  double res1 = 0.0;
160  // Highland formula
161  // Parameters given at Particle Data Group https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-passage-particles-matter.pdf
162  double beta = std::sqrt(1 - ((m_muon * m_muon) / (p * p + m_muon * m_muon)));
163  if (xx > 0 && p > 0)
164  res1 = (13.6 / (p * beta)) * (1.0 + 0.038 * TMath::Log(l0 / cet::square(beta))) *
165  std::sqrt(l0);
166  res1 *= correction_;
167 
168  res1 = std::hypot(res1, theta0);
169 
170  double const diff = yy - res1;
171  result += cet::square(diff / ey);
172  }
173 
174  if (std::isnan(result) || std::isinf(result)) {
175  MF_LOG_DEBUG("TrackMomentumCalculator") << " Is nan in my_mcs_chi2 ! ";
176  return -1;
177  }
178 
179  return result;
180  }
181 
182  private:
183  std::vector<double> const xmeas_;
184  std::vector<double> const ymeas_;
185  std::vector<double> const eymeas_;
186  double const correction_;
187  };
188 
189  class FcnWrapperLLHD {
190  public:
191  explicit FcnWrapperLLHD(std::vector<double>&& dEi,
192  std::vector<double>&& dEj,
193  std::vector<double>&& dthij,
194  std::vector<double>&& ind,
195  std::vector<bool>&& dthij_valid,
196  double stepsize,
197  double correction)
198  : dEi_{std::move(dEi)}
199  , dEj_{std::move(dEj)}
200  , dthij_{std::move(dthij)}
201  , ind_{std::move(ind)}
202  , dthij_valid_{std::move(dthij_valid)}
203  , stepsize_{stepsize}
204  , correction_{correction}
205  {}
206 
207  double my_mcs_llhd(double const* x) const
208  {
209  double result = 0.0;
210 
211  double red_length = stepsize_ / rad_length;
212 
213  double p = x[0];
214  double theta0 = x[1];
215 
216  // Total initial energy of the muon (converting the input "p" into energy with muon mass)
217  double Etot = std::hypot(p, m_muon);
218 
219  double Ei{Etot};
220 
221  double dEi{0};
222  auto const n = dEi_.size(); // number of segments of energy
223 
224  bool addpenality = false;
225  for (std::size_t i = 0; i < n; ++i) {
226  Ei -= dEi;
227  if (Ei >= E_GeV[0]) { //otherwise keep dEi the same (In any case, this is close to p = 0)
228  dEi = dEdx_vs_E_spline3.Eval(Ei) * stepsize_;
229  }
230  else {
231  dEi = dEdx_GeV_per_cm[0] * stepsize_;
232  }
233  if (Ei <= m_muon) {
234  Ei =
235  m_muon +
236  0.010; // Reached zero energy, keep this value constant to not evaluate `nan` in pij nor tH0
237  addpenality = true;
238  }
239 
240  if (dthij_valid_.at(i) == false) continue;
241 
242  // Total momentum of the muon including momentum lost upstream of this segment (converting Eij to momentum)
243  double pij = std::sqrt(Ei * Ei - m_muon * m_muon);
244 
245  // Highland formula
246  // Parameters given at Particle Data Group https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-passage-particles-matter.pdf
247  // Modified with uboone studies: https://iopscience.iop.org/article/10.1088/1748-0221/12/10/P10010
248  double tH0 = ComputeExpetecteRMS(pij, red_length);
249 
250  // The tH0 (theta rms) is calculated for projected angles
251  // If space angles are used instead, tH0 needs to be multiplied by sqrt(2)
252  tH0 *= correction_;
253 
254  double rms_square = -1.0;
255 
256  double prob = 0;
257  double DT = 0;
258  // Computes the rms of angle (no need to evaluate sqrt)
259  rms_square = cet::sum_of_squares(tH0, theta0);
260 
261  DT = dthij_.at(i);
262 
263  // Formula is modified so we don't compute sqrt(rms), use factor in log instead
264  prob = -0.5 * std::log(2.0 * TMath::Pi()) - 0.5 * std::log(rms_square) -
265  0.5 * DT * DT / rms_square;
266 
267  if (addpenality) { prob -= 2 * rms_square; }
268 
269  result = result - 2.0 * prob; // Adds for each segment
270  }
271 
272  return result;
273  }
274 
275  private:
276  std::vector<double> const dEi_;
277  std::vector<double> const dEj_;
278  std::vector<double> const dthij_;
279  std::vector<double> const ind_;
280  std::vector<bool> const dthij_valid_;
281  double const stepsize_;
282  double const correction_;
283  };
284 
285 }
286 
287 namespace trkf {
288 
289  TrackMomentumCalculator::TrackMomentumCalculator(double const min,
290  double const max,
291  double const stepsize,
292  int const angleMethod,
293  int const nsteps)
294  : minLength{min}
295  , maxLength{max}
296  , steps_size{stepsize}
297  , fMCSAngleMethod{static_cast<ScatterAngleMethods>(angleMethod)}
298  , n_steps{nsteps}
299  {
300  for (int i = 1; i <= n_steps; i++) {
301  steps.push_back(steps_size * i);
302  }
303  }
304 
305  double TrackMomentumCalculator::GetTrackMomentum(double trkrange, int pdg) const
306  {
307  /* Muon range-momentum tables from CSDA (Argon density = 1.4 g/cm^3)
308  website:
309  http://pdg.lbl.gov/2012/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_289.pdf
310 
311  CSDA table values:
312  double Range_grampercm[30] = {9.833E-1, 1.786E0, 3.321E0,
313  6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2,
314  2.385E2, 4.934E2, 6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
315  4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4,
316  4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5}; double KE_MeV[30] = {10, 14,
317  20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000,
318  4000, 8000, 10000, 14000, 20000, 30000, 40000, 80000, 100000, 140000,
319  200000, 300000, 400000};
320 
321  Functions below are obtained by fitting polynomial fits to KE_MeV vs
322  Range (cm) graph. A better fit was obtained by splitting the graph into
323  two: Below range<=200cm,a polynomial of power 4 was a good fit; above
324  200cm, a polynomial of power 6 was a good fit
325 
326  Fit errors for future purposes:
327  Below 200cm, Forpoly4 fit: p0 err=1.38533;p1 err=0.209626; p2
328  err=0.00650077; p3 err=6.42207E-5; p4 err=1.94893E-7; Above 200cm,
329  Forpoly6 fit: p0 err=5.24743;p1 err=0.0176229; p2 err=1.6263E-5; p3
330  err=5.9155E-9; p4 err=9.71709E-13; p5 err=7.22381E-17;p6
331  err=1.9709E-21;*/
332 
334  //*********For muon, the calculations are valid up to 1.91E4 cm range
335  //corresponding to a Muon KE of 40 GeV**********//
337 
338  /*Proton range-momentum tables from CSDA (Argon density = 1.4 g/cm^3):
339  website: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html
340 
341  CSDA values:
342  double KE_MeV_P_Nist[31]={10, 15, 20, 30, 40, 80, 100, 150, 200, 250, 300,
343  350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
344  1500, 2000, 2500, 3000, 4000, 5000};
345 
346  double Range_gpercm_P_Nist[31]={1.887E-1,3.823E-1, 6.335E-1, 1.296,
347  2.159, 7.375, 1.092E1, 2.215E1, 3.627E1, 5.282E1, 7.144E1,
348  9.184E1, 1.138E2, 1.370E2, 1.614E2, 1.869E2, 2.132E2, 2.403E2,
349  2.681E2, 2.965E2, 3.254E2, 3.548E2, 3.846E2, 4.148E2, 4.454E2,
350  7.626E2, 1.090E3, 1.418E3, 1.745E3, 2.391E3, 3.022E3};
351 
352  Functions below are obtained by fitting power and polynomial fits to
353  KE_MeV vs Range (cm) graph. A better fit was obtained by splitting the
354  graph into two: Below range<=80cm,a a*(x^b) was a good fit; above 80cm, a
355  polynomial of power 6 was a good fit
356 
357  Fit errors for future purposes:
358  For power function fit: a=0.388873; and b=0.00347075
359  Forpoly6 fit: p0 err=3.49729;p1 err=0.0487859; p2 err=0.000225834; p3
360  err=4.45542E-7; p4 err=4.16428E-10; p5 err=1.81679E-13;p6
361  err=2.96958E-17;*/
362 
364  //*********For proton, the calculations are valid up to 3.022E3 cm range
365  //corresponding to a Muon KE of 5 GeV**********//
367 
368  if (trkrange < 0 || std::isnan(trkrange)) {
369  mf::LogError("TrackMomentumCalculator") << "Invalid track range " << trkrange << " return -1";
370  return -1.;
371  }
372 
373  double KE, Momentum, M;
374  constexpr double Muon_M = m_muon * 1e3, Proton_M = 938.272;
375 
376  if (abs(pdg) == 13) {
377  M = Muon_M;
378  KE = KEvsR_spline3.Eval(trkrange);
379  }
380  else if (abs(pdg) == 2212) {
381  M = Proton_M;
382  if (trkrange > 0 && trkrange <= 80)
383  KE = 29.9317 * std::pow(trkrange, 0.586304);
384  else if (trkrange > 80 && trkrange <= 3.022E3)
385  KE = 149.904 + (3.34146 * trkrange) + (-0.00318856 * trkrange * trkrange) +
386  (4.34587E-6 * trkrange * trkrange * trkrange) +
387  (-3.18146E-9 * trkrange * trkrange * trkrange * trkrange) +
388  (1.17854E-12 * trkrange * trkrange * trkrange * trkrange * trkrange) +
389  (-1.71763E-16 * trkrange * trkrange * trkrange * trkrange * trkrange * trkrange);
390  else
391  KE = -999;
392  }
393  else
394  KE = -999;
395 
396  if (KE < 0)
397  Momentum = -999;
398  else
399  Momentum = std::sqrt((KE * KE) + (2 * M * KE));
400 
401  Momentum = Momentum / 1000;
402 
403  return Momentum;
404  }
405 
406  // Momentum measurement via Multiple Coulomb Scattering (MCS)
407 
408  // author: Leonidas N. Kalousis (July 2015)
409 
410  // email: kalousis@vt.edu
411 
412  // Updated by: Henrique Vieira de Souza (June 2024)
413  // email: hvsouza@apc.in2p3.fr
414 
416  const bool checkValidPoints,
417  const int maxMomentum_MeV,
418  const double min_resolution,
419  const double max_resolution,
420  const bool check_valid_scattered,
421  const double angle_correction)
422  {
423 
424  std::vector<double> recoX;
425  std::vector<double> recoY;
426  std::vector<double> recoZ;
427 
428  int n_points = trk->NumberTrajectoryPoints();
429 
430  for (int i = 0; i < n_points; i++) {
431  if (checkValidPoints && !trk->HasValidPoint(i)) continue;
432  auto const& pos = trk->LocationAtPoint(i);
433  recoX.push_back(pos.X());
434  recoY.push_back(pos.Y());
435  recoZ.push_back(pos.Z());
436  }
437 
438  if (recoX.size() < 2) return -1.0;
439 
440  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
441 
442  double const seg_size{steps_size};
443 
444  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size, check_valid_scattered);
445  if (!segments.has_value()) return -1.0;
446 
447  auto const seg_steps = segments->x.size();
448  if (seg_steps < 2) return -1;
449 
450  double const recoSegmentLength = segments->L.at(seg_steps - 1);
451  if (recoSegmentLength < minLength || recoSegmentLength > maxLength) return -1;
452 
453  std::vector<double> dEi;
454  std::vector<double> dEj;
455  std::vector<double> dthij;
456  std::vector<double> ind;
457  std::vector<bool> dthij_valid = segments->nvalid;
458  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
459  return -1;
460 
461  auto const ndEi = dEi.size();
462  if (ndEi < 1) return -1;
463 
464  double correction = 1.;
465  if (fMCSAngleMethod == kAngleCombined) { correction = std::sqrt(2.); }
466 
467  // Gets energy evaluated by range
468  auto recoL = trk->Length();
469  double const minP = this->GetTrackMomentum(recoL, 13);
470 
471  ROOT::Minuit2::Minuit2Minimizer mP{};
472  FcnWrapperLLHD const wrapper{std::move(dEi),
473  std::move(dEj),
474  std::move(dthij),
475  std::move(ind),
476  std::move(dthij_valid),
477  seg_size,
478  correction};
479  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_llhd(xs); }, 2);
480 
481  mP.SetFunction(FCA);
482 
483  // Start point for resolution
484  // If max_resolution is zero, startpoint is equal min_resolution
485  double startpoint = 2;
486  if (startpoint < min_resolution)
487  startpoint = (max_resolution - min_resolution) / 2.; // prevent wrong values
488  if (max_resolution == 0) startpoint = min_resolution;
489 
490  // Starting energy as double of the energy by range
491  // Assumes that the smallet possible energy is given by 60% with CSDA
492  // Step as 10 %
493  mP.SetLimitedVariable(0, "p_{MCS}", minP * 2, minP * 0.1, minP * 0.6, maxMomentum_MeV / 1.e3);
494  mP.SetLimitedVariable(
495  1, "#delta#theta", startpoint, startpoint / 2., min_resolution, max_resolution);
496  if (max_resolution == 0) { mP.FixVariable(1); }
497  mP.SetMaxFunctionCalls(1.E9);
498  mP.SetMaxIterations(1.E9);
499  mP.SetTolerance(0.01);
500  mP.SetStrategy(2);
501  mP.SetErrorDef(1);
502 
503  bool const mstatus = mP.Minimize();
504 
505  mP.Hesse();
506 
507  const double* pars = mP.X();
508 
509  double const p_mcs = pars[0];
510 
511  return mstatus ? p_mcs : -1.0;
512  }
513 
515  {
516  double const LLHDp = GetMuMultiScatterLLHD3(trk, true);
517  double const LLHDm = GetMuMultiScatterLLHD3(trk, false);
518 
519  if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
520  int const n_points = trk->NumberTrajectoryPoints();
521  return trk->LocationAtPoint<TVector3>(n_points - 1);
522  }
523  else {
524  return trk->LocationAtPoint<TVector3>(0);
525  }
526 
527  // Should never get here
528  return TVector3{};
529  }
530 
532  bool const dir)
533  {
534  std::vector<double> recoX;
535  std::vector<double> recoY;
536  std::vector<double> recoZ;
537 
538  int const n_points = trk->NumberTrajectoryPoints();
539  for (int i = 0; i < n_points; ++i) {
540  auto const index = dir ? i : n_points - 1 - i;
541  auto const& pos = trk->LocationAtPoint(index);
542  recoX.push_back(pos.X());
543  recoY.push_back(pos.Y());
544  recoZ.push_back(pos.Z());
545  }
546 
547  if (recoX.size() < 2) return -1.0;
548 
549  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
550 
551  constexpr double seg_size{5.0};
552  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
553  if (!segments.has_value()) return -1.0;
554 
555  auto const seg_steps = segments->x.size();
556  if (seg_steps < 2) return -1;
557 
558  double const recoSegmentLength = segments->L.at(seg_steps - 1);
559  if (recoSegmentLength < 15.0 || recoSegmentLength > maxLength) return -1;
560 
561  std::vector<double> dEi;
562  std::vector<double> dEj;
563  std::vector<double> dthij;
564  std::vector<double> ind;
565  double angle_correction = 1; // never tested in this function
566  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
567  return -1;
568 
569  auto const recoL = trk->Length();
570  double const p_range = recoL * kcal;
571  double const logL = my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
572 
573  return logL;
574  }
575 
576  int TrackMomentumCalculator::getDeltaThetaij_(std::vector<double>& ei,
577  std::vector<double>& ej,
578  std::vector<double>& th,
579  std::vector<double>& ind,
580  Segments const& segments,
581  double const thick,
582  double const angle_correction) const
583  {
584  auto const& segnx = segments.nx;
585  auto const& segny = segments.ny;
586  auto const& segnz = segments.nz;
587  auto const& segL = segments.L;
588 
589  int const a1 = segnx.size();
590  int const a2 = segny.size();
591  int const a3 = segnz.size();
592 
593  if (a1 != a2 || a1 != a3) {
594  std::cout << " ( Get thij ) Error ! " << std::endl;
595  return -1.0;
596  }
597 
598  int tot = a1;
599 
600  for (int i = 0; i < tot - 1; i++) {
601  double const dx = segnx.at(i);
602  double const dy = segny.at(i);
603  double const dz = segnz.at(i);
604 
605  // Assumes z as propagation angle
606  TVector3 const vec_z{dx, dy, dz};
607  TVector3 vec_x;
608  TVector3 vec_y;
609 
610  double const switcher = basex.Dot(vec_z);
611  if (std::abs(switcher) <= 0.995) {
612  vec_y = vec_z.Cross(basex).Unit();
613  vec_x = vec_y.Cross(vec_z);
614  }
615  else {
616  // cout << " It switched ! Isn't this lovely !!! " << endl; //
617  vec_y = basez.Cross(vec_z).Unit();
618  vec_x = vec_y.Cross(vec_z);
619  }
620 
621  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
622  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
623  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
624 
625  double const here_dx = segnx.at(i + 1);
626  double const here_dy = segny.at(i + 1);
627  double const here_dz = segnz.at(i + 1);
628 
629  TVector3 const here_vec{here_dx, here_dy, here_dz};
630  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
631 
632  double const scx = rot_here.X();
633  double const scy = rot_here.Y();
634  double const scz = rot_here.Z();
635 
636  double const azy = find_angle(scz, scy);
637  double const azx = find_angle(scz, scx);
638 
639  constexpr double ULim = 10000.0; // Avoid huge (wrong) angles
640  constexpr double LLim = -10000.0;
641 
642  double const Li = segL.at(i);
643  double const Lj = segL.at(i + 1);
644 
645  if (azy <= ULim && azy >= LLim) { // safe scatter in the yz plane
646 
647  if (azx <= ULim && azx >= LLim) { // safe scatter in the za plane
648  ei.push_back(Li); // Energy deposited at i
649  ej.push_back(Lj); // Energy deposited at i+1
650  if (fMCSAngleMethod == kAnglezx) {
651  th.push_back(azx); // scattered angle z-x
652  }
653  else if (fMCSAngleMethod == kAnglezy) {
654  th.push_back(azy);
655  }
656  else if (fMCSAngleMethod == kAngleCombined) {
657  th.push_back(std::hypot(azx, azy) /
658  angle_correction); // space angle (applying correction)
659  }
660  }
661  }
662  }
663 
664  return 0;
665  }
666 
668  const bool checkValidPoints,
669  const int maxMomentum_MeV,
670  const double min_resolution,
671  const double max_resolution)
672  {
673  std::vector<double> recoX;
674  std::vector<double> recoY;
675  std::vector<double> recoZ;
676 
677  int n_points = trk->NumberTrajectoryPoints();
678 
679  for (int i = 0; i < n_points; i++) {
680  if (checkValidPoints && !trk->HasValidPoint(i)) continue;
681  auto const& pos = trk->LocationAtPoint(i);
682  recoX.push_back(pos.X());
683  recoY.push_back(pos.Y());
684  recoZ.push_back(pos.Z());
685  }
686 
687  if (recoX.size() < 2) return -1.0;
688 
689  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
690 
691  double const seg_size{steps_size};
692  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
693  if (!segments.has_value()) return -1.0;
694 
695  auto const seg_steps = segments->x.size();
696  if (seg_steps < 2) return -1;
697 
698  double const recoSegmentLength = segments->L.at(seg_steps - 1);
699  if (recoSegmentLength < minLength || recoSegmentLength > maxLength) return -1;
700 
701  double ymax = -999.0;
702  double ymin = +999.0;
703 
704  std::vector<double> xmeas;
705  std::vector<double> ymeas;
706  std::vector<double> eymeas;
707  xmeas.reserve(n_steps);
708  ymeas.reserve(n_steps);
709  eymeas.reserve(n_steps);
710  for (int j = 0; j < n_steps; j++) {
711  double const trial = steps.at(j);
712  // computes the rms by groups of trial, if seg_size was chosen as 10, trials will be 10, 20, etc.. until 10 * n_steps
713  auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
714 
715  if (std::isnan(mean) || std::isinf(mean)) {
716  mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
717  continue;
718  }
719  if (std::isnan(rms) || std::isinf(rms)) {
720  mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
721  continue;
722  }
723  if (std::isnan(rmse) || std::isinf(rmse)) {
724  mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
725  continue;
726  }
727 
728  if (mean == -1 && rms == -1 && rmse == -1) continue;
729  xmeas.push_back(trial); // x values are different steps length, ex: 10, 20, 30 cm
730  ymeas.push_back(rms); // y values are the RMS of the scattered angle for each step defined
731  eymeas.push_back(
732  std::hypot(rmse, 0.05 * rms)); // <--- conservative syst. error to fix chi^{2} behaviour !!!
733 
734  if (ymin > rms) ymin = rms;
735  if (ymax < rms) ymax = rms;
736  }
737 
738  assert(xmeas.size() == ymeas.size());
739  assert(xmeas.size() == eymeas.size());
740  if (xmeas.empty()) { return -1.0; }
741 
742  TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
743 
744  gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
745  "thickness in cm; (#Delta#theta)_{rms} in mrad");
746 
747  gr_meas.SetLineColor(kBlack);
748  gr_meas.SetMarkerColor(kBlack);
749  gr_meas.SetMarkerStyle(20);
750  gr_meas.SetMarkerSize(1.2);
751 
752  gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0), steps.at(n_steps - 1) + steps.at(0));
753  gr_meas.SetMinimum(0.0);
754  gr_meas.SetMaximum(1.80 * ymax);
755 
756  double correction = 1.;
757  if (fMCSAngleMethod == kAngleCombined) { correction = std::sqrt(2.); }
758 
759  ROOT::Minuit2::Minuit2Minimizer mP{};
760  FcnWrapper const wrapper{std::move(xmeas), std::move(ymeas), std::move(eymeas), correction};
761  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
762 
763  // Start point for resolution
764  double startpoint = 2;
765  if (startpoint < min_resolution) startpoint = (max_resolution - min_resolution) / 2.;
766  if (max_resolution == 0) startpoint = min_resolution;
767 
768  mP.SetFunction(FCA);
769  mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, maxMomentum_MeV / 1.e3);
770  mP.SetLimitedVariable(
771  1, "#delta#theta", startpoint, startpoint / 2., min_resolution, max_resolution);
772  if (max_resolution == 0) { mP.FixVariable(1); }
773  mP.SetMaxFunctionCalls(1.E9);
774  mP.SetMaxIterations(1.E9);
775  mP.SetTolerance(0.01);
776  mP.SetStrategy(2);
777  mP.SetErrorDef(1);
778 
779  bool const mstatus = mP.Minimize();
780 
781  mP.Hesse();
782 
783  const double* pars = mP.X();
784 
785  auto const recoL = trk->Length();
786  double const deltap = (recoL * kcal) / 2.0;
787 
788  double const p_mcs = pars[0] + deltap;
789 
790  return mstatus ? p_mcs : -1.0;
791  }
792 
793  bool TrackMomentumCalculator::plotRecoTracks_(std::vector<double> const& xxx,
794  std::vector<double> const& yyy,
795  std::vector<double> const& zzz)
796  {
797  auto const n = xxx.size();
798  auto const y_size = yyy.size();
799  auto const z_size = zzz.size();
800 
801  if (n != y_size || n != z_size) {
802  cout << " ( Get reco tacks ) Error ! " << endl;
803  return false;
804  }
805 
806  // Here, we perform a const-cast to double* because, sadly,
807  // TPolyLine3D requires a pointer to a non-const object. We will
808  // trust that ROOT does not mess around with the underlying data.
809  auto xs = const_cast<double*>(xxx.data());
810  auto ys = const_cast<double*>(yyy.data());
811  auto zs = const_cast<double*>(zzz.data());
812 
813  auto const narrowed_size = static_cast<int>(n); // ROOT requires a signed integral type
814  delete gr_reco_xyz;
815  gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
816  gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
817  gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
818  gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
819 
820  return true;
821  }
822 
823  // Compute the deviation in `segx, ...` of the segment stores at `segnx, ...`
824  // vx, vy, vz are used and cleared afterwards
825  // The `maximum fluctuation` is given by the direction of the give by the Principal Component Analysis (PCA)
826  void TrackMomentumCalculator::compute_max_fluctuation_vector(const std::vector<double> segx,
827  const std::vector<double> segy,
828  const std::vector<double> segz,
829  std::vector<double>& segnx,
830  std::vector<double>& segny,
831  std::vector<double>& segnz,
832  std::vector<bool>& segn_isvalid,
833  std::vector<double>& vx,
834  std::vector<double>& vy,
835  std::vector<double>& vz,
836  bool check_valid_scattered)
837  {
838  auto const na = vx.size();
839 
840  double sumx = 0.0;
841  double sumy = 0.0;
842  double sumz = 0.0;
843 
844  bool isvalid = true;
845  // In case vx, vy, etc have 3 points, probably two are just "linear"
846  // interpolations. In this case, the angle of scattering will be zero.
847  if (check_valid_scattered && na <= 2) { isvalid = false; }
848 
849  // computes the average in x, y, z
850  for (std::size_t i = 0; i < na; ++i) {
851  sumx += vx.at(i);
852  sumy += vy.at(i);
853  sumz += vz.at(i);
854  }
855 
856  sumx /= na;
857  sumy /= na;
858  sumz /= na;
859 
860  std::vector<double> mx;
861  std::vector<double> my;
862  std::vector<double> mz;
863 
864  TMatrixDSym m(3);
865 
866  // Computes the Covariance matrix (Principal Component Analysis (PCA)
867  for (std::size_t i = 0; i < na; ++i) {
868  double const xxw1 = vx.at(i);
869  double const yyw1 = vy.at(i);
870  double const zzw1 = vz.at(i);
871 
872  mx.push_back(xxw1 - sumx);
873  my.push_back(yyw1 - sumy);
874  mz.push_back(zzw1 - sumz);
875 
876  double const xxw0 = mx.at(i);
877  double const yyw0 = my.at(i);
878  double const zzw0 = mz.at(i);
879 
880  m(0, 0) += xxw0 * xxw0 / na;
881  m(0, 1) += xxw0 * yyw0 / na;
882  m(0, 2) += xxw0 * zzw0 / na;
883 
884  m(1, 0) += yyw0 * xxw0 / na;
885  m(1, 1) += yyw0 * yyw0 / na;
886  m(1, 2) += yyw0 * zzw0 / na;
887 
888  m(2, 0) += zzw0 * xxw0 / na;
889  m(2, 1) += zzw0 * yyw0 / na;
890  m(2, 2) += zzw0 * zzw0 / na;
891  }
892 
893  TMatrixDSymEigen me(m);
894 
895  // retrieve eigenvalues and vectors
896  TVectorD eigenval = me.GetEigenValues();
897  TMatrixD eigenvec = me.GetEigenVectors();
898 
899  double max1 = -666.0;
900 
901  double ind1 = 0;
902 
903  // get maximum eingevalue
904  for (int i = 0; i < 3; ++i) {
905  double const p1 = eigenval(i);
906 
907  if (p1 > max1) {
908  max1 = p1;
909  ind1 = i;
910  }
911  }
912 
913  // set the `direction` vector that points to the maximum fluctuation
914  double ax = eigenvec(0, ind1);
915  double ay = eigenvec(1, ind1);
916  double az = eigenvec(2, ind1);
917 
918  if (n_seg > 1) {
919  // for x, y and z, check if the last point is bigger then the previous point.
920  // Ensures that the computed fluctation follows the trend of the track
921  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
922  ax = std::abs(ax);
923  else
924  ax = -1.0 * std::abs(ax);
925 
926  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
927  ay = std::abs(ay);
928  else
929  ay = -1.0 * std::abs(ay);
930 
931  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
932  az = std::abs(az);
933  else
934  az = -1.0 * std::abs(az);
935 
936  segnx.push_back(ax);
937  segny.push_back(ay);
938  segnz.push_back(az);
939  segn_isvalid.push_back(isvalid);
940  }
941 
942  // clear the vectors
943  vx.clear();
944  vy.clear();
945  vz.clear();
946  }
947 
948  /* This function will group each point of the track inside segments with
949  * fixed size `seg_size`.
950  * It returns computed new points `segx, segy, segz` separated by `seg_size`,
951  * the deviation `segnx, ...` between this points and the distance `segL` in
952  * steps of `seg_size` between these points
953  */
954  std::optional<TrackMomentumCalculator::Segments> TrackMomentumCalculator::getSegTracks_(
955  std::vector<double> const& xxx,
956  std::vector<double> const& yyy,
957  std::vector<double> const& zzz,
958  double const seg_size,
959  bool const check_valid_scattered)
960  {
961  double stag = 0.0;
962 
963  int a1 = xxx.size();
964  int a2 = yyy.size();
965  int a3 = zzz.size();
966 
967  if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
968  cout << " ( Digitize reco tacks ) Error ! " << endl;
969  return std::nullopt;
970  }
971 
972  int const stopper = seg_stop / seg_size;
973 
974  // values to be filled and returned
975  std::vector<double> segx, segnx;
976  std::vector<double> segy, segny;
977  std::vector<double> segz, segnz;
978  std::vector<bool> segn_isvalid;
979  std::vector<double> segL;
980 
981  n_seg = 0;
982 
983  double x0{};
984  double y0{};
985  double z0{};
986 
987  double x00 = xxx.at(0);
988  double y00 = yyy.at(0);
989  double z00 = zzz.at(0);
990 
991  int indC = 0;
992 
993  // These vectors will keep the points inside reach segments
994  // They are cleared inside the function `compute_max_fluctuation_vector`
995  // each time it finishes with a segment
996  std::vector<double> vx;
997  std::vector<double> vy;
998  std::vector<double> vz;
999 
1000  for (int i = 0; i < a1; i++) {
1001  x0 = xxx.at(i);
1002  y0 = yyy.at(i);
1003  z0 = zzz.at(i);
1004 
1005  double const RR0 = std::hypot(x00 - x0, y00 - y0, z00 - z0);
1006 
1007  if (RR0 >= stag) { // stag is aways set to zero, this is always true
1008 
1009  segx.push_back(x0);
1010  segy.push_back(y0);
1011  segz.push_back(z0);
1012 
1013  segL.push_back(0);
1014 
1015  // TGraph
1016  x_seg[n_seg] = x0;
1017  y_seg[n_seg] = y0;
1018  z_seg[n_seg] = z0;
1019 
1020  n_seg++;
1021 
1022  vx.push_back(x0);
1023  vy.push_back(y0);
1024  vz.push_back(z0);
1025 
1026  indC = i + 1;
1027 
1028  break;
1029  }
1030  }
1031 
1032  for (int i = indC; i < a1 - 1; i++) { // starting at second point (i=1) if stag set to zero
1033  // current point
1034  double const x1 = xxx.at(i);
1035  double const y1 = yyy.at(i);
1036  double const z1 = zzz.at(i);
1037 
1038  // distante from previous point
1039  double const dr1 = std::hypot(x1 - x0, y1 - y0, z1 - z0);
1040 
1041  // next point
1042  double const x2 = xxx.at(i + 1);
1043  double const y2 = yyy.at(i + 1);
1044  double const z2 = zzz.at(i + 1);
1045 
1046  // distant of next point to previous point
1047  double const dr2 = std::hypot(x2 - x0, y2 - y0, z2 - z0);
1048 
1049  if (dr1 < seg_size) {
1050  vx.push_back(x1);
1051  vy.push_back(y1);
1052  vz.push_back(z1);
1053  }
1054 
1055  /* If current point is inside segment length w.r.t. the first point of
1056  * the segment (x0,y0,z0), but the next point is outsize: create a new
1057  * point in between (x1,y1,z1) and (x2,y2,z2) in which will be exacly at
1058  * the segment length (w.r.t. to the first point) This is done using the
1059  * cosines law for a given factor `t` times dr (x2-x1, ...), and so:
1060  *
1061  * (ds)^2 = (dr1)^2 + (t*dr)^2 + 2*dot_product(dr1, t*dr)
1062  * Using cos(180-theta) = -cos(theta))
1063  *
1064  * Solve for `t` the second degree equation: t^2 + beta*t + gamma = 0
1065  */
1066  if (dr1 <= seg_size && dr2 > seg_size) {
1067  double const dx = x2 - x1;
1068  double const dy = y2 - y1;
1069  double const dz = z2 - z1;
1070  double const dr = std::hypot(dx, dy, dz);
1071 
1072  if (dr == 0) {
1073  cout << " ( Zero ) Error ! " << endl;
1074  return std::nullopt;
1075  }
1076 
1077  double const beta = 2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
1078 
1079  double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
1080  double const delta = beta * beta - 4.0 * gamma;
1081 
1082  if (delta < 0.0) {
1083  cout << " ( Discriminant ) Error ! " << endl;
1084  return std::nullopt;
1085  }
1086 
1087  // solves for t
1088  double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
1089  double const t = lysi1;
1090 
1091  // find next points in that will exactly at the segment length from x0
1092  double const xp = x1 + t * dx;
1093  double const yp = y1 + t * dy;
1094  double const zp = z1 + t * dz;
1095 
1096  // Add points to be returned
1097  segx.push_back(xp);
1098  segy.push_back(yp);
1099  segz.push_back(zp);
1100 
1101  segL.push_back(n_seg * seg_size);
1102 
1103  // for TGraph
1104  x_seg[n_seg] = xp;
1105  y_seg[n_seg] = yp;
1106  z_seg[n_seg] = zp;
1107 
1108  n_seg++;
1109 
1110  // This are the new `x0` points (for next segment)
1111  x0 = xp;
1112  y0 = yp;
1113  z0 = zp;
1114 
1115  // Add points to segment
1116  vx.push_back(x0);
1117  vy.push_back(y0);
1118  vz.push_back(z0);
1119 
1120  if (n_seg <= 1) // This should never happen
1121  return std::nullopt;
1122 
1123  // Now, compute the deviation in `segx, ...` of the segment
1124  // vx, vy, vz are used and cleared afterwards
1126  segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
1127 
1128  // Starting over
1129  vx.push_back(x0);
1130  vy.push_back(y0);
1131  vz.push_back(z0);
1132  }
1133  else if (dr1 > seg_size) { // in this case, just interpolate until reach `seg_size`
1134 
1135  // Rolling `i` back for the next iteration.
1136  // Because the current point does not belong to this segment
1137  i = (i - 1);
1138 
1139  double const dx = x1 - x0;
1140  double const dy = y1 - y0;
1141  double const dz = z1 - z0;
1142  double const dr = std::hypot(dx, dy, dz);
1143 
1144  if (dr == 0) {
1145  cout << " ( Zero ) Error ! " << endl;
1146  return std::nullopt;
1147  }
1148 
1149  // computes the point by simple interpolation
1150  double const t = seg_size / dr;
1151  double const xp = x0 + t * dx;
1152  double const yp = y0 + t * dy;
1153  double const zp = z0 + t * dz;
1154 
1155  segx.push_back(xp);
1156  segy.push_back(yp);
1157  segz.push_back(zp);
1158  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
1159 
1160  // for TGraph
1161  x_seg[n_seg] = xp;
1162  y_seg[n_seg] = yp;
1163  z_seg[n_seg] = zp;
1164 
1165  n_seg++;
1166 
1167  x0 = xp;
1168  y0 = yp;
1169  z0 = zp;
1170 
1171  vx.push_back(x0);
1172  vy.push_back(y0);
1173  vz.push_back(z0);
1174 
1175  if (n_seg <= 1) // This should never happen
1176  return std::nullopt;
1177 
1178  // Now, compute the deviation in `segx, ...` of the segment
1179  // vx, vy, vz are used and cleared afterwards
1181  segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
1182 
1183  // vectors are cleared in previous step
1184  vx.push_back(x0);
1185  vy.push_back(y0);
1186  vz.push_back(z0);
1187  }
1188 
1189  if (n_seg >= (stopper + 1.0) && seg_stop != -1) break;
1190  }
1191 
1192  delete gr_seg_xyz;
1193  gr_seg_xyz = new TPolyLine3D{n_seg, z_seg, x_seg, y_seg};
1194  gr_seg_yz = TGraph{n_seg, z_seg, y_seg};
1195  gr_seg_xz = TGraph{n_seg, z_seg, x_seg};
1196  gr_seg_xy = TGraph{n_seg, x_seg, y_seg};
1197 
1198  return std::make_optional<Segments>(
1199  Segments{segx, segnx, segy, segny, segz, segnz, segL, segn_isvalid});
1200  }
1201 
1202  /* Computes the rms by groups of `thick`
1203  *
1204  */
1205  std::tuple<double, double, double> TrackMomentumCalculator::getDeltaThetaRMS_(
1206  Segments const& segments,
1207  double const thick) const
1208  {
1209  auto const& segnx = segments.nx;
1210  auto const& segny = segments.ny;
1211  auto const& segnz = segments.nz;
1212  auto const& segL = segments.L;
1213 
1214  int const a1 = segnx.size();
1215  int const a2 = segny.size();
1216  int const a3 = segnz.size();
1217 
1218  if (a1 != a2 || a1 != a3) {
1219  cout << " ( Get RMS ) Error ! " << endl;
1220  return std::make_tuple(0., 0., 0.);
1221  }
1222 
1223  int const tot = a1;
1224 
1225  std::vector<double> buf0;
1226 
1227  for (int i = 0; i < tot; i++) {
1228  double const dx = segnx.at(i);
1229  double const dy = segny.at(i);
1230  double const dz = segnz.at(i);
1231 
1232  TVector3 const vec_z{dx, dy, dz};
1233  TVector3 vec_x;
1234  TVector3 vec_y;
1235 
1236  double const switcher = basex.Dot(vec_z);
1237 
1238  if (switcher <= 0.995) {
1239  vec_y = vec_z.Cross(basex).Unit();
1240  vec_x = vec_y.Cross(vec_z);
1241  }
1242  else {
1243  // cout << " It switched ! Isn't this lovely !!! " << endl;
1244  vec_y = basez.Cross(vec_z).Unit();
1245  vec_x = vec_y.Cross(vec_z);
1246  }
1247 
1248  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
1249  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
1250  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
1251 
1252  double const refL = segL.at(i);
1253 
1254  for (int j = i; j < tot; j++) {
1255  double const L1 = segL.at(j);
1256 
1257  double const dz1 = L1 - refL;
1258 
1259  if (dz1 >= thick) {
1260  double const here_dx = segnx.at(j);
1261  double const here_dy = segny.at(j);
1262  double const here_dz = segnz.at(j);
1263 
1264  TVector3 const here_vec{here_dx, here_dy, here_dz};
1265  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
1266 
1267  double const scx = rot_here.X();
1268  double const scy = rot_here.Y();
1269  double const scz = rot_here.Z();
1270 
1271  double const azy = find_angle(scz, scy);
1272  double const azx = find_angle(scz, scx);
1273 
1274  constexpr double ULim = 10000.0; // Avoid huge (wrong) angles
1275  constexpr double LLim = -10000.0;
1276 
1277  if (azy <= ULim && azy >= LLim) { // safe scatter in the yz plane
1278 
1279  if (azx <= ULim && azx >= LLim) { // safe scatter in the za plane
1280  if (fMCSAngleMethod == kAnglezx) {
1281  buf0.push_back(azx); // scattered angle z-x
1282  }
1283  else if (fMCSAngleMethod == kAnglezy) {
1284  buf0.push_back(azy);
1285  }
1286  else if (fMCSAngleMethod == kAngleCombined) {
1287  buf0.push_back(std::hypot(azx, azy));
1288  }
1289  }
1290  }
1291  break; // of course !
1292  }
1293  }
1294  }
1295 
1296  int const nmeas = buf0.size();
1297  double nnn = 0.0;
1298 
1299  double mean = 0.0;
1300  double rms = 0.0;
1301  double rmse = 0.0;
1302 
1303  for (int i = 0; i < nmeas; i++) {
1304  mean += buf0.at(i);
1305  nnn++;
1306  }
1307 
1308  mean = mean / nnn;
1309 
1310  for (int i = 0; i < nmeas; i++)
1311  rms += ((buf0.at(i)) * (buf0.at(i)));
1312 
1313  rms = rms / (nnn);
1314  rms = std::sqrt(rms);
1315  rmse = rms / std::sqrt(2.0 * tot);
1316 
1317  double rms1 = rms;
1318 
1319  rms = 0.0;
1320 
1321  double ntot1 = 0.0;
1322  double const lev1 = 2.50;
1323 
1324  for (int i = 0; i < nmeas; i++) {
1325  double const amp = buf0.at(i);
1326  if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
1327  ++ntot1;
1328  rms += amp * amp;
1329  }
1330  }
1331 
1332  rms = rms / (ntot1);
1333  rms = std::sqrt(rms);
1334  rmse = rms / std::sqrt(2.0 * ntot1);
1335  return std::make_tuple(mean, rms, rmse);
1336  }
1337 
1338  double TrackMomentumCalculator::find_angle(double vz, double vy) const
1339  {
1340  double thetayz = -999.0;
1341 
1342  if (vz > 0 && vy > 0) {
1343  double ratio = std::abs(vy / vz);
1344  thetayz = std::atan(ratio);
1345  }
1346 
1347  else if (vz < 0 && vy > 0) {
1348  double ratio = std::abs(vy / vz);
1349  thetayz = std::atan(ratio);
1350  thetayz = TMath::Pi() - thetayz;
1351  }
1352 
1353  else if (vz < 0 && vy < 0) {
1354  double ratio = std::abs(vy / vz);
1355  thetayz = std::atan(ratio);
1356  thetayz = thetayz + TMath::Pi();
1357  }
1358 
1359  else if (vz > 0 && vy < 0) {
1360  double ratio = std::abs(vy / vz);
1361  thetayz = std::atan(ratio);
1362  thetayz = 2.0 * TMath::Pi() - thetayz;
1363  }
1364 
1365  else if (vz == 0 && vy > 0) {
1366  thetayz = TMath::Pi() / 2.0;
1367  }
1368 
1369  else if (vz == 0 && vy < 0) {
1370  thetayz = 3.0 * TMath::Pi() / 2.0;
1371  }
1372  else if (vz > 0 && vy == 0) {
1373  thetayz = 0;
1374  }
1375  else if (vz < 0 && vy == 0) {
1376  thetayz = TMath::Pi();
1377  }
1378 
1379  if (thetayz > TMath::Pi()) { thetayz = thetayz - 2.0 * TMath::Pi(); }
1380 
1381  return 1000.0 * thetayz;
1382  }
1383 
1384  double TrackMomentumCalculator::my_g(double xx, double Q, double s) const
1385  {
1386  if (s == 0.) {
1387  cout << " Error : The code tries to divide by zero ! " << endl;
1388  return 0.;
1389  }
1390 
1391  double const arg = (xx - Q) / s;
1392  double const result = -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1393 
1394  if (std::isnan(result) || std::isinf(result)) {
1395  cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
1396  }
1397 
1398  return result;
1399  }
1400 
1401  double TrackMomentumCalculator::my_mcs_llhd(std::vector<double> const& dEi,
1402  std::vector<double> const& dEj,
1403  std::vector<double> const& dthij,
1404  std::vector<double> const& ind,
1405  double const x0,
1406  double const x1) const
1407  {
1408  double p = x0;
1409  double theta0x = x1;
1410 
1411  double result = 0.0;
1412  double nnn1 = dEi.size(); // number of segments of energy
1413 
1414  double red_length = (steps_size) / rad_length;
1415  double addth = 0;
1416 
1417  for (int i = 0; i < nnn1; i++) {
1418  double Ei = p - dEi.at(i); // Estimated energy at point i
1419  double Ej = p - dEj.at(i); // Estimated enery at point j
1420 
1421  // If the momentum p choosen allows that the muon stopped inside, add 1 rad to the change in scatter angle (as the particle stops)
1422  if (Ei > 0 && Ej < 0) addth = 3.14 * 1000.0;
1423 
1424  Ei = std::abs(Ei);
1425  Ej = std::abs(Ej);
1426 
1427  // Highland formula
1428  // Parameters given at Particle Data Group https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-passage-particles-matter.pdf
1429  double tH0 =
1430  (13.6 / std::sqrt(Ei * Ej)) * (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1431 
1432  double rms = -1.0;
1433 
1434  if (ind.at(i) == 1) {
1435  // Computes the rms of angle
1436  rms = std::hypot(tH0, theta0x);
1437 
1438  double const DT = dthij.at(i) + addth;
1439  double const prob = my_g(DT, 0.0, rms); // Computes log likelihood
1440 
1441  result = result - 2.0 * prob; // Adds for each segment
1442  }
1443  }
1444 
1445  if (std::isnan(result) || std::isinf(result)) {
1446  std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
1447  }
1448  return result;
1449  }
1450 
1451 } // namespace track
Float_t x
Definition: compare.C:6
double GetMomentumMultiScatterLLHD(art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const double min_resolution=0.001, const double max_resolution=800, const bool check_valid_scattered=false, const double angle_correction=0.757)
Calculate muon momentum (GeV) using multiple coulomb scattering by log likelihood.
Double_t xx
Definition: macro.C:12
double GetMomentumMultiScatterChi2(art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const double min_resolution=0, const double max_resolution=45)
Calculate muon momentum (GeV) using multiple coulomb scattering. Chi2 minimization of the Highland fo...
double my_g(double xx, double Q, double s) const
chi square minizer using Minuit2, it will minize (xx-Q)/s
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:145
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
constexpr auto abs(T v)
Returns the absolute value of the argument.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::optional< Segments > getSegTracks_(std::vector< double > const &xxx, std::vector< double > const &yyy, std::vector< double > const &zzz, double seg_size, bool check_valid_scattered=false)
Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/...
Float_t y2[n_points_geant4]
Definition: compare.C:26
int getDeltaThetaij_(std::vector< double > &ei, std::vector< double > &ej, std::vector< double > &th, std::vector< double > &ind, Segments const &segments, double thick, double const angle_correction) const
Gets the scatterd angle for all the segments.
#define a2
#define a3
Use scattered angle z-x (z is along the particle&#39;s direction)
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
double find_angle(double vz, double vy) const
Gets angle between two vy and vz.
Float_t E
Definition: plot.C:20
Float_t thick
Definition: plot.C:23
double my_mcs_llhd(std::vector< double > const &dEi, std::vector< double > const &dEj, std::vector< double > const &dthij, std::vector< double > const &ind, double x0, double x1) const
Minimizer of log likelihood for scattered angle.
bool plotRecoTracks_(std::vector< double > const &xxx, std::vector< double > const &yyy, std::vector< double > const &zzz)
Struct to store segments. x, y and z are the 3D points of the segment nx, ny, nz forms the vector tha...
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Gets the scattered angle RMS for a all segments.
Use space angle: sqrt( zx^2 + zy^2 )/sqrt(2)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2671
double value
Definition: spectrum.C:18
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
void compute_max_fluctuation_vector(const std::vector< double > segx, const std::vector< double > segy, const std::vector< double > segz, std::vector< double > &segnx, std::vector< double > &segny, std::vector< double > &segnz, std::vector< bool > &segn_isvalid, std::vector< double > &vx, std::vector< double > &vy, std::vector< double > &vz, bool check_valid_scattered)
Computes the vector with most scattering inside a segment with size steps_size.
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Provides recob::Track data product.
TDirectory * dir
Definition: macro.C:5
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
#define MF_LOG_DEBUG(id)
TVector3 GetMultiScatterStartingPoint(art::Ptr< recob::Track > const &trk)
Char_t n[5]
Float_t x2[n_points_geant4]
Definition: compare.C:26
Float_t e
Definition: plot.C:35
double GetTrackMomentum(double trkrange, int pdg) const
#define a1