QeFinder_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////           
/// \brief  A module for finding numu CC QE interactions
/// \author  Nicholas Raddatz - raddatz@physics.umn.edu
///  \date    February 2013
//////////////////////////////////////////////////////////////////// 

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

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "fhiclcpp/ParameterSet.h"
#include "canvas/Persistency/Common/FindManyP.h"

//NOvA includes
#include "QEEventFinder/QePId.h"
#include "ReMId/ReMId.h"
#include "RecoBase/Track.h"
#include "NumuEnergy/NumuE.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "Geometry/Geometry.h"


namespace qeef{
  /// A QE Event Finder for muons 
  class QeFinder : public art::EDProducer{
  public:
    explicit QeFinder(fhicl::ParameterSet const& pset);
    virtual ~QeFinder();

    void beginJob();
    void produce(art::Event& evt);

  private:
    void Init();

    std::string fSliceLabel;
    std::string fTrackLabel;
    std::string fPIDLabel;
    std::string fEnergyLabel;
    std::string fKnnLibPath;
    std::string fknn1trkfile;
    std::string fknn2trkfile;
    bool trainfileloaded;
    TMVA::Reader *knn1trk, *knn2trk;
    Float_t offEK, ediff, ediffZ, dedx;

  };
}

namespace qeef{

  //.......................................................................
  QeFinder::QeFinder(fhicl::ParameterSet const& pset):
    knn1trk(nullptr),
    knn2trk(nullptr)
  {
    fSliceLabel    = (pset.get< std::string >("SliceLabel"));
    fTrackLabel    = (pset.get< std::string >("TrackLabel"));
    fPIDLabel      = (pset.get< std::string >("PIDLabel"));
    fEnergyLabel   = (pset.get< std::string >("EnergyLabel"));
    fKnnLibPath    = (pset.get< std::string >("KnnLibPath"));
    produces< std::vector<qeef::QePId> >();
    produces< art::Assns<qeef::QePId, rb::Cluster> >();
  }

  //.......................................................................
  QeFinder::~QeFinder()
  {
    if (knn1trk) delete knn1trk;
    if (knn2trk) delete knn2trk;
  }

  //.......................................................................
  void QeFinder::beginJob()
  {
    trainfileloaded = false;    
  }

  //.......................................................................
  void QeFinder::Init()
  {
    if(trainfileloaded) return;

    art::ServiceHandle<geo::Geometry> geom;
    // Expand fcl input lib path
    std::string knnlibpath = util::EnvExpansion(fKnnLibPath);
    // set up FD training file
    if(geom->DetId() == novadaq::cnv::kFARDET){
      fknn1trkfile = knnlibpath+"TMVAClassification_KNN_qe_1track.weights.xml"; 
      fknn2trkfile = knnlibpath+"TMVAClassification_KNN_qe_2track.weights.xml";
    }
    // set up ND training file
    else if(geom->DetId() == novadaq::cnv::kNEARDET){
      fknn1trkfile = knnlibpath+"TMVAClassification_KNN_qe_1track.nd.weights.xml"; 
      fknn2trkfile = knnlibpath+"TMVAClassification_KNN_qe_2track.nd.weights.xml";
    }
    else{ std::cerr << "This file isn't near or far detector, algorithm doesn't work on any other" << std::endl;}

    knn1trk = new TMVA::Reader();
    knn1trk->AddVariable("offEK",&offEK);
    knn1trk->AddVariable("ediff",&ediff);
    knn1trk->AddVariable("ediffZ",&ediffZ);    

    knn1trk->BookMVA("KNN method",fknn1trkfile);
      
    knn2trk = new TMVA::Reader();
    knn2trk->AddVariable("offEK",&offEK);
    knn2trk->AddVariable("ediff",&ediff);
    knn2trk->AddVariable("ediffZ",&ediffZ);
    knn2trk->AddVariable("dedx",&dedx);    

    knn2trk->BookMVA("KNN method",fknn2trkfile);
     
    trainfileloaded = true;
  }
  
  //.......................................................................
  void QeFinder::produce(art::Event& evt)
  {
    Init();

    // going to produce qe pid objects and associations with slices
    std::unique_ptr< std::vector<qeef::QePId> > qeefOut(new std::vector<qeef::QePId>);
    std::unique_ptr< art::Assns<qeef::QePId, rb::Cluster> > qeAss(new art::Assns<qeef::QePId, rb::Cluster>);
    
    
    // get all of slices out of the event
    art::Handle<std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSliceLabel,slices);
    
    // set up associations between slices and tracks 
    art::FindManyP<rb::Track> trkAss(slices,evt,fTrackLabel);
    // set up associations between slices and energy
    art::FindManyP<numue::NumuE> eAss(slices,evt,fEnergyLabel);

    for(size_t isl =0; isl < slices->size(); ++isl){

      // Get the slice
      art::Ptr<rb::Cluster> slice(slices,isl);
      // skip if it is noise
      if(slice->IsNoise() || !trkAss.isValid()){ continue; }

      // make a qepid object with default pid value of -5
      qeef::QePId qepid(14,-5);
      
      // get the tracks that are associated with this slice
      const std::vector<art::Ptr<rb::Track> > tracks = trkAss.at(isl);

      // if more than 3 or no tracks no qe pid will be made so skip this
      if(tracks.size() < 3 && tracks.size() > 0){
        
        // set up associations between tracks and remid
        art::FindOneP<remid::ReMId> remAss(tracks,evt,fPIDLabel);
        
        //Find index of highest pid track
        unsigned int bestId = remid::HighestPIDTrack(tracks, fPIDLabel, evt);

        int n3dtrk = 0; // keep track of the number of 3d tracks     
        std::vector<art::Ptr<rb::CellHit> > trkHits; // cellhits from every reco track
        
        // loop over tracks and find the highest pid one
        for(size_t itrk = 0; itrk < tracks.size(); ++itrk){
          // get the track
          const art::Ptr<rb::Track> track = tracks[itrk];
          
          if(track->Is3D()){ ++n3dtrk; }
          
          // get all of the cell hits associated with this track
          for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
            trkHits.push_back(track->Cell(ihit));
          }
          
        }

        // continue if highest pid track found
        if(bestId != 999){

          // check number of 3d trks
          if((n3dtrk == 1 && tracks.size() == 1) ||
             (n3dtrk >=1 && tracks.size() == 2)){

            // get all the hits not assocaited with a track in the event
            rb::Cluster nontrkHits;
            for(unsigned int slhit = 0; slhit < slice->NCell(); ++slhit){
              bool intrk = false;
              // check if this hit is in any track
              for(unsigned int trhit = 0; trhit < trkHits.size(); ++trhit){
                if(slice->Cell(slhit) == trkHits[trhit]){
                  intrk = true;
                  break;
                }
              }
              if(!intrk){ nontrkHits.Add(slice->Cell(slhit)); }
            }

            // knn variables
            float offE;
            float ediff;
            float ediffZ;
            float dedx;

            // get muon pid and track
            art::Ptr<remid::ReMId> murem = remAss.at(bestId);
            art::Ptr<rb::Track> mutrk = tracks[bestId];

            // get numu energy estimates
            if(eAss.isValid()){
              std::vector<art::Ptr<numue::NumuE> > numuEnergy = eAss.at(isl);
              if(numuEnergy.size() > 0){
                
                // get the energy variables for the qe pid
                float angleQEE = numuEnergy[0]->AngleQEE();
                float errAngleQEE = numuEnergy[0]->AngleQEError();
                float trkQEE = numuEnergy[0]->TrkQEE();

                // Don't pid for cases where angleQEE or trkqeE is nonsense
                if( angleQEE < 0.0 || angleQEE > 40 || errAngleQEE < 0.0 || errAngleQEE > 40 || trkQEE > 20){
                  qepid.SetNtrk(tracks.size());
                  qepid.SetVal(-1);
                }
                else{
                  
                  ediff = (angleQEE - trkQEE)/trkQEE;
                  ediffZ = (angleQEE - trkQEE)/errAngleQEE;
                  
                  if(tracks.size() == 1){
                    double vertE = 0.0;
                    if(nontrkHits.NCell() > 0){
                      vertE = nontrkHits.TotalGeV();
                    }
                    offE = vertE/mutrk->TotalGeV();
                    
                    // gather all pid variables 
                    std::vector<float> knnInputs(3);
                    knnInputs[0] = offE;
                    knnInputs[1] = ediff;
                    knnInputs[2] = ediffZ;
                    double kNN = knn1trk->EvaluateMVA(knnInputs,"KNN method");
                    
                    // update qe pid info
                    qepid.SetNtrk(1);
                    qepid.SetOffTrkE(offE);
                    qepid.SetEdiff(ediff);
                    qepid.SetEdiffZ(ediffZ);
                    qepid.SetVal(kNN);
                  } // end of 1 track case
                  else{
                    int nonmuId = 1;
                    if(bestId == 1){ nonmuId = 0; }
                    art::Ptr<remid::ReMId> nonmurem = remAss.at(nonmuId);
                    art::Ptr<rb::Track> nonmutrk = tracks[nonmuId];
                    double vertE = 0.0;
                    if(nontrkHits.NCell() > 0){
                      vertE = nontrkHits.TotalGeV();
                    }
                    offE = vertE/(mutrk->TotalGeV()+nonmutrk->TotalGeV());
                    dedx = nonmurem->AvgDedx()/murem->AvgDedx();
                    
                    // gather all pid variables
                    std::vector<float> knnInputs(4);
                    knnInputs[0] = offE;
                    knnInputs[1] = ediff;
                    knnInputs[2] = ediffZ;
                    knnInputs[3] = dedx;
                    double kNN = knn2trk->EvaluateMVA(knnInputs,"KNN method");
                    
                    // update qe pid info
                    qepid.SetNtrk(2);
                    qepid.SetOffTrkE(offE);
                    qepid.SetEdiff(ediff);
                    qepid.SetEdiffZ(ediffZ);
                    qepid.SetDedx(dedx);
                    qepid.SetVal(kNN);
                  } // end of 2 track case

                } // end of else nonsense variables
              } // check if we have numu energy for this slice
            } // check that numu energy association is valid
            
          } // check if we have the correct number of tracks for qe piding
                  
        } // check if we found the best muon
        


      } // end of if < 3 total tracks
      
      // add this pid object to the vector of ones to add to the event
      qeefOut->push_back(qepid);
      // make an association to the slice it came from
      util::CreateAssn(*this,evt,*(qeefOut.get()),slice,*(qeAss.get()));

    } // end loop over slices
    
    evt.put(std::move(qeefOut));
    evt.put(std::move(qeAss));

  } // end of produce

} // end of qeef namespace

namespace qeef
{
  DEFINE_ART_MODULE(QeFinder)
}