CCHitFinderAlg.h

Go to the documentation of this file.

// ClusterCrawlerAlg.h
//
// ClusterCrawlerAlg class
//
// Bruce Baller
//
#ifndef CCHITFINDERALG_H
#define CCHITFINDERALG_H

// C/C++ standard libraries
#include <string>
#include <vector>
#include <memory> // std::unique_ptr<>
#include <ostream> // std::endl

// framework libraries
#include "fhiclcpp/ParameterSet.h" 
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h" 
#include "canvas/Persistency/Common/Ptr.h" 
#include "canvas/Persistency/Common/PtrVector.h" 

// LArSoft libraries
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larreco/RecoAlg/GausFitCache.h"


namespace hit {
  
  class CCHitFinderAlg {
  
  public:
    
    std::vector<recob::Hit> allhits;

    CCHitFinderAlg(fhicl::ParameterSet const& pset);
    virtual ~CCHitFinderAlg() = default;

    virtual void reconfigure(fhicl::ParameterSet const& pset);

    void RunCCHitFinder(std::vector<recob::Wire> const& Wires);
    
    std::vector<recob::Hit>&& YieldHits() { return std::move(allhits); }
    
    template <typename Stream>
    void PrintStats(Stream& out) const;
    
  private:
    
    std::vector<float> fMinPeak;
    std::vector<float> fMinRMS;
    unsigned short fMaxBumps; // make a crude hit if > MaxBumps are found in the RAT
    unsigned short fMaxXtraHits; // max num of hits in Region Above Threshold
    float fChiSplit;      
  //  float ChgNorm;     // Area norm for the wire we are working on

    std::vector<float> fChiNorms;
    std::vector<float> fTimeOffsets;
    std::vector<float> fChgNorms;

    raw::ChannelID_t theChannel;
    unsigned short theWireNum;
    unsigned short thePlane;

    float chinorm;
  //  float timeoff;
    static constexpr float Sqrt2Pi = 2.5066;
    static constexpr float SqrtPi  = 1.7725;
    
    bool fUseChannelFilter;

//    bool prt;
    
    art::ServiceHandle<geo::Geometry> geom;

    // fit n Gaussians possibly with bounds setting (parmin, parmax)
    void FitNG(unsigned short nGaus, unsigned short npt, float *ticks,
       float *signl);
    // parameters, errors, lower limit, upper limits for FitNG
    std::vector<double> par;
    std::vector<double> parerr;
    std::vector<double> parmin;
    std::vector<double> parmax;
    float chidof;
    int dof;
    std::vector<unsigned short> bumps;
    
    class HitChannelInfo_t {
        public:
      recob::Wire const* wire;
      geo::WireID wireID;
      geo::SigType_t sigType;
      
      HitChannelInfo_t
        (recob::Wire const* w, geo::WireID wid, geo::Geometry const& geom);
    }; // HitChannelInfo_t
    
    // make a cruddy hit if fitting fails
    void MakeCrudeHit(unsigned short npt, float *ticks, float *signl);
    // store the hits
    void StoreHits(unsigned short TStart, unsigned short npt, 
      HitChannelInfo_t info, float adcsum
      );

    // study hit finding and fitting
    bool fStudyHits;
    std::vector< short > fUWireRange, fUTickRange;
    std::vector< short > fVWireRange, fVTickRange;
    std::vector< short > fWWireRange, fWTickRange;
    void StudyHits(unsigned short flag, unsigned short npt = 0,
      float *ticks = 0, float *signl = 0, unsigned short tstart = 0);
    std::vector<int> bumpCnt;
    std::vector<int> RATCnt;
    std::vector<float> bumpChi;
    std::vector<float> bumpRMS;
    std::vector<int> hitCnt;
    std::vector<float> hitRMS;
    // use to determine the slope of protons
    std::vector<float> loWire;
    std::vector<float> loTime;
    std::vector<float> hiWire;
    std::vector<float> hiTime;
    bool SelRAT; // set true if a Region Above Threshold should be studied
    
    bool fUseFastFit; 
    
    std::unique_ptr<GausFitCache> FitCache; 
    
    
    typedef struct {
      unsigned int FastFits; 
      std::vector<unsigned int> MultiGausFits; 
      
      void Reset(unsigned int nGaus);
      
      void AddMultiGaus(unsigned int nGaus);
      
      void AddFast() { ++FastFits; }
      
    } FitStats_t;
    
    FitStats_t FinalFitStats; 
    FitStats_t TriedFitStats; 
    
    static bool FastGaussianFit(
      unsigned short npt, float const*ticks, float const*signl,
      std::array<double, 3>& params,
      std::array<double, 3>& paramerrors,
      float& chidof
      );
    
    static constexpr unsigned int MaxGaussians = 20;
    
  }; // class CCHitFinderAlg
  
} // namespace hit


//==============================================================================
//===  Template implementation
//===
template <typename Stream>
void hit::CCHitFinderAlg::PrintStats(Stream& out) const {
  
  out << "CCHitFinderAlg fit statistics:";
  if (fUseFastFit) {
    out << "\n  fast 1-Gaussian fits: " << FinalFitStats.FastFits
      << " succeeded (" << TriedFitStats.FastFits << " tried)";
  }
  else
    out << "\n  fast 1-Gaussian fits: disabled";
  
  for (unsigned int nGaus = 1; nGaus < MaxGaussians; ++nGaus) {
    if (TriedFitStats.MultiGausFits[nGaus-1] == 0) continue;
    out << "\n  " << nGaus << "-Gaussian fits: "
      << FinalFitStats.MultiGausFits[nGaus-1]
      << " accepted (" << TriedFitStats.MultiGausFits[nGaus-1] << " tried)";
  } // for nGaus
  if (TriedFitStats.MultiGausFits.back() > 0) {
    out << "\n  " << FinalFitStats.MultiGausFits.size()
      << "-Gaussian fits or higher: " << FinalFitStats.MultiGausFits.back()
      << " accepted (" << TriedFitStats.MultiGausFits.back() << " tried)";
  }
  out << std::endl;
  
} // CCHitFinderAlg::FitStats_t::Print()




#endif // ifndef CCHITFINDERALG_H