SkimmingUtils.h

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
// \version $ $
//
// \brief Utility object to make skimming easier
//
// \author brebel@fnal.gov
//
////////////////////////////////////////////////////////////////////////////

#ifndef SKIMMINGUTILS_H
#define SKIMMINGUTILS_H

#include <map>

#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Assns.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/EDFilter.h"
#include "canvas/Persistency/Provenance/ProductID.h" 
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindOne.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Utilities/AssociationUtil.h"
#include "GeometryObjects/CellUniqueId.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Track.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Vertex.h"
#include "ShowerLID/ShowerLID.h"
#include "MCCheater/BackTracker.h"
#include "MCReweight/FluxWeights.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/GTruth.h"

namespace skim {
  class SkimmingUtils;
}

namespace skim  {
  class SkimmingUtils {
    
  public:
    
    SkimmingUtils();
    SkimmingUtils(art::Handle<std::vector<rb::CellHit> > const& cellHandle,
                  std::string                            const& instance=std::string());
    virtual ~SkimmingUtils();
    
    // method that fills a map of std::vector< T const* > to an index for a slice
    template<class T> void FillSliceMap(size_t                                  const& s,
                                        art::FindMany<T>                             & fmp,
                                        std::map<size_t, std::vector<T const*> >     & ptrMap);

    // method that fills a std::vector< T const* > for a slice
    template<class T> void FillSliceVector(size_t           const& s,
                                           art::FindMany<T>      & fmp,
                                           std::vector<T const*> & ptrMap);

    // method that fills a std::vector< T const* > for a slice based on another type
    // U where there is a one to one association between T and U.  Use this sparingly
    // as this method creates a FindOne which is expensive - can't seem to find a
    // better way of doing it though
    template<class T, class U> void FillSliceVector(art::Event       const& e,
                                                    std::string      const& label,
                                                    std::vector<U const*> & otherPtrs,
                                                    std::vector<T const*> & fillPtrs);

    template<class T, class U> void FillSliceVector(art::Event       const& e,
                                                    art::InputTag    const& tag,
                                                    std::vector<U const*> & otherPtrs,
                                                    std::vector<T const*> & fillPtrs);

    
    // method that copies the requested data products and creates an association between the
    // product and the slice from which it came.  The data product must not contain any kind
    // of art::Ptrs or similar external references because we are using a copy constructor
    // and those references will not copy correctly
    template<class T> bool CopyProductAndSliceAssn(art::EDProducer           const& prod,
                                                   art::Event                     & e,
                                                   std::vector<art::Ptr<T> > const& objIn,
                                                   std::vector<T>                     & objOut,
                                                   std::vector<rb::Cluster>       & sliceOut,
                                                   art::Assns<rb::Cluster, T>     & slcCROAssn);
    
    bool CopySlice(art::Ptr<rb::Cluster>     const& sliceIn,
                   std::vector<rb::Cluster>       & sliceOut,
                   size_t                         & sliceNum);
    
    bool CopyVertex(art::EDProducer                     const& prod,
                    art::Event                               & e,
                    art::Ptr<rb::Vertex>                const& vertexIn,
                    std::vector<rb::Vertex>                  & vertexOut,
                    std::vector<rb::Cluster>                 & sliceOut,
                    art::Assns<rb::Cluster, rb::Vertex>      & slcVtxAssn);
    
    bool CopyShowers(art::EDProducer                      const& prod,
                     art::Event                                & e,
                     std::vector< art::Ptr<rb::Shower> >        const& showerIn,
                     std::vector<rb::Shower>                         & showerOut,
                     std::vector<rb::Cluster>                        & sliceOut,
                     art::Assns<rb::Cluster, rb::Shower>             & slcShwAssn);
    
    bool CopyTracks(art::EDProducer                      const& prod,
                    art::Event                                & e,
                    std::vector< art::Ptr<rb::Track> >   const& trackIn,
                    std::vector<rb::Track>                    & trackOut,
                    std::vector<rb::Cluster>                  & sliceOut,
                    art::Assns<rb::Cluster, rb::Track>        & slcTrkAssn,
                    std::string                          const& slcInstance="",
                    std::string                          const& trkInstance="");
    
    bool CopyProngs(art::EDProducer                      const& prod,
                     art::Event                               & e,
                     std::vector< art::Ptr<rb::Prong> >  const& prongIn,
                     std::vector<rb::Prong>                          & prongOut,
                     std::vector<rb::Cluster>                       & sliceOut,
                     art::Assns<rb::Cluster, rb::Prong>              & slcPrgAssn);

  bool SliceMCTruthAssociation(art::EDProducer                        const& prod,
                               art::Event                                  & e,
                               size_t                                 const& slcIdx,
                               std::vector<rb::Cluster>                    & sliceOut,
                               art::Assns<rb::Cluster, simb::MCTruth>      & slcTruAssn,
                               art::Assns<rb::Cluster, simb::MCFlux>       & slcFlxAssn,
                               art::Assns<rb::Cluster, fxwgt::FluxWeights> & slcFlxWgtAssn,
                               art::Assns<rb::Cluster, simb::GTruth>       & slcGTruAssn,
                               std::string                            const& generatorLabel,
                               std::string                            const& fluxWeightLabel);
  
  private:
  
    std::map<std::pair<geo::CellUniqueId, float>, art::Ptr<rb::CellHit> > fCellIdToHitPtr;
    art::ServiceHandle<geo::Geometry>                                     fGeom;
    std::string                                                           fInstance;
  };
  
} // end namespace

//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::CopySlice(art::Ptr<rb::Cluster>     const& sliceIn,
                                           std::vector<rb::Cluster>       & sliceOut,
                                           size_t                         & sliceNum)
{
  if(fCellIdToHitPtr.size() < 1)
    throw cet::exception("SkimmingUtils")
    << "CellId/TNS to art::Ptr<rb::CellHit> map is empty, no way to copy the slice";
  
  art::ServiceHandle<geo::Geometry> geom;
  // make a art::PtrVector<rb::CellHit> for all the hits in the slice
  auto allCells = sliceIn->AllCells();

  std::set<std::pair<geo::CellUniqueId, float> > cuid;
  
  art::PtrVector<rb::CellHit> slcCells;
  std::vector<double> slcWeights;

  // figure out which cells belong in the slice and the weight of each
  for(size_t ch = 0; ch < allCells.size(); ++ch){
    auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
                              allCells[ch]->TNS());
    
    // don't duplicate cells in the output list
    if( !cuid.insert(key).second ){
      LOG_WARNING("SkimmingUtils")
      << "trying to insert a duplicate cell " << (*allCells[ch]);
      continue;
    }
    
    if( fCellIdToHitPtr.count(key) > 0 ){
      slcCells.push_back( fCellIdToHitPtr.find(key)->second );
      slcWeights.emplace_back( sliceIn->Weight(ch) );
    }
  }

  // make the new slice on the output collection
  sliceOut.emplace_back(sliceIn->View(), sliceNum);
  sliceOut.back().Add(slcCells, slcWeights);

  return true;
}

//------------------------------------------------------------------------------
template<class T> inline void skim::SkimmingUtils::FillSliceMap(size_t                                 const& s,
                                                                art::FindMany<T>                            & fmp,
                                                                std::map<size_t, std::vector<T const*> >    & ptrMap)
{  
  std::vector<T const*> ptrs;
  if( fmp.isValid() ){
    fmp.get(s, ptrs);
    ptrMap[s] = ptrs;
  }

  return;
}

//------------------------------------------------------------------------------
template<class T, class U> inline void skim::SkimmingUtils::FillSliceVector(art::Event       const& e,
                                                                            std::string      const& label,
                                                                            std::vector<U const*> & otherPtrs,
                                                                            std::vector<T const*> & fillPtrs)
{
  // make the FindOne
  art::InputTag tag(label);

  this->FillSliceVector(e, tag, otherPtrs, fillPtrs);
  
  return;
}

//------------------------------------------------------------------------------
template<class T, class U> inline void skim::SkimmingUtils::FillSliceVector(art::Event       const& e,
                                                                            art::InputTag    const& tag,
                                                                            std::vector<U const*> & otherPtrs,
                                                                            std::vector<T const*> & fillPtrs)
{
  
  // make the FindOne
  art::FindOne<T> fo(otherPtrs, e, tag);
  
  //clear the input vector
  fillPtrs.clear();
  fillPtrs.reserve(otherPtrs.size());
                   
  // loop over the thesePtrs and get the associated others
  for(size_t t = 0; t < otherPtrs.size(); ++t){
    if(fo.at(t).isValid()) fillPtrs.emplace_back(&(fo.at(t).ref()));
  }

  fillPtrs.shrink_to_fit();
                   
  return;
}


//------------------------------------------------------------------------------
template<class T> inline void skim::SkimmingUtils::FillSliceVector(size_t           const& s,
                                                                   art::FindMany<T>      & fmp,
                                                                   std::vector<T const*> & ptrVec)
{
  if( fmp.isValid() ) fmp.get(s, ptrVec);

  return;
}

//------------------------------------------------------------------------------
template<class T> inline bool skim::SkimmingUtils::CopyProductAndSliceAssn(art::EDProducer           const& prod,
                                                                           art::Event                     & e,
                                                                           std::vector<art::Ptr<T> > const& objIn,
                                                                           std::vector<T>                     & objOut,
                                                                           std::vector<rb::Cluster>       & sliceOut,
                                                                           art::Assns<rb::Cluster, T>     & slcObjAssn)
{

  for(auto in : objIn){
    objOut.push_back( *in );
    util::CreateAssn(prod,
                     e,
                     sliceOut,
                     objOut,
                     slcObjAssn,
                     objOut.size()-1,
                     objOut.size(),
                     UINT_MAX,
                     fInstance,
                     fInstance);
 }    

  return true;
}
    
//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::CopyVertex(art::EDProducer                     const& prod,
                                            art::Event                               & e,
                                            art::Ptr<rb::Vertex>                const& vertexIn,
                                            std::vector<rb::Vertex>                  & vertexOut,
                                            std::vector<rb::Cluster>                 & sliceOut,
                                            art::Assns<rb::Cluster, rb::Vertex>      & slcVtxAssn)
{
  // I am assured there is only one vertex per slice by the elastic arms module
  vertexOut.emplace_back(vertexIn->GetXYZ(), vertexIn->GetT() );
  util::CreateAssn(prod,
                   e,
                   sliceOut,
                   vertexOut,
                   slcVtxAssn,
                   vertexOut.size()-1,
                   vertexOut.size(),
                   UINT_MAX,
                   fInstance,
                   fInstance);
  
  return true;
}

//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::CopyShowers(art::EDProducer                     const& prod,
                                             art::Event                               & e,
                                             std::vector< art::Ptr<rb::Shower> > const& showerIn,
                                             std::vector<rb::Shower>                  & showerOut,
                                             std::vector<rb::Cluster>                 & sliceOut,
                                             art::Assns<rb::Cluster, rb::Shower>      & slcShwAssn)
{
  if(fCellIdToHitPtr.size() < 1)
    throw cet::exception("SkimmingUtils")
    << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the showers";

  // loop over all the showers
  art::ServiceHandle<geo::Geometry> geom;
  art::PtrVector<rb::CellHit> allCells;
  for(auto shw : showerIn){
    
    allCells = shw->AllCells();
    
    std::vector<double> shwCellWt(shw->NCell(), 1.);
    art::PtrVector<rb::CellHit> shwHitPtrs;
    
    // loop over all the hits in the shower to get their weights for the 
    // shower and to determine which of the Ptrs of the new slice collection
    // are needed
    for(size_t ch = 0; ch < allCells.size(); ++ch){
      auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
                                allCells[ch]->TNS());

      if( fCellIdToHitPtr.count(key) > 0 ){
        shwHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
        shwCellWt[ch] = shw->Weight(ch);
      }
    } // end loop over cells in the current shower

    // make the copy of the shower
    showerOut.emplace_back(shwHitPtrs,
                           shwCellWt,
                           shw->TotalLength(),
                           shw->Start(),
                           shw->Dir(),
                           shw->ID() );

    util::CreateAssn(prod,
                     e,
                     sliceOut,
                     showerOut,
                     slcShwAssn,
                     showerOut.size()-1,
                     showerOut.size(),
                     UINT_MAX,
                     fInstance,
                     fInstance);

  } // end loop over showers

  return true;
}

//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::CopyTracks(art::EDProducer                    const& prod,
                                            art::Event                              & e,
                                            std::vector< art::Ptr<rb::Track> > const& trackIn,
                                            std::vector<rb::Track>                  & trackOut,
                                            std::vector<rb::Cluster>                & sliceOut,
                                            art::Assns<rb::Cluster, rb::Track>      & slcTrkAssn,
                                            std::string                        const& slcInstance,
                                            std::string                        const& trkInstance)
{
  
  std::string slcInst = (slcInstance.empty()) ? fInstance : slcInstance;
  std::string trkInst = (trkInstance.empty()) ? fInstance : trkInstance;
  
  if(fCellIdToHitPtr.size() < 1)
    throw cet::exception("SkimmingUtils")
    << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the tracks";

  // loop over all the showers
  art::ServiceHandle<geo::Geometry> geom;
  art::PtrVector<rb::CellHit> allCells;
  for(auto trk : trackIn){

    //////
    //    if( trk->View() != geo::kXorY) continue; // Copy Only 3D Tracks
    /////    
    allCells = trk->AllCells();
    
    std::vector<double> trkCellWt(trk->NCell(), 1.);
    art::PtrVector<rb::CellHit> trkHitPtrs;
    
    // loop over all the hits in the shower to get their weights for the 
    // shower and to determine which of the Ptrs of the new slice collection
    // are needed
    for(size_t ch = 0; ch < allCells.size(); ++ch){
      auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
                                allCells[ch]->TNS());

      if( fCellIdToHitPtr.count(key) > 0 ){
        trkHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
        trkCellWt[ch] = trk->Weight(ch);
      }
    } // end loop over cells in the current shower

    // make the copy of the track
    trackOut.emplace_back( trkHitPtrs,
                          trkCellWt,
                          trk->Trajectory(),
                          trk->Start(),
                          trk->Dir(),
                          trk->ID(),
                          trk->View());

    util::CreateAssn(prod,
                     e,
                     sliceOut,
                     trackOut,
                     slcTrkAssn,
                     trackOut.size()-1,
                     trackOut.size(),
                     UINT_MAX,
                     slcInst,
                     trkInst);

  } // end loop over tracks

  return true;
}


//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::CopyProngs(art::EDProducer                    const& prod,
                                            art::Event                              & e,
                                            std::vector< art::Ptr<rb::Prong> > const& prongIn,
                                            std::vector<rb::Prong>                  & prongOut,
                                            std::vector<rb::Cluster>                & sliceOut,
                                            art::Assns<rb::Cluster, rb::Prong>      & slcPrgAssn)
{
  if(fCellIdToHitPtr.size() < 1)
    throw cet::exception("SkimmingUtils")
    << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the prongs";
  
  // loop over all the prongs
  art::ServiceHandle<geo::Geometry> geom;
  art::PtrVector<rb::CellHit> allCells;
  for(auto png : prongIn){
    
    allCells = png->AllCells();
    
    std::vector<double> pngCellWt(png->NCell(), 1.);
    art::PtrVector<rb::CellHit> pngHitPtrs;
    
    // loop over all the hits in the prong to get their weights for the
    // prong and to determine which of the Ptrs of the new slice collection
    // are needed
    for(size_t ch = 0; ch < allCells.size(); ++ch){
      auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
                                allCells[ch]->TNS());
      
      if( fCellIdToHitPtr.count(key) > 0 ){
        pngHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
        pngCellWt[ch] = png->Weight(ch);
      }
    } // end loop over cells in the current prong
    
    // make the copy of the prong
    prongOut.emplace_back(pngHitPtrs,
                          png->Start(),
                          png->Dir(),
                          png->ID() );
    
    util::CreateAssn(prod,
                     e,
                     sliceOut,
                     prongOut,
                     slcPrgAssn,
                     prongOut.size()-1,
                     prongOut.size(),
                     UINT_MAX,
                     fInstance,
                     fInstance);
    
  } // end loop over showers
  
  return true;
}

//------------------------------------------------------------------------------
inline bool skim::SkimmingUtils::SliceMCTruthAssociation(art::EDProducer                        const& prod,
                                                         art::Event                                  & e,
                                                         size_t                                 const& slcIdx,
                                                         std::vector<rb::Cluster>                    & sliceOut,
                                                         art::Assns<rb::Cluster, simb::MCTruth>      & slcTruAssn,
                                                         art::Assns<rb::Cluster, simb::MCFlux>       & slcFlxAssn,
                                                         art::Assns<rb::Cluster, fxwgt::FluxWeights> & slcFlxWgtAssn,
                                                         art::Assns<rb::Cluster, simb::GTruth>       & slcGTruAssn,
                                                         std::string                            const& generatorLabel,
                                                         std::string                            const& fluxWeightLabel)
{
  
  // have to get the BackTracker service to be able to make the truth to slice
  // association
  art::ServiceHandle<cheat::BackTracker> bt;
  
  // We want to get the slice truth, and to match the CAF analysis we use the
  // energy matching metric
  std::vector<rb::Cluster> slice({sliceOut[slcIdx]});
  auto sliceTruthTable = bt->SlicesToMCTruthsTable(slice);
  
  if(sliceTruthTable[0].size() < 1){
    LOG_WARNING("SkimmingUtils")
    << "Could not find MCTruth objects for this slice in backtracker, so "
    << "return without setting MCTruth or MCFlux associations for this slice";
    return false;
  }
  
  auto nuIdx = bt->SliceToOrderedNuIds(bt->allMCTruth(),
                                       sliceTruthTable,
                                       cheat::EnergyMetric,
                                       cheat::EnergyMetric);

  // the first entry in nuIdx is supposed to be the best match of this slice
  // to a true interaction there is only one slice in sliceTruthTable because
  // we only passed one in the vector that made it
  auto trueInt = sliceTruthTable[0][nuIdx[0]].neutrinoInt;
  
  util::CreateAssn(prod,
                   e,
                   sliceOut,
                   trueInt,
                   slcTruAssn,
                   slcIdx,
                   fInstance);
  
  // find the MCTruth to MCFlux association and associate the flux to the slice
  art::PtrVector<simb::MCTruth> mctp;
  mctp.push_back(trueInt);
  art::FindOneP<simb::MCFlux> fofx(mctp, e, generatorLabel);
  if(fofx.isValid()) util::CreateAssn(prod,
                                      e,
                                      sliceOut,
                                      fofx.at(0),
                                      slcFlxAssn,
                                      slcIdx,
                                      fInstance);

  art::FindOneP<fxwgt::FluxWeights> fofw(mctp, e, fluxWeightLabel);
  
  if(fofw.isValid()) util::CreateAssn(prod,
                                      e,
                                      sliceOut,
                                      fofw.at(0),
                                      slcFlxWgtAssn,
                                      slcIdx,
                                      fInstance);

  LOG_VERBATIM("SkimmingUtils")
  << "The FindOne for the FluxWeights is valid: "
  << fofw.isValid()
  << " and there are "
  << slcFlxWgtAssn.size()
  << " associations "
  << fofw.at(0)->GetHadronProductionCentralValue();
  
  art::FindOneP<simb::GTruth> fogt(mctp, e, generatorLabel);
  if(fogt.isValid()) util::CreateAssn(prod,
                                      e,
                                      sliceOut,
                                      fogt.at(0),
                                      slcGTruAssn,
                                      slcIdx,
                                      fInstance);


  return true;
}


#endif // SkimmingUtils