PFParticleTrackAna_module.cc

Go to the documentation of this file.

 
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
 
#include "TTree.h"
 
#include <string>
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
namespace lar_pandora {
 
  class PFParticleTrackAna : public art::EDAnalyzer {
  public:
    PFParticleTrackAna(fhicl::ParameterSet const& pset);
 
    virtual ~PFParticleTrackAna();
 
    void beginJob();
    void endJob();
    void analyze(const art::Event& evt);
    void reconfigure(fhicl::ParameterSet const& pset);
 
  private:
    TTree* m_pCaloTree; 
 
    int m_run;     
    int m_event;   
    int m_index;   
    int m_ntracks; 
    int m_trkid;   
    int m_plane;   
 
    double m_length;        
    double m_dEdx;          
    double m_dNdx;          
    double m_dQdx;          
    double m_residualRange; 
 
    double m_x;  
    double m_y;  
    double m_z;  
    double m_px; 
    double m_py; 
    double m_pz; 
 
    bool m_useModBox; 
    bool m_isCheated; 
 
    std::string m_trackModuleLabel; 
  };
 
  DEFINE_ART_MODULE(PFParticleTrackAna)
 
} // namespace lar_pandora
 
//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows
 
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
 
#include "larcore/Geometry/Geometry.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardataobj/RecoBase/Track.h"
 
#include "larpandora/LArPandoraInterface/LArPandoraHelper.h"
 
#include <iostream>
 
namespace lar_pandora {
 
  PFParticleTrackAna::PFParticleTrackAna(fhicl::ParameterSet const& pset) : art::EDAnalyzer(pset)
  {
    this->reconfigure(pset);
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  PFParticleTrackAna::~PFParticleTrackAna() {}
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleTrackAna::reconfigure(fhicl::ParameterSet const& pset)
  {
    m_useModBox = pset.get<bool>("UeModBox", true);
    m_isCheated = pset.get<bool>("IsCheated", false);
    m_trackModuleLabel = pset.get<std::string>("TrackModule", "pandora");
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleTrackAna::beginJob()
  {
    //
    art::ServiceHandle<art::TFileService const> tfs;
 
    m_pCaloTree = tfs->make<TTree>("calorimetry", "LAr Track Calo Tree");
    m_pCaloTree->Branch("run", &m_run, "run/I");
    m_pCaloTree->Branch("event", &m_event, "event/I");
    m_pCaloTree->Branch("index", &m_index, "index/I");
    m_pCaloTree->Branch("ntracks", &m_ntracks, "ntracks/I");
    m_pCaloTree->Branch("trkid", &m_trkid, "trkid/I");
    m_pCaloTree->Branch("plane", &m_plane, "plane/I");
    m_pCaloTree->Branch("length", &m_length, "length/D");
    m_pCaloTree->Branch("dEdx", &m_dEdx, "dEdx/D");
    m_pCaloTree->Branch("dNdx", &m_dNdx, "dNdx/D");
    m_pCaloTree->Branch("dQdx", &m_dQdx, "dQdx/D");
    m_pCaloTree->Branch("residualRange", &m_residualRange, "residualRange/D");
    m_pCaloTree->Branch("x", &m_x, "x/D");
    m_pCaloTree->Branch("y", &m_y, "y/D");
    m_pCaloTree->Branch("z", &m_z, "z/D");
    m_pCaloTree->Branch("px", &m_px, "px/D");
    m_pCaloTree->Branch("py", &m_py, "py/D");
    m_pCaloTree->Branch("pz", &m_pz, "pz/D");
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleTrackAna::endJob() {}
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleTrackAna::analyze(const art::Event& evt)
  {
    std::cout << " *** PFParticleTrackAna::analyze(...) *** " << std::endl;
 
    m_run = evt.run();
    m_event = evt.id().event();
    m_index = 0;
 
    m_ntracks = 0;
    m_trkid = 0;
    m_plane = 0;
    m_length = 0.0;
    m_dEdx = 0.0;
    m_dNdx = 0.0;
    m_dQdx = 0.0;
    m_residualRange = 0.0;
 
    m_x = 0.0;
    m_y = 0.0;
    m_z = 0.0;
    m_px = 0.0;
    m_py = 0.0;
    m_pz = 0.0;
 
    std::cout << "  Run: " << m_run << std::endl;
    std::cout << "  Event: " << m_event << std::endl;
 
    TrackVector trackVector;
    TracksToHits tracksToHits;
    LArPandoraHelper::CollectTracks(evt, m_trackModuleLabel, trackVector, tracksToHits);
 
    std::cout << "  Tracks: " << trackVector.size() << std::endl;
 
    // art::ServiceHandle<geo::Geometry const> theGeometry;
    // auto const* theDetector = lar::providerFrom<detinfo::DetectorPropertiesService>();
 
 
    // const double adc2eU(5.1e-3);
    // const double adc2eV(5.2e-3);
    // const double adc2eW(5.4e-3);
    // const double adc2eCheat(theDetector->ElectronsToADC());
 
    // const double tau(theDetector->ElectronLifetime());
 
    m_ntracks = trackVector.size();
 
    for (TrackVector::const_iterator iter = trackVector.begin(), iterEnd = trackVector.end();
         iter != iterEnd;
         ++iter) {
      const art::Ptr<recob::Track> track = *iter;
 
      m_trkid = track->ID();
      m_length = track->Length();
 
      m_plane = 0;
      m_dEdx = 0.0;
      m_dNdx = 0.0;
      m_dQdx = 0.0;
      m_residualRange = 0.0;
 
      m_x = 0.0;
      m_y = 0.0;
      m_z = 0.0;
      m_px = 0.0;
      m_py = 0.0;
      m_pz = 0.0;
 
      for (unsigned int p = 0; p < track->NumberTrajectoryPoints(); ++p) {
        auto pos = track->LocationAtPoint(p);
        auto dir = track->DirectionAtPoint(p);
 
        m_residualRange = track->Length(p);
 
        m_x = pos.x();
        m_y = pos.y();
        m_z = pos.z();
        m_px = dir.x();
        m_py = dir.y();
        m_pz = dir.z();
 
        /*************************************************************/
        /*                          WARNING                          */
        /*************************************************************/
        /* The dQdx information in recob::Track has been deprecated  */
        /* since 2016 and in 11/2018 the recob::Track interface was  */
        /* changed and DQdxAtPoint and NumberdQdx were removed.      */
        /* Therefore the code below is now commented out             */
        /* (note that it was most likely not functional anyways).    */
        /* For any issue please contact: larsoft-team@fnal.gov       */
        /*************************************************************/
        /*
            const double dQdxU(track->DQdxAtPoint(p, geo::kU)); // plane 0
            const double dQdxV(track->DQdxAtPoint(p, geo::kV)); // plane 1
            const double dQdxW(track->DQdxAtPoint(p, geo::kW)); // plane 2
 
            m_plane = ((dQdxU > 0.0) ? geo::kU : (dQdxV > 0.0) ? geo::kV : geo::kW);
 
            const double adc2e(m_isCheated ? adc2eCheat : (geo::kU == m_plane) ? adc2eU : (geo::kV == m_plane) ? adc2eV : adc2eW);
 
            m_dQdx = ((geo::kU == m_plane) ? dQdxU : (geo::kV == m_plane) ? dQdxV : dQdxW);
 
            // TODO: Need to include T0 information (currently assume T0 = 0)
 
            m_dNdx = ((m_dQdx / adc2e) * exp((m_x / theDetector->GetXTicksCoefficient()) * theDetector->SamplingRate() * 1.e-3 / tau));
 
            m_dEdx = (m_useModBox ? theDetector->ModBoxCorrection(m_dNdx) : theDetector->BirksCorrection(m_dNdx));
        */
        /*************************************************************/
 
        m_pCaloTree->Fill();
        ++m_index;
      }
    }
  }
 
} //namespace lar_pandora