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)
 
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
 

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 22 of file TrackMomentumCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 133 of file 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  }
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55

Member Function Documentation

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

Definition at line 1180 of file TrackMomentumCalculator.cxx.

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

Definition at line 1052 of file TrackMomentumCalculator.cxx.

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

Definition at line 493 of file 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  auto const& pos = trk->LocationAtPoint(i);
504  recoX.push_back(pos.X());
505  recoY.push_back(pos.Y());
506  recoZ.push_back(pos.Z());
507  }
508 
509  if (recoX.size() < 2)
510  return -1.0;
511 
512  if (!plotRecoTracks_(recoX, recoY, recoZ))
513  return -1.0;
514 
515  double const seg_size{steps_size};
516  auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
517  if (!segments.has_value())
518  return -1.0;
519 
520  auto const seg_steps = segments->x.size();
521  if (seg_steps < 2)
522  return -1;
523 
524  double const recoL = segments->L.at(seg_steps - 1);
525  if (recoL < minLength || recoL > maxLength)
526  return -1;
527 
528  double ymax = -999.0;
529  double ymin = +999.0;
530 
531  std::vector<double> xmeas;
532  std::vector<double> ymeas;
533  std::vector<double> eymeas;
534  xmeas.reserve(n_steps);
535  ymeas.reserve(n_steps);
536  eymeas.reserve(n_steps);
537  for (int j = 0; j < n_steps; j++) {
538  double const trial = steps.at(j);
539  auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
540 
541  if (std::isnan(mean) || std::isinf(mean)) {
542  mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
543  continue;
544  }
545  if (std::isnan(rms) || std::isinf(rms)) {
546  mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
547  continue;
548  }
549  if (std::isnan(rmse) || std::isinf(rmse)) {
550  mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
551  continue;
552  }
553 
554  xmeas.push_back(trial); // Is this what is intended?
555  ymeas.push_back(rms);
556  eymeas.push_back(std::sqrt(cet::sum_of_squares(rmse, 0.05 * rms))); // <--- conservative syst. error to fix chi^{2} behaviour !!!
557 
558  if (ymin > rms)
559  ymin = rms;
560  if (ymax < rms)
561  ymax = rms;
562  }
563 
564  assert(xmeas.size() == ymeas.size());
565  assert(xmeas.size() == eymeas.size());
566  if (xmeas.empty()) {
567  return -1.0;
568  }
569 
570  TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
571 
572  gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
573  "thickness in cm; (#Delta#theta)_{rms} in mrad");
574 
575  gr_meas.SetLineColor(kBlack);
576  gr_meas.SetMarkerColor(kBlack);
577  gr_meas.SetMarkerStyle(20);
578  gr_meas.SetMarkerSize(1.2);
579 
580  gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0),
581  steps.at(n_steps - 1) + steps.at(0));
582  gr_meas.SetMinimum(0.0);
583  gr_meas.SetMaximum(1.80 * ymax);
584 
585  ROOT::Minuit2::Minuit2Minimizer mP{};
586  FcnWrapper const wrapper{move(xmeas), move(ymeas), move(eymeas)};
587  ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
588 
589  mP.SetFunction(FCA);
590  mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, 7.5);
591  mP.SetLimitedVariable(1, "#delta#theta", 0.0, 1.0, 0.0, 45.0);
592  mP.SetMaxFunctionCalls(1.E9);
593  mP.SetMaxIterations(1.E9);
594  mP.SetTolerance(0.01);
595  mP.SetStrategy(2);
596  mP.SetErrorDef(1.0);
597 
598  bool const mstatus = mP.Minimize();
599 
600  mP.Hesse();
601 
602  const double* pars = mP.X();
603  const double* erpars = mP.Errors();
604 
605  double const deltap = (recoL * kcal) / 2.0;
606 
607  double const p_mcs = pars[0] + deltap;
608  double const p_mcs_e [[maybe_unused]] = erpars[0];
609  return mstatus ? p_mcs : -1.0;
610  }
cet::sqlite::rms
double rms(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:40
cet::sum_of_squares
constexpr T sum_of_squares(T x, T y)
Definition: pow.h:139
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
trkf::TrackMomentumCalculator::getSegTracks_
std::optional< Segments > getSegTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
Definition: TrackMomentumCalculator.cxx:645
submit_ppfxjobs.wrapper
wrapper
Definition: submit_ppfxjobs.py:165
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
trkf::TrackMomentumCalculator::getDeltaThetaRMS_
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Definition: TrackMomentumCalculator.cxx:1052
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:208
cet::sqlite::mean
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:613
double trkf::TrackMomentumCalculator::GetMomentumMultiScatterLLHD ( art::Ptr< recob::Track > const &  trk)

Definition at line 253 of file 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: TrackMomentumCalculator.cxx:390
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
l
static QStrList * l
Definition: config.cpp:1044
trkf::TrackMomentumCalculator::getSegTracks_
std::optional< Segments > getSegTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
Definition: TrackMomentumCalculator.cxx:645
e
const double e
Definition: gUpMuFluxGen.cxx:165
trkf::TrackMomentumCalculator::my_mcs_llhd
double my_mcs_llhd(std::vector< float > const &dEi, std::vector< float > const &dEj, std::vector< float > const &dthij, std::vector< float > const &ind, double x0, double x1) const
Definition: TrackMomentumCalculator.cxx:1242
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:613
TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const &  trk)

Definition at line 322 of file 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  }
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Definition: TrackMomentumCalculator.cxx:340
double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3 ( art::Ptr< recob::Track > const &  trk,
bool  dir 
)

Definition at line 340 of file 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: TrackMomentumCalculator.cxx:390
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
dir
string dir
Definition: INukeNucleonCorr.cxx:61
trkf::TrackMomentumCalculator::getSegTracks_
std::optional< Segments > getSegTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz, double seg_size)
Definition: TrackMomentumCalculator.cxx:645
trkf::TrackMomentumCalculator::my_mcs_llhd
double my_mcs_llhd(std::vector< float > const &dEi, std::vector< float > const &dEj, std::vector< float > const &dthij, std::vector< float > const &ind, double x0, double x1) const
Definition: TrackMomentumCalculator.cxx:1242
trkf::TrackMomentumCalculator::plotRecoTracks_
bool plotRecoTracks_(std::vector< float > const &xxx, std::vector< float > const &yyy, std::vector< float > const &zzz)
Definition: TrackMomentumCalculator.cxx:613
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 645 of file TrackMomentumCalculator.cxx.

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

Definition at line 144 of file 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  }
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
make_THn_beam_input.Momentum
Momentum
Definition: make_THn_beam_input.py:70
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
double trkf::TrackMomentumCalculator::my_g ( double  xx,
double  Q,
double  s 
) const
private

Definition at line 1223 of file TrackMomentumCalculator.cxx.

1224  {
1225  if (s == 0.) {
1226  cout << " Error : The code tries to divide by zero ! " << endl;
1227  return 0.;
1228  }
1229 
1230  double const arg = (xx - Q) / s;
1231  double const result =
1232  -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1233 
1234  if (std::isnan(result) || std::isinf(result)) {
1235  cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
1236  }
1237 
1238  return result;
1239  }
result
static QCString result
Definition: fortranscanner.cpp:56614
s
static QCString * s
Definition: config.cpp:1042
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
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 1242 of file TrackMomentumCalculator.cxx.

1248  {
1249  double p = x0;
1250  double theta0x = x1;
1251 
1252  double result = 0.0;
1253  double nnn1 = dEi.size();
1254 
1255  double red_length = (10.0) / 14.0;
1256  double addth = 0;
1257 
1258  for (int i = 0; i < nnn1; i++) {
1259  double Ei = p - dEi.at(i);
1260  double Ej = p - dEj.at(i);
1261 
1262  if (Ei > 0 && Ej < 0)
1263  addth = 3.14 * 1000.0;
1264 
1265  Ei = std::abs(Ei);
1266  Ej = std::abs(Ej);
1267 
1268  double tH0 = (13.6 / std::sqrt(Ei * Ej)) *
1269  (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1270 
1271  double rms = -1.0;
1272 
1273  if (ind.at(i) == 1) {
1274  rms = std::sqrt(tH0 * tH0 + pow(theta0x, 2.0));
1275 
1276  double const DT = dthij.at(i) + addth;
1277  double const prob = my_g(DT, 0.0, rms);
1278 
1279  result = result - 2.0 * prob;
1280  }
1281  }
1282 
1283  if (std::isnan(result) || std::isinf(result)) {
1284  std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
1285  }
1286  return result;
1287  }
result
static QCString result
Definition: fortranscanner.cpp:56614
cet::sqlite::rms
double rms(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:40
trkf::TrackMomentumCalculator::my_g
double my_g(double xx, double Q, double s) const
Definition: TrackMomentumCalculator.cxx:1223
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
test.p
p
Definition: test.py:223
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
bool trkf::TrackMomentumCalculator::plotRecoTracks_ ( std::vector< float > const &  xxx,
std::vector< float > const &  yyy,
std::vector< float > const &  zzz 
)
private

Definition at line 613 of file TrackMomentumCalculator.cxx.

616  {
617  auto const n = xxx.size();
618  auto const y_size = yyy.size();
619  auto const z_size = zzz.size();
620 
621  if (n != y_size || n != z_size) {
622  cout << " ( Get reco tacks ) Error ! " << endl;
623  return false;
624  }
625 
626  // Here, we perform a const-cast to float* because, sadly,
627  // TPolyLine3D requires a pointer to a non-const object. We will
628  // trust that ROOT does not mess around with the underlying data.
629  auto xs = const_cast<float*>(xxx.data());
630  auto ys = const_cast<float*>(yyy.data());
631  auto zs = const_cast<float*>(zzz.data());
632 
633  auto const narrowed_size =
634  static_cast<int>(n); // ROOT requires a signed integral type
635  delete gr_reco_xyz;
636  gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
637  gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
638  gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
639  gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
640 
641  return true;
642  }
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
genie::units::zs
static constexpr double zs
Definition: Units.h:102
genie::units::ys
static constexpr double ys
Definition: Units.h:103
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030

Member Data Documentation

TGraph trkf::TrackMomentumCalculator::gr_reco_xy {}
private

Definition at line 95 of file TrackMomentumCalculator.h.

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

Definition at line 94 of file TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_reco_xz {}
private

Definition at line 97 of file TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_reco_yz {}
private

Definition at line 96 of file TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_xy {}
private

Definition at line 100 of file TrackMomentumCalculator.h.

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

Definition at line 99 of file TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_xz {}
private

Definition at line 102 of file TrackMomentumCalculator.h.

TGraph trkf::TrackMomentumCalculator::gr_seg_yz {}
private

Definition at line 101 of file TrackMomentumCalculator.h.

double trkf::TrackMomentumCalculator::maxLength
private

Definition at line 82 of file TrackMomentumCalculator.h.

double trkf::TrackMomentumCalculator::minLength
private

Definition at line 81 of file TrackMomentumCalculator.h.

int trkf::TrackMomentumCalculator::n_seg {}
private

Definition at line 69 of file TrackMomentumCalculator.h.

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

Definition at line 78 of file TrackMomentumCalculator.h.

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

Definition at line 68 of file TrackMomentumCalculator.h.

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

Definition at line 79 of file TrackMomentumCalculator.h.

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

Definition at line 77 of file TrackMomentumCalculator.h.

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

Definition at line 71 of file TrackMomentumCalculator.h.

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

Definition at line 72 of file TrackMomentumCalculator.h.

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

Definition at line 73 of file TrackMomentumCalculator.h.

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