BPFCVNAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    BPFCVNAna_module.cc
// \brief   Analyzer module for CVN prongs with BPF goodness 
// \author  psihas@fnal.gov essmith@fnal.gov
////////////////////////////////////////////////////////////////////////

// C/C++ includes
#include <iostream>
#include <string>
#include <sstream>
#include <math.h>


// ROOT includes
#include "TTree.h"
#include "TFile.h"

// Framework includes
#include "art/Framework/Core/FileBlock.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"

// NOvASoft includes
#include "ReMId/ReMId.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FitSum.h"
#include "RecoBase/PID.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Shower.h"
#include "RecoBase/FilterList.h"
#include "Utilities/AssociationUtil.h"
#include "ShowerLID/ShowerLID.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"

#include "BreakPointFitter/DataProducts/BPFPId.h"
#include "BreakPointFitter/func/SliceSummary.h"
#include "CVN/func/TrainingData.h"
#include "CVN/func/InteractionType.h"
#include "CVN/func/ProngType.h"
#include "CVN/func/AssignLabels.h"
#include "CVN/func/Result.h"

#include "NumuEnergy/NumuE.h"
#include "NumuEnergy/NumuEAlg.h"

// -------------------------------------------------------------
//  TO DO
//
//  (fill out my to do list...)
// 
// -------------------------------------------------------------

namespace bpf {

  class BPFCVNAna : public art::EDAnalyzer {
  public:
    explicit BPFCVNAna(fhicl::ParameterSet const& pset);
    ~BPFCVNAna();
    void analyze(const art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void endJob();
    void beginRun(art::Run const&);
    static bool CompareE( const art::Ptr<rb::Prong>& a, const art::Ptr<rb::Prong>& b);

  private:
    const art::ProductToken<std::vector<rb::Cluster>> fClusterToken;

    std::string fVertexLabel;
    std::string fPixelMapInput;
    std::string fProngModLabel;
    std::string fProng3DLabel;
    std::string fProng2DLabel;
    std::string fCVNLabel;
    std::string fTrackLabel;
    std::string fFitSumLabel;
    std::string fBPFPidLabel;
    std::string fEnergyModuleLabel;  
    std::string fCVNSliceLabel;
    std::string fShowerLabel;
    std::string fShowerLIDLabel;
    std::vector<std::string> fPreselectionLabels;
    std::string fCosRejLabel;

    art::ServiceHandle<cheat::BackTracker> bt;
    
    TTree *fCVNBPFslice;
    SliceSummary *fInfo;

    unsigned int Nprongs;
    unsigned int Nprongs2D;

    //TFile* fFile;

    static bool byScore (std::pair<int,float> i,std::pair<int,float> j){
      return ( (i.second > j.second) );
    }
    // Will use to order pairs of < PDG, score >
    static bool scoreCompare(const art::Ptr<rb::PID>& pidA, const art::Ptr<rb::PID>& pidB) {
      return pidA->Value() > pidB->Value();
    }
    
  };


  //.......................................................................
  BPFCVNAna::BPFCVNAna(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
    , fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel")))
//     , fFile(nullptr)
  {
    this->reconfigure(pset);

  }

  //......................................................................
  BPFCVNAna::~BPFCVNAna()
  {
  }

  //......................................................................
  void BPFCVNAna::reconfigure(const fhicl::ParameterSet& pset)
  {
    fPixelMapInput          = pset.get<std::string>("PixelMapInput");
    fVertexLabel            = pset.get<std::string>("VertexLabel");
    fProngModLabel          = pset.get<std::string>("ProngModLabel");
    fProng3DLabel           = pset.get<std::string>("Prong3DLabel");
    fProng2DLabel           = pset.get<std::string>("Prong2DLabel");
    fCVNLabel               = pset.get<std::string>("CVNLabel");
    fCVNSliceLabel          = pset.get<std::string>("CVNSliceLabel");
    fFitSumLabel            = pset.get<std::string>("FitSumLabel");
    fBPFPidLabel            = pset.get<std::string>("BPFPidLabel");
    fTrackLabel             = pset.get<std::string>("TrackLabel");
    fEnergyModuleLabel      = pset.get<std::string>("EnergyModuleLabel");
    fShowerLabel            = pset.get<std::string>("ShowerLabel");
    fShowerLIDLabel         = pset.get<std::string>("ShowerLIDLabel");
    fPreselectionLabels     = pset.get<std::vector<std::string>>("PreselectionLabels");
    fCosRejLabel            = pset.get<std::string>("CosRejLabel");
  }

  //......................................................................
  void BPFCVNAna::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;

    fCVNBPFslice    = tfs->make<TTree>("CVNBPFnu","CVNBPFnu");
    fCVNBPFslice->Branch("fSlice", "bpf::SliceSummary", &fInfo);

  }

  //......................................................................
  void BPFCVNAna::endJob()
  {

  }
  
  //......................................................................
  
  void BPFCVNAna::beginRun(art::Run const&)
  {
  }

  bool BPFCVNAna::CompareE( const art::Ptr<rb::Prong>& a, const art::Ptr<rb::Prong>& b)
  {
    return a->TotalGeV() > b->TotalGeV();
  }

  //......................................................................
  void BPFCVNAna::analyze(const art::Event& evt)
  {

    art::Ptr<simb::MCTruth> truth;
    
    art::Handle< std::vector< rb::Cluster > > sliceHandle;
    evt.getByToken(fClusterToken, sliceHandle);

    const std::vector<rb::Cluster> & clusters = *sliceHandle;
    art::PtrVector<rb::Cluster> slices;
    for ( unsigned int i=0; i<sliceHandle->size(); ++i )
      slices.push_back(art::Ptr<rb::Cluster>(sliceHandle, i));

    art::FindOneP<cvn::Result>    foSliceCVN(sliceHandle, evt, fCVNSliceLabel);
    if ( !foSliceCVN.isValid() ) return;

    // Get Numu Energy objects
    art::FindOneP<numue::NumuE> foNumuE(sliceHandle, evt, fEnergyModuleLabel);
    if ( !foNumuE.isValid() ) return;
   
    // Get MCNeutrino Idx vector matched to the slices
    std::vector<cheat::NeutrinoWithIndex>               nus = bt->allMCTruth();
    std::vector<std::vector<cheat::NeutrinoEffPur>> slTruth = bt->SlicesToMCTruthsTable(slices);
    std::vector<int> sliceNuId = bt->SliceToOrderedNuIds(nus, slTruth, cheat::EnergyMetric, cheat::EnergyMetric);
    art::FindManyP<rb::Track> fmSliceKalman(sliceHandle, evt, "kalmantrackmerge");
    for ( size_t iSlice = 0; iSlice < clusters.size(); ++iSlice ) {

      fInfo = new SliceSummary;

      const rb::Cluster& slice = clusters[iSlice];

      if ( rb::IsFiltered(evt,sliceHandle,iSlice,fPreselectionLabels) ) continue;
      if ( sliceNuId[iSlice] == -1 ) continue;
      if ( slice.IsNoise() ) continue;

      // Fill slice info
      fInfo->slice          = iSlice;
      fInfo->slice_totalGeV = slice.TotalGeV();
      fInfo->slice_nhits    = slice.NCell(geo::kXorY);
      truth = slTruth[iSlice][sliceNuId[iSlice]].neutrinoInt;
      if ( !truth->NeutrinoSet() ) continue;

      // Fill true nu info
      simb::MCNeutrino thisNu = truth->GetNeutrino();
      fInfo->trueE_nu         = thisNu.Nu().E();
      fInfo->trueE_lep        = thisNu.Lepton().E();
      fInfo->trueE_had        = thisNu.Nu().E() -  thisNu.Lepton().E();
      fInfo->truePDG_lep      = thisNu.Lepton().PdgCode();
      fInfo->true_CCNC        = thisNu.CCNC();

      // Fill CVN ID
      art::Ptr<cvn::Result> cvnSliceResult = foSliceCVN.at(iSlice);
      if (!cvnSliceResult) continue;

      // from CAFMaker::FillCVNResultVars
      fInfo->cvnID_NUMU = cvnSliceResult->fOutput[0] + cvnSliceResult->fOutput[1]
                        + cvnSliceResult->fOutput[2] + cvnSliceResult->fOutput[3];
      fInfo->cvnID_NUE  = cvnSliceResult->fOutput[4] + cvnSliceResult->fOutput[5]
                        + cvnSliceResult->fOutput[6] + cvnSliceResult->fOutput[7];

      // Get 3D Prongs
      std::vector<art::Ptr<rb::Prong>> prongs3D;
      art::FindManyP<rb::Prong>        fmProng3D(sliceHandle, evt, 
                                                 art::InputTag(fProngModLabel, fProng3DLabel));
      if ( !fmProng3D.isValid() ) continue;
      if ( fmProng3D.isValid() )  prongs3D = fmProng3D.at(iSlice);

      // Get prongs by energy
      std::sort(prongs3D.begin(),prongs3D.end(), CompareE);
      Nprongs = std::min((int)prongs3D.size(), 20);     

      // Find CVN result and BPF track
      art::FindManyP<rb::PID>    fmCVNID(prongs3D, evt, fCVNLabel);
      art::FindOneP<cvn::Result>   foCVN(prongs3D, evt, fCVNLabel);
      art::FindManyP<rb::Track>  fmTrack(prongs3D, evt, fTrackLabel);
      art::FindOneP<rb::Shower>     foSh(prongs3D, evt, fShowerLabel);

      // Make container to contain all hits
      rb::Cluster hitsMPP; //hits for muon, proton, pions
      rb::Cluster hits3DProngsOther;

      for ( unsigned int iProng = 0; iProng < Nprongs; ++iProng ){              

        if (fmCVNID.isValid()) {
          std::vector< art::Ptr<rb::PID> > cvnPIDs = fmCVNID.at(iProng);
          
          // Get cvn results by PID value
          std::sort(cvnPIDs.begin(), cvnPIDs.end(), scoreCompare);
          for( unsigned int iPID = 0; iPID < cvnPIDs.size(); ++iPID ){
            fInfo->prng_cvnscore[iProng][iPID] = cvnPIDs[iPID]->Value();
            fInfo->prng_cvnpPDG[iProng][iPID]  = cvnPIDs[iPID]->Pdg();
          }
          
        } // end if(fmCVNID.isValid())

        // Get 2D and 3D truth and max purity 
        cvn::ProngType ptype3D, ptypeX, ptypeY;
        double purity3D, purityX, purityY, recE;
        unsigned int ncellX, ncellY;
        bool isprimary;
        
        ptype3D = ProngClassify(*prongs3D[iProng], &ptype3D, &ptypeX, &ptypeY, &isprimary, 
                                &purity3D, &purityX, &purityY, &recE, &ncellX, &ncellY);

        // Get particle contrubutions by purity
        std::vector< std::pair<int, double> > purPDG  = cvn::GetProngPurityByPDG(*prongs3D[iProng]);
        std::vector< std::pair<int, double> > enePDG  = cvn::GetProngEnergyByPDG(*prongs3D[iProng]);
        std::vector< std::pair<int, double> > effPDG  = cvn::GetProngEfficiencyByPDG(*prongs3D[iProng], slice);

        if(purPDG.size()==0) continue;
        for( int iPur = 0; iPur < std::min((int)purPDG.size(),20); ++iPur ){
          fInfo->prng_truePDG[iProng][iPur]  = (int)purPDG[iPur].first;
          fInfo->prng_truepur[iProng][iPur]  = (double)purPDG[iPur].second;
          fInfo->prng_trueeff[iProng][iPur]  = (double)effPDG[iPur].second;
          fInfo->prng_trueE[iProng][iPur]    = (double)enePDG[iPur].second;
        }

        // account for hits that are in prongs that are muon, proton, or pion (by truth)
        // TODO - this, but for prongs that are muon, proton, or pion by CVN...?
        
        if ((int)purPDG[0].first == 13 || (int)purPDG[0].first == 2212 || (int)purPDG[0].first == 211){

          // loop over all hits in prong
          for(unsigned int h = 0; h < prongs3D[iProng]->NCell(); ++h){

            const rb::CellHit *pronghit = prongs3D[iProng]->Cell(h).get();

            //check to see if this hit is already in the list
            bool onlist = false;
            for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
              const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
              if((*pronghit) == (*prongHitInList)) {
                onlist = true;
                break;
              }
            } // end loop over hits in list (hpr)

            // Add the hit if it's not already on the MPP list
            if(!onlist) hitsMPP.Add(prongs3D[iProng]->Cell(h));
          } // end loop over hits 
        }

        // hits that are in 3D prongs that are *not* muon, proton, or pion by truth
        else {
 
          // loop over all hits in prong
          for(unsigned int h = 0; h < prongs3D[iProng]->NCell(); ++h){

            const rb::CellHit *pronghit = prongs3D[iProng]->Cell(h).get();

            // check to see if this hit is already in the MPP list
            bool onlist = false;
            for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
              const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
              if((*pronghit) == (*prongHitInList)) {
                onlist = true;
                break;
              }
            } // end loop over hits in MPP list (hpr)

            // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
            if (!onlist){
              for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
                const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
                if((*pronghit) == (*prongHitInList)) {
                  onlist = true;
                  break;
                }
              }
            }

            // Add the hit if its not already on the "other" list
            if(!onlist) hits3DProngsOther.Add(prongs3D[iProng]->Cell(h));
          } // end loop over hits 
        } //end else for non-{muon, proton, pion} 

        // get  LID info
        if ( foSh.isValid() ){
          cet::maybe_ref<art::Ptr<rb::Shower> const> roSh(foSh.at(iProng));
          art::Ptr<rb::Shower> shower = roSh.ref();

          std::vector<art::Ptr<rb::Shower>> showers;
          showers.push_back(shower);

          art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
          if ( foShLID.isValid() ){
            cet::maybe_ref<art::Ptr<slid::ShowerLID> const> roShLID(foShLID.at(0));
            art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
            if ( roShLID.ref() ) {
              fInfo->shwLID_ePID[iProng] = showerLID->Value();
              fInfo->shwLID_shwE[iProng] = showerLID->ShowerEnergy();

              // Get particle contrubutions by purity
              std::vector< std::pair<int,double> > purPDGs=cvn::GetProngPurityByPDG(*showers[0]);
              std::vector< std::pair<int,double> > enePDGs=cvn::GetProngEnergyByPDG(*showers[0]);
              std::vector< std::pair<int,double> > effPDGs=cvn::GetProngEfficiencyByPDG(*showers[0], slice);
            
              if(purPDGs.size()==0) continue;
              for( int iPur = 0; iPur < std::min((int)purPDGs.size(),20); ++iPur ){
                fInfo->shwLID_truePDG[iProng][iPur] = (int)purPDGs[iPur].first;
                fInfo->shwLID_truepur[iProng][iPur] = (double)purPDGs[iPur].second;
                fInfo->shwLID_trueeff[iProng][iPur] = (double)effPDGs[iPur].second;
              }    
            }
          }
        } // end shower LID

        art::Ptr<rb::Prong> prong = prongs3D[iProng];
        // Get the tracks associated with this prong
        std::vector< art::Ptr< rb::Track > > tracks;
        if (!fmTrack.isValid()) continue;
        tracks = fmTrack.at(iProng);

        art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
        // loop over BPF tracks
        art::FindManyP<bpfit::BPFPId> fmBPFTrackBPFPid(tracks, evt, fBPFPidLabel);
        for ( unsigned int itrack = 0; itrack < tracks.size(); ++itrack ){
          std::vector< art::Ptr< rb::FitSum > > fitsums;
          if ( fitsumFMP.isValid() ) fitsums = fitsumFMP.at(itrack);
          if ( fitsums.size() != 1 ) continue;

          art::Ptr< rb::Track > track = tracks.at(itrack);

          double phi = atan2(track->Dir().Y(), track->Dir().X())*180.0/M_PI;
          if(phi<0.0) phi += 360.0;

          fInfo->trkBPF_PDG[iProng][itrack] = fitsums[0]->PDG();
          fInfo->trkBPF_E[iProng][itrack]   = fitsums[0]->FourMom().E();
          fInfo->trkBPF_len[iProng][itrack] = track->TotalLength();
          fInfo->trkBPF_theta[iProng][itrack]       = acos(track->Dir().Z())*180.0/M_PI;
          fInfo->trkBPF_phi[iProng][itrack]         = phi;

          std::vector< std::pair<int, double> > purPDG = cvn::GetProngPurityByPDG(*(tracks[itrack]));
          std::vector< std::pair<int, double> > effPDG = cvn::GetProngEfficiencyByPDG(*(tracks[itrack]), slice);

          // loop over particle contributors to this track
          if(purPDG.size() == 0) continue;
          for(int iMatch = 0; iMatch < std::min((int)purPDG.size(),20); ++iMatch) {

            fInfo->trkBPF_truePDG[iProng][itrack][iMatch] = (int)purPDG[iMatch].first;
            fInfo->trkBPF_truepur[iProng][itrack][iMatch] = (double)purPDG[iMatch].second;
            fInfo->trkBPF_trueeff[iProng][itrack][iMatch] = (double)effPDG[iMatch].second;          

          } // end loop over matches (iMatch)

          // get the BPFpid value for this track...
          if(fmBPFTrackBPFPid.isValid()) {
            std::vector<art::Ptr<bpfit::BPFPId>> bpfpids = fmBPFTrackBPFPid.at(itrack);
            if(bpfpids.size() < 1) continue;
            // NOTE: There should only be one PID associated with each track. PIDs for different particle assumptions
            //       are expected to be associated with the appropriate BPF track
            fInfo->trkBPF_pID[iProng][itrack] = bpfpids[0]->Value();
          }
        }//bpf tracks     

        // Fill other prong information
        fInfo->prng_trueIDx[iProng]    = ptypeX;
        fInfo->prng_trueIDy[iProng]    = ptypeY;
        fInfo->prng_trueID[iProng]     = ptype3D;
        fInfo->prng_nhitsX[iProng]     = ncellX;
        fInfo->prng_nhitsY[iProng]     = ncellY;
        fInfo->prng_nhits3D[iProng]    = ncellX+ncellY;
        fInfo->prng_purX[iProng]       = purityX;
        fInfo->prng_purY[iProng]       = purityY;
        fInfo->prng_pur3D[iProng]      = purity3D;
        fInfo->prng_length[iProng]     = prong->TotalLength();
        fInfo->prng_totalGeV[iProng]   = prong->TotalGeV();     
        fInfo->prng_missingpln[iProng] = prong->NMissingPlanes(geo::kXorY);

      }// over prongs

      // Get nu, muon, hadronic energies as calculated by SA numu analysis       
      const art::Ptr<numue::NumuE> sliceEnergy = foNumuE.at(iSlice); 
      if ( sliceEnergy ){
        fInfo->numuE_nu     = sliceEnergy->E();
        fInfo->numuE_lep    = sliceEnergy->RecoMuonE();
        fInfo->numuE_had    = sliceEnergy->RecoTrkCCHadE();
      }
      
      // Fill info for Kalman tracks
      std::vector< art::Ptr<rb::Track> > sliceKalmanTracks;
      if ( fmSliceKalman.isValid() ){

        sliceKalmanTracks = fmSliceKalman.at(iSlice);
        art::FindOneP<remid::ReMId> foKalmanRemid(sliceKalmanTracks, evt, "remid");
        // loop over all tracks
        for(int iTrk = 0; iTrk < std::min((int)sliceKalmanTracks.size(),20); ++iTrk) {
          art::Ptr< rb::Track > track = sliceKalmanTracks.at(iTrk);

          std::vector< std::pair<int, double> > purPDG = cvn::GetProngPurityByPDG(*track);
          std::vector< std::pair<int, double> > effPDG = cvn::GetProngEfficiencyByPDG(*track, slice);

          double phi = atan2(track->Dir().Y(), track->Dir().X())*180.0/M_PI;
          if(phi<0.0) phi += 360.0;

          fInfo->trkKal_len[iTrk]   = track->TotalLength();
          fInfo->trkKal_theta[iTrk] = acos(track->Dir().Z())*180.0/M_PI;
          fInfo->trkKal_phi[iTrk]   = phi;

          // loop over particle contributors to this track
          if(purPDG.size() == 0) continue; // don't need this line
          for(int iMatch = 0; iMatch < std::min((int)purPDG.size(),20); ++iMatch) {

            fInfo->trkKal_truePDG[iTrk][iMatch] = (int)purPDG[iMatch].first;
            fInfo->trkKal_truepur[iTrk][iMatch] = (double)purPDG[iMatch].second;
            fInfo->trkKal_trueeff[iTrk][iMatch] = (double)effPDG[iMatch].second;            

          } // end loop over matches (iMatch)

          // Get the remid values for each track and fill...
          if(foKalmanRemid.isValid()) {
            fInfo->trkKal_pID[iTrk] = (foKalmanRemid.at(iTrk))->Value();
          }

        } // end loop over tracks (iTrk)
        
      } //kalmantrk

      // Get 2D Prongs
      std::vector<art::Ptr<rb::Prong>> prongs2D;
      art::FindManyP<rb::Prong>        fmProng2D(sliceHandle, evt,
                                                 art::InputTag(fProngModLabel, fProng2DLabel));

      rb::Cluster hits2DProngs;
      if ( fmProng2D.isValid() ) {
        prongs2D = fmProng2D.at(iSlice);

        std::sort(prongs2D.begin(),prongs2D.end(),CompareE);
        Nprongs2D = std::min((int)prongs2D.size(), 20);

        for (unsigned int iProng = 0; iProng < Nprongs2D; ++iProng){
          std::vector< std::pair<int, double> > purPDG = cvn::Get2DProngPurityByPDG(*prongs2D[iProng]);
          std::vector< std::pair<int, double> > enePDG = cvn::Get2DProngEnergyByPDG(*prongs2D[iProng]);
          std::vector< std::pair<int, double> > effPDG = cvn::Get2DProngEfficiencyByPDG(*prongs2D[iProng], slice); 

          if(purPDG.size()==0) continue; 
          for (int iPur = 0; iPur < std::min((int)purPDG.size(),20); ++iPur ){
            fInfo->prng2D_truePDG[iProng][iPur]  = (int)purPDG[iPur].first;
            fInfo->prng2D_truepur[iProng][iPur]  = (double)purPDG[iPur].second;
            fInfo->prng2D_trueeff[iProng][iPur]  = (double)effPDG[iPur].second;
            fInfo->prng2D_trueE[iProng][iPur]    = (double)enePDG[iPur].second;
          }

          for(unsigned int h2D = 0; h2D < prongs2D[iProng]->NCell(); ++ h2D){
            const rb::CellHit *pronghit = prongs2D[iProng]->Cell(h2D).get();
            // check to see if this hit is already in the 3D lists
            bool onlist = false;
            for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
              const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
              if((*pronghit) == (*prongHitInList)) {
                onlist = true;
                break;
              }
            } // end loop over hits in MPP list (hpr)

            // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
            if (!onlist){
              for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
                const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
                if((*pronghit) == (*prongHitInList)) {
                  onlist = true;
                  break;
                }
              }
            }

            // if not on either of those lists, check to see if it's already on 2D prong list
            if (!onlist){
              for(unsigned int hpr = 0; hpr < hits2DProngs.NCell(); ++hpr) {
                const rb::CellHit *prongHitInList = hits2DProngs.Cell(hpr).get();
                if((*pronghit) == (*prongHitInList)) {
                  onlist = true;
                  break;
                }
              }
            }
            if (!onlist) hits2DProngs.Add(prongs2D[iProng]->Cell(h2D));
          }//end loop over hits

          art::Ptr<rb::Prong> prong2D = prongs2D[iProng];
          // Misc 2D prong info
          fInfo->prng2D_length[iProng]     = prong2D->TotalLength();
          fInfo->prng2D_totalGeV[iProng]   = prong2D->TotalGeV();       
          fInfo->prng2D_missingpln[iProng] = prong2D->NMissingPlanes(geo::kXorY);
          fInfo->prng2D_nhits[iProng]      = prong2D->NCell(geo::kXorY);
        }//end 2D prong loop
      }// end 2D prongs isValid

      // orphaned hits - adapted from BPFEnergyEstimator
      rb::Cluster orphanedHits;
      int nOrphHits = 0;
      // Loop over all hits in the slice
      for (unsigned int hsl = 0; hsl < slice.NCell(); ++hsl) {
        const rb::CellHit *sliceHit = slice.Cell(hsl).get();

        // check to see if this hit is in any of the lists
        bool onlist = false;
        for(unsigned int hpr = 0; hpr < hitsMPP.NCell(); ++hpr) {
          const rb::CellHit *prongHitInList = hitsMPP.Cell(hpr).get();
          if((*sliceHit) == (*prongHitInList)) {
            onlist = true;
            break;
          }
        } // end loop over hits in MPP list (hpr)

        // if not already in the MPP list, check to see if this is already in the "other" 3D prong list
        if (!onlist){
          for(unsigned int hpr = 0; hpr < hits3DProngsOther.NCell(); ++hpr) {
            const rb::CellHit *prongHitInList = hits3DProngsOther.Cell(hpr).get();
            if((*sliceHit) == (*prongHitInList)) {
              onlist = true;
              break;
            }
          }
        }

        // if not on either of those lists, check to see if it's already on 2D prong list
        if (!onlist){
          for(unsigned int hpr = 0; hpr < hits2DProngs.NCell(); ++hpr) {
            const rb::CellHit *prongHitInList = hits2DProngs.Cell(hpr).get();
            if((*sliceHit) == (*prongHitInList)) {
              onlist = true;
              break;
            }
          }
        }
        // if it's not on a 3D prong or a 2D prong, put it in the list of orphaned hits
        if(!onlist) {
          orphanedHits.Add(slice.Cell(hsl));
          nOrphHits++;
        }
      }
      fInfo->orph_nhits = nOrphHits;
      if(orphanedHits.NCell()>0) fInfo->orph_totalGeV = orphanedHits.TotalGeV();

      // for cosmic rejection cuts
      art::FindManyP<cosrej::CosRejObj> fmCosRej(sliceHandle,evt,fCosRejLabel);
      if (fmCosRej.isValid()) {
        const std::vector<art::Ptr<cosrej::CosRejObj> > cosrej = fmCosRej.at(iSlice);
        if (cosrej.size() > 0){
          fInfo->anglekal  = cosrej[0]->AngleKal();
          fInfo->concospid = cosrej[0]->ConCosPID();
        }
        
      }
 
      fInfo->Nprongs   = Nprongs;
      fInfo->Nprongs2D = Nprongs2D;
      fInfo->Npids     = 10; //arbitrarily larger than expected contributions to one prong
      fInfo->run       = evt.run();
      fInfo->subrun    = evt.subRun();
      fInfo->event     = evt.id().event();

      fCVNBPFslice->Fill();
    }// slices
    
  }

  //----------------------------------------------------------------------
  DEFINE_ART_MODULE(bpf::BPFCVNAna)
} // end namespace bpf
///////////////////////////////////////////////////////////////////