MCParticleTrackMatching_module.cc

Go to the documentation of this file.

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h" 
#include "canvas/Persistency/Common/Ptr.h" 
#include "canvas/Persistency/Common/PtrVector.h" 
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h" 
#include "art/Framework/Services/Optional/TFileDirectory.h"

#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include <memory>
#include <iostream>
#include <map>
#include <iterator>

// LArSoft
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"

#include "lardata/Utilities/AssociationUtil.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "lardataobj/AnalysisBase/BackTrackerMatchingData.h"

namespace t0 {
  class MCParticleTrackMatching;
}

class t0::MCParticleTrackMatching : public art::EDProducer {
public:
  explicit MCParticleTrackMatching(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  MCParticleTrackMatching(MCParticleTrackMatching const &) = delete;
  MCParticleTrackMatching(MCParticleTrackMatching &&) = delete;
  MCParticleTrackMatching & operator = (MCParticleTrackMatching const &) = delete; 
 MCParticleTrackMatching & operator = (MCParticleTrackMatching &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;

  // Selected optional functions.
  void beginJob() override;
  void reconfigure(fhicl::ParameterSet const & p) ;


private:
  
  art::InputTag fTrackModuleLabel;
  art::InputTag fTrackHitAssnLabel;
  art::InputTag fHitModuleLabel;
  art::InputTag fHitParticleAssnLabel;
  
};


t0::MCParticleTrackMatching::MCParticleTrackMatching(fhicl::ParameterSet const & p)
{
  reconfigure(p);
  produces< art::Assns<recob::Track , simb::MCParticle, anab::BackTrackerMatchingData > > ();
}

void t0::MCParticleTrackMatching::reconfigure(fhicl::ParameterSet const & p)
{
  fTrackModuleLabel = p.get<art::InputTag>("TrackModuleLabel");
  fTrackHitAssnLabel = p.get<art::InputTag>("TrackHitAssnLabel",fTrackModuleLabel);
  fHitModuleLabel = p.get<art::InputTag>("HitModuleLabel");
  fHitParticleAssnLabel = p.get<art::InputTag>("HitParticleAssnLabel");
  
}

void t0::MCParticleTrackMatching::beginJob()
{
}

void t0::MCParticleTrackMatching::produce(art::Event & evt)
{
  if (evt.isRealData()) return;

  //auto mcpartHandle = evt.getValidHandle< std::vector<simb::MCParticle> >("largeant");
  std::unique_ptr< art::Assns<recob::Track, simb::MCParticle, anab::BackTrackerMatchingData > > MCPartTrackassn( new art::Assns<recob::Track, simb::MCParticle, anab::BackTrackerMatchingData >);
  
  
  double maxe = -1;
  double tote = 0;
  // int    trkid = -1;
  //int    maxtrkid = -1;
  //double maxn = -1;
  //double totn = 0;
  //int maxntrkid = -1;

  anab::BackTrackerMatchingData btdata;
  std::unordered_map<int,double> trkide;
  
  art::Handle< std::vector<recob::Track> > trackListHandle;
  evt.getByLabel(fTrackModuleLabel,trackListHandle);
  
  art::Handle< std::vector<recob::Hit> > hitListHandle;
  evt.getByLabel(fHitModuleLabel,hitListHandle);

  if(!trackListHandle.isValid()){
    std::cerr << "Track handle is not valid!" << std::endl;
    return;
  }
  
  if(!hitListHandle.isValid()){
    std::cerr << "Hit handle is not valid!" << std::endl;
    return;
  }

  auto const& trackList(*trackListHandle);
  art::FindManyP<recob::Hit> fmtht(trackListHandle, evt, fTrackHitAssnLabel);
  //auto const& mcpartList(*mcpartHandle);

  for(size_t i_t=0; i_t<trackList.size(); ++i_t){
    art::Ptr<recob::Track> trkPtr(trackListHandle,i_t);
    trkide.clear();
    tote = 0; maxe=-1; art::Ptr<simb::MCParticle> maxp;
    
    std::vector< art::Ptr<recob::Hit> > allHits = fmtht.at(i_t);

    std::vector<anab::BackTrackerHitMatchingData const*> bthmd_vec;
    std::vector< art::Ptr<simb::MCParticle> > matchedParticlePtrs;
    
    art::FindManyP<simb::MCParticle,anab::BackTrackerHitMatchingData>
      particles_per_hit(hitListHandle, evt, fHitParticleAssnLabel);

    for(size_t i_h=0; i_h<allHits.size(); ++i_h){
      bthmd_vec.clear(); matchedParticlePtrs.clear();
      particles_per_hit.get(allHits[i_h].key(),matchedParticlePtrs,bthmd_vec);

      for(size_t i_p=0; i_p<matchedParticlePtrs.size(); ++i_p){
        trkide[ matchedParticlePtrs[i_p]->TrackId() ] += bthmd_vec[i_p]->energy;
        tote += bthmd_vec[i_p]->energy;
        if( trkide[ matchedParticlePtrs[i_p]->TrackId() ] > maxe ){
          maxe = trkide[ matchedParticlePtrs[i_p]->TrackId() ];
          maxp = matchedParticlePtrs[i_p];
        }
      }//end loop over particles per hit

    }//end loop over hits

    btdata.cleanliness = maxe/tote;
    if(maxe>0)
      MCPartTrackassn->addSingle(trkPtr, maxp, btdata);

  }//end loop over tracks

  evt.put(std::move(MCPartTrackassn));
} // Produce

DEFINE_ART_MODULE(t0::MCParticleTrackMatching)