MEFinder_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//  Michel Electron Finder package
//  Dan Pershey
//  pershey@caltech.edu
//
//  Jul 9 2015.
//  MEFinder is a Michel E tagging package designed to meet all of NOvA's 
//  Michel reconstruction needs.  It is a revamp of MichelEFilter and can
//  satisfy both our calibration and current analysis requirements for Michels.
//
///////////////////////////////////////////////////////////////////////////////

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Assns.h"
#include "cetlib/search_path.h"


// NOvASoft includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "MEFinder/MEClusters.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"

#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TMath.h"

namespace me
{

  // These guys come up a lot...
  typedef art::Ptr<rb::CellHit>     PtrCHit;
  typedef std::vector<PtrCHit>      VecCHits;
  typedef std::vector<me::VecCHits> VecVecCHits;


  // Holds hists for computing MID, the Michel E LL Identifier
  struct LikeHists
  {
    TH2F* fMEDeltaT_Dist;
    TH2F* fMECalE_NCells;

    TH2F* fHFDeltaT_Dist;
    TH2F* fHFCalE_NCells;
  };

  // Collection of relevant variables to ME cluster
  struct MEVars
  {
    double distToSlc;
    art::Ptr<rb::Cluster> parentSlc;
    // Only valid for TrkMEs
    art::Ptr<rb::Track>   parentTrk;
  };

  // Holds fhicl parameters
  struct MEFinderParams
  {
    template<class T> using Atom = fhicl::Atom<T>;
    template<class T> using Table = fhicl::Table<T>;
    using Comment = fhicl::Comment;
    using Name = fhicl::Name;


    // Grab necessary input from other reco modules
    Atom<std::string> slicerInput
    {
      Name("SlicerInput"),
      Comment("Input for physics slices")
    };
    Atom<std::string> trackInput
    {
      Name("TrackInput"),
      Comment("Input for tracks with which to associate MEs")
    };

    // Params for deciding if noise hit is good enough to physics slice
    Atom<double> maxDist
    {
      Name("MaxDist"),
      Comment("Max dist from phsyics slice to selected noise hit")
    };
    Atom<double> minHitADC
    {
      Name("MinHitADC"),
      Comment("Min ADC for a hit to be considered - needed for FD speed ")
    };

    // Should we calculate MID values for this job
    Atom<bool> calibAvailable
    {
      Name("CalibAvailable"),
      Comment("Are there calibration constants available for getting CalE")
    };

    // Preselection parameters
    Atom<double> maxClosestApproach
    {
      Name("MaxClosestApproach"),
      Comment("ME clust must be at least this close to a hit in phys slice")
    };
    Atom<double> trkClosestApproach
    {
      Name("TrkClosestApproach"),
      Comment("Dist cut to be considered TrkME")
    };
    Atom<unsigned int> maxNCells
    {
      Name("MaxNCells"),
      Comment("Presel Cut")
    };
    Atom<unsigned int> minNCells
    {
      Name("MinNCells"),
      Comment("Presel Cut")
    };
    Atom<double> maxDeltaT
    {
      Name("MaxDeltaT"),
      Comment("Presel Cut")
    };
    Atom<double> minDeltaT
    {
      Name("MinDeltaT"),
      Comment("Presel Cut")
    };
    Atom<double> maxCalE
    {
      Name("MaxCalE"),
      Comment("Presel Cut")
    };
    Atom<double> minCalE
    {
      Name("MinCalE"),
      Comment("Presel Cut")
    };
    Atom<double> minDeltaTReTrigger
    {
      Name("MinDeltaTReTrigger"),
      Comment("See DocDB 12295, apd dead time not well modeled in MC")
    };

    // DBSCAN parameters
    Atom<double> eps
    {
      Name("EpsForClustering"),
      Comment("Dist to search out to in the DBSCAN algorithm")
    };
    Atom<double> planeScale
    {
      Name("PlaneScale"),
      Comment("For distances in DBSCAN")
    };
    Atom<double> cellScale
    {
      Name("CellScale"),
      Comment("For distances in DBSCAN")
    };
    Atom<double> timeScale
    {
      Name("TimeScale"),
      Comment("For distances in DBSCAN")
    };

    Atom<std::string> llFileName
    {
      Name("LLFileName"),
      Comment("Contains hist files for calculating Delta LL ME PID, MID")
    };

    Atom<bool> matchMEtoSlcByDist
    {
      Name("MatchMEtoSlcByDist"),
      Comment("Match MEs to Slices by distance only")
    };

  };

  class MEFinder : public art::EDProducer
  {
  public:

    // Allows 'nova --print-description' to work
    using Parameters = art::EDProducer::Table<MEFinderParams>;

    explicit MEFinder(const Parameters& params);
    virtual  ~MEFinder();

    void produce(art::Event& evt);
    void beginJob();

  protected:

    // Following functions work to making a set of ME clusters of rb::CellHits

    // Sorts through noise and small physics hits to find candidate Michel hits
    void GetPhysNoiseSlcs(art::Handle<std::vector<rb::Cluster> > slcs,
                          std::vector<rb::Cluster>               &noiseSlcs,
                          std::vector<rb::Cluster>               &physSlcs);
    me::VecCHits MECandHits(std::vector<rb::Cluster> noiseSlcs,
                            std::vector<rb::Cluster> physSlcs);
    void FindNoiseHitRange(me::VecCHits hits, double physTNS,
                           int& firstIdx, int& latsIdx);
    // Adds noise hit to list of Michel hits if it is near a physics slice
    void AddMEHit(me::VecCHits &meHits,
                  me::PtrCHit  noiseCell,
                  rb::Cluster  physSlc);
    // These two find min distance between a hit and a physics slice
    double MinHitSlcDist(me::PtrCHit hit, rb::Cluster slc);
    double HitToHitDistance(const PtrCHit &hit1,
                            const PtrCHit &hit2);
    // Sorts these ME hits into clusters of rb::CellHits
    me::VecVecCHits ClusterMEHits(VecCHits);
    // Holds the distance function between points for clustering
    double MEDistanceMetric(me::PtrCHit hit1, me::PtrCHit hit2);


    // Following functions make art SlcME/TrkME products and compute variables

    me::SlcME MakeSlcME(me::VecCHits michel,
                        art::Handle<std::vector<rb::Cluster> > slcs,
                        int& parentSlcIdx);
    me::TrkME MakeTrkME(me::SlcME michel,
                        art::Handle<std::vector<rb::Cluster> > slcs,
                        art::Ptr<rb::Track> parentTrack);
    void BestSlcMatch(me::SlcME michel,
                      me::MEVars& vars,
                      art::Handle<std::vector<rb::Cluster> > slcs,
                      int& parentSlcIdx);
    bool PassesMEPresel(me::SlcME michel);
    double GetMID(double deltaT, double distToSlc, double calE, int nCells);


    // Asks if ME is a TrkME, and fills out relevant extra variables

    bool IsTrkME(me::SlcME michel,
                 const std::vector< art::Ptr<rb::Track> > trks,
                 art::Ptr<rb::Track>& parentTrk);
    double TrkMEDist(me::SlcME slcME, TVector3 trkEnd);
    std::vector<double> MEMeanXYZ(me::VecCHits me, bool&, bool&);


    // Need a hold of the geometry to calculated cell-cell distances
    art::ServiceHandle<geo::Geometry> fgeom;

    MEFinderParams fParams;

    LikeHists      fLikeHists;
  };
}

namespace me
{
  MEFinder::MEFinder(const Parameters& params)
    : fParams(params())
  {
    produces<std::vector<me::SlcME> >();
    produces<art::Assns<me::SlcME, rb::Cluster> >();
    produces<std::vector<me::TrkME> >();
    produces<art::Assns<me::TrkME, rb::Track> >();
    produces<art::Assns<me::TrkME, rb::Cluster> >();
  }

  MEFinder::~MEFinder() { }

  // beginJob pulls out hist files for MID
  void MEFinder::beginJob()
  {
    cet::search_path sp("FW_SEARCH_PATH");
    std::string LikeFileName;
    sp.find_file(fParams.llFileName(), LikeFileName);
    TFile* like = TFile::Open(LikeFileName.c_str());
    like->cd();

    fLikeHists.fMEDeltaT_Dist  = (TH2F*)gDirectory->Get("MEDeltaT_Dist");
    fLikeHists.fMECalE_NCells  = (TH2F*)gDirectory->Get("MECalE_NCells");

    fLikeHists.fHFDeltaT_Dist  = (TH2F*)gDirectory->Get("HFDeltaT_Dist");
    fLikeHists.fHFCalE_NCells  = (TH2F*)gDirectory->Get("HFCalE_NCells");
  }

  // Main produce function.  Finds ME's and puts them in ART
  void MEFinder::produce(art::Event& evt)
  {
    // Grab ahold of slc and trk info from the art event, set up Assns
    art::Handle<std::vector<rb::Cluster> > slcs;
    evt.getByLabel(fParams.slicerInput(), slcs);
    art::FindManyP<rb::Track> trkInfo(slcs, evt, fParams.trackInput());

    std::unique_ptr< std::vector<me::SlcME> >
      slcMEs(new std::vector<me::SlcME>);
    std::unique_ptr< art::Assns<me::SlcME, rb::Cluster> >
      slcAssn(new art::Assns<me::SlcME, rb::Cluster>);

    std::unique_ptr< std::vector<me::TrkME> >
      trkMEs(new std::vector<me::TrkME>);
    std::unique_ptr< art::Assns<me::TrkME, rb::Track> >
      trkAssn(new art::Assns<me::TrkME, rb::Track>);
    std::unique_ptr< art::Assns<me::TrkME, rb::Cluster> >
      trkMEToSlcAssn(new art::Assns<me::TrkME, rb::Cluster>);

    // Vector of all CellHits close to any parent slice in time and space
    std::vector<rb::Cluster> noiseSlcs, physSlcs;
    GetPhysNoiseSlcs(slcs, noiseSlcs, physSlcs);
    me::VecCHits meHits = MECandHits(noiseSlcs, physSlcs);

    // Then sort this vec of CellHits into ME clusters
    me::VecVecCHits meClusters = ClusterMEHits(meHits);

    for (unsigned int i = 0; i < meClusters.size(); i++){
      art::Ptr<rb::Track> parentTrk;
      int parentSlcIdx;
      me::SlcME slcME =  MakeSlcME(meClusters[i], slcs, parentSlcIdx);
      if (!PassesMEPresel(slcME)) continue;
      const std::vector< art::Ptr<rb::Track> > trks = trkInfo.at(parentSlcIdx);
      if (IsTrkME(slcME, trks, parentTrk)){
        me::TrkME trkME = MakeTrkME(slcME, slcs, parentTrk);
        trkMEs->push_back(trkME);
        util::CreateAssn(*this,evt,*trkMEs,trkME.ParentTrk(),*(trkAssn.get()));
        util::CreateAssn(*this,evt,*trkMEs,trkME.ParentSlc(),
                         *(trkMEToSlcAssn.get()));
      }
      else {
        slcMEs->push_back(slcME);
        util::CreateAssn(*this,evt,*slcMEs,slcME.ParentSlc(),*(slcAssn.get()));
      }
    }
    evt.put(std::move(slcMEs));
    evt.put(std::move(slcAssn));
    evt.put(std::move(trkMEs));
    evt.put(std::move(trkAssn));
    evt.put(std::move(trkMEToSlcAssn));
  }

  // Outputs a list of noise cells close enough to physics to be a ME hit
  void MEFinder::GetPhysNoiseSlcs(art::Handle< std::vector<rb::Cluster> > slcs,
                                  std::vector<rb::Cluster> &noiseSlcs,
                                  std::vector<rb::Cluster> &physSlcs)
  {
    // Label each slice either noise or physics.  A noise slice can be a slice 
    // flagged as noise, or slices with only a few hits  ( <= 12 )
    for (unsigned int i = 0; i < slcs->size(); i++){
      rb::Cluster slice = slcs->at(i);
      if (slice.IsNoise())
        noiseSlcs.push_back(slice);
      else if (slice.NCell() <= fParams.maxNCells())
        noiseSlcs.push_back(slice);
      else
        physSlcs.push_back(slice);
    }
  }

  me::VecCHits MEFinder::MECandHits(std::vector<rb::Cluster> noiseSlcs,
                                    std::vector<rb::Cluster> physSlcs){
    // Loop over all noise hits, then ask each physics slice if the noise
    // hit is close in time / space, then add the hit to the list of candidate
    // hits appropriately.
    me::VecCHits meHits;

    for (unsigned int nIdx = 0; nIdx < noiseSlcs.size(); nIdx++){
      me::VecCHits noiseHits;
      for (unsigned int cIdx = 0; cIdx < noiseSlcs[nIdx].NCell(); cIdx++)
        if (noiseSlcs[nIdx].Cell(cIdx)->ADC() > fParams.minHitADC())
          noiseHits.push_back(noiseSlcs[nIdx].Cell(cIdx));
      rb::SortByTime(noiseHits);
      for (unsigned int pIdx = 0; pIdx < physSlcs.size(); pIdx++){
        int fIdx = 0; int lIdx = 0;
        FindNoiseHitRange(noiseHits, physSlcs[pIdx].MeanTNS(), fIdx, lIdx);
        for (int cIdx = fIdx; cIdx <= lIdx; cIdx++){
          AddMEHit(meHits, noiseHits[cIdx], physSlcs[pIdx]);
        }
      }
    }
    return meHits;
  }
  
  // Calculates the first / last noise hit with a delta T in the allowed 
  // region for MEs for a given slice
  void MEFinder::FindNoiseHitRange(me::VecCHits nHits, double physTNS,
                                   int& fIdx, int& lIdx)
  {
    fIdx = 0; lIdx = 0;
    double fMETime = physTNS+fParams.minDeltaT();
    double lMETime  = physTNS+fParams.maxDeltaT();
    int nIter = log2(nHits.size())+1;
    for (int i = 1; i <= nIter; i++){
      int fmidpoint = fIdx + (nHits.size()-1)/pow(2,i);
      fIdx = (nHits[fmidpoint]->TNS() > fMETime) ? fIdx : fmidpoint;
      int lmidpoint = lIdx + nHits.size()/pow(2,i);
      lIdx = (nHits[lmidpoint]->TNS() > lMETime) ? lIdx : lmidpoint;
    }
    fIdx++;
  }

  // Adds current noiseCell to list of ME hits if close to physSlc
  void MEFinder::AddMEHit(me::VecCHits &meHits,
                          me::PtrCHit  noiseCell,
                          rb::Cluster  physSlc){
    double minDist = MinHitSlcDist(noiseCell, physSlc);
    double deltaT  = noiseCell->TNS() - physSlc.MeanTNS();
    if (noiseCell->ADC() < fParams.minHitADC())
      return;
    // Treat Retriggers differently because of APD dead-time (DocDB 12295)
    // Additionally require deltaT > 0 for tagging retrigger hits to
    // accomodate negative search time for measuring uncorrelated background 
    double curMinDeltaT = ( minDist == 0 && deltaT > 0 ) ?
      fParams.minDeltaTReTrigger() : fParams.minDeltaT();
    if (deltaT  > curMinDeltaT &&
        deltaT  < fParams.maxDeltaT() &&
        minDist < fParams.maxDist()){
      bool addHit = true;
      for (int i = 0; i < (int)meHits.size(); i++){
        if (noiseCell->Plane() == meHits[i]->Plane() &&
            noiseCell->Cell()  == meHits[i]->Cell()  &&
            noiseCell->ADC()   == meHits[i]->ADC())
          addHit = false;
      }
      if (addHit) meHits.push_back(noiseCell);
    }
    return;
  }

  // Sorts these hits into clusters
  me::VecVecCHits MEFinder::ClusterMEHits(me::VecCHits meHits){
    me::VecVecCHits meClusters;
    if (meHits.size() == 0)
      return meClusters;

    // This will keep track of wheter hit is already clustered
    std::vector<bool> isUsed(meHits.size(), false);

    for (unsigned int i = 0; i < meHits.size(); i++){
      if (isUsed[i]) continue;
      me::VecCHits curCluster;
      curCluster.push_back(meHits[i]);
      isUsed[i] = true;
      for (unsigned int j = i; j < meHits.size(); j++){
        if (isUsed[j]) continue;
        if (MEDistanceMetric(meHits[i], meHits[j]) < fParams.eps()){
          isUsed[j] = true;
          curCluster.push_back(meHits[j]);
          // Check if any of the hits we just passed over are close to jth idx
          for (unsigned int k = i; k < j; k++){
            if (isUsed[k]) continue;
            if (MEDistanceMetric(meHits[j], meHits[k]) < fParams.eps()){
              isUsed[k] = true;
              curCluster.push_back(meHits[k]);
            }
          }
        }
      }
      meClusters.push_back(curCluster);
    }
    return meClusters;
  }

  double MEFinder::MEDistanceMetric(me::PtrCHit hit1, me::PtrCHit hit2)
  {
    double deltaC = (hit1->View() == hit2->View()) ?
      double(hit1->Cell()  - hit2->Cell())  / fParams.cellScale() : 0;
    double deltaP = double(hit1->Plane() - hit2->Plane());
    deltaP /= fParams.planeScale();
    double deltaT = double(hit1->TNS()   - hit2->TNS());
    deltaT /= fParams.timeScale();
    return sqrt(deltaC*deltaC + deltaP*deltaP + deltaT*deltaT);
  }

  // Min dist between a ME hit and a physics slice
  double MEFinder::MinHitSlcDist(PtrCHit hit, rb::Cluster slc){
    double minDist = 10000;
    for (unsigned int i = 0; i < slc.NCell(); i++){
      double curDist = HitToHitDistance(hit, slc.Cell(i));
      if (curDist < minDist){
        minDist = curDist;
      }
    }
    return minDist;
  }

  // Calculates dist between two cells
  double MEFinder::HitToHitDistance(const PtrCHit &hit1, const PtrCHit &hit2){
    // Only calculate if both hits are in the same view
    if (hit1->View() != hit2->View()) return 100000.;
    else{
      double xyz1[3], xyz2[3];

      unsigned short plane1_idx = hit1->Plane();
      unsigned short cell1_idx = hit1->Cell();
      const geo::PlaneGeo* plane1 = fgeom->Plane(plane1_idx);
      const geo::CellGeo* cell1 = plane1->Cell(cell1_idx);
      cell1->GetCenter(xyz1);

      unsigned short plane2_idx = hit2->Plane();
      unsigned short cell2_idx = hit2->Cell();
      const geo::PlaneGeo* plane2 = fgeom->Plane(plane2_idx);
      const geo::CellGeo* cell2 = plane2->Cell(cell2_idx);
      cell2->GetCenter(xyz2);

      // Agnostic to whether hits are in X or Y views
      return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
    }
  }


  me::SlcME MEFinder::MakeSlcME(me::VecCHits michel,
                                art::Handle<std::vector<rb::Cluster> > slcs,
                                int& parentSlcIdx){
    me::SlcME slcME;
    me::MEVars vars;

    for (unsigned int i = 0; i < michel.size(); i++)
      slcME.Add(michel[i]);
    BestSlcMatch(slcME, vars, slcs, parentSlcIdx);
    slcME.SetDistToSlc    (vars.distToSlc);
    slcME.SetParentSlc    (vars.parentSlc);
    slcME.SetDeltaT       (slcME.MeanTNS() - vars.parentSlc->MeanTNS());
    double mid = -10;
    if (fParams.calibAvailable()){
      mid = GetMID(slcME.DeltaT(),
                   slcME.DistToSlc(),
                   slcME.CalorimetricEnergy(),
                   slcME.NCell());
    }
    slcME.SetMID(mid);
    return slcME;
  }

  me::TrkME MEFinder::MakeTrkME(me::SlcME slcME,
                                art::Handle<std::vector<rb::Cluster> > slcs,
                                art::Ptr<rb::Track> parentTrk){
    me::TrkME trkME;
    me::MEVars vars;

    for (unsigned int i = 0; i < slcME.NCell(); i++)
      trkME.Add(slcME.Cell(i));

    trkME.SetDistToSlc    (slcME.DistToSlc());
    trkME.SetParentSlc    (parentTrk);
    trkME.SetDeltaT       (slcME.MeanTNS() - parentTrk->MeanTNS());
    double mid = -10;
    if (fParams.calibAvailable()){
      mid = GetMID(slcME.DeltaT(),
                   slcME.DistToSlc(),
                   slcME.CalorimetricEnergy(),
                   slcME.NCell());
    }
    trkME.SetMID(mid);
    trkME.SetParentTrk(parentTrk);
    return trkME;
  }

  // Takes a ME Cluster, and finds out what physics slice it is most likely
  // associated with using deltaT and distToSlc info
  void MEFinder::BestSlcMatch(me::SlcME slcME,
                              me::MEVars& vars,
                              art::Handle<std::vector<rb::Cluster> > slcs,
                              int& parentSlcIdx){
    // Grab the MID building histogram
    TH2F* hAssn = fLikeHists.fMEDeltaT_Dist;
    // The assnScore is larger for slices that are appropriately near the ME
    double assnScore = -1;
    parentSlcIdx = 0;
    double allSlcMinDist = 1.0e6;
    for (unsigned int i = 0; i < slcs->size(); i++){
      rb::Cluster slc = slcs->at(i);
      if (slc.IsNoise() || slc.NCell() <= fParams.maxNCells())
        continue;
      art::Ptr<rb::Cluster> slcPtr(slcs, i);

      double deltaT = slcME.MeanTNS() - slc.MeanTNS();
      double minDist = 10000;

      for (unsigned int j = 0; j < slcME.NCell(); j++){
        double dist = MinHitSlcDist(slcME.Cell(j), slcs->at(i));
        if (dist < minDist)
          minDist = dist;
      } // Loop over ME hits

      // Match by distance only
      if ( fParams.matchMEtoSlcByDist() )
      {
        if ( minDist < allSlcMinDist )
        {
          parentSlcIdx = i;
          allSlcMinDist = minDist;
          vars.parentSlc     = slcPtr;
          vars.distToSlc     = minDist;
        }
      }

      // Match by distance and time (standard method)
      else
      {
        double curScore = hAssn->GetBinContent(hAssn->FindBin(deltaT, minDist));
        if (curScore > assnScore){
          parentSlcIdx = i;
          assnScore = curScore;
          vars.parentSlc     = slcPtr;
          vars.distToSlc     = minDist;
        }
      }

    } // Loop over slices
  }

  bool MEFinder::PassesMEPresel(me::SlcME slcME){
    if (slcME.DistToSlc() > fParams.maxClosestApproach())
      return false;
    if (slcME.NCell() > fParams.maxNCells())
      return false;
    if (slcME.NCell() < fParams.minNCells())
      return false;
    if (slcME.DeltaT() > fParams.maxDeltaT())
      return false;
    if (slcME.DeltaT() < fParams.minDeltaT())
      return false;
    // Only apply calE presel if a calibration is available
    if (fParams.calibAvailable()){
      if (slcME.CalorimetricEnergy() > fParams.maxCalE())
        return false;
      if (slcME.CalorimetricEnergy() < fParams.minCalE())
        return false;
    }
    return true;
  }

  double MEFinder::GetMID(double DeltaT, double Dist, double CalE, int NCells){
    TH2* meDeltaT_Dist = fLikeHists.fMEDeltaT_Dist;
    TH2* meCalE_NCells = fLikeHists.fMECalE_NCells;

    int meDeltaT_Dist_Bin = meDeltaT_Dist->FindBin(DeltaT, Dist);
    int meCalE_NCells_Bin = meCalE_NCells->FindBin(CalE, NCells);
    double meLike = meDeltaT_Dist->GetBinContent(meDeltaT_Dist_Bin) *
                    meCalE_NCells->GetBinContent(meCalE_NCells_Bin);

    TH2* hfDeltaT_Dist = fLikeHists.fHFDeltaT_Dist;
    TH2* hfCalE_NCells = fLikeHists.fHFCalE_NCells;

    int hfDeltaT_Dist_Bin = hfDeltaT_Dist->FindBin(DeltaT, Dist);
    int hfCalE_NCells_Bin = hfCalE_NCells->FindBin(CalE, NCells);
    double hfLike = hfDeltaT_Dist->GetBinContent(hfDeltaT_Dist_Bin) *
                    hfCalE_NCells->GetBinContent(hfCalE_NCells_Bin);

    if (meLike == 0 && hfLike == 0) return -10;
    if (meLike == 0)
      meLike = 1e-10;
    if (hfLike == 0)
      hfLike = 1e-10;

    return TMath::Log(meLike) - TMath::Log(hfLike);

  }

  bool MEFinder::IsTrkME(me::SlcME slcME,
                         const std::vector< art::Ptr<rb::Track> > trks,
                         art::Ptr<rb::Track>& parentTrk){
    for (unsigned int i = 0; i < trks.size(); i++){
      if (trks[i]->NXCell() == 0 || trks[i]->NYCell() == 0)
        continue;
      double dist = TrkMEDist(slcME, trks[i]->Stop());
      if (dist < fParams.trkClosestApproach()){
        parentTrk = trks[i];
        return true;
      }
    }
    return false;
  }

  double MEFinder::TrkMEDist(me::SlcME me, TVector3 trkEnd){
    double deltaX = ( (me.NXCell() == 0) ? 0 : me.MeanX() - trkEnd.X() );
    double deltaY = ( (me.NYCell() == 0) ? 0 : me.MeanY() - trkEnd.Y() );
    double deltaZ = me.MeanZ() - trkEnd.Z();
    return util::pythag( deltaX, deltaY, deltaZ );
  }
}

namespace me
{
  DEFINE_ART_MODULE(MEFinder)
}