LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::TrackMomentumCalculator Class Reference

#include "TrackMomentumCalculator.h"

Classes

struct  Segments
 Struct to store segments. x, y and z are the 3D points of the segment nx, ny, nz forms the vector that point to the direction of most scattering L is the length of the segment, using steps_size. More...
 

Public Member Functions

 TrackMomentumCalculator (double minLength=100.0, double maxLength=1350.0, double steps_size=10.)
 Constructor. More...
 
double GetTrackMomentum (double trkrange, int pdg) const
 
double GetMomentumMultiScatterChi2 (art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500)
 Calculate muon momentum (GeV) using multiple coulomb scattering. Chi2 minimization of the Highland formula. More...
 
double GetMomentumMultiScatterLLHD (art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const int MomentumStep_MeV=10, const int max_resolution=0)
 Calculate muon momentum (GeV) using multiple coulomb scattering by log likelihood. More...
 
double GetMuMultiScatterLLHD3 (art::Ptr< recob::Track > const &trk, bool dir)
 
TVector3 GetMultiScatterStartingPoint (art::Ptr< recob::Track > const &trk)
 

Private Member Functions

bool plotRecoTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
 
void compute_max_fluctuation_vector (const std::vector< float > segx, const std::vector< float > segy, const std::vector< float > segz, std::vector< float > &segnx, std::vector< float > &segny, std::vector< float > &segnz, std::vector< float > &vx, std::vector< float > &vy, std::vector< float > &vz)
 Computes the vector with most scattering inside a segment with size steps_size. More...
 
std::optional< Segments > getSegTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
 Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/P10010. More...
 
std::tuple< double, double, double > getDeltaThetaRMS_ (Segments const &segments, double thick) const
 Gets the scattered angle RMS for a all segments. More...
 
int getDeltaThetaij_ (std::vector< float > &ei, std::vector< float > &ej, std::vector< float > &th, std::vector< float > &ind, Segments const &segments, double thick) const
 Gets the scatterd angle for all the segments. More...
 
double my_g (double xx, double Q, double s) const
 chi square minizer using Minuit2, it will minize (xx-Q)/s More...
 
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
 Minimizer of log likelihood for scattered angle. More...
 
double find_angle (double vz, double vy) const
 Gets angle between two vy and vz. More...
 

Private Attributes

float seg_stop {-1.}
 
int n_seg {}
 
float x_seg [50000]
 
float y_seg [50000]
 
float z_seg [50000]
 
int n_steps {6}
 
std::vector< float > steps
 
double minLength
 
double maxLength
 
double steps_size
 
double rad_length {14.0}
 
TPolyLine3D * gr_reco_xyz {nullptr}
 
TGraph gr_reco_xy {}
 
TGraph gr_reco_yz {}
 
TGraph gr_reco_xz {}
 
TPolyLine3D * gr_seg_xyz {nullptr}
 
TGraph gr_seg_xy {}
 
TGraph gr_seg_yz {}
 
TGraph gr_seg_xz {}
 

Detailed Description

Definition at line 21 of file TrackMomentumCalculator.h.

Constructor & Destructor Documentation

trkf::TrackMomentumCalculator::TrackMomentumCalculator ( double  minLength = 100.0,
double  maxLength = 1350.0,
double  steps_size = 10. 
)

Constructor.

Parameters are relevant for Multiple Coulomb Scattering

Parameters
minLengthminimum length in cm of tracks (length here is based on the amout of segments)
maxLengthmaximum length in cm of tracks (length here is based on the amout of segments)
steps_sizesize in cm of each segment to compute scattering

Definition at line 129 of file TrackMomentumCalculator.cxx.

References maxLength, n_steps, steps, and steps_size.

132  : minLength{min}, maxLength{max}, steps_size{stepsize}
133  {
134  for (int i = 1; i <= n_steps; i++) {
135  steps.push_back(steps_size * i);
136  }
137  }

Member Function Documentation

void trkf::TrackMomentumCalculator::compute_max_fluctuation_vector ( const std::vector< float >  segx,
const std::vector< float >  segy,
const std::vector< float >  segz,
std::vector< float > &  segnx,
std::vector< float > &  segny,
std::vector< float > &  segnz,
std::vector< float > &  vx,
std::vector< float > &  vy,
std::vector< float > &  vz 
)
private

Computes the vector with most scattering inside a segment with size steps_size.

Parameters
segx,segy,segzsegments points
segnx,segny,segnzvector components to be filled
vectorused to control points to be used at segments

Definition at line 616 of file TrackMomentumCalculator.cxx.

References util::abs(), and n_seg.

Referenced by getSegTracks_().

625  {
626  auto const na = vx.size();
627 
628  double sumx = 0.0;
629  double sumy = 0.0;
630  double sumz = 0.0;
631 
632  for (std::size_t i = 0; i < na; ++i) {
633  sumx += vx.at(i);
634  sumy += vy.at(i);
635  sumz += vz.at(i);
636  }
637 
638  sumx /= na;
639  sumy /= na;
640  sumz /= na;
641 
642  std::vector<double> mx;
643  std::vector<double> my;
644  std::vector<double> mz;
645 
646  TMatrixDSym m(3);
647 
648  for (std::size_t i = 0; i < na; ++i) {
649  double const xxw1 = vx.at(i);
650  double const yyw1 = vy.at(i);
651  double const zzw1 = vz.at(i);
652 
653  mx.push_back(xxw1 - sumx);
654  my.push_back(yyw1 - sumy);
655  mz.push_back(zzw1 - sumz);
656 
657  double const xxw0 = mx.at(i);
658  double const yyw0 = my.at(i);
659  double const zzw0 = mz.at(i);
660 
661  m(0, 0) += xxw0 * xxw0 / na;
662  m(0, 1) += xxw0 * yyw0 / na;
663  m(0, 2) += xxw0 * zzw0 / na;
664 
665  m(1, 0) += yyw0 * xxw0 / na;
666  m(1, 1) += yyw0 * yyw0 / na;
667  m(1, 2) += yyw0 * zzw0 / na;
668 
669  m(2, 0) += zzw0 * xxw0 / na;
670  m(2, 1) += zzw0 * yyw0 / na;
671  m(2, 2) += zzw0 * zzw0 / na;
672  }
673 
674  TMatrixDSymEigen me(m);
675 
676  TVectorD eigenval = me.GetEigenValues();
677  TMatrixD eigenvec = me.GetEigenVectors();
678 
679  double max1 = -666.0;
680 
681  double ind1 = 0;
682 
683  for (int i = 0; i < 3; ++i) {
684  double const p1 = eigenval(i);
685 
686  if (p1 > max1) {
687  max1 = p1;
688  ind1 = i;
689  }
690  }
691 
692  double ax = eigenvec(0, ind1);
693  double ay = eigenvec(1, ind1);
694  double az = eigenvec(2, ind1);
695 
696  if (n_seg > 1) {
697  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
698  ax = std::abs(ax);
699  else
700  ax = -1.0 * std::abs(ax);
701 
702  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
703  ay = std::abs(ay);
704  else
705  ay = -1.0 * std::abs(ay);
706 
707  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
708  az = std::abs(az);
709  else
710  az = -1.0 * std::abs(az);
711 
712  segnx.push_back(ax);
713  segny.push_back(ay);
714  segnz.push_back(az);
715  }
716  vx.clear();
717  vy.clear();
718  vz.clear();
719  }
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
double trkf::TrackMomentumCalculator::find_angle ( double  vz,
double  vy 
) const
private

Gets angle between two vy and vz.

Parameters
vz
vy
Returns
angle in mrad

Definition at line 1057 of file TrackMomentumCalculator.cxx.

References util::abs().

Referenced by getDeltaThetaij_(), and getDeltaThetaRMS_().

1058  {
1059  double thetayz = -999.0;
1060 
1061  if (vz > 0 && vy > 0) {
1062  double ratio = std::abs(vy / vz);
1063  thetayz = std::atan(ratio);
1064  }
1065 
1066  else if (vz < 0 && vy > 0) {
1067  double ratio = std::abs(vy / vz);
1068  thetayz = std::atan(ratio);
1069  thetayz = TMath::Pi() - thetayz;
1070  }
1071 
1072  else if (vz < 0 && vy < 0) {
1073  double ratio = std::abs(vy / vz);
1074  thetayz = std::atan(ratio);
1075  thetayz = thetayz + TMath::Pi();
1076  }
1077 
1078  else if (vz > 0 && vy < 0) {
1079  double ratio = std::abs(vy / vz);
1080  thetayz = std::atan(ratio);
1081  thetayz = 2.0 * TMath::Pi() - thetayz;
1082  }
1083 
1084  else if (vz == 0 && vy > 0) {
1085  thetayz = TMath::Pi() / 2.0;
1086  }
1087 
1088  else if (vz == 0 && vy < 0) {
1089  thetayz = 3.0 * TMath::Pi() / 2.0;
1090  }
1091 
1092  if (thetayz > TMath::Pi()) { thetayz = thetayz - 2.0 * TMath::Pi(); }
1093 
1094  return 1000.0 * thetayz;
1095  }
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
int trkf::TrackMomentumCalculator::getDeltaThetaij_ ( std::vector< float > &  ei,
std::vector< float > &  ej,
std::vector< float > &  th,
std::vector< float > &  ind,
Segments const &  segments,
double  thick 
) const
private

Gets the scatterd angle for all the segments.

Parameters
ei,ejwill be filled with energy lost at point i and j
thdelta theta to be filled
indselection of scattering plane
segments
thickis the steps_size
Returns
sucess or failure

Definition at line 371 of file TrackMomentumCalculator.cxx.

References a1, a2, a3, util::abs(), find_angle(), trkf::TrackMomentumCalculator::Segments::L, trkf::TrackMomentumCalculator::Segments::nx, trkf::TrackMomentumCalculator::Segments::ny, trkf::TrackMomentumCalculator::Segments::nz, trkf::TrackMomentumCalculator::Segments::x, trkf::TrackMomentumCalculator::Segments::y, and trkf::TrackMomentumCalculator::Segments::z.

Referenced by GetMomentumMultiScatterLLHD(), and GetMuMultiScatterLLHD3().

377  {
378  int const a1 = segments.x.size();
379  int const a2 = segments.y.size();
380  int const a3 = segments.z.size();
381 
382  if (a1 != a2 || a1 != a3) {
383  std::cout << " ( Get thij ) Error ! " << std::endl;
384  return -1.0;
385  }
386 
387  auto const& segnx = segments.nx;
388  auto const& segny = segments.ny;
389  auto const& segnz = segments.nz;
390  auto const& segL = segments.L;
391 
392  int tot = a1 - 1;
393  double thick1 = thick + 0.13; // adds a small offset to the 10 cm segment
394 
395  for (int i = 0; i < tot; i++) {
396  double const dx = segnx.at(i);
397  double const dy = segny.at(i);
398  double const dz = segnz.at(i);
399 
400  // Assumes z as propagation angle
401  TVector3 const vec_z{dx, dy, dz};
402  TVector3 vec_x;
403  TVector3 vec_y;
404 
405  double const switcher = basex.Dot(vec_z);
406  if (std::abs(switcher) <= 0.995) {
407  vec_y = vec_z.Cross(basex).Unit();
408  vec_x = vec_y.Cross(vec_z);
409  }
410  else {
411  // cout << " It switched ! Isn't this lovely !!! " << endl; //
412  vec_y = basez.Cross(vec_z).Unit();
413  vec_x = vec_y.Cross(vec_z);
414  }
415 
416  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
417  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
418  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
419 
420  double const refL = segL.at(i);
421 
422  // Now loop over next segments until thick1 is reached in the next segment
423  for (int j = i; j < tot; j++) {
424  double const L1 = segL.at(j);
425  double const L2 = segL.at(j + 1);
426 
427  double const dz1 = L1 - refL;
428  double const dz2 = L2 - refL;
429 
430  if (dz1 <= thick1 && dz2 > thick1) {
431  double const here_dx = segnx.at(j);
432  double const here_dy = segny.at(j);
433  double const here_dz = segnz.at(j);
434 
435  TVector3 const here_vec{here_dx, here_dy, here_dz};
436  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
437 
438  double const scx = rot_here.X();
439  double const scy = rot_here.Y();
440  double const scz = rot_here.Z();
441 
442  double const azy = find_angle(scz, scy);
443  double const azx = find_angle(scz, scx);
444 
445  constexpr double ULim = 10000.0; // Avoid huge (wrong) angles
446  constexpr double LLim = -10000.0;
447 
448  double const cL = kcal;
449  double const Li = segL.at(i);
450  double const Lj = segL.at(j);
451 
452  if (azy <= ULim && azy >= LLim) { // save scatter in the yz plane
453  ei.push_back(Li * cL); // Energy deposited at i
454  ej.push_back(Lj * cL); // Energy deposited at j
455  th.push_back(azy); // scattered angle
456  ind.push_back(2);
457  }
458 
459  if (azx <= ULim && azx >= LLim) { // save scatter in the za plane
460  ei.push_back(Li * cL);
461  ej.push_back(Lj * cL);
462  th.push_back(azx);
463  ind.push_back(1);
464  }
465 
466  break; // of course !
467  }
468  }
469  }
470 
471  return 0;
472  }
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
a2
#define a2
Definition: HoughBaseAlg.cxx:64
a3
#define a3
Definition: HoughBaseAlg.cxx:65
trkf::TrackMomentumCalculator::find_angle
double find_angle(double vz, double vy) const
Gets angle between two vy and vz.
Definition: TrackMomentumCalculator.cxx:1057
thick
Float_t thick
Definition: plot.C:23
a1
#define a1
Definition: HoughBaseAlg.cxx:63
std::tuple< double, double, double > trkf::TrackMomentumCalculator::getDeltaThetaRMS_ ( Segments const &  segments,
double  thick 
) const
private

Gets the scattered angle RMS for a all segments.

Parameters
segmentssegments computed
thickis the steps_size

TODO: Add better description of steps

Returns
tuple with mean value, rms and error of rms

Definition at line 935 of file TrackMomentumCalculator.cxx.

References a1, a2, a3, find_angle(), trkf::TrackMomentumCalculator::Segments::L, pmtana::mean(), trkf::TrackMomentumCalculator::Segments::nx, trkf::TrackMomentumCalculator::Segments::ny, and trkf::TrackMomentumCalculator::Segments::nz.

Referenced by GetMomentumMultiScatterChi2().

938  {
939  auto const& segnx = segments.nx;
940  auto const& segny = segments.ny;
941  auto const& segnz = segments.nz;
942  auto const& segL = segments.L;
943 
944  int const a1 = segnx.size();
945  int const a2 = segny.size();
946  int const a3 = segnz.size();
947 
948  if (a1 != a2 || a1 != a3) {
949  cout << " ( Get RMS ) Error ! " << endl;
950  return std::make_tuple(0., 0., 0.);
951  }
952 
953  int const tot = a1 - 1;
954 
955  double const thick1 = thick + 0.13;
956 
957  std::vector<float> buf0;
958 
959  for (int i = 0; i < tot; i++) {
960  double const dx = segnx.at(i);
961  double const dy = segny.at(i);
962  double const dz = segnz.at(i);
963 
964  TVector3 const vec_z{dx, dy, dz};
965  TVector3 vec_x;
966  TVector3 vec_y;
967 
968  double const switcher = basex.Dot(vec_z);
969 
970  if (switcher <= 0.995) {
971  vec_y = vec_z.Cross(basex).Unit();
972  vec_x = vec_y.Cross(vec_z);
973  }
974  else {
975  // cout << " It switched ! Isn't this lovely !!! " << endl;
976  vec_y = basez.Cross(vec_z).Unit();
977  vec_x = vec_y.Cross(vec_z);
978  }
979 
980  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
981  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
982  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
983 
984  double const refL = segL.at(i);
985 
986  for (int j = i; j < tot; j++) {
987  double const L1 = segL.at(j);
988  double const L2 = segL.at(j + 1);
989 
990  double const dz1 = L1 - refL;
991  double const dz2 = L2 - refL;
992 
993  if (dz1 <= thick1 && dz2 > thick1) {
994  double const here_dx = segnx.at(j);
995  double const here_dy = segny.at(j);
996  double const here_dz = segnz.at(j);
997 
998  TVector3 const here_vec{here_dx, here_dy, here_dz};
999  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
1000 
1001  double const scx = rot_here.X();
1002  double const scz = rot_here.Z();
1003 
1004  double const azx = find_angle(scz, scx);
1005 
1006  double const ULim = 10000.0;
1007  double const LLim = -10000.0;
1008 
1009  if (azx <= ULim && azx >= LLim) { buf0.push_back(azx); }
1010 
1011  break; // of course !
1012  }
1013  }
1014  }
1015 
1016  int const nmeas = buf0.size();
1017  double nnn = 0.0;
1018 
1019  double mean = 0.0;
1020  double rms = 0.0;
1021  double rmse = 0.0;
1022 
1023  for (int i = 0; i < nmeas; i++) {
1024  mean += buf0.at(i);
1025  nnn++;
1026  }
1027  mean = mean / nnn;
1028 
1029  for (int i = 0; i < nmeas; i++)
1030  rms += ((buf0.at(i)) * (buf0.at(i)));
1031 
1032  rms = rms / (nnn);
1033  rms = std::sqrt(rms);
1034  rmse = rms / std::sqrt(2.0 * tot);
1035 
1036  double rms1 = rms;
1037 
1038  rms = 0.0;
1039 
1040  double ntot1 = 0.0;
1041  double const lev1 = 2.50;
1042 
1043  for (int i = 0; i < nmeas; i++) {
1044  double const amp = buf0.at(i);
1045  if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
1046  ++ntot1;
1047  rms += amp * amp;
1048  }
1049  }
1050 
1051  rms = rms / (ntot1);
1052  rms = std::sqrt(rms);
1053  rmse = rms / std::sqrt(2.0 * ntot1);
1054  return std::make_tuple(mean, rms, rmse);
1055  }
a2
#define a2
Definition: HoughBaseAlg.cxx:64
a3
#define a3
Definition: HoughBaseAlg.cxx:65
trkf::TrackMomentumCalculator::find_angle
double find_angle(double vz, double vy) const
Gets angle between two vy and vz.
Definition: TrackMomentumCalculator.cxx:1057
thick
Float_t thick
Definition: plot.C:23
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
a1
#define a1
Definition: HoughBaseAlg.cxx:63
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2 ( art::Ptr< recob::Track > const &  trk,
const bool  checkValidPoints = false,
const int  maxMomentum_MeV = 7500 
)

Calculate muon momentum (GeV) using multiple coulomb scattering. Chi2 minimization of the Highland formula.

Parameters
trkthe muon track
checkValidPointsrather take into account only valid points or not
maxMomentum_MeVmaximum momentum in MeV for the minimization

TODO: Add better description of the steps done.

Returns
momentum in GeV

Definition at line 474 of file TrackMomentumCalculator.cxx.

References getDeltaThetaRMS_(), getSegTracks_(), recob::Track::HasValidPoint(), recob::Track::LocationAtPoint(), maxLength, pmtana::mean(), n_steps, recob::Track::NumberTrajectoryPoints(), plotRecoTracks_(), steps, and steps_size.

Referenced by trkf::KalmanFilterFinalTrackFitter::setMomValue().

477  {
478  std::vector<float> recoX;
479  std::vector<float> recoY;
480  std::vector<float> recoZ;
481 
482  int n_points = trk->NumberTrajectoryPoints();
483 
484  for (int i = 0; i < n_points; i++) {
485  if (checkValidPoints && !trk->HasValidPoint(i)) continue;
486  auto const& pos = trk->LocationAtPoint(i);
487  recoX.push_back(pos.X());
488  recoY.push_back(pos.Y());
489  recoZ.push_back(pos.Z());
490  }
491 
492  if (recoX.size() < 2) return -1.0;
493 
494  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
495 
496  double const seg_size{steps_size};
497  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
498  if (!segments.has_value()) return -1.0;
499 
500  auto const seg_steps = segments->x.size();
501  if (seg_steps < 2) return -1;
502 
503  double const recoL = segments->L.at(seg_steps - 1);
504  if (recoL < minLength || recoL > maxLength) return -1;
505 
506  double ymax = -999.0;
507  double ymin = +999.0;
508 
509  std::vector<double> xmeas;
510  std::vector<double> ymeas;
511  std::vector<double> eymeas;
512  xmeas.reserve(n_steps);
513  ymeas.reserve(n_steps);
514  eymeas.reserve(n_steps);
515  for (int j = 0; j < n_steps; j++) {
516  double const trial = steps.at(j);
517  auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
518 
519  if (std::isnan(mean) || std::isinf(mean)) {
520  mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
521  continue;
522  }
523  if (std::isnan(rms) || std::isinf(rms)) {
524  mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
525  continue;
526  }
527  if (std::isnan(rmse) || std::isinf(rmse)) {
528  mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
529  continue;
530  }
531 
532  xmeas.push_back(trial); // Is this what is intended?
533  ymeas.push_back(rms);
534  eymeas.push_back(std::sqrt(cet::sum_of_squares(
535  rmse, 0.05 * rms))); // <--- conservative syst. error to fix chi^{2} behaviour !!!
536 
537  if (ymin > rms) ymin = rms;
538  if (ymax < rms) ymax = rms;
539  }
540 
541  assert(xmeas.size() == ymeas.size());
542  assert(xmeas.size() == eymeas.size());
543  if (xmeas.empty()) { return -1.0; }
544 
545  TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
546 
547  gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
548  "thickness in cm; (#Delta#theta)_{rms} in mrad");
549 
550  gr_meas.SetLineColor(kBlack);
551  gr_meas.SetMarkerColor(kBlack);
552  gr_meas.SetMarkerStyle(20);
553  gr_meas.SetMarkerSize(1.2);
554 
555  gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0), steps.at(n_steps - 1) + steps.at(0));
556  gr_meas.SetMinimum(0.0);
557  gr_meas.SetMaximum(1.80 * ymax);
558 
559  ROOT::Minuit2::Minuit2Minimizer mP{};
560  FcnWrapper const wrapper{move(xmeas), move(ymeas), move(eymeas)};
561  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
562 
563  mP.SetFunction(FCA);
564  mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, maxMomentum_MeV / 1.e3);
565  mP.SetLimitedVariable(1, "#delta#theta", 0.0, 1.0, 0.0, 45.0);
566  mP.SetMaxFunctionCalls(1.E9);
567  mP.SetMaxIterations(1.E9);
568  mP.SetTolerance(0.01);
569  mP.SetStrategy(2);
570  mP.SetErrorDef(1.0);
571 
572  bool const mstatus = mP.Minimize();
573 
574  mP.Hesse();
575 
576  const double* pars = mP.X();
577  const double* erpars = mP.Errors();
578 
579  double const deltap = (recoL * kcal) / 2.0;
580 
581  double const p_mcs = pars[0] + deltap;
582  double const p_mcs_e [[maybe_unused]] = erpars[0];
583  return mstatus ? p_mcs : -1.0;
584  }
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
recob::Track::HasValidPoint
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:145
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
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)
Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/...
Definition: TrackMomentumCalculator.cxx:721
trkf::TrackMomentumCalculator::getDeltaThetaRMS_
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Gets the scattered angle RMS for a all segments.
Definition: TrackMomentumCalculator.cxx:935
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:208
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:586
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD ( art::Ptr< recob::Track > const &  trk,
const bool  checkValidPoints = false,
const int  maxMomentum_MeV = 7500,
const int  MomentumStep_MeV = 10,
const int  max_resolution = 0 
)

Calculate muon momentum (GeV) using multiple coulomb scattering by log likelihood.

Parameters
trkthe muon track
checkValidPointsrather take into account only valid points or not
maxMomentum_MeVmaximum momentum in MeV for the minimization
MomentumStep_MeVenergy steps for minimization
max_resolutionmaximum angular resolution for fit. Setting to zero will cause the fit only over momentum and fixed resolution of 2 mrad

TODO: Add better description of the steps done

Returns
momentum in GeV

Definition at line 246 of file TrackMomentumCalculator.cxx.

References e, getDeltaThetaij_(), getSegTracks_(), recob::Track::HasValidPoint(), recob::Track::LocationAtPoint(), maxLength, my_mcs_llhd(), recob::Track::NumberTrajectoryPoints(), plotRecoTracks_(), and steps_size.

251  {
252  std::vector<float> recoX;
253  std::vector<float> recoY;
254  std::vector<float> recoZ;
255 
256  int n_points = trk->NumberTrajectoryPoints();
257 
258  for (int i = 0; i < n_points; i++) {
259  if (checkValidPoints && !trk->HasValidPoint(i)) continue;
260  auto const& pos = trk->LocationAtPoint(i);
261  recoX.push_back(pos.X());
262  recoY.push_back(pos.Y());
263  recoZ.push_back(pos.Z());
264  }
265 
266  if (recoX.size() < 2) return -1.0;
267 
268  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
269 
270  double const seg_size{steps_size};
271 
272  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
273  if (!segments.has_value()) return -1.0;
274 
275  auto const seg_steps = segments->x.size();
276  if (seg_steps < 2) return -1;
277 
278  double const recoL = segments->L.at(seg_steps - 1);
279  if (recoL < minLength || recoL > maxLength) return -1;
280 
281  std::vector<float> dEi;
282  std::vector<float> dEj;
283  std::vector<float> dthij;
284  std::vector<float> ind;
285  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0) return -1.0;
286 
287  double logL = 1e+16;
288  double bf = -666.0; // double errs = -666.0;
289 
290  int const start1{};
291  int const end1{static_cast<int>(maxMomentum_MeV / MomentumStep_MeV)};
292  int const start2{};
293  int const end2{max_resolution}; // 800.0;
294 
295  for (int k = start1; k <= end1; ++k) {
296  double const p_test = 0.001 + k * 0.01;
297 
298  for (int l = start2; l <= end2; ++l) {
299  double const res_test = (start2 == end2) ? 2.0 : 0.001 + l * 1.0; // 0.001+l*1.0;
300  double const fv = my_mcs_llhd(dEi, dEj, dthij, ind, p_test, res_test);
301 
302  if (fv < logL) {
303  bf = p_test;
304  logL = fv;
305  // errs = res_test;
306  }
307  }
308  }
309  return bf;
310  }
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
Gets the scatterd angle for all the segments.
Definition: TrackMomentumCalculator.cxx:371
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
recob::Track::HasValidPoint
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:145
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
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)
Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/...
Definition: TrackMomentumCalculator.cxx:721
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
Minimizer of log likelihood for scattered angle.
Definition: TrackMomentumCalculator.cxx:1114
e
Float_t e
Definition: plot.C:35
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:586
TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const &  trk)

Definition at line 312 of file TrackMomentumCalculator.cxx.

References GetMuMultiScatterLLHD3(), recob::Track::LocationAtPoint(), and recob::Track::NumberTrajectoryPoints().

313  {
314  double const LLHDp = GetMuMultiScatterLLHD3(trk, true);
315  double const LLHDm = GetMuMultiScatterLLHD3(trk, false);
316 
317  if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
318  int const n_points = trk->NumberTrajectoryPoints();
319  return trk->LocationAtPoint<TVector3>(n_points - 1);
320  }
321  else {
322  return trk->LocationAtPoint<TVector3>(0);
323  }
324 
325  // Should never get here
326  return TVector3{};
327  }
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Definition: TrackMomentumCalculator.cxx:329
double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3 ( art::Ptr< recob::Track > const &  trk,
bool  dir 
)

Definition at line 329 of file TrackMomentumCalculator.cxx.

References getDeltaThetaij_(), getSegTracks_(), recob::Track::LocationAtPoint(), maxLength, my_mcs_llhd(), recob::Track::NumberTrajectoryPoints(), and plotRecoTracks_().

Referenced by GetMultiScatterStartingPoint().

331  {
332  std::vector<float> recoX;
333  std::vector<float> recoY;
334  std::vector<float> recoZ;
335 
336  int const n_points = trk->NumberTrajectoryPoints();
337  for (int i = 0; i < n_points; ++i) {
338  auto const index = dir ? i : n_points - 1 - i;
339  auto const& pos = trk->LocationAtPoint(index);
340  recoX.push_back(pos.X());
341  recoY.push_back(pos.Y());
342  recoZ.push_back(pos.Z());
343  }
344 
345  if (recoX.size() < 2) return -1.0;
346 
347  if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
348 
349  constexpr double seg_size{5.0};
350  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
351  if (!segments.has_value()) return -1.0;
352 
353  auto const seg_steps = segments->x.size();
354  if (seg_steps < 2) return -1;
355 
356  double const recoL = segments->L.at(seg_steps - 1);
357  if (recoL < 15.0 || recoL > maxLength) return -1;
358 
359  std::vector<float> dEi;
360  std::vector<float> dEj;
361  std::vector<float> dthij;
362  std::vector<float> ind;
363  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0) return -1.0;
364 
365  double const p_range = recoL * kcal;
366  double const logL = my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
367 
368  return logL;
369  }
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
Gets the scatterd angle for all the segments.
Definition: TrackMomentumCalculator.cxx:371
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
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)
Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/...
Definition: TrackMomentumCalculator.cxx:721
dir
TDirectory * dir
Definition: macro.C:5
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
Minimizer of log likelihood for scattered angle.
Definition: TrackMomentumCalculator.cxx:1114
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:586
std::optional< TrackMomentumCalculator::Segments > trkf::TrackMomentumCalculator::getSegTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz,
double  seg_size 
)
private

Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/P10010.

Parameters
xxx3D reconstructed points x-axis
yyy3D reconstructed points y-axiy
zzz3D reconstructed points z-axiz
seg_sizeSegments size defined in class constructor

TODO: Add better description of steps

Returns
Segments

Definition at line 721 of file TrackMomentumCalculator.cxx.

References a1, a2, a3, compute_max_fluctuation_vector(), gr_seg_xy, gr_seg_xyz, gr_seg_xz, gr_seg_yz, n_seg, seg_stop, x1, x2, x_seg, y1, y2, y_seg, and z_seg.

Referenced by GetMomentumMultiScatterChi2(), GetMomentumMultiScatterLLHD(), and GetMuMultiScatterLLHD3().

726  {
727  double stag = 0.0;
728 
729  int a1 = xxx.size();
730  int a2 = yyy.size();
731  int a3 = zzz.size();
732 
733  if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
734  cout << " ( Digitize reco tacks ) Error ! " << endl;
735  return std::nullopt;
736  }
737 
738  int const stopper = seg_stop / seg_size;
739 
740  std::vector<float> segx, segnx;
741  std::vector<float> segy, segny;
742  std::vector<float> segz, segnz;
743  std::vector<float> segL;
744 
745  int ntot = 0;
746 
747  n_seg = 0;
748 
749  double x0{};
750  double y0{};
751  double z0{};
752 
753  double x00 = xxx.at(0);
754  double y00 = yyy.at(0);
755  double z00 = zzz.at(0);
756 
757  int indC = 0;
758 
759  std::vector<float> vx;
760  std::vector<float> vy;
761  std::vector<float> vz;
762 
763  for (int i = 0; i < a1; i++) {
764  x0 = xxx.at(i);
765  y0 = yyy.at(i);
766  z0 = zzz.at(i);
767 
768  double const RR0 = std::sqrt(cet::sum_of_squares(x00 - x0, y00 - y0, z00 - z0));
769 
770  if (RR0 >= stag) {
771  segx.push_back(x0);
772  segy.push_back(y0);
773  segz.push_back(z0);
774 
775  segL.push_back(stag);
776 
777  x_seg[n_seg] = x0;
778  y_seg[n_seg] = y0;
779  z_seg[n_seg] = z0;
780 
781  n_seg++;
782 
783  vx.push_back(x0);
784  vy.push_back(y0);
785  vz.push_back(z0);
786 
787  ntot++;
788 
789  indC = i + 1;
790 
791  break;
792  }
793  }
794 
795  for (int i = indC; i < a1 - 1; i++) {
796  double const x1 = xxx.at(i);
797  double const y1 = yyy.at(i);
798  double const z1 = zzz.at(i);
799 
800  double const dr1 = std::sqrt(cet::sum_of_squares(x1 - x0, y1 - y0, z1 - z0));
801 
802  double const x2 = xxx.at(i + 1);
803  double const y2 = yyy.at(i + 1);
804  double const z2 = zzz.at(i + 1);
805 
806  double const dr2 = std::sqrt(cet::sum_of_squares(x2 - x0, y2 - y0, z2 - z0));
807 
808  if (dr1 < seg_size) {
809  vx.push_back(x1);
810  vy.push_back(y1);
811  vz.push_back(z1);
812 
813  ntot++;
814  }
815 
816  if (dr1 <= seg_size && dr2 > seg_size) {
817  double const dx = x2 - x1;
818  double const dy = y2 - y1;
819  double const dz = z2 - z1;
820  double const dr = std::sqrt(dx * dx + dy * dy + dz * dz);
821 
822  if (dr == 0) {
823  cout << " ( Zero ) Error ! " << endl;
824  return std::nullopt;
825  }
826 
827  double const beta = 2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
828 
829  double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
830  double const delta = beta * beta - 4.0 * gamma;
831 
832  if (delta < 0.0) {
833  cout << " ( Discriminant ) Error ! " << endl;
834  return std::nullopt;
835  }
836 
837  double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
838  double const t = lysi1;
839 
840  double const xp = x1 + t * dx;
841  double const yp = y1 + t * dy;
842  double const zp = z1 + t * dz;
843 
844  segx.push_back(xp);
845  segy.push_back(yp);
846  segz.push_back(zp);
847 
848  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
849 
850  x_seg[n_seg] = xp;
851  y_seg[n_seg] = yp;
852  z_seg[n_seg] = zp;
853  n_seg++;
854 
855  x0 = xp;
856  y0 = yp;
857  z0 = zp;
858 
859  vx.push_back(x0);
860  vy.push_back(y0);
861  vz.push_back(z0);
862 
863  ntot++;
864  if (n_seg <= 1) // This should never happen
865  return std::nullopt;
866 
867  compute_max_fluctuation_vector(segx, segy, segz, segnx, segny, segnz, vx, vy, vz);
868 
869  ntot = 1;
870  vx.push_back(x0);
871  vy.push_back(y0);
872  vz.push_back(z0);
873  }
874  else if (dr1 > seg_size) {
875 
876  i = (i - 1); // In this case, we keep at the same point until seg_size is reached
877  double const dx = x1 - x0;
878  double const dy = y1 - y0;
879  double const dz = z1 - z0;
880  double const dr = std::sqrt(cet::sum_of_squares(dx, dy, dz));
881 
882  if (dr == 0) {
883  cout << " ( Zero ) Error ! " << endl;
884  return std::nullopt;
885  }
886 
887  double const t = seg_size / dr;
888  double const xp = x0 + t * dx;
889  double const yp = y0 + t * dy;
890  double const zp = z0 + t * dz;
891 
892  segx.push_back(xp);
893  segy.push_back(yp);
894  segz.push_back(zp);
895  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
896 
897  x_seg[n_seg] = xp;
898  y_seg[n_seg] = yp;
899  z_seg[n_seg] = zp;
900  n_seg++;
901 
902  x0 = xp;
903  y0 = yp;
904  z0 = zp;
905 
906  vx.push_back(x0);
907  vy.push_back(y0);
908  vz.push_back(z0);
909 
910  ntot++;
911  if (n_seg <= 1) // This should never happen
912  return std::nullopt;
913 
914  compute_max_fluctuation_vector(segx, segy, segz, segnx, segny, segnz, vx, vy, vz);
915 
916  // vectors are cleared in previous step
917  ntot = 1;
918  vx.push_back(x0);
919  vy.push_back(y0);
920  vz.push_back(z0);
921  }
922 
923  if (n_seg >= (stopper + 1.0) && seg_stop != -1) break;
924  }
925 
926  delete gr_seg_xyz;
927  gr_seg_xyz = new TPolyLine3D{n_seg, z_seg, x_seg, y_seg};
928  gr_seg_yz = TGraph{n_seg, z_seg, y_seg};
929  gr_seg_xz = TGraph{n_seg, z_seg, x_seg};
930  gr_seg_xy = TGraph{n_seg, x_seg, y_seg};
931 
932  return std::make_optional<Segments>(Segments{segx, segnx, segy, segny, segz, segnz, segL});
933  }
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
trkf::TrackMomentumCalculator::compute_max_fluctuation_vector
void compute_max_fluctuation_vector(const std::vector< float > segx, const std::vector< float > segy, const std::vector< float > segz, std::vector< float > &segnx, std::vector< float > &segny, std::vector< float > &segnz, std::vector< float > &vx, std::vector< float > &vy, std::vector< float > &vz)
Computes the vector with most scattering inside a segment with size steps_size.
Definition: TrackMomentumCalculator.cxx:616
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
a2
#define a2
Definition: HoughBaseAlg.cxx:64
a3
#define a3
Definition: HoughBaseAlg.cxx:65
x2
Float_t x2[n_points_geant4]
Definition: compare.C:26
a1
#define a1
Definition: HoughBaseAlg.cxx:63
double trkf::TrackMomentumCalculator::GetTrackMomentum ( double  trkrange,
int  pdg 
) const

Definition at line 139 of file TrackMomentumCalculator.cxx.

References util::abs(), and E.

Referenced by trkmkr::KalmanFilterFitTrackMaker::getMomentum(), trkf::KalmanFilterTrajectoryFitter::setMomValue(), and trkf::KalmanFilterFinalTrackFitter::setMomValue().

140  {
141  /* Muon range-momentum tables from CSDA (Argon density = 1.4 g/cm^3)
142  website:
143  http://pdg.lbl.gov/2012/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_289.pdf
144 
145  CSDA table values:
146  float Range_grampercm[30] = {9.833E-1, 1.786E0, 3.321E0,
147  6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2,
148  2.385E2, 4.934E2, 6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
149  4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4,
150  4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5}; float KE_MeV[30] = {10, 14,
151  20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000,
152  4000, 8000, 10000, 14000, 20000, 30000, 40000, 80000, 100000, 140000,
153  200000, 300000, 400000};
154 
155  Functions below are obtained by fitting polynomial fits to KE_MeV vs
156  Range (cm) graph. A better fit was obtained by splitting the graph into
157  two: Below range<=200cm,a polynomial of power 4 was a good fit; above
158  200cm, a polynomial of power 6 was a good fit
159 
160  Fit errors for future purposes:
161  Below 200cm, Forpoly4 fit: p0 err=1.38533;p1 err=0.209626; p2
162  err=0.00650077; p3 err=6.42207E-5; p4 err=1.94893E-7; Above 200cm,
163  Forpoly6 fit: p0 err=5.24743;p1 err=0.0176229; p2 err=1.6263E-5; p3
164  err=5.9155E-9; p4 err=9.71709E-13; p5 err=7.22381E-17;p6
165  err=1.9709E-21;*/
166 
168  //*********For muon, the calculations are valid up to 1.91E4 cm range
169  //corresponding to a Muon KE of 40 GeV**********//
171 
172  /*Proton range-momentum tables from CSDA (Argon density = 1.4 g/cm^3):
173  website: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html
174 
175  CSDA values:
176  double KE_MeV_P_Nist[31]={10, 15, 20, 30, 40, 80, 100, 150, 200, 250, 300,
177  350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
178  1500, 2000, 2500, 3000, 4000, 5000};
179 
180  double Range_gpercm_P_Nist[31]={1.887E-1,3.823E-1, 6.335E-1, 1.296,
181  2.159, 7.375, 1.092E1, 2.215E1, 3.627E1, 5.282E1, 7.144E1,
182  9.184E1, 1.138E2, 1.370E2, 1.614E2, 1.869E2, 2.132E2, 2.403E2,
183  2.681E2, 2.965E2, 3.254E2, 3.548E2, 3.846E2, 4.148E2, 4.454E2,
184  7.626E2, 1.090E3, 1.418E3, 1.745E3, 2.391E3, 3.022E3};
185 
186  Functions below are obtained by fitting power and polynomial fits to
187  KE_MeV vs Range (cm) graph. A better fit was obtained by splitting the
188  graph into two: Below range<=80cm,a a*(x^b) was a good fit; above 80cm, a
189  polynomial of power 6 was a good fit
190 
191  Fit errors for future purposes:
192  For power function fit: a=0.388873; and b=0.00347075
193  Forpoly6 fit: p0 err=3.49729;p1 err=0.0487859; p2 err=0.000225834; p3
194  err=4.45542E-7; p4 err=4.16428E-10; p5 err=1.81679E-13;p6
195  err=2.96958E-17;*/
196 
198  //*********For proton, the calculations are valid up to 3.022E3 cm range
199  //corresponding to a Muon KE of 5 GeV**********//
201 
202  if (trkrange < 0 || std::isnan(trkrange)) {
203  mf::LogError("TrackMomentumCalculator") << "Invalid track range " << trkrange << " return -1";
204  return -1.;
205  }
206 
207  double KE, Momentum, M;
208  constexpr double Muon_M = 105.7, Proton_M = 938.272;
209 
210  if (abs(pdg) == 13) {
211  M = Muon_M;
212  KE = KEvsR_spline3.Eval(trkrange);
213  }
214  else if (abs(pdg) == 2212) {
215  M = Proton_M;
216  if (trkrange > 0 && trkrange <= 80)
217  KE = 29.9317 * std::pow(trkrange, 0.586304);
218  else if (trkrange > 80 && trkrange <= 3.022E3)
219  KE = 149.904 + (3.34146 * trkrange) + (-0.00318856 * trkrange * trkrange) +
220  (4.34587E-6 * trkrange * trkrange * trkrange) +
221  (-3.18146E-9 * trkrange * trkrange * trkrange * trkrange) +
222  (1.17854E-12 * trkrange * trkrange * trkrange * trkrange * trkrange) +
223  (-1.71763E-16 * trkrange * trkrange * trkrange * trkrange * trkrange * trkrange);
224  else
225  KE = -999;
226  }
227  else
228  KE = -999;
229 
230  if (KE < 0)
231  Momentum = -999;
232  else
233  Momentum = std::sqrt((KE * KE) + (2 * M * KE));
234 
235  Momentum = Momentum / 1000;
236 
237  return Momentum;
238  }
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
E
Float_t E
Definition: plot.C:20
double trkf::TrackMomentumCalculator::my_g ( double  xx,
double  Q,
double  s 
) const
private

chi square minizer using Minuit2, it will minize (xx-Q)/s

Parameters
xx
Q
s
Returns
Momentum in GeV

Definition at line 1097 of file TrackMomentumCalculator.cxx.

Referenced by my_mcs_llhd().

1098  {
1099  if (s == 0.) {
1100  cout << " Error : The code tries to divide by zero ! " << endl;
1101  return 0.;
1102  }
1103 
1104  double const arg = (xx - Q) / s;
1105  double const result = -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1106 
1107  if (std::isnan(result) || std::isinf(result)) {
1108  cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
1109  }
1110 
1111  return result;
1112  }
xx
Double_t xx
Definition: macro.C:12
double trkf::TrackMomentumCalculator::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
private

Minimizer of log likelihood for scattered angle.

Parameters
dEienergy at step i
dEjenergy at step j
dthijscattered angle between points i and j
indselection of scattering plane
x0momentum to be fitted
x1resolution to be fitted

TODO: Add better description of steps

Returns
momentum in GeV

Definition at line 1114 of file TrackMomentumCalculator.cxx.

References util::abs(), my_g(), rad_length, steps_size, and x1.

Referenced by GetMomentumMultiScatterLLHD(), and GetMuMultiScatterLLHD3().

1120  {
1121  double p = x0;
1122  double theta0x = x1;
1123 
1124  double result = 0.0;
1125  double nnn1 = dEi.size(); // number of segments of energy
1126 
1127  double red_length = (steps_size) / rad_length;
1128  double addth = 0;
1129 
1130  for (int i = 0; i < nnn1; i++) {
1131  double Ei = p - dEi.at(i); // Estimated energy at point i
1132  double Ej = p - dEj.at(i); // Estimated enery at point j
1133 
1134  // If the momentum p choosen allows that the muon stopped inside, add 1 rad to the change in scatter angle (as the particle stops)
1135  if (Ei > 0 && Ej < 0) addth = 3.14 * 1000.0;
1136 
1137  Ei = std::abs(Ei);
1138  Ej = std::abs(Ej);
1139 
1140  // Highland formula
1141  // Parameters given at Particle Data Group https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-passage-particles-matter.pdf
1142  double tH0 =
1143  (13.6 / std::sqrt(Ei * Ej)) * (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1144 
1145  double rms = -1.0;
1146 
1147  if (ind.at(i) == 1) {
1148  // Computes the rms of angle
1149  rms = std::sqrt(tH0 * tH0 + cet::square(theta0x));
1150 
1151  double const DT = dthij.at(i) + addth;
1152  double const prob = my_g(DT, 0.0, rms); // Computes log likelihood
1153 
1154  result = result - 2.0 * prob; // Adds for each segment
1155  }
1156  }
1157 
1158  if (std::isnan(result) || std::isinf(result)) {
1159  std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
1160  }
1161  return result;
1162  }
trkf::TrackMomentumCalculator::my_g
double my_g(double xx, double Q, double s) const
chi square minizer using Minuit2, it will minize (xx-Q)/s
Definition: TrackMomentumCalculator.cxx:1097
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
bool trkf::TrackMomentumCalculator::plotRecoTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz 
)
private

Definition at line 586 of file TrackMomentumCalculator.cxx.

References gr_reco_xy, gr_reco_xyz, gr_reco_xz, gr_reco_yz, and n.

Referenced by GetMomentumMultiScatterChi2(), GetMomentumMultiScatterLLHD(), and GetMuMultiScatterLLHD3().

589  {
590  auto const n = xxx.size();
591  auto const y_size = yyy.size();
592  auto const z_size = zzz.size();
593 
594  if (n != y_size || n != z_size) {
595  cout << " ( Get reco tacks ) Error ! " << endl;
596  return false;
597  }
598 
599  // Here, we perform a const-cast to float* because, sadly,
600  // TPolyLine3D requires a pointer to a non-const object. We will
601  // trust that ROOT does not mess around with the underlying data.
602  auto xs = const_cast<float*>(xxx.data());
603  auto ys = const_cast<float*>(yyy.data());
604  auto zs = const_cast<float*>(zzz.data());
605 
606  auto const narrowed_size = static_cast<int>(n); // ROOT requires a signed integral type
607  delete gr_reco_xyz;
608  gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
609  gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
610  gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
611  gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
612 
613  return true;
614  }
n
Char_t n[5]
Definition: comparison_ascii.C:55

Member Data Documentation

TGraph trkf::TrackMomentumCalculator::gr_reco_xy {}
private

Definition at line 224 of file TrackMomentumCalculator.h.

Referenced by plotRecoTracks_().

TPolyLine3D* trkf::TrackMomentumCalculator::gr_reco_xyz {nullptr}
private

Definition at line 223 of file TrackMomentumCalculator.h.

Referenced by plotRecoTracks_().

TGraph trkf::TrackMomentumCalculator::gr_reco_xz {}
private

Definition at line 226 of file TrackMomentumCalculator.h.

Referenced by plotRecoTracks_().

TGraph trkf::TrackMomentumCalculator::gr_reco_yz {}
private

Definition at line 225 of file TrackMomentumCalculator.h.

Referenced by plotRecoTracks_().

TGraph trkf::TrackMomentumCalculator::gr_seg_xy {}
private

Definition at line 229 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

TPolyLine3D* trkf::TrackMomentumCalculator::gr_seg_xyz {nullptr}
private

Definition at line 228 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

TGraph trkf::TrackMomentumCalculator::gr_seg_xz {}
private

Definition at line 231 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

TGraph trkf::TrackMomentumCalculator::gr_seg_yz {}
private

Definition at line 230 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

double trkf::TrackMomentumCalculator::maxLength
private

Definition at line 209 of file TrackMomentumCalculator.h.

Referenced by GetMomentumMultiScatterChi2(), GetMomentumMultiScatterLLHD(), GetMuMultiScatterLLHD3(), and TrackMomentumCalculator().

double trkf::TrackMomentumCalculator::minLength
private

Definition at line 208 of file TrackMomentumCalculator.h.

int trkf::TrackMomentumCalculator::n_seg {}
private

Definition at line 189 of file TrackMomentumCalculator.h.

Referenced by compute_max_fluctuation_vector(), and getSegTracks_().

int trkf::TrackMomentumCalculator::n_steps {6}
private

Definition at line 205 of file TrackMomentumCalculator.h.

Referenced by GetMomentumMultiScatterChi2(), and TrackMomentumCalculator().

double trkf::TrackMomentumCalculator::rad_length {14.0}
private

Definition at line 211 of file TrackMomentumCalculator.h.

Referenced by my_mcs_llhd().

float trkf::TrackMomentumCalculator::seg_stop {-1.}
private

Definition at line 188 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

std::vector<float> trkf::TrackMomentumCalculator::steps
private

Definition at line 206 of file TrackMomentumCalculator.h.

Referenced by GetMomentumMultiScatterChi2(), and TrackMomentumCalculator().

double trkf::TrackMomentumCalculator::steps_size
private

Definition at line 210 of file TrackMomentumCalculator.h.

Referenced by GetMomentumMultiScatterChi2(), GetMomentumMultiScatterLLHD(), my_mcs_llhd(), and TrackMomentumCalculator().

float trkf::TrackMomentumCalculator::x_seg[50000]
private

Definition at line 191 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

float trkf::TrackMomentumCalculator::y_seg[50000]
private

Definition at line 192 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

float trkf::TrackMomentumCalculator::z_seg[50000]
private

Definition at line 193 of file TrackMomentumCalculator.h.

Referenced by getSegTracks_().

The documentation for this class was generated from the following files: