RecoValidationTutorial_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file   
// \brief  Simple demonstration of grabbing validation from art::Event
// \author Fernanda Psihas - psihas@fnal.gov
//
////////////////////////////////////////////////////////////////////////

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

// NOvA includes
#include "MCCheater/BackTracker.h"
#include "RecoBaseHit/CellHit.h"
#include "RecoBaseHit/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "TH1F.h"
#include "TVector3.h"

namespace tut
{
  class RecoValidationTutorial : public art::EDAnalyzer
  {
    
  public:
    explicit RecoValidationTutorial(const fhicl::ParameterSet& p);
    ~RecoValidationTutorial();
    
    // Required functions.
    void analyze(const art::Event& evt);

    // Selected optional functions.
    void reconfigure(const fhicl::ParameterSet& p);

    void beginJob();
    
  protected:
    
    std::string fCellHitLabel;   ///< module label for CellHit
    std::string fClusterLabel;   ///< module label for clusters
    std::string fVertexLabel;   ///< module label for vertices

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

    TH1* fVertexX;
    TH1* fVertexY;
    TH1* fVertexZ;

  };
}
  
namespace tut {

  //.......................................................................
  RecoValidationTutorial::RecoValidationTutorial(const fhicl::ParameterSet& pset)
    : EDAnalyzer(pset)
  {
    reconfigure(pset);
  }

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

  //.......................................................................
  void RecoValidationTutorial::reconfigure(const fhicl::ParameterSet& pset)
  {
    fClusterLabel = pset.get<std::string>("ClusterLabel");
    fVertexLabel  = pset.get<std::string>("VertexLabel");
  }

  //.......................................................................
  void RecoValidationTutorial::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fVertexX = tfs->make<TH1F>("vertexXdiff", ";Reco - True Vertex X (cm)",80, -200, 200);
    fVertexY = tfs->make<TH1F>("vertexYdiff", ";Reco - True Vertex Y (cm)",80, -200, 200);
    fVertexZ = tfs->make<TH1F>("vertexZdiff", ";Reco - True Vertex Z (cm)",80, -200, 200);
  }

  //-----------------------------------------------------------------------
  void RecoValidationTutorial::analyze(const art::Event& evt)
  {

    // ---------- Get data products ----------

    // Get all Slices in the event
    art::Handle< std::vector<rb::Cluster> > sliceHandle;
    evt.getByLabel(fClusterLabel, sliceHandle);
    
    art::PtrVector<rb::Cluster> slices;
    for ( unsigned int i = 0; i < sliceHandle->size(); ++i ) {
      art::Ptr<rb::Cluster> sliceptr(sliceHandle, i);
      slices.push_back(sliceptr);
    }

    
    // Get vertices associated with the slices
    art::FindManyP<rb::Vertex> fmVert(sliceHandle, evt, fVertexLabel);
    

    // Get all truths related to this slice
    std::vector<std::vector<cheat::NeutrinoEffPur>> slTruth = bt->SlicesToMCTruthsTable(slices);
    // Get all neutrinos from the event
    std::vector<cheat::NeutrinoWithIndex>               nus = bt->allMCTruth();
    std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur = bt->SlicesToMCTruthsTable(slices);
    // and order them by purity 
    std::vector<int> sliceNuId = bt->SliceToOrderedNuIds(nus, slTruth, 
                                                         cheat::EffMetric, 
                                                         cheat::EffMetric);


    // ---------- Loop Over Slices ----------
    for ( unsigned int iSlice = 0; iSlice < slices.size(); ++iSlice ){
      
      // Simple example. Print number of cells in each slice.
      const art::Ptr<rb::Cluster> slice = slices[iSlice];
      unsigned int nCells = slice->NCell();
      
      mf::LogInfo("RecoValidationTutorial") << " Slice # " << iSlice
                                             << " has "     << nCells
                                             << " cells."   << std::endl;      
      
      if ( slice->IsNoise() ) continue;
      
      
      // ---------- Reco Vertex Info ----------

      // Always check the assns are ok
      if( !fmVert.isValid() ) continue;
      
      const std::vector< art::Ptr<rb::Vertex> > vtx = fmVert.at(iSlice);
      if ( vtx.size() != 1) continue;
      
      TVector3 vxyz = vtx[0]->GetXYZ();
      mf::LogInfo("RecoValidationTutorial") << "  It has Reco Vertex at "
                                            << "( "<< vxyz.X()
                                            << ", "<< vxyz.Y()
                                            << ", "<< vxyz.Z()
                                            << " ) " <<std::endl;
      

      // ---------- True Vertex Info ----------

      // Check there are related nus
      if ( sliceNuId[iSlice] == -1 ) return;
      
      // Get this slice's highest eff. neutrino
      art::Ptr<simb::MCTruth> truth = slTruth[iSlice][sliceNuId[iSlice]].neutrinoInt;      
      simb::MCNeutrino thisNu = truth->GetNeutrino();
      
      //now get vertex info and make comparisons
      fVertexX->Fill( vxyz.X() - thisNu.Nu().Vx() );
      fVertexY->Fill( vxyz.Y() - thisNu.Nu().Vy() );
      fVertexZ->Fill( vxyz.Z() - thisNu.Nu().Vz() );
      
    } //iSlice    
  }// analyze
  //-----------------------------------------------------------------------

} // end namespce tut

// Magic. Do not touch this
DEFINE_ART_MODULE(tut::RecoValidationTutorial)