WCTrackRecoAnalysis_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
/// \file    WCTrackRecoAnalysis_module.cc
/// \module  analyzer
/// \brief   Analyze and benchmark various reconstructions for the
///          ToF reconstruction.
/// \author  Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
///          David Duenas (University of Cincinnati)
/// \date    May 2019
////////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/FindManyP.h"

// nova
#include "RawData/RawBeamline.h"
#include "BeamlineRecoBase/BeamlineDigit.h"
#include "BeamlineRecoBase/BeamlinePID.h"
#include "BeamlineRecoBase/MuonStackTrack.h"
#include "BeamlineRecoBase/ToF.h"
#include "BeamlineRecoBase/WCTrack.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// root
#include "TH1F.h"

// stl
#include <iostream>

// -----------------------------------------------------------------------
namespace beamlinerecoana {

  class WCTrackRecoAnalysis : public art::EDAnalyzer {

  public:

    WCTrackRecoAnalysis(const fhicl::ParameterSet& pset);

    void reconfigure(const fhicl::ParameterSet& pset);
    void analyze(const art::Event& evt);

  private:

    std::string fWCTrackModuleLabel;
    std::string fMCTruthModuleLabel;

    TH1F* hRecoMomentum;

    art::ServiceHandle<art::TFileService> tfs;

  };

  DEFINE_ART_MODULE(WCTrackRecoAnalysis)

}

// -----------------------------------------------------------------------
beamlinerecoana::WCTrackRecoAnalysis::WCTrackRecoAnalysis(const fhicl::ParameterSet& pset) : EDAnalyzer(pset) {
  this->reconfigure(pset);
  hRecoMomentum = tfs->make<TH1F>("RecoMomentum", ";Reconstructed WC Momentum (MeV);", 100, 0, 2000);
}

// -----------------------------------------------------------------------
void beamlinerecoana::WCTrackRecoAnalysis::reconfigure(const fhicl::ParameterSet& pset) {
  fWCTrackModuleLabel = pset.get<std::string>("WCTrackModuleLabel");
  fMCTruthModuleLabel = pset.get<std::string>("fMCTruthModuleLabel");
}

// -----------------------------------------------------------------------
void beamlinerecoana::WCTrackRecoAnalysis::analyze(const art::Event& evt) {

  std::cout << "Processing Event " << evt.event() << std::endl;

  // get the wc tracks from the event
  art::Handle<std::vector<brb::WCTrack> > wcTrackHandle;
  std::vector<art::Ptr<brb::WCTrack> > wcTracks;
  if (evt.getByLabel(fWCTrackModuleLabel, wcTrackHandle))
    art::fill_ptr_vector(wcTracks, wcTrackHandle);

  // print stuff
  std::cout << "There are " << wcTracks.size() << " WC tracks in this event" << std::endl;

  if (wcTracks.size()) {
    std::cout << "Track momentum is " << wcTracks[0]->Momentum() << std::endl;
    hRecoMomentum->Fill(wcTracks[0]->Momentum());
  }

  // get the generator-level information from the event
  art::Handle<std::vector<simb::MCTruth> > mcTruthHandle;
  std::vector<art::Ptr<simb::MCTruth> > mcTruths;
  if (evt.getByLabel(fMCTruthModuleLabel, mcTruthHandle))
    art::fill_ptr_vector(mcTruths, mcTruthHandle);

  // print more stuff
  std::cout << "Thers is " << mcTruths.size() << " MCTruth simulation-level objects in this event" << std::endl;

  if (mcTruths.size()) {
    std::cout << "There are " << mcTruths[0]->NParticles() << " particles in this simulation" << std::endl;
    if (mcTruths[0]->NParticles()) {
      const simb::MCParticle& particle = mcTruths[0]->GetParticle(0);
      std::cout << "Particle momentum " << particle.Momentum().Vect().Mag() << std::endl;
      std::cout << "Particle momentum " << particle.P() << std::endl;
    }
  }

  // make plots -- 2D reco vs true momentum, any others?
  // be aware there are sometimes multiple particles in the event (there will (I think) always be one MCTruth, but it may have NParticles() == 2)

  return;

}