#include "TrackMomentumCalculator.h"

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

Public Member Functions
	TrackMomentumCalculator (double minLength=100.0, double maxLength=1350.0, double steps_size=10., int angleMethod=1, int nsteps=6)
	Constructor. More...

double	GetTrackMomentum (double trkrange, int pdg) const

double	GetMomentumMultiScatterChi2 (art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const double min_resolution=0, const double max_resolution=45)
	Calculate muon momentum (GeV) using multiple coulomb scattering. Chi2 minimization of the Highland formula. More...

double	GetMomentumMultiScatterLLHD (art::Ptr< recob::Track > const &trk, const bool checkValidPoints=false, const int maxMomentum_MeV=7500, const 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. More...

double	GetMuMultiScatterLLHD3 (art::Ptr< recob::Track > const &trk, bool dir)

TVector3	GetMultiScatterStartingPoint (art::Ptr< recob::Track > const &trk)

Private Types
enum	ScatterAngleMethods { kAnglezx = 1, kAnglezy, kAngleCombined }

Private Member Functions
bool	plotRecoTracks_ (std::vector< double > const &xxx, std::vector< double > const &yyy, std::vector< double > const &zzz)

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. More...

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/P10010. More...

std::tuple< double, double, double >	getDeltaThetaRMS_ (Segments const &segments, double thick) const
	Gets the scattered angle RMS for a all segments. More...

int	getDeltaThetaij_ (std::vector< 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. More...

double	my_g (double xx, double Q, double s) const
	chi square minizer using Minuit2, it will minize (xx-Q)/s More...

double	my_mcs_llhd (std::vector< 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. More...

double	find_angle (double vz, double vy) const
	Gets angle between two vy and vz. More...

Private Attributes
double	seg_stop {-1.}

int	n_seg {}

double	x_seg [50000]

double	y_seg [50000]

double	z_seg [50000]

std::vector< double >	steps

double	minLength

double	maxLength

double	steps_size

double	rad_length {14.0}

ScatterAngleMethods	fMCSAngleMethod

int	n_steps

TPolyLine3D *	gr_reco_xyz {nullptr}

TGraph	gr_reco_xy {}

TGraph	gr_reco_yz {}

TGraph	gr_reco_xz {}

TPolyLine3D *	gr_seg_xyz {nullptr}

TGraph	gr_seg_xy {}

TGraph	gr_seg_yz {}

TGraph	gr_seg_xz {}

Detailed Description

Definition at line 21 of file TrackMomentumCalculator.h.

Member Enumeration Documentation

enum trkf::TrackMomentumCalculator::ScatterAngleMethods

private

Enumerator
kAnglezx	Use scattered angle z-x (z is along the particle's direction)
kAnglezy	Use scattered angle z-y.
kAngleCombined	Use space angle: sqrt( zx^2 + zy^2 )/sqrt(2)

Definition at line 231 of file TrackMomentumCalculator.h.

                              {
       kAnglezx = 1,   
       kAnglezy,       
       kAngleCombined, 
     };

Constructor & Destructor Documentation

trkf::TrackMomentumCalculator::TrackMomentumCalculator	(	double	minLength = `100.0`,
		double	maxLength = `1350.0`,
		double	steps_size = `10.`,
		int	angleMethod = `1`,
		int	nsteps = `6`
	)

Constructor.

Parameters are relevant for Multiple Coulomb Scattering

Parameters

minLength	minimum length in cm of tracks (length here is based on the amout of segments)
maxLength	maximum length in cm of tracks (length here is based on the amout of segments)
steps_size	size in cm of each segment to compute scattering

Definition at line 302 of file TrackMomentumCalculator.cxx.

References fMCSAngleMethod, maxLength, n_steps, steps, and steps_size.

     : minLength{min}
     , maxLength{max}
     , steps_size{stepsize}
     , fMCSAngleMethod{static_cast<ScatterAngleMethods>(angleMethod)}
     , n_steps{nsteps}
   {
     for (int i = 1; i <= n_steps; i++) {
       steps.push_back(steps_size * i);
     }
   }

Member Function Documentation

void trkf::TrackMomentumCalculator::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
	)

private

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

Parameters

segx,segy,segz	segments points
segnx,segny,segnz	vector components to be filled
vector	used to control points to be used at segments

Definition at line 866 of file TrackMomentumCalculator.cxx.

References util::abs(), and n_seg.

Referenced by getSegTracks_().

   {
     auto const na = vx.size();
 
     double sumx = 0.0;
     double sumy = 0.0;
     double sumz = 0.0;
 
     bool isvalid = true;
     // In case vx, vy, etc have 3 points, probably two are just "linear"
     // interpolations. In this case, the angle of scattering will be zero.
     if (check_valid_scattered && na <= 2) { isvalid = false; }
 
     // computes the average in x, y, z
     for (std::size_t i = 0; i < na; ++i) {
       sumx += vx.at(i);
       sumy += vy.at(i);
       sumz += vz.at(i);
     }
 
     sumx /= na;
     sumy /= na;
     sumz /= na;
 
     std::vector<double> mx;
     std::vector<double> my;
     std::vector<double> mz;
 
     TMatrixDSym m(3);
 
     // Computes the Covariance matrix (Principal Component Analysis (PCA)
     for (std::size_t i = 0; i < na; ++i) {
       double const xxw1 = vx.at(i);
       double const yyw1 = vy.at(i);
       double const zzw1 = vz.at(i);
 
       mx.push_back(xxw1 - sumx);
       my.push_back(yyw1 - sumy);
       mz.push_back(zzw1 - sumz);
 
       double const xxw0 = mx.at(i);
       double const yyw0 = my.at(i);
       double const zzw0 = mz.at(i);
 
       m(0, 0) += xxw0 * xxw0 / na;
       m(0, 1) += xxw0 * yyw0 / na;
       m(0, 2) += xxw0 * zzw0 / na;
 
       m(1, 0) += yyw0 * xxw0 / na;
       m(1, 1) += yyw0 * yyw0 / na;
       m(1, 2) += yyw0 * zzw0 / na;
 
       m(2, 0) += zzw0 * xxw0 / na;
       m(2, 1) += zzw0 * yyw0 / na;
       m(2, 2) += zzw0 * zzw0 / na;
     }
 
     TMatrixDSymEigen me(m);
 
     // retrieve eigenvalues and vectors
     TVectorD eigenval = me.GetEigenValues();
     TMatrixD eigenvec = me.GetEigenVectors();
 
     double max1 = -666.0;
 
     double ind1 = 0;
 
     // get maximum eingevalue
     for (int i = 0; i < 3; ++i) {
       double const p1 = eigenval(i);
 
       if (p1 > max1) {
         max1 = p1;
         ind1 = i;
       }
     }
 
     // set the `direction` vector that points to the maximum fluctuation
     double ax = eigenvec(0, ind1);
     double ay = eigenvec(1, ind1);
     double az = eigenvec(2, ind1);
 
     if (n_seg > 1) {
       // for x, y and z, check if the last point is bigger then the previous point.
       // Ensures that the computed fluctation follows the trend of the track
       if (segx.at(n_seg - 1) - segx.at(n_seg - 2) > 0)
         ax = std::abs(ax);
       else
         ax = -1.0 * std::abs(ax);
 
       if (segy.at(n_seg - 1) - segy.at(n_seg - 2) > 0)
         ay = std::abs(ay);
       else
         ay = -1.0 * std::abs(ay);
 
       if (segz.at(n_seg - 1) - segz.at(n_seg - 2) > 0)
         az = std::abs(az);
       else
         az = -1.0 * std::abs(az);
 
       segnx.push_back(ax);
       segny.push_back(ay);
       segnz.push_back(az);
       segn_isvalid.push_back(isvalid);
     }
 
     // clear the vectors
     vx.clear();
     vy.clear();
     vz.clear();
   }

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

private

Gets angle between two vy and vz.

Parameters

vz
vy

Returns: angle in mrad

Definition at line 1378 of file TrackMomentumCalculator.cxx.

References util::abs().

Referenced by getDeltaThetaij_(), and getDeltaThetaRMS_().

   {
     double thetayz = -999.0;
 
     if (vz > 0 && vy > 0) {
       double ratio = std::abs(vy / vz);
       thetayz = std::atan(ratio);
     }
 
     else if (vz < 0 && vy > 0) {
       double ratio = std::abs(vy / vz);
       thetayz = std::atan(ratio);
       thetayz = TMath::Pi() - thetayz;
     }
 
     else if (vz < 0 && vy < 0) {
       double ratio = std::abs(vy / vz);
       thetayz = std::atan(ratio);
       thetayz = thetayz + TMath::Pi();
     }
 
     else if (vz > 0 && vy < 0) {
       double ratio = std::abs(vy / vz);
       thetayz = std::atan(ratio);
       thetayz = 2.0 * TMath::Pi() - thetayz;
     }
 
     else if (vz == 0 && vy > 0) {
       thetayz = TMath::Pi() / 2.0;
     }
 
     else if (vz == 0 && vy < 0) {
       thetayz = 3.0 * TMath::Pi() / 2.0;
     }
     else if (vz > 0 && vy == 0) {
       thetayz = 0;
     }
     else if (vz < 0 && vy == 0) {
       thetayz = TMath::Pi();
     }
 
     if (thetayz > TMath::Pi()) { thetayz = thetayz - 2.0 * TMath::Pi(); }
 
     return 1000.0 * thetayz;
   }

int trkf::TrackMomentumCalculator::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

private

Gets the scatterd angle for all the segments.

Parameters

ei,ej	will be filled with energy lost at point i and j
th	delta theta to be filled
ind	selection of scattering plane
segments
thick	is the steps_size
angle_correction	Correction due to oversmoothing, applied for space angle only

Returns: sucess or failure

Definition at line 611 of file TrackMomentumCalculator.cxx.

References a1, a2, a3, util::abs(), find_angle(), fMCSAngleMethod, kAngleCombined, kAnglezx, kAnglezy, trkf::TrackMomentumCalculator::Segments::L, trkf::TrackMomentumCalculator::Segments::nx, trkf::TrackMomentumCalculator::Segments::ny, and trkf::TrackMomentumCalculator::Segments::nz.

Referenced by GetMomentumMultiScatterLLHD(), and GetMuMultiScatterLLHD3().

   {
     auto const& segnx = segments.nx;
     auto const& segny = segments.ny;
     auto const& segnz = segments.nz;
     auto const& segL = segments.L;
 
     int const a1 = segnx.size();
     int const a2 = segny.size();
     int const a3 = segnz.size();
 
     if (a1 != a2 || a1 != a3) {
       std::cout << " ( Get thij ) Error ! " << std::endl;
       return -1.0;
     }
 
     int tot = a1;
 
     for (int i = 0; i < tot - 1; i++) {
       double const dx = segnx.at(i);
       double const dy = segny.at(i);
       double const dz = segnz.at(i);
 
       // Assumes z as propagation angle
       TVector3 const vec_z{dx, dy, dz};
       TVector3 vec_x;
       TVector3 vec_y;
 
       double const switcher = basex.Dot(vec_z);
       if (std::abs(switcher) <= 0.995) {
         vec_y = vec_z.Cross(basex).Unit();
         vec_x = vec_y.Cross(vec_z);
       }
       else {
         // cout << " It switched ! Isn't this lovely !!! " << endl; //
         vec_y = basez.Cross(vec_z).Unit();
         vec_x = vec_y.Cross(vec_z);
       }
 
       TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
       TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
       TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
 
       double const here_dx = segnx.at(i + 1);
       double const here_dy = segny.at(i + 1);
       double const here_dz = segnz.at(i + 1);
 
       TVector3 const here_vec{here_dx, here_dy, here_dz};
       TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
 
       double const scx = rot_here.X();
       double const scy = rot_here.Y();
       double const scz = rot_here.Z();
 
       double const azy = find_angle(scz, scy);
       double const azx = find_angle(scz, scx);
 
       constexpr double ULim = 10000.0; // Avoid huge (wrong) angles
       constexpr double LLim = -10000.0;
 
       double const Li = segL.at(i);
       double const Lj = segL.at(i + 1);
 
       if (azy <= ULim && azy >= LLim) { // safe scatter in the yz plane
 
         if (azx <= ULim && azx >= LLim) { // safe scatter in the za plane
           ei.push_back(Li);               // Energy deposited at i
           ej.push_back(Lj);               // Energy deposited at i+1
           if (fMCSAngleMethod == kAnglezx) {
             th.push_back(azx); // scattered angle z-x
           }
           else if (fMCSAngleMethod == kAnglezy) {
             th.push_back(azy);
           }
           else if (fMCSAngleMethod == kAngleCombined) {
             th.push_back(std::hypot(azx, azy) /
                          angle_correction); // space angle (applying correction)
           }
         }
       }
     }
 
     return 0;
   }

std::tuple< double, double, double > trkf::TrackMomentumCalculator::getDeltaThetaRMS_	(	Segments const &	segments,
		double	thick
	)		const

private

Gets the scattered angle RMS for a all segments.

Parameters

segments	segments computed
thick	is the steps_size

TODO: Add better description of steps

Returns: tuple with mean value, rms and error of rms

Definition at line 1245 of file TrackMomentumCalculator.cxx.

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

Referenced by GetMomentumMultiScatterChi2().

   {
     auto const& segnx = segments.nx;
     auto const& segny = segments.ny;
     auto const& segnz = segments.nz;
     auto const& segL = segments.L;
 
     int const a1 = segnx.size();
     int const a2 = segny.size();
     int const a3 = segnz.size();
 
     if (a1 != a2 || a1 != a3) {
       cout << " ( Get RMS ) Error ! " << endl;
       return std::make_tuple(0., 0., 0.);
     }
 
     int const tot = a1;
 
     std::vector<double> buf0;
 
     for (int i = 0; i < tot; i++) {
       double const dx = segnx.at(i);
       double const dy = segny.at(i);
       double const dz = segnz.at(i);
 
       TVector3 const vec_z{dx, dy, dz};
       TVector3 vec_x;
       TVector3 vec_y;
 
       double const switcher = basex.Dot(vec_z);
 
       if (switcher <= 0.995) {
         vec_y = vec_z.Cross(basex).Unit();
         vec_x = vec_y.Cross(vec_z);
       }
       else {
         // cout << " It switched ! Isn't this lovely !!! " << endl;
         vec_y = basez.Cross(vec_z).Unit();
         vec_x = vec_y.Cross(vec_z);
       }
 
       TVector3 const Rx{vec_x.Dot(basex), vec_x.Dot(basey), vec_x.Dot(basez)};
       TVector3 const Ry{vec_y.Dot(basex), vec_y.Dot(basey), vec_y.Dot(basez)};
       TVector3 const Rz{vec_z.Dot(basex), vec_z.Dot(basey), vec_z.Dot(basez)};
 
       double const refL = segL.at(i);
 
       for (int j = i; j < tot; j++) {
         double const L1 = segL.at(j);
 
         double const dz1 = L1 - refL;
 
         if (dz1 >= thick) {
           double const here_dx = segnx.at(j);
           double const here_dy = segny.at(j);
           double const here_dz = segnz.at(j);
 
           TVector3 const here_vec{here_dx, here_dy, here_dz};
           TVector3 const rot_here{Rx.Dot(here_vec), Ry.Dot(here_vec), Rz.Dot(here_vec)};
 
           double const scx = rot_here.X();
           double const scy = rot_here.Y();
           double const scz = rot_here.Z();
 
           double const azy = find_angle(scz, scy);
           double const azx = find_angle(scz, scx);
 
           constexpr double ULim = 10000.0; // Avoid huge (wrong) angles
           constexpr double LLim = -10000.0;
 
           if (azy <= ULim && azy >= LLim) { // safe scatter in the yz plane
 
             if (azx <= ULim && azx >= LLim) { // safe scatter in the za plane
               if (fMCSAngleMethod == kAnglezx) {
                 buf0.push_back(azx); // scattered angle z-x
               }
               else if (fMCSAngleMethod == kAnglezy) {
                 buf0.push_back(azy);
               }
               else if (fMCSAngleMethod == kAngleCombined) {
                 buf0.push_back(std::hypot(azx, azy));
               }
             }
           }
           break; // of course !
         }
       }
     }
 
     int const nmeas = buf0.size();
     double nnn = 0.0;
 
     double mean = 0.0;
     double rms = 0.0;
     double rmse = 0.0;
 
     for (int i = 0; i < nmeas; i++) {
       mean += buf0.at(i);
       nnn++;
     }
 
     mean = mean / nnn;
 
     for (int i = 0; i < nmeas; i++)
       rms += ((buf0.at(i)) * (buf0.at(i)));
 
     rms = rms / (nnn);
     rms = std::sqrt(rms);
     rmse = rms / std::sqrt(2.0 * tot);
 
     double rms1 = rms;
 
     rms = 0.0;
 
     double ntot1 = 0.0;
     double const lev1 = 2.50;
 
     for (int i = 0; i < nmeas; i++) {
       double const amp = buf0.at(i);
       if (amp < (mean + lev1 * rms1) && amp > (mean - lev1 * rms1)) {
         ++ntot1;
         rms += amp * amp;
       }
     }
 
     rms = rms / (ntot1);
     rms = std::sqrt(rms);
     rmse = rms / std::sqrt(2.0 * ntot1);
     return std::make_tuple(mean, rms, rmse);
   }

double trkf::TrackMomentumCalculator::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 formula.

Parameters

trk	the muon track
checkValidPoints	rather take into account only valid points or not
maxMomentum_MeV	maximum momentum in MeV for the minimization

TODO: Add better description of the steps done.

Returns: momentum in GeV

Definition at line 702 of file TrackMomentumCalculator.cxx.

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

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

   {
     std::vector<double> recoX;
     std::vector<double> recoY;
     std::vector<double> recoZ;
 
     int n_points = trk->NumberTrajectoryPoints();
 
     for (int i = 0; i < n_points; i++) {
       if (checkValidPoints && !trk->HasValidPoint(i)) continue;
       auto const& pos = trk->LocationAtPoint(i);
       recoX.push_back(pos.X());
       recoY.push_back(pos.Y());
       recoZ.push_back(pos.Z());
     }
 
     if (recoX.size() < 2) return -1.0;
 
     if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
 
     double const seg_size{steps_size};
     auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
     if (!segments.has_value()) return -1.0;
 
     auto const seg_steps = segments->x.size();
     if (seg_steps < 2) return -1;
 
     double const recoSegmentLength = segments->L.at(seg_steps - 1);
     if (recoSegmentLength < minLength || recoSegmentLength > maxLength) return -1;
 
     double ymax = -999.0;
     double ymin = +999.0;
 
     std::vector<double> xmeas;
     std::vector<double> ymeas;
     std::vector<double> eymeas;
     xmeas.reserve(n_steps);
     ymeas.reserve(n_steps);
     eymeas.reserve(n_steps);
     for (int j = 0; j < n_steps; j++) {
       double const trial = steps.at(j);
       // computes the rms by groups of trial, if seg_size was chosen as 10, trials will be 10, 20, etc.. until 10 * n_steps
       auto const [mean, rms, rmse] = getDeltaThetaRMS_(*segments, trial);
 
       if (std::isnan(mean) || std::isinf(mean)) {
         mf::LogDebug("TrackMomentumCalculator") << "Returned mean is either nan or infinity.";
         continue;
       }
       if (std::isnan(rms) || std::isinf(rms)) {
         mf::LogDebug("TrackMomentumCalculator") << "Returned rms is either nan or infinity.";
         continue;
       }
       if (std::isnan(rmse) || std::isinf(rmse)) {
         mf::LogDebug("TrackMomentumCalculator") << "Returned rmse is either nan or infinity.";
         continue;
       }
 
       if (mean == -1 && rms == -1 && rmse == -1) continue;
       xmeas.push_back(trial); // x values are different steps length, ex: 10, 20, 30 cm
       ymeas.push_back(rms);   // y values are the RMS of the scattered angle for each step defined
       eymeas.push_back(
         std::hypot(rmse, 0.05 * rms)); // <--- conservative syst. error to fix chi^{2} behaviour !!!
 
       if (ymin > rms) ymin = rms;
       if (ymax < rms) ymax = rms;
     }
 
     assert(xmeas.size() == ymeas.size());
     assert(xmeas.size() == eymeas.size());
     if (xmeas.empty()) { return -1.0; }
 
     TGraphErrors gr_meas{n_steps, xmeas.data(), ymeas.data(), nullptr, eymeas.data()};
 
     gr_meas.SetTitle("(#Delta#theta)_{rms} versus material thickness; Material "
                      "thickness in cm; (#Delta#theta)_{rms} in mrad");
 
     gr_meas.SetLineColor(kBlack);
     gr_meas.SetMarkerColor(kBlack);
     gr_meas.SetMarkerStyle(20);
     gr_meas.SetMarkerSize(1.2);
 
     gr_meas.GetXaxis()->SetLimits(steps.at(0) - steps.at(0), steps.at(n_steps - 1) + steps.at(0));
     gr_meas.SetMinimum(0.0);
     gr_meas.SetMaximum(1.80 * ymax);
 
     double correction = 1.;
     if (fMCSAngleMethod == kAngleCombined) { correction = std::sqrt(2.); }
 
     ROOT::Minuit2::Minuit2Minimizer mP{};
     FcnWrapper const wrapper{std::move(xmeas), std::move(ymeas), std::move(eymeas), correction};
     ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_chi2(xs); }, 2);
 
     // Start point for resolution
     double startpoint = 2;
     if (startpoint < min_resolution) startpoint = (max_resolution - min_resolution) / 2.;
     bool fixresolution = false;
     double maxres = max_resolution;
     if (max_resolution == 0 || max_resolution <= min_resolution) {
       fixresolution = true;
       startpoint = min_resolution;
       maxres += 1;
     }
 
     mP.SetFunction(FCA);
     mP.SetLimitedVariable(0, "p_{MCS}", 1.0, 0.01, 0.001, maxMomentum_MeV / 1.e3);
     mP.SetLimitedVariable(1, "#delta#theta", startpoint, startpoint / 2., min_resolution, maxres);
     if (fixresolution) { mP.FixVariable(1); }
     mP.SetMaxFunctionCalls(1.E9);
     mP.SetMaxIterations(1.E9);
     mP.SetTolerance(0.01);
     mP.SetStrategy(2);
     mP.SetErrorDef(1);
 
     bool const mstatus = mP.Minimize();
 
     mP.Hesse();
 
     const double* pars = mP.X();
 
     auto const recoL = trk->Length();
     double const deltap = (recoL * kcal) / 2.0;
 
     double const p_mcs = pars[0] + deltap;
 
     return mstatus ? p_mcs : -1.0;
   }

double trkf::TrackMomentumCalculator::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.

Parameters

trk	the muon track
checkValidPoints	rather take into account only valid points or not
maxMomentum_MeV	maximum momentum in MeV for the minimization
MomentumStep_MeV	energy steps for minimization
max_resolution	maximum angular resolution for fit. Setting to zero will cause the fit only over momentum and fixed resolution of 2 mrad
check_valid_scattered	set to true will check if scatter angles are valid. Angles are invalid if there is only two points in one segment.
angle_correction	Correction for space angle due to possible oversmoothing. The value (0.757) was set based on studies with MC. Change this value through the fhcl file

TODO: Add better description of the steps done

Returns: momentum in GeV

Definition at line 445 of file TrackMomentumCalculator.cxx.

References fMCSAngleMethod, getDeltaThetaij_(), getSegTracks_(), GetTrackMomentum(), recob::Track::HasValidPoint(), kAngleCombined, recob::Track::Length(), recob::Track::LocationAtPoint(), maxLength, recob::Track::NumberTrajectoryPoints(), plotRecoTracks_(), and steps_size.

   {
 
     std::vector<double> recoX;
     std::vector<double> recoY;
     std::vector<double> recoZ;
 
     int n_points = trk->NumberTrajectoryPoints();
 
     for (int i = 0; i < n_points; i++) {
       if (checkValidPoints && !trk->HasValidPoint(i)) continue;
       auto const& pos = trk->LocationAtPoint(i);
       recoX.push_back(pos.X());
       recoY.push_back(pos.Y());
       recoZ.push_back(pos.Z());
     }
 
     if (recoX.size() < 2) return -1.0;
 
     if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
 
     double const seg_size{steps_size};
 
     auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size, check_valid_scattered);
     if (!segments.has_value()) return -1.0;
 
     auto const seg_steps = segments->x.size();
     if (seg_steps < 2) return -1;
 
     double const recoSegmentLength = segments->L.at(seg_steps - 1);
     if (recoSegmentLength < minLength || recoSegmentLength > maxLength) return -1;
 
     std::vector<double> dEi;
     std::vector<double> dEj;
     std::vector<double> dthij;
     std::vector<double> ind;
     std::vector<bool> dthij_valid = segments->nvalid;
     if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
       return -1;
 
     auto const ndEi = dEi.size();
     if (ndEi < 1) return -1;
 
     double correction = 1.;
     if (fMCSAngleMethod == kAngleCombined) { correction = std::sqrt(2.); }
 
     // Gets energy evaluated by range
     auto recoL = trk->Length();
     double const minP = this->GetTrackMomentum(recoL, 13);
 
     ROOT::Minuit2::Minuit2Minimizer mP{};
     FcnWrapperLLHD const wrapper{std::move(dEi),
                                  std::move(dEj),
                                  std::move(dthij),
                                  std::move(ind),
                                  std::move(dthij_valid),
                                  seg_size,
                                  correction};
     ROOT::Math::Functor FCA([&wrapper](double const* xs) { return wrapper.my_mcs_llhd(xs); }, 2);
 
     mP.SetFunction(FCA);
 
     // Start point for resolution
     // If max_resolution is zero, startpoint is equal min_resolution
     double startpoint = 2;
     if (startpoint < min_resolution)
       startpoint = (max_resolution - min_resolution) / 2.; // prevent wrong values
     bool fixresolution = false;
     double maxres = max_resolution;
     if (max_resolution == 0 || max_resolution <= min_resolution) {
       fixresolution = true;
       startpoint = min_resolution;
       maxres += 1;
     }
 
     // Starting energy as double of the energy by range
     // Assumes that the smallet possible energy is given by 60%  with CSDA
     // Step as 10 %
     mP.SetLimitedVariable(0, "p_{MCS}", minP * 2, minP * 0.1, minP * 0.6, maxMomentum_MeV / 1.e3);
     mP.SetLimitedVariable(1, "#delta#theta", startpoint, startpoint / 2., min_resolution, maxres);
     if (fixresolution) { mP.FixVariable(1); }
     mP.SetMaxFunctionCalls(1.E9);
     mP.SetMaxIterations(1.E9);
     mP.SetTolerance(0.01);
     mP.SetStrategy(2);
     mP.SetErrorDef(1);
 
     bool const mstatus = mP.Minimize();
 
     mP.Hesse();
 
     const double* pars = mP.X();
 
     double const p_mcs = pars[0];
 
     return mstatus ? p_mcs : -1.0;
   }

TVector3 trkf::TrackMomentumCalculator::GetMultiScatterStartingPoint ( art::Ptr< recob::Track > const & trk )

Definition at line 549 of file TrackMomentumCalculator.cxx.

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

   {
     double const LLHDp = GetMuMultiScatterLLHD3(trk, true);
     double const LLHDm = GetMuMultiScatterLLHD3(trk, false);
 
     if (LLHDp != -1 && LLHDm != -1 && LLHDp > LLHDm) {
       int const n_points = trk->NumberTrajectoryPoints();
       return trk->LocationAtPoint<TVector3>(n_points - 1);
     }
     else {
       return trk->LocationAtPoint<TVector3>(0);
     }
 
     // Should never get here
     return TVector3{};
   }

double trkf::TrackMomentumCalculator::GetMuMultiScatterLLHD3	(	art::Ptr< recob::Track > const &	trk,
		bool	dir
	)

Definition at line 566 of file TrackMomentumCalculator.cxx.

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

Referenced by GetMultiScatterStartingPoint().

   {
     std::vector<double> recoX;
     std::vector<double> recoY;
     std::vector<double> recoZ;
 
     int const n_points = trk->NumberTrajectoryPoints();
     for (int i = 0; i < n_points; ++i) {
       auto const index = dir ? i : n_points - 1 - i;
       auto const& pos = trk->LocationAtPoint(index);
       recoX.push_back(pos.X());
       recoY.push_back(pos.Y());
       recoZ.push_back(pos.Z());
     }
 
     if (recoX.size() < 2) return -1.0;
 
     if (!plotRecoTracks_(recoX, recoY, recoZ)) return -1.0;
 
     constexpr double seg_size{5.0};
     auto const segments = getSegTracks_(recoX, recoY, recoZ, seg_size);
     if (!segments.has_value()) return -1.0;
 
     auto const seg_steps = segments->x.size();
     if (seg_steps < 2) return -1;
 
     double const recoSegmentLength = segments->L.at(seg_steps - 1);
     if (recoSegmentLength < 15.0 || recoSegmentLength > maxLength) return -1;
 
     std::vector<double> dEi;
     std::vector<double> dEj;
     std::vector<double> dthij;
     std::vector<double> ind;
     double angle_correction = 1; // never tested in this function
     if (getDeltaThetaij_(dEi, dEj, dthij, ind, *segments, seg_size, angle_correction) != 0)
       return -1;
 
     auto const recoL = trk->Length();
     double const p_range = recoL * kcal;
     double const logL = my_mcs_llhd(dEi, dEj, dthij, ind, p_range, 5.65);
 
     return logL;
   }

std::optional< TrackMomentumCalculator::Segments > trkf::TrackMomentumCalculator::getSegTracks_	(	std::vector< double > const &	xxx,
		std::vector< double > const &	yyy,
		std::vector< double > const &	zzz,
		double	seg_size,
		bool	check_valid_scattered = `false`
	)

private

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

Parameters

xxx	3D reconstructed points x-axis
yyy	3D reconstructed points y-axiy
zzz	3D reconstructed points z-axiz
seg_size	Segments size defined in class constructor
check_valid_scattered	set to true will check if scatter angles are valid. Angles are invalid if there is only two points in one segment. Set true only for LLHD! TODO: Add better description of steps

Returns: Segments

Definition at line 994 of file TrackMomentumCalculator.cxx.

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

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

   {
     double stag = 0.0;
 
     int a1 = xxx.size();
     int a2 = yyy.size();
     int a3 = zzz.size();
 
     if ((a1 != a2) || (a1 != a3) || (a2 != a3)) {
       cout << " ( Digitize reco tacks ) Error ! " << endl;
       return std::nullopt;
     }
 
     int const stopper = seg_stop / seg_size;
 
     // values to be filled and returned
     std::vector<double> segx, segnx;
     std::vector<double> segy, segny;
     std::vector<double> segz, segnz;
     std::vector<bool> segn_isvalid;
     std::vector<double> segL;
 
     n_seg = 0;
 
     double x0{};
     double y0{};
     double z0{};
 
     double x00 = xxx.at(0);
     double y00 = yyy.at(0);
     double z00 = zzz.at(0);
 
     int indC = 0;
 
     // These vectors will keep the points inside reach segments
     // They are cleared inside the function `compute_max_fluctuation_vector`
     // each time it finishes with a segment
     std::vector<double> vx;
     std::vector<double> vy;
     std::vector<double> vz;
 
     for (int i = 0; i < a1; i++) {
       x0 = xxx.at(i);
       y0 = yyy.at(i);
       z0 = zzz.at(i);
 
       double const RR0 = std::hypot(x00 - x0, y00 - y0, z00 - z0);
 
       if (RR0 >= stag) { // stag is aways set to zero, this is always true
 
         segx.push_back(x0);
         segy.push_back(y0);
         segz.push_back(z0);
 
         segL.push_back(0);
 
         // TGraph
         x_seg[n_seg] = x0;
         y_seg[n_seg] = y0;
         z_seg[n_seg] = z0;
 
         n_seg++;
 
         vx.push_back(x0);
         vy.push_back(y0);
         vz.push_back(z0);
 
         indC = i + 1;
 
         break;
       }
     }
 
     for (int i = indC; i < a1 - 1; i++) { // starting at second point (i=1) if stag set to zero
       // current point
       double const x1 = xxx.at(i);
       double const y1 = yyy.at(i);
       double const z1 = zzz.at(i);
 
       // distante from previous point
       double const dr1 = std::hypot(x1 - x0, y1 - y0, z1 - z0);
 
       // next point
       double const x2 = xxx.at(i + 1);
       double const y2 = yyy.at(i + 1);
       double const z2 = zzz.at(i + 1);
 
       // distant of next point to previous point
       double const dr2 = std::hypot(x2 - x0, y2 - y0, z2 - z0);
 
       if (dr1 < seg_size) {
         vx.push_back(x1);
         vy.push_back(y1);
         vz.push_back(z1);
       }
 
       /* If current point is inside segment length w.r.t. the first point of
        * the segment (x0,y0,z0), but the next point is outsize: create a new
        * point in between (x1,y1,z1) and (x2,y2,z2) in which will be exacly at
        * the segment length (w.r.t. to the first point) This is done using the
        * cosines law for a given factor `t` times dr (x2-x1, ...), and so:
        *
        * (ds)^2 = (dr1)^2 + (t*dr)^2 + 2*dot_product(dr1, t*dr)
        * Using cos(180-theta) = -cos(theta))
        *
        * Solve for `t` the second degree equation: t^2 + beta*t + gamma = 0
        */
       if (dr1 <= seg_size && dr2 > seg_size) {
         double const dx = x2 - x1;
         double const dy = y2 - y1;
         double const dz = z2 - z1;
         double const dr = std::hypot(dx, dy, dz);
 
         if (dr == 0) {
           cout << " ( Zero ) Error ! " << endl;
           return std::nullopt;
         }
 
         double const beta = 2.0 * ((x1 - x0) * dx + (y1 - y0) * dy + (z1 - z0) * dz) / (dr * dr);
 
         double const gamma = (dr1 * dr1 - seg_size * seg_size) / (dr * dr);
         double const delta = beta * beta - 4.0 * gamma;
 
         if (delta < 0.0) {
           cout << " ( Discriminant ) Error ! " << endl;
           return std::nullopt;
         }
 
         // solves for t
         double const lysi1 = (-beta + std::sqrt(delta)) / 2.0;
         double const t = lysi1;
 
         // find next points in that will exactly at the segment length from x0
         double const xp = x1 + t * dx;
         double const yp = y1 + t * dy;
         double const zp = z1 + t * dz;
 
         // Add points to be returned
         segx.push_back(xp);
         segy.push_back(yp);
         segz.push_back(zp);
 
         segL.push_back(n_seg * seg_size);
 
         // for TGraph
         x_seg[n_seg] = xp;
         y_seg[n_seg] = yp;
         z_seg[n_seg] = zp;
 
         n_seg++;
 
         // This are the new `x0` points (for next segment)
         x0 = xp;
         y0 = yp;
         z0 = zp;
 
         // Add points to segment
         vx.push_back(x0);
         vy.push_back(y0);
         vz.push_back(z0);
 
         if (n_seg <= 1) // This should never happen
           return std::nullopt;
 
         // Now, compute the deviation in `segx, ...` of the segment
         // vx, vy, vz are used and cleared afterwards
         compute_max_fluctuation_vector(
           segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
 
         // Starting over
         vx.push_back(x0);
         vy.push_back(y0);
         vz.push_back(z0);
       }
       else if (dr1 > seg_size) { // in this case, just interpolate until reach `seg_size`
 
         // Rolling `i` back for the next iteration.
         // Because the current point does not belong to this segment
         i = (i - 1);
 
         double const dx = x1 - x0;
         double const dy = y1 - y0;
         double const dz = z1 - z0;
         double const dr = std::hypot(dx, dy, dz);
 
         if (dr == 0) {
           cout << " ( Zero ) Error ! " << endl;
           return std::nullopt;
         }
 
         // computes the point by simple interpolation
         double const t = seg_size / dr;
         double const xp = x0 + t * dx;
         double const yp = y0 + t * dy;
         double const zp = z0 + t * dz;
 
         segx.push_back(xp);
         segy.push_back(yp);
         segz.push_back(zp);
         segL.push_back(1.0 * n_seg * 1.0 * seg_size + stag);
 
         // for TGraph
         x_seg[n_seg] = xp;
         y_seg[n_seg] = yp;
         z_seg[n_seg] = zp;
 
         n_seg++;
 
         x0 = xp;
         y0 = yp;
         z0 = zp;
 
         vx.push_back(x0);
         vy.push_back(y0);
         vz.push_back(z0);
 
         if (n_seg <= 1) // This should never happen
           return std::nullopt;
 
         // Now, compute the deviation in `segx, ...` of the segment
         // vx, vy, vz are used and cleared afterwards
         compute_max_fluctuation_vector(
           segx, segy, segz, segnx, segny, segnz, segn_isvalid, vx, vy, vz, check_valid_scattered);
 
         // vectors are cleared in previous step
         vx.push_back(x0);
         vy.push_back(y0);
         vz.push_back(z0);
       }
 
       if (n_seg >= (stopper + 1.0) && seg_stop != -1) break;
     }
 
     delete gr_seg_xyz;
     gr_seg_xyz = new TPolyLine3D{n_seg, z_seg, x_seg, y_seg};
     gr_seg_yz = TGraph{n_seg, z_seg, y_seg};
     gr_seg_xz = TGraph{n_seg, z_seg, x_seg};
     gr_seg_xy = TGraph{n_seg, x_seg, y_seg};
 
     return std::make_optional<Segments>(
       Segments{segx, segnx, segy, segny, segz, segnz, segL, segn_isvalid});
   }

double trkf::TrackMomentumCalculator::GetTrackMomentum	(	double	trkrange,
		int	pdg
	)		const

Definition at line 318 of file TrackMomentumCalculator.cxx.

References util::abs(), and E.

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

   {
     /* Muon range-momentum tables from CSDA (Argon density = 1.4 g/cm^3)
        website:
        http://pdg.lbl.gov/2012/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_289.pdf
 
        CSDA table values:
        double Range_grampercm[30] = {9.833E-1, 1.786E0, 3.321E0,
        6.598E0, 1.058E1, 3.084E1, 4.250E1, 6.732E1, 1.063E2, 1.725E2,
        2.385E2, 4.934E2, 6.163E2, 8.552E2, 1.202E3, 1.758E3, 2.297E3,
        4.359E3, 5.354E3, 7.298E3, 1.013E4, 1.469E4, 1.910E4, 3.558E4,
        4.326E4, 5.768E4, 7.734E4, 1.060E5, 1.307E5}; double KE_MeV[30] = {10, 14,
        20, 30, 40, 80, 100, 140, 200, 300, 400, 800, 1000, 1400, 2000, 3000,
        4000, 8000, 10000, 14000, 20000, 30000, 40000, 80000, 100000, 140000,
        200000, 300000, 400000};
 
        Functions below are obtained by fitting polynomial fits to KE_MeV vs
        Range (cm) graph. A better fit was obtained by splitting the graph into
        two: Below range<=200cm,a polynomial of power 4 was a good fit; above
        200cm, a polynomial of power 6 was a good fit
 
        Fit errors for future purposes:
        Below 200cm, Forpoly4 fit: p0 err=1.38533;p1 err=0.209626; p2
        err=0.00650077; p3 err=6.42207E-5; p4 err=1.94893E-7; Above 200cm,
        Forpoly6 fit: p0 err=5.24743;p1 err=0.0176229; p2 err=1.6263E-5; p3
        err=5.9155E-9; p4 err=9.71709E-13; p5 err=7.22381E-17;p6
        err=1.9709E-21;*/
 
     //*********For muon, the calculations are valid up to 1.91E4 cm range
     //corresponding to a Muon KE of 40 GeV**********//
 
     /*Proton range-momentum tables from CSDA (Argon density = 1.4 g/cm^3):
       website: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html
 
       CSDA values:
       double KE_MeV_P_Nist[31]={10, 15, 20, 30, 40, 80, 100, 150, 200, 250, 300,
       350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
       1500, 2000, 2500, 3000, 4000, 5000};
 
       double Range_gpercm_P_Nist[31]={1.887E-1,3.823E-1, 6.335E-1, 1.296,
       2.159, 7.375, 1.092E1, 2.215E1, 3.627E1, 5.282E1, 7.144E1,
       9.184E1, 1.138E2, 1.370E2, 1.614E2, 1.869E2, 2.132E2, 2.403E2,
       2.681E2, 2.965E2, 3.254E2, 3.548E2, 3.846E2, 4.148E2, 4.454E2,
       7.626E2, 1.090E3, 1.418E3, 1.745E3, 2.391E3, 3.022E3};
 
       Functions below are obtained by fitting power and polynomial fits to
       KE_MeV vs Range (cm) graph. A better fit was obtained by splitting the
       graph into two: Below range<=80cm,a a*(x^b) was a good fit; above 80cm, a
       polynomial of power 6 was a good fit
 
       Fit errors for future purposes:
       For power function fit: a=0.388873; and b=0.00347075
       Forpoly6 fit: p0 err=3.49729;p1 err=0.0487859; p2 err=0.000225834; p3
       err=4.45542E-7; p4 err=4.16428E-10; p5 err=1.81679E-13;p6
       err=2.96958E-17;*/
 
     //*********For proton, the calculations are valid up to 3.022E3 cm range
     //corresponding to a Muon KE of 5 GeV**********//
 
     /* Pion range-momentum extracted from Bethe-Bloch equation.  The
        computation are all one using
        https://github.com/sungbinoh/PDSPTreeAnalyzer.git (last commit 4a3d1d2)
 
        The computation is done with the spline of the Range_grampercm and
        KE_MeV_Pion.*/
 
     //*********For pion, the calculations are valid only when there is no
     //inelastic scattering present. So this method should be used
     //carefully.**********//
 
     if (trkrange < 0 || std::isnan(trkrange)) {
       mf::LogError("TrackMomentumCalculator") << "Invalid track range " << trkrange << " return -1";
       return -1.;
     }
 
     double KE, Momentum, M;
     constexpr double Muon_M = m_muon * 1e3, Proton_M = 938.272, Pion_M = 139.57039;
 
     if (abs(pdg) == 13) {
       M = Muon_M;
       KE = KEvsR_spline3.Eval(trkrange);
     }
     else if (abs(pdg) == 2212) {
       M = Proton_M;
       if (trkrange > 0 && trkrange <= 80)
         KE = 29.9317 * std::pow(trkrange, 0.586304);
       else if (trkrange > 80 && trkrange <= 3.022E3)
         KE = 149.904 + (3.34146 * trkrange) + (-0.00318856 * trkrange * trkrange) +
              (4.34587E-6 * trkrange * trkrange * trkrange) +
              (-3.18146E-9 * trkrange * trkrange * trkrange * trkrange) +
              (1.17854E-12 * trkrange * trkrange * trkrange * trkrange * trkrange) +
              (-1.71763E-16 * trkrange * trkrange * trkrange * trkrange * trkrange * trkrange);
       else
         KE = -999;
     }
     else if (abs(pdg) == 211) {
       M = Pion_M;
       KE = KEvsRPi_spline3.Eval(trkrange);
     }
     else
       KE = -999;
 
     if (KE < 0)
       Momentum = -999;
     else
       Momentum = std::sqrt((KE * KE) + (2 * M * KE));
 
     Momentum = Momentum / 1000;
 
     return Momentum;
   }

double trkf::TrackMomentumCalculator::my_g	(	double	xx,
		double	Q,
		double	s
	)		const

private

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

Parameters

xx
Q
s

Returns: Momentum in GeV

Definition at line 1424 of file TrackMomentumCalculator.cxx.

Referenced by my_mcs_llhd().

   {
     if (s == 0.) {
       cout << " Error : The code tries to divide by zero ! " << endl;
       return 0.;
     }
 
     double const arg = (xx - Q) / s;
     double const result = -0.5 * std::log(2.0 * TMath::Pi()) - std::log(s) - 0.5 * arg * arg;
 
     if (std::isnan(result) || std::isinf(result)) {
       cout << " Is nan ! my_g ! " << -std::log(s) << ", " << s << endl;
     }
 
     return result;
   }

double trkf::TrackMomentumCalculator::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

private

Minimizer of log likelihood for scattered angle.

Parameters

dEi	energy at step i
dEj	energy at step j
dthij	scattered angle between points i and j
ind	selection of scattering plane
x0	momentum to be fitted
x1	resolution to be fitted

TODO: Add better description of steps

Returns: momentum in GeV

Definition at line 1441 of file TrackMomentumCalculator.cxx.

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

Referenced by GetMuMultiScatterLLHD3().

   {
     double p = x0;
     double theta0x = x1;
 
     double result = 0.0;
     double nnn1 = dEi.size(); // number of segments of energy
 
     double red_length = (steps_size) / rad_length;
     double addth = 0;
 
     for (int i = 0; i < nnn1; i++) {
       double Ei = p - dEi.at(i); // Estimated energy at point i
       double Ej = p - dEj.at(i); // Estimated enery at point j
 
       // If the momentum p choosen allows that the muon stopped inside, add 1 rad to the change in scatter angle (as the particle stops)
       if (Ei > 0 && Ej < 0) addth = 3.14 * 1000.0;
 
       Ei = std::abs(Ei);
       Ej = std::abs(Ej);
 
       // Highland formula
       // Parameters given at Particle Data Group https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-passage-particles-matter.pdf
       double tH0 =
         (13.6 / std::sqrt(Ei * Ej)) * (1.0 + 0.038 * std::log(red_length)) * std::sqrt(red_length);
 
       double rms = -1.0;
 
       if (ind.at(i) == 1) {
         // Computes the rms of angle
         rms = std::hypot(tH0, theta0x);
 
         double const DT = dthij.at(i) + addth;
         double const prob = my_g(DT, 0.0, rms); // Computes log likelihood
 
         result = result - 2.0 * prob; // Adds for each segment
       }
     }
 
     if (std::isnan(result) || std::isinf(result)) {
       std::cout << " Is nan ! my_mcs_llhd ( 1 ) ! " << std::endl;
     }
     return result;
   }

bool trkf::TrackMomentumCalculator::plotRecoTracks_	(	std::vector< double > const &	xxx,
		std::vector< double > const &	yyy,
		std::vector< double > const &	zzz
	)

private

Definition at line 833 of file TrackMomentumCalculator.cxx.

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

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

   {
     auto const n = xxx.size();
     auto const y_size = yyy.size();
     auto const z_size = zzz.size();
 
     if (n != y_size || n != z_size) {
       cout << " ( Get reco tacks ) Error ! " << endl;
       return false;
     }
 
     // Here, we perform a const-cast to double* because, sadly,
     // TPolyLine3D requires a pointer to a non-const object.  We will
     // trust that ROOT does not mess around with the underlying data.
     auto xs = const_cast<double*>(xxx.data());
     auto ys = const_cast<double*>(yyy.data());
     auto zs = const_cast<double*>(zzz.data());
 
     auto const narrowed_size = static_cast<int>(n); // ROOT requires a signed integral type
     delete gr_reco_xyz;
     gr_reco_xyz = new TPolyLine3D{narrowed_size, zs, xs, ys};
     gr_reco_yz = TGraph{narrowed_size, zzz.data(), yyy.data()};
     gr_reco_xz = TGraph{narrowed_size, zzz.data(), xxx.data()};
     gr_reco_xy = TGraph{narrowed_size, xxx.data(), yyy.data()};
 
     return true;
   }