PCHitsList_module.cc

Go to the documentation of this file.

// \file    PCHitsList_module.cc
// \brief   Module to fill vector of PCHits, calibration hit info from cosmic files
// \version $Id: PCHitsList_module.cc,v 1.0 2013/01/13 18:20:04 rtoner Exp $
// \author  rtoner
// \date    $Date: 2013/01/13 22:06:10 $

//C++ and C:
#include <cmath>
#include <iostream>
#include <vector>
#include <utility>

//Framework:
#include "art/Framework/Core/EDProducer.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h" 

//NovaSoft:
#include "CalibrationDataProducts/PCHit.h"
#include "ChannelInfo/BadChanList.h"
#include "DAQChannelMap/DAQChannelMap.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/CellGeo.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "Utilities/func/MathUtil.h"
#include "CalibrationUtils/CalibUtil.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/HitMap.h"
#include "RecoBase/Track.h"
#include "RecoBase/RecoHit.h"
#include "Utilities/AssociationUtil.h"
#include "CMap/service/DetectorService.h"
#include "Simulation/FLSHit.h"

#include "TVector3.h"

namespace caldp { class PCHit; }

namespace calib { 
  
  class PCHitsList : public art::EDProducer {
    
  public:
    explicit PCHitsList(fhicl::ParameterSet const & pset);
    virtual ~PCHitsList();
    
    void produce(art::Event & e) override;
    
    void reconfigure(fhicl::ParameterSet const & p);
    
  private:
    bool IsGoodTrack(rb::Track   const& track,
                     rb::Cluster const& slice);

    bool GoodSteps(rb::Track const& track);
    
    void GetTrueEnergyPathAndLightForCell(const rb::Track*   track,
                                          size_t      const& ihit,
                                          double           & mev,
                                          double           & truePath,
                                          double           & trueW,
                                          double           & truePE,
                                          double           & poissonLambda) const;

    double TrueW( const sim::FLSHit& fls ) const;

    void ProcessTrack(const rb::Track*              track,
                      daqchannelmap::DAQChannelMap* cmap,
                      std::vector<caldp::PCHit>*    pcHitsXY,
                      std::vector<caldp::PCHit>*    pcHitsZ,
                      std::vector<caldp::PCHit>*    pcHitsAvg) const;

    void ProcessTrackTrajectory(const rb::Track*              track,
                                daqchannelmap::DAQChannelMap* cmap,
                                std::set<double>       const& zPlaneBounds,
                                std::vector<caldp::PCHit>*    pcHitsTraj) const;

    void ProcessTrackForBelowThresholdHits(const rb::Track*              track,
                                           daqchannelmap::DAQChannelMap* cmap,
                                           std::vector<caldp::PCHit>*    pcHitsXYThresh) const;

    void ProcessTrackForFLSHits(const rb::Track*              track2,
                                daqchannelmap::DAQChannelMap* cmap,
                                std::vector<caldp::PCHit>*    pcHitsXY) const;

    caldp::PCHit CellHitToPCHit(const rb::CellHit*               chit,
                                const rb::Track*                     trk,
                                double                        const& path,
                                double                        const& mev,
                                double                        const& truePE,
                                double                        const& truePath,
                                double                        const& trueW,
                                double                        const& poissonLambda,
                                daqchannelmap::DAQChannelMap*    cmap) const;

    std::string fCosmicLabel;
    std::string fQualXYName;
    std::string fQualZName;
    std::string fQualAvgName;
    std::string fQualTrajName;
    std::string fQualXYThreshName;
    std::string fQualXYFLSName;

    bool   fIsBeam;
    double fPathCut;            ///< Only keep hits with pathlength less than this
    double fMinPath;            ///< Only keep hits with pathlength greater than this
    int    fExtentZCut;         ///< Only keep tracks crossing this distance in z
    double fDirZCut;            ///< Only keep tracks shallower than this
    double fCompletenessCut;    ///< Only keep tracks with more than this fraction of slice hits included in each view
    double fMaxVtxDistFromEdge; ///< Only keep tracks with an initial point within this dist (cm) from the edge
    double fMaxEndDistFromEdge; ///< Only keep tracks whose final point is not more than this distance from the edge
    double fBadStepSizeLimit;   ///< Only keep tracks without steps outside this limit, see GoodSteps comments
    double fMaxCellsPerPlane;   ///< Only keep tracks with fewer <cells/plane> than this limit
    double fMaxDeltaEndPlane;   ///< Only keep tracks where the difference in end plane between the views is less 
                                ///< than this number
    double fMaxPlaneAsymmetry;  ///< Only keep tracks where the asymmetry in the planes in the views is less than
                                ///< this number

    art::ServiceHandle<cheat::BackTracker>  fBT;
    art::ServiceHandle<geo::Geometry>       fGeom;
    art::ServiceHandle<geo::LiveGeometry>   fLiveGeom;
    art::ServiceHandle<ds::DetectorService> fDS;
  };
}

namespace calib { 

  //---------------------------------------------------------------------//    

  PCHitsList::PCHitsList(fhicl::ParameterSet const & pset)
    : EDProducer(pset)
  {
    reconfigure(pset);
      
    produces< std::vector<caldp::PCHit> >(fQualXYName);
    produces< std::vector<caldp::PCHit> >(fQualZName);
    produces< std::vector<caldp::PCHit> >(fQualAvgName);
    produces< std::vector<caldp::PCHit> >(fQualTrajName);
    produces< std::vector<caldp::PCHit> >(fQualXYThreshName);
    produces< std::vector<caldp::PCHit> >(fQualXYFLSName);

    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYName);
    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualZName);
    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualAvgName);
    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualTrajName);
    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYThreshName);
    produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYFLSName);
  }

  //---------------------------------------------------------------------//
    
  PCHitsList::~PCHitsList()
  {
  }

  //---------------------------------------------------------------------//    

  void PCHitsList::reconfigure(fhicl::ParameterSet const & pset)
  {
    fCosmicLabel        = pset.get< std::string >("CosmicLabel");      // Label of tracker module
    fQualXYName         = pset.get< std::string >("QualXYName");       // Instance label, "XY" quality hits
    fQualZName          = pset.get< std::string >("QualZName");        // Instance label, "Z" quality hits
    fQualAvgName        = pset.get< std::string >("QualAvgName");      // Instance label, "Avg" quality hits
    fQualTrajName       = pset.get< std::string >("QualTrajName");     // Instance label, "Traj" quality hits
    fQualXYThreshName   = pset.get< std::string >("QualXYThreshName"); // Instance label, "XY" quality hits
                                                                       // falling below threshold
    fQualXYFLSName      = pset.get< std::string >("QualXYFLSName");    // Instance label, "XY" quality hits from FLS

    fIsBeam             = pset.get< bool        >("IsBeam");
    fPathCut            = pset.get< double      >("PathCut");
    fMinPath            = pset.get< double      >("MinPathlength");
    fExtentZCut         = pset.get< int         >("ExtentZCut");
    fDirZCut            = pset.get< double      >("DirZCut");
    fCompletenessCut    = pset.get< double      >("CompletenessCut");
    fMaxVtxDistFromEdge = pset.get< double      >("MaxVtxDistFromEdge", 10.);
    fMaxEndDistFromEdge = pset.get< double      >("MaxEndDistFromEdge", 10.);
    fBadStepSizeLimit   = pset.get< double      >("BadStepSizeLimit",   2. );
    fMaxCellsPerPlane   = pset.get< double      >("MaxCellsPerPlane",   6. );
    fMaxDeltaEndPlane   = pset.get< double      >("MaxDeltaEndPlane",   3. );
    fMaxPlaneAsymmetry  = pset.get< double      >("MaxPlaneAsymmetry",  0.1);
    
    return;
  }
    
  //---------------------------------------------------------------------//    
    
  void PCHitsList::produce(art::Event & e)
  {
    //Ptrs to vectors to be filled (one instance for each degree of hit quality)
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsXY      (new std::vector<caldp::PCHit>);
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsZ       (new std::vector<caldp::PCHit>);
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsAvg     (new std::vector<caldp::PCHit>);
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsTraj    (new std::vector<caldp::PCHit>);
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsXYThresh(new std::vector<caldp::PCHit>);
    std::unique_ptr< std::vector<caldp::PCHit>  > pcHitsXYFLS   (new std::vector<caldp::PCHit>);
      
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXY      (new art::Assns<caldp::PCHit, rb::Track>);
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsZ       (new art::Assns<caldp::PCHit, rb::Track>);
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsAvg     (new art::Assns<caldp::PCHit, rb::Track>);
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsTraj    (new art::Assns<caldp::PCHit, rb::Track>);
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXYThresh(new art::Assns<caldp::PCHit, rb::Track>);
    std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXYFLS   (new art::Assns<caldp::PCHit, rb::Track>);

    //Get cosmic tracks to use for calibration:
    art::Handle<std::vector<rb::Track> > trackcol;
    e.getByLabel(fCosmicLabel, trackcol);
      
    std::vector< art::Ptr<rb::Track> > tracks;
    art::fill_ptr_vector(tracks, trackcol);

    art::FindManyP<rb::Cluster> fmpslc(trackcol, e, fCosmicLabel);

    daqchannelmap::DAQChannelMap* cmap = daqchannelmap::DAQChannelMap::getInstance(fDS->DetId());
    
    // prepare to find the path lengths in a cell based on the track trajectory
    std::set<double> planeZBounds;
    calib::FindZBoundaries(planeZBounds);
    
    for(size_t itrk = 0; itrk < tracks.size(); ++itrk){

      // Find the slice associated with this track
      // There should be exactly one
      const std::vector<art::Ptr<rb::Cluster> > slices = fmpslc.at(itrk);
      if(slices.size() != 1)
        throw cet::exception("PCHitsList") << "incorrect number of slices associated to track";
      
      const rb::Cluster& slice = *slices[0];

      // Remember how many hits of each type we had so far, so we know what
      // index to start associating from when we're done with this track.
      const size_t xyBeginIdx       = pcHitsXY->size();
      const size_t zBeginIdx        = pcHitsZ->size();
      const size_t avgBeginIdx      = pcHitsAvg->size();
      const size_t trajBeginIdx     = pcHitsTraj->size();
      const size_t xyThreshBeginIdx = pcHitsXYThresh->size();
      const size_t xyFLSBeginIdx    = pcHitsXYFLS->size();

      const art::Ptr<rb::Track> track = tracks[itrk];

      MF_LOG_DEBUG("PCHitsList") << "TRACK " << itrk << ": " << track->NCell() << " cells to process.";
        
      //-------Track cuts:----------------//
      if( !this->IsGoodTrack(*track, slice) ) continue;

      // track passes the cuts, so process it
      ProcessTrack(track.get(), 
                   cmap,
                   pcHitsXY.get(),
                   pcHitsZ.get(),
                   pcHitsAvg.get());

      ProcessTrackTrajectory(track.get(),
                             cmap,
                             planeZBounds,
                             pcHitsTraj.get());

      ProcessTrackForBelowThresholdHits(track.get(),
                                        cmap,
                                        pcHitsXYThresh.get());

      ProcessTrackForFLSHits(track.get(), 
                             cmap,
                             pcHitsXYFLS.get());


      // Associate all the PCHits we created for this track. We do this in
      // bulk because making lots of individual calls as we go is slow.
      util::CreateAssn(e, *pcHitsXY,       track, *hitAssnsXY,       xyBeginIdx,       UINT_MAX, fQualXYName);
      util::CreateAssn(e, *pcHitsZ,        track, *hitAssnsZ,        zBeginIdx,        UINT_MAX, fQualZName);
      util::CreateAssn(e, *pcHitsAvg,      track, *hitAssnsAvg,      avgBeginIdx,      UINT_MAX, fQualAvgName);
      util::CreateAssn(e, *pcHitsTraj,     track, *hitAssnsTraj,     trajBeginIdx,     UINT_MAX, fQualTrajName);
      util::CreateAssn(e, *pcHitsXYThresh, track, *hitAssnsXYThresh, xyThreshBeginIdx, UINT_MAX, fQualXYThreshName);
      util::CreateAssn(e, *pcHitsXYFLS,    track, *hitAssnsXYFLS,    xyFLSBeginIdx,    UINT_MAX, fQualXYFLSName);

    }//End track loop
      
    MF_LOG_DEBUG("PCHitsList")
      << "SUMMARY: Size of XY vector  = " << pcHitsXY->size() 
      << ", Size of Z vector   = " << pcHitsZ->size() 
      << ", Size of Avg vector = " << pcHitsAvg->size();
      
    e.put(std::move(pcHitsXY),         fQualXYName);
    e.put(std::move(pcHitsZ),          fQualZName);
    e.put(std::move(pcHitsAvg),        fQualAvgName);
    e.put(std::move(pcHitsTraj),       fQualTrajName);
    e.put(std::move(pcHitsXYThresh),   fQualXYThreshName);
    e.put(std::move(pcHitsXYFLS),      fQualXYFLSName);
    e.put(std::move(hitAssnsXY),       fQualXYName);
    e.put(std::move(hitAssnsZ),        fQualZName);
    e.put(std::move(hitAssnsAvg),      fQualAvgName);
    e.put(std::move(hitAssnsTraj),     fQualTrajName);
    e.put(std::move(hitAssnsXYThresh), fQualXYThreshName);
    e.put(std::move(hitAssnsXYFLS),    fQualXYFLSName);
     
    MF_LOG_DEBUG("PCHitsList") << "Event finished successfully.";

  }//end produce


  //---------------------------------------------------------------------//
  bool PCHitsList::IsGoodTrack(rb::Track   const& track,
                               rb::Cluster const& slice)
  {
    MF_LOG_DEBUG("PCHitsList") << "NCellsX=" << track.NXCell() << " NCellsY=" << track.NYCell();
    
    //Track must have at least two X cells and 2 Y cells:
    if(track.NXCell() < 2 || track.NYCell() < 2) { 
      MF_LOG_DEBUG("PCHitsList") << "!---->Failed NCell Cut!";
      return false; 
    }
        
    //Check that direction vector is well-defined:
    const TVector3 dr = track.Dir().Unit();
    const double   dx = std::abs(dr.X());
    const double   dy = std::abs(dr.Y());
    const double   dz = std::abs(dr.Z());

    //No undefined directions:
    if(std::isnan(dx) || 
       std::isnan(dy) || 
       std::isnan(dz)) return false;
    
    //At least one direction must be non-zero:
    if(dx == 0 && 
       dy == 0 && 
       dz == 0) {
      MF_LOG_DEBUG("PCHitsList") << "!---->Failed Direction Cut!";
      return false;
    }
    
    //Must be a unit vector (otherwise, abort event, something went very wrong):
    if ( std::abs( util::pythag(dx, dy, dz)-1 ) > 1e-4){
      throw cet::exception("SkipEvent") << "ERROR: Direction was not a unit vector.  Skipping Event.\n";
    }

    // Additional quality cuts. These remove tracks that often have
    // badly-reconstructed W positions.

    // Cut very steep tracks two ways. These can be misreconstructed, or
    // they can simply make it difficult for the W estimation to do a good
    // job.  OK to use the rb::Track::Dir() method here because the initial
    // direction of the track is what we want to cut on
    if(std::abs(track.Start().Z() - track.Stop().Z() )  <= fExtentZCut || 
       std::abs(track.Dir().Z()) <  fDirZCut) return false;

    // Cut out actual reconstruction failures, where the track doesn't use
    // a significant amount of the slice hits.  If either view is not
    // complete enough discard the whole thing
    if(track.NXCell() < fCompletenessCut*slice.NXCell() || 
       track.NYCell() < fCompletenessCut*slice.NYCell()) return false;

    // Check that the average number of cells per plane is similar in 
    // each view.  If the value is very high in either view, or larger in 
    // one view than the other then you might have some form of FEB flash
    // or cosmic ray brem that will mess up the calibration
    std::set<unsigned short> xPlaneList;
    std::set<unsigned short> yPlaneList;
    for(size_t c = 0; c < track.NCell(); ++c){
      if     (track.Cell(c)->View() == geo::kX) xPlaneList.insert(track.Cell(c)->Plane());
      else if(track.Cell(c)->View() == geo::kY) yPlaneList.insert(track.Cell(c)->Plane());
    }// end loop over cells in track
    
    double xCells  = track.NXCell();
    double yCells  = track.NYCell();
    double xPlanes = 1. * xPlaneList.size();
    double yPlanes = 1. * yPlaneList.size();
    double cellsPerPlaneX = (xPlanes > 0) ? xCells/xPlanes : 999.;
    double cellsPerPlaneY = (yPlanes > 0) ? yCells/yPlanes : 999.;
    if( std::max(cellsPerPlaneX, cellsPerPlaneY) > fMaxCellsPerPlane ) return false;

    // remove tracks where there is an asymmetry between the number 
    // of planes hit in each view.  One class of asymmetry is where 
    // there are extra planes off either end of the track in one view.
    // Since there is no information in the other view for those planes
    // you can't be sure you got the length/containment correct
    double xStart = track.MaxPlane(geo::kX);
    double yStart = track.MaxPlane(geo::kY);
    double xEnd   = track.MinPlane(geo::kX);
    double yEnd   = track.MinPlane(geo::kY);
    if(std::abs(xStart - yStart)   > fMaxDeltaEndPlane ||
       std::abs(xEnd   - yEnd  )   > fMaxDeltaEndPlane ||
       std::abs(xPlanes - yPlanes) > fMaxPlaneAsymmetry*(xPlanes + yPlanes) ) return false;
    

    // Check that the track doesn't take one or two big steps
    // compared to the rest, that is usually a reconstruction failure.
    if(!this->GoodSteps(track) ) return false;

    int junk = -999;

    // Remove tracks whose vertex position is well away from the edge of the
    // detector. Muon catcher is not counted as part of the detector proper in
    // ND so have to check for that as well.
    if( !fIsBeam && std::abs(fLiveGeom->DistToEdgeXY(track.Start())) > fMaxVtxDistFromEdge &&
       std::abs(fLiveGeom->DistToEdgeZ (track.Start())) > fMaxVtxDistFromEdge &&
       (fDS->DetId() == novadaq::cnv::kFARDET ||
        std::abs(fLiveGeom->DistanceToEdgeInMC(track.Start(), junk)) > fMaxVtxDistFromEdge) ) {
         return false;
    }

    // Remove tracks whose end position is more than fMaxEndDistFromEdge
    // outside of the detector The minus sign is because LiveGeometry returns a
    // negative value if the point is outside of the active detector.  A track
    // can stop anywhere inside the detector and still be OK.
    if( !fIsBeam && (fLiveGeom->DistToEdgeXY(track.Stop()) < -fMaxEndDistFromEdge ||
         fLiveGeom->DistToEdgeZ (track.Stop()) < -fMaxEndDistFromEdge) &&
       (fDS->DetId() == novadaq::cnv::kFARDET ||
        fLiveGeom->DistanceToEdgeInMC(track.Stop(), junk) < -fMaxEndDistFromEdge) ) {
         return false;
       }

    return true;
  }


  //......................................................................  
  // tracks that have huge changes in the step size from one trajectory
  // point to another are probably poorly reconstructed, and should be removed
  // from the analysis.  Usually the biggest steps are at either the start 
  // or end of the track in cases of really bad reconstruction
  bool PCHitsList::GoodSteps(rb::Track const& track)
  {

    // first find the median step size
    size_t size = track.NTrajectoryPoints() - 1;

    // track with only two distances between trajectory points, 
    // probably not very useful for the analysis so return false
    if(size < 3) return false;

    std::vector<float> stepsize(size, 0.);
    float largestStep = -1.;
    for(size_t t = 1; t < size+1; ++t){
      TVector3 const& tp1 = track.TrajectoryPoint(t-1);
      TVector3 const& tp2 = track.TrajectoryPoint(t);
      stepsize[t-1] = (tp1 - tp2).Mag();
      largestStep = std::max(stepsize[t-1], largestStep);
    }

    // sort the vector of step sizes and find the median
    std::sort(stepsize.begin(), stepsize.end());

    float median = 0.;
    if(size%2 == 0) 
      median = 0.5*(stepsize[size/2] + stepsize[(size/2) - 1]);
    else
      median = stepsize[size/2];

    // the median represents the value of the step size 
    // distribution where half of the distribution is
    // below that value and half is above it.  If the 
    // largest step is more than twice the median, then 
    // it is very far outside the distribution and this
    // track was not well reconstructed
    if(largestStep > median*fBadStepSizeLimit) return false;

    return true;
  }

  //---------------------------------------------------------------------//
  // all input references are assumed to have been set to zero by the calling method
  void PCHitsList::GetTrueEnergyPathAndLightForCell(const rb::Track*   track,
                                                    size_t      const& ihit,
                                                    double           & mev,
                                                    double           & truePath,
                                                    double           & trueW,
                                                    double           & truePE,
                                                    double           & poissonLambda) const
  {
    
    if (fBT->HaveTruthInfo()){
      const std::vector<cheat::TrackIDE> ides = fBT->HitToTrackIDE(track->Cell(ihit));
      
      for(size_t n = 0; n < ides.size(); ++n){
        mev += ides[n].energy*1000;

        // the ides are sorted to have the most contributing track id be first in the list
        // so grab that track id and get the particle it corresponds to
        if(n < 1){
          auto particle = fBT->ParticleNavigator()[ides[n].trackID];
          
          std::vector<sim::FLSHit> flsHits = fBT->ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
          for(auto fls : flsHits){
            if(fls.GetPlaneID() == track->Cell(ihit)->Plane() &&
               fls.GetCellID()  == track->Cell(ihit)->Cell() ){
              truePath = fls.GetTotalPathLength();
              trueW = TrueW(fls);
              break;
            }
          }
        } // end if on the first ide
      } // end loop over ides
      
      const std::vector<sim::PhotonSignal> pes = fBT->HitToPhotonSignal(track->Cell(ihit));
      for(size_t n = 0; n < pes.size(); ++n){
        truePE += pes[n].NPhoton();
        poissonLambda += pes[n].PoissonLambda();
      }
      
    }

    return;
  }



  //---------------------------------------------------------------------//
  double PCHitsList::TrueW(const sim::FLSHit& fls) const
  {
    double local_enter[3] = {fls.GetEntryX(), fls.GetEntryY(), fls.GetEntryZ()};
    double local_exit[3]  = {fls.GetExitX(),  fls.GetExitY(),  fls.GetExitZ()};
    double enter[3], exit[3];
    int ip = fls.GetPlaneID();
    int ic = fls.GetCellID();

    (fGeom->Plane(ip)->Cell(ic))->LocalToWorld(local_enter,enter);
    (fGeom->Plane(ip)->Cell(ic))->LocalToWorld(local_exit ,exit);

    if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kX)
      return ( enter[1] + exit[1] )/2;
    else if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kY){
      // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
      if(fGeom->DetId() == novadaq::cnv::kTESTBEAM) 
        return -( enter[0] + exit[0] )/2;
      return ( enter[0] + exit[0] )/2;
    }
    else
      mf::LogError("TrueW") << "bad view for fls hit";
    return 9999.;
  }

  
  //---------------------------------------------------------------------//
  void PCHitsList::ProcessTrack(const rb::Track*              track,
                                daqchannelmap::DAQChannelMap* cmap,
                                std::vector<caldp::PCHit>*    pcHitsXY,
                                std::vector<caldp::PCHit>*    pcHitsZ,
                                std::vector<caldp::PCHit>*    pcHitsAvg) const
  {
    std::vector<rb::CellHit> chits;
    for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
      chits.push_back(*track->Cell(ihit).get());
    }

    //Calculate path quality:
    std::vector<EPathQuality> quals;
    const std::vector<double> paths = BestPathEstimates(track->OfflineChans(), track->PlaneDirMap(), quals);

    //Check there is a path and quality measure for each cell; if not, something wrong, abort event:
    if(chits.size() != quals.size() ||
       chits.size() != paths.size()){
      throw cet::exception("SkipEvent") << "ERROR: Number of pathlengths or qualities != number of cells.  Skipping Event.\n";
    }

    //Loop through cells in track:
    for(size_t ihit = 0; ihit < chits.size(); ++ihit){
          
      //Get hit
      const rb::CellHit& chit = chits[ihit];

      //Get path:
      const double path = paths[ihit]; //RBT: Way of ensuring cell and path line up??

      //Check path not too long:
      if(path > fPathCut) continue;

      //Get the correct MeV and true PE in this hit
      double mev           = 0.;
      double truePE        = 0.;
      double poissonLambda = 0.;
      double truePath      = 0.;
      double trueW         = -9999.;
      this->GetTrueEnergyPathAndLightForCell(track, ihit, mev, truePath, trueW, truePE, poissonLambda);

      caldp::PCHit pchit = CellHitToPCHit(&chit, track,
                                          path, mev, truePE, truePath, trueW, poissonLambda,
                                          cmap);

      switch(quals[ihit]){

        case kXYAdjacents:
          MF_LOG_DEBUG("PCHitsList")
          << "Quality XY!  - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
          << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
          << " W=" << track->W(&chit) << " mev=" << mev;
          pcHitsXY->push_back(pchit);
          break;
          
        case kZAdjacents:
          MF_LOG_DEBUG("PCHitsList")
          << "Quality Z!   - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
          << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
          << " W=" << track->W(&chit) << " mev=" << mev;
          pcHitsZ->push_back(pchit);
          break;
          
        case kAverage:
          MF_LOG_DEBUG("PCHitsList")
          << "Quality Avg! - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
          << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
          << " W=" << track->W(&chit) << " mev=" << mev;
          pcHitsAvg->push_back(pchit);
          break;
          
        case kBadTrack:
          MF_LOG_DEBUG("PCHitsList") << "!!!!!!BAD hit/track!!!!!! - hit=" << ihit;
            // Something wrong with this hit/track, so do nothing
          break;
          
        default:
          throw cet::exception("SkipEvent") << "ERROR: Unknown Path Quality type.  Skipping Event.\n";
            //Unknown path type; something went very wrong; abort this event.
      }
      
    }//End cell loop

  }//end ProcessTrack

  //---------------------------------------------------------------------//
  void PCHitsList::ProcessTrackTrajectory(const rb::Track*              track,
                                          daqchannelmap::DAQChannelMap* cmap,
                                          std::set<double>       const& zPlaneBounds,
                                          std::vector<caldp::PCHit>*    pcHitsTraj) const
  {
    auto boundaryMap = calib::MakeZBoundaryMap(zPlaneBounds, track->Trajectory());

    // don't bother with this track if no z boundaries were found, that implies
    // that all of the trajectory points are outside of the live geometry
    if(boundaryMap.empty()) return;

    for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
      const rb::CellHit& chit = *track->Cell(ihit);

      std::pair<uint32_t, uint32_t> pc = {chit.Plane(), chit.Cell()};

      double xyz[3];
      fGeom->Plane(chit.Plane())->Cell(chit.Cell())->GetCenter(xyz);
      xyz[1-chit.View()] = track->W(&chit);
      // TODO - wouldn't this be a better estimate of the hit position?
      // track->InterpolateXY(xyz[2], xyz[0], xyz[1]);
      TVector3 point(xyz[0], xyz[1], xyz[2]);

      MF_LOG_DEBUG("PCHitsList")
      << ihit << "/" << track->NCell()
      << " recohit: " << point.X()
      << " " << point.Y() << " " << point.Z();
      
      // don't bother with the point if it is reconstructed outside of the detector
      // silly but apparently true, IsPointLive cannot handle the muon catcher
      // IsMuonCatcher returns false if in the FD, so in the FD this check
      // should be determined by IsPointLive
      if( !fLiveGeom->IsPointLive(point) && !fLiveGeom->InMuonCatcher(point) ){
        MF_LOG_DEBUG("PCHitsList")
        << ihit
        << " is not valid "
        << "(x,y,z) = "
        << point.X() << "," << point.Y() << "," << point.Z();
        continue;
      }

      try{
        float path = calib::PathLengthInCell(boundaryMap, point, pc);
        
        //Check path not too long:
        if(path > fPathCut || path < fMinPath) continue;
        
        //Get the correct MeV and true PE in this hit
        double mev           = 0.;
        double truePE        = 0.;
        double trueW         = -9999.;
        double truePath      = 0.;
        double poissonLambda = 0.;
        this->GetTrueEnergyPathAndLightForCell(track, ihit, mev, truePath, trueW, truePE, poissonLambda);

        pcHitsTraj->push_back(CellHitToPCHit(&chit, track,
                                             path, mev, truePE, truePath, trueW, poissonLambda,
                                             cmap));
      }
      catch(cet::exception &e){
        MF_LOG_WARNING("PCHitsList")
        << "caught exception \n" << e
        << "\n cell " << pc.first << " " << pc.second << " not added to PCHitsList";
      }
      
    }//End cell loop
    
  }//end ProcessTrackTrajectory

  //---------------------------------------------------------------------//
  void PCHitsList::ProcessTrackForBelowThresholdHits(const rb::Track*              track,
                                                     daqchannelmap::DAQChannelMap* cmap,
                                                     std::vector<caldp::PCHit>*    pcHitsXYThresh) const
  {
    std::map<geo::OfflineChan, double> trueE;
    std::map<geo::OfflineChan, double> truePE;
    std::map<geo::OfflineChan, double> truePath;
    std::map<geo::OfflineChan, double> trueW;
    std::map<geo::OfflineChan, double> poissonLambda;
    if(fBT->HaveTruthInfo()){
      const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track->AllCells());
      if(!parts.empty()){
        const int id = parts[0]->TrackId(); // Biggest contributor
        std::set<geo::OfflineChan> flsChans;
        const std::vector<sim::FLSHit> flss = fBT->ParticleToFLSHit(id);
        for(const sim::FLSHit& fls: flss){
          const geo::OfflineChan chan(fls.GetPlaneID(), fls.GetCellID());
          trueE[chan] += 1000*fls.GetEdep();
          truePath[chan] = fls.GetTotalPathLength();
          trueW[chan] = TrueW(fls);
          flsChans.insert(chan);
        }
        
        for(geo::OfflineChan chan: flsChans){
          // We won't have an actual CellHit for this cell, so have to use
          // CellTo not HitTo. Could try to do something clever with matching
          // the times up. But this must be a small effect.
          const std::vector<sim::PhotonSignal> pes = fBT->CellToPhotonSignal(chan.Plane(), chan.Cell());
          for(sim::PhotonSignal pe: pes){
            truePE[chan] += pe.NPhoton();
            poissonLambda[chan] += pe.PoissonLambda();
          }
        }
      }
    }

    // Find cells that would have made it into the sample if they'd been
    // hit. Downstream analyzers can interpret these as a threshold effect.
    std::vector<EPathQuality> threshQuals;
    std::vector<double> threshPaths;
    std::vector<geo::OfflineChan> threshChans;

    FindBelowThresholdCalibCandidates(track->OfflineChans(), track->PlaneDirMap(), threshQuals, threshPaths, threshChans);

    for(unsigned int threshIdx = 0; threshIdx < threshQuals.size(); ++threshIdx){
      const double path = threshPaths[threshIdx];
      // Check path not too long. Apply same cut as normal candidates so as
      // not to bias ourselves.
      if(path > 10) continue;

      const int plane = threshChans[threshIdx].Plane();
      const int cell = threshChans[threshIdx].Cell();
      const geo::OfflineChan chan(plane, cell);
      const geo::View_t view = fGeom->Plane(plane)->View();
      rb::CellHit dummy;
      dummy.SetPlane(plane);
      dummy.SetCell(cell);
      dummy.SetView(view);

      caldp::PCHit pchit = CellHitToPCHit(&dummy, track,
                                          path, trueE[chan],
                                          truePE[chan], truePath[chan], trueW[chan], poissonLambda[chan],
                                          cmap);

      if(threshQuals[threshIdx] == kXYAdjacents)
      pcHitsXYThresh->push_back(pchit);
    } // end for threshIdx

  } // end ProcessTrackForBelowThresholdHits

  //---------------------------------------------------------------------//
  void PCHitsList::ProcessTrackForFLSHits(const rb::Track*              track,
                                          daqchannelmap::DAQChannelMap* cmap,
                                          std::vector<caldp::PCHit>*    pcHitsXY) const
  {
    if(!fBT->HaveTruthInfo()) return;

    rb::HitMap hmap(track->AllCells());

    // BackTracker won't be able to figure out the true energy, since
    // there's no CellHit to work back from. Accumulate it here.
    std::map<geo::OfflineChan, double> extrasTruthMap;
    std::map<geo::OfflineChan, double> extrasPEMap;
    std::map<geo::OfflineChan, double> extrasTruePath;
    std::map<geo::OfflineChan, double> extrasTrueW;
    std::map<geo::OfflineChan, double> extrasPoissonLambda;
    std::vector<geo::OfflineChan> extras;

    const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track->AllCells());
    if(!parts.empty()){
      const int id = parts[0]->TrackId(); // Biggest contributor
      const std::vector<sim::FLSHit> flss = fBT->ParticleToFLSHit(id);
      const unsigned int nFLS = flss.size();
      std::set<geo::OfflineChan> extraset;
      std::set<geo::OfflineChan> flsCells;
      for(unsigned int flsIdx = 0; flsIdx < nFLS; ++flsIdx){
        const sim::FLSHit& fls = flss[flsIdx];
        const geo::OfflineChan chan(fls.GetPlaneID(), fls.GetCellID());
        extraset.insert(chan);
        extrasTruthMap[chan] += 1000*fls.GetEdep();
        extrasTruePath[chan] = fls.GetTotalPathLength();
        if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kX)
          extrasTrueW[chan] = fls.GetYAverage();
        else if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kY)
          extrasTrueW[chan] = fls.GetXAverage();
        flsCells.insert(chan);
      }

      for(geo::OfflineChan chan: flsCells){
        // We might not have an actual CellHit for this cell, so have to use
        // CellTo not HitTo. Could try to do something clever with matching
        // the times up. But this must be a small effect.
        const std::vector<sim::PhotonSignal> pes = fBT->CellToPhotonSignal(chan.Plane(), chan.Cell());
        for(sim::PhotonSignal pe: pes){
          extrasPEMap[chan] += pe.NPhoton();
          extrasPoissonLambda[chan] += pe.PoissonLambda();
        }
      }

      // extraset is a set to get a collection of distinct cells, but
      // later on we need it to be a vector.
      extras.insert(extras.end(), extraset.begin(), extraset.end());
    }

    const double meant = track->MeanTNS();

    rb::Track trackNoHits(*track);
    trackNoHits.Clear();

    //Calculate path quality
    std::vector<EPathQuality> quals;
    const std::vector<double> paths = BestPathEstimates(trackNoHits.OfflineChans(), trackNoHits.PlaneDirMap(), quals, extras);

    //Loop through cells in track:
    for(unsigned int ihit = 0; ihit < extras.size(); ++ihit){
      const geo::OfflineChan chan = extras[ihit];

      //Get path:
      const double path = paths[ihit];
      //Check path not too long:
      if(path > 10) continue;

      if(quals[ihit] == kXYAdjacents){
          //Get the correct MeV in this hit
        const double mev = extrasTruthMap[chan];
        const double truePE = extrasPEMap[chan];
        const double truePath = extrasTruePath[chan];
        const double trueW    = extrasTrueW[chan];
        const double poissonLambda = extrasPoissonLambda[chan];
        
        rb::CellHit chit;
        chit.SetPlane(chan.Plane());
        chit.SetCell(chan.Cell());
        chit.SetView(fGeom->Plane(chan.Plane())->View());
        if(hmap.CellExists(chan.Plane(), chan.Cell())){
          chit.SetPE(hmap.Cell(chan.Plane(), chan.Cell())->PE());
        }
        else{
          chit.SetPE(0);
        }
        chit.SetTNS(meant, true);
        
        caldp::PCHit pchit = CellHitToPCHit(&chit, track,
                                            path, mev, truePE, truePath, trueW, poissonLambda,
                                            cmap);
        
        pcHitsXY->push_back(pchit);
      }
    }//End cell loop

  }//end ProcessTrackForFLSHit

  //---------------------------------------------------------------------//
  caldp::PCHit PCHitsList::CellHitToPCHit(const rb::CellHit*            chit,
                                          const rb::Track*                      trk,
                                          double                        const&  path,
                                          double                        const&  mev,
                                          double                        const&  truePE,
                                          double                        const&  truePath,
                                          double                        const&  trueW,
                                          double                        const&  poissonLambda,
                                          daqchannelmap::DAQChannelMap* cmap) const
  {
    //get more information on hit location
    uint32_t lchan   = cmap->encodeLChan(fGeom->DetId(), chit->Plane(), chit->Cell());
    uint32_t dchan   = cmap->encodeDChan(lchan);
    uint32_t pixel   = cmap->getPixel(dchan);
    uint32_t apd     = cmap->getFEB(dchan);
    uint32_t dcm     = cmap->getDCM(dchan);
    uint32_t diblock = cmap->getDiBlock(dchan);

    double w = trk->W(chit);

    //get distance to readout electronics
    double xyzd[4];
    GetXYZD(chit->OfflineChan(), w, xyzd);
    //make sure the reconstructed hit coordinate is inside the physical
    //detector if not then change the distance to the electronics to -1.0
    if(std::abs(xyzd[0]) > fGeom->DetHalfWidth()  || 
       std::abs(xyzd[1]) > fGeom->DetHalfHeight() || 
                xyzd[2]  > fGeom->DetLength()     || 
                xyzd[2]  < 0.0) xyzd[3] = -1.0;

    double flightlen = trk->DistanceFromStart(xyzd[2]); 

    //Define and fill a PCHit for this cell hit
    caldp::PCHit pchit;
    pchit.SetPlane        ( chit->Plane()     ); //Set plane number                        
    pchit.SetCell         ( chit->Cell()      ); //Set cell number                         
    pchit.SetPixel        ( pixel             ); //Set pixel number                        
    pchit.SetAPD          ( apd               ); //Set apd number                          
    pchit.SetDCM          ( dcm               ); //Set dcm number                          
    pchit.SetDiblock      ( diblock           ); //set diblock number              
    pchit.SetView         ( chit->View()      ); //set hit view                    
    pchit.SetPE           ( chit->PE()        ); //Set uncalibrated PE             
    pchit.SetPath         ( path              ); //Set path length                         
    pchit.SetW            ( w                 ); //Set W (unmeasured coordinate) value 
    pchit.SetTrueMeV      ( mev               ); //Set true energy of hit                  
    pchit.SetTruePE       ( truePE            ); //Set true PE of hit
    pchit.SetTruePath     ( truePath          ); //Set true pathlength through cell
    pchit.SetTrueW        ( trueW             );
    pchit.SetPoissonLambda(poissonLambda      );
    pchit.SetTNS          ( chit->TNS()       ); //Set time of hit                                 
    pchit.SetFlightLen    ( flightlen         ); //Set path from start of track to hit     
    pchit.SetReadDist     ( xyzd[3]           ); //Set distance from hit to readout electronics
    pchit.SetGoodTime     ( chit->GoodTiming()); //Set if cell had a proper time fit           

    return pchit;
  } // end CellHitToPCHit
  
} //end namespace

//---------------------------------------------------------------------//    

DEFINE_ART_MODULE(calib::PCHitsList)