TrackKalmanCheater_module.cc

Go to the documentation of this file.

// Class:       TrackKalmanCheater
// Module Type: producer
// File:        TrackKalmanCheater.h
//
// This class produces RecoBase/Track objects using KalmanFilterService.
// MC truth information is used to associate Hits used as input to
// the Kalman filter.
//
// Configuration parameters:
//
// Hist               - Histogram flag (generate histograms if true).
// UseClusterHits     - Use clustered hits if true, use all hits if false.
// HitModuleLabel     - Module label for unclustered Hits.
// ClusterModuleLabel - Module label for Clusters.
// MaxTcut            - Maximum delta ray energy in Mev for dE/dx.
// KalmanFilterAlg    - Parameter set for KalmanFilterAlg.
// SpacePointAlg      - Parmaeter set for space points.
//

#include <vector>
#include <deque>

#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h" 
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"

#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larsim/MCCheater/BackTrackerService.h"
#include "larsim/MCCheater/ParticleInventoryService.h"
#include "larreco/RecoAlg/KalmanFilterAlg.h"
#include "lardata/RecoObjects/KHitContainerWireX.h"
#include "lardata/RecoObjects/SurfYZPlane.h"
#include "lardata/RecoObjects/PropYZPlane.h"
#include "lardata/RecoObjects/KHit.h"
#include "larreco/RecoAlg/SpacePointAlg.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "nusimdata/SimulationBase/MCParticle.h"

#include "TH1F.h"

namespace trkf {

  class Propagator;

  class TrackKalmanCheater : public art::EDProducer {
  public:

    // Copnstructors, destructor.

    explicit TrackKalmanCheater(fhicl::ParameterSet const & pset);
    virtual ~TrackKalmanCheater();

    // Overrides.

    virtual void reconfigure(fhicl::ParameterSet const & pset);
    virtual void produce(art::Event & e);
    virtual void beginJob();
    virtual void endJob();

  private:

    // Fcl parameters.

    bool fHist;                        
    KalmanFilterAlg fKFAlg;            
    SpacePointAlg fSpacePointAlg;      
    bool fUseClusterHits;              
    std::string fHitModuleLabel;       
    std::string fClusterModuleLabel;   
    double fMaxTcut;                   

    const Propagator* fProp;

    // Histograms.

    TH1F* fHIncChisq;   
    TH1F* fHPull;       

    // Statistics.

    int fNumEvent;    
    int fNumTrack;    

  };

  DEFINE_ART_MODULE(TrackKalmanCheater)

} // namespace trkf

//------------------------------------------------------------------------------
trkf::TrackKalmanCheater::TrackKalmanCheater(fhicl::ParameterSet const & pset)
  : fHist(false)
  , fKFAlg(pset.get<fhicl::ParameterSet>("KalmanFilterAlg"))
  , fSpacePointAlg(pset.get<fhicl::ParameterSet>("SpacePointAlg"))
  , fUseClusterHits(false)
  , fMaxTcut(0.)
  , fProp(0)
  , fHIncChisq(0)
  , fHPull(0)
  , fNumEvent(0)
  , fNumTrack(0)
{
  reconfigure(pset);
  produces<std::vector<recob::Track>                   >();
  produces<std::vector<recob::SpacePoint>              >();
  produces<art::Assns<recob::Track, recob::Hit>        >();
  produces<art::Assns<recob::Track, recob::SpacePoint> >();
  produces<art::Assns<recob::SpacePoint, recob::Hit>   >();

  // Report.

  mf::LogInfo("TrackKalmanCheater") 
    << "TrackKalmanCheater configured with the following parameters:\n"
    << "  UseClusterHits = " << fUseClusterHits << "\n"
    << "  HitModuleLabel = " << fHitModuleLabel << "\n"
    << "  ClusterModuleLabel = " << fClusterModuleLabel;
}

//------------------------------------------------------------------------------
trkf::TrackKalmanCheater::~TrackKalmanCheater()
{}

//------------------------------------------------------------------------------
void trkf::TrackKalmanCheater::reconfigure(fhicl::ParameterSet const & pset)
{
  fHist = pset.get<bool>("Hist");
  fKFAlg.reconfigure(pset.get<fhicl::ParameterSet>("KalmanFilterAlg"));
  fSpacePointAlg.reconfigure(pset.get<fhicl::ParameterSet>("SpacePointAlg"));
  fUseClusterHits = pset.get<bool>("UseClusterHits");
  fHitModuleLabel = pset.get<std::string>("HitModuleLabel");
  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
  fMaxTcut = pset.get<double>("MaxTcut");
  if(fProp != 0)
    delete fProp;
  fProp = new PropYZPlane(fMaxTcut, true);
}

//------------------------------------------------------------------------------
void trkf::TrackKalmanCheater::beginJob()
{
  if(fHist) {

    // Book histograms.

    art::ServiceHandle<art::TFileService> tfs;
    art::TFileDirectory dir = tfs->mkdir("hitkalman", "TrackKalmanCheater histograms");

    fHIncChisq = dir.make<TH1F>("IncChisq", "Incremental Chisquare", 100, 0., 20.);
    fHPull = dir.make<TH1F>("Pull", "Hit Pull", 100, -10., 10.);
  }
}

//------------------------------------------------------------------------------
void trkf::TrackKalmanCheater::produce(art::Event & evt)
{
  ++fNumEvent;

  // Make a collection of tracks, plus associations, that will
  // eventually be inserted into the event.

  std::unique_ptr<std::vector<recob::Track> > tracks(new std::vector<recob::Track>);
  std::unique_ptr< art::Assns<recob::Track, recob::Hit> > th_assn(new art::Assns<recob::Track, recob::Hit>);
  std::unique_ptr< art::Assns<recob::Track, recob::SpacePoint> > tsp_assn(new art::Assns<recob::Track, recob::SpacePoint>);

  // Make a collection of space points, plus associations, that will
  // be inserted into the event.

  std::unique_ptr<std::vector<recob::SpacePoint> > spts(new std::vector<recob::SpacePoint>);
  std::unique_ptr< art::Assns<recob::SpacePoint, recob::Hit> > sph_assn(new art::Assns<recob::SpacePoint, recob::Hit>);

  // Make a collection of KGTracks where we will save our results.

  std::deque<KGTrack> kalman_tracks;

  // Get Services.

  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;

  // Reset space point algorithm.

  fSpacePointAlg.clearHitMap();


  // Get Hits.

  art::PtrVector<recob::Hit> hits;

  if(fUseClusterHits) {

    // Get clusters.

    art::Handle< std::vector<recob::Cluster> > clusterh;
    evt.getByLabel(fClusterModuleLabel, clusterh);

    // Get hits from all clusters.
    art::FindManyP<recob::Hit> fm(clusterh, evt, fClusterModuleLabel);

    if(clusterh.isValid()) {
      int nclus = clusterh->size();

      for(int i = 0; i < nclus; ++i) {
        art::Ptr<recob::Cluster> pclus(clusterh, i);
        std::vector< art::Ptr<recob::Hit> > clushits = fm.at(i);
        int nhits = clushits.size();
        hits.reserve(hits.size() + nhits);

        for(std::vector< art::Ptr<recob::Hit> >::const_iterator ihit = clushits.begin();
            ihit != clushits.end(); ++ihit) {
          hits.push_back(*ihit);
        }
      }
    }
  }
  else {

    // Get unclustered hits.

    art::Handle< std::vector<recob::Hit> > hith;
    evt.getByLabel(fHitModuleLabel, hith);
    if(hith.isValid()) {
      int nhits = hith->size();
      hits.reserve(nhits);

      for(int i = 0; i < nhits; ++i)
        hits.push_back(art::Ptr<recob::Hit>(hith, i));
    }
  }

  // Sort hits into separate PtrVectors based on track id.

  std::map<int, art::PtrVector<recob::Hit> > hitmap;

  // Loop over hits.

  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
        ihit != hits.end(); ++ihit) {
    //const recob::Hit& hit = **ihit;

    // Get track ids for this hit.

    std::vector<sim::TrackIDE> tids = bt_serv->HitToTrackIDEs(*ihit);

    // Loop over track ids.

    for(std::vector<sim::TrackIDE>::const_iterator itid = tids.begin();
        itid != tids.end(); ++itid) {
      int trackID = itid->trackID;

      // Add hit to PtrVector corresponding to this track id.

      hitmap[trackID].push_back(*ihit);
    }
  }

  // Extract geant mc particles.

  sim::ParticleList const& plist = pi_serv->ParticleList();

  // Loop over geant particles.

  {
    // use mf::LogDebug instead of LOG_DEBUG because we reuse it in many lines;
    // insertions are protected by mf::isDebugEnabled()
    mf::LogDebug log("TrackKalmanCheater");
    for(sim::ParticleList::const_iterator ipart = plist.begin();
        ipart != plist.end(); ++ipart) {
      const simb::MCParticle* part = (*ipart).second;
      if (!part) throw cet::exception("TrackKalmanCheater") << "no particle!\n";
      int pdg = part->PdgCode();

      // Ignore everything except stable charged nonshowering particles.

      int apdg = std::abs(pdg);
      if(apdg == 13 ||     // Muon
         apdg == 211 ||    // Charged pion
         apdg == 321 ||    // Charged kaon
         apdg == 2212) {   // (Anti)proton

        int trackid = part->TrackId();
        int stat = part->StatusCode();
        int nhit = 0;
        if(hitmap.count(trackid) != 0)
          nhit = hitmap[trackid].size();
        const TLorentzVector& pos = part->Position();
        const TLorentzVector& mom = part->Momentum();
        
        if (mf::isDebugEnabled()) {
          log << "Trackid=" << trackid 
              << ", pdgid=" << pdg 
              << ", status=" << stat
              << ", NHit=" << nhit << "\n"
              << "  x = " << pos.X()
              << ", y = " << pos.Y()
              << ", z = " << pos.Z() << "\n"
              << "  px= " << mom.Px()
              << ", py= " << mom.Py()
              << ", pz= " << mom.Pz() << "\n";
        } // if debugging
        
        double x = pos.X();
        double y = pos.Y();
        double z = pos.Z();
        double px = mom.Px();
        double py = mom.Py();
        double pz = mom.Pz();
        double p = std::sqrt(px*px + py*py + pz*pz);

        if(nhit > 0 && pz != 0.) {
          const art::PtrVector<recob::Hit>& trackhits = hitmap[trackid];
          if (trackhits.empty())
            throw cet::exception("TrackKalmanCheater") << "No hits in track\n";

          // Make a seed track (KTrack).

          std::shared_ptr<const Surface> psurf(new SurfYZPlane(0., y, z, 0.));
          TrackVector vec(5);
          vec(0) = x;
          vec(1) = 0.;
          vec(2) = px / pz;
          vec(3) = py / pz;
          vec(4) = 1. / p;
          Surface::TrackDirection dir = (pz > 0. ? Surface::FORWARD : Surface::BACKWARD);
          KTrack trk(psurf, vec, dir, pdg);

          // Fill KHitContainer with hits.

          KHitContainerWireX cont;
          cont.fill(trackhits, -1);

          // Count hits in each plane.  Set the preferred plane to be
          // the one with the most hits.

          std::vector<unsigned int> planehits(3, 0);
          for(art::PtrVector<recob::Hit>::const_iterator ih = trackhits.begin();
              ih != trackhits.end(); ++ih) {
            const recob::Hit& hit = **ih;
            unsigned int plane = hit.WireID().Plane;

            if (plane >= planehits.size()) {
              throw cet::exception("TrackKalmanCheater") << "plane " << plane << "...\n"; }
            ++planehits[plane];
          }
          unsigned int prefplane = 0;
          for(unsigned int i=0; i<planehits.size(); ++i) {
            if(planehits[i] > planehits[prefplane])
              prefplane = i;
          }
          if (mf::isDebugEnabled())
            log << "Preferred plane = " << prefplane << "\n";

          // Build and smooth track.

          KGTrack trg(prefplane);
          fKFAlg.setPlane(prefplane);
          bool ok = fKFAlg.buildTrack(trk, trg, fProp, Propagator::FORWARD, cont, false);
          if(ok) {
            ok = fKFAlg.smoothTrackIter(5, trg, fProp);
            if(ok) {
              KETrack tremom;
              bool pok = fKFAlg.fitMomentum(trg, fProp, tremom);
              if(pok)
                fKFAlg.updateMomentum(tremom, fProp, trg);
              ++fNumTrack;
              kalman_tracks.push_back(trg);
            }
          }
          if (mf::isDebugEnabled())
            log << (ok? "Build track succeeded.": "Build track failed.") << "\n";
          
        }
      }
    }
  }

  // Fill histograms.

  if(fHist) {

    // First loop over tracks.

    for(std::deque<KGTrack>::const_iterator k = kalman_tracks.begin();
        k != kalman_tracks.end(); ++k) {
      const KGTrack& trg = *k;

      // Loop over measurements in this track.

      const std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
      for(std::multimap<double, KHitTrack>::const_iterator ih = trackmap.begin();
          ih != trackmap.end(); ++ih) {
        const KHitTrack& trh = (*ih).second;
        const std::shared_ptr<const KHitBase>& hit = trh.getHit();
        double chisq = hit->getChisq();
        fHIncChisq->Fill(chisq);
        const KHit<1>* ph1 = dynamic_cast<const KHit<1>*>(&*hit);
        if(ph1 != 0) {
          double pull = ph1->getResVector()(0) / std::sqrt(ph1->getResError()(0, 0));
          fHPull->Fill(pull);
        }
      }
    }
  }

  // Process Kalman filter tracks into persistent objects.

  tracks->reserve(kalman_tracks.size());
  for(std::deque<KGTrack>::const_iterator k = kalman_tracks.begin();
      k != kalman_tracks.end(); ++k) {
    const KGTrack& kalman_track = *k;

    // Add Track object to collection.

    tracks->push_back(recob::Track());
    kalman_track.fillTrack(tracks->back(), tracks->size() - 1);

    // Make Track to Hit associations.  

    art::PtrVector<recob::Hit> trhits;
    std::vector<unsigned int> hittpindex;
    kalman_track.fillHits(hits, hittpindex);
    util::CreateAssn(*this, evt, *tracks, trhits, *th_assn, tracks->size()-1);

    // Make space points from this track.

    int nspt = spts->size();
    fSpacePointAlg.fillSpacePoints(*spts, kalman_track.TrackMap());

    // Associate newly created space points with hits.
    // Also associate track with newly created space points.

    art::PtrVector<recob::SpacePoint> sptvec;

    // Loop over newly created space points.
      
    for(unsigned int ispt = nspt; ispt < spts->size(); ++ispt) {
      const recob::SpacePoint& spt = (*spts)[ispt];
      art::ProductID sptid = getProductID<std::vector<recob::SpacePoint> >();
      art::Ptr<recob::SpacePoint> sptptr(sptid, ispt, evt.productGetter(sptid));
      sptvec.push_back(sptptr);

      // Make space point to hit associations.

      const art::PtrVector<recob::Hit>& sphits = 
        fSpacePointAlg.getAssociatedHits(spt);
      util::CreateAssn(*this, evt, *spts, sphits, *sph_assn, ispt);
    }

    // Make track to space point associations.

    util::CreateAssn(*this, evt, *tracks, sptvec, *tsp_assn, tracks->size()-1);
  }

  // Add tracks and associations to event.

  evt.put(std::move(tracks));
  evt.put(std::move(spts));
  evt.put(std::move(th_assn));
  evt.put(std::move(tsp_assn));
  evt.put(std::move(sph_assn));
}

//------------------------------------------------------------------------------
void trkf::TrackKalmanCheater::endJob()
{
  mf::LogInfo("TrackKalmanCheater") 
    << "TrackKalmanCheater statistics:\n"
    << "  Number of events = " << fNumEvent << "\n"
    << "  Number of tracks created = " << fNumTrack;
}