WSNumu_module.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////
/// \brief  Wrong-sign BDT estimator for Numu analysis
/// \author Bruno Zamorano B.Zamorano@sussex.ac.uk
////////////////////////////////////////////////////////
#include <string>
#include "TMath.h"

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "WrongSignNumu/WSNumuObj.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "ReMId/ReMId.h"
#include "NumuEnergy/NumuE.h"
#include "TimingFit/TimingFitAlg.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "CVN/func/Result.h"

#include "TFile.h"
#include "TMVA/Reader.h"

namespace wsnumu {

  class WSNumu : public art::EDProducer {
  public:
    explicit WSNumu(fhicl::ParameterSet const & pset);
    virtual ~WSNumu();

    void produce (art::Event & evt);
    void beginRun(art::Run& run);

  protected:

    std::string        fSlicerLabel;
    std::string        fTrackLabel;
    std::string        fPIDLabel;
    std::string        fEnergyLabel;
    std::string        fTMVApath;
    bool               fUseLongestTrack;
    TMVA::Reader*      fReader;
    float              TMVAvars[6];
  };
}


namespace wsnumu
{

  //----------------------------------------------------------------------
  WSNumu::WSNumu(fhicl::ParameterSet const & pset):
    fSlicerLabel      (pset.get< std::string >("SlicerLabel")),
    fTrackLabel       (pset.get< std::string >("TrackLabel")),
    fPIDLabel         (pset.get< std::string >("PIDLabel")),
    fEnergyLabel      (pset.get< std::string >("EnergyLabel")),
    fTMVApath         (pset.get< std::string >("tmvapath")),
    fUseLongestTrack  (pset.get< bool        >("UseLongestTrack")),
    fReader           (nullptr)
  {
    produces< std::vector<wsnumu::WSNumuObj> >();
    produces< art::Assns<wsnumu::WSNumuObj, rb::Cluster> >();

  }


  //----------------------------------------------------------------------
  WSNumu::~WSNumu()
  {
    if(fReader)  delete fReader;
  }

  //-------------------------------------------------------------------
  void WSNumu::beginRun(art::Run& run)
  {

    art::ServiceHandle<geo::Geometry>  geom;

    if(geom->DetId() == novadaq::cnv::kNEARDET) {

      fReader = new TMVA::Reader;
      fReader->AddVariable("theta",&TMVAvars[0]);
      fReader->AddVariable("hadEhit",&TMVAvars[1]);
      fReader->AddVariable("hadEfrac",&TMVAvars[2]);
      fReader->AddVariable("muonE",&TMVAvars[3]);
      fReader->AddVariable("nPlaneGap",&TMVAvars[4]);
      fReader->AddSpectator("PDGcode",&TMVAvars[5]);
      std::string pidpath = util::EnvExpansion(fTMVApath);
      fReader->BookMVA("BDTG",pidpath+"WSNumu_BDTG_prod3_NDtraining.weights.xml");

    }
    
  } //end beginRun
  
  //-------------------------------------------------------------------
  void WSNumu::produce(art::Event& evt)
  {
      
    art::ServiceHandle<geo::Geometry>  geom;
    art::Handle<std::vector<rb::Cluster> > slicevec;
    evt.getByLabel(fSlicerLabel,slicevec);
    
    if(slicevec->empty()) {
      mf::LogWarning ("No Slices")<<"No Slices in the input file";
      return;
    }
    
    art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);
    art::FindOneP<numue::NumuE> sliceEnergyList(slicevec, evt, fEnergyLabel);

    std::unique_ptr< std::vector<wsnumu::WSNumuObj> > outputObjects(new std::vector<wsnumu::WSNumuObj>);
    std::unique_ptr< art::Assns<wsnumu::WSNumuObj, rb::Cluster> > assoc(new art::Assns<wsnumu::WSNumuObj, rb::Cluster>);

    for(size_t i = 0; i < slicevec->size(); ++i){

      art::Ptr<rb::Cluster> slice(slicevec, i);
      
      if(slice->IsNoise()) {
        continue;
      }

      wsnumu::WSNumuObj wsnumuobj; // initialization in constructor
      
      unsigned int bestIdx = 999;  // for index in list of the highest ReMId PID track

      if(!trackAssnList.isValid()) {
        mf::LogWarning("WSNumu") << "WSNumu: No Kalman Tracks!";
        outputObjects->push_back(wsnumuobj);
        util::CreateAssn(*this, evt, *(outputObjects.get()), slice, *(assoc.get()));
        continue;
      }

      const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);
      art::FindOneP<remid::ReMId> remidVec(sliceTracks, evt, fPIDLabel);
      
      if(fUseLongestTrack) {
        bestIdx = 0;
        if(sliceTracks.size() == 0) {
          outputObjects->push_back(wsnumuobj);
          continue;
        } 
        const art::Ptr<rb::Track> track = sliceTracks[0];
        if(!track->Is3D()) {
          outputObjects->push_back(wsnumuobj);
          continue;
        }
      }
      else if(!remidVec.isValid()) {
        mf::LogWarning("WSNumu") << "WSNumu: No ReMId information!";
        outputObjects->push_back(wsnumuobj);
        util::CreateAssn(*this, evt, *(outputObjects.get()), slice, *(assoc.get()));
        continue;
      }
      else {
        bestIdx = remid::HighestPIDTrack(sliceTracks, fPIDLabel, evt);
      }
      

      if(geom->DetId() == novadaq::cnv::kNEARDET) {  // Only for ND
        // Main condition
        if(bestIdx!=999) {

          float theta     = 0.;
          float hadEhit   = 0.;
          float hadEfrac  = 0.;
          float muonE     = 0.;
          int nPlaneGap   = 0;
        
          const art::Ptr<rb::Track> track = sliceTracks[bestIdx];
          const art::Ptr<numue::NumuE> energy = sliceEnergyList.at(i);
          if (!energy) continue;
        
          double tdx = track->Dir().X();
          double tdy = track->Dir().Y();
          double tdz = track->Dir().Z();

          if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)) continue; // mis-reconstructed track

          theta = TMath::ACos(tdz)*180./TMath::Pi();
          wsnumuobj.SetTheta(theta);

          if( energy->HadCluster().NCell() > 0 ) {
            hadEhit = energy->HadCalE() / energy->HadCluster().NCell();
          }
          wsnumuobj.SetHadEhit(hadEhit);

          if( energy->TrkCCE() > 0 ) {
            hadEfrac = energy->RecoTrkCCHadE() / energy->TrkCCE();
          }
          wsnumuobj.SetHadEfrac(hadEfrac);

          if( energy->RecoMuonE() > 0 ) {
            muonE = energy->RecoMuonE();
          }
          wsnumuobj.SetMuonE(muonE);

          nPlaneGap = track->NMissingPlanes(geo::kXorY);
          wsnumuobj.SetNPlaneGap(nPlaneGap);    
                  
          //---
          TMVAvars[0] = wsnumuobj.Theta();
          TMVAvars[1] = wsnumuobj.HadEhit();
          TMVAvars[2] = wsnumuobj.HadEfrac();
          TMVAvars[3] = wsnumuobj.MuonE();
          TMVAvars[4] = wsnumuobj.NPlaneGap();
          TMVAvars[5] = 0; // Spectator(s)

          wsnumuobj.SetWSPID(fReader->EvaluateMVA("BDTG"));
        }       
      } // end if good track

      outputObjects->push_back(wsnumuobj);  // write out object with all information
      util::CreateAssn(*this, evt, *(outputObjects.get()), slice, *(assoc.get()));

    } //end loop over slices

    evt.put(std::move(outputObjects));
    evt.put(std::move(assoc));

  } //end produce

} //end namespace

DEFINE_ART_MODULE(wsnumu::WSNumu)