NumuEAna_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
// Class:       NumuEAna
// Module Type: analyzer
// File:        NumuEAna_module.cc
////////////////////////////////////////////////////////////////////////

#include "TH1.h"
#include "TH2.h"
#include "TTree.h"

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art_root_io/TFileService.h"

#include "MCCheater/BackTracker.h"
#include "RunHistory/service/RunHistoryService.h"
#include "CosRej/CosRejObj.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Energy.h"
#include "Simulation/Simulation.h"
#include "Simulation/Particle.h"
#include "Geometry/Geometry.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"
#include "ReMId/ReMId.h"
#include "QEEventFinder/QePId.h"
#include "NumuEnergy/NumuE.h"
#include "NovaDAQConventions/DAQConventions.h"

class TH1D;
class TH2D;
class TTree;

namespace numue {
  
  class NumuEAna : public art::EDAnalyzer {
  public:
    explicit NumuEAna(fhicl::ParameterSet const &p);
    virtual ~NumuEAna();
    
    virtual void analyze(art::Event const &e);
    
    virtual void reconfigure(fhicl::ParameterSet const & p);
    virtual void beginJob();
    virtual void beginRun(art::Run const &r);
    double TotalLengthInDetector(const sim::Particle* part, geo::View_t view);
    
  private:
    
    std::string fHitModuleLabel;         ///< lable for module making the hits
    std::string fTrackModuleLabel;       ///< label for module making the tracks
    std::string fCosmicTrackModuleLabel; ///< label for module making cosmic tracks
    std::string fSlicerModuleLabel;      ///< label for module making the slices
    std::string fEnergyModuleLabel;      ///< label for module writing energy objects
    std::string fPIDModuleLabel;         ///< label for module writing pid objects
    std::string fCosRejLabel;            ///< label for module writing cosmic rejection
    std::string fQepidLabel;             ///< label for module writing qe pid objects

    art::ServiceHandle<nova::dbi::RunHistoryService> rh;
    
    TTree*      fNumuETree;           ///< n-tuple used for further analyses (pure population)
    TH1D*       fCalCCEnergy;         ///< histogram of calorimetric neutrino energy for general CC population
    TH1D*       fTrkQEEnergy;         ///< histogram of muon track length neutrino energy for QE super pure population
    TH1D*       fTrkNonQEEnergy;      ///< histogram of muon track length neutrino energy for NonQE super pure population
    TH1D*       fAngleQEEnergy;       ///< histogram of neutrino energy for QE population using angle formula
    
    TH1D*       fAngleQEError;        ///< histogram of error on the neutrino energy for QE population using angle formula
    
    TH1D*       fCalCCEnergyDiff;     ///< histogram of the diff. between calorimetric neutrino E and true neutrino E (for pure pop.)
    TH1D*       fTrkQEEnergyDiff;     ///< histogram of the diff. between muon track length neutrino E and true neutrino E (for pure QE pop.)
    TH1D*       fTrkNonQEEnergyDiff;  ///< histogram of the diff. between muon track length neutrino E and true neutrino E (for pure NonQE pop.)
    TH1D*       fAngleQEEnergyDiff;   ///< histogram of the diff. between neutrino E using angle formula and true neutrino E (for pure QE pop.)
    
    TH2D*       fCalCCEnergy2D;       ///< histogram of calorimetric neutrino E and true neutrino E (for pure pop.)
    TH2D*       fTrkQEEnergy2D;       ///< histogram of muon track length neutrino E and true neutrino E (for pure QE pop.)
    TH2D*       fTrkNonQEEnergy2D;    ///< histogram of muon track length neutrino E and true neutrino E (for pure NonQE pop.)
    TH2D*       fAngleQEEnergy2D;     ///< histogram of neutrino E using angle formula and true neutrino E (for pure QE pop.)
    
    TH1D*       fHadronCalEnergy;     ///< histogram of hadronic calorimetric energy NOT on the muon track
    TH1D*       fHadronTrkEnergy;     ///< histogram of hadronic calorimetric energy on the muon track
    
    TH1D*       fNDTrkLenActive;      ///< histogram of ND muon track length in the active region
    TH1D*       fNDTrkLenCatcher;     ///< histogram of ND muon track length in the muon catcher region
    TH1D*       fNDTrkCalActive;      ///< histogram of ND muon calorimetric energy in the active region
    TH1D*       fNDTrkCalTransition;  ///< histogram of ND muon calorimetric energy in the transition region
    TH1D*       fNDTrkCalCatcher;     ///< histogram of ND muon calorimetric energy in the muon catcher region
    TH1D*       fNDHadCalActive;      ///< histogram of ND hadron calorimetric energy in the active region
    TH1D*       fNDHadCalTransition;  ///< histogram of ND hadron calorimetric energy in the transition region
    TH1D*       fNDHadCalCatcher;     ///< histogram of ND hadron calorimetric energy in the muon catcher region
    
    // Main energy estimators
    float CalCCEnergy   ; 
    float QEEnergy      ; 
    float NonQEEnergy   ; 
    float CCEnergy      ; 
    float AngleQEEnergy ; 
    float AngleQEError  ; 
    float HadCalEnergy  ; 
    float HadTrkEnergy  ; 
    float ndTrkLenAct   ; 
    float ndTrkLenCat   ; 
    float ndTrkCalAct   ; 
    float ndTrkCalTran  ; 
    float ndTrkCalCat   ; 
    float ndHadCalAct   ; 
    float ndHadCalTran  ; 
    float ndHadCalCat   ;
    float ndHadTrkAct   ; 
    float ndHadTrkTran  ; 
    float ndHadTrkCat   ;
    float ndTrkTranX    ;
    float ndTrkTranY    ;
      
    bool fND;                  ///< is detector ND?
    bool fFD;                  ///< is detector FD?

    bool isCC;                 ///< is the interaction CC?
    bool TrueQEInt;            ///< is the interaction QE?
    bool trackMatchesMuon;     ///< track is muon-like
    bool contained;            ///< does it pass numu containment cuts?

    int Run;
    int SubRun;
    int Evt;
    int SubEvt;
    int goodMuon;
    int bestTrack;
    int PDGCode;
    int nKalman;
    int nCosmic;
    int qepidntrk;              ///< Number of tracks used by QePId

    // All of these have to do with containment
    int firstplane;
    int lastplane;
    int nplanestofront;
    int nplanestoback;
    unsigned int mincellsfromedge;
    int nslicehit;
    int nslicecontplanes;

    double NuPurity;             ///< Neutrino selection purity
    double NuEfficiency;         ///< Neutrino selection efficiency
    double trackLength;          ///< Reconstructed tracklength
    double catcherE;             ///< Energy before entering the muon catcher (ND)
    double bestPID;              ///< Highest RemID
    double qepidPID;             ///< PID value of output kNN from QePId
    double cosrejanglekal;       ///< Cosine of angle of best ReMId Kalman Track from CosRej
    double cosrejnumucontpid;    ///< Cosmic rejection PID for contained events from CosRej
    double trackMuonEffHits;
    double trackMuonPurHits;
    double start_x_reco;
    double start_y_reco;
    double start_z_reco;
    double end_x_reco;
    double end_y_reco;
    double end_z_reco;
    double TrueMuonE;            ///< True energy of the muon
    double RecoMuonE;            ///< Reco energy of the muon
    double TrueNuE;              ///< True energy of the neutrino
    double TrueTrackLength;      ///< True tracklength of the muon
    double fMaxTrkLen;           ///< length of detector from corner to corner
  };
  
  //-------------------------------------------------------------------
  NumuEAna::NumuEAna(fhicl::ParameterSet const &p):
    EDAnalyzer(p),
    fND(false),
    fFD(false),
    fMaxTrkLen(-1)
  {
    this->reconfigure(p);
  }
  
  //-------------------------------------------------------------------
  NumuEAna::~NumuEAna() 
  {
  }
  
  //-------------------------------------------------------------------
  void NumuEAna::analyze(art::Event const &e) 
  { 
    if(!(fND || fFD)){
      mf::LogWarning("NumuEAna") << "attempting to run on detector that is not "
                                 << "nd or fd, bail";
      return;
    }

    art::ServiceHandle<geo::Geometry> geom;
    
    Run = e.run();
    SubRun = e.subRun();
    Evt    = e.event();


    // get all hits in the event 
    art::Handle< std::vector<rb::CellHit> > hithdl;
    e.getByLabel(fHitModuleLabel, hithdl);
    
    
    art::PtrVector<rb::CellHit> allhits;
    for(size_t h = 0; h < hithdl->size(); ++h){
      art::Ptr<rb::CellHit> hit(hithdl, h);
      allhits.push_back(hit);
    }
    
    // get the tracks out of the eventpdg
    art::Handle<std::vector<rb::Track> > trkcol;
    e.getByLabel(fTrackModuleLabel,trkcol);
    
    //Load the slicer list from the event
    art::Handle<std::vector<rb::Cluster> > slicecol;
    e.getByLabel(fSlicerModuleLabel,slicecol);
    if(slicecol->empty()){
      mf::LogWarning ("No Slices")<<"No Slices in the input file";
      return;
    }
    
    art::PtrVector<rb::Cluster>  slicelist;
    for(unsigned int i = 0; i<slicecol->size();++i){
      art::Ptr<rb::Cluster> clust(slicecol,i);
      slicelist.push_back(clust);
    }
    
    art::ServiceHandle<cheat::BackTracker> bt;
    if(!bt->HaveTruthInfo()){
      mf::LogWarning("NumuEAna") << "attempting to run MC truth check on "
                                 << "things without truth, bail";
      return;
    }
    
    //Associating slices to tracks and energies
    art::FindManyP<rb::Track>   sliceToTracks(slicecol, e, fTrackModuleLabel);
    art::FindManyP<rb::Track>   sliceToCosmicTracks(slicecol, e, fCosmicTrackModuleLabel);
    art::FindOneP<cosrej::CosRejObj>  sliceToCosRej(slicecol, e, fCosRejLabel);
    art::FindOneP<NumuE> sliceToEnergy(slicecol, e, fEnergyModuleLabel);
    art::FindOneP<qeef::QePId> sliceToQepid(slicecol, e, fQepidLabel);
    
    //Loop over slices
    for(size_t s = 0; s < slicelist.size(); ++s){

      SubEvt = s;
  
      const art::Ptr<NumuE> sliceEnergy = sliceToEnergy.at(s);
      if (!sliceEnergy) continue;
      
      const art::Ptr<cosrej::CosRejObj> sliceCosrej = sliceToCosRej.at(s);
      if (!sliceCosrej) continue;

      cosrejanglekal = sliceCosrej->AngleKal();
      cosrejnumucontpid = sliceCosrej->ConCosPID();

      const art::Ptr<qeef::QePId> sliceQepid = sliceToQepid.at(s);
      qepidPID  = -5;
      qepidntrk = -5;
      if (sliceQepid){
        qepidPID  = sliceQepid->Value();
        qepidntrk = sliceQepid->Ntrk();
      }

      // A great exercise of copy and paste starts here
      unsigned int dibmask = (rh->GoodDiBlockMask()&rh->GetConfiguration());
      if(dibmask == ((1<<16)-1)) {
        dibmask = rh->GetConfiguration();
      }
      if(dibmask == 0) {
        dibmask = rh->GetConfiguration();
      }

      unsigned int dibfirst = 0;
      unsigned int diblast = 0;
      if (dibmask) {
        int iD;
        iD=0; while (!((dibmask>>iD)&1)) iD++; dibfirst=iD+1;
        iD=0; while (dibmask>>iD)        iD++; diblast=iD;
      }

      unsigned int cellsfromedge = 99999999;
      mincellsfromedge = 99999999;

      nslicehit = slicelist.at(s)->NCell();
      nslicecontplanes = slicelist.at(s)->MostContiguousPlanes(geo::kXorY);

      //Loop over slice in question here:
     for(unsigned int hitIdx = 0; hitIdx < slicelist.at(s)->NCell(); ++hitIdx)
     {
       const art::Ptr<rb::CellHit>& chit = slicelist.at(s)->Cell(hitIdx);
       const int planeNum = chit->Plane();
       cellsfromedge = std::min((unsigned int)chit->Cell(), geom->Plane(planeNum)->Ncells() - 1 - chit->Cell());

       if(cellsfromedge < mincellsfromedge)
         mincellsfromedge = cellsfromedge;

     } //END OF SLICE LOOP HERE

     // And ends here (lo and behold)
     
      //Pulling out all info from the energy object and plotting it - these plots are without any cuts
      CalCCEnergy   = sliceEnergy->CalCCE();
      QEEnergy      = sliceEnergy->TrkQEE();
      NonQEEnergy   = sliceEnergy->TrkNonQEE();
      CCEnergy      = sliceEnergy->TrkCCE();
      AngleQEEnergy = sliceEnergy->AngleQEE();
      AngleQEError  = sliceEnergy->AngleQEError();
      HadCalEnergy  = sliceEnergy->HadCalE();
      HadTrkEnergy  = sliceEnergy->HadTrkE();
      ndTrkLenAct   = sliceEnergy->NDTrkLenAct();
      ndTrkLenCat   = sliceEnergy->NDTrkLenCat();
      ndTrkCalAct   = sliceEnergy->NDTrkCalAct();
      ndTrkCalTran  = sliceEnergy->NDTrkCalTran();
      ndTrkCalCat   = sliceEnergy->NDTrkCalCat();
      ndHadCalAct   = sliceEnergy->NDHadCalAct();
      ndHadCalTran  = sliceEnergy->NDHadCalTran();
      ndHadCalCat   = sliceEnergy->NDHadCalCat();
      ndHadTrkAct   = sliceEnergy->NDHadTrkAct();
      ndHadTrkTran  = sliceEnergy->NDHadTrkTran();
      ndHadTrkCat   = sliceEnergy->NDHadTrkCat();
      ndTrkTranX    = sliceEnergy->NDTrkTranX();
      ndTrkTranY    = sliceEnergy->NDTrkTranY();
   
      fCalCCEnergy->Fill(CalCCEnergy);
      fTrkQEEnergy->Fill(QEEnergy);
      fTrkNonQEEnergy->Fill(NonQEEnergy);
      fAngleQEEnergy->Fill(AngleQEEnergy);
      fAngleQEError->Fill(AngleQEError);
      fHadronCalEnergy->Fill(HadCalEnergy);
      fHadronTrkEnergy->Fill(HadTrkEnergy);
      fNDTrkLenActive->Fill(ndTrkLenAct);
      fNDTrkLenCatcher->Fill(ndTrkLenCat);
      fNDTrkCalActive->Fill(ndTrkCalAct);
      fNDTrkCalTransition->Fill(ndTrkCalTran);
      fNDTrkCalCatcher->Fill(ndTrkCalCat);
      fNDHadCalActive->Fill(ndHadCalAct);
      fNDHadCalTransition->Fill(ndHadCalTran);
      fNDHadCalCatcher->Fill(ndHadCalCat);
            
      //Now start trying to make super pure samples with truth to fill resolution plots for neutrino energy  
      art::PtrVector< rb::CellHit > sliceHits;
      sliceHits = slicelist[s]->AllCells();
      
      std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(sliceHits,allhits);
      if (funcReturn.empty()) continue;
      
      NuPurity = funcReturn[0].purity;
      NuEfficiency = funcReturn[0].efficiency;
      
       //Getting CC neutrinos     
      const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
      isCC = (test_neutrino.CCNC()==simb::kCC);
      PDGCode = test_neutrino.Nu().PdgCode();
      
      //Checking if QE or not
      TrueQEInt = false;
      if ((test_neutrino.InteractionType()==simb::kCCQE)) TrueQEInt = true;
      
      //Determining information about the muon - does it have enough hits to make this a decent interaction?
      std::vector<const sim::Particle*> particleList = bt->MCTruthToParticles(funcReturn[0].neutrinoInt);
      std::set<int> neutrinoTrackIDs;
      for (unsigned int i = 0; i < particleList.size(); ++i){
        neutrinoTrackIDs.insert(particleList[i]->TrackId());
      }
      
      std::set<int>::iterator itr = neutrinoTrackIDs.begin();
      
      goodMuon        = -1;
      bestTrack       = -1;
      trackLength     = 0.0;
      catcherE        = 0.0;
      TrueTrackLength = 0.0;
      TrueMuonE       = 0.0;

      while( itr != neutrinoTrackIDs.end() ){
        
        int id_int = *itr;
        const sim::Particle* particle = bt->TrackIDToParticle(id_int); 
        int PDG = particle->PdgCode();
        
        if ((PDG != 13)&&(PDG != -13)){  
          ++itr;
          continue;
        }//End of loop over non-muons
        
        //Trying to find main muon for the event
        const sim::Particle* mother = bt->TrackIDToMotherParticle(id_int);
        if (mother->TrackId() != id_int){
          ++itr;
          continue;
        }//End of loop over muons who don't come from the neutrino
        
        bool passCuts = bt->PassMuonCuts(id_int,sliceHits);
        
        if (!passCuts){
          ++itr;
          continue;
        }//End of loop over muons who don't pass cuts
        
        goodMuon = id_int;
        TrueTrackLength = TotalLengthInDetector(particle,geo::kXorY);
        TrueMuonE = particle->E();

        double closeToCentre = 50.0;

        std::vector<sim::FLSHit> flshits = bt->ParticleToFLSHit(particle->TrackId());

        for(unsigned int q = 0; q < flshits.size(); ++q){
          if((flshits[q].GetPlaneID() == 191) && (std::abs(flshits[q].GetPDG())==13)){
            if(fabs(flshits[q].GetEntryY()) < closeToCentre){
              catcherE = flshits[q].GetEntryEnergy() + 0.105658;
              closeToCentre = fabs(flshits[q].GetEntryY());
            }
            if(fabs(flshits[q].GetExitY()) < closeToCentre){
              catcherE = flshits[q].GetExitEnergy() + 0.105658;
              closeToCentre = fabs(flshits[q].GetExitY());
            }
          } // End of loop over plane 191
        } // End of loop over flshits

        break;      
      }//End of while loop over neutrino particles
      
      if (goodMuon == -1) continue;
      
      //Now we know we have an interaction of the correct type, passes muon cuts. We now loop over tracks.
      
    
      // Get all of the tracks in the slice
      const std::vector< art::Ptr<rb::Track> > sliceTracks = sliceToTracks.at(s);    
      if(sliceTracks.empty()) continue; 

      // Get all of the cosmic tracks in the slice
      const std::vector< art::Ptr<rb::Track> > sliceCosmicTracks = sliceToCosmicTracks.at(s);  
      nCosmic = sliceCosmicTracks.size();

      bestPID    = 0.0;
      RecoMuonE  = -5.0;
      
      art::FindOneP<remid::ReMId> trackRemid(sliceTracks, e, fPIDModuleLabel);  
      art::FindOneP<rb::Energy>   trackEnergy(sliceTracks, e, fEnergyModuleLabel); 
      
      // Get Reco values for the best track
      bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
      
      if(bestTrack == -1) continue;
      if(bestTrack == 999) continue;

      trackLength = sliceTracks[bestTrack]->TotalLength();
      bestPID = trackRemid.at(bestTrack)->Value();

      if(trackEnergy.isValid()){
        RecoMuonE = trackEnergy.at(bestTrack)->E();
      }
      
      //Checking if reco muon track matches truth muon well enough
      trackMatchesMuon = false;
      std::set<int> muonTrackID;
      muonTrackID.insert(goodMuon);
      art::PtrVector< rb::CellHit > bestTrackHits;
      bestTrackHits = sliceTracks[bestTrack]->AllCells();  
      trackMuonEffHits = bt->HitCollectionEfficiency(muonTrackID, bestTrackHits, allhits, geo::kXorY, 0, false);  
      trackMuonPurHits = bt->HitCollectionPurity(muonTrackID, bestTrackHits, 0, 0, false); 
      
      if ((trackMuonEffHits > 0.8)&&(trackMuonPurHits > 0.8)) trackMatchesMuon = true;  
      
      //Determining if track is contained enough
      contained = false;  

      firstplane = -1;
      lastplane = -1;
      nplanestofront = -1;
      nplanestoback = -1;

      if (sliceTracks[bestTrack]->NCell()>0){
        
        firstplane = sliceTracks[bestTrack]->MinPlane();
        lastplane  = sliceTracks[bestTrack]->MaxPlane();
        if(fFD)
        {
          nplanestofront = firstplane - (dibfirst-1)*64;
          nplanestoback  = (diblast*64) - 1 - lastplane;
        }
        if(fND)
        {
          nplanestofront = firstplane;
          nplanestoback  = 214 - lastplane;
        }
        start_x_reco = sliceTracks[bestTrack]->Start().X();
        start_y_reco = sliceTracks[bestTrack]->Start().Y();
        start_z_reco = sliceTracks[bestTrack]->Start().Z();
        
        end_x_reco   = sliceTracks[bestTrack]->Stop().X();
        end_y_reco   = sliceTracks[bestTrack]->Stop().Y();
        end_z_reco   = sliceTracks[bestTrack]->Stop().Z();
        
        if (fFD){
          if(mincellsfromedge > 1){
            if (sliceCosrej->CosFwdCell() > 0 && sliceCosrej->CosBakCell() > 0){
              if (sliceCosrej->KalFwdCell() > 10 && sliceCosrej->KalBakCell() > 10){
                if (nplanestofront > 1 && nplanestoback > 1){
                  contained = true;             
                }
              }
            }
          } 
        }//End of FD containment logic chain
        
        if (fND){
          if(mincellsfromedge > 1){
            if(sliceCosrej->KalFwdCellND() > 4 && sliceCosrej->KalBakCellND() > 8){
              if(firstplane > 1 && lastplane < 212){
                if(start_z_reco < 1150 && (end_z_reco < 1275 || sliceCosrej->KalYPosAtTrans() < 55)){
                  if( (ndHadCalTran + ndHadCalCat) < 0.03){
                    contained = true;
                  }
                }
              }
            }
          }
        }//End of ND containemnt logic chain       
      }//End of loop over best tracks to find start/stop
      
      //Filling remaining histograms with relevant super pure populations      
      TrueNuE = test_neutrino.Nu().E();
      
      fCalCCEnergyDiff->Fill((CalCCEnergy-TrueNuE)/TrueNuE);
      fCalCCEnergy2D->Fill(CalCCEnergy,TrueNuE);
      
      if (TrueQEInt){
        fTrkQEEnergyDiff  ->Fill((QEEnergy-TrueNuE)/TrueNuE);
        fAngleQEEnergyDiff->Fill((AngleQEEnergy-TrueNuE)/TrueNuE);  
        fTrkQEEnergy2D  ->Fill(QEEnergy,TrueNuE);
        fAngleQEEnergy2D->Fill(AngleQEEnergy,TrueNuE);
      }//End of pure QE slices
      else{
        fTrkNonQEEnergyDiff->Fill((NonQEEnergy-TrueNuE)/TrueNuE);
        fTrkNonQEEnergy2D->Fill(NonQEEnergy,TrueNuE);
      }//End of pure NonQE slices
    
    nKalman = sliceTracks.size();

    fNumuETree->Fill(); // Fill the tree for the pure sample  
    }//End of loop over slices
    return;
  }//End of analyze
  
  //-------------------------------------------------------------------
  void NumuEAna::reconfigure(fhicl::ParameterSet const & p) 
  {
    fHitModuleLabel         = p.get< std::string >("HitModuleLabel"        );
    fTrackModuleLabel       = p.get< std::string >("TrackModuleLabel"      );
    fCosmicTrackModuleLabel = p.get< std::string >("CosmicTrackModuleLabel");
    fSlicerModuleLabel      = p.get< std::string >("SlicerModuleLabel"     );
    fEnergyModuleLabel      = p.get< std::string >("EnergyModuleLabel"     );
    fPIDModuleLabel         = p.get< std::string >("PIDModuleLabel"        );
    fCosRejLabel            = p.get< std::string >("CosRejLabel"           );
    fQepidLabel             = p.get< std::string >("QepidLabel"            );
  }
  
  //-------------------------------------------------------------------
  void NumuEAna::beginJob() 
  {
    art::ServiceHandle<art::TFileService> tfs;
    
    fNumuETree = tfs->make<TTree>("NumuETree","Numu analysis Energy tree"); // Output tree
    fNumuETree->Branch("CalCCEnergy",&CalCCEnergy); 
    fNumuETree->Branch("QEEnergy",&QEEnergy); 
    fNumuETree->Branch("NonQEEnergy",&NonQEEnergy);
    fNumuETree->Branch("CCEnergy",&CCEnergy); 
    fNumuETree->Branch("AngleQEEnergy",&AngleQEEnergy); 
    fNumuETree->Branch("AngleQEError",&AngleQEError); 
    fNumuETree->Branch("HadCalEnergy",&HadCalEnergy); 
    fNumuETree->Branch("HadTrkEnergy",&HadTrkEnergy); 
    fNumuETree->Branch("ndTrkLenAct",&ndTrkLenAct); 
    fNumuETree->Branch("ndTrkLenCat",&ndTrkLenCat); 
    fNumuETree->Branch("ndTrkCalAct",&ndTrkCalAct); 
    fNumuETree->Branch("ndTrkCalTran",&ndTrkCalTran); 
    fNumuETree->Branch("ndTrkCalCat",&ndTrkCalCat); 
    fNumuETree->Branch("ndHadCalAct",&ndHadCalAct); 
    fNumuETree->Branch("ndHadCalTran",&ndHadCalTran); 
    fNumuETree->Branch("ndHadCalCat",&ndHadCalCat);
    fNumuETree->Branch("ndHadTrkAct",&ndHadTrkAct); 
    fNumuETree->Branch("ndHadTrkTran",&ndHadTrkTran); 
    fNumuETree->Branch("ndHadTrkCat",&ndHadTrkCat);
    fNumuETree->Branch("ndTrkTranX",&ndTrkTranX);
    fNumuETree->Branch("ndTrkTranY",&ndTrkTranY);
    fNumuETree->Branch("ND",&fND);
    fNumuETree->Branch("FD",&fFD);
    // Booleans
    fNumuETree->Branch("isCC",&isCC);
    fNumuETree->Branch("TrueQEInt",&TrueQEInt); 
    fNumuETree->Branch("trackMatchesMuon",&trackMatchesMuon);
    fNumuETree->Branch("contained",&contained);
    // Integers
    fNumuETree->Branch("Run",&Run);
    fNumuETree->Branch("SubRun",&SubRun);
    fNumuETree->Branch("Evt",&Evt);
    fNumuETree->Branch("SubEvt",&SubEvt);
    fNumuETree->Branch("PDGCode",&PDGCode);
    fNumuETree->Branch("goodMuon",&goodMuon);
    fNumuETree->Branch("bestTrack",&bestTrack);
    fNumuETree->Branch("nKalman",&nKalman);
    fNumuETree->Branch("nCosmic",&nCosmic);
    // Containment integers
    fNumuETree->Branch("firstplane",&firstplane);
    fNumuETree->Branch("lastplane",&lastplane);
    fNumuETree->Branch("nplanestofront",&nplanestofront);
    fNumuETree->Branch("nplanestoback",&nplanestoback);
    fNumuETree->Branch("mincellsfromedge",&mincellsfromedge);
    fNumuETree->Branch("nslicehit",&nslicehit);
    fNumuETree->Branch("nslicecontplanes",&nslicecontplanes);
    // Doubles
    fNumuETree->Branch("NuPurity",&NuPurity);
    fNumuETree->Branch("NuEfficiency",&NuEfficiency);
    fNumuETree->Branch("trackLength",&trackLength);
    fNumuETree->Branch("catcherE",&catcherE);
    fNumuETree->Branch("bestPID",&bestPID);
    fNumuETree->Branch("qepidPID",&qepidPID);
    fNumuETree->Branch("qepidntrk",&qepidntrk);
    fNumuETree->Branch("cosrejanglekal",&cosrejanglekal);
    fNumuETree->Branch("cosrejnumucontpid",&cosrejnumucontpid);
    fNumuETree->Branch("trackMuonEffHits",&trackMuonEffHits);
    fNumuETree->Branch("trackMuonPurHits",&trackMuonPurHits);
    fNumuETree->Branch("start_x_reco",&start_x_reco);
    fNumuETree->Branch("start_y_reco",&start_y_reco);
    fNumuETree->Branch("start_z_reco",&start_z_reco);
    fNumuETree->Branch("end_x_reco",&end_x_reco);
    fNumuETree->Branch("end_y_reco",&end_y_reco);
    fNumuETree->Branch("end_z_reco",&end_z_reco);
    fNumuETree->Branch("TrueNuE",&TrueNuE);
    fNumuETree->Branch("TrueMuonE",&TrueMuonE);
    fNumuETree->Branch("RecoMuonE",&RecoMuonE);
    fNumuETree->Branch("TrueTrackLength",&TrueTrackLength);

    // And finally the histograms
    fCalCCEnergy        = tfs->make<TH1D>("CalCCEnergy",        ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
    fTrkQEEnergy        = tfs->make<TH1D>("TrkQEEnergy",        ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
    fTrkNonQEEnergy     = tfs->make<TH1D>("TrkNonQEEnergy",     ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
    fAngleQEEnergy      = tfs->make<TH1D>("AngleQEEnergy",      ";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
    
    fAngleQEError       = tfs->make<TH1D>("AngleQEError",       ";Error on Neutrino Energy;Slices",200, -0.1, 5.0);
    
    fCalCCEnergyDiff    = tfs->make<TH1D>("CalCCEnergyDiff",    ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
    fTrkQEEnergyDiff    = tfs->make<TH1D>("TrkQEEnergyDiff",    ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
    fTrkNonQEEnergyDiff = tfs->make<TH1D>("TrkNonQEEnergyDiff", ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
    fAngleQEEnergyDiff  = tfs->make<TH1D>("AngleQEEnergyDiff",  ";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
    
    fCalCCEnergy2D      = tfs->make<TH2D>("CalCCEnergy2D",      ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
    fTrkQEEnergy2D      = tfs->make<TH2D>("TrkQEEnergy2D",      ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
    fTrkNonQEEnergy2D   = tfs->make<TH2D>("TrkNonQEEnergy2D",   ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
    fAngleQEEnergy2D    = tfs->make<TH2D>("AngleQEEnergy2D",    ";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
    
    fHadronCalEnergy    = tfs->make<TH1D>("HadronCalEnergy",    ";Reco Hadron Energy (GeV);Slices",200, -0.1, 10);
    fHadronTrkEnergy    = tfs->make<TH1D>("HadronTrkEnergy",    ";Track Hadron Energy (GeV);Slices",200, -0.1, 5);
    
    fNDTrkLenActive     = tfs->make<TH1D>("NDTrkLenActive",     ";Track Length (cm) in Active Region;Slices",200, 0., 1400);  
    fNDTrkLenCatcher    = tfs->make<TH1D>("NDTrkLenCatcher",    ";Track Length (cm) in Catcher Region;Slices",200, 0., 800);
    fNDTrkCalActive     = tfs->make<TH1D>("NDTrkCalActive",     ";Track Energy (GeV) in Active Region;Slices",200, 0., 10);
    fNDTrkCalTransition = tfs->make<TH1D>("NDTrkCalTransition", ";Track Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
    fNDTrkCalCatcher    = tfs->make<TH1D>("NDTrkCalCatcher",    ";Track Energy (GeV) in Catcher Region;Slices",200, 0., 2);
    fNDHadCalActive     = tfs->make<TH1D>("NDHadCalActive",     ";Hadron Energy (GeV) in Active Region;Slices",200, 0., 10);
    fNDHadCalTransition = tfs->make<TH1D>("NDHadCalTransition", ";Hadron Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
    fNDHadCalCatcher    = tfs->make<TH1D>("NDHadCalCatcher",    ";Hadron Energy (GeV) in Catcher Region;Slices",200, 0., 2);
    
  }//End of begin Run
  
  //-------------------------------------------------------------------
  void NumuEAna::beginRun(art::Run const &r) 
  {
    art::ServiceHandle<geo::Geometry>  geom;    
    double detHalfWidth   =  geom->DetHalfWidth();
    double detHalfHeight  =  geom->DetHalfHeight();
    double detLength      =  geom->DetLength();
    
    fMaxTrkLen = util::pythag(detHalfWidth*2, detHalfHeight*2, detLength);
    
    fND = false;
    fFD = false;
    
    novadaq::cnv::DetId detID = geom->DetId();
    if (detID == novadaq::cnv::kNEARDET) fND = true;
    if (detID == novadaq::cnv::kFARDET)  fFD = true;
    
  }//End of beginRun

  double NumuEAna::TotalLengthInDetector(const sim::Particle* p, geo::View_t view)
  {
  const unsigned int N = p->NumberTrajectoryPoints();

  if(N < 2){
    mf::LogWarning("RecoCheckAna") << "Particle has no trajectory points";
    return 0;
  }

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

  double dist = 0;

  for(unsigned int n = 0; n < N-1; ++n){
    TVector3 p0 = p->Position(n).Vect();
    TVector3 p1 = p->Position(n+1).Vect();
    if(p0.Z() < 0 || p1.Z() < 0 ||
       p0.Z() > geom->DetLength() || p1.Z() > geom->DetLength() ||
       fabs(p0.X()) > geom->DetHalfWidth() || fabs(p1.X()) > geom->DetHalfWidth() ||
       fabs(p0.Y()) > geom->DetHalfHeight() || fabs(p1.Y()) > geom->DetHalfHeight()) continue;

    if(view == geo::kX){
      p0.SetY(0);
      p1.SetY(0);
    }
    if(view == geo::kY){
      p0.SetX(0);
      p1.SetX(0);
    }
    dist += (p1-p0).Mag();
  }

  return dist;
} // end of TotalLengthInDetector
  
  //----------------------------------------------------------
  DEFINE_ART_MODULE(NumuEAna)
  
}//End of namespace