MakeNCCosRej_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       MakeNCCosRej
// Module Type: producer
// \file        MakeNCCosRej_module.cc
// \brief       A module to write out the NCCosRej object
//
// \author      Shaokai Yang - yangs9uc@gmail.com
////////////////////////////////////////////////////////////////////////
// system includes
#include <memory>
#include <string>
#include <vector>

// art includes
#include "art/Framework/Core/EDProducer.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 "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"

// ROOT includes
#include "TMath.h"
#include "TMVA/Reader.h"
#include "TVector3.h"

// novasoft includes
#include "CosRej/NueCosRej.h"
#include "CVN/func/Result.h"
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "NCID/NCCosRej.h"
#include "ShowerLID/ParticleIDAlg.h"
#include "ShowerLID/ShowerLID.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Shower.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/AssociationUtil.h"


namespace ncid
{
  class MakeNCCosRej : public art::EDProducer
  {
  public:
    
    explicit MakeNCCosRej(fhicl::ParameterSet const & p);
    ~MakeNCCosRej();

    void beginRun(art::Run &evt) override;
    void reconfigure(const fhicl::ParameterSet& p);
    void produce(art::Event & e) override;
 
  private:

    std::string fSliceLabel;
    std::string fProngLabel;
    std::string fProngInstance;
    std::string fCVNLabel;
    std::string fShowerLIDLabel;
    std::string fShowerLabel;
    std::string fCosmicTrackLabel;
    std::string fCellHitLabel;
    std::string fCosRejLabel;
    std::string fPIDLibPath;
    std::string fNCCosRejPIDFile;
    std::string fNCCosRejPIDGFile;

    TMVA::Reader* fReaderDT;
    TMVA::Reader* fReaderGBDT;

    // Variables need to send to TMVA::Reader
    float TMVAvars[13];

    /// 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
{
  MakeNCCosRej::MakeNCCosRej(fhicl::ParameterSet const & p) :
    EDProducer(p),
    fReaderDT(nullptr),
    fReaderGBDT(nullptr),
    //fParticleIDAlg(0),
    fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
  {
    produces< std::vector<ncid::NCCosRej> > ();
    produces< art::Assns<ncid::NCCosRej, rb::Cluster> > ();
    this->reconfigure(p);
  }
  
  MakeNCCosRej::~MakeNCCosRej()
  {
    if(fReaderDT)   delete fReaderDT;
    if(fReaderGBDT) delete fReaderGBDT;
  }
  
  void MakeNCCosRej::beginRun(art::Run& run)
  {
    std::string pidlibpath = util::EnvExpansion(fPIDLibPath);
    
    // Make sure we can find the Real Adaptive Boosted Decision Trees weight file before we set up MVA
    if(!cet::file_exists(pidlibpath+fNCCosRejPIDFile)){
      throw cet::exception("MakeNCCosRej")
        << "Couldn't find file " << pidlibpath+fNCCosRejPIDFile <<std::endl;
    }
    // Make sure we can find the Gradient Boosted Decision Trees weight file before we set up MVA
    if(!cet::file_exists(pidlibpath+fNCCosRejPIDGFile)){
      throw cet::exception("MakeNCCosRej")
        << "Couldn't find file " << pidlibpath+fNCCosRejPIDGFile <<std::endl;
    }
    
    art::ServiceHandle<geo::Geometry>  geom;
    
    fReaderDT = new TMVA::Reader;
    fReaderDT->AddVariable("cosmicid",   &TMVAvars[0]);  //cosmicid
    fReaderDT->AddVariable("partptp",    &TMVAvars[1]);  //partptp
    fReaderDT->AddVariable("shwnhit",    &TMVAvars[2]);  //shwnhit
    fReaderDT->AddVariable("shwxminusy", &TMVAvars[3]);  //shwminusy
    fReaderDT->AddVariable("shwxplusy",  &TMVAvars[4]);  //shwxplusy
    fReaderDT->AddVariable("shwxovery",  &TMVAvars[5]);  //shwxpvery
    fReaderDT->AddVariable("shwcalE",    &TMVAvars[6]);  //shwcalE
    fReaderDT->AddVariable("shwdirY",    &TMVAvars[7]);  //shwdirY
    fReaderDT->AddVariable("shwlen",     &TMVAvars[8]);  //shwlen
    fReaderDT->AddVariable("shwwwidth",  &TMVAvars[9]);  //shwwwidth
    fReaderDT->AddVariable("shwGap",     &TMVAvars[10]); //shwGap
    fReaderDT->AddVariable("nshwlid",    &TMVAvars[11]); //nshwlid
    fReaderDT->AddVariable("nmiphit",    &TMVAvars[12]); //nmiphit    
    fReaderDT->BookMVA("BDTA",pidlibpath+fNCCosRejPIDFile);

    fReaderGBDT = new TMVA::Reader;
    fReaderGBDT->AddVariable("cosmicid",   &TMVAvars[0]);  //cosmicid
    fReaderGBDT->AddVariable("partptp",    &TMVAvars[1]);  //partptp
    fReaderGBDT->AddVariable("shwnhit",    &TMVAvars[2]);  //shwnhit
    fReaderGBDT->AddVariable("shwxminusy", &TMVAvars[3]);  //shwminusy
    fReaderGBDT->AddVariable("shwxplusy",  &TMVAvars[4]);  //shwxplusy
    fReaderGBDT->AddVariable("shwxovery",  &TMVAvars[5]);  //shwxpvery
    fReaderGBDT->AddVariable("shwcalE",    &TMVAvars[6]);  //shwcalE
    fReaderGBDT->AddVariable("shwdirY",    &TMVAvars[7]);  //shwdirY
    fReaderGBDT->AddVariable("shwlen",     &TMVAvars[8]);  //shwlen
    fReaderGBDT->AddVariable("shwwwidth",  &TMVAvars[9]);  //shwwwidth
    fReaderGBDT->AddVariable("shwGap",     &TMVAvars[10]); //shwGap
    fReaderGBDT->AddVariable("nshwlid",    &TMVAvars[11]); //nshwlid
    fReaderGBDT->AddVariable("nmiphit",    &TMVAvars[12]); //nmiphit    
    fReaderGBDT->BookMVA("BDTG", pidlibpath+fNCCosRejPIDGFile);
}
  
  void MakeNCCosRej::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");
    fCosRejLabel          = p.get< std::string >("CosRejLabel");
    fNCCosRejPIDFile      = p.get< std::string >("NCCosRejPIDFile");
    fNCCosRejPIDGFile     = p.get< std::string >("NCCosRejPIDGFile");
    fParticleIDAlgPSet    = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");  
  }
  
  //////////////////////////////////////////////////////////////////////////
  void MakeNCCosRej::produce(art::Event & evt)
  {
    art::ServiceHandle<geo::Geometry>  geom;
    art::ServiceHandle<cheat::BackTracker> bt;
    
    std::unique_ptr< std::vector<ncid::NCCosRej> >
      ncCosRejjies(new std::vector< ncid::NCCosRej >);
    std::unique_ptr< art::Assns<ncid::NCCosRej, rb::Cluster> >
      sliceAssn(new art::Assns<ncid::NCCosRej, 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);
    
    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);

    const int nslices = slices.size();     
    for(int isli = 0; isli < nslices; ++isli){
      // Loop over slices
      art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
      
      if(thisSlice->IsNoise())  continue;
      if(!prongAssn.isValid())  continue;
      
      ncid::NCCosRej nccosrej;
      std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
      
      if(prongs.empty()) continue;
      
      nccosrej.SetProngTransMom(TransMomFraction(prongs));
      std::vector< art::Ptr<cosrej::NueCosRej> > cosrejs = fmcr.at(isli);
      
      if(cosrejs.empty()) continue;
      
      // Need CVN, input to Real Adaptive BDT
      art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);

      if(!fmCVN.isValid()) continue;
      
      std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);   
      // CVN values
      float cosmicid  = 0.;
      if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
      
      // LIDShower values 
      float shwGap    = -5.;
      float shwnhit   = -5.; 
      float shwnhitx  = -5.;
      float shwnhity  = -5.;
      float shwwwidth = -5.;
      float shwcalE   = -5.;    
      float shwdirY   = -5.;
      float shwlen    = -5.;
      float partptp   = -5.;
      /////////////////////// Fill the PIDs ////////////////////////////  
      // Only execute when there are showers associated with the shwlids
      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);
      float nshwlid = shwlids.size();
      
      if(!foSh.isValid()) continue;
      
      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;
      
      cet::maybe_ref<art::Ptr<slid::ShowerLID> const> roShLID(foShLID.at(0));
      art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();

      int nMipHit = 0;
      for(size_t hitIdx = 0; hitIdx < slices[isli]->NCell(); ++hitIdx){
        const art::Ptr<rb::CellHit>&  chptr = slices[isli]->Cell(hitIdx);
        const rb::RecoHit rhit = slices[isli]->RecoHit(chptr);
        if(!rhit.IsCalibrated()) continue;
        if((rhit.PECorr() > 100.0) && (rhit.PECorr() < 245.0)) nMipHit++;
      } // end of hits loop 
   
      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();
      partptp   = cosrejs[0]->ParticleShowerTransMom();
      float shwxminusy = shwnhitx - shwnhity; 
      float shwxplusy  = shwnhitx + shwnhity; 
      float shwxovery  = shwxminusy/shwxplusy;        
      float nmiphit    = nMipHit;


      TMVAvars[0]  = cosmicid;                     //ncid
      TMVAvars[1]  = partptp;                      //partptp
      //TMVAvars[1]  = nccosrej.ProngTransMom();     //partptp
      TMVAvars[2]  = shwnhit;                      //shwnhit;
      TMVAvars[3]  = shwxminusy;
      TMVAvars[4]  = shwxplusy;
      TMVAvars[5]  = shwxovery;
      TMVAvars[6]  = shwcalE;                     //shwcalE
      TMVAvars[7]  = shwdirY;                     //shwdirY
      TMVAvars[8]  = shwlen;                      //shwlen
      TMVAvars[9]  = shwwwidth;                   //shwwidth
      TMVAvars[10] = shwGap;                      //shwgap
      TMVAvars[11] = nshwlid;
      TMVAvars[12] = nmiphit;                     //nmiphit;
      
      nccosrej.SetCosPIDDT(fReaderDT->EvaluateMVA("BDTA"));
      nccosrej.SetCosPIDDTG(fReaderGBDT->EvaluateMVA("BDTG"));    
      ncCosRejjies->push_back(nccosrej);
      util::CreateAssn(evt, *ncCosRejjies, thisSlice ,*sliceAssn);
    } // end loop over slices
    
    evt.put(std::move(sliceAssn));
    evt.put(std::move(ncCosRejjies));
  } // end of producer
  
  ///////////////////////////////////////////////////////////////////////
  double MakeNCCosRej::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::MakeNCCosRej) 
////////////////////////////////////////////////////////////////////////