trkf::TrackAna Class Reference

Inheritance diagram for trkf::TrackAna:

Classes
struct	MCHists
struct	RecoHists

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

Private Member Functions
void	analyze (const art::Event &evt) override
void	endJob () override
void	anaStitch (const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp)
std::vector< size_t >	fsort_indexes (const std::vector< double > &v)

Private Attributes
std::string	fTrackModuleLabel
std::string	fMCTrackModuleLabel
std::string	fSpacepointModuleLabel
std::string	fStitchModuleLabel
std::string	fTrkSpptAssocModuleLabel
std::string	fHitSpptAssocModuleLabel
std::string	fHitModuleLabel
int	fDump
double	fMinMCKE
double	fMinMCLen
double	fMatchColinearity
double	fMatchDisp
double	fWMatchDisp
double	fMatchLength
bool	fIgnoreSign
bool	fStitchedAnalysis
std::string	fOrigin
bool	fCheckOrigin
simb::Origin_t	fOriginValue
int	fPrintLevel
std::map< int, MCHists >	fMCHistMap
std::map< int, RecoHists >	fRecoHistMap
int	fNumEvent

Detailed Description

Definition at line 260 of file TrackAna_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 654 of file TrackAna_module.cc.

    : EDAnalyzer(pset)
    , fTrackModuleLabel(pset.get<std::string>("TrackModuleLabel"))
    , fMCTrackModuleLabel(pset.get<std::string>("MCTrackModuleLabel"))
    , fSpacepointModuleLabel(pset.get<std::string>("SpacepointModuleLabel"))
    , fStitchModuleLabel(pset.get<std::string>("StitchModuleLabel"))
    , fTrkSpptAssocModuleLabel(pset.get<std::string>("TrkSpptAssocModuleLabel"))
    , fHitSpptAssocModuleLabel(pset.get<std::string>("HitSpptAssocModuleLabel"))
    , fHitModuleLabel(pset.get<std::string>("HitModuleLabel"))
    , fDump(pset.get<int>("Dump"))
    , fMinMCKE(pset.get<double>("MinMCKE"))
    , fMinMCLen(pset.get<double>("MinMCLen"))
    , fMatchColinearity(pset.get<double>("MatchColinearity"))
    , fMatchDisp(pset.get<double>("MatchDisp"))
    , fWMatchDisp(pset.get<double>("WMatchDisp"))
    , fMatchLength(pset.get<double>("MatchLength"))
    , fIgnoreSign(pset.get<bool>("IgnoreSign"))
    , fStitchedAnalysis(pset.get<bool>("StitchedAnalysis", false))
    , fOrigin(pset.get<std::string>("MCTrackOrigin", "Any"))
    , fPrintLevel(pset.get<int>("PrintLevel", 0))
    , fNumEvent(0)
  {

    // Decide whether to check MCTrack origin
    fCheckOrigin = false;
    fOriginValue = simb::kUnknown;
    if (fOrigin.find("Beam") != std::string::npos) {
      fCheckOrigin = true;
      fOriginValue = simb::kBeamNeutrino;
    }
    else if (fOrigin.find("Cosmic") != std::string::npos) {
      fCheckOrigin = true;
      fOriginValue = simb::kCosmicRay;
    }
    else if (fOrigin.find("Super") != std::string::npos) {
      fCheckOrigin = true;
      fOriginValue = simb::kSuperNovaNeutrino;
    }
    else if (fOrigin.find("Single") != std::string::npos) {
      fCheckOrigin = true;
      fOriginValue = simb::kSingleParticle;
    }

    // Report.

    mf::LogInfo("TrackAna") << "TrackAna configured with the following parameters:\n"
                            << "  TrackModuleLabel = " << fTrackModuleLabel << "\n"
                            << "  MCTrackModuleLabel = " << fMCTrackModuleLabel << "\n"
                            << "  StitchModuleLabel = " << fStitchModuleLabel << "\n"
                            << "  TrkSpptAssocModuleLabel = " << fTrkSpptAssocModuleLabel << "\n"
                            << "  HitSpptAssocModuleLabel = " << fHitSpptAssocModuleLabel << "\n"
                            << "  HitModuleLabel = " << fHitModuleLabel << "\n"
                            << "  Dump = " << fDump << "\n"
                            << "  MinMCKE = " << fMinMCKE << "\n"
                            << "  MinMCLen = " << fMinMCLen << "  Origin = " << fOrigin
                            << " Origin value " << fOriginValue;
  }

Member Function Documentation

void trkf::TrackAna::analyze ( const art::Event &  evt )

overrideprivate

Definition at line 718 of file TrackAna_module.cc.

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

    ++fNumEvent;

    // Optional dump stream.

    std::unique_ptr<mf::LogInfo> pdump;
    if (fDump > 0) {
      --fDump;
      pdump = std::make_unique<mf::LogInfo>("TrackAna");
    }

    // Make sure histograms are booked.

    bool mc = !evt.isRealData();

    // Get MCTracks.

    art::Handle<std::vector<sim::MCTrack>> mctrackh;
    evt.getByLabel(fMCTrackModuleLabel, mctrackh);
    // pair of MCTrack and index of matched reco track
    std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;

    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);

    if (mc && mctrackh.isValid()) {

      selected_mctracks.reserve(mctrackh->size());

      // Dump MCTracks.

      if (pdump) {
        *pdump << "MC Tracks\n"
               << "       Id   pdg           x         y         z          dx        dy        dz "
                  "          p\n"
               << "--------------------------------------------------------------------------------"
                  "-----------\n";
      }

      // Loop over mc tracks, and fill histograms that depend only
      // on mc information.  Also, fill a secondary list of mc tracks
      // that pass various selection criteria.

      for (std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();
           imctrk != mctrackh->end();
           ++imctrk) {
        const sim::MCTrack& mctrk = *imctrk;
        int pdg = mctrk.PdgCode();
        if (fIgnoreSign) pdg = std::abs(pdg);

        // 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

          // check MC track origin?
          if (fCheckOrigin && mctrk.Origin() != fOriginValue) continue;

          // Apply minimum energy cut.

          if (mctrk.Start().E() >= mctrk.Start().Momentum().Mag() + 1000. * fMinMCKE) {

            // Calculate the x offset due to nonzero mc particle time.

            double mctime = mctrk.Start().T();                      // nsec
            double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm

            // Calculate the length of this mc particle inside the fiducial volume.

            TVector3 mcstart;
            TVector3 mcend;
            TVector3 mcstartmom;
            TVector3 mcendmom;
            double plen =
              length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);

            // Apply minimum fiducial length cut.  Always reject particles that have
            // zero length in the tpc regardless of the configured cut.

            if (plen > 0. && plen > fMinMCLen) {

              // This is a good mc particle (capable of making a track).

              selected_mctracks.push_back(std::make_pair(&mctrk, -1));

              // Dump MC particle information here.

              if (pdump) {
                double pstart = mcstartmom.Mag();
                double pend = mcendmom.Mag();
                *pdump << "\nOffset" << std::setw(3) << mctrk.TrackID() << std::setw(6)
                       << mctrk.PdgCode() << "  " << std::fixed << std::setprecision(2)
                       << std::setw(10) << mcdx << "\nStart " << std::setw(3) << mctrk.TrackID()
                       << std::setw(6) << mctrk.PdgCode() << "  " << std::fixed
                       << std::setprecision(2) << std::setw(10) << mcstart[0] << std::setw(10)
                       << mcstart[1] << std::setw(10) << mcstart[2];
                if (pstart > 0.) {
                  *pdump << "  " << std::fixed << std::setprecision(3) << std::setw(10)
                         << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
                         << std::setw(10) << mcstartmom[2] / pstart;
                }
                else
                  *pdump << std::setw(32) << " ";
                *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
                *pdump << "\nEnd   " << std::setw(3) << mctrk.TrackID() << std::setw(6)
                       << mctrk.PdgCode() << "  " << std::fixed << std::setprecision(2)
                       << std::setw(10) << mcend[0] << std::setw(10) << mcend[1] << std::setw(10)
                       << mcend[2];
                if (pend > 0.01) {
                  *pdump << "  " << std::fixed << std::setprecision(3) << std::setw(10)
                         << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend
                         << std::setw(10) << mcendmom[2] / pend;
                }
                else
                  *pdump << std::setw(32) << " ";
                *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
                       << "\nLength: " << plen << "\n";
              }

              // Fill histograms.

              if (fMCHistMap.count(pdg) == 0) {
                std::ostringstream ostr;
                ostr << "MC" << (fIgnoreSign ? "All" : (pdg > 0 ? "Pos" : "Neg")) << std::abs(pdg);
                fMCHistMap.emplace(pdg, MCHists{ostr.str()});
              }
              const MCHists& mchists = fMCHistMap.at(pdg);

              double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
              double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
              double mcmom = mcstartmom.Mag();
              double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
              double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);

              mchists.fHmcstartx->Fill(mcstart.X());
              mchists.fHmcstarty->Fill(mcstart.Y());
              mchists.fHmcstartz->Fill(mcstart.Z());
              mchists.fHmcendx->Fill(mcend.X());
              mchists.fHmcendy->Fill(mcend.Y());
              mchists.fHmcendz->Fill(mcend.Z());
              mchists.fHmctheta->Fill(mcstartmom.Theta());
              mchists.fHmcphi->Fill(mcstartmom.Phi());
              mchists.fHmctheta_xz->Fill(mctheta_xz);
              mchists.fHmctheta_yz->Fill(mctheta_yz);
              mchists.fHmcmom->Fill(mcmom);
              mchists.fHmcmoml->Fill(mcmom);
              mchists.fHmcke->Fill(mcke);
              mchists.fHmckel->Fill(mcke);
              mchists.fHmclen->Fill(plen);
              mchists.fHmclens->Fill(plen);
            }
          }
        }
      }
    } //mc

    // Get tracks and spacepoints and hits
    art::Handle<std::vector<recob::Track>> trackh;
    art::Handle<std::vector<art::PtrVector<recob::Track>>> trackvh;
    art::Handle<std::vector<recob::Hit>> hith;

    evt.getByLabel(fTrackModuleLabel, trackh);
    evt.getByLabel(fStitchModuleLabel, trackvh);
    evt.getByLabel(fHitModuleLabel, hith);

    // Extract all hits into a vector of art::Ptrs (the format used by back tracker).

    std::vector<art::Ptr<recob::Hit>> allhits;
    if (hith.isValid()) {
      allhits.reserve(hith->size());
      for (unsigned int i = 0; i < hith->size(); ++i) {
        allhits.emplace_back(hith, i);
      }
    }

    // Construct FindManyP object to be used for finding track-hit associations.

    art::FindManyP<recob::Hit> tkhit_find(trackh, evt, fTrackModuleLabel);

    // This new top part of TrackAna between two long lines of ************s
    // is particular to analyzing Stitched Tracks.
    // *******************************************************************//

    if (trackvh.isValid() && fStitchedAnalysis) {
      mf::LogDebug("TrackAna") << "TrackAna read " << trackvh->size()
                               << "  vectors of Stitched PtrVectorsof tracks.";
      anaStitch(evt, clockData, detProp);
    }

    if (trackh.isValid()) {

      if (pdump) {
        *pdump << "\nReconstructed Tracks\n"
               << "       Id  MCid           x         y         z          dx        dy        dz "
                  "          p\n"
               << "--------------------------------------------------------------------------------"
                  "-----------\n";
      }

      // Loop over tracks.

      int ntrack = trackh->size();
      for (int i = 0; i < ntrack; ++i) {
        art::Ptr<recob::Track> ptrack(trackh, i);
        const recob::Track& track = *ptrack;

        // Extract hits associated with this track.

        std::vector<art::Ptr<recob::Hit>> trackhits;
        tkhit_find.get(i, trackhits);

        // Calculate the x offset due to nonzero reconstructed time.

        double trackdx = 0;

        // Fill histograms involving reco tracks only.

        int ntraj = track.NumberTrajectoryPoints();
        if (ntraj > 0) {
          TVector3 pos = track.Vertex<TVector3>();
          TVector3 dir = track.VertexDirection<TVector3>();
          TVector3 end = track.End<TVector3>();
          pos[0] += trackdx;
          end[0] += trackdx;

          double dpos = bdist(pos);
          double dend = bdist(end);
          double tlen = length(track);
          double theta_xz = std::atan2(dir.X(), dir.Z());
          double theta_yz = std::atan2(dir.Y(), dir.Z());

          if (fRecoHistMap.count(0) == 0) fRecoHistMap.emplace(0, RecoHists{"Reco"});
          const RecoHists& rhists = fRecoHistMap.at(0);

          rhists.fHstartx->Fill(pos.X());
          rhists.fHstarty->Fill(pos.Y());
          rhists.fHstartz->Fill(pos.Z());
          rhists.fHstartd->Fill(dpos);
          rhists.fHendx->Fill(end.X());
          rhists.fHendy->Fill(end.Y());
          rhists.fHendz->Fill(end.Z());
          rhists.fHendd->Fill(dend);
          rhists.fHtheta->Fill(dir.Theta());
          rhists.fHphi->Fill(dir.Phi());
          rhists.fHtheta_xz->Fill(theta_xz);
          rhists.fHtheta_yz->Fill(theta_yz);

          double mom = 0.;
          if (track.HasMomentum()) mom = track.VertexMomentum();
          rhists.fHmom->Fill(mom);
          rhists.fHmoml->Fill(mom);
          rhists.fHlen->Fill(tlen);
          rhists.fHlens->Fill(tlen);

          // Id of matching mc particle.

          int mcid = -1;

          // Loop over direction.

          for (int swap = 0; swap < 2; ++swap) {

            // Analyze reversed tracks only if start momentum = end momentum.

            if (swap != 0 && track.HasMomentum() &&
                std::abs(track.VertexMomentum() - track.EndMomentum()) > 1.e-3)
              continue;

            // Calculate the global-to-local rotation matrix.

            int start_point = (swap == 0 ? 0 : ntraj - 1);
            TMatrixD rot = track.GlobalToLocalRotationAtPoint<TMatrixD>(start_point);

            // Update track data for reversed case.

            if (swap != 0) {
              rot(1, 0) = -rot(1, 0);
              rot(2, 0) = -rot(2, 0);
              rot(1, 1) = -rot(1, 1);
              rot(2, 1) = -rot(2, 1);
              rot(1, 2) = -rot(1, 2);
              rot(2, 2) = -rot(2, 2);

              pos = track.End<TVector3>();
              dir = -track.EndDirection<TVector3>();
              end = track.Vertex<TVector3>();
              pos[0] += trackdx;
              end[0] += trackdx;

              dpos = bdist(pos);
              dend = bdist(end);
              theta_xz = std::atan2(dir.X(), dir.Z());
              theta_yz = std::atan2(dir.Y(), dir.Z());

              if (track.HasMomentum()) mom = track.EndMomentum();
            }

            // Get covariance matrix.

            const auto& cov = (swap == 0 ? track.VertexCovariance() : track.EndCovariance());

            // Loop over track-like mc particles.

            for (unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
              const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
              int pdg = mctrk.PdgCode();
              if (fIgnoreSign) pdg = std::abs(pdg);
              auto iMCHistMap = fMCHistMap.find(pdg);
              if (iMCHistMap == fMCHistMap.end())
                throw cet::exception("TrackAna") << "no particle with ID=" << pdg << "\n";
              const MCHists& mchists = iMCHistMap->second;

              // Calculate the x offset due to nonzero mc particle time.

              double mctime = mctrk.Start().T();                      // nsec
              double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm

              // Calculate the points where this mc particle enters and leaves the
              // fiducial volume, and the length in the fiducial volume.

              TVector3 mcstart;
              TVector3 mcend;
              TVector3 mcstartmom;
              TVector3 mcendmom;
              double plen =
                length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);

              // Get the displacement of this mc particle in the global coordinate system.

              TVector3 mcpos = mcstart - pos;

              // Rotate the momentum and position to the
              // track-local coordinate system.

              TVector3 mcmoml = rot * mcstartmom;
              TVector3 mcposl = rot * mcpos;

              double colinearity = mcmoml.Z() / mcmoml.Mag();

              double u = mcposl.X();
              double v = mcposl.Y();
              double w = mcposl.Z();

              double pu = mcmoml.X();
              double pv = mcmoml.Y();
              double pw = mcmoml.Z();

              double dudw = pu / pw;
              double dvdw = pv / pw;

              double u0 = u - w * dudw;
              double v0 = v - w * dvdw;
              double uv0 = std::hypot(u0, v0);

              mchists.fHduvcosth->Fill(colinearity, uv0);

              if (std::abs(uv0) < fMatchDisp) {

                // Fill slope matching histograms.

                mchists.fHmcdudw->Fill(dudw);
                mchists.fHmcdvdw->Fill(dvdw);
                mchists.fHdudwpull->Fill(dudw / std::sqrt(cov(2, 2)));
                mchists.fHdvdwpull->Fill(dvdw / std::sqrt(cov(3, 3)));
              }
              mchists.fHcosth->Fill(colinearity);
              if (colinearity > fMatchColinearity) {

                // Fill displacement matching histograms.

                mchists.fHmcu->Fill(u0);
                mchists.fHmcv->Fill(v0);
                mchists.fHmcw->Fill(w);
                mchists.fHupull->Fill(u0 / std::sqrt(cov(0, 0)));
                mchists.fHvpull->Fill(v0 / std::sqrt(cov(1, 1)));

                if (std::abs(uv0) < fMatchDisp) {

                  // Fill matching histograms.

                  double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
                  double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
                  double mcmom = mcstartmom.Mag();
                  double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
                  double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);

                  mchists.fHstartdx->Fill(pos.X() - mcstart.X());
                  mchists.fHstartdy->Fill(pos.Y() - mcstart.Y());
                  mchists.fHstartdz->Fill(pos.Z() - mcstart.Z());
                  mchists.fHenddx->Fill(end.X() - mcend.X());
                  mchists.fHenddy->Fill(end.Y() - mcend.Y());
                  mchists.fHenddz->Fill(end.Z() - mcend.Z());
                  mchists.fHlvsl->Fill(plen, tlen);
                  mchists.fHdl->Fill(tlen - plen);
                  mchists.fHpvsp->Fill(mcmom, mom);
                  double dp = mom - mcmom;
                  mchists.fHdp->Fill(dp);
                  mchists.fHppull->Fill(dp / std::sqrt(cov(4, 4)));
                  if (std::abs(dpos) >= 5. && std::abs(dend) >= 5.) {
                    mchists.fHpvspc->Fill(mcmom, mom);
                    mchists.fHdpc->Fill(dp);
                    mchists.fHppullc->Fill(dp / std::sqrt(cov(4, 4)));
                  }

                  // Count this track as well-reconstructed if it is matched to an
                  // mc particle (which it is if get here), and if in addition the
                  // starting position (w) matches and the reconstructed track length
                  // is more than 0.5 of the mc particle trajectory length.

                  bool good = std::abs(w) <= fWMatchDisp && tlen > fMatchLength * plen;
                  if (good) {
                    mcid = mctrk.TrackID();

                    // Calculate and fill hit efficiency and purity.

                    std::set<int> tkidset;
                    tkidset.insert(mcid);
                    double hiteff = bt_serv->HitCollectionEfficiency(
                      clockData, tkidset, trackhits, allhits, geo::k3D);
                    double hitpurity = bt_serv->HitCollectionPurity(clockData, tkidset, trackhits);
                    mchists.fHHitEff->Fill(hiteff);
                    mchists.fHHitPurity->Fill(hitpurity);

                    // Fill efficiency numerator histograms.

                    mchists.fHgstartx->Fill(mcstart.X());
                    mchists.fHgstarty->Fill(mcstart.Y());
                    mchists.fHgstartz->Fill(mcstart.Z());
                    mchists.fHgendx->Fill(mcend.X());
                    mchists.fHgendy->Fill(mcend.Y());
                    mchists.fHgendz->Fill(mcend.Z());
                    mchists.fHgtheta->Fill(mcstartmom.Theta());
                    mchists.fHgphi->Fill(mcstartmom.Phi());
                    mchists.fHgtheta_xz->Fill(mctheta_xz);
                    mchists.fHgtheta_yz->Fill(mctheta_yz);
                    mchists.fHgmom->Fill(mcmom);
                    mchists.fHgmoml->Fill(mcmom);
                    mchists.fHgke->Fill(mcke);
                    mchists.fHgkel->Fill(mcke);
                    mchists.fHglen->Fill(plen);
                    mchists.fHglens->Fill(plen);

                    // set the match flag
                    selected_mctracks[imc].second = i;

                    if (fPrintLevel > 0) {
                      const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mcid);
                      float KE = ptkl->E() - ptkl->Mass();
                      std::string KEUnits = " GeV";
                      if (mctrk.Origin() != simb::kCosmicRay) {
                        // MeV for low energy particles
                        KE *= 1000;
                        KEUnits = " MeV";
                      }
                      mf::LogVerbatim("TrackAna")
                        << evt.run() << "." << evt.event() << " Match MCTkID " << std::setw(6)
                        << mctrk.TrackID() << " Origin " << mctrk.Origin() << " PDG" << std::setw(5)
                        << mctrk.PdgCode() << " KE" << std::setw(4) << (int)KE << KEUnits
                        << " RecoTrkID " << track.ID() << " hitEff " << std::setprecision(2)
                        << hiteff << " hitPur " << hitpurity;
                      int sWire, sTick, eWire, eTick;
                      // this won't work for DUNE
                      for (unsigned short ipl = 0; ipl < geom->Nplanes(); ++ipl) {
                        sWire = geom->NearestWire(mcstart, ipl, 0, 0);
                        sTick = detProp.ConvertXToTicks(mcstart[0], ipl, 0, 0);
                        eWire = geom->NearestWire(mcend, ipl, 0, 0);
                        eTick = detProp.ConvertXToTicks(mcend[0], ipl, 0, 0);
                        mf::LogVerbatim("TrackAna")
                          << "   Wire:Tick in Pln " << ipl << " W:T " << sWire << ":" << sTick
                          << " - " << eWire << ":" << eTick;
                      } // ipl
                    }   // fPrintLevel == 2
                  }     // good
                }
              }
            }
          }

          // Dump track information here.

          if (pdump) {
            TVector3 pos = track.Vertex<TVector3>();
            TVector3 dir = track.VertexDirection<TVector3>();
            TVector3 end = track.End<TVector3>();
            pos[0] += trackdx;
            end[0] += trackdx;
            TVector3 enddir = track.EndDirection<TVector3>();
            double pstart = track.VertexMomentum();
            double pend = track.EndMomentum();
            *pdump << "\nOffset" << std::setw(3) << track.ID() << std::setw(6) << mcid << "  "
                   << std::fixed << std::setprecision(2) << std::setw(10) << trackdx << "\nStart "
                   << std::setw(3) << track.ID() << std::setw(6) << mcid << "  " << std::fixed
                   << std::setprecision(2) << std::setw(10) << pos[0] << std::setw(10) << pos[1]
                   << std::setw(10) << pos[2];
            if (pstart > 0.) {
              *pdump << "  " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
                     << std::setw(10) << dir[1] << std::setw(10) << dir[2];
            }
            else
              *pdump << std::setw(32) << " ";
            *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
            *pdump << "\nEnd   " << std::setw(3) << track.ID() << std::setw(6) << mcid << "  "
                   << std::fixed << std::setprecision(2) << std::setw(10) << end[0] << std::setw(10)
                   << end[1] << std::setw(10) << end[2];
            if (pend > 0.01) {
              *pdump << "  " << std::fixed << std::setprecision(3) << std::setw(10) << enddir[0]
                     << std::setw(10) << enddir[1] << std::setw(10) << enddir[2];
            }
            else
              *pdump << std::setw(32) << " ";
            *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
                   << "\nLength: " << tlen << "\n";
          }
        }
      }
    } // i

    // print out un-matched MC tracks
    if (fPrintLevel > 0) {
      for (unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
        if (selected_mctracks[imc].second >= 0) continue;
        const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
        const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mctrk.TrackID());
        float KE = ptkl->E() - ptkl->Mass();
        std::string KEUnits = " GeV";
        if (mctrk.Origin() != simb::kCosmicRay) {
          // MeV for low energy particles
          KE *= 1000;
          KEUnits = " MeV";
        }
        // find the start/end wire:time in each plane
        TVector3 mcstart, mcend, mcstartmom, mcendmom;
        double mcdx = mctrk.Start().T() * 1.e-3 * detProp.DriftVelocity(); // cm
        double plen = length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
        mf::LogVerbatim("TrackAna")
          << evt.run() << "." << evt.event() << " NoMat MCTkID " << std::setw(6) << mctrk.TrackID()
          << " Origin " << mctrk.Origin() << " PDG" << std::setw(5) << mctrk.PdgCode() << " KE"
          << std::setw(4) << (int)KE << KEUnits << " Length " << std::fixed << std::setprecision(1)
          << plen << " cm";
        if (fPrintLevel > 1) {
          int sWire, sTick, eWire, eTick;
          // this won't work for DUNE
          for (unsigned short ipl = 0; ipl < geom->Nplanes(); ++ipl) {
            sWire = geom->NearestWire(mcstart, ipl, 0, 0);
            sTick = detProp.ConvertXToTicks(mcstart[0], ipl, 0, 0);
            eWire = geom->NearestWire(mcend, ipl, 0, 0);
            eTick = detProp.ConvertXToTicks(mcend[0], ipl, 0, 0);
            mf::LogVerbatim("TrackAna") << "   Wire:Tick in Pln " << ipl << " W:T " << sWire << ":"
                                        << sTick << " - " << eWire << ":" << eTick;
          } // ipl
        }   // fPrintLevel > 1
      }     // imc
    }       // fPrintLevel > 0
  }

void trkf::TrackAna::anaStitch ( const art::Event &  evt, detinfo::DetectorClocksData const &  clockData, detinfo::DetectorPropertiesData const &  detProp  )

private

Definition at line 1283 of file TrackAna_module.cc.

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

    std::map<int, std::map<int, art::PtrVector<recob::Hit>>> hitmap; // trkID, otrk, hitvec
    std::map<int, int> KEmap; // length traveled in det [cm]?, trkID want to sort by KE
    bool mc = !evt.isRealData();
    art::Handle<std::vector<recob::Track>> trackh;
    art::Handle<std::vector<recob::SpacePoint>> sppth;
    art::Handle<std::vector<art::PtrVector<recob::Track>>> trackvh;

    evt.getByLabel(fTrackModuleLabel, trackh);
    evt.getByLabel(fSpacepointModuleLabel, sppth);
    evt.getByLabel(fStitchModuleLabel, trackvh);
    int ntv(trackvh->size());

    std::vector<art::PtrVector<recob::Track>>::const_iterator cti = trackvh->begin();

    if (trackh.isValid()) {
      art::FindManyP<recob::SpacePoint> fswhole(trackh, evt, fTrkSpptAssocModuleLabel);
      int nsppts_assnwhole = fswhole.size();
      std::cout << "TrackAna: Number of clumps of Spacepoints from Assn for all Tracks: "
                << nsppts_assnwhole << std::endl;
    }

    if (fRecoHistMap.count(0) == 0) {
      fRecoHistMap.emplace(0, RecoHists{"Reco"});
      std::cout << "\n"
                << "\t\t  TrkAna: Fresh fRecoHistMap[0] ******* \n"
                << std::endl;
    }
    const RecoHists& rhistsStitched = fRecoHistMap.at(0);

    std::vector<std::vector<unsigned int>> NtrkIdsAll;
    std::vector<double> ntvsorted;
    hitmap.clear();
    KEmap.clear();

    // Look at the components of the stitched tracks. Grab their sppts/hits from Assns.
    for (int o = 0; o < ntv; ++o) // o for outer
    {
      const art::PtrVector<recob::Track> pvtrack(*(cti++));
      int ntrack = pvtrack.size();
      std::vector<std::vector<unsigned int>> NtrkId_Hit; // hit IDs in inner tracks
      std::vector<unsigned int> vecMode;
      art::FindManyP<recob::SpacePoint> fs(pvtrack, evt, fTrkSpptAssocModuleLabel);

      for (int i = 0; i < ntrack; ++i) {

        // From gdb> ptype fs, the vector of Ptr<SpacePoint>s it appears is
        // grabbed after fs.at(0)
        bool assns(true);
        try {
          // Got Spacepoints from this Track; now get Hits from those
          // Spacepoints.

          int nsppts_assn = fs.at(i).size();

          const auto& sppt = fs.at(i); //.at(0);
          // since we're in a try and worried about failure, we won't pull the following
          // FindManyP out of the loop.
          art::FindManyP<recob::Hit> fh(sppt, evt, fHitSpptAssocModuleLabel);

          // Importantly, loop on all sppts, though they don't all contribute to the track.
          // As opposed to looping on the trajectory pts, which is a lower number.
          // Also, important, in job in whch this runs I set TrackKal3DSPS parameter MaxPass=1,
          // cuz I don't want merely the sparse set of sppts as follows from the uncontained
          // MS-measurement in 2nd pass.
          std::vector<unsigned int> vecNtrkIds;
          for (int is = 0; is < nsppts_assn; ++is) {
            int nhits = fh.at(is).size(); // should be 2 or 3: number of planes.
            for (int ih = 0; ih < nhits; ++ih) {
              const auto& hit = fh.at(is).at(ih); // Our vector is after the .at(is) this time.
              if (hit->SignalType() != geo::kCollection) continue;
              rhistsStitched.fHHitChg->Fill(hit->Integral());
              rhistsStitched.fHHitWidth->Fill(2. * hit->RMS());
              if (mc) {
                std::vector<sim::TrackIDE> tids = bt_serv->HitToTrackIDEs(clockData, hit);
                // more here.
                // Loop over track ids.
                bool justOne(true); // Only take first trk that contributed to this hit
                for (std::vector<sim::TrackIDE>::const_iterator itid = tids.begin();
                     itid != tids.end();
                     ++itid) {
                  int trackID = std::abs(itid->trackID);
                  hitmap[trackID][o].push_back(hit);

                  if (justOne) {
                    vecNtrkIds.push_back(trackID);
                    justOne = false;
                  }
                  // Add hit to PtrVector corresponding to this track id.
                  rhistsStitched.fHHitTrkId->Fill(trackID);
                  const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(trackID);
                  if (!part) break;

                  rhistsStitched.fHHitPdg->Fill(part->PdgCode());
                  // This really needs to be indexed as KE deposited in volTPC, not just KE. EC, 24-July-2014.

                  TVector3 mcstart;
                  TVector3 mcend;
                  TVector3 mcstartmom;
                  TVector3 mcendmom;
                  double mctime = part->T();                              // nsec
                  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm

                  double plen =
                    length(clockData, detProp, *part, mcdx, mcstart, mcend, mcstartmom, mcendmom);

                  KEmap[(int)(1e6 * plen)] = trackID; // multiple assignment but always the same, so
                                                      // fine.
                }

              } // mc
            }   //  hits

          } //    spacepoints

          if (mc) {
            NtrkId_Hit.push_back(vecNtrkIds);
            // Find the trkID mode for this i^th track
            unsigned int ii(1);
            int max(-12), n(1), ind(0);
            std::sort(vecNtrkIds.begin(), vecNtrkIds.end());
            std::vector<unsigned int> strkIds(vecNtrkIds);
            while (ii < vecNtrkIds.size()) {
              if (strkIds.at(ii) != strkIds.at(ii - 1)) { n = 1; }
              else {
                n++;
              }
              if (n > max) {
                max = n;
                ind = ii;
              }
              ii++;
            }
            unsigned int mode(sim::NoParticleId);
            if (strkIds.begin() != strkIds.end()) mode = strkIds.at(ind);
            vecMode.push_back(mode);

            if (strkIds.size() != 0)
              rhistsStitched.fModeFrac->Fill((double)max / (double)strkIds.size());
            else
              rhistsStitched.fModeFrac->Fill(-1.0);
          } // mc

        } // end try
        catch (cet::exception& x) {
          assns = false;
        }
        if (!assns) throw cet::exception("TrackAna") << "Bad Associations. \n";

      } // i

      if (mc) {
        // one vector per o trk, for all modes of stitched i trks
        NtrkIdsAll.push_back(vecMode);

        std::unique(NtrkIdsAll.back().begin(), NtrkIdsAll.back().end());
        double sum(0.0);
        for (auto const val : NtrkIdsAll.back()) {
          sum += hitmap[val][o].size();
        }
        ntvsorted.push_back(sum);
      }

      //
    } // o

    int vtmp(0);
    // get KEmap indices by most energetic first, least last.
    for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
      vtmp++;
    }

    // get o trk indices by their hits. Most populous track first, least last.
    for (auto const o : fsort_indexes(ntvsorted)) {
      int v(0);
      // get KEmap indices by longest trajectory first, least last.
      for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
        int val = it->second; // grab trkIDs in order, since they're sorted by KE
        rhistsStitched.fNTrkIdTrks3->Fill(o, v, hitmap[val][o].size());
        v++;
      }
    }

    // In how many o tracks did each trkId appear? Histo it. Would like it to be
    // precisely 1. Histo it vs. particle KE.
    flattener flat(NtrkIdsAll);
    std::vector<unsigned int>& v = flat;
    //  auto const it ( std::unique(v.begin(),v.end()) ); // never use this it,
    //  perhaps.
    for (auto const val : v) {
      if (val != (unsigned int)sim::NoParticleId) {
        const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(val);
        double T(part->E() - 0.001 * part->Mass());
        rhistsStitched.fNTrkIdTrks->Fill(std::count(v.begin(), v.end(), val));
        rhistsStitched.fNTrkIdTrks2->Fill(std::count(v.begin(), v.end(), val), T);
      }
      else {
        rhistsStitched.fNTrkIdTrks2->Fill(-1.0, 0.0);
      }
    }
  }

void trkf::TrackAna::endJob ( )

overrideprivate

Definition at line 1493 of file TrackAna_module.cc.

  {
    // Print summary.

    mf::LogInfo("TrackAna") << "TrackAna statistics:\n"
                            << "  Number of events = " << fNumEvent;

    // Fill efficiency histograms.

    for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
         ++i) {
      const MCHists& mchists = i->second;
      effcalc(mchists.fHgstartx, mchists.fHmcstartx, mchists.fHestartx);
      effcalc(mchists.fHgstarty, mchists.fHmcstarty, mchists.fHestarty);
      effcalc(mchists.fHgstartz, mchists.fHmcstartz, mchists.fHestartz);
      effcalc(mchists.fHgendx, mchists.fHmcendx, mchists.fHeendx);
      effcalc(mchists.fHgendy, mchists.fHmcendy, mchists.fHeendy);
      effcalc(mchists.fHgendz, mchists.fHmcendz, mchists.fHeendz);
      effcalc(mchists.fHgtheta, mchists.fHmctheta, mchists.fHetheta);
      effcalc(mchists.fHgphi, mchists.fHmcphi, mchists.fHephi);
      effcalc(mchists.fHgtheta_xz, mchists.fHmctheta_xz, mchists.fHetheta_xz);
      effcalc(mchists.fHgtheta_yz, mchists.fHmctheta_yz, mchists.fHetheta_yz);
      effcalc(mchists.fHgmom, mchists.fHmcmom, mchists.fHemom);
      effcalc(mchists.fHgmoml, mchists.fHmcmoml, mchists.fHemoml);
      effcalc(mchists.fHgke, mchists.fHmcke, mchists.fHeke);
      effcalc(mchists.fHgkel, mchists.fHmckel, mchists.fHekel);
      effcalc(mchists.fHglen, mchists.fHmclen, mchists.fHelen);
      effcalc(mchists.fHglens, mchists.fHmclens, mchists.fHelens);
    }
  }

std::vector< size_t > trkf::TrackAna::fsort_indexes ( const std::vector< double > &  v )

private

Definition at line 1526 of file TrackAna_module.cc.

  {
    // initialize original index locations
    std::vector<size_t> idx(v.size());
    for (size_t i = 0; i != idx.size(); ++i)
      idx[i] = i;
    // sort indexes based on comparing values in v
    std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) {
      return v[i1] > v[i2];
    }); // Want most occupied trks first. was <, EC, 21-July.
    return idx;
  }

Member Data Documentation

bool trkf::TrackAna::fCheckOrigin

private

Definition at line 431 of file TrackAna_module.cc.

int trkf::TrackAna::fDump

private

Definition at line 420 of file TrackAna_module.cc.

std::string trkf::TrackAna::fHitModuleLabel

private

Definition at line 418 of file TrackAna_module.cc.

std::string trkf::TrackAna::fHitSpptAssocModuleLabel

private

Definition at line 417 of file TrackAna_module.cc.

bool trkf::TrackAna::fIgnoreSign

private

Definition at line 427 of file TrackAna_module.cc.

double trkf::TrackAna::fMatchColinearity

private

Definition at line 423 of file TrackAna_module.cc.

double trkf::TrackAna::fMatchDisp

private

Definition at line 424 of file TrackAna_module.cc.

double trkf::TrackAna::fMatchLength

private

Definition at line 426 of file TrackAna_module.cc.

std::map<int, MCHists> trkf::TrackAna::fMCHistMap

private

Definition at line 437 of file TrackAna_module.cc.

std::string trkf::TrackAna::fMCTrackModuleLabel

private

Definition at line 413 of file TrackAna_module.cc.

double trkf::TrackAna::fMinMCKE

private

Definition at line 421 of file TrackAna_module.cc.

double trkf::TrackAna::fMinMCLen

private

Definition at line 422 of file TrackAna_module.cc.

int trkf::TrackAna::fNumEvent

private

Definition at line 442 of file TrackAna_module.cc.

std::string trkf::TrackAna::fOrigin

private

Definition at line 430 of file TrackAna_module.cc.

simb::Origin_t trkf::TrackAna::fOriginValue

private

Definition at line 432 of file TrackAna_module.cc.

int trkf::TrackAna::fPrintLevel

private

Definition at line 433 of file TrackAna_module.cc.

std::map<int, RecoHists> trkf::TrackAna::fRecoHistMap

private

Definition at line 438 of file TrackAna_module.cc.

std::string trkf::TrackAna::fSpacepointModuleLabel

private

Definition at line 414 of file TrackAna_module.cc.

bool trkf::TrackAna::fStitchedAnalysis

private

Definition at line 428 of file TrackAna_module.cc.

std::string trkf::TrackAna::fStitchModuleLabel

private

Definition at line 415 of file TrackAna_module.cc.

std::string trkf::TrackAna::fTrackModuleLabel

private

Definition at line 412 of file TrackAna_module.cc.

std::string trkf::TrackAna::fTrkSpptAssocModuleLabel

private

Definition at line 416 of file TrackAna_module.cc.

double trkf::TrackAna::fWMatchDisp

private

Definition at line 425 of file TrackAna_module.cc.

---

The documentation for this class was generated from the following file:

• TrackAna_module.cc