All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::TrackMomentumCalculator Class Reference

#include <TrackMomentumCalculator.h>

Classes

struct  Segments
 

Public Member Functions

 TrackMomentumCalculator (double minLength=100.0, double maxLength=1350.0)
 
double GetTrackMomentum (double trkrange, int pdg) const
 
double GetMomentumMultiScatterChi2 (art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false)
 
double GetMomentumMultiScatterLLHD (art::Ptr< recob::Track > const &trk)
 
double GetMuMultiScatterLLHD3 (art::Ptr< recob::Track > const &trk, bool dir)
 
TVector3 GetMultiScatterStartingPoint (art::Ptr< recob::Track > const &trk)
 
 TrackMomentumCalculator (double minLength=100.0, double maxLength=1350.0)
 
double GetTrackMomentum (double trkrange, int pdg) const
 
double GetMomentumMultiScatterChi2 (art::Ptr< recob::Track > const &trk)
 
double GetMomentumMultiScatterLLHD (art::Ptr< recob::Track > const &trk)
 
double GetMuMultiScatterLLHD3 (art::Ptr< recob::Track > const &trk, bool dir)
 
TVector3 GetMultiScatterStartingPoint (art::Ptr< recob::Track > const &trk)
 
 TrackMomentumCalculator (double minLength=100.0, double maxLength=1350.0)
 
double GetTrackMomentum (double trkrange, int pdg) const
 
double GetMomentumMultiScatterChi2 (art::Ptr< recob::Track > const &trk)
 
double GetMomentumMultiScatterLLHD (art::Ptr< recob::Track > const &trk)
 
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)
 
std::optional< Segments > getSegTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
 
std::tuple< double, double,
double > 		getDeltaThetaRMS_ (Segments const &segments, double thick) const
 
int getDeltaThetaij_ (std::vector< float > &ei, std::vector< float > &ej, std::vector< float > &th, std::vector< float > &ind, Segments const &segments, double thick) const
 
double my_g (double xx, double Q, double s) const
 
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
 
double find_angle (double vz, double vy) const
 
bool plotRecoTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
 
std::optional< Segments > getSegTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
 
std::tuple< double, double,
double > 		getDeltaThetaRMS_ (Segments const &segments, double thick) const
 
int getDeltaThetaij_ (std::vector< float > &ei, std::vector< float > &ej, std::vector< float > &th, std::vector< float > &ind, Segments const &segments, double thick) const
 
double my_g (double xx, double Q, double s) const
 
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
 
double find_angle (double vz, double vy) const
 
bool plotRecoTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
 
std::optional< Segments > getSegTracks_ (std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
 
std::tuple< double, double,
double > 		getDeltaThetaRMS_ (Segments const &segments, double thick) const
 
int getDeltaThetaij_ (std::vector< float > &ei, std::vector< float > &ej, std::vector< float > &th, std::vector< float > &ind, Segments const &segments, double thick) const
 
double my_g (double xx, double Q, double s) const
 
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
 
double find_angle (double vz, double vy) const
 

Private Attributes

float seg_stop {-1.}
 
int n_seg {}
 
float x_seg [50000]
 
float y_seg [50000]
 
float z_seg [50000]
 
float steps_size {10.}
 
int n_steps {6}
 
std::vector< float > steps
 
double minLength
 
double maxLength
 
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 larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 133 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

135  : minLength{min}
136  , maxLength{max}
137  {
138  for (int i = 1; i <= n_steps; i++) {
139  steps.push_back(steps_size * i);
140  }
141  }
trkf::TrackMomentumCalculator::TrackMomentumCalculator ( double  minLength = 100.0,
double  maxLength = 1350.0 
)
trkf::TrackMomentumCalculator::TrackMomentumCalculator ( double  minLength = 100.0,
double  maxLength = 1350.0 
)

Member Function Documentation

double trkf::TrackMomentumCalculator::find_angle ( double  vz,
double  vy 
) const
private

Definition at line 1181 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

1182  {
1183  double thetayz = -999.0;
1184 
1185  if (vz > 0 && vy > 0) {
1186  double ratio = std::abs(vy / vz);
1187  thetayz = std::atan(ratio);
1188  }
1189 
1190  else if (vz < 0 && vy > 0) {
1191  double ratio = std::abs(vy / vz);
1192  thetayz = std::atan(ratio);
1193  thetayz = TMath::Pi() - thetayz;
1194  }
1195 
1196  else if (vz < 0 && vy < 0) {
1197  double ratio = std::abs(vy / vz);
1198  thetayz = std::atan(ratio);
1199  thetayz = thetayz + TMath::Pi();
1200  }
1201 
1202  else if (vz > 0 && vy < 0) {
1203  double ratio = std::abs(vy / vz);
1204  thetayz = std::atan(ratio);
1205  thetayz = 2.0 * TMath::Pi() - thetayz;
1206  }
1207 
1208  else if (vz == 0 && vy > 0) {
1209  thetayz = TMath::Pi() / 2.0;
1210  }
1211 
1212  else if (vz == 0 && vy < 0) {
1213  thetayz = 3.0 * TMath::Pi() / 2.0;
1214  }
1215 
1216  if (thetayz > TMath::Pi()) {
1217  thetayz = thetayz - 2.0 * TMath::Pi();
1218  }
1219 
1220  return 1000.0 * thetayz;
1221  }
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
double trkf::TrackMomentumCalculator::find_angle ( double  vz,
double  vy 
) const
private
double trkf::TrackMomentumCalculator::find_angle ( double  vz,
double  vy 
) const
private
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

Definition at line 390 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

396  {
397  int const a1 = segments.x.size();
398  int const a2 = segments.y.size();
399  int const a3 = segments.z.size();
400 
401  if (a1 != a2 || a1 != a3) {
402  std::cout << " ( Get thij ) Error ! " << std::endl;
403  return -1.0;
404  }
405 
406  auto const& segnx = segments.nx;
407  auto const& segny = segments.ny;
408  auto const& segnz = segments.nz;
409  auto const& segL = segments.L;
410 
411  int tot = a1 - 1;
412  double thick1 = thick + 0.13;
413 
414  for (int i = 0; i < tot; i++) {
415  double const dx = segnx.at(i);
416  double const dy = segny.at(i);
417  double const dz = segnz.at(i);
418 
419  TVector3 const vec_z{dx, dy, dz};
420  TVector3 vec_x;
421  TVector3 vec_y;
422 
423  double const switcher = basex.Dot(vec_z);
424  if (std::abs(switcher) <= 0.995) {
425  vec_y = vec_z.Cross(basex).Unit();
426  vec_x = vec_y.Cross(vec_z);
427  }
428  else {
429  // cout << " It switched ! Isn't this lovely !!! " << endl; //
430  vec_y = basez.Cross(vec_z).Unit();
431  vec_x = vec_y.Cross(vec_z);
432  }
433 
434  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
435  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
436  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
437 
438  double const refL = segL.at(i);
439 
440  for (int j = i; j < tot; j++) {
441  double const L1 = segL.at(j);
442  double const L2 = segL.at(j + 1);
443 
444  double const dz1 = L1 - refL;
445  double const dz2 = L2 - refL;
446 
447  if (dz1 <= thick1 && dz2 > thick1) {
448  double const here_dx = segnx.at(j);
449  double const here_dy = segny.at(j);
450  double const here_dz = segnz.at(j);
451 
452  TVector3 const here_vec{here_dx, here_dy, here_dz};
453  TVector3 const rot_here{
454  Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
455 
456  double const scx = rot_here.X();
457  double const scy = rot_here.Y();
458  double const scz = rot_here.Z();
459 
460  double const azy = find_angle(scz, scy);
461  double const azx = find_angle(scz, scx);
462 
463  constexpr double ULim = 10000.0;
464  constexpr double LLim = -10000.0;
465 
466  double const cL = kcal;
467  double const Li = segL.at(i);
468  double const Lj = segL.at(j);
469 
470  if (azy <= ULim && azy >= LLim) {
471  ei.push_back(Li * cL);
472  ej.push_back(Lj * cL);
473  th.push_back(azy);
474  ind.push_back(2);
475  }
476 
477  if (azx <= ULim && azx >= LLim) {
478  ei.push_back(Li * cL);
479  ej.push_back(Lj * cL);
480  th.push_back(azx);
481  ind.push_back(1);
482  }
483 
484  break; // of course !
485  }
486  }
487  }
488 
489  return 0;
490  }
a2
#define a2
Definition: HoughBaseAlg.cxx:64
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
a3
#define a3
Definition: HoughBaseAlg.cxx:65
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
a1
#define a1
Definition: HoughBaseAlg.cxx:63
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
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
std::tuple< double, double, double > trkf::TrackMomentumCalculator::getDeltaThetaRMS_ ( Segments const &  segments,
double  thick 
) const
private

Definition at line 1053 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

1055  {
1056  auto const& segnx = segments.nx;
1057  auto const& segny = segments.ny;
1058  auto const& segnz = segments.nz;
1059  auto const& segL = segments.L;
1060 
1061  int const a1 = segnx.size();
1062  int const a2 = segny.size();
1063  int const a3 = segnz.size();
1064 
1065  if (a1 != a2 || a1 != a3) {
1066  cout << " ( Get RMS ) Error ! " << endl;
1067  return std::make_tuple(0., 0., 0.);
1068  }
1069 
1070  int const tot = a1 - 1;
1071 
1072  double const thick1 = thick + 0.13;
1073 
1074  std::vector<float> buf0;
1075 
1076  for (int i = 0; i < tot; i++) {
1077  double const dx = segnx.at(i);
1078  double const dy = segny.at(i);
1079  double const dz = segnz.at(i);
1080 
1081  TVector3 const vec_z{dx, dy, dz};
1082  TVector3 vec_x;
1083  TVector3 vec_y;
1084 
1085  double const switcher = basex.Dot(vec_z);
1086 
1087  if (switcher <= 0.995) {
1088  vec_y = vec_z.Cross(basex).Unit();
1089  vec_x = vec_y.Cross(vec_z);
1090  }
1091  else {
1092  // cout << " It switched ! Isn't this lovely !!! " << endl;
1093  vec_y = basez.Cross(vec_z).Unit();
1094  vec_x = vec_y.Cross(vec_z);
1095  }
1096 
1097  TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
1098  TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
1099  TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
1100 
1101  double const refL = segL.at(i);
1102 
1103  for (int j = i; j < tot; j++) {
1104  double const L1 = segL.at(j);
1105  double const L2 = segL.at(j + 1);
1106 
1107  double const dz1 = L1 - refL;
1108  double const dz2 = L2 - refL;
1109 
1110  if (dz1 <= thick1 && dz2 > thick1) {
1111  double const here_dx = segnx.at(j);
1112  double const here_dy = segny.at(j);
1113  double const here_dz = segnz.at(j);
1114 
1115  TVector3 const here_vec{here_dx, here_dy, here_dz};
1116  TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
1117 
1118  double const scx = rot_here.X();
1119  double const scy = rot_here.Y();
1120  double const scz = rot_here.Z();
1121 
1122  double azy = find_angle(scz, scy);
1123  azy *= 1.0;
1124 
1125  double const azx = find_angle(scz, scx);
1126 
1127  double const ULim = 10000.0;
1128  double const LLim = -10000.0;
1129 
1130  if (azx <= ULim && azx >= LLim) {
1131  buf0.push_back(azx);
1132  }
1133 
1134  break; // of course !
1135  }
1136  }
1137  }
1138 
1139  int const nmeas = buf0.size();
1140  double nnn = 0.0;
1141 
1142  double mean = 0.0;
1143  double rms = 0.0;
1144  double rmse = 0.0;
1145 
1146  for (int i = 0; i < nmeas; i++) {
1147  mean += buf0.at(i);
1148  nnn++;
1149  }
1150  mean = mean / nnn;
1151 
1152  for (int i = 0; i < nmeas; i++)
1153  rms += ((buf0.at(i)) * (buf0.at(i)));
1154 
1155  rms = rms / (nnn);
1156  rms = std::sqrt(rms);
1157  rmse = rms / std::sqrt(2.0 * tot);
1158 
1159  double rms1 = rms;
1160 
1161  rms = 0.0;
1162 
1163  double ntot1 = 0.0;
1164  double const lev1 = 2.50;
1165 
1166  for (int i = 0; i < nmeas; i++) {
1167  double const amp = buf0.at(i);
1168  if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
1169  ++ntot1;
1170  rms += amp * amp;
1171  }
1172  }
1173 
1174  rms = rms / (ntot1);
1175  rms = std::sqrt(rms);
1176  rmse = rms / std::sqrt(2.0 * ntot1);
1177  return std::make_tuple(mean, rms, rmse);
1178  }
a2
#define a2
Definition: HoughBaseAlg.cxx:64
a3
#define a3
Definition: HoughBaseAlg.cxx:65
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
a1
#define a1
Definition: HoughBaseAlg.cxx:63
std::tuple<double, double, double> trkf::TrackMomentumCalculator::getDeltaThetaRMS_ ( Segments const &  segments,
double  thick 
) const
private
std::tuple<double, double, double> trkf::TrackMomentumCalculator::getDeltaThetaRMS_ ( Segments const &  segments,
double  thick 
) const
private
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2 ( art::Ptr< recob::Track > const &  trk,
const bool  checkValidPoints = false 
)

Definition at line 493 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

495  {
496  std::vector<float> recoX;
497  std::vector<float> recoY;
498  std::vector<float> recoZ;
499 
500  int n_points = trk->NumberTrajectoryPoints();
501 
502  for (int i = 0; i < n_points; i++) {
503  if (checkValidPoints && !trk->HasValidPoint(i)) continue;
504  auto const& pos = trk->LocationAtPoint(i);
505  recoX.push_back(pos.X());
506  recoY.push_back(pos.Y());
507  recoZ.push_back(pos.Z());
508  }
509 
510  if (recoX.size() < 2)
511  return -1.0;
512 
513  if (!plotRecoTracks_(recoX, recoY, recoZ))
514  return -1.0;
515 
516  double const seg_size{steps_size};
517  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
518  if (!segments.has_value())
519  return -1.0;
520 
521  auto const seg_steps = segments->x.size();
522  if (seg_steps < 2)
523  return -1;
524 
525  double const recoL = segments->L.at(seg_steps - 1);
526  if (recoL < minLength || recoL > maxLength)
527  return -1;
528 
529  double ymax = -999.0;
530  double ymin = +999.0;
531 
532  std::vector<double> xmeas;
533  std::vector<double> ymeas;
534  std::vector<double> eymeas;
535  xmeas.reserve(n_steps);
536  ymeas.reserve(n_steps);
537  eymeas.reserve(n_steps);
538  for (int j = 0; j < n_steps; j++) {
539  double const trial = steps.at(j);
540  auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
541 
542  if (std::isnan(mean) || std::isinf(mean)) {
543  mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
544  continue;
545  }
546  if (std::isnan(rms) || std::isinf(rms)) {
547  mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
548  continue;
549  }
550  if (std::isnan(rmse) || std::isinf(rmse)) {
551  mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
552  continue;
553  }
554 
555  xmeas.push_back(trial); // Is this what is intended?
556  ymeas.push_back(rms);
557  eymeas.push_back(std::sqrt(cet::sum_of_squares(rmse, 0.05 * rms))); // <--- conservative syst. error to fix chi^{2} behaviour !!!
558 
559  if (ymin > rms)
560  ymin = rms;
561  if (ymax < rms)
562  ymax = rms;
563  }
564 
565  assert(xmeas.size() == ymeas.size());
566  assert(xmeas.size() == eymeas.size());
567  if (xmeas.empty()) {
568  return -1.0;
569  }
570 
571  TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
572 
573  gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
574  "thickness in cm; (#Delta#theta)_{rms} in mrad");
575 
576  gr_meas.SetLineColor(kBlack);
577  gr_meas.SetMarkerColor(kBlack);
578  gr_meas.SetMarkerStyle(20);
579  gr_meas.SetMarkerSize(1.2);
580 
581  gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0),
582  steps.at(n_steps - 1) + steps.at(0));
583  gr_meas.SetMinimum(0.0);
584  gr_meas.SetMaximum(1.80 * ymax);
585 
586  ROOT::Minuit2::Minuit2Minimizer mP{};
587  FcnWrapper const wrapper{move(xmeas), move(ymeas), move(eymeas)};
588  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
589 
590  mP.SetFunction(FCA);
591  mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, 7.5);
592  mP.SetLimitedVariable(1, "#delta#theta", 0.0, 1.0, 0.0, 45.0);
593  mP.SetMaxFunctionCalls(1.E9);
594  mP.SetMaxIterations(1.E9);
595  mP.SetTolerance(0.01);
596  mP.SetStrategy(2);
597  mP.SetErrorDef(1.0);
598 
599  bool const mstatus = mP.Minimize();
600 
601  mP.Hesse();
602 
603  const double* pars = mP.X();
604  const double* erpars = mP.Errors();
605 
606  double const deltap = (recoL * kcal) / 2.0;
607 
608  double const p_mcs = pars[0] + deltap;
609  double const p_mcs_e [[maybe_unused]] = erpars[0];
610  return mstatus ? p_mcs : -1.0;
611  }
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:646
trkf::TrackMomentumCalculator::getDeltaThetaRMS_
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:1053
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:614
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2 ( art::Ptr< recob::Track > const &  trk)
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2 ( art::Ptr< recob::Track > const &  trk)

Definition at line 469 of file sbnana/sbnanalysis/ana/SBNOscReco/LArReco/TrackMomentumCalculator.cxx.

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  }
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:646
trkf::TrackMomentumCalculator::getDeltaThetaRMS_
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:1053
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:614
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD ( art::Ptr< recob::Track > const &  trk)

Definition at line 253 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

255  {
256  std::vector<float> recoX;
257  std::vector<float> recoY;
258  std::vector<float> recoZ;
259 
260  int n_points = trk->NumberTrajectoryPoints();
261 
262  for (int i = 0; i < n_points; i++) {
263  auto const& pos = trk->LocationAtPoint(i);
264  recoX.push_back(pos.X());
265  recoY.push_back(pos.Y());
266  recoZ.push_back(pos.Z());
267  }
268 
269  if (recoX.size() < 2)
270  return -1.0;
271 
272  if (!plotRecoTracks_(recoX, recoY, recoZ))
273  return -1.0;
274 
275  constexpr double seg_size{10.};
276 
277  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
278  if (!segments.has_value())
279  return -1.0;
280 
281  auto const seg_steps = segments->x.size();
282  if (seg_steps < 2)
283  return -1;
284 
285  double const recoL = segments->L.at(seg_steps - 1);
286  if (recoL < minLength || recoL > maxLength)
287  return -1;
288 
289  std::vector<float> dEi;
290  std::vector<float> dEj;
291  std::vector<float> dthij;
292  std::vector<float> ind;
293  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0)
294  return -1.0;
295 
296  double logL = 1e+16;
297  double bf = -666.0; // double errs = -666.0;
298 
299  int const start1{};
300  int const end1{750};
301  int const start2{};
302  int const end2{}; // 800.0;
303 
304  for (int k = start1; k <= end1; ++k) {
305  double const p_test = 0.001 + k * 0.01;
306 
307  for (int l = start2; l <= end2; ++l) {
308  double const res_test = 2.0; // 0.001+l*1.0;
309  double const fv = my_mcs_llhd(dEi, dEj, dthij, ind, p_test, res_test);
310 
311  if (fv < logL) {
312  bf = p_test;
313  logL = fv;
314  // errs = res_test;
315  }
316  }
317  }
318  return bf;
319  }
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:390
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:646
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:1243
e
do i e
Definition: updateVersionList.sh:68
k
pdgs k
Definition: selectors.fcl:22
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:614
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD ( art::Ptr< recob::Track > const &  trk)
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD ( art::Ptr< recob::Track > const &  trk)
TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const &  trk)

Definition at line 322 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

324  {
325  double const LLHDp = GetMuMultiScatterLLHD3(trk, true);
326  double const LLHDm = GetMuMultiScatterLLHD3(trk, false);
327 
328  if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
329  int const n_points = trk->NumberTrajectoryPoints();
330  return trk->LocationAtPoint<TVector3>(n_points - 1);
331  } else {
332  return trk->LocationAtPoint<TVector3>(0);
333  }
334 
335  // Should never get here
336  return TVector3{};
337  }
trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:340
TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const &  trk)
TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const &  trk)
double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3 ( art::Ptr< recob::Track > const &  trk,
bool  dir 
)

Definition at line 340 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

343  {
344  std::vector<float> recoX;
345  std::vector<float> recoY;
346  std::vector<float> recoZ;
347 
348  int const n_points = trk->NumberTrajectoryPoints();
349  for (int i = 0; i < n_points; ++i) {
350  auto const index = dir ? i : n_points - 1 - i;
351  auto const& pos = trk->LocationAtPoint(index);
352  recoX.push_back(pos.X());
353  recoY.push_back(pos.Y());
354  recoZ.push_back(pos.Z());
355  }
356 
357  if (recoX.size() < 2)
358  return -1.0;
359 
360  if (!plotRecoTracks_(recoX, recoY, recoZ))
361  return -1.0;
362 
363  constexpr double seg_size{5.0};
364  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
365  if (!segments.has_value())
366  return -1.0;
367 
368  auto const seg_steps = segments->x.size();
369  if (seg_steps < 2)
370  return -1;
371 
372  double const recoL = segments->L.at(seg_steps - 1);
373  if (recoL < 15.0 || recoL > maxLength)
374  return -1;
375 
376  std::vector<float> dEi;
377  std::vector<float> dEj;
378  std::vector<float> dthij;
379  std::vector<float> ind;
380  if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size) != 0)
381  return -1.0;
382 
383  double const p_range = recoL * kcal;
384  double const logL = my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
385 
386  return logL;
387  }
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:390
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:646
dropbox.dir
tuple dir
Definition: dropbox.py:28
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: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:1243
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:614
double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3 ( art::Ptr< recob::Track > const &  trk,
bool  dir 
)
double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3 ( art::Ptr< recob::Track > const &  trk,
bool  dir 
)
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

Definition at line 646 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

650  {
651  double stag = 0.0;
652 
653  int a1 = xxx.size();
654  int a2 = yyy.size();
655  int a3 = zzz.size();
656 
657  if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
658  cout << " ( Digitize reco tacks ) Error ! " << endl;
659  return std::nullopt;
660  }
661 
662  int const stopper = seg_stop / seg_size;
663 
664  std::vector<float> segx, segnx;
665  std::vector<float> segy, segny;
666  std::vector<float> segz, segnz;
667  std::vector<float> segL;
668 
669  int ntot = 0;
670 
671  n_seg = 0;
672 
673  double x0{};
674  double y0{};
675  double z0{};
676 
677  double x00 = xxx.at(0);
678  double y00 = yyy.at(0);
679  double z00 = zzz.at(0);
680 
681  int indC = 0;
682 
683  std::vector<float> vx;
684  std::vector<float> vy;
685  std::vector<float> vz;
686 
687  for (int i = 0; i < a1; i++) {
688  x0 = xxx.at(i);
689  y0 = yyy.at(i);
690  z0 = zzz.at(i);
691 
692  double const RR0 =
693  std::sqrt(cet::sum_of_squares(x00 - x0, y00 - y0, z00 - z0));
694 
695  if (RR0 >= stag) {
696  segx.push_back(x0);
697  segy.push_back(y0);
698  segz.push_back(z0);
699 
700  segL.push_back(stag);
701 
702  x_seg[n_seg] = x0;
703  y_seg[n_seg] = y0;
704  z_seg[n_seg] = z0;
705 
706  n_seg++;
707 
708  vx.push_back(x0);
709  vy.push_back(y0);
710  vz.push_back(z0);
711 
712  ntot++;
713 
714  indC = i + 1;
715 
716  break;
717  }
718  }
719 
720  for (int i = indC; i < a1 - 1; i++) {
721  double const x1 = xxx.at(i);
722  double const y1 = yyy.at(i);
723  double const z1 = zzz.at(i);
724 
725  double const dr1 = std::sqrt(cet::sum_of_squares(x1 - x0, y1 - y0, z1 - z0));
726 
727  double const x2 = xxx.at(i + 1);
728  double const y2 = yyy.at(i + 1);
729  double const z2 = zzz.at(i + 1);
730 
731  double const dr2 = std::sqrt(cet::sum_of_squares(x2 - x0, y2 - y0, z2 - z0));
732 
733  if (dr1 < seg_size) {
734  vx.push_back(x1);
735  vy.push_back(y1);
736  vz.push_back(z1);
737 
738  ntot++;
739  }
740 
741  if (dr1 <= seg_size && dr2 > seg_size) {
742  double const dx = x2 - x1;
743  double const dy = y2 - y1;
744  double const dz = z2 - z1;
745  double const dr = std::sqrt(dx * dx + dy * dy + dz * dz);
746 
747  if (dr == 0) {
748  cout << " ( Zero ) Error ! " << endl;
749  return std::nullopt;
750  }
751 
752  double const beta =
753  2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
754 
755  double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
756  double const delta = beta * beta - 4.0 * gamma;
757 
758  if (delta < 0.0) {
759  cout << " ( Discriminant ) Error ! " << endl;
760  return std::nullopt;
761  }
762 
763  double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
764  double const t = lysi1;
765 
766  double const xp = x1 + t * dx;
767  double const yp = y1 + t * dy;
768  double const zp = z1 + t * dz;
769 
770  segx.push_back(xp);
771  segy.push_back(yp);
772  segz.push_back(zp);
773 
774  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
775 
776  x_seg[n_seg] = xp;
777  y_seg[n_seg] = yp;
778  z_seg[n_seg] = zp;
779  n_seg++;
780 
781  x0 = xp;
782  y0 = yp;
783  z0 = zp;
784 
785  vx.push_back(x0);
786  vy.push_back(y0);
787  vz.push_back(z0);
788 
789  ntot++;
790 
791  auto const na = vx.size();
792  double sumx = 0.0;
793  double sumy = 0.0;
794  double sumz = 0.0;
795 
796  for (std::size_t i = 0; i < na; ++i) {
797  sumx += vx.at(i);
798  sumy += vy.at(i);
799  sumz += vz.at(i);
800  }
801 
802  sumx /= na;
803  sumy /= na;
804  sumz /= na;
805 
806  std::vector<double> mx;
807  std::vector<double> my;
808  std::vector<double> mz;
809 
810  TMatrixDSym m(3);
811 
812  for (std::size_t i = 0; i < na; ++i) {
813  double const xxw1 = vx.at(i);
814  double const yyw1 = vy.at(i);
815  double const zzw1 = vz.at(i);
816 
817  mx.push_back(xxw1 - sumx);
818  my.push_back(yyw1 - sumy);
819  mz.push_back(zzw1 - sumz);
820 
821  double const xxw0 = mx.at(i);
822  double const yyw0 = my.at(i);
823  double const zzw0 = mz.at(i);
824 
825  m(0, 0) += xxw0 * xxw0 / na;
826  m(0, 1) += xxw0 * yyw0 / na;
827  m(0, 2) += xxw0 * zzw0 / na;
828 
829  m(1, 0) += yyw0 * xxw0 / na;
830  m(1, 1) += yyw0 * yyw0 / na;
831  m(1, 2) += yyw0 * zzw0 / na;
832 
833  m(2, 0) += zzw0 * xxw0 / na;
834  m(2, 1) += zzw0 * yyw0 / na;
835  m(2, 2) += zzw0 * zzw0 / na;
836  }
837 
838  TMatrixDSymEigen me(m);
839 
840  TVectorD eigenval = me.GetEigenValues();
841  TMatrixD eigenvec = me.GetEigenVectors();
842 
843  double max1 = -666.0;
844 
845  double ind1 = 0;
846 
847  for (int i = 0; i < 3; ++i) {
848  double const p1 = eigenval(i);
849 
850  if (p1 > max1) {
851  max1 = p1;
852  ind1 = i;
853  }
854  }
855 
856  double ax = eigenvec(0, ind1);
857  double ay = eigenvec(1, ind1);
858  double az = eigenvec(2, ind1);
859 
860  if (n_seg > 1) {
861  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
862  ax = std::abs(ax);
863  else
864  ax = -1.0 * std::abs(ax);
865 
866  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
867  ay = std::abs(ay);
868  else
869  ay = -1.0 * std::abs(ay);
870 
871  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
872  az = std::abs(az);
873  else
874  az = -1.0 * std::abs(az);
875 
876  segnx.push_back(ax);
877  segny.push_back(ay);
878  segnz.push_back(az);
879  }
880  else
881  return std::nullopt;
882 
883  ntot = 0;
884 
885  vx.clear();
886  vy.clear();
887  vz.clear();
888 
889  vx.push_back(x0);
890  vy.push_back(y0);
891  vz.push_back(z0);
892 
893  ntot++;
894  }
895  else if (dr1 > seg_size) {
896  double const dx = x1 - x0;
897  double const dy = y1 - y0;
898  double const dz = z1 - z0;
899  double const dr = std::sqrt(cet::sum_of_squares(dx, dy, dz));
900 
901  if (dr == 0) {
902  cout << " ( Zero ) Error ! " << endl;
903  return std::nullopt;
904  }
905 
906  double const t = seg_size / dr;
907  double const xp = x0 + t * dx;
908  double const yp = y0 + t * dy;
909  double const zp = z0 + t * dz;
910 
911  segx.push_back(xp);
912  segy.push_back(yp);
913  segz.push_back(zp);
914  segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
915 
916  x_seg[n_seg] = xp;
917  y_seg[n_seg] = yp;
918  z_seg[n_seg] = zp;
919  n_seg++;
920 
921  x0 = xp;
922  y0 = yp;
923  z0 = zp;
924 
925  i = (i - 1);
926 
927  vx.push_back(x0);
928  vy.push_back(y0);
929  vz.push_back(z0);
930 
931  ntot++;
932 
933  double na = vx.size();
934 
935  double sumx = 0.0;
936  double sumy = 0.0;
937  double sumz = 0.0;
938 
939  for (std::size_t i = 0; i < na; ++i) {
940  sumx += vx.at(i);
941  sumy += vy.at(i);
942  sumz += vz.at(i);
943  }
944 
945  sumx /= na;
946  sumy /= na;
947  sumz /= na;
948 
949  std::vector<double> mx;
950  std::vector<double> my;
951  std::vector<double> mz;
952 
953  TMatrixDSym m(3);
954 
955  for (int i = 0; i < na; ++i) {
956  double const xxw1 = vx.at(i);
957  double const yyw1 = vy.at(i);
958  double const zzw1 = vz.at(i);
959 
960  mx.push_back(xxw1 - sumx);
961  my.push_back(yyw1 - sumy);
962  mz.push_back(zzw1 - sumz);
963 
964  double const xxw0 = mx.at(i);
965  double const yyw0 = my.at(i);
966  double const zzw0 = mz.at(i);
967 
968  m(0, 0) += xxw0 * xxw0 / na;
969  m(0, 1) += xxw0 * yyw0 / na;
970  m(0, 2) += xxw0 * zzw0 / na;
971 
972  m(1, 0) += yyw0 * xxw0 / na;
973  m(1, 1) += yyw0 * yyw0 / na;
974  m(1, 2) += yyw0 * zzw0 / na;
975 
976  m(2, 0) += zzw0 * xxw0 / na;
977  m(2, 1) += zzw0 * yyw0 / na;
978  m(2, 2) += zzw0 * zzw0 / na;
979  }
980 
981  TMatrixDSymEigen me(m);
982 
983  TVectorD eigenval = me.GetEigenValues();
984  TMatrixD eigenvec = me.GetEigenVectors();
985 
986  double max1 = -666.0;
987  double ind1 = 0;
988 
989  for (int i = 0; i < 3; ++i) {
990  double const p1 = eigenval(i);
991 
992  if (p1 > max1) {
993  max1 = p1;
994  ind1 = i;
995  }
996  }
997 
998  double ax = eigenvec(0, ind1);
999  double ay = eigenvec(1, ind1);
1000  double az = eigenvec(2, ind1);
1001 
1002  if (n_seg > 1) {
1003  if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
1004  ax = std::abs(ax);
1005  else
1006  ax = -1.0 * std::abs(ax);
1007 
1008  if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
1009  ay = std::abs(ay);
1010  else
1011  ay = -1.0 * std::abs(ay);
1012 
1013  if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
1014  az = std::abs(az);
1015  else
1016  az = -1.0 * std::abs(az);
1017 
1018  segnx.push_back(ax);
1019  segny.push_back(ay);
1020  segnz.push_back(az);
1021  }
1022 
1023  else
1024  return std::nullopt;
1025 
1026  ntot = 0;
1027 
1028  vx.clear();
1029  vy.clear();
1030  vz.clear();
1031 
1032  vx.push_back(x0);
1033  vy.push_back(y0);
1034  vz.push_back(z0);
1035 
1036  ntot++;
1037  }
1038 
1039  if (n_seg >= (stopper + 1.0) && seg_stop != -1)
1040  break;
1041  }
1042 
1043  delete gr_seg_xyz;
1044  gr_seg_xyz = new TPolyLine3D{n_seg, z_seg, x_seg, y_seg};
1045  gr_seg_yz = TGraph{n_seg, z_seg, y_seg};
1046  gr_seg_xz = TGraph{n_seg, z_seg, x_seg};
1047  gr_seg_xy = TGraph{n_seg, x_seg, y_seg};
1048 
1049  return std::make_optional<Segments>(Segments{segx, segnx, segy, segny, segz, segnz, segL});
1050  }
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
my
my($xml, $fcl, $workdir, $check, $merge)
a2
#define a2
Definition: HoughBaseAlg.cxx:64
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
a3
#define a3
Definition: HoughBaseAlg.cxx:65
p1
physics associatedGroupsWithLeft p1
Definition: test_associatedgroupswithleft_full.fcl:24
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
a1
#define a1
Definition: HoughBaseAlg.cxx:63
std::optional<Segments> trkf::TrackMomentumCalculator::getSegTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz,
double  seg_size 
)
private
std::optional<Segments> trkf::TrackMomentumCalculator::getSegTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz,
double  seg_size 
)
private
double trkf::TrackMomentumCalculator::GetTrackMomentum ( double  trkrange,
int  pdg 
) const

Definition at line 144 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

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

Definition at line 1224 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

1225  {
1226  if (s == 0.) {
1227  cout << " Error : The code tries to divide by zero ! " << endl;
1228  return 0.;
1229  }
1230 
1231  double const arg = (xx - Q) / s;
1232  double const result =
1233  -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1234 
1235  if (std::isnan(result) || std::isinf(result)) {
1236  cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
1237  }
1238 
1239  return result;
1240  }
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
double trkf::TrackMomentumCalculator::my_g ( double  xx,
double  Q,
double  s 
) const
private
double trkf::TrackMomentumCalculator::my_g ( double  xx,
double  Q,
double  s 
) const
private
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

Definition at line 1243 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

1249  {
1250  double p = x0;
1251  double theta0x = x1;
1252 
1253  double result = 0.0;
1254  double nnn1 = dEi.size();
1255 
1256  double red_length = (10.0) / 14.0;
1257  double addth = 0;
1258 
1259  for (int i = 0; i < nnn1; i++) {
1260  double Ei = p - dEi.at(i);
1261  double Ej = p - dEj.at(i);
1262 
1263  if (Ei > 0 && Ej < 0)
1264  addth = 3.14 * 1000.0;
1265 
1266  Ei = std::abs(Ei);
1267  Ej = std::abs(Ej);
1268 
1269  double tH0 = (13.6 / std::sqrt(Ei * Ej)) *
1270  (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1271 
1272  double rms = -1.0;
1273 
1274  if (ind.at(i) == 1) {
1275  rms = std::sqrt(tH0 * tH0 + pow(theta0x, 2.0));
1276 
1277  double const DT = dthij.at(i) + addth;
1278  double const prob = my_g(DT, 0.0, rms);
1279 
1280  result = result - 2.0 * prob;
1281  }
1282  }
1283 
1284  if (std::isnan(result) || std::isinf(result)) {
1285  std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
1286  }
1287  return result;
1288  }
trkf::TrackMomentumCalculator::my_g
double my_g(double xx, double Q, double s) const
Definition: larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx:1224
p
pdgs p
Definition: selectors.fcl:22
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
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
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
bool trkf::TrackMomentumCalculator::plotRecoTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz 
)
private

Definition at line 614 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.cxx.

617  {
618  auto const n = xxx.size();
619  auto const y_size = yyy.size();
620  auto const z_size = zzz.size();
621 
622  if (n != y_size || n != z_size) {
623  cout << " ( Get reco tacks ) Error ! " << endl;
624  return false;
625  }
626 
627  // Here, we perform a const-cast to float* because, sadly,
628  // TPolyLine3D requires a pointer to a non-const object. We will
629  // trust that ROOT does not mess around with the underlying data.
630  auto xs = const_cast<float*>(xxx.data());
631  auto ys = const_cast<float*>(yyy.data());
632  auto zs = const_cast<float*>(zzz.data());
633 
634  auto const narrowed_size =
635  static_cast<int>(n); // ROOT requires a signed integral type
636  delete gr_reco_xyz;
637  gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
638  gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
639  gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
640  gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
641 
642  return true;
643  }
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
bool trkf::TrackMomentumCalculator::plotRecoTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz 
)
private
bool trkf::TrackMomentumCalculator::plotRecoTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz 
)
private

Member Data Documentation

TGraph trkf::TrackMomentumCalculator::gr_reco_xy {}
private

Definition at line 94 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

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

Definition at line 93 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_reco_xz {}
private

Definition at line 96 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_reco_yz {}
private

Definition at line 95 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_xy {}
private

Definition at line 99 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

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

Definition at line 98 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_xz {}
private

Definition at line 101 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_yz {}
private

Definition at line 100 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

double trkf::TrackMomentumCalculator::maxLength
private

Definition at line 81 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

double trkf::TrackMomentumCalculator::minLength
private

Definition at line 80 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

int trkf::TrackMomentumCalculator::n_seg {}
private

Definition at line 68 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

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

Definition at line 77 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

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

Definition at line 67 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

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

Definition at line 78 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

float trkf::TrackMomentumCalculator::steps_size {10.}
private

Definition at line 76 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

float trkf::TrackMomentumCalculator::x_seg
private

Definition at line 70 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

float trkf::TrackMomentumCalculator::y_seg
private

Definition at line 71 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

float trkf::TrackMomentumCalculator::z_seg
private

Definition at line 72 of file larreco/larreco/RecoAlg/TrackMomentumCalculator.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5