ReMIdValidate_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////           
/// \brief  A module for analyzing the ReMId objects in 
///         neutrino interactions    
/// \author Nicholas Raddatz - raddatz@physics.umn.edu
/// \date   September 2012 
//////////////////////////////////////////////////////////////////// 

#include <string>

// ROOT includes
#include "TVector3.h"
#include "TH1F.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"

// NOvA includes
#include "ReMId/ReMId.h"
#include "MCCheater/BackTracker.h"
#include "Calibrator/Calibrator.h"
#include "RecoBase/Track.h"
#include "Utilities/AssociationUtil.h"


namespace remid {
  /// A module to analyze remid objects produced by the RecoMuon module
  class ReMIdValidate : public art::EDAnalyzer {
  public:
    
    explicit ReMIdValidate(fhicl::ParameterSet const& pset);
    virtual ~ReMIdValidate();
    
    void beginJob();
    void analyze(art::Event const& evt);
    void endJob();
    static bool SortByHits(const art::Ptr<rb::Track>&  a, const art::Ptr<rb::Track>&  b);
    
  private:
    std::string fSliceLabel;         ///<Where to find reconstructed slices
    std::string fTrackLabel;         ///<Where to find recondtructed tracks
    std::string fPidLabel;           ///<Where to find the pids
    
    TH1F* fpid;                      ///<Pid value of every track
    TH1F* fpidnc;
    TH1F* fpidcc;
    TH1F* fpidcccut;
    TH1F* fnscat;
    TH1F* fndedx;
    TH1F* fcosscat;
    TH1F* favgdedx;
    TH1F* fdedx;
    TH1F* fpidhigh;
    TH1F* fpidnchigh;
    TH1F* fpidcchigh;
    TH1F* fnscathigh;
    TH1F* fndedxhigh;
    TH1F* fcosscathigh;
    TH1F* favgdedxhigh;
    TH1F* fdedxhigh;
    TH1F* fpidmu;
    TH1F* fnscatmu;
    TH1F* fndedxmu;
    TH1F* fcosscatmu;
    TH1F* favgdedxmu;
    TH1F* fdedxmu;
    TH1F* fpidmucont;
    TH1F* fnscatmucont;
    TH1F* fndedxmucont;
    TH1F* fcosscatmucont;
    TH1F* favgdedxmucont;
    TH1F* fdedxmucont;
    
    bool IsContained(art::Ptr<rb::Track> track);   ///< Function to determine whether or not a track is contained
    
  };
}


namespace remid{
  
  //.......................................................................
  
  ReMIdValidate::ReMIdValidate(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    fSliceLabel = (pset.get< std::string >("SliceLabel"));
    fTrackLabel = (pset.get< std::string >("TrackLabel"));
    fPidLabel   = (pset.get< std::string >("PidLabel"  ));
  }
  
  //.......................................................................
  
  ReMIdValidate::~ReMIdValidate()
  {
  }
  
  //.......................................................................

  void ReMIdValidate::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fpid             = tfs->make<TH1F>("fpid",";pid;NTracks",95,-0.1,1.1);
    fpidnc           = tfs->make<TH1F>("fpidnc",";pid;NTracks",95,-0.1,1.1);
    fpidcc           = tfs->make<TH1F>("fpidcc",";pid;NTracks",95,-0.1,1.1);
    fpidcccut        = tfs->make<TH1F>("fpidcccut",";pid;NTracks",95,-0.1,1.1);
    fnscat           = tfs->make<TH1F>("fnscat",";nscat;NTracks",301,-0.5,301);
    fndedx           = tfs->make<TH1F>("fndedx",";ndedx;NTracks",301,-0.5,301);
    fcosscat         = tfs->make<TH1F>("fcosscat",";cos #theta_{scat};",300,-1.1,1.1);
    favgdedx         = tfs->make<TH1F>("favgdedx",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
    fdedx            = tfs->make<TH1F>("fdedx",";dE/dx (GeV/cm);",500,0.0,0.01);
    fpidnchigh       = tfs->make<TH1F>("fpidnchigh",";pid;NTracks",95,-0.1,1.1);
    fpidcchigh       = tfs->make<TH1F>("fpidcchigh",";pid;NTracks",95,-0.1,1.1);
    fpidhigh         = tfs->make<TH1F>("fpidhigh",";pid;NTracks",95,-0.1,1.1);
    fnscathigh       = tfs->make<TH1F>("fnscathigh",";nscat;NTracks",301,-0.5,301);
    fndedxhigh       = tfs->make<TH1F>("fndedxhigh",";ndedx;NTracks",301,-0.5,301);
    fcosscathigh     = tfs->make<TH1F>("fcosscathigh",";cos #theta_{scat};",300,-1.1,1.1);
    favgdedxhigh     = tfs->make<TH1F>("favgdedxhigh",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
    fdedxhigh        = tfs->make<TH1F>("fdedxhigh",";dE/dx (GeV/cm);",500,0.0,0.01);
    fpidmu           = tfs->make<TH1F>("fpidmu",";pid;NTracks",95,-0.1,1.1);
    fnscatmu         = tfs->make<TH1F>("fnscatmu",";nscat;NTracks",301,-0.5,301);
    fndedxmu         = tfs->make<TH1F>("fndedxmu",";ndedx;NTracks",301,-0.5,301);
    fcosscatmu       = tfs->make<TH1F>("fcosscatmu",";cos #theta_{scat};",300,-1.1,1.1);
    favgdedxmu       = tfs->make<TH1F>("favgdedxmu",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
    fdedxmu          = tfs->make<TH1F>("fdedxmu",";dE/dx (GeV/cm);",500,0.0,0.01);
    fpidmucont       = tfs->make<TH1F>("fpidmucont",";pid;NTracks",95,-0.1,1.1);
    fnscatmucont     = tfs->make<TH1F>("fnscatmucont",";nscat;NTracks",301,-0.5,301);
    fndedxmucont     = tfs->make<TH1F>("fndedxmucont",";ndedx;NTracks",301,-0.5,301);
    fcosscatmucont   = tfs->make<TH1F>("fcosscatmucont",";cos #theta_{scat};",300,-1.1,1.1);
    favgdedxmucont   = tfs->make<TH1F>("favgdedxmucont",";average dE/dx (GeV/cm);NTracks",500,0.0,0.01);
    fdedxmucont      = tfs->make<TH1F>("fdedxmucont",";dE/dx (GeV/cm);",500,0.0,0.01);
  }
  
  //.......................................................................
  
  void ReMIdValidate::analyze(art::Event const& evt)
  {
    // get the slices out of the event
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSliceLabel,slicecol);
    
    // get the tracks out of the event
    art::Handle<std::vector<rb::Track> > trkcol;
    evt.getByLabel(fTrackLabel,trkcol);
    
    art::PtrVector<rb::Cluster> slicelist;
    art::PtrVector<rb::CellHit> Allevthits; // PtrVector of all of the hits in the event
    
    for(unsigned int i = 0; i<slicecol->size();++i){
      art::Ptr<rb::Cluster>slice(slicecol,i);
      slicelist.push_back(slice);
      for(unsigned int ihit = 0; ihit<slice->NCell();++ihit){
        Allevthits.push_back(slice->Cell(ihit));
      }
    }
    
    // get  assosciations between slices and tracks
    art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
    art::ServiceHandle<cheat::BackTracker> bt;
    
    // loop over all slices and find the most effiecient slice for a neutrino interaction
    double bestEff = -1;
    unsigned int mostEffSlice = 0;  // Index of the most efficient neutrino slice
    int intType = 6;  // This is a index to keep track of what type of interaction this is
    // get the neutrino interaction type and slice it came from 
    for(unsigned int iSlice = 0;iSlice<slicelist.size();++iSlice){
      // get the neutrino interaction from this slice
      std::vector<cheat::NeutrinoEffPur> nuint = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),Allevthits);
      if(nuint.size()>0 && nuint[0].efficiency>bestEff){
        bestEff = nuint[0].efficiency;
        mostEffSlice = iSlice;
        // get the neutrino
        const simb::MCNeutrino nu = nuint[0].neutrinoInt->GetNeutrino();
        if(nu.CCNC() == simb::kNC){ intType = 5; }
        else{
          if(nu.Nu().PdgCode() == 12 || nu.Nu().PdgCode() == -12){ intType = 4; }
          else if(nu.Mode() == 0){intType = 0;}
          else if(nu.Mode() == 1){intType = 1;}
          else if(nu.Mode() == 2){intType = 2;}
          else if(nu.Mode() == 3){intType = 3;}
        }
      }
    }
    
    // loop over the slices
    for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
      if(intType <0){ continue; }
      if(mostEffSlice != iSlice){ continue; }
      // get all of the tracks associated to this slice
      std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
      if(tracks.size() == 0){ continue; }
      art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
      double bestpid = -1.0; 
      int bestInd = -1;
      // find the true primary muon track id if a CC interaction
      std::set<int> muID;
      if(intType < 4){
        std::set<int> trackIDs = bt->GetTrackIDList();
        for(std::set<int>::iterator it= trackIDs.begin(); it != trackIDs.end(); ++it){
          // get the particle
          const sim::Particle* part = bt->TrackIDToParticle(*it);
          // check if this is the muon
          // std::cout<<"On particle id: "<<*it<<" with pdg: "<<part->PdgCode()<<" and mother: "<<part->Mother()<<std::endl;
          if((part->PdgCode() == 13 || part->PdgCode() == -13) && part->Mother() == 0){
            muID.insert(*it);
            break;
          }
        }
      }
      double bestmueff = -1.0;
      int muInd = -1;
      for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
        art::Ptr<remid::ReMId> pidobj = foids.at(itrack);
        fpid->Fill(pidobj->Value());
        if(intType < 4){
          fpidcc->Fill(pidobj->Value());
          // find if this track is the most efficiently reconstructed muon
          double eff = bt->HitCollectionEfficiency(muID,tracks[itrack]->AllCells(),slicelist[iSlice]->AllCells(),geo::kXorY);
          if(eff > bestmueff){
            bestmueff = eff;
            muInd = itrack;
          }
        }
        else{ fpidnc->Fill(pidobj->Value());}
        fnscat->Fill(pidobj->NScatters());
        for(unsigned int i = 0; i<pidobj->NScatters(); ++i){
          fcosscat->Fill(pidobj->CosScat(i));
        }
        favgdedx->Fill(pidobj->AvgDedx());
        fndedx->Fill(pidobj->NDedx());
        for(unsigned int i = 0; i<pidobj->NDedx(); ++i){
          fdedx->Fill(pidobj->Dedx(i));
        }
        if(pidobj->Value() > bestpid){
          bestInd = itrack;
          bestpid = pidobj->Value();
        }

      }

      // make plots for highest pid track
      if(bestInd <0){ continue; }
      art::Ptr<remid::ReMId> pidobj = foids.at(bestInd);
      fpidhigh->Fill(pidobj->Value());
      if(intType < 4){fpidcchigh->Fill(pidobj->Value()); }
      else{ fpidnchigh->Fill(pidobj->Value()); }
      fnscathigh->Fill(pidobj->NScatters());
      for(unsigned int i = 0; i<pidobj->NScatters(); ++i){
        fcosscathigh->Fill(pidobj->CosScat(i));
      }
      favgdedxhigh->Fill(pidobj->AvgDedx());
      fndedxhigh->Fill(pidobj->NDedx());
      for(unsigned int i = 0; i<pidobj->NDedx(); ++i){
        fdedxhigh->Fill(pidobj->Dedx(i));
      }


      // make plots for true mu track by efficiency
      if(intType >3){ continue; }
      if(muInd < 0){continue; }
      art::Ptr<remid::ReMId> mupidobj = foids.at(muInd);
      fpidmu->Fill(mupidobj->Value());
      fnscatmu->Fill(mupidobj->NScatters());
      for(unsigned int i = 0; i<mupidobj->NScatters(); ++i){
        fcosscatmu->Fill(mupidobj->CosScat(i));
      }
      favgdedxmu->Fill(mupidobj->AvgDedx());
      fndedxmu->Fill(mupidobj->NDedx());
      for(unsigned int i = 0; i<mupidobj->NDedx(); ++i){
        fdedxmu->Fill(mupidobj->Dedx(i));
      }

      if(IsContained(tracks[muInd])){
        fpidmucont->Fill(mupidobj->Value());
        fnscatmucont->Fill(mupidobj->NScatters());
        for(unsigned int i = 0; i<mupidobj->NScatters(); ++i){
          fcosscatmucont->Fill(mupidobj->CosScat(i));
        }
        favgdedxmucont->Fill(mupidobj->AvgDedx());
        fndedxmucont->Fill(mupidobj->NDedx());
        for(unsigned int i = 0; i<mupidobj->NDedx(); ++i){
          fdedxmucont->Fill(mupidobj->Dedx(i));
        }
      }

      // make plots for contained qe sample
      if(!IsContained(tracks[bestInd]) && !IsContained(tracks[0])){ continue; }
      if(tracks[bestInd]->View() != geo::kXorY){ continue; }
      if(intType != 0){continue; }
      fpidcccut->Fill(pidobj->Value());

    } // end loop over slices


  }

  //.......................................................................

  void ReMIdValidate::endJob()
  {
  }

  //.......................................................................

  bool ReMIdValidate::SortByHits(const art::Ptr<rb::Track>& a, const art::Ptr<rb::Track>& b)
  {
    return a->NCell()>b->NCell();
  }

  //.......................................................................

  bool ReMIdValidate::IsContained(art::Ptr<rb::Track> track)
  {
    double xmin,xmax,ymin,ymax,zmin,zmax;
    // // define the boundaries of containment
    // int concells = 3;
    // art::ServiceHandle<geo::Geometry>  geom;    
    // double detHalfWidth   =  geom->DetHalfWidth();
    // double detHalfHeight  =  geom->DetHalfHeight();
    // double detLength =  geom->DetLength();
    // // get the boundaries of the containment boundary
    // const geo::PlaneGeo* pi = geom->Plane(concells-1);
    // const geo::PlaneGeo* pf = geom->Plane(geom->NPlanes()-concells);
    // // get zmin and zmax first
    // const geo::CellGeo* ci = pi->Cell(concells-1);
    // const geo::CellGeo* cl = pi->Cell(pi->Ncells()-1-concells);
    // const geo::CellGeo* cf = pf->Cell(0);
    // // double xyzf[3];
    // double xyzi[3], xyzf[3], xyzl[3];
    // ci->GetCenter(xyzi);
    // cl->GetCenter(xyzl);
    // cf->GetCenter(xyzf);
    // std::cout<<"3 cell z center position: "<<xyzf[2]<<std::endl;
    // double dzi = ci->HalfD();
    // zmin = xyzi[2]+dzi;

    // double dzf = cf->HalfD();
    // zmax = xyzf[2]-dzf;

    // if(pi->View() == geo::kX){
    //   xmin = xyzi[0]+ci->HalfW();
    //   xmax = xyzl[0]-cl->HalfW();
    //   int nextplane = geom->NextPlaneOtherView(concells-1);
    //   pi = geom->Plane(nextplane);
    //   ci = pi->Cell(concells-1);
    //   cl = pi->Cell(pi->Ncells()-1-concells);
    //   ci->GetCenter(xyzi);
    //   cl->GetCenter(xyzl);
    //   ymin = xyzi[1]+ci->HalfW();
    //   ymax = xyzl[1]-ci->HalfW();
    // }
    // else if(pi->View() == geo::kY){
    //   ymin = xyzi[1]+ci->HalfW();
    //   ymax = xyzl[1]-cl->HalfW();
    //   int nextplane = geom->NextPlaneOtherView(concells-1);
    //   pi = geom->Plane(nextplane);
    //   ci = pi->Cell(concells-1);
    //   cl = pi->Cell(pi->Ncells()-1-concells);
    //   ci->GetCenter(xyzi);
    //   cl->GetCenter(xyzl);
    //   xmin = xyzi[0]+ci->HalfW();
    //   xmax = xyzl[0]-ci->HalfW();
    // }

    // contianement cuts for a 28 block far detector. Track ends must be 3 cells + 100 cm from detector edge
    xmin = -750.213;
    xmax = 750.213;
    ymin = -750.213;
    ymax = 750.213;
    zmin = 19.9156;
    zmax = 6201.02; // 30 blocks
    zmax = 5945.29; // 28 blocks
    // zmax = detLength-zmin;
    TVector3 start = track->Start();
    TVector3 stop = track->Stop();
    if(start.X() > xmin+100.00 && start.X() < xmax-100.00 && start.Y() > ymin+100.00 && start.Y() < ymax-100.00 && start.Z() > zmin+100.00 && start.Z() < zmax-100.00 &&
       stop.X() > xmin && stop.X() < xmax && stop.Y() > ymin && stop.Y() < ymax && stop.Z() > zmin && stop.Z() < zmax){ return true; }
    else{ return false; }

  }


} // end of remid namespace

namespace remid
{
  DEFINE_ART_MODULE(ReMIdValidate)
}