TrackMomentumCalculator.h

Go to the documentation of this file.

//  \author  sowjanyag@phys.ksu.edu

#ifndef TrackMomentumCalculator_H
#define TrackMomentumCalculator_H

#include "canvas/Persistency/Common/Ptr.h"
#include "lardataobj/RecoBase/Track.h"

#include "TGraph.h"
#include "TVector3.h"

#include <optional>
#include <tuple>
#include <vector>

class TPolyLine3D;

namespace trkf {

  class TrackMomentumCalculator {
  public:
    TrackMomentumCalculator(double minLength = 100.0,
                            double maxLength = 1350.0,
                            double steps_size = 10.,
                            int angleMethod = 1,
                            int nsteps = 6);

    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);
    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);
    double GetMuMultiScatterLLHD3(art::Ptr<recob::Track> const& trk, bool dir);
    TVector3 GetMultiScatterStartingPoint(art::Ptr<recob::Track> const& trk);

  private:
    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);
    struct Segments {
      std::vector<double> x, nx;
      std::vector<double> y, ny;
      std::vector<double> z, nz;
      std::vector<double> L;
      std::vector<bool> nvalid;
    };

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

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

    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;

    double my_g(double xx, double Q, double s) const;

    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;

    double seg_stop{-1.};
    int n_seg{};

    double x_seg[50000];
    double y_seg[50000];
    double z_seg[50000];

    double find_angle(double vz, double vy) const;

    std::vector<double> steps;

    double minLength;
    double maxLength;
    double steps_size;
    double rad_length{14.0};

    enum ScatterAngleMethods {
      kAnglezx = 1,   
      kAnglezy,       
      kAngleCombined, 
    };

    ScatterAngleMethods fMCSAngleMethod;
    int n_steps;

    // The following are objects that are created but not drawn or
    // saved.  This class should consider accepting a "debug"
    // parameter where if it is specified, then the graphs will be
    // created; otherwise, their creation is unnecessary and impedes
    // efficiency.
    //
    // N.B. TPolyLine3D objects are owned by ROOT, and we thus refer
    // to them by pointer.  It is important that 'delete' is not
    // called on the TPolyLine3D pointers during destruction of a
    // TrackMomentumCalculator object.
    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{};
  };

} // namespace trkf

#endif // TrackMomentumCalculator_H