LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
TrackMomentumCalculator.cxx
Go to the documentation of this file.
1 // \file TrackMomentumCalculator.cxx
2 //
3 // \author sowjanyag@phys.ksu.edu
4 
5 #include "cetlib/pow.h"
6 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
7 #include "messagefacility/MessageLogger/MessageLogger.h"
8 
9 #include <array>
10 #include <cassert>
11 
12 using std::cout;
13 using std::endl;
14 
15 namespace {
16 
17  constexpr auto range_gramper_cm()
18  {
19  std::array<float, 29> Range_grampercm{{
20  9.833E-1, 1.786E0, 3.321E0, 6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1,
21  1.063E2, 1.725E2, 2.385E2, 4.934E2, 6.163E2, 8.552E2, 1.202E3, 1.758E3,
22  2.297E3, 4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4,
23  4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5}};
24  for (float& value : Range_grampercm) {
25  value /= 1.396; // convert to cm
26  }
27  return Range_grampercm;
28  }
29 
30  constexpr auto Range_grampercm = range_gramper_cm();
31  constexpr std::array<float, 29> KE_MeV{{
32  10, 14, 20, 30, 40, 80, 100, 140, 200, 300,
33  400, 800, 1000, 1400, 2000, 3000, 4000, 8000, 10000, 14000,
34  20000, 30000, 40000, 80000, 100000, 140000, 200000, 300000, 400000}};
35  TGraph const KEvsR{29, Range_grampercm.data(), KE_MeV.data()};
36  TSpline3 const KEvsR_spline3{"KEvsRS", &KEvsR};
37 
38 
39  TVector3 const basex{1, 0, 0};
40  TVector3 const basey{0, 1, 0};
41  TVector3 const basez{0, 0, 1};
42  constexpr float kcal{0.0024};
43 
44  class FcnWrapper {
45  public:
46  explicit FcnWrapper(std::vector<double>&& xmeas,
47  std::vector<double>&& ymeas,
48  std::vector<double>&& eymeas)
49  : xmeas_{xmeas}
50  , ymeas_{ymeas}
51  , eymeas_{eymeas}
52  {}
53 
54  double
55  my_mcs_chi2(double const* x) const
56  {
57  double result = 0.0;
58 
59  double p = x[0];
60  double theta0 = x[1];
61 
62  auto const n = xmeas_.size();
63  assert(n == ymeas_.size());
64  assert(n == eymeas_.size());
65 
66  for (std::size_t i = 0; i < n; ++i) {
67  double const xx = xmeas_[i];
68  double const yy = ymeas_[i];
69  double const ey = eymeas_[i];
70 
71  if (std::abs(ey) < std::numeric_limits<double>::epsilon()) {
72  std::cout << " Zero denominator in my_mcs_chi2 ! " << std::endl;
73  return -1;
74  }
75 
76  constexpr double rad_length{14.0};
77  double const l0 = xx / rad_length;
78  double res1 = 0.0;
79 
80  if (xx > 0 && p > 0)
81  res1 = (13.6 / p) * std::sqrt(l0) * (1.0 + 0.038 * std::log(l0));
82 
83  res1 = std::sqrt(res1 * res1 + theta0 * theta0);
84 
85  double const diff = yy - res1;
86  result += cet::square(diff / ey);
87  }
88 
89  result += 2.0 / (4.6) * theta0; // *std::log( 1.0/14.0 );
90 
91  if (std::isnan(result) || std::isinf(result)) {
92  LOG_DEBUG("TrackMomentumCalculator") << " Is nan in my_mcs_chi2 ! ";
93  return -1;
94  }
95 
96  return result;
97  }
98 
99  private:
100  std::vector<double> const xmeas_;
101  std::vector<double> const ymeas_;
102  std::vector<double> const eymeas_;
103  };
104 
105 }
106 
107 namespace trkf {
108 
109  TrackMomentumCalculator::TrackMomentumCalculator(double const min,
110  double const max)
111  : minLength{min}
112  , maxLength{max}
113  {
114  for (int i = 1; i <= n_steps; i++) {
115  steps.push_back(steps_size * i);
116  }
117  }
118 
119  double
120  TrackMomentumCalculator::GetTrackMomentum(double trkrange, int pdg) const
121  {
122  /* Muon range-momentum tables from CSDA (Argon density = 1.4 g/cm^3)
123  website:
124  http://pdg.lbl.gov/2012/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_289.pdf
125 
126  CSDA table values:
127  float Range_grampercm[30] = {9.833E-1, 1.786E0, 3.321E0,
128  6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2,
129  2.385E2, 4.934E2, 6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
130  4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4,
131  4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5}; float KE_MeV[30] = {10, 14,
132  20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000,
133  4000, 8000, 10000, 14000, 20000, 30000, 40000, 80000, 100000, 140000,
134  200000, 300000, 400000};
135 
136  Functions below are obtained by fitting polynomial fits to KE_MeV vs
137  Range (cm) graph. A better fit was obtained by splitting the graph into
138  two: Below range<=200cm,a polynomial of power 4 was a good fit; above
139  200cm, a polynomial of power 6 was a good fit
140 
141  Fit errors for future purposes:
142  Below 200cm, Forpoly4 fit: p0 err=1.38533;p1 err=0.209626; p2
143  err=0.00650077; p3 err=6.42207E-5; p4 err=1.94893E-7; Above 200cm,
144  Forpoly6 fit: p0 err=5.24743;p1 err=0.0176229; p2 err=1.6263E-5; p3
145  err=5.9155E-9; p4 err=9.71709E-13; p5 err=7.22381E-17;p6
146  err=1.9709E-21;*/
147 
149  //*********For muon, the calculations are valid up to 1.91E4 cm range
150  //corresponding to a Muon KE of 40 GeV**********//
152 
153  /*Proton range-momentum tables from CSDA (Argon density = 1.4 g/cm^3):
154  website: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html
155 
156  CSDA values:
157  double KE_MeV_P_Nist[31]={10, 15, 20, 30, 40, 80, 100, 150, 200, 250, 300,
158  350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
159  1500, 2000, 2500, 3000, 4000, 5000};
160 
161  double Range_gpercm_P_Nist[31]={1.887E-1,3.823E-1, 6.335E-1, 1.296,
162  2.159, 7.375, 1.092E1, 2.215E1, 3.627E1, 5.282E1, 7.144E1,
163  9.184E1, 1.138E2, 1.370E2, 1.614E2, 1.869E2, 2.132E2, 2.403E2,
164  2.681E2, 2.965E2, 3.254E2, 3.548E2, 3.846E2, 4.148E2, 4.454E2,
165  7.626E2, 1.090E3, 1.418E3, 1.745E3, 2.391E3, 3.022E3};
166 
167  Functions below are obtained by fitting power and polynomial fits to
168  KE_MeV vs Range (cm) graph. A better fit was obtained by splitting the
169  graph into two: Below range<=80cm,a a*(x^b) was a good fit; above 80cm, a
170  polynomial of power 6 was a good fit
171 
172  Fit errors for future purposes:
173  For power function fit: a=0.388873; and b=0.00347075
174  Forpoly6 fit: p0 err=3.49729;p1 err=0.0487859; p2 err=0.000225834; p3
175  err=4.45542E-7; p4 err=4.16428E-10; p5 err=1.81679E-13;p6
176  err=2.96958E-17;*/
177 
179  //*********For proton, the calculations are valid up to 3.022E3 cm range
180  //corresponding to a Muon KE of 5 GeV**********//
182 
183  if (trkrange < 0 || std::isnan(trkrange)) {
184  mf::LogError("TrackMomentumCalculator")
185  << "Invalid track range " << trkrange << " return -1";
186  return -1.;
187  }
188 
189  double KE, Momentum, M;
190  constexpr double Muon_M = 105.7, Proton_M = 938.272;
191 
192  if (abs(pdg) == 13) {
193  M = Muon_M;
194  KE = KEvsR_spline3.Eval(trkrange);
195  } else if (abs(pdg) == 2212) {
196  M = Proton_M;
197  if (trkrange > 0 && trkrange <= 80)
198  KE = 29.9317 * std::pow(trkrange, 0.586304);
199  else if (trkrange > 80 && trkrange <= 3.022E3)
200  KE =
201  149.904 + (3.34146 * trkrange) + (-0.00318856 * trkrange * trkrange) +
202  (4.34587E-6 * trkrange * trkrange * trkrange) +
203  (-3.18146E-9 * trkrange * trkrange * trkrange * trkrange) +
204  (1.17854E-12 * trkrange * trkrange * trkrange * trkrange * trkrange) +
205  (-1.71763E-16 * trkrange * trkrange * trkrange * trkrange * trkrange *
206  trkrange);
207  else
208  KE = -999;
209  } else
210  KE = -999;
211 
212  if (KE < 0)
213  Momentum = -999;
214  else
215  Momentum = std::sqrt((KE * KE) + (2 * M * KE));
216 
217  Momentum = Momentum / 1000;
218 
219  return Momentum;
220  }
221 
222  // Momentum measurement via Multiple Coulomb Scattering (MCS)
223 
224  // author: Leonidas N. Kalousis (July 2015)
225 
226  // email: kalousis@vt.edu
227 
228  double
229  TrackMomentumCalculator::GetMomentumMultiScatterLLHD(
230  const art::Ptr<recob::Track>& trk)
231  {
232  std::vector<float> recoX;
233  std::vector<float> recoY;
234  std::vector<float> recoZ;
235 
236  int n_points = trk->NumberTrajectoryPoints();
237 
238  for (int i = 0; i < n_points; i++) {
239  auto const& pos = trk->LocationAtPoint(i);
240  recoX.push_back(pos.X());
241  recoY.push_back(pos.Y());
242  recoZ.push_back(pos.Z());
243  }
244 
245  if (recoX.size() < 2)
246  return -1.0;
247 
248  if (!plotRecoTracks_(recoX, recoY, recoZ))
249  return -1.0;
250 
251  constexpr double seg_size{10.};
252 
253  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
254  if (!segments.has_value())
255  return -1.0;
256 
257  auto const seg_steps = segments->x.size();
258  if (seg_steps < 2)
259  return -1;
260 
261  double const recoL = segments->L.at(seg_steps - 1);
262  if (recoL < minLength || recoL > maxLength)
263  return -1;
264 
265  std::vector<float> dEi;
266  std::vector<float> dEj;
267  std::vector<float> dthij;
268  std::vector<float> ind;
269  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0)
270  return -1.0;
271 
272  double logL = 1e+16;
273  double bf = -666.0; // double errs = -666.0;
274 
275  int const start1{};
276  int const end1{750};
277  int const start2{};
278  int const end2{}; // 800.0;
279 
280  for (int k = start1; k <= end1; ++k) {
281  double const p_test = 0.001 + k * 0.01;
282 
283  for (int l = start2; l <= end2; ++l) {
284  double const res_test = 2.0; // 0.001+l*1.0;
285  double const fv = my_mcs_llhd(dEi, dEj, dthij, ind, p_test, res_test);
286 
287  if (fv < logL) {
288  bf = p_test;
289  logL = fv;
290  // errs = res_test;
291  }
292  }
293  }
294  return bf;
295  }
296 
297  TVector3
298  TrackMomentumCalculator::GetMultiScatterStartingPoint(
299  const art::Ptr<recob::Track>& trk)
300  {
301  double const LLHDp = GetMuMultiScatterLLHD3(trk, true);
302  double const LLHDm = GetMuMultiScatterLLHD3(trk, false);
303 
304  if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
305  int const n_points = trk->NumberTrajectoryPoints();
306  return trk->LocationAtPoint<TVector3>(n_points - 1);
307  } else {
308  return trk->LocationAtPoint<TVector3>(0);
309  }
310 
311  // Should never get here
312  return TVector3{};
313  }
314 
315  double
316  TrackMomentumCalculator::GetMuMultiScatterLLHD3(
317  art::Ptr<recob::Track> const& trk,
318  bool const dir)
319  {
320  std::vector<float> recoX;
321  std::vector<float> recoY;
322  std::vector<float> recoZ;
323 
324  int const n_points = trk->NumberTrajectoryPoints();
325  for (int i = 0; i < n_points; ++i) {
326  auto const index = dir ? i : n_points - 1 - i;
327  auto const& pos = trk->LocationAtPoint(index);
328  recoX.push_back(pos.X());
329  recoY.push_back(pos.Y());
330  recoZ.push_back(pos.Z());
331  }
332 
333  if (recoX.size() < 2)
334  return -1.0;
335 
336  if (!plotRecoTracks_(recoX, recoY, recoZ))
337  return -1.0;
338 
339  constexpr double seg_size{5.0};
340  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
341  if (!segments.has_value())
342  return -1.0;
343 
344  auto const seg_steps = segments->x.size();
345  if (seg_steps < 2)
346  return -1;
347 
348  double const recoL = segments->L.at(seg_steps - 1);
349  if (recoL < 15.0 || recoL > maxLength)
350  return -1;
351 
352  std::vector<float> dEi;
353  std::vector<float> dEj;
354  std::vector<float> dthij;
355  std::vector<float> ind;
356  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0)
357  return -1.0;
358 
359  double const p_range = recoL * kcal;
360  double const logL = my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
361 
362  return logL;
363  }
364 
365  int
366  TrackMomentumCalculator::getDeltaThetaij_(std::vector<float>& ei,
367  std::vector<float>& ej,
368  std::vector<float>& th,
369  std::vector<float>& ind,
370  Segments const& segments,
371  double const thick) const
372  {
373  int const a1 = segments.x.size();
374  int const a2 = segments.y.size();
375  int const a3 = segments.z.size();
376 
377  if (a1 != a2 || a1 != a3) {
378  std::cout << " ( Get thij ) Error ! " << std::endl;
379  return -1.0;
380  }
381 
382  auto const& segnx = segments.nx;
383  auto const& segny = segments.ny;
384  auto const& segnz = segments.nz;
385  auto const& segL = segments.L;
386 
387  int tot = a1 - 1;
388  double thick1 = thick + 0.13;
389 
390  for (int i = 0; i < tot; i++) {
391  double const dx = segnx.at(i);
392  double const dy = segny.at(i);
393  double const dz = segnz.at(i);
394 
395  TVector3 const vec_z{dx, dy, dz};
396  TVector3 vec_x;
397  TVector3 vec_y;
398 
399  double const switcher = basex.Dot(vec_z);
400  if (std::abs(switcher) <= 0.995) {
401  vec_y = vec_z.Cross(basex).Unit();
402  vec_x = vec_y.Cross(vec_z);
403  }
404  else {
405  // cout << " It switched ! Isn't this lovely !!! " << endl; //
406  vec_y = basez.Cross(vec_z).Unit();
407  vec_x = vec_y.Cross(vec_z);
408  }
409 
410  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
411  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
412  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
413 
414  double const refL = segL.at(i);
415 
416  for (int j = i; j < tot; j++) {
417  double const L1 = segL.at(j);
418  double const L2 = segL.at(j + 1);
419 
420  double const dz1 = L1 - refL;
421  double const dz2 = L2 - refL;
422 
423  if (dz1 <= thick1 && dz2 > thick1) {
424  double const here_dx = segnx.at(j);
425  double const here_dy = segny.at(j);
426  double const here_dz = segnz.at(j);
427 
428  TVector3 const here_vec{here_dx, here_dy, here_dz};
429  TVector3 const rot_here{
430  Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
431 
432  double const scx = rot_here.X();
433  double const scy = rot_here.Y();
434  double const scz = rot_here.Z();
435 
436  double const azy = find_angle(scz, scy);
437  double const azx = find_angle(scz, scx);
438 
439  constexpr double ULim = 10000.0;
440  constexpr double LLim = -10000.0;
441 
442  double const cL = kcal;
443  double const Li = segL.at(i);
444  double const Lj = segL.at(j);
445 
446  if (azy <= ULim && azy >= LLim) {
447  ei.push_back(Li * cL);
448  ej.push_back(Lj * cL);
449  th.push_back(azy);
450  ind.push_back(2);
451  }
452 
453  if (azx <= ULim && azx >= LLim) {
454  ei.push_back(Li * cL);
455  ej.push_back(Lj * cL);
456  th.push_back(azx);
457  ind.push_back(1);
458  }
459 
460  break; // of course !
461  }
462  }
463  }
464 
465  return 0;
466  }
467 
468  double
469  TrackMomentumCalculator::GetMomentumMultiScatterChi2(
470  const art::Ptr<recob::Track>& trk)
471  {
472  std::vector<float> recoX;
473  std::vector<float> recoY;
474  std::vector<float> recoZ;
475 
476  int n_points = trk->NumberTrajectoryPoints();
477 
478  for (int i = 0; i < n_points; i++) {
479  auto const& pos = trk->LocationAtPoint(i);
480  recoX.push_back(pos.X());
481  recoY.push_back(pos.Y());
482  recoZ.push_back(pos.Z());
483  }
484 
485  if (recoX.size() < 2)
486  return -1.0;
487 
488  if (!plotRecoTracks_(recoX, recoY, recoZ))
489  return -1.0;
490 
491  double const seg_size{steps_size};
492  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
493  if (!segments.has_value())
494  return -1.0;
495 
496  auto const seg_steps = segments->x.size();
497  if (seg_steps < 2)
498  return -1;
499 
500  double const recoL = segments->L.at(seg_steps - 1);
501  if (recoL < minLength || recoL > maxLength)
502  return -1;
503 
504  double ymax = -999.0;
505  double ymin = +999.0;
506 
507  std::vector<double> xmeas;
508  std::vector<double> ymeas;
509  std::vector<double> eymeas;
510  xmeas.reserve(n_steps);
511  ymeas.reserve(n_steps);
512  eymeas.reserve(n_steps);
513  for (int j = 0; j < n_steps; j++) {
514  double const trial = steps.at(j);
515  auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
516 
517  if (std::isnan(mean) || std::isinf(mean)) {
518  mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
519  continue;
520  }
521  if (std::isnan(rms) || std::isinf(rms)) {
522  mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
523  continue;
524  }
525  if (std::isnan(rmse) || std::isinf(rmse)) {
526  mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
527  continue;
528  }
529 
530  xmeas.push_back(trial); // Is this what is intended?
531  ymeas.push_back(rms);
532  eymeas.push_back(std::sqrt(cet::sum_of_squares(rmse, 0.05 * rms))); // <--- conservative syst. error to fix chi^{2} behaviour !!!
533 
534  if (ymin > rms)
535  ymin = rms;
536  if (ymax < rms)
537  ymax = rms;
538  }
539 
540  assert(xmeas.size() == ymeas.size());
541  assert(xmeas.size() == eymeas.size());
542  if (xmeas.empty()) {
543  return -1.0;
544  }
545 
546  TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
547 
548  gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
549  "thickness in cm; (#Delta#theta)_{rms} in mrad");
550 
551  gr_meas.SetLineColor(kBlack);
552  gr_meas.SetMarkerColor(kBlack);
553  gr_meas.SetMarkerStyle(20);
554  gr_meas.SetMarkerSize(1.2);
555 
556  gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0),
557  steps.at(n_steps - 1) + steps.at(0));
558  gr_meas.SetMinimum(0.0);
559  gr_meas.SetMaximum(1.80 * ymax);
560 
561  ROOT::Minuit2::Minuit2Minimizer mP{};
562  FcnWrapper const wrapper{move(xmeas), move(ymeas), move(eymeas)};
563  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
564 
565  mP.SetFunction(FCA);
566  mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, 7.5);
567  mP.SetLimitedVariable(1, "#delta#theta", 0.0, 1.0, 0.0, 45.0);
568  mP.SetMaxFunctionCalls(1.E9);
569  mP.SetMaxIterations(1.E9);
570  mP.SetTolerance(0.01);
571  mP.SetStrategy(2);
572  mP.SetErrorDef(1.0);
573 
574  bool const mstatus = mP.Minimize();
575 
576  mP.Hesse();
577 
578  const double* pars = mP.X();
579  const double* erpars = mP.Errors();
580 
581  double const deltap = (recoL * kcal) / 2.0;
582 
583  double const p_mcs = pars[0] + deltap;
584  double const p_mcs_e [[maybe_unused]] = erpars[0];
585  return mstatus ? p_mcs : -1.0;
586  }
587 
588  bool
589  TrackMomentumCalculator::plotRecoTracks_(std::vector<float> const& xxx,
590  std::vector<float> const& yyy,
591  std::vector<float> const& zzz)
592  {
593  auto const n = xxx.size();
594  auto const y_size = yyy.size();
595  auto const z_size = zzz.size();
596 
597  if (n != y_size || n != z_size) {
598  cout << " ( Get reco tacks ) Error ! " << endl;
599  return false;
600  }
601 
602  // Here, we perform a const-cast to float* because, sadly,
603  // TPolyLine3D requires a pointer to a non-const object. We will
604  // trust that ROOT does not mess around with the underlying data.
605  auto xs = const_cast<float*>(xxx.data());
606  auto ys = const_cast<float*>(yyy.data());
607  auto zs = const_cast<float*>(zzz.data());
608 
609  auto const narrowed_size =
610  static_cast<int>(n); // ROOT requires a signed integral type
611  delete gr_reco_xyz;
612  gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
613  gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
614  gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
615  gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
616 
617  return true;
618  }
619 
620  std::optional<TrackMomentumCalculator::Segments>
621  TrackMomentumCalculator::getSegTracks_(std::vector<float> const& xxx,
622  std::vector<float> const& yyy,
623  std::vector<float> const& zzz,
624  double const seg_size)
625  {
626  double stag = 0.0;
627 
628  int a1 = xxx.size();
629  int a2 = yyy.size();
630  int a3 = zzz.size();
631 
632  if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
633  cout << " ( Digitize reco tacks ) Error ! " << endl;
634  return std::nullopt;
635  }
636 
637  int const stopper = seg_stop / seg_size;
638 
639  std::vector<float> segx, segnx;
640  std::vector<float> segy, segny;
641  std::vector<float> segz, segnz;
642  std::vector<float> segL;
643 
644  int ntot = 0;
645 
646  n_seg = 0;
647 
648  double x0{};
649  double y0{};
650  double z0{};
651 
652  double x00 = xxx.at(0);
653  double y00 = yyy.at(0);
654  double z00 = zzz.at(0);
655 
656  int indC = 0;
657 
658  std::vector<float> vx;
659  std::vector<float> vy;
660  std::vector<float> vz;
661 
662  for (int i = 0; i < a1; i++) {
663  x0 = xxx.at(i);
664  y0 = yyy.at(i);
665  z0 = zzz.at(i);
666 
667  double const RR0 =
668  std::sqrt(cet::sum_of_squares(x00 - x0, y00 - y0, z00 - z0));
669 
670  if (RR0 >= stag) {
671  segx.push_back(x0);
672  segy.push_back(y0);
673  segz.push_back(z0);
674 
675  segL.push_back(stag);
676 
677  x_seg[n_seg] = x0;
678  y_seg[n_seg] = y0;
679  z_seg[n_seg] = z0;
680 
681  n_seg++;
682 
683  vx.push_back(x0);
684  vy.push_back(y0);
685  vz.push_back(z0);
686 
687  ntot++;
688 
689  indC = i + 1;
690 
691  break;
692  }
693  }
694 
695  for (int i = indC; i < a1 - 1; i++) {
696  double const x1 = xxx.at(i);
697  double const y1 = yyy.at(i);
698  double const z1 = zzz.at(i);
699 
700  double const dr1 = std::sqrt(cet::sum_of_squares(x1 - x0, y1 - y0, z1 - z0));
701 
702  double const x2 = xxx.at(i + 1);
703  double const y2 = yyy.at(i + 1);
704  double const z2 = zzz.at(i + 1);
705 
706  double const dr2 = std::sqrt(cet::sum_of_squares(x2 - x0, y2 - y0, z2 - z0));
707 
708  if (dr1 < seg_size) {
709  vx.push_back(x1);
710  vy.push_back(y1);
711  vz.push_back(z1);
712 
713  ntot++;
714  }
715 
716  if (dr1 <= seg_size && dr2 > seg_size) {
717  double const dx = x2 - x1;
718  double const dy = y2 - y1;
719  double const dz = z2 - z1;
720  double const dr = std::sqrt(dx * dx + dy * dy + dz * dz);
721 
722  if (dr == 0) {
723  cout << " ( Zero ) Error ! " << endl;
724  return std::nullopt;
725  }
726 
727  double const beta =
728  2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
729 
730  double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
731  double const delta = beta * beta - 4.0 * gamma;
732 
733  if (delta < 0.0) {
734  cout << " ( Discriminant ) Error ! " << endl;
735  return std::nullopt;
736  }
737 
738  double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
739  double const t = lysi1;
740 
741  double const xp = x1 + t * dx;
742  double const yp = y1 + t * dy;
743  double const zp = z1 + t * dz;
744 
745  segx.push_back(xp);
746  segy.push_back(yp);
747  segz.push_back(zp);
748 
749  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
750 
751  x_seg[n_seg] = xp;
752  y_seg[n_seg] = yp;
753  z_seg[n_seg] = zp;
754  n_seg++;
755 
756  x0 = xp;
757  y0 = yp;
758  z0 = zp;
759 
760  vx.push_back(x0);
761  vy.push_back(y0);
762  vz.push_back(z0);
763 
764  ntot++;
765 
766  auto const na = vx.size();
767  double sumx = 0.0;
768  double sumy = 0.0;
769  double sumz = 0.0;
770 
771  for (std::size_t i = 0; i < na; ++i) {
772  sumx += vx.at(i);
773  sumy += vy.at(i);
774  sumz += vz.at(i);
775  }
776 
777  sumx /= na;
778  sumy /= na;
779  sumz /= na;
780 
781  std::vector<double> mx;
782  std::vector<double> my;
783  std::vector<double> mz;
784 
785  TMatrixDSym m(3);
786 
787  for (std::size_t i = 0; i < na; ++i) {
788  double const xxw1 = vx.at(i);
789  double const yyw1 = vy.at(i);
790  double const zzw1 = vz.at(i);
791 
792  mx.push_back(xxw1 - sumx);
793  my.push_back(yyw1 - sumy);
794  mz.push_back(zzw1 - sumz);
795 
796  double const xxw0 = mx.at(i);
797  double const yyw0 = my.at(i);
798  double const zzw0 = mz.at(i);
799 
800  m(0, 0) += xxw0 * xxw0 / na;
801  m(0, 1) += xxw0 * yyw0 / na;
802  m(0, 2) += xxw0 * zzw0 / na;
803 
804  m(1, 0) += yyw0 * xxw0 / na;
805  m(1, 1) += yyw0 * yyw0 / na;
806  m(1, 2) += yyw0 * zzw0 / na;
807 
808  m(2, 0) += zzw0 * xxw0 / na;
809  m(2, 1) += zzw0 * yyw0 / na;
810  m(2, 2) += zzw0 * zzw0 / na;
811  }
812 
813  TMatrixDSymEigen me(m);
814 
815  TVectorD eigenval = me.GetEigenValues();
816  TMatrixD eigenvec = me.GetEigenVectors();
817 
818  double max1 = -666.0;
819 
820  double ind1 = 0;
821 
822  for (int i = 0; i < 3; ++i) {
823  double const p1 = eigenval(i);
824 
825  if (p1 > max1) {
826  max1 = p1;
827  ind1 = i;
828  }
829  }
830 
831  double ax = eigenvec(0, ind1);
832  double ay = eigenvec(1, ind1);
833  double az = eigenvec(2, ind1);
834 
835  if (n_seg > 1) {
836  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
837  ax = std::abs(ax);
838  else
839  ax = -1.0 * std::abs(ax);
840 
841  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
842  ay = std::abs(ay);
843  else
844  ay = -1.0 * std::abs(ay);
845 
846  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
847  az = std::abs(az);
848  else
849  az = -1.0 * std::abs(az);
850 
851  segnx.push_back(ax);
852  segny.push_back(ay);
853  segnz.push_back(az);
854  }
855  else
856  return std::nullopt;
857 
858  ntot = 0;
859 
860  vx.clear();
861  vy.clear();
862  vz.clear();
863 
864  vx.push_back(x0);
865  vy.push_back(y0);
866  vz.push_back(z0);
867 
868  ntot++;
869  }
870  else if (dr1 > seg_size) {
871  double const dx = x1 - x0;
872  double const dy = y1 - y0;
873  double const dz = z1 - z0;
874  double const dr = std::sqrt(cet::sum_of_squares(dx, dy, dz));
875 
876  if (dr == 0) {
877  cout << " ( Zero ) Error ! " << endl;
878  return std::nullopt;
879  }
880 
881  double const t = seg_size / dr;
882  double const xp = x0 + t * dx;
883  double const yp = y0 + t * dy;
884  double const zp = z0 + t * dz;
885 
886  segx.push_back(xp);
887  segy.push_back(yp);
888  segz.push_back(zp);
889  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
890 
891  x_seg[n_seg] = xp;
892  y_seg[n_seg] = yp;
893  z_seg[n_seg] = zp;
894  n_seg++;
895 
896  x0 = xp;
897  y0 = yp;
898  z0 = zp;
899 
900  i = (i - 1);
901 
902  vx.push_back(x0);
903  vy.push_back(y0);
904  vz.push_back(z0);
905 
906  ntot++;
907 
908  double na = vx.size();
909 
910  double sumx = 0.0;
911  double sumy = 0.0;
912  double sumz = 0.0;
913 
914  for (std::size_t i = 0; i < na; ++i) {
915  sumx += vx.at(i);
916  sumy += vy.at(i);
917  sumz += vz.at(i);
918  }
919 
920  sumx /= na;
921  sumy /= na;
922  sumz /= na;
923 
924  std::vector<double> mx;
925  std::vector<double> my;
926  std::vector<double> mz;
927 
928  TMatrixDSym m(3);
929 
930  for (int i = 0; i < na; ++i) {
931  double const xxw1 = vx.at(i);
932  double const yyw1 = vy.at(i);
933  double const zzw1 = vz.at(i);
934 
935  mx.push_back(xxw1 - sumx);
936  my.push_back(yyw1 - sumy);
937  mz.push_back(zzw1 - sumz);
938 
939  double const xxw0 = mx.at(i);
940  double const yyw0 = my.at(i);
941  double const zzw0 = mz.at(i);
942 
943  m(0, 0) += xxw0 * xxw0 / na;
944  m(0, 1) += xxw0 * yyw0 / na;
945  m(0, 2) += xxw0 * zzw0 / na;
946 
947  m(1, 0) += yyw0 * xxw0 / na;
948  m(1, 1) += yyw0 * yyw0 / na;
949  m(1, 2) += yyw0 * zzw0 / na;
950 
951  m(2, 0) += zzw0 * xxw0 / na;
952  m(2, 1) += zzw0 * yyw0 / na;
953  m(2, 2) += zzw0 * zzw0 / na;
954  }
955 
956  TMatrixDSymEigen me(m);
957 
958  TVectorD eigenval = me.GetEigenValues();
959  TMatrixD eigenvec = me.GetEigenVectors();
960 
961  double max1 = -666.0;
962  double ind1 = 0;
963 
964  for (int i = 0; i < 3; ++i) {
965  double const p1 = eigenval(i);
966 
967  if (p1 > max1) {
968  max1 = p1;
969  ind1 = i;
970  }
971  }
972 
973  double ax = eigenvec(0, ind1);
974  double ay = eigenvec(1, ind1);
975  double az = eigenvec(2, ind1);
976 
977  if (n_seg > 1) {
978  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
979  ax = std::abs(ax);
980  else
981  ax = -1.0 * std::abs(ax);
982 
983  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
984  ay = std::abs(ay);
985  else
986  ay = -1.0 * std::abs(ay);
987 
988  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
989  az = std::abs(az);
990  else
991  az = -1.0 * std::abs(az);
992 
993  segnx.push_back(ax);
994  segny.push_back(ay);
995  segnz.push_back(az);
996  }
997 
998  else
999  return std::nullopt;
1000 
1001  ntot = 0;
1002 
1003  vx.clear();
1004  vy.clear();
1005  vz.clear();
1006 
1007  vx.push_back(x0);
1008  vy.push_back(y0);
1009  vz.push_back(z0);
1010 
1011  ntot++;
1012  }
1013 
1014  if (n_seg >= (stopper + 1.0) && seg_stop != -1)
1015  break;
1016  }
1017 
1018  delete gr_seg_xyz;
1019  gr_seg_xyz = new TPolyLine3D{n_seg, z_seg, x_seg, y_seg};
1020  gr_seg_yz = TGraph{n_seg, z_seg, y_seg};
1021  gr_seg_xz = TGraph{n_seg, z_seg, x_seg};
1022  gr_seg_xy = TGraph{n_seg, x_seg, y_seg};
1023 
1024  return std::make_optional<Segments>(Segments{segx, segnx, segy, segny, segz, segnz, segL});
1025  }
1026 
1027  std::tuple<double, double, double>
1028  TrackMomentumCalculator::getDeltaThetaRMS_(Segments const& segments,
1029  double const thick) const
1030  {
1031  auto const& segnx = segments.nx;
1032  auto const& segny = segments.ny;
1033  auto const& segnz = segments.nz;
1034  auto const& segL = segments.L;
1035 
1036  int const a1 = segnx.size();
1037  int const a2 = segny.size();
1038  int const a3 = segnz.size();
1039 
1040  if (a1 != a2 || a1 != a3) {
1041  cout << " ( Get RMS ) Error ! " << endl;
1042  return std::make_tuple(0., 0., 0.);
1043  }
1044 
1045  int const tot = a1 - 1;
1046 
1047  double const thick1 = thick + 0.13;
1048 
1049  std::vector<float> buf0;
1050 
1051  for (int i = 0; i < tot; i++) {
1052  double const dx = segnx.at(i);
1053  double const dy = segny.at(i);
1054  double const dz = segnz.at(i);
1055 
1056  TVector3 const vec_z{dx, dy, dz};
1057  TVector3 vec_x;
1058  TVector3 vec_y;
1059 
1060  double const switcher = basex.Dot(vec_z);
1061 
1062  if (switcher <= 0.995) {
1063  vec_y = vec_z.Cross(basex).Unit();
1064  vec_x = vec_y.Cross(vec_z);
1065  }
1066  else {
1067  // cout << " It switched ! Isn't this lovely !!! " << endl;
1068  vec_y = basez.Cross(vec_z).Unit();
1069  vec_x = vec_y.Cross(vec_z);
1070  }
1071 
1072  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
1073  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
1074  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
1075 
1076  double const refL = segL.at(i);
1077 
1078  for (int j = i; j < tot; j++) {
1079  double const L1 = segL.at(j);
1080  double const L2 = segL.at(j + 1);
1081 
1082  double const dz1 = L1 - refL;
1083  double const dz2 = L2 - refL;
1084 
1085  if (dz1 <= thick1 && dz2 > thick1) {
1086  double const here_dx = segnx.at(j);
1087  double const here_dy = segny.at(j);
1088  double const here_dz = segnz.at(j);
1089 
1090  TVector3 const here_vec{here_dx, here_dy, here_dz};
1091  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
1092 
1093  double const scx = rot_here.X();
1094  double const scy = rot_here.Y();
1095  double const scz = rot_here.Z();
1096 
1097  double azy = find_angle(scz, scy);
1098  azy *= 1.0;
1099 
1100  double const azx = find_angle(scz, scx);
1101 
1102  double const ULim = 10000.0;
1103  double const LLim = -10000.0;
1104 
1105  if (azx <= ULim && azx >= LLim) {
1106  buf0.push_back(azx);
1107  }
1108 
1109  break; // of course !
1110  }
1111  }
1112  }
1113 
1114  int const nmeas = buf0.size();
1115  double nnn = 0.0;
1116 
1117  double mean = 0.0;
1118  double rms = 0.0;
1119  double rmse = 0.0;
1120 
1121  for (int i = 0; i < nmeas; i++) {
1122  mean += buf0.at(i);
1123  nnn++;
1124  }
1125  mean = mean / nnn;
1126 
1127  for (int i = 0; i < nmeas; i++)
1128  rms += ((buf0.at(i)) * (buf0.at(i)));
1129 
1130  rms = rms / (nnn);
1131  rms = std::sqrt(rms);
1132  rmse = rms / std::sqrt(2.0 * tot);
1133 
1134  double rms1 = rms;
1135 
1136  rms = 0.0;
1137 
1138  double ntot1 = 0.0;
1139  double const lev1 = 2.50;
1140 
1141  for (int i = 0; i < nmeas; i++) {
1142  double const amp = buf0.at(i);
1143  if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
1144  ++ntot1;
1145  rms += amp * amp;
1146  }
1147  }
1148 
1149  rms = rms / (ntot1);
1150  rms = std::sqrt(rms);
1151  rmse = rms / std::sqrt(2.0 * ntot1);
1152  return std::make_tuple(mean, rms, rmse);
1153  }
1154 
1155  double
1156  TrackMomentumCalculator::find_angle(double vz, double vy) const
1157  {
1158  double thetayz = -999.0;
1159 
1160  if (vz > 0 && vy > 0) {
1161  double ratio = std::abs(vy / vz);
1162  thetayz = std::atan(ratio);
1163  }
1164 
1165  else if (vz < 0 && vy > 0) {
1166  double ratio = std::abs(vy / vz);
1167  thetayz = std::atan(ratio);
1168  thetayz = TMath::Pi() - thetayz;
1169  }
1170 
1171  else if (vz < 0 && vy < 0) {
1172  double ratio = std::abs(vy / vz);
1173  thetayz = std::atan(ratio);
1174  thetayz = thetayz + TMath::Pi();
1175  }
1176 
1177  else if (vz > 0 && vy < 0) {
1178  double ratio = std::abs(vy / vz);
1179  thetayz = std::atan(ratio);
1180  thetayz = 2.0 * TMath::Pi() - thetayz;
1181  }
1182 
1183  else if (vz == 0 && vy > 0) {
1184  thetayz = TMath::Pi() / 2.0;
1185  }
1186 
1187  else if (vz == 0 && vy < 0) {
1188  thetayz = 3.0 * TMath::Pi() / 2.0;
1189  }
1190 
1191  if (thetayz > TMath::Pi()) {
1192  thetayz = thetayz - 2.0 * TMath::Pi();
1193  }
1194 
1195  return 1000.0 * thetayz;
1196  }
1197 
1198  double
1199  TrackMomentumCalculator::my_g(double xx, double Q, double s) const
1200  {
1201  if (s == 0.) {
1202  cout << " Error : The code tries to divide by zero ! " << endl;
1203  return 0.;
1204  }
1205 
1206  double const arg = (xx - Q) / s;
1207  double const result =
1208  -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1209 
1210  if (std::isnan(result) || std::isinf(result)) {
1211  cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
1212  }
1213 
1214  return result;
1215  }
1216 
1217  double
1218  TrackMomentumCalculator::my_mcs_llhd(std::vector<float> const& dEi,
1219  std::vector<float> const& dEj,
1220  std::vector<float> const& dthij,
1221  std::vector<float> const& ind,
1222  double const x0,
1223  double const x1) const
1224  {
1225  double p = x0;
1226  double theta0x = x1;
1227 
1228  double result = 0.0;
1229  double nnn1 = dEi.size();
1230 
1231  double red_length = (10.0) / 14.0;
1232  double addth = 0;
1233 
1234  for (int i = 0; i < nnn1; i++) {
1235  double Ei = p - dEi.at(i);
1236  double Ej = p - dEj.at(i);
1237 
1238  if (Ei > 0 && Ej < 0)
1239  addth = 3.14 * 1000.0;
1240 
1241  Ei = std::abs(Ei);
1242  Ej = std::abs(Ej);
1243 
1244  double tH0 = (13.6 / std::sqrt(Ei * Ej)) *
1245  (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1246 
1247  double rms = -1.0;
1248 
1249  if (ind.at(i) == 1) {
1250  rms = std::sqrt(tH0 * tH0 + pow(theta0x, 2.0));
1251 
1252  double const DT = dthij.at(i) + addth;
1253  double const prob = my_g(DT, 0.0, rms);
1254 
1255  result = result - 2.0 * prob;
1256  }
1257  }
1258 
1259  if (std::isnan(result) || std::isinf(result)) {
1260  std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
1261  }
1262  return result;
1263  }
1264 
1265 } // namespace track
x
Float_t x
Definition: compare.C:6
s
Float_t s
Definition: plot.C:23
xx
Double_t xx
Definition: macro.C:12
E
Float_t E
Definition: plot.C:23
trkf::TrackMomentumCalculator::my_g
double my_g(double xx, double Q, double s) const
Definition: TrackMomentumCalculator.cxx:1199
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2
double GetMomentumMultiScatterChi2(art::Ptr< recob::Track > const &trk)
Definition: TrackMomentumCalculator.cxx:469
trkf::TrackMomentumCalculator::getDeltaThetaij_
int getDeltaThetaij_(std::vector< float > &ei, std::vector< float > &ej, std::vector< float > &th, std::vector< float > &ind, Segments const &segments, double thick) const
Definition: TrackMomentumCalculator.cxx:366
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:129
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:105
beta
Double_t beta
Definition: TG43_relative_dose.C:15
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:214
trkf::TrackMomentumCalculator::getSegTracks_
std::optional< Segments > getSegTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
Definition: TrackMomentumCalculator.cxx:621
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
a2
#define a2
Definition: HoughBaseAlg.cxx:67
max
Int_t max
Definition: plot.C:27
a3
#define a3
Definition: HoughBaseAlg.cxx:68
trkf::TrackMomentumCalculator::find_angle
double find_angle(double vz, double vy) const
Definition: TrackMomentumCalculator.cxx:1156
trkf::TrackMomentumCalculator::getDeltaThetaRMS_
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Definition: TrackMomentumCalculator.cxx:1028
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:15
fhicl::detail::atom::value
std::string value(boost::any const &)
Definition: printing_helpers.cc:97
trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Definition: TrackMomentumCalculator.cxx:316
dir
TDirectory * dir
Definition: macro.C:5
min
Int_t min
Definition: plot.C:26
trkf::TrackMomentumCalculator::my_mcs_llhd
double my_mcs_llhd(std::vector< float > const &dEi, std::vector< float > const &dEj, std::vector< float > const &dthij, std::vector< float > const &ind, double x0, double x1) const
Definition: TrackMomentumCalculator.cxx:1218
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:211
trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint
TVector3 GetMultiScatterStartingPoint(art::Ptr< recob::Track > const &trk)
Definition: TrackMomentumCalculator.cxx:298
LOG_DEBUG
#define LOG_DEBUG(id)
Definition: MessageLogger.h:263
n
Char_t n[5]
Definition: comparison_ascii.C:55
x2
Float_t x2[n_points_geant4]
Definition: compare.C:26
e
Float_t e
Definition: plot.C:34
trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD
double GetMomentumMultiScatterLLHD(art::Ptr< recob::Track > const &trk)
Definition: TrackMomentumCalculator.cxx:229
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg) const
Definition: TrackMomentumCalculator.cxx:120
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:589
a1
#define a1
Definition: HoughBaseAlg.cxx:66