Inheritance diagram for trkf::TrackKalmanCheater:

Public Member Functions
	TrackKalmanCheater (fhicl::ParameterSet const &pset)

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	produce (art::Event &e) override

void	beginJob () override
	Begin job method. More...

void	endJob () override
	End job method. More...

Private Attributes
bool	fHist
	Make histograms. More...

KalmanFilterAlg	fKFAlg
	Kalman filter algorithm. More...

SpacePointAlg	fSpacePointAlg
	Space point algorithm. More...

bool	fUseClusterHits
	Use cluster hits or all hits? More...

std::string	fHitModuleLabel
	Unclustered Hits. More...

std::string	fClusterModuleLabel
	Clustered Hits. More...

double	fMaxTcut
	Maximum delta ray energy in MeV for restricted dE/dx. More...

TH1F *	fHIncChisq {nullptr}
	Incremental chisquare. More...

TH1F *	fHPull {nullptr}
	Hit pull. More...

int	fNumEvent {0}
	Number of events seen. More...

int	fNumTrack {0}
	Number of tracks produced. More...

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 63 of file TrackKalmanCheater_module.cc.

Constructor & Destructor Documentation

trkf::TrackKalmanCheater::TrackKalmanCheater ( fhicl::ParameterSet const & pset )

explicit

Constructor.

Arguments:

p - Fcl parameters.

Definition at line 104 of file TrackKalmanCheater_module.cc.

   : EDProducer{pset}
   , fHist(pset.get<bool>("Hist"))
   , fKFAlg(pset.get<fhicl::ParameterSet>("KalmanFilterAlg"))
   , fSpacePointAlg(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"))
 {
   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;
 }

Member Function Documentation

void trkf::TrackKalmanCheater::beginJob ( )

overrideprivatevirtual

Begin job method.

Reimplemented from art::EDProducer.

Definition at line 132 of file TrackKalmanCheater_module.cc.

 {
   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::endJob ( )

overrideprivatevirtual

End job method.

Reimplemented from art::EDProducer.

Definition at line 450 of file TrackKalmanCheater_module.cc.

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

void trkf::TrackKalmanCheater::produce ( art::Event & evt )

overrideprivatevirtual

Produce method.

Arguments:

e - Art event.

This method extracts Hit from the event and produces and adds Track objects.

Implements art::EDProducer.

Definition at line 157 of file TrackKalmanCheater_module.cc.

 {
   ++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 const> geom;
   art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
   art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
   // 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(clockData, *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();
 
   auto const det_prop =
     art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
 
   PropYZPlane const prop{det_prop, fMaxTcut, true};
 
   // Loop over geant particles.
 
   {
     // use mf::LogDebug instead of MF_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 (not accepted_particle(apdg)) { continue; }
 
       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::hypot(px, py, 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(det_prop, 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);
         if (fKFAlg.buildTrack(trk, trg, prop, Propagator::FORWARD, cont, false)) {
           if (fKFAlg.smoothTrackIter(5, trg, prop)) {
             KETrack tremom;
             if (fKFAlg.fitMomentum(trg, prop, tremom)) { fKFAlg.updateMomentum(tremom, prop, trg); }
             ++fNumTrack;
             kalman_tracks.push_back(trg);
           }
         }
       }
     }
   }
 
   // Fill histograms.
 
   if (fHist) {
 
     // First loop over tracks.
 
     for (auto const& trg : kalman_tracks) {
 
       // 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 (auto const& kalman_track : kalman_tracks) {
     tracks->push_back(recob::Track());
     kalman_track.fillTrack(det_prop, 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(evt, *tracks, trhits, *th_assn, tracks->size() - 1);
 
     // Make space points from this track.
 
     int nspt = spts->size();
     fSpacePointAlg.fillSpacePoints(det_prop, *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 = evt.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(evt, *spts, sphits, *sph_assn, ispt);
     }
 
     // Make track to space point associations.
 
     util::CreateAssn(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));
 }