CVNEventProngDump_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
// \file    CVNClassify_module.cc
// \brief   Analyzer for creating both slice and prong by prong CVN
//          PixelMap objects 
// \author  Fernanda Psihas - psihas@fnal.gov
////////////////////////////////////////////////////////////////////////

// C/C++ includes
#include <iostream>
#include <sstream>

// ROOT includes
#include "TTree.h"
#include "TH2F.h"
#include "TVector3.h"

// Framework includes
#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 "Geometry/Geometry.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Vertex.h"
#include "MCCheater/BackTracker.h"
#include "ShowerLID/ShowerLID.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "Utilities/AssociationUtil.h"

#include "CVN/func/AssignLabels.h"
#include "CVN/func/TrainingData.h"
#include "CVN/func/InteractionType.h"
#include "CVN/func/FinalStateType.h"
#include "CVN/func/FinalStateProngsType.h"
#include "CVN/func/EventLabeledPMaps.h"
#include "CVN/func/ProngType.h"
#include "CVN/func/ProngTrainingData.h"

#include "RecoValidation/CAFCutter.h"
#include "StandardRecord/StandardRecord.h"


namespace cvn {
  class CVNEventProngDump : public art::EDAnalyzer {
  public:
    explicit CVNEventProngDump(fhicl::ParameterSet const& pset);    
    ~CVNEventProngDump();                        
    
    void analyze(const art::Event& evt); 
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void endJob();


    
  private:

    std::string fClusterLabel;

    /// Module label for input prongs
    std::string fProngModLabel;

    /// Instance label for input 3D prongs
    std::string fProng3DLabel;

    /// Instance label for input 2D prongs
    std::string fProng2DLabel;

    /// Instance label for cluster pixelmaps
    std::string fClusterPMLabel;

    /// Instance label for prong pixelmaps
    std::string fProngPMLabel;

    /// Instance label for vertex
    std::string fVertexLabel;

    /// Instance label for showers
    std::string fShowerLabel;

    /// Instance label for showers
    std::string fShowerLIDLabel;

    /// Option to save maps for prongs
    bool        fMapProngs;

    std::string fPixelMapInput;
    bool        fWriteMapTH2;

    bool fObeyPreselection;
    std::vector<std::string> fPreselectionLabels;

    EventLabeledPMaps* fEventPMapsData;

    TTree*        fTrainTree;

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

    /// Function to extract TH2 from PixelMap and write to TFile
    void WriteMapTH2(const art::Event& evt, const PixelMap& pm,
                     bool isProng, unsigned int slice, unsigned int id);

    std::string fCAFLabel;     ///< label for process that made standard records
    bool        fApplyCAFCuts; ///< do you want to apply CAF-level cuts?
    int         fCutType;      ///< what cuts do you want to apply? (see CAFCutter.h)

    recovalid::CAFCutter fCAFCutter; ///< CAFCutter helper class
  };



  //.......................................................................
  CVNEventProngDump::CVNEventProngDump(fhicl::ParameterSet const& pset) 
    : EDAnalyzer(pset)
  {
    this->reconfigure(pset);
  }

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

  //......................................................................
  void CVNEventProngDump::reconfigure(const fhicl::ParameterSet& pset)
  {
    fClusterLabel   = pset.get<std::string>("ClusterLabel");
    fPixelMapInput  = pset.get<std::string>("PixelMapInput");
    fProngModLabel  = pset.get<std::string>("ProngModLabel");
    fProng3DLabel   = pset.get<std::string>("Prong3DLabel");
    fProng2DLabel   = pset.get<std::string>("Prong2DLabel");
    fClusterPMLabel = pset.get<std::string>("ClusterPMLabel");
    fProngPMLabel   = pset.get<std::string>("ProngPMLabel");
    fVertexLabel    = pset.get<std::string>("VertexLabel");
    fShowerLabel    = pset.get<std::string>("ShowerLabel");
    fShowerLIDLabel = pset.get<std::string>("ShowerLIDLabel");
    fMapProngs      = pset.get<bool>       ("MapProngs");
    fWriteMapTH2    = pset.get<bool>       ("WriteMapTH2");
    fObeyPreselection = pset.get<bool>     ("ObeyPreselection");
    fPreselectionLabels = pset.get< std::vector<std::string> >("PreselectionLabels");
    fCAFLabel       = pset.get<std::string>("CAFLabel");
    fApplyCAFCuts   = pset.get<bool>       ("ApplyCAFCuts");
    fCutType        = pset.get<int>        ("CutType");

  }
      
  //......................................................................
  void CVNEventProngDump::beginJob()
  {

    art::ServiceHandle<art::TFileService> tfs;

    fEventPMapsData = new EventLabeledPMaps;
    fTrainTree = tfs->make<TTree>("CVNTrainTree", "Training records");
    fTrainTree->Branch("EventPMapsData", "cvn::EventLabeledPMaps", &fEventPMapsData);
    
  }

  //......................................................................
  void CVNEventProngDump::endJob()
  {
  }

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

    art::ServiceHandle<geo::Geometry> geom;
    
    InteractionType           interaction = kOther;
    FinalStateType             finalstate = kOtherFS;
    FinalStateProngsType finalstateprongs = kOtherFSp;

    // get the slices
    art::Handle< std::vector< rb::Cluster > > clusterHandle;
    evt.getByLabel(fClusterLabel, clusterHandle);
    const std::vector<rb::Cluster> & clusters = *clusterHandle;
    art::PtrVector<rb::Cluster> slices;
    for(unsigned int i=0; i<clusterHandle->size(); ++i){
      slices.push_back(art::Ptr<rb::Cluster>(clusterHandle, i));
    }

    // Get the Standard Records
    art::FindManyP<caf::StandardRecord> fmpSR(slices, evt, fCAFLabel);

    art::FindManyP<PixelMap> fmPixelMap(clusterHandle, evt, fPixelMapInput);

    std::vector<int> bestNuId;
    std::vector<cheat::NeutrinoWithIndex> nus;
    std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur;
    nus = fBT->allMCTruth();
    sEffPur = fBT->SlicesToMCTruthsTable(slices);
    bestNuId =fBT->SliceToOrderedNuIdsByEnergy(nus, sEffPur);


    for(size_t iClust = 0; iClust < clusters.size(); ++iClust){

      if(fObeyPreselection && rb::IsFiltered(evt, clusterHandle, iClust,
                                             fPreselectionLabels)) continue;
      // get record associated with the slice
      std::vector< art::Ptr<caf::StandardRecord> > srs;
      bool pass = true;
      // Apply the CAF-level cuts
      if(fApplyCAFCuts && fmpSR.isValid()){
        srs = fmpSR.at(iClust);
        if(srs.size() > 0) pass = fCAFCutter.passCuts(fCutType, srs[0].get());
      }

      // if we don't pass the CAF cuts (when applying), skip this slice
      if(!pass) continue;
      
      const rb::Cluster& cluster = clusters[iClust];
      if(!fmPixelMap.isValid()) continue;

      const std::vector<art::Ptr<PixelMap>> pixelMaps = fmPixelMap.at(iClust);
      if(pixelMaps.empty()) continue;

      art::Ptr<simb::MCTruth> truth;
      if(bestNuId[iClust] != -1) truth = sEffPur[iClust][bestNuId[iClust]].neutrinoInt;
      else continue;

      TrainingData train;
      float nuEnergy, lepEnergy;
      double vtxx, vtxy, vtxz;
      interaction       = SliceClassify(truth, cluster, nuEnergy, lepEnergy, vtxx, vtxy, vtxz);
      finalstate        = GetFinalStateType(truth);
      finalstateprongs  = GetFinalStateProngsType(truth);
       
      train = TrainingData(interaction, finalstate, finalstateprongs, nuEnergy, lepEnergy, *pixelMaps[0]);
  
      if(fWriteMapTH2) WriteMapTH2(evt, train.fPMap, 0, (unsigned int)iClust, 0);


      if(fMapProngs){

        art::FindManyP<rb::Vertex> fmv(clusterHandle, evt, fVertexLabel);
        art::FindManyP<rb::Prong> fmProng3D(clusterHandle, evt, art::InputTag(fProngModLabel, fProng3DLabel));
        art::FindManyP<rb::Prong> fmProng2D(clusterHandle, evt, art::InputTag(fProngModLabel, fProng2DLabel));

        std::vector< art::Ptr<rb::Prong>  > prongs3D;
        std::vector< art::Ptr<rb::Prong>  > prongs2D;
        std::vector< art::Ptr<rb::Vertex> > vert;

        if(fmProng2D.isValid() ) prongs2D = fmProng2D.at(iClust);
        if(fmProng3D.isValid() ) prongs3D = fmProng3D.at(iClust);
        if(fmv.isValid()       ) vert     = fmv.at(iClust);
        if(!fmProng2D.isValid() && !fmProng3D.isValid()) continue;

        art::FindManyP<PixelMap> fmPixelMap3D(prongs3D, evt, fPixelMapInput);
        art::FindManyP<PixelMap> fmPixelMap2D(prongs2D, evt, fPixelMapInput);
        art::FindOneP<rb::Shower>        foSh(prongs3D, evt, fShowerLabel);

        TVector3 v;
        if(!vert.empty() ) v = vert[0]->GetXYZ();
        else v.SetXYZ(-999.,-999.,-999.);


        double   length;
        TVector3 prongdir;
        TVector3 beamdir  = geom->NuMIBeamDirection();
        float    angle;


        for(unsigned int iProng = 0; iProng < prongs3D.size(); ++iProng){

          length = -5.0;
          angle  = -5.0;

          if(!fmPixelMap3D.isValid()) continue;

          const std::vector< art::Ptr<PixelMap> > pixelMaps3D = fmPixelMap3D.at(iProng);
          if(pixelMaps3D.empty()) continue;

          // Get gap to vertex from showerlid
          double gap = -5.0;
          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();
              gap = showerLID->Gap();
            }
          }

          // Currently not stored, for pm drawing only
          ProngTrainingData p_train;
          ProngType ptype3D, ptypeX, ptypeY;
          double purity3D, purityX, purityY, recE;
          unsigned int ncellX, ncellY;
          bool isprimary;
          
          TrainingData pmtrain(interaction, finalstate, finalstateprongs,
                               nuEnergy, lepEnergy, *pixelMaps3D[0]);

          length = prongs3D[iProng]->TotalLength();
          angle  = prongdir*beamdir;

          ptype3D = ProngClassify(*prongs3D[iProng], &ptype3D, &ptypeX, &ptypeY,
                                  &isprimary, &purity3D, &purityX, &purityY,
                                  &recE, &ncellX, &ncellY);       

          p_train = ProngTrainingData(ptype3D, ptypeX, ptypeY, isprimary,
                                      purity3D, purityX, purityY, recE,
                                      ncellX, ncellY, v.X(), v.Y(), v.Z(),
                                      length, (double)angle, gap, *pixelMaps3D[0]);



          if(fWriteMapTH2) WriteMapTH2(evt, p_train.fProngPMap, 1, (unsigned int)iClust, iProng+1);

          fEventPMapsData->fSliceLabeledPM = train;
          fEventPMapsData->fProngLabeledPMs = p_train;
          fTrainTree->Fill();
          
        } // loop over 3D prongs
        
        for(unsigned int iProng2 = 0; iProng2 < prongs2D.size(); ++iProng2){
          
          if(!fmPixelMap2D.isValid()) continue;

          const std::vector< art::Ptr<PixelMap> > pixelMaps2D = fmPixelMap2D.at(iProng2);
          if(pixelMaps2D.empty()) continue;

          ProngTrainingData p2_train;
          ProngType ptype2D, ptypeX2D, ptypeY2D;
          double purity2D, purityX2D, purityY2D, recE;
          unsigned int ncellX, ncellY;
          bool isprimary;

          ptype2D  = ProngClassify(*prongs2D[iProng2], &ptype2D, &ptypeX2D, &ptypeY2D,
                                   &isprimary, &purity2D, &purityX2D, &purityY2D,
                                   &recE, &ncellX, &ncellY);
          
          p2_train = ProngTrainingData(ptype2D, ptypeX2D, ptypeY2D, isprimary,
                                       purity2D, purityX2D, purityY2D, recE,
                                       ncellX, ncellY, v.X(), v.Y(), v.Z(),
                                       0.0, 0.0, 0.0, *pixelMaps2D[0]);

          fEventPMapsData->fSliceLabeledPM = train;
          fEventPMapsData->fProngLabeledPMs = p2_train;
          fTrainTree->Fill();
        
        } // loop over 2D prongs
      } // if(fMapProngs)


    } // loop over slices
  } // analyze

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

  void CVNEventProngDump::WriteMapTH2(const art::Event& evt, const PixelMap& pm,
                                      bool isProng, unsigned int slice, unsigned int id)
  {
    
    TString mapid = "_r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"
      +std::to_string(evt.event())+"_sl"+std::to_string(slice)+"_p"+std::to_string(id);
    TString data  = "ProngPixelMap";
    if (!isProng) data  = "SlicePixelMap";
    
    TH2F* hist  = pm.ToTH2(data+mapid);
    TH2F* histX = pm.SingleViewToTH2(0, data+"ViewX"+mapid);
    TH2F* histY = pm.SingleViewToTH2(1, data+"ViewY"+mapid);
    
    art::ServiceHandle<art::TFileService> tfs;
    
    TH2F* histWrite = tfs->make<TH2F>(*hist);
    histWrite->Write();
    TH2F* histWriteX = tfs->make<TH2F>(*histX);
    histWriteX->Write();
    TH2F* histWriteY = tfs->make<TH2F>(*histY);
    histWriteY->Write();
    
    delete hist;
    delete histWrite;
    delete histX;
    delete histWriteX;
    delete histY;
    delete histWriteY;
    
  }
  
  DEFINE_ART_MODULE(CVNEventProngDump)
} // end namespace cvn
////////////////////////////////////////////////////////////////////////