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 <array>
#include <memory>  // std::unique_ptr<>
#include <ostream> // std::endl
#include <vector>

// framework libraries
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/fwd.h"

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

namespace hit {

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

    explicit CCHitFinderAlg(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;
    static constexpr float Sqrt2Pi = 2.5066;
    static constexpr float SqrtPi = 1.7725;

    bool fUseChannelFilter;

    geo::WireReadoutGeom const* wireReadoutGeom{
      &art::ServiceHandle<geo::WireReadout const>()->Get()};

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

    struct HitChannelInfo_t {
      recob::Wire const* wire;
      geo::WireID wireID;
      geo::SigType_t sigType;
    };

    // 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* signl = nullptr,
                   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; 

    GausFitCache FitCache; 

    struct FitStats_t {
      unsigned int FastFits;                   
      std::vector<unsigned int> MultiGausFits; 

      void Reset(unsigned int nGaus);
      void AddMultiGaus(unsigned int nGaus);
      void AddFast() { ++FastFits; }
    };

    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