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

    double GetTrackMomentum(double trkrange, int pdg) const;
    double GetMomentumMultiScatterChi2(art::Ptr<recob::Track> const& trk,
                                       const bool checkValidPoints = false,
                                       const int maxMomentum_MeV = 7500);
    double GetMomentumMultiScatterLLHD(art::Ptr<recob::Track> const& trk,
                                       const bool checkValidPoints = false,
                                       const int maxMomentum_MeV = 7500,
                                       const int MomentumStep_MeV = 10,
                                       const int max_resolution = 0);
    double GetMuMultiScatterLLHD3(art::Ptr<recob::Track> const& trk, bool dir);
    TVector3 GetMultiScatterStartingPoint(art::Ptr<recob::Track> const& trk);

  private:
    bool plotRecoTracks_(std::vector<float> const& xxx,
                         std::vector<float> const& yyy,
                         std::vector<float> const& zzz);

    void compute_max_fluctuation_vector(const std::vector<float> segx,
                                        const std::vector<float> segy,
                                        const std::vector<float> segz,
                                        std::vector<float>& segnx,
                                        std::vector<float>& segny,
                                        std::vector<float>& segnz,
                                        std::vector<float>& vx,
                                        std::vector<float>& vy,
                                        std::vector<float>& vz);
    struct Segments {
      std::vector<float> x, nx;
      std::vector<float> y, ny;
      std::vector<float> z, nz;
      std::vector<float> L;
    };

    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;

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

    float x_seg[50000];
    float y_seg[50000];
    float z_seg[50000];

    double find_angle(double vz, double vy) const;

    int n_steps{6};
    std::vector<float> steps;

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

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