NCNNKeras_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       NCNNKeras
// Module Type: producer
// \file        NCNNKeras_module.cc
// \brief       A module to write out the NC Cosmic Rejection NN (Keras)
//
// \author      Shaokai Yang - yangs9uc@gmail.com
////////////////////////////////////////////////////////////////////////
// system includes
#include <memory>
#include <string>
#include <vector>

// art includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindManyP.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 "CVN/func/Result.h"
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "NCID/func/KerasModel.h"
#include "NCID/NCKeras.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "nusimdata/SimulationBase/MCTruth.h"
//#include "ShowerLID/ParticleIDAlg.h"
#include "ShowerLID/ShowerLID.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/AssociationUtil.h"

// ROOT includes
#include "TFile.h"
#include "TH1D.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TSystem.h"
#include "TTree.h"
#include "TVector3.h"
#include "TVectorD.h"

namespace ncid
{
  class NCNNKeras : public art::EDProducer
  {
  public:
    
    explicit NCNNKeras(fhicl::ParameterSet const & p);
    ~NCNNKeras();
    
    void beginJob() override;
    void beginRun(art::Run &evt) override;
    void reconfigure(const fhicl::ParameterSet& pset);
    void produce(art::Event & e) override;
    
  private:

    std::string fSliceLabel;
    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;
    
    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;

    /// Particle ID alorithm (loglikelihoods and dE/dx)
    //slid::ParticleIDAlg* fParticleIDAlg;
    /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
    fhicl::ParameterSet  fParticleIDAlgPSet;
    /// Return transverse momentum fraction of the event. Calculation based on reconstructed prongs.
    double TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const;
  };
  
  bool CompareByCalEnergy(const art::Ptr<rb::Prong> a, 
                          const art::Ptr<rb::Prong> b);
}
//////////////////////////////////////////////////////////////////////////
namespace ncid
{
  NCNNKeras::NCNNKeras(fhicl::ParameterSet const & p) :
      EDProducer(p),
      fMCTruthModuleLabel(p.get<std::string>("MCTruthModuleLabel")),
    //fParticleIDAlg(0),
    fParticleIDAlgPSet(p.get<fhicl::ParameterSet>("ParticleIDAlgPSet"))
  {
    produces< std::vector<ncid::NCKeras> > ();
    produces< art::Assns<ncid::NCKeras, rb::Cluster> > ();
    this->reconfigure(p);
  }

  //////////////////////////////////////////////////////////////////////////
  NCNNKeras::~NCNNKeras()
  {
  }

  //////////////////////////////////////////////////////////////////////////
  void NCNNKeras::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");    
  }

  //////////////////////////////////////////////////////////////////////////
  void NCNNKeras::beginRun(art::Run& 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 NCNNKeras::reconfigure(fhicl::ParameterSet const & p)
  {
    fMCTruthModuleLabel   = p.get< std::string >("MCTruthModuleLabel");
    fSliceLabel           = p.get< std::string >("SliceLabel");
    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");
    fParticleIDAlgPSet    = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");  
  }
  
  //////////////////////////////////////////////////////////////////////////
  void NCNNKeras::produce(art::Event & evt)
  {

    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;
    }
    
    art::ServiceHandle<geo::Geometry>      geom;
    art::ServiceHandle<cheat::BackTracker> bt;
        
    std::unique_ptr< std::vector<ncid::NCKeras> >
      ncCosKeras(new std::vector<ncid::NCKeras>);
    std::unique_ptr< art::Assns<ncid::NCKeras, rb::Cluster> >
      sliceAssn(new art::Assns<ncid::NCKeras, rb::Cluster>);
    
    // 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);
    
    // Get things associated with slices
    art::FindManyP<rb::Prong>      prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
    art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
    
    const int nslices = slices.size();     
    for(int isli = 0; isli < nslices; ++isli){
      
      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;
      const simb::MCNeutrino& nu = SliceTruth->GetNeutrino();
      art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
      
      if(thisSlice->IsNoise())  continue;
      if(!prongAssn.isValid())  continue;
      
      int intCurrent  = nu.CCNC();
      
      if(intCurrent != simb::kNC) continue; // select only CC
      
      ncid::NCKeras nckeras;
      std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
      if(prongs.empty()) continue;
      
      nckeras.SetProngTransMom(TransMomFraction(prongs));
      
      art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
      if(!fmCVN.isValid()) continue; // CVN Validation
      std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);   
      cosmicid    = -5.; // CVN values
      if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
      shwGap      = -5.; // LIDShower values
      shwnhit     = -5.;
      shwnhitx    = -5.;
      shwnhity    = -5.;
      shwwwidth   = -5.;
      shwcalE     = -5.;        
      partptp     = nckeras.ProngTransMom();
      shwdirY     = -5.;
      shwlen      = -5.;
      kerasValue  = -5.;
      
      if(!showerLidAssn.isValid()) continue; // showerLidAssn Validation
      
      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);
      if(!foShLID.isValid()) continue; // foShLID validation
      
      cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
      art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
      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    = thisSlice->NCell();

      // 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;
      ncCosKeras->push_back(nckeras);
      
      util::CreateAssn(evt, *ncCosKeras, thisSlice ,*sliceAssn);
      fTree->Fill();
      } // end loop over slices
    evt.put(std::move(sliceAssn));
    evt.put(std::move(ncCosKeras));
  } // end of producer
  
  ///////////////////////////////////////////////////////////////////////
  double NCNNKeras::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
  {
    if(prongs.empty()) return 0.;
    
    TVector3 ret;
    for(const art::Ptr<rb::Prong>& prong: prongs){
      if(!prong->Is3D()) continue;
      const double w = prong->CalorimetricEnergy();
      ret += w*prong->Dir();
    }
    
    ret = ret.Unit();
    
    art::ServiceHandle<geo::Geometry> geom;
    const TVector3 beamDir = geom->NuMIBeamDirection();
    const double beamProj  = ret.Dot(beamDir);
    
    return (ret-beamProj*beamDir).Mag();
  }
  
  bool CompareByCalEnergy(const art::Ptr<rb::Prong> a,
                          const art::Ptr<rb::Prong> b)
  {
    return a->CalorimetricEnergy() > b->CalorimetricEnergy();
  }
  
} // end namespace ncid

DEFINE_ART_MODULE(ncid::NCNNKeras) 
////////////////////////////////////////////////////////////////////////