ParticleIDAlg.h

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  ParticleIDAlg.h
/// \brief   Calculate dE/dx and log likelihoods for different particle hypotheses.
///     This information will be of use for identification of any particle that
///     uses a Shower as a starting point.
///

///
///     Upon instantiation the histograms are loaded.  Ideally one would reuse an instance of this
///     class to save reloading of the histograms.
///
/// \author  smith@physics.umn.edu
////////////////////////////////////////////////////////////////////////


#ifndef ParticleIDAlg_h
#define ParticleIDAlg_h
#include <vector>
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Provenance/EventID.h"
#include "Geometry/Geometry.h"
#include "ShowerLID/NuEEnergyAlg.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "ShowerLID/DedxDistribution.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/WeightedHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "RecoBase/Shower.h"
#include "Utilities/func/MathUtil.h"
#include "TGeoManager.h"
#include "TGeoNode.h"



namespace rb {
  class Cluster;
  class Track;
  class CellHit;
  struct WeightedHit;
}


namespace slid {

  class ParticleIDAlg {
  public:
    explicit ParticleIDAlg(fhicl::ParameterSet const& pset, 
                           novadaq::cnv::DetId detId,
                           NuEEnergyAlg* NuEEnergyAlgIn,
                           bool loadLibs = true);

    ~ParticleIDAlg();

    /// \brief Reinitialize  parameters read from FCL file
    void reconfigure(const fhicl::ParameterSet & pset, 
                     novadaq::cnv::DetId detId);

    /// \brief  Calculate transverse dE/dx log-likelihood for specific particle hypothesis
    double DedxTransLL(const slid::DedxParticleType partHyp,
                       const rb::Shower* vShower,
                       const std::vector< const rb::Shower* > showercol,
                       const art::Event& evt);


    /// \brief  Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
    double DedxLongLL(const slid::DedxParticleType partHyp, 
                      const rb::Shower* vShower,
                      const std::vector< const rb::Shower* > showercol,
                      const art::Event& evt);

    /// \brief  Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
    /// assuming that the particle direction if opposite to the reconstructed direction
    double DedxInverseLongLL(const slid::DedxParticleType partHyp, 
                             const rb::Shower* vShower,
                             const std::vector< const rb::Shower* > showercol,
                             const art::Event& evt);

    /// \brief  Set rb::Prong to be analyzed.  This must be set before any calculations can be done.
    void SetShower(const rb::Shower* vShower, 
                   const std::vector< const rb::Shower* > showercol,
                   const art::Event& evt);



    ///\brief Calculate the gap between the vertex and the start of the shower
    double GetGapVertexToShowerStart(const rb::Shower* vShower, TVector3 evtVtx,
                                     const art::Event& evt);

    ///\brief Calculate the mass of candidate and any other shower in the slice that is closest to the pi0 mass.
    double Pi0Mass(const std::vector< const rb::Shower* > showerCollection,
                   unsigned int iShower, const art::Event& evt);

    ///\brief Index in showerCollection of the shower for which the reconstructed 
    /// invariant mass was closest to pi0 mass for the shower at index iShower.
    int Pi0SecondPhoton(const std::vector< const rb::Shower* > showerCollection,
                           unsigned int iShower, const art::Event& evt);
    
    ///\brief Calculate the shower radius
    double Radius(const std::vector<const rb::Shower*> showerCollection,
                                unsigned int iShower, const art::Event& evt);

    ///\brief return longitudinal dedx for a specified plane index
    double PlaneLongDedx(unsigned int pIdx);

    ///\brief return transverse dedx for specified transverse plane
    double PlaneTransDedx(unsigned int tpIdx);

    ///\brief return transverse dedx for specified transverse plane
    double CellPlaneDedx(int tpIdx, unsigned int pIdx);

    ///\brief return shower radius for a plane index
    double PlaneRadius(unsigned int pIdx);

    ///\brief Return point of intersection between shower core and plane
    TVector3 PlaneHitXYZ(unsigned int pIdx);

    ///\brief Return point of intersection cell between shower core and plane
    double PlaneHitCell(unsigned int pIdx);
    ///\brief Return the gloabel plane index for an index in the shower coordinate
    unsigned int PlaneToGlobalPlane(unsigned int pIdx);

    ///\brief Return the energy weighted mean shower positon for the plane in cm
    double PlaneCentroid(unsigned int pIdx);

    ///\brief Return the asymmetry of the shower
    double ShowerAsym() const {return fAsym;}

    ///\brief Return the inertia calculation
    TVector3 InertiaXYZ() const {return fIxyz;}

    /// \brief Get distance of closest approach between shower and vertex
    double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop);

    /// \brief Number of MIP planes within the first 20 planes of the start of the shower
    int GetMIPPlanes(const rb::Shower* shower,
                     const std::vector< const rb::Shower* > showercol,
                     const art::Event& evt);

    /// Map of the longitudinal ll by paricle type
    std::map<int, float> fPartLongLL;

    /// Map of the transverse ll by particle type
    std::map<int, float> fPartTransLL;

  private:
    
    /// \brief Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector.
    double CalcPlaneLongitudinalDedx(unsigned int plane);

    /// \brief Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index.
    double CalcPlaneTransverseDedx(unsigned int iTransversePlane);


    /// \brief Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number.
    double CalcCellPlaneTransverseDedx(int iTransversePlane, unsigned int iPlane);

    /// \brief  Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane
    /// for shower already read into object by SetShower
    std::map< int, std::map< int, int > > CalcPlaneHits(const rb::Shower* vShower );
    
    /// \brief  Return a map of mean cell posisiton for each plane
    /// for shower already read into object by SetShower
    std::map< int, double > CalcPlaneCentroid( const rb::Shower* vShower);

    /// \brief Sum the energy in a longitudinal plane
    double CalcEnergyInLongitudinalPlane(int iPlane);

    /// \brief Determine detector XY region
    void CalcDetectorXYRegionIndex();

    /// \brief  Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesis
    double CalcPlaneLongDedxProb(DedxParticleType particleType, int igev, unsigned int iPlane, unsigned int plane);

    /// \brief Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis
    double CalcPlaneTransDedxProb(DedxParticleType particleType, int igev, unsigned int iTransversePlane);

    /// \brief From Jianming's code, still working out changes.
    double CalcInterPlaneDedxProb(DedxParticleType type, unsigned int iLongPlane, unsigned int plane);

    /// \brief From Jianming's code, still working out changes.
    double CalcInterCellTransDedxProb(const DedxParticleType type, unsigned int iTransPlane);

    /// \brief Load all of the histograms from files into arrays of TH1D* objects.  The details of handling the
    ///  histograms are taken care of within object of the helper class slid::DedxDistribution.
    ///
    bool LoadDedxHistogramFiles(const fhicl::ParameterSet & pset, novadaq::cnv::DetId geom);

    /// \brief Calculate path length of shower through this plane.
    double CalcTrkHitPath(unsigned int plane);

    /// \brief  Calculate energy bin indices and energy values from the histograms
    ///  corresponding energy of shower to be used to interpolate between bins in energy
    void CalculateEnergyBinIndex();

    /// \brief Calculate the asymmetry and inertia of the shower 
    void CalcAsymIneria();

    /// \brief Helper function for interia and asymmetry calculation
    TVector3 GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2);

    /// data/MC flag
    //bool isData;

    /// event ID
    art::EventID eventID;
 
    /// Closest mass to pi0 mass formed by this shower and any other shower in the slice.
    double fClosestPi0Mass;
                
    /// dE/dx distributions (histogram files and data) for each detector type and particle category
    std::vector<slid::DedxDistribution> fDedxDistributions;

    /// rb::Shower object for which particle ID is to be calculated.  This is cached.
    const rb::Shower* fShower;
    std::vector< const rb::Shower* > fShowerCol;

    /// Shower energy
    double fShowerEnergyGev;

    /// geometry
    art::ServiceHandle<geo::Geometry> fGeom;

    /// Nue energy algorithms
    NuEEnergyAlg* fNuEEnergyAlg;


    /// Map of longitudinal dedx by plane. The first element is
    /// plane number and the second is the dedx for that plane.
    /// This is used for lazy caching of the dedx plane values
    /// in order to save recalculating for the same shower.
    std::map< int, double> fDedxLongByPlane;

    /// Map of transverse dedx by plane.  The first element is
    /// plane number and the second is the dedx for that plane.
    /// This is used for lazy caching of the dedx plane values
    /// in order to save recalculating for the same shower.
    std::map< int, double> fDedxTransByPlane;


    /// Map where the first element is
    /// plane number and the second is a map with the key being the
    /// cell number in the plane and the second value is the global
    /// cell index in a shower 
    std::map< int, std::map< int, int > > fPlaneHits;

    /// Map where the key is the index of the plane in the shower and the second is the global plane
    /// This map saves looping in some of the utility functions
    std::map<unsigned int, unsigned int> fMapIdxToGlobalPlane;

    /// Map where the first index is 
    /// plane number and the second is the energy weighted center position for that plane
    std::map< int, double > fPlaneCentroid;

    /// Map of transverse energy for each transverse plane
    std::map<int, double > fTransEnergy;

    /// Map of path length for each transverse plane
    std::map<int, double > fTransPathLength;

    /// Index of the other photon shower used for pi0 mass calculation 
    /// for a given shower
    int fPi0Index;
    
    /// Detector sub-region bin.   Calculated once and cached for a given shower/event combo.
    int iReg;

    /// Shower start position
    TVector3 pos;

    /// Shower adjacent energy bins.   Calculated once and cached for a given shower/event combo.
    int igev0;
    int igev1;

    /// Energy values of energy bins adjacent to shower.   Calculated once and cached for a given shower/event combo.
    double e0;
    double e1;

    /// shower moment of ineria
    TVector3 fIxyz;
    /// shower asymmetry
    double fAsym;
    
    /// Upper bound of MIP Range
    double fMipRangeHigh;

    /// Lower bound of MIP Range
    double fMipRangeLow;

    /// Detector ID
    novadaq::cnv::DetId fDetId;

  };
}
#endif // ParticleIDAlg_h

////////////////////////////////////////////////////////////////////////