NCNNKerasVal_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       NCCosRej
// Module Type: analyzer
// \file        NCNNKerasVal_module.cc
// \brief       A module to for Validation 
// \author      Shaokai Yang - yangs9uc@gmail.com
////////////////////////////////////////////////////////////////////////
// system includes
#include <memory>
#include <string>
#include <vector>

// art includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "canvas/Persistency/Common/FindMany.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_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Utilities/InputTag.h"
#include "cetlib_except/exception.h"
#include "cetlib/filesystem.h"
#include "cetlib/search_path.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// novasoft includes
#include "CosRej/NueCosRej.h"
#include "CVN/func/Result.h"
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "NCID/func/KerasModel.h"
#include "NCID/NCKeras.h"
//#include "ShowerLID/ParticleIDAlg.h"
#include "ShowerLID/ShowerLID.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/AssociationUtil.h"

// ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TMath.h"
#include "TSystem.h"
#include "TTree.h"
#include "TVector3.h"


namespace ncid
{
  class NCNNKerasVal : public art::EDAnalyzer
  {
  public:
    explicit NCNNKerasVal(fhicl::ParameterSet const& p);
    virtual ~NCNNKerasVal();
    void analyze(const art::Event& evt);   
    void beginJob();   
    void beginRun(art::Run const& run);
    void reconfigure(fhicl::ParameterSet const& p);
 
  private:
    std::string fSliceLabel;
    std::string fCellHitLabel;  
    std::string fVertexLabel;
    std::string fProngLabel;
    std::string fProngInstance;
    std::string fCVNLabel;
    std::string fShowerLIDLabel;
    std::string fShowerLabel;
    std::string fCosmicTrackLabel;
    std::string fMCTruthModuleLabel;
    std::string fPIDLibPath;
    std::string fNCKerasFile;
    std::string fNCCosRejPIDFile;
    std::string fCosRejLabel;
    keras::KerasModel* keras_5nn;
    TTree* fTree;

    float cosmicid;
    float shwGap;
    float shwnhit;
    float shwnhitx;
    float shwnhity;
    float shwwwidth;
    float shwcalE;
    float partptp;
    float shwdirY;
    float shwlen;
    float shwxminusy;
    float shwxplusy;
    float shwxovery;
    float nmiphit;
    float nshwlid;
    float kerasValue;
    float abdtValue;
    float gbdtValue;
    //slid::ParticleIDAlg* fParticleIDAlg;
    fhicl::ParameterSet  fParticleIDAlgPSet;
  };
  
  bool CompareByCalEnergy(const art::Ptr<rb::Prong> a, 
                          const art::Ptr<rb::Prong> b);
}

//////////////////////////////////////////////////////////////////////////
namespace ncid
{
  //////////////////////////////////////////////////////////////////////////
  NCNNKerasVal::NCNNKerasVal(fhicl::ParameterSet const& p) :
    EDAnalyzer(p),
    fMCTruthModuleLabel(p.get< std::string>("MCTruthModuleLabel")),
    //fParticleIDAlg(0),
    fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
  { 
    this->reconfigure(p);
  }

  //////////////////////////////////////////////////////////////////////////
  NCNNKerasVal::~NCNNKerasVal()
  {
  }
  
  //////////////////////////////////////////////////////////////////////////
  void NCNNKerasVal::reconfigure(fhicl::ParameterSet const& p)
  {
    fMCTruthModuleLabel   = p.get< std::string >("MCTruthModuleLabel");
    fSliceLabel           = p.get< std::string >("SliceLabel");
    fCellHitLabel         = p.get< std::string >("CellHitLabel");
    fProngLabel           = p.get< std::string >("ProngLabel");
    fProngInstance        = p.get< std::string >("ProngInstance");
    fCVNLabel             = p.get< std::string >("CVNLabel");
    fShowerLabel          = p.get< std::string >("ShowerLabel");
    fShowerLIDLabel       = p.get< std::string >("ShowerLIDLabel");
    fPIDLibPath           = p.get< std::string >("PIDLibPath");
    fNCKerasFile          = p.get< std::string >("NCKerasFile");
    fCosRejLabel          = p.get< std::string >("CosRejLabel");
    fParticleIDAlgPSet    = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");  
  }  

  //////////////////////////////////////////////////////////////////////////
  void NCNNKerasVal::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fTree = tfs->make<TTree>("Testtree", "Testtree");
    
    fTree->Branch("cosmicid", &cosmicid, "cosmicid/F");    
    fTree->Branch("shwGap", &shwGap, "shwGap/F");    
    fTree->Branch("shwnhit", &shwnhit, "shwnhit/F");    
    fTree->Branch("shwnhitx", &shwnhitx, "shwnhitx/F");    
    fTree->Branch("shwnhity", &shwnhity, "shwnhity/F");    
    fTree->Branch("shwwwidth", &shwwwidth, "shwwwidth/F");    
    fTree->Branch("shwcalE", &shwcalE, "shwcalE/F");    
    fTree->Branch("partptp", &partptp, "partptp/F");    
    fTree->Branch("shwdirY", &shwdirY, "shwdirY/F");    
    fTree->Branch("shwlen", &shwlen, "shwlen/F");    
    fTree->Branch("shwxminusy", &shwxminusy, "shwxminusy/F");    
    fTree->Branch("shwxplusy", &shwxplusy, "shwxplusy/F");    
    fTree->Branch("shwxovery", &shwxovery, "shwxovery/F");    
    fTree->Branch("nmiphit", &nmiphit, "nmiphit/F");    
    fTree->Branch("nshwlid", &nshwlid, "nshwlid/F");    
    fTree->Branch("kerasValue", &kerasValue, "kerasValue/F");    
    fTree->Branch("abdtValue", &abdtValue, "abdtValue/F");    
    fTree->Branch("gbdtValue", &gbdtValue, "gbdtValue/F");    
  }

  //////////////////////////////////////////////////////////////////////////
  void NCNNKerasVal::beginRun(art::Run const& run)
  {
    std::string pidlibpath = util::EnvExpansion(fPIDLibPath);
    // Make sure we can find the Neural Networks weight file       
    if(!cet::file_exists(pidlibpath+fNCKerasFile)){
      throw cet::exception("NCKeras")
        << "Couldn't find file " << pidlibpath+fNCKerasFile <<std::endl;
    }
    keras_5nn = new keras::KerasModel(pidlibpath+fNCKerasFile, false);     
  }

  //////////////////////////////////////////////////////////////////////////
  void NCNNKerasVal::analyze(const art::Event& evt)
  {
    art::ServiceHandle<cheat::BackTracker> bt;
    bool isMCTruthListEmpty = false;
    art::Handle< std::vector<simb::MCTruth> > mclist;
    try {
      evt.getByLabel(fMCTruthModuleLabel, mclist);
      if(mclist->empty()) isMCTruthListEmpty = true;
    }
    catch(...) { isMCTruthListEmpty = true; }
    if(isMCTruthListEmpty){
      mf::LogError("MCTruth Error") << "Error retrieving MCTruth list." << std::endl;
      return;
    }   
   
    // Get the slices in the event
    art::Handle<std::vector<rb::Cluster> > sliceHandle;
    evt.getByLabel( fSliceLabel, sliceHandle);
    
    std::vector<art::Ptr<rb::Cluster> > slices;
    art::fill_ptr_vector(slices, sliceHandle);
    
    art::Handle<std::vector<rb::CellHit>> htHandle;
    evt.getByLabel(fCellHitLabel, htHandle);
    
    art::PtrVector<rb::CellHit> allHits;
    for(unsigned int m = 0; m < htHandle->size(); ++m){
      allHits.push_back(art::Ptr<rb::CellHit>(htHandle,m));
    }
    // Get things associated with slices
    art::FindManyP<rb::Prong>         prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
    art::FindMany<slid::ShowerLID>    showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
    art::FindManyP<cosrej::NueCosRej> fmcr(slices, evt, fCosRejLabel);    
    art::FindManyP<rb::Vertex>        vertexAssn(sliceHandle, evt, fVertexLabel);
    const int nslices = slices.size();     
    cosmicid    = -5.; // CVN values
    shwGap      = -5.; // LIDShower values 
    shwnhit     = -5.; 
    shwnhitx    = -5.;
    shwnhity    = -5.;
    shwwwidth   = -5.;
    shwcalE     = -5.;          
    partptp     = -5.; //nckeras.ProngTransMom();
    shwdirY     = -5.;
    shwlen      = -5.;
    kerasValue  = -5.;

    // Loop over slices
    for(int isli = 0; isli < nslices; ++isli){
      art::Ptr<rb::Cluster> thisSlice = slices.at(isli); 
      if(thisSlice->IsNoise())  continue;       
      if(!prongAssn.isValid())  continue;
      //Get True neutrino info
      std::vector<cheat::NeutrinoEffPur> SlicePurr = bt->SliceToMCTruth(slices[isli]->AllCells(), slices[isli]->AllCells());
      
      if(SlicePurr.size() == 0) continue;
      
      const art::Ptr<simb::MCTruth> &SliceTruth = SlicePurr[0].neutrinoInt;
      if(!SliceTruth->NeutrinoSet()) continue;
      
      ncid::NCKeras nckeras;
      std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
      if(prongs.empty()) continue;
      
      art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
      // CVN Validation
      if(!fmCVN.isValid()) continue;
      
      std::vector< art::Ptr<cosrej::NueCosRej> > cosrejs = fmcr.at(isli);
      if(cosrejs.empty()) continue;
      
      std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
      //  showerLidAssn Validation
      if(!showerLidAssn.isValid()) continue;
      
      std::vector<const slid::ShowerLID*> shwlids = showerLidAssn.at(isli);
      if(shwlids.empty()) continue;
      
      sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
      art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
      nshwlid = shwlids.size();   
      cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(0));
      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);
      // foShLID validation
      if(!foShLID.isValid()) continue;
      
      cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
      art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
      /////////////OK, Cuts///////////////
      /// Event Quality Cuts
      if(cosrejs[0]->HitsPerPlane() >= 8) continue;
      if(showers.size() == 0) continue;
      if(showerLID->Gap() >= 100.) continue;
      if(thisSlice->MostContiguousPlanes(geo::kXorY) <= 2) continue;
      /// NC Official Containment Cuts
      if(cosrejs[0]->StartDistToFront() < 10) continue;
      if(cosrejs[0]->StopDistToFront()  < 10) continue;
      if(cosrejs[0]->StartDistToBack()  < 10) continue;
      if(cosrejs[0]->StopDistToBack()   < 10) continue;
      if(cosrejs[0]->StartDistToTop()   < 10) continue;
      if(cosrejs[0]->StopDistToTop()    < 10) continue;
      if(cosrejs[0]->StartDistToWest()  < 10) continue;
      if(cosrejs[0]->StopDistToWest()   < 10) continue;
      if(cosrejs[0]->StartDistToEast()  < 10) continue;
      if(cosrejs[0]->StopDistToEast()   < 10) continue;

      int nMipHitt = 0;
      for(size_t hitIdx = 0; hitIdx < slices[isli]->NCell(); ++hitIdx){
        const art::Ptr<rb::CellHit>&  chptr = slices[isli]->Cell(hitIdx);
        const rb::RecoHit Reco = slices[isli]->RecoHit(chptr);
        if(!Reco.IsCalibrated()) continue;
        if((Reco.PECorr() > 100.0) && (Reco.PECorr() < 245.0)) nMipHitt++;
      } // end of hits loop
      
      if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
      shwGap     = showerLID->Gap();
      shwnhit    = shower->NCell();
      shwwwidth  = showerLID->Radius();
      shwcalE    = showerLID->ShowerEnergy();
      shwnhitx   = shower->NXCell();
      shwdirY    = (float)shower->Dir().Y();
      shwnhity   = shower->NYCell();
      shwlen     = shower->TotalLength();         
      shwxminusy = shwnhitx - shwnhity; 
      shwxplusy  = shwnhitx + shwnhity; 
      shwxovery  = shwxminusy/shwxplusy;              
      nmiphit    = nMipHitt;
      partptp    = cosrejs[0]->ParticleShowerTransMom();
      
      // Keras model
      keras::DataChunk *sample = new keras::DataChunkFlat();          
      float myfloats[] = {partptp,
                          cosmicid,
                          nmiphit,
                          shwcalE,
                          shwnhit,
                          nshwlid,
                          shwnhitx,
                          shwnhity,
                          shwxminusy,
                          shwxplusy,
                          shwxovery,
                          shwdirY,
                          shwlen,
                          shwwwidth,
                          shwGap};
      
      std::vector<float> v (myfloats, myfloats + sizeof(myfloats) / sizeof(float) );
      sample->set_data(v);   
      kerasValue = keras_5nn->compute_output(sample)[0];
      nckeras.SetCosPIDKeras(kerasValue);
      delete sample;
      
      fTree->Fill();
    } // end loop over slices
  } // end of analyzer  
} // end namespace ncid

DEFINE_ART_MODULE(ncid::NCNNKerasVal)