MuonRemove_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \file    MuonRemove.cxx
/// \brief   Remove raw digits corresponding to muon hits.
/// \author  Kanika Sachdev (ksachdev@physics.umn.edu)
////////////////////////////////////////////////////////////////////////



// C/C++ includes
#include <cmath>
#include <iostream>
#include <map>
#include <vector>
#include <algorithm>
#include <list>
// ROOT includes
#include "TMath.h"
#include "TVector3.h"
#include "TDatabasePDG.h"

// Framework includes
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "CVN/func/Result.h"
#include "StandardRecord/CVNLabels.h"
#include "RecoBase/FilterList.h"

// NOvASoft includes
#include "NovaDAQConventions/DAQConventions.h"
#include "DAQDataFormats/NanoSliceVersionConvention.h"
#include "MuonRemove/art/TrackCleanUpAlg.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "CMap/service/CMapService.h"
#include "MCCheater/BackTracker.h"
#include "Simulation/Particle.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/FLSHit.h"
#include "Simulation/PhotonSignal.h"
#include "RawData/RawDigit.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/WeightedHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/FitSum.h"
#include "Calibrator/Calibrator.h"
#include "Utilities/AssociationUtil.h"
#include "ReMId/ReMId.h"
#include "NumuEnergy/NumuEAlg.h"
#include "NumuEnergy/NumuE.h"

namespace murem
{
  class MuonRemove : public art::EDProducer
  {
  public:

    explicit MuonRemove(fhicl::ParameterSet const& pset);
    ~MuonRemove();

    void beginRun(art::Run& run);

    void produce(art::Event& evt);

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

    ///\brief Returns a vector of Cellhits attributed to the muon by checking with truth.
    /// For implementing muon removal by truth
    std::vector< rb::WeightedHit > RemoveByTruth(const art::Ptr<rb::Cluster>& slice);

    ///\brief A function that returns the longest track amongst a vector of tracks
    rb::Track LongestTrack( std::vector< art::Ptr<rb::Track> >& tracks, art::Event& evt, double& muonEnergy);

    ///\brief Returns the highest efficiency (by truth) muon track for removal;
    rb::Track RemoveByEfficiency( std::vector< art::Ptr<rb::Track> >& tracks,
                                  const rb::Cluster &slice, art::Event& evt, double& muonEnergy);

    ///\brief Returns the track with the highest PID value returned by the ReMId package
    rb::Track RemoveByReMId( std::vector< art::Ptr<rb::Track> >& tracks,
                             art::Event& evt, double& muonEnergy);
    ///\brief Returns the BPF muon track with the highest Prong CVN value
    rb::Track RemoveByProngCVN( std::vector< art::Ptr<rb::Track> >& tracks, const art::Ptr<rb::Cluster>& slice,
                                art::Event& evt, double& muonEnergy);

    ///\brief Function for preselecting the slices for the further mrcc process
    bool SliceSelection ( const art::Event& evt, const art::Ptr<rb::Cluster>& slice);

    ///\brief Saves the information about the removed muon in a simb::MCParticle object
    /// The energy of the muon is obtained from NumuEAlg.
    simb::MCParticle MuonInfo(rb::Track muonTrack, double muonEnergy);

    ///\brief Saves muon removed FLSHitLists, PhotonSignals, Particle list and
    /// information about the removed muon in a simb::MCParticle.
    void FillTruthInfo( const art::Event& evt,
                        std::unique_ptr< std::vector<sim::FLSHitList>    >&     mrccFLSHits,
                        std::unique_ptr< std::vector<sim::PhotonSignal>  >&     mrccPhotSig,
                        std::unique_ptr< std::vector<sim::Particle>      >&     mrccPartList,
                        std::unique_ptr< std::vector<simb::MCParticle>   >&     muonInfo) ;

  private:

    /// An object of TrackCleanUpAlg
    TrackCleanUpAlg*    fTrkCleanUpAlg;

    /// ParameterSet of TrackCleanUpAlg
    fhicl::ParameterSet fTrkCleanUpPSet;

    /// An object of NumuEAlg
    numue::NumuEAlg*    fNumuEAlg;

    /// ParameterSet of NumuEAlg
    fhicl::ParameterSet fNumuEAlgPSet;

    /// A FHiCL flag to turn of track clean up. Default if true.
    bool                fTrackCleanUp;

    /// ADC thresholds

    float fADCThreshold_Far;
    float fADCThreshold_NDOS;
    float fADCThreshold_Near;
    float fADCThreshold;

    // Configurable module parameters
    /// Label of module that produced verteces
    std::string         fElasticInput;
    /// Label of module that produced prongs
    std::string         fFuzzy3DInput;
    /// Label of module that produced Prong CVN PIDs
    std::string         fCVNParticleLabel;
    /// Label of module that produced BPF tracks
    std::string         fBPFTrackModuleLabel;
    /// Label of module that produced slice CVN PID
    std::string         fCVNLabel;

    /// Label of module which will be used for track loading
    std::string         fWhichTrackModule;

    /// Label of module that produced raw digits
    std::string         fRawDataLabel;
    /// Label of module that produced tracks
    std::string         fKalmanTrackModuleLabel;
    /// Label of module that produced slice clusters
    std::string         fClusterModuleLabel;
    /// Label of module that produced muon pid object
    std::string         fMuonID;
    /// Label of module that produced MC truth information
    std::string         fMCTruthModuleLabel;
    /// Label of module that produced Numu PID
    std::string         fNumuReMIdLabel;
    /// Label of module that produced Energy object (in case of Kalman tracks)
    std::string         fEnergyLabel;
    /// Label of module that produced Energy object (in case of BPF tracks)
    std::string         fEnergyLabelBPF;
    
    /// Numu PID cut
    double              fNumuPIDCut;
    /// DCS Distance: multipoint readout look back
    int                 fDCSDistance;

    /// Handle to Geometry
    art::ServiceHandle<geo::Geometry> fgeom;
  };

  /// \brief Takes in a vector of raw digits
  /// and sorts it by plane.
  void SortDigitsByPlaneAndCell( std::list< rawdata::RawDigit >& digitcol );

  /// \brief Compare two rawdigits by plane. Return true if
  /// plane of ra is greater then the plane of rb.
  bool CompareDigitsByPlaneAndCell( rawdata::RawDigit ra,
                                    rawdata::RawDigit rb);

  void SortByPlaneAndCell( std::vector< rb::WeightedHit >& whits );

  /// \brief Compare two rawdigits by plane. Return true if
  /// plane of ra is greater then the plane of rb.
  bool CompareByPlaneAndCell( rb::WeightedHit whita,
                              rb::WeightedHit whitb);

}


namespace murem
{

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  MuonRemove::MuonRemove(fhicl::ParameterSet const& pset)
    : fTrkCleanUpAlg(0)
    , fTrkCleanUpPSet(pset.get< fhicl::ParameterSet >("TrkCleanUpPSet") )
    , fNumuEAlg(0)
    , fNumuEAlgPSet(pset.get< fhicl::ParameterSet >("NumuEAlgPSet") )
  {
    this->reconfigure(pset);

    // Produces a vector of muon removed raw digit
    produces< std::vector <rawdata::RawDigit> >();
    produces< std::vector< simb::MCParticle> >();
    produces< art::Assns<simb::MCParticle, rb::Cluster> >();

    // If removing muon by truth, also want to store the mrcc
    // FLSHit list to make BackTracking possible
    if(fMuonID == "removeByTruth"){
      produces< std::vector <sim::FLSHitList>   >();
      produces< std::vector <sim::PhotonSignal> >();
      produces< std::vector <sim::Particle>     >();
    }

    mf::LogInfo("MuonRemove") << " MuonRemove::MuonRemove()\n";
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  MuonRemove::~MuonRemove()
  {
    if(fTrkCleanUpAlg) delete fTrkCleanUpAlg;
    if(fNumuEAlg)      delete fNumuEAlg;
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void MuonRemove::reconfigure(const fhicl::ParameterSet& pset)
  {
    fMCTruthModuleLabel     = pset.get< std::string      >("MCTruthModuleLabel");
    fClusterModuleLabel     = pset.get< std::string      >("ClusterModuleLabel");
    fKalmanTrackModuleLabel = pset.get< std::string      >("KalmanTrackModuleLabel");
    fRawDataLabel           = pset.get< std::string      >("RawDataLabel", "daq");
    fMuonID                 = pset.get< std::string      >("MuonID");
    fNumuReMIdLabel         = pset.get< std::string      >("NumuReMIdLabel");
    fEnergyLabel            = pset.get< std::string      >("EnergyLabel");
    fEnergyLabelBPF         = pset.get< std::string      >("EnergyLabelBPF");
    fElasticInput           = pset.get< std::string      >("ElasticInput");
    fFuzzy3DInput           = pset.get< std::string      >("Fuzzy3DInput");
    fCVNParticleLabel       = pset.get< std::string      >("CVNParticleLabel");
    fBPFTrackModuleLabel    = pset.get< std::string      >("BPFTrackModuleLabel");
    fWhichTrackModule       = pset.get< std::string      >("WhichTrackModule");
    fCVNLabel               = pset.get< std::string      >("CVNLabel");

    fTrackCleanUp       = pset.get< bool             >("TrackCleanUp");
    fADCThreshold_Far   = pset.get< float            >("ADCThreshold_Far");
    fADCThreshold_NDOS  = pset.get< float            >("ADCThreshold_NDOS");
    fADCThreshold_Near  = pset.get< float            >("ADCThreshold_Near");
    fNumuPIDCut         = pset.get< double           >("NumuPIDCut");
    fDCSDistance        = pset.get< int              >("DCSDistance");
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void MuonRemove::beginRun(art::Run& run)
  {

    fTrkCleanUpAlg = new TrackCleanUpAlg(fTrkCleanUpPSet);
    fNumuEAlg      = new numue::NumuEAlg(fNumuEAlgPSet);


    novadaq::cnv::DetId detID = fgeom->DetId();
    if (detID == novadaq::cnv::kNEARDET)
      fADCThreshold   =  fADCThreshold_Near;
    if (detID == novadaq::cnv::kFARDET)
      fADCThreshold   =  fADCThreshold_Far;
    if (detID == novadaq::cnv::kNDOS)
      fADCThreshold   =  fADCThreshold_NDOS;

    return;
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void MuonRemove::produce(art::Event& evt)
  {


    std::unique_ptr< std::vector<rawdata::RawDigit> >
      muonRemovedDigits(new std::vector<rawdata::RawDigit>);
    std::unique_ptr< std::vector<simb::MCParticle > >
      muonInfo(new std::vector< simb::MCParticle>);
    std::unique_ptr< art::Assns<simb::MCParticle, rb::Cluster> >
      assns(new art::Assns<simb::MCParticle, rb::Cluster>);

    std::unique_ptr< std::vector<sim::Particle >    >
      mrccPartList(new std::vector< sim::Particle>);
    std::unique_ptr< std::vector<sim::FLSHitList>   >
      mrccFLSHits(new std::vector<sim::FLSHitList>);
    std::unique_ptr< std::vector<sim::PhotonSignal> >
      mrccPhotSig(new std::vector<sim::PhotonSignal>);

    std::vector< rawdata::RawDigit>                 newDigits;
    std::vector< rb::WeightedHit>                   allTrackHits;
    art::ServiceHandle< calib::Calibrator>          cal;

    art::Handle< std::vector<rawdata::RawDigit> > digitcol;
    evt.getByLabel(fRawDataLabel, digitcol);
    int nRaw = digitcol->size();
    std::list< rawdata::RawDigit >digs;
    for( int iRaw = 0; iRaw < nRaw; ++iRaw ){
      art::Ptr< rawdata::RawDigit > iDig( digitcol, iRaw);
      digs.push_back( *(iDig));
    }// end loop over raw digits

    murem::SortDigitsByPlaneAndCell( digs );
    std::list< rawdata::RawDigit >::iterator digIter, digIterBegin, digIterEnd;

    // Get the slices in the event
    art::Handle< std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fClusterModuleLabel, slicecol);
    art::PtrVector< rb::Cluster> slices;
    for(unsigned int i = 0; i < slicecol->size(); ++i){
        bool filt = true;
        art::Ptr< rb::Cluster> sli(slicecol, i);
        // Use only slices that passed preselection procedure and rockpresel filter
        if(rb::IsFiltered(evt, slicecol, i, "rockpresel")) filt = false;
        bool sel = SliceSelection(evt, sli);
        if (sel && filt) slices.push_back(sli);
    }

    rb::Track muonTrack;

    /// Specify lable which will be used for getting all tracks
    std::string CallSpecificModule = "";
    if (fWhichTrackModule == "kalman") CallSpecificModule = fKalmanTrackModuleLabel;
    if (fWhichTrackModule == "bpf") CallSpecificModule = fBPFTrackModuleLabel;

    art::FindManyP<rb::Track> fmt(slices, evt, CallSpecificModule);

    int nSlices = slices.size();

    // Find the muon track in each slice. Longest track is assumed
    // to be from the Muon for now.
    for( int iSli = 0; iSli < nSlices; iSli++ ){

      // Find all the tracks in kth slice
      std::vector< art::Ptr<rb::Track> > alltracks    = fmt.at(iSli);
      std::vector< art::Ptr<rb::Track> > tracks;
      /// In case of BPF tracks there a re three tracks for each physical track
      /// with different fit assumptions, loop tracks, check the fit sum and choose the muon ones
      if (fWhichTrackModule == "bpf") {
          art::FindManyP<rb::FitSum> fmBPFFitSums(alltracks, evt, fBPFTrackModuleLabel);
          for (unsigned int trackId = 0; trackId < alltracks.size(); ++trackId){
               if (fmBPFFitSums.isValid()){
                   std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId); 
                   if ( abs(bpfFitSums[0]->PDG()) == 13 ) tracks.push_back(alltracks[trackId]);
               }
          }
      }
      /// No special procedures for Kalman tracks
      if (fWhichTrackModule == "kalman")  tracks    =  alltracks;
      
      std::vector< rb::WeightedHit>      trackHits;

      // Select the muon track based on the configuration in fhicl file.
      //Muon energy estimation to be extract from Energy objects. 
      double MuEnergy = -5;

      if(fMuonID == "removeByLength")
        muonTrack = LongestTrack( tracks, evt, MuEnergy );

      if(fMuonID == "removeByEfficiency")
        muonTrack = RemoveByEfficiency( tracks, *slices[iSli], evt, MuEnergy );

      art::ServiceHandle<cheat::BackTracker> bt;
      if(fMuonID == "removeByTruth"){
        if( bt->HaveTruthInfo() )
            trackHits = RemoveByTruth( slices[iSli]);
        else
          std::cerr<<" No MC truth info in this file.\n"
                   <<" Choose a different method for muon removal.\n";
      }

      /// ReMId works only for Kalman tracks
      if ( fMuonID == "removeByReMId" && fWhichTrackModule == "kalman" )
           muonTrack = RemoveByReMId( tracks, evt, MuEnergy );

      // Works only for BPF tracks
      if ( fMuonID == "removeByProngCVN" && fWhichTrackModule == "bpf" ){ 
           muonTrack = RemoveByProngCVN( tracks, slices[iSli], evt, MuEnergy ); 
      }

      if(fMuonID != "removeByTruth" && muonTrack.NCell() > 0){
        muonInfo->push_back( MuonInfo(muonTrack, MuEnergy) );
        util::CreateAssn(*this, evt, *muonInfo, 
                         slices[iSli], *assns);
      }

      if( muonTrack.NCell() != 0 && fMuonID != "removeByTruth" )
        trackHits = fTrkCleanUpAlg->CleanUpTrack(muonTrack, *slices[iSli]);

      allTrackHits.insert( allTrackHits.end(), trackHits.begin(), trackHits.end() );

    }// End loop over slices

    SortByPlaneAndCell( allTrackHits );

    std::vector< rb::WeightedHit >::iterator trkIterStart = allTrackHits.begin(), trkIter;
    for( auto digIter = digs.begin(); digIter != digs.end(); ++digIter){

      rawdata::RawDigit dig = *digIter;
      rb::CellHit cHit(cal->MakeCellHit( &dig ));
      int digPlane = cHit.Plane();
      int digCell  = cHit.Cell();
      int digTDC   = cHit.TDC();

      bool matched = false;

      for( trkIter = trkIterStart; trkIter != allTrackHits.end(); ++trkIter ){

        art::Ptr<rb::CellHit> hit = trkIter->hit;

        if( digPlane < hit->Plane() ){
          trkIterStart = trkIter;
          break;
        }
        else if( digPlane == hit->Plane() &&
                 digCell  <  hit->Cell() ){
          trkIterStart = trkIter;
          break;
        }

        if(digPlane == hit->Plane() &&
           digCell  == hit->Cell()  &&
           digTDC   == hit->TDC() ){

          float muonWeight = trkIter->weight;
          float newADC     = (1-muonWeight)*cHit.ADC();

          if( muonWeight == 1 || newADC < fADCThreshold ){
            trkIterStart = trkIter;
            matched = true;
            break;
          }
          else if(muonWeight > 0 && muonWeight < 1){
            // Need the following steps to be able to handle
            // data/mc from before the multi-point readout era

            if( dig.IsMC() && dig.Version() == 0 && dig.NADC() == 3){
              dig.SetADC(0, newADC );
              dig.SetADC(1, (1-muonWeight)*cHit.ADC(1));
              dig.SetADC(2, (1-muonWeight)*cHit.ADC(2));
            }
            else if( dig.Version() == 0 )
              dig.SetADC(0, newADC );
            else{

              static daqdataformats::NanoSliceVersionConvention conv;
              unsigned int nPretrig = conv.getNPretriggeredSamples( dig.Version());
              int16_t baseline = dig.ADC(nPretrig - fDCSDistance );
              int nadc = dig.NADC();
              for( int iadc = 0; iadc < nadc; ++iadc){
                int16_t newadc = ((dig.ADC(iadc) - baseline )*( 1-muonWeight) ) + baseline;
                dig.SetADC( iadc, newadc);
              }
            }
            matched = true;
            muonRemovedDigits->push_back( dig  );
            trkIterStart = trkIter;
            break;
          }// end if 0 < muonweight < 1    OC
        }
      }// end loop over all track hits
      if( !matched)
        muonRemovedDigits->push_back( *digIter  );
    }// end loop over all digits


    std::cout<<"raw digits before and after removal: "
             <<nRaw<<" "<<muonRemovedDigits->size()<<std::endl;


    if( fMuonID == "removeByTruth"){
      FillTruthInfo(evt, mrccFLSHits, mrccPhotSig, mrccPartList, muonInfo);


      evt.put(std::move(mrccPartList));
      evt.put(std::move(mrccPhotSig));
      evt.put(std::move(mrccFLSHits));

    }

    evt.put(std::move(muonRemovedDigits));
    evt.put(std::move(muonInfo));
    evt.put(std::move(assns));

    return;

  }// end of producer



  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////
  /// Preslection procedure
  bool MuonRemove::SliceSelection (const art::Event& evt, const art::Ptr<rb::Cluster>& slice)
  {

   /// Get slice CVN PID
   /// NB it can cause the bias in the further CVNe/presel ratio
   /*art::FindManyP<cvn::Result> fmCVN ({slice}, evt, fCVNLabel);

   cvn::Result cvnResult;
   float numuid = -1;
   float nueid = -1;
   float nutauid = -1;
   float ncid = -1;
   if (fmCVN.isValid()){
       std::vector<art::Ptr<cvn::Result>> prods = fmCVN.at(0);
       if (!prods.empty()){
           cvnResult = *prods[0];
           numuid = cvnResult.fOutput[0];
           nueid = cvnResult.fOutput[1];
           nutauid = -5;
           ncid = cvnResult.fOutput[2];
       }
   }*/

   /// Loose quality
   /// Vertex exists
   bool vtxValid = false;
   /// Number of prongs > 0
   bool NProngValid = false;
   /// at least one track with E > 0
   double muonEmax = -1;

   art::FindManyP<rb::Vertex> fmElastic ({slice}, evt, fElasticInput);

   if (fmElastic.isValid()){
       std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
       if(elastics.size() > 0) {
          vtxValid = true; /// vertex exists
          /// N Prongs > 0
          art::FindManyP<rb::Prong> fmFuzzyProng3D ( elastics, evt, fFuzzy3DInput);
          std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
          if (prongs.size() > 0) NProngValid = true; /// N prongs > 0 is valid
          /// Check parent muon energy > 0
          /// Check all tracks in this slice
          std::string CallSpecificModule = "";
          if (fWhichTrackModule == "kalman") CallSpecificModule = fKalmanTrackModuleLabel;
          if (fWhichTrackModule == "bpf") CallSpecificModule = fBPFTrackModuleLabel;
          /// Get all tracks
          art::FindManyP<rb::Track> fmt({slice}, evt, CallSpecificModule); 
          std::vector< art::Ptr<rb::Track> > tracks    = fmt.at(0);
          for (unsigned int i = 0; i<tracks.size(); i++){
               /// Use them as they are for Kalman tracks
               if (fWhichTrackModule == "kalman") {
                   art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
                   double muonE = muE.at(0)->E();
                   if (muonE > muonEmax) muonEmax = muonE;
               }       
               if (fWhichTrackModule == "bpf") {
                   /// Check BPF track's fit sum pdg (use only muons)
                   art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[i]}, evt, fBPFTrackModuleLabel);
                   if (fmBPFFitSums.isValid()) {
                       std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
                       if (abs(bpfFitSums[0]->PDG()) == 13) {
                           double muonE = bpfFitSums[0]->FourMom().E();
                           if (muonE > muonEmax) muonEmax = muonE;
                       }
                   }
               }
          } // end of loop over the tracks
       } // end of elastics.size() > 0
   } /// end of vtx is valid 
  /// resulting decision - all these things should exist and CVNm > 0.5 (?) 
  bool ret = vtxValid && NProngValid && (muonEmax > 0); //&& (numuid > 0.5);
 
  return ret;

  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  /// Returns the longest track amongst the given vector of tracks.

  rb::Track MuonRemove::LongestTrack( std::vector< art::Ptr<rb::Track> >& tracks, art::Event& evt, double& muonEnergy)
  {
    rb::Track longestTrack;
    double maxLength = 0;
    double length;
    for(unsigned int i = 0; i < tracks.size(); i++){
        length = tracks[i]->TotalLength();
        if( length > maxLength){
            maxLength = length;
            longestTrack = *tracks[i];
            std::string SpecificEnergyLabel = "";
            /// check the muon energy
            if (fWhichTrackModule == "kalman") {
                SpecificEnergyLabel = fEnergyLabel;
                art::FindOneP<rb::Energy> muE({tracks[i]}, evt, SpecificEnergyLabel);
                muonEnergy = muE.at(0)->E();
            }
            if (fWhichTrackModule == "bpf") {
                art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[i]}, evt, fBPFTrackModuleLabel);
                if (fmBPFFitSums.isValid()) {
                    std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
                    muonEnergy = bpfFitSums[0]->FourMom().E();
                }
            }
        } // end of if
    } // end of loop 

    return longestTrack;
  }


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  /// Finds muon hits by using truth
  /// Loops through all the hits in a given slice and checks if
  /// the photon signal was produced by a muon primary.  If all
  /// the light in a hit is not produced by the muon alone, then
  /// only a fraction of the energy in the hit is attributed to the muon.

  std::vector<rb::WeightedHit> MuonRemove::RemoveByTruth( const art::Ptr<rb::Cluster>& slice)
  {

    std::vector<rb::WeightedHit> muonHits;
    art::ServiceHandle<cheat::BackTracker> bt;

    const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
    int nPrimaries = pnav.NumberOfPrimaries();
    std::vector< int> primMuTrackIDs;

    for(int p = 0; p < nPrimaries; ++p){
      const sim::Particle* prim = pnav.Primary(p);
      if( abs(prim->PdgCode()) == 13 ){

        int muMotherID = prim->TrackId();
        primMuTrackIDs.push_back( muMotherID);

        for( int i = 0; i < (int)pnav.size(); i++ ){

          const sim::Particle* part = pnav.Particle(i);
          if(part->Mother() == muMotherID){
            primMuTrackIDs.push_back( part->TrackId() );
          }
        }//end loop over particles
      }
    }// end loop over primaries

    for(int s = 0; s< (int)slice->NCell(); s++){

      if( !(bt->IsNoise( slice->Cell(s) ) ) ){

        const std::vector< sim::PhotonSignal> pSignal = bt->HitToPhotonSignal( slice->Cell(s));

        int  muonPhotons = 0;
        int  totPhotons  = 0;

        for( int p = 0; p < (int)pSignal.size(); p++){
          totPhotons += pSignal[p].NPhoton();
          int partID = pSignal[p].TrackId();

          // Find the photon signal produced by the muon
          for(int i = 0; i < (int)primMuTrackIDs.size(); i++ ){
            if( partID == primMuTrackIDs[i] ){
              muonPhotons += pSignal[p].NPhoton();
            }
          }
        }

        if(muonPhotons > 0){
          if( muonPhotons == totPhotons ){
            rb::WeightedHit newHit(slice->Cell(s), 1);
            muonHits.push_back(newHit);
          }
          else {
            rb::WeightedHit newHit(slice->Cell(s), (float)muonPhotons/totPhotons);
            muonHits.push_back(newHit);
          }
        }
      }// endif is noise
    }// End loop over slice hits.

    return(muonHits);
  }// End of RemoveByTruth


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  /// Finds the highest efficiency muon track amongst a
  /// vector of tracks. The efficiency is calculated based
  /// on all hits in the slice that the track belongs to,
  /// not on all hits in the event. BackTracker is used
  /// to calculate efficiencies. Since there may be multiple
  /// numuCC interactions in a slice, the process is twofold-
  /// First find the highest efficiency track for every primary
  /// muon, then select the track with highest effciency.


  rb::Track MuonRemove::RemoveByEfficiency( std::vector< art::Ptr<rb::Track> >& tracks,
                                            const rb::Cluster &slice, art::Event& evt, double& muonEnergy)
  {

    rb::Track muonTrack;
    art::ServiceHandle<cheat::BackTracker> bt;

    art::PtrVector< rb::CellHit > allHits   = slice.AllCells();

    const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
    int nPrimaries = pnav.NumberOfPrimaries();
    std::vector< std::set <int> > primMuTrackIDs;

    for(int n = 0; n < nPrimaries; ++n){
      std::set <int> thisMuIDs;
      const sim::Particle* prim = pnav.Primary(n);

      if( abs( prim->PdgCode() ) == 13 ){

        int muMotherID = prim->TrackId();
        thisMuIDs.insert(muMotherID);

        // Find the track IDs of primary muons and its daughters
        for(int p = 0; p < (int)pnav.size(); ++p){
          const sim::Particle* part = pnav.Particle(p);
          //Check if the primary mother is a Muon
          if (part->Mother() == muMotherID )
            thisMuIDs.insert( part->TrackId() );

        }// end loop over particles
      }
      primMuTrackIDs.push_back(thisMuIDs);
    }// end loop over primaries

    double maxEff = -1;
    double eff;

    for( int p = 0; p< (int)primMuTrackIDs.size(); p++){

      // The following loop finds the highest efficiency track for every primary muon.

      for( int i = 0; i< (int)tracks.size(); i++){

        art::PtrVector< rb::CellHit > trackHits = tracks[i]->AllCells();

        eff = bt->HitCollectionEfficiency( primMuTrackIDs[p], trackHits,
                                           allHits,  geo::kXorY);

        if( eff > maxEff){
          maxEff = eff;
          muonTrack = *tracks[i];
          art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
          muonEnergy = muE.at(0)->E();
        }
      }// end loop over tracks
    }// end loop over primary muons

    return muonTrack;
  } // End of RemoveByEfficiency

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////


  /// A function to save the removed muon information in a simb::MCTrajectory
  /// object. Note that only the first trajectory poin is filled put.

  simb::MCParticle MuonRemove::MuonInfo(rb::Track muonTrack, double muonEnergy)
  {

    art::PtrVector< rb::CellHit> muonHits = muonTrack.AllCells();
    rb::SortByPlane(muonHits);

    static TDatabasePDG*  pdgt = TDatabasePDG::Instance();
    TParticlePDG* pdgpm = pdgt->GetParticle(13);
    assert(pdgpm);
    const double muonMass = pdgpm->Mass();


    // trackID, pdgCode, process, motherTrackID, status 1 = stableFinalState

    double muonMom        = TMath::Sqrt((muonEnergy*muonEnergy) - (muonMass*muonMass));

    TLorentzVector startPosition( muonTrack.Start().X(), muonTrack.Start().Y(),
                                  muonTrack.Start().Z(), muonTrack.MinTNS() );

    TLorentzVector startMomentum( muonTrack.Dir().X() * muonMom,
                                  muonTrack.Dir().Y() * muonMom ,
                                  muonTrack.Dir().Z() * muonMom, muonEnergy );

    TLorentzVector stopPosition( muonTrack.Stop().X(),
                                 muonTrack.Stop().Y(),
                                 muonTrack.Stop().Z(),
                                 muonTrack.Cell( muonTrack.NCell() -1 )->TNS() );
    TLorentzVector stopMomentum( 0,0,0,0 );

    simb::MCParticle muonPart( -1000, 13, "primary", -1000, muonMass, 1);
    muonPart.AddTrajectoryPoint( startPosition, startMomentum );
    muonPart.AddTrajectoryPoint( stopPosition,  stopMomentum );
    return muonPart;
  }


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////


  /// A function to write muon removed truth information to the file.
  /// It fills the unique_ptrs for FLSHitList, PhotonSignal, Particles
  /// Also saves the removed muon's information in a vector of simb::MCTrajectory


  void MuonRemove::FillTruthInfo( const art::Event &evt,
                                  std::unique_ptr< std::vector<sim::FLSHitList>     >&     mrccFLSHits,
                                  std::unique_ptr< std::vector<sim::PhotonSignal>   >&     mrccPhotSig,
                                  std::unique_ptr< std::vector<sim::Particle>       >&     mrccPartList,
                                  std::unique_ptr< std::vector<simb::MCParticle>    >&     muonInfo)
  {

    art::ServiceHandle<cheat::BackTracker> bt;

    art::Handle< std::vector<simb::MCTruth> > MCTruthList;
    evt.getByLabel(fMCTruthModuleLabel, MCTruthList);
    if (MCTruthList->empty())
      return;

    const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
    std::vector< int> primMuTrackIDs;

    // Get the TrackId's of Primary Muon and its daughters
    for(int p = 0; p < (int)pnav.NumberOfPrimaries(); ++p){
      const sim::Particle* prim = pnav.Primary(p);

      if( abs(prim->PdgCode()) == 13){

        muonInfo->push_back( *prim);
        int muMotherID = prim->TrackId();
        primMuTrackIDs.push_back(muMotherID);

        for( int i = 0; i < (int)pnav.size(); i++){
          const sim::Particle* part = pnav.Particle(i);
          if( part->Mother() == muMotherID ){
            primMuTrackIDs.push_back( part->TrackId() );
          }
        }// end loop over particles
      }
    }// end loop over primaries

    // Make  FLSHitLists for FLSHits not originating from Primary muon and its daughters.

    art::Handle< std::vector<sim::FLSHitList> > flslistHandle;
    evt.getByLabel("geantgen", flslistHandle);


    for( size_t l = 0; l < flslistHandle->size(); ++l){
      sim::FLSHitList newFLSList;

      const std::vector<sim::FLSHit> &hits = flslistHandle->at(l).fHits;

      for( size_t f = 0; f < hits.size(); ++f){
        bool hitIsFromMuon = false;

        for( int id = 0; id< (int)primMuTrackIDs.size(); id++){
          if(hits[f].GetTrackID() == primMuTrackIDs[id])
            hitIsFromMuon = true;
        }
        if(!hitIsFromMuon)
          newFLSList.fHits.push_back(hits[f]);
      }// end loop over flshits in this flshitlist


      mrccFLSHits->push_back(newFLSList);
    }// end loop over flshit lists


    // Make vector of photon  signals not originating from primary muon and its daughters.

    art::Handle<std::vector<sim::PhotonSignal> > phots;
    evt.getByLabel("photrans", phots);

    for(size_t i = 0; i < phots->size(); ++i){
      const sim::PhotonSignal& phot = (*phots)[i];
      bool sigIsFromMuon = false;

      for(int id = 0; id<(int)primMuTrackIDs.size(); id++){
        if( phot.TrackId() == primMuTrackIDs[id])
          sigIsFromMuon = true;
      }

      if(!sigIsFromMuon)
        mrccPhotSig->push_back(phot);
    } // end of loop over Photon Signals

    // Make particle lists for all things besides muon and its daughters.
    // Also, the muonInfo vector get s filled at this point.

    for (unsigned int s = 0; s < MCTruthList->size(); s++){

      art::Ptr<simb::MCTruth> mc  (MCTruthList, s);
      const art::Ptr<simb::MCTruth> neutrinoInteraction = mc;

      std::vector<const sim::Particle*> particleList = bt->MCTruthToParticles(neutrinoInteraction);

      for (unsigned int i = 0; i < particleList.size(); ++i){
        bool partIsFromMuon = false;
        for(int id = 0; id<(int)primMuTrackIDs.size(); id++){
          if(particleList[i]->TrackId() == primMuTrackIDs[id])
            partIsFromMuon = true;
        }
        if(!partIsFromMuon)
          mrccPartList->push_back(*particleList[i]);
      }
    }// end loop over MCTruth

  }// End of FillTruthInfo

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  rb::Track MuonRemove::RemoveByReMId( std::vector< art::Ptr<rb::Track> >& tracks,
                                       art::Event &evt, double& muonEnergy)
  {
    float maxPidVal = -1;
    rb::Track mostMuLikeTrack;
    art::FindOneP<remid::ReMId> fo(tracks, evt, fNumuReMIdLabel);
    if(fo.isValid()){
      for(int i = 0; i < (int)tracks.size(); i++){
        // Find the ReMId object associated with this track
        art::Ptr<remid::ReMId> pid = fo.at(i);
        float pidVal = pid->Value();
        if( pidVal > maxPidVal ){
          maxPidVal = pidVal;
          mostMuLikeTrack = *tracks[i];
          art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
          muonEnergy = muE.at(0)->E();
        }
      }//end loop over tracks
    }
    else{
      std::cerr<<" There are no remid objects in this file.\n"
               <<" Choose another method of muon removal.\n";
    }
    
    if( maxPidVal > fNumuPIDCut )
      return mostMuLikeTrack;
    else{
      rb::Track voidTrack;
      return voidTrack;
    }
  }// End of ReMId

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  rb::Track MuonRemove::RemoveByProngCVN( std::vector< art::Ptr<rb::Track> >& tracks, const art::Ptr<rb::Cluster>& slice,
                                          art::Event &evt, double& muonEnergy)
  {
     /// Prong CVN works for prongs only, so need to find the most muon like prong first
     /// and then use respective track for the mrcc purposes
     float maxPidVal = -1;
     art::Ptr<rb::Prong> mostMuLikeProng;
     rb::Track mostMuLikeTrack; 
     /// Get the slice's vrtx
     art::FindManyP<rb::Vertex> fmElastic ({slice}, evt, fElasticInput);  
     std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
     if( elastics.size() > 0 ){
         /// Get all prongs for this vtx
         art::FindManyP<rb::Prong> fmFuzzyProng3D ( elastics, evt, fFuzzy3DInput);
         std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
         /// Get the collection of CVNs for these prongs
         art::FindManyP<rb::PID> fmCVNParticleResult(prongs, evt, fCVNParticleLabel);
         for ( unsigned int iPng = 0; iPng < prongs.size(); ++iPng ) {
               art::Ptr<rb::Prong> prong = prongs[iPng];
               std::vector<art::Ptr<rb::PID>> cvnparts = fmCVNParticleResult.at(iPng);
               for ( unsigned int j=0; j<cvnparts.size(); ++j ){
                     if ( abs(cvnparts[j]->Pdg()) == 13 ){
                          double thisPid = cvnparts[j]->Value();
                          if ( thisPid > maxPidVal) {
                               maxPidVal = thisPid;
                               mostMuLikeProng = prongs[iPng];
                          }
                      }
               }
               /// Prong CVN is empty for the prongs with L > 500 cm, assume these are definitely muons
               if ( cvnparts.size() == 0 && fmCVNParticleResult.isValid() && prong->TotalLength() > 500){
                    maxPidVal = 1; ///?
                    mostMuLikeProng = prongs[iPng];
               }
         } // end of loop over the prongs
     }
     /// Done with searhing the most muon like prong 
     /// Find respective tracks (again x3 more tracks in BPF since there are three fit assumptions)
     art::FindManyP<rb::Track> fmMuLikeTracks ( {mostMuLikeProng}, evt, fBPFTrackModuleLabel);
     std::vector< art::Ptr<rb::Track> > MuLikeTracks = fmMuLikeTracks.at(0);
     art::FindManyP<rb::FitSum> fmBPFFitSums ( MuLikeTracks, evt, fBPFTrackModuleLabel );
     /// Loop over all BPF tracks, look for the muon one
     for (unsigned int trackId = 0; trackId < MuLikeTracks.size(); ++trackId){
          if (fmBPFFitSums.isValid()){
              std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId);
              if ( abs(bpfFitSums[0]->PDG()) == 13 ) {
                   muonEnergy = bpfFitSums[0]->FourMom().E();
                   mostMuLikeTrack=*MuLikeTracks[trackId];
              }
          }
     } 
     if ( maxPidVal > fNumuPIDCut )
          return mostMuLikeTrack;
     else{
       rb::Track voidTrack;
       return voidTrack;
     }

  } // End remove by prong CVN 

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  // Compares two raw digits by plane. Same as CompareByPlane
  // function for CellHits in CellHit.h/.cxx. Helper to
  // SortDigitsByPlane.
  bool CompareDigitsByPlaneAndCell( rawdata::RawDigit ra,
                                    rawdata::RawDigit rb)
  {
    art::ServiceHandle<cmap::CMap> cmap;
    if( cmap->GetPlane( &ra ) < cmap->GetPlane( &rb ) )
      return true;
    if( cmap->GetPlane( &ra ) > cmap->GetPlane( &rb ) )
      return false;
    return cmap->GetCell( &ra ) < cmap->GetCell( &rb );

  }// End of CompareByPlane

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  // Sorts a vector of raw digits in order of increasing plane
  // number. Similar to SortByPlane function for CellHits.
  void SortDigitsByPlaneAndCell( std::list< rawdata::RawDigit >& digs)
  {
    digs.sort( CompareDigitsByPlaneAndCell );
  }// End of SortDigitsByPlane


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  bool CompareByPlaneAndCell( rb::WeightedHit whita,
                       rb::WeightedHit whitb)
  {
    if( whita.hit->Plane() < whitb.hit->Plane() )
      return true;
    if( whita.hit->Plane() > whitb.hit->Plane() )
      return false;
    return  whita.hit->Cell() < whitb.hit->Cell();
  }// End of CompareByPlane

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void SortByPlaneAndCell( std::vector< rb::WeightedHit >& whits)
  {
    std::sort( whits.begin(), whits.end(), CompareByPlaneAndCell );
  }// End of SortDigitsByPlane


  DEFINE_ART_MODULE(MuonRemove)

} // End namespace murem

/**********************************************************************/
////////////////////////////////////////////////////////////////////////