6 #include "cetlib/pow.h" 17 #include "Math/Functor.h" 18 #include "Math/GenVector/PositionVector3D.h" 19 #include "Minuit2/Minuit2Minimizer.h" 22 #include "TGraphErrors.h" 24 #include "TMatrixDSymEigen.h" 25 #include "TMatrixDSymfwd.h" 26 #include "TMatrixDfwd.h" 28 #include "TMatrixTSym.h" 29 #include "TPolyLine3D.h" 31 #include "TVectorDfwd.h" 42 constexpr
double LAr_density{1.396};
43 constexpr
double rad_length{14.0};
44 constexpr
double m_muon{0.1057};
45 constexpr
auto range_gramper_cm()
47 std::array<double, 73> Range_grampercm{
48 {0.9833, 1.36, 1.786, 2.507, 3.321, 4.859, 6.598, 8.512, 10.58, 12.78,
49 15.1, 17.52, 20.04, 25.31, 30.84, 36.59, 42.5, 54.73, 67.32, 86.66,
50 106.3, 139.4, 172.5, 205.6, 238.5, 271.1, 303.5, 335.7, 367.7, 431.0,
51 493.4, 555.2, 616.3, 736.8, 855.2, 1030.0, 1202.0, 1482.0, 1758.0, 2029.0,
52 2297.0, 2562.0, 2825.0, 3085.0, 3343.0, 3854.0, 4359.0, 4859.0, 5354.0, 6333.0,
53 7298.0, 8726.0, 10130.0, 12430.0, 14690.0, 16920.0, 19100.0, 21260.0, 23380.0, 25480.0,
54 27550.0, 31610.0, 35580.0, 39460.0, 43260.0, 50620.0, 57680.0, 67780.0, 77340.0, 92220.0,
55 1.06e+05, 1.188e+05, 1.307e+05}};
56 for (
double&
value : Range_grampercm) {
59 return Range_grampercm;
62 constexpr
auto Range_grampercm = range_gramper_cm();
63 constexpr std::array<double, 73> KE_MeV{
64 {10.0, 12.0, 14.0, 17.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0,
65 50.0, 55.0, 60.0, 70.0, 80.0, 90.0, 100.0, 120.0, 140.0, 170.0,
66 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 700.0,
67 800.0, 900.0, 1000.0, 1200.0, 1400.0, 1700.0, 2000.0, 2500.0, 3000.0, 3500.0,
68 4000.0, 4500.0, 5000.0, 5500.0, 6000.0, 7000.0, 8000.0, 9000.0, 10000.0, 12000.0,
69 14000.0, 17000.0, 20000.0, 25000.0, 30000.0, 35000.0, 40000.0, 45000.0, 50000.0, 55000.0,
70 60000.0, 70000.0, 80000.0, 90000.0, 1
e+05, 1.2e+05, 1.4e+05, 1.7e+05, 2
e+05, 2.5e+05,
71 3
e+05, 3.5e+05, 4
e+05}};
72 TGraph
const KEvsR{73, Range_grampercm.data(), KE_MeV.data()};
73 TSpline3
const KEvsR_spline3{
"KEvsRS", &KEvsR};
76 std::vector<double> dedx_GeV_per_cm()
79 std::vector<double>
dEdx{
80 {5.687, 4.979, 4.461, 3.902, 3.502, 3.042, 2.731, 2.508, 2.34, 2.21, 2.107, 2.023, 1.954,
81 1.848, 1.771, 1.713, 1.67, 1.609, 1.57, 1.536, 1.519, 1.508, 1.51, 1.517, 1.526, 1.537,
82 1.548, 1.559, 1.57, 1.591, 1.61, 1.628, 1.645, 1.675, 1.7, 1.733, 1.761, 1.799, 1.829,
83 1.855, 1.877, 1.897, 1.914, 1.93, 1.944, 1.969, 1.991, 2.01, 2.028, 2.058, 2.084, 2.119,
84 2.149, 2.193, 2.232, 2.269, 2.304, 2.337, 2.369, 2.4, 2.43, 2.49, 2.548, 2.606, 2.663,
85 2.776, 2.888, 3.055, 3.224, 3.498, 3.774, 4.052, 4.332}};
87 value *= LAr_density * 1
e-3;
93 const std::vector<double> dEdx_GeV_per_cm = dedx_GeV_per_cm();
94 const std::vector<double> E_GeV{
95 {0.115695, 0.117697, 0.119693, 0.122694, 0.125695, 0.13069, 0.135694, 0.14069, 0.145714,
96 0.150689, 0.155682, 0.160666, 0.165693, 0.17566, 0.185714, 0.1957, 0.205644, 0.225683,
97 0.245698, 0.275669, 0.305658, 0.355669, 0.405711, 0.45563, 0.505671, 0.555646, 0.605694,
98 0.655676, 0.705661, 0.805664, 0.905689, 1.00557, 1.10606, 1.30529, 1.50571, 1.8061,
99 2.10565, 2.60614, 3.1058, 3.60555, 4.10536, 4.60521, 5.10609, 5.606, 6.10591,
100 7.10579, 8.10569, 9.10561, 10.1106, 12.1105, 14.1104, 17.1103, 20.1103, 25.1102,
101 30.1102, 35.1102, 40.1101, 45.1101, 50.1101, 55.1101, 60.1101, 70.1101, 80.1101,
102 90.1101, 100.1, 120.1, 140.1, 170.1, 200.1, 250.1, 300.1, 350.1,
104 TGraph
const dEdx_vs_E{ndedx, E_GeV.data(), dEdx_GeV_per_cm.data()};
105 TSpline3
const dEdx_vs_E_spline3{
"dEdx_vs_E", &dEdx_vs_E};
107 TVector3
const basex{1, 0, 0};
108 TVector3
const basey{0, 1, 0};
109 TVector3
const basez{0, 0, 1};
110 constexpr
double kcal{0.0024};
112 constexpr
double MomentumDependentConstant(
const double p)
117 return (a / (p * p)) + c;
119 double ComputeExpetecteRMS(
const double p,
const double red_length)
121 double beta = std::sqrt(1 - ((m_muon * m_muon) / (p * p + m_muon * m_muon)));
122 return (MomentumDependentConstant(p) / (p * beta)) *
123 (1.0 + 0.038 * TMath::Log(red_length / cet::square(beta))) * std::sqrt(red_length);
127 explicit FcnWrapper(std::vector<double>&& xmeas,
128 std::vector<double>&& ymeas,
129 std::vector<double>&& eymeas,
131 : xmeas_{std::move(xmeas)}
132 , ymeas_{std::move(ymeas)}
133 , eymeas_{std::move(eymeas)}
134 , correction_{correction}
137 double my_mcs_chi2(
double const*
x)
const 142 double theta0 = x[1];
144 auto const n = xmeas_.size();
145 assert(
n == ymeas_.size());
146 assert(
n == eymeas_.size());
148 for (std::size_t i = 0; i <
n; ++i) {
149 double const xx = xmeas_[i];
150 double const yy = ymeas_[i];
151 double const ey = eymeas_[i];
153 if (
std::abs(ey) < std::numeric_limits<double>::epsilon()) {
154 std::cout <<
" Zero denominator in my_mcs_chi2 ! " << std::endl;
158 double const l0 = xx / rad_length;
162 double beta = std::sqrt(1 - ((m_muon * m_muon) / (p * p + m_muon * m_muon)));
164 res1 = (13.6 / (p * beta)) * (1.0 + 0.038 * TMath::Log(l0 / cet::square(beta))) *
168 res1 = std::hypot(res1, theta0);
170 double const diff = yy - res1;
171 result += cet::square(diff / ey);
174 if (std::isnan(result) || std::isinf(result)) {
175 MF_LOG_DEBUG(
"TrackMomentumCalculator") <<
" Is nan in my_mcs_chi2 ! ";
183 std::vector<double>
const xmeas_;
184 std::vector<double>
const ymeas_;
185 std::vector<double>
const eymeas_;
186 double const correction_;
189 class FcnWrapperLLHD {
191 explicit FcnWrapperLLHD(std::vector<double>&& dEi,
192 std::vector<double>&& dEj,
193 std::vector<double>&& dthij,
194 std::vector<double>&& ind,
195 std::vector<bool>&& dthij_valid,
198 : dEi_{std::move(dEi)}
199 , dEj_{std::move(dEj)}
200 , dthij_{std::move(dthij)}
201 , ind_{std::move(ind)}
202 , dthij_valid_{std::move(dthij_valid)}
203 , stepsize_{stepsize}
204 , correction_{correction}
207 double my_mcs_llhd(
double const* x)
const 211 double red_length = stepsize_ / rad_length;
214 double theta0 = x[1];
217 double Etot = std::hypot(p, m_muon);
222 auto const n = dEi_.size();
224 bool addpenality =
false;
225 for (std::size_t i = 0; i <
n; ++i) {
227 if (Ei >= E_GeV[0]) {
228 dEi = dEdx_vs_E_spline3.Eval(Ei) * stepsize_;
231 dEi = dEdx_GeV_per_cm[0] * stepsize_;
240 if (dthij_valid_.at(i) ==
false)
continue;
243 double pij = std::sqrt(Ei * Ei - m_muon * m_muon);
248 double tH0 = ComputeExpetecteRMS(pij, red_length);
254 double rms_square = -1.0;
259 rms_square = cet::sum_of_squares(tH0, theta0);
264 prob = -0.5 * std::log(2.0 * TMath::Pi()) - 0.5 * std::log(rms_square) -
265 0.5 * DT * DT / rms_square;
267 if (addpenality) { prob -= 2 * rms_square; }
269 result = result - 2.0 * prob;
276 std::vector<double>
const dEi_;
277 std::vector<double>
const dEj_;
278 std::vector<double>
const dthij_;
279 std::vector<double>
const ind_;
280 std::vector<bool>
const dthij_valid_;
281 double const stepsize_;
282 double const correction_;
289 TrackMomentumCalculator::TrackMomentumCalculator(
double const min,
291 double const stepsize,
292 int const angleMethod,
300 for (
int i = 1; i <=
n_steps; i++) {
368 if (trkrange < 0 || std::isnan(trkrange)) {
369 mf::LogError(
"TrackMomentumCalculator") <<
"Invalid track range " << trkrange <<
" return -1";
373 double KE, Momentum, M;
374 constexpr
double Muon_M = m_muon * 1e3, Proton_M = 938.272;
376 if (
abs(pdg) == 13) {
378 KE = KEvsR_spline3.Eval(trkrange);
380 else if (
abs(pdg) == 2212) {
382 if (trkrange > 0 && trkrange <= 80)
383 KE = 29.9317 * std::pow(trkrange, 0.586304);
384 else if (trkrange > 80 && trkrange <= 3.022E3)
385 KE = 149.904 + (3.34146 * trkrange) + (-0.00318856 * trkrange * trkrange) +
386 (4.34587E-6 * trkrange * trkrange * trkrange) +
387 (-3.18146
E-9 * trkrange * trkrange * trkrange * trkrange) +
388 (1.17854E-12 * trkrange * trkrange * trkrange * trkrange * trkrange) +
389 (-1.71763
E-16 * trkrange * trkrange * trkrange * trkrange * trkrange * trkrange);
399 Momentum = std::sqrt((KE * KE) + (2 * M * KE));
401 Momentum = Momentum / 1000;
416 const bool checkValidPoints,
417 const int maxMomentum_MeV,
418 const double min_resolution,
419 const double max_resolution,
420 const bool check_valid_scattered,
421 const double angle_correction)
424 std::vector<double> recoX;
425 std::vector<double> recoY;
426 std::vector<double> recoZ;
430 for (
int i = 0; i < n_points; i++) {
433 recoX.push_back(pos.X());
434 recoY.push_back(pos.Y());
435 recoZ.push_back(pos.Z());
438 if (recoX.size() < 2)
return -1.0;
444 auto const segments =
getSegTracks_(recoX, recoY, recoZ, seg_size, check_valid_scattered);
445 if (!segments.has_value())
return -1.0;
447 auto const seg_steps = segments->x.size();
448 if (seg_steps < 2)
return -1;
450 double const recoSegmentLength = segments->L.at(seg_steps - 1);
451 if (recoSegmentLength < minLength || recoSegmentLength >
maxLength)
return -1;
453 std::vector<double> dEi;
454 std::vector<double> dEj;
455 std::vector<double> dthij;
456 std::vector<double> ind;
457 std::vector<bool> dthij_valid = segments->nvalid;
458 if (
getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
461 auto const ndEi = dEi.size();
462 if (ndEi < 1)
return -1;
464 double correction = 1.;
468 auto recoL = trk->
Length();
471 ROOT::Minuit2::Minuit2Minimizer mP{};
472 FcnWrapperLLHD
const wrapper{std::move(dEi),
476 std::move(dthij_valid),
479 ROOT::Math::Functor FCA([&wrapper](
double const* xs) {
return wrapper.my_mcs_llhd(xs); }, 2);
485 double startpoint = 2;
486 if (startpoint < min_resolution)
487 startpoint = (max_resolution - min_resolution) / 2.;
488 if (max_resolution == 0) startpoint = min_resolution;
493 mP.SetLimitedVariable(0,
"p_{MCS}", minP * 2, minP * 0.1, minP * 0.6, maxMomentum_MeV / 1.e3);
494 mP.SetLimitedVariable(
495 1,
"#delta#theta", startpoint, startpoint / 2., min_resolution, max_resolution);
496 if (max_resolution == 0) { mP.FixVariable(1); }
497 mP.SetMaxFunctionCalls(1.E9);
498 mP.SetMaxIterations(1.E9);
499 mP.SetTolerance(0.01);
503 bool const mstatus = mP.Minimize();
507 const double* pars = mP.X();
509 double const p_mcs = pars[0];
511 return mstatus ? p_mcs : -1.0;
519 if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
534 std::vector<double> recoX;
535 std::vector<double> recoY;
536 std::vector<double> recoZ;
539 for (
int i = 0; i < n_points; ++i) {
540 auto const index = dir ? i : n_points - 1 - i;
542 recoX.push_back(pos.X());
543 recoY.push_back(pos.Y());
544 recoZ.push_back(pos.Z());
547 if (recoX.size() < 2)
return -1.0;
551 constexpr
double seg_size{5.0};
552 auto const segments =
getSegTracks_(recoX, recoY, recoZ, seg_size);
553 if (!segments.has_value())
return -1.0;
555 auto const seg_steps = segments->x.size();
556 if (seg_steps < 2)
return -1;
558 double const recoSegmentLength = segments->L.at(seg_steps - 1);
559 if (recoSegmentLength < 15.0 || recoSegmentLength >
maxLength)
return -1;
561 std::vector<double> dEi;
562 std::vector<double> dEj;
563 std::vector<double> dthij;
564 std::vector<double> ind;
565 double angle_correction = 1;
566 if (
getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
569 auto const recoL = trk->
Length();
570 double const p_range = recoL * kcal;
571 double const logL =
my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
577 std::vector<double>& ej,
578 std::vector<double>& th,
579 std::vector<double>& ind,
582 double const angle_correction)
const 584 auto const& segnx = segments.
nx;
585 auto const& segny = segments.
ny;
586 auto const& segnz = segments.
nz;
587 auto const& segL = segments.
L;
589 int const a1 = segnx.size();
590 int const a2 = segny.size();
591 int const a3 = segnz.size();
593 if (a1 != a2 || a1 != a3) {
594 std::cout <<
" ( Get thij ) Error ! " << std::endl;
600 for (
int i = 0; i < tot - 1; i++) {
601 double const dx = segnx.at(i);
602 double const dy = segny.at(i);
603 double const dz = segnz.at(i);
606 TVector3
const vec_z{dx, dy, dz};
610 double const switcher = basex.Dot(vec_z);
612 vec_y = vec_z.Cross(basex).Unit();
613 vec_x = vec_y.Cross(vec_z);
617 vec_y = basez.Cross(vec_z).Unit();
618 vec_x = vec_y.Cross(vec_z);
621 TVector3
const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
622 TVector3
const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
623 TVector3
const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
625 double const here_dx = segnx.at(i + 1);
626 double const here_dy = segny.at(i + 1);
627 double const here_dz = segnz.at(i + 1);
629 TVector3
const here_vec{here_dx, here_dy, here_dz};
630 TVector3
const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
632 double const scx = rot_here.X();
633 double const scy = rot_here.Y();
634 double const scz = rot_here.Z();
639 constexpr
double ULim = 10000.0;
640 constexpr
double LLim = -10000.0;
642 double const Li = segL.at(i);
643 double const Lj = segL.at(i + 1);
645 if (azy <= ULim && azy >= LLim) {
647 if (azx <= ULim && azx >= LLim) {
657 th.push_back(std::hypot(azx, azy) /
668 const bool checkValidPoints,
669 const int maxMomentum_MeV,
670 const double min_resolution,
671 const double max_resolution)
673 std::vector<double> recoX;
674 std::vector<double> recoY;
675 std::vector<double> recoZ;
679 for (
int i = 0; i < n_points; i++) {
682 recoX.push_back(pos.X());
683 recoY.push_back(pos.Y());
684 recoZ.push_back(pos.Z());
687 if (recoX.size() < 2)
return -1.0;
692 auto const segments =
getSegTracks_(recoX, recoY, recoZ, seg_size);
693 if (!segments.has_value())
return -1.0;
695 auto const seg_steps = segments->x.size();
696 if (seg_steps < 2)
return -1;
698 double const recoSegmentLength = segments->L.at(seg_steps - 1);
699 if (recoSegmentLength < minLength || recoSegmentLength >
maxLength)
return -1;
701 double ymax = -999.0;
702 double ymin = +999.0;
704 std::vector<double> xmeas;
705 std::vector<double> ymeas;
706 std::vector<double> eymeas;
710 for (
int j = 0; j <
n_steps; j++) {
711 double const trial =
steps.at(j);
715 if (std::isnan(
mean) || std::isinf(
mean)) {
716 mf::LogDebug(
"TrackMomentumCalculator") <<
"Returned mean is either nan or infinity.";
719 if (std::isnan(rms) || std::isinf(rms)) {
720 mf::LogDebug(
"TrackMomentumCalculator") <<
"Returned rms is either nan or infinity.";
723 if (std::isnan(rmse) || std::isinf(rmse)) {
724 mf::LogDebug(
"TrackMomentumCalculator") <<
"Returned rmse is either nan or infinity.";
728 if (
mean == -1 && rms == -1 && rmse == -1)
continue;
729 xmeas.push_back(trial);
730 ymeas.push_back(rms);
732 std::hypot(rmse, 0.05 * rms));
734 if (ymin > rms) ymin = rms;
735 if (ymax < rms) ymax = rms;
738 assert(xmeas.size() == ymeas.size());
739 assert(xmeas.size() == eymeas.size());
740 if (xmeas.empty()) {
return -1.0; }
742 TGraphErrors gr_meas{
n_steps, xmeas.data(), ymeas.data(),
nullptr, eymeas.data()};
744 gr_meas.SetTitle(
"(#Delta#theta)_{rms} versus material thickness; Material " 745 "thickness in cm; (#Delta#theta)_{rms} in mrad");
747 gr_meas.SetLineColor(kBlack);
748 gr_meas.SetMarkerColor(kBlack);
749 gr_meas.SetMarkerStyle(20);
750 gr_meas.SetMarkerSize(1.2);
753 gr_meas.SetMinimum(0.0);
754 gr_meas.SetMaximum(1.80 * ymax);
756 double correction = 1.;
759 ROOT::Minuit2::Minuit2Minimizer mP{};
760 FcnWrapper
const wrapper{std::move(xmeas), std::move(ymeas), std::move(eymeas), correction};
761 ROOT::Math::Functor FCA([&wrapper](
double const* xs) {
return wrapper.my_mcs_chi2(xs); }, 2);
764 double startpoint = 2;
765 if (startpoint < min_resolution) startpoint = (max_resolution - min_resolution) / 2.;
766 if (max_resolution == 0) startpoint = min_resolution;
769 mP.SetLimitedVariable(0,
"p_{MCS}", 1.0, 0.01, 0.001, maxMomentum_MeV / 1.e3);
770 mP.SetLimitedVariable(
771 1,
"#delta#theta", startpoint, startpoint / 2., min_resolution, max_resolution);
772 if (max_resolution == 0) { mP.FixVariable(1); }
773 mP.SetMaxFunctionCalls(1.E9);
774 mP.SetMaxIterations(1.E9);
775 mP.SetTolerance(0.01);
779 bool const mstatus = mP.Minimize();
783 const double* pars = mP.X();
785 auto const recoL = trk->
Length();
786 double const deltap = (recoL * kcal) / 2.0;
788 double const p_mcs = pars[0] + deltap;
790 return mstatus ? p_mcs : -1.0;
794 std::vector<double>
const& yyy,
795 std::vector<double>
const& zzz)
797 auto const n = xxx.size();
798 auto const y_size = yyy.size();
799 auto const z_size = zzz.size();
801 if (n != y_size || n != z_size) {
802 cout <<
" ( Get reco tacks ) Error ! " << endl;
809 auto xs =
const_cast<double*
>(xxx.data());
810 auto ys =
const_cast<double*
>(yyy.data());
811 auto zs =
const_cast<double*
>(zzz.data());
813 auto const narrowed_size =
static_cast<int>(
n);
815 gr_reco_xyz =
new TPolyLine3D{narrowed_size, zs, xs, ys};
816 gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
817 gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
818 gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
827 const std::vector<double> segy,
828 const std::vector<double> segz,
829 std::vector<double>& segnx,
830 std::vector<double>& segny,
831 std::vector<double>& segnz,
832 std::vector<bool>& segn_isvalid,
833 std::vector<double>& vx,
834 std::vector<double>& vy,
835 std::vector<double>& vz,
836 bool check_valid_scattered)
838 auto const na = vx.size();
847 if (check_valid_scattered && na <= 2) { isvalid =
false; }
850 for (std::size_t i = 0; i < na; ++i) {
860 std::vector<double> mx;
861 std::vector<double> my;
862 std::vector<double> mz;
867 for (std::size_t i = 0; i < na; ++i) {
868 double const xxw1 = vx.at(i);
869 double const yyw1 = vy.at(i);
870 double const zzw1 = vz.at(i);
872 mx.push_back(xxw1 - sumx);
873 my.push_back(yyw1 - sumy);
874 mz.push_back(zzw1 - sumz);
876 double const xxw0 = mx.at(i);
877 double const yyw0 = my.at(i);
878 double const zzw0 = mz.at(i);
880 m(0, 0) += xxw0 * xxw0 / na;
881 m(0, 1) += xxw0 * yyw0 / na;
882 m(0, 2) += xxw0 * zzw0 / na;
884 m(1, 0) += yyw0 * xxw0 / na;
885 m(1, 1) += yyw0 * yyw0 / na;
886 m(1, 2) += yyw0 * zzw0 / na;
888 m(2, 0) += zzw0 * xxw0 / na;
889 m(2, 1) += zzw0 * yyw0 / na;
890 m(2, 2) += zzw0 * zzw0 / na;
893 TMatrixDSymEigen me(m);
896 TVectorD eigenval = me.GetEigenValues();
897 TMatrixD eigenvec = me.GetEigenVectors();
899 double max1 = -666.0;
904 for (
int i = 0; i < 3; ++i) {
905 double const p1 = eigenval(i);
914 double ax = eigenvec(0, ind1);
915 double ay = eigenvec(1, ind1);
916 double az = eigenvec(2, ind1);
921 if (segx.at(
n_seg - 1) - segx.at(
n_seg - 2) > 0)
926 if (segy.at(
n_seg - 1) - segy.at(
n_seg - 2) > 0)
931 if (segz.at(
n_seg - 1) - segz.at(
n_seg - 2) > 0)
939 segn_isvalid.push_back(isvalid);
955 std::vector<double>
const& xxx,
956 std::vector<double>
const& yyy,
957 std::vector<double>
const& zzz,
958 double const seg_size,
959 bool const check_valid_scattered)
967 if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
968 cout <<
" ( Digitize reco tacks ) Error ! " << endl;
972 int const stopper =
seg_stop / seg_size;
975 std::vector<double> segx, segnx;
976 std::vector<double> segy, segny;
977 std::vector<double> segz, segnz;
978 std::vector<bool> segn_isvalid;
979 std::vector<double> segL;
987 double x00 = xxx.at(0);
988 double y00 = yyy.at(0);
989 double z00 = zzz.at(0);
996 std::vector<double> vx;
997 std::vector<double> vy;
998 std::vector<double> vz;
1000 for (
int i = 0; i <
a1; i++) {
1005 double const RR0 = std::hypot(x00 - x0, y00 - y0, z00 - z0);
1032 for (
int i = indC; i < a1 - 1; i++) {
1034 double const x1 = xxx.at(i);
1035 double const y1 = yyy.at(i);
1036 double const z1 = zzz.at(i);
1039 double const dr1 = std::hypot(x1 - x0, y1 - y0, z1 - z0);
1042 double const x2 = xxx.at(i + 1);
1043 double const y2 = yyy.at(i + 1);
1044 double const z2 = zzz.at(i + 1);
1047 double const dr2 = std::hypot(x2 - x0, y2 - y0, z2 - z0);
1049 if (dr1 < seg_size) {
1066 if (dr1 <= seg_size && dr2 > seg_size) {
1067 double const dx = x2 -
x1;
1068 double const dy = y2 -
y1;
1069 double const dz = z2 - z1;
1070 double const dr = std::hypot(dx, dy, dz);
1073 cout <<
" ( Zero ) Error ! " << endl;
1074 return std::nullopt;
1077 double const beta = 2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
1079 double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
1080 double const delta = beta * beta - 4.0 * gamma;
1083 cout <<
" ( Discriminant ) Error ! " << endl;
1084 return std::nullopt;
1088 double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
1089 double const t = lysi1;
1092 double const xp = x1 + t * dx;
1093 double const yp = y1 + t * dy;
1094 double const zp = z1 + t * dz;
1101 segL.push_back(
n_seg * seg_size);
1121 return std::nullopt;
1126 segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
1133 else if (dr1 > seg_size) {
1139 double const dx = x1 - x0;
1140 double const dy = y1 - y0;
1141 double const dz = z1 - z0;
1142 double const dr = std::hypot(dx, dy, dz);
1145 cout <<
" ( Zero ) Error ! " << endl;
1146 return std::nullopt;
1150 double const t = seg_size / dr;
1151 double const xp = x0 + t * dx;
1152 double const yp = y0 + t * dy;
1153 double const zp = z0 + t * dz;
1158 segL.push_back(1.0 *
n_seg * 1.0 * seg_size + stag);
1176 return std::nullopt;
1181 segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
1198 return std::make_optional<Segments>(
1199 Segments{segx, segnx, segy, segny, segz, segnz, segL, segn_isvalid});
1207 double const thick)
const 1209 auto const& segnx = segments.
nx;
1210 auto const& segny = segments.
ny;
1211 auto const& segnz = segments.
nz;
1212 auto const& segL = segments.
L;
1214 int const a1 = segnx.size();
1215 int const a2 = segny.size();
1216 int const a3 = segnz.size();
1218 if (a1 != a2 || a1 != a3) {
1219 cout <<
" ( Get RMS ) Error ! " << endl;
1220 return std::make_tuple(0., 0., 0.);
1225 std::vector<double> buf0;
1227 for (
int i = 0; i < tot; i++) {
1228 double const dx = segnx.at(i);
1229 double const dy = segny.at(i);
1230 double const dz = segnz.at(i);
1232 TVector3
const vec_z{dx, dy, dz};
1236 double const switcher = basex.Dot(vec_z);
1238 if (switcher <= 0.995) {
1239 vec_y = vec_z.Cross(basex).Unit();
1240 vec_x = vec_y.Cross(vec_z);
1244 vec_y = basez.Cross(vec_z).Unit();
1245 vec_x = vec_y.Cross(vec_z);
1248 TVector3
const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
1249 TVector3
const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
1250 TVector3
const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
1252 double const refL = segL.at(i);
1254 for (
int j = i; j < tot; j++) {
1255 double const L1 = segL.at(j);
1257 double const dz1 = L1 - refL;
1260 double const here_dx = segnx.at(j);
1261 double const here_dy = segny.at(j);
1262 double const here_dz = segnz.at(j);
1264 TVector3
const here_vec{here_dx, here_dy, here_dz};
1265 TVector3
const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
1267 double const scx = rot_here.X();
1268 double const scy = rot_here.Y();
1269 double const scz = rot_here.Z();
1274 constexpr
double ULim = 10000.0;
1275 constexpr
double LLim = -10000.0;
1277 if (azy <= ULim && azy >= LLim) {
1279 if (azx <= ULim && azx >= LLim) {
1281 buf0.push_back(azx);
1284 buf0.push_back(azy);
1287 buf0.push_back(std::hypot(azx, azy));
1296 int const nmeas = buf0.size();
1303 for (
int i = 0; i < nmeas; i++) {
1310 for (
int i = 0; i < nmeas; i++)
1311 rms += ((buf0.at(i)) * (buf0.at(i)));
1314 rms = std::sqrt(rms);
1315 rmse = rms / std::sqrt(2.0 * tot);
1322 double const lev1 = 2.50;
1324 for (
int i = 0; i < nmeas; i++) {
1325 double const amp = buf0.at(i);
1326 if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
1332 rms = rms / (ntot1);
1333 rms = std::sqrt(rms);
1334 rmse = rms / std::sqrt(2.0 * ntot1);
1335 return std::make_tuple(mean, rms, rmse);
1340 double thetayz = -999.0;
1342 if (vz > 0 && vy > 0) {
1344 thetayz = std::atan(ratio);
1347 else if (vz < 0 && vy > 0) {
1349 thetayz = std::atan(ratio);
1350 thetayz = TMath::Pi() - thetayz;
1353 else if (vz < 0 && vy < 0) {
1355 thetayz = std::atan(ratio);
1356 thetayz = thetayz + TMath::Pi();
1359 else if (vz > 0 && vy < 0) {
1361 thetayz = std::atan(ratio);
1362 thetayz = 2.0 * TMath::Pi() - thetayz;
1365 else if (vz == 0 && vy > 0) {
1366 thetayz = TMath::Pi() / 2.0;
1369 else if (vz == 0 && vy < 0) {
1370 thetayz = 3.0 * TMath::Pi() / 2.0;
1372 else if (vz > 0 && vy == 0) {
1375 else if (vz < 0 && vy == 0) {
1376 thetayz = TMath::Pi();
1379 if (thetayz > TMath::Pi()) { thetayz = thetayz - 2.0 * TMath::Pi(); }
1381 return 1000.0 * thetayz;
1387 cout <<
" Error : The code tries to divide by zero ! " << endl;
1391 double const arg = (xx - Q) / s;
1392 double const result = -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
1394 if (std::isnan(result) || std::isinf(result)) {
1395 cout <<
" Is nan ! my_g ! " << -std::log(s) <<
", " << s << endl;
1402 std::vector<double>
const& dEj,
1403 std::vector<double>
const& dthij,
1404 std::vector<double>
const& ind,
1406 double const x1)
const 1409 double theta0x =
x1;
1411 double result = 0.0;
1412 double nnn1 = dEi.size();
1417 for (
int i = 0; i < nnn1; i++) {
1418 double Ei = p - dEi.at(i);
1419 double Ej = p - dEj.at(i);
1422 if (Ei > 0 && Ej < 0) addth = 3.14 * 1000.0;
1430 (13.6 / std::sqrt(Ei * Ej)) * (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
1434 if (ind.at(i) == 1) {
1436 rms = std::hypot(tH0, theta0x);
1438 double const DT = dthij.at(i) + addth;
1439 double const prob =
my_g(DT, 0.0, rms);
1441 result = result - 2.0 * prob;
1445 if (std::isnan(result) || std::isinf(result)) {
1446 std::cout <<
" Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
double GetMomentumMultiScatterLLHD(art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const double min_resolution=0.001, const double max_resolution=800, const bool check_valid_scattered=false, const double angle_correction=0.757)
Calculate muon momentum (GeV) using multiple coulomb scattering by log likelihood.
double GetMomentumMultiScatterChi2(art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const double min_resolution=0, const double max_resolution=45)
Calculate muon momentum (GeV) using multiple coulomb scattering. Chi2 minimization of the Highland fo...
double my_g(double xx, double Q, double s) const
chi square minizer using Minuit2, it will minize (xx-Q)/s
Float_t y1[n_points_granero]
Float_t x1[n_points_granero]
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
constexpr auto abs(T v)
Returns the absolute value of the argument.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::optional< Segments > getSegTracks_(std::vector< double > const &xxx, std::vector< double > const &yyy, std::vector< double > const &zzz, double seg_size, bool check_valid_scattered=false)
Split tracks into segments to calculate the scattered angle later. Check DOI 10.1088/1748-0221/12/10/...
Float_t y2[n_points_geant4]
int getDeltaThetaij_(std::vector< double > &ei, std::vector< double > &ej, std::vector< double > &th, std::vector< double > &ind, Segments const &segments, double thick, double const angle_correction) const
Gets the scatterd angle for all the segments.
ScatterAngleMethods fMCSAngleMethod
Use scattered angle z-x (z is along the particle's direction)
double Length(size_t p=0) const
Access to various track properties.
double find_angle(double vz, double vy) const
Gets angle between two vy and vz.
double my_mcs_llhd(std::vector< double > const &dEi, std::vector< double > const &dEj, std::vector< double > const &dthij, std::vector< double > const &ind, double x0, double x1) const
Minimizer of log likelihood for scattered angle.
bool plotRecoTracks_(std::vector< double > const &xxx, std::vector< double > const &yyy, std::vector< double > const &zzz)
Struct to store segments. x, y and z are the 3D points of the segment nx, ny, nz forms the vector tha...
std::tuple< double, double, double > getDeltaThetaRMS_(Segments const &segments, double thick) const
Gets the scattered angle RMS for a all segments.
Use space angle: sqrt( zx^2 + zy^2 )/sqrt(2)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
void compute_max_fluctuation_vector(const std::vector< double > segx, const std::vector< double > segy, const std::vector< double > segz, std::vector< double > &segnx, std::vector< double > &segny, std::vector< double > &segnz, std::vector< bool > &segn_isvalid, std::vector< double > &vx, std::vector< double > &vy, std::vector< double > &vz, bool check_valid_scattered)
Computes the vector with most scattering inside a segment with size steps_size.
double GetMuMultiScatterLLHD3(art::Ptr< recob::Track > const &trk, bool dir)
Provides recob::Track data product.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
TVector3 GetMultiScatterStartingPoint(art::Ptr< recob::Track > const &trk)
Float_t x2[n_points_geant4]
TPolyLine3D * gr_reco_xyz
std::vector< double > steps
double GetTrackMomentum(double trkrange, int pdg) const