trkf::SeedAna Class Reference

Inheritance diagram for trkf::SeedAna:

Classes
struct	MCHists
struct	RecoHists

Public Member Functions
	SeedAna (fhicl::ParameterSet const &pset)
Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDAnalyzer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept
	Analyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	Analyzer (Table< Config > const &config)
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::Observer
	~Observer () noexcept
	Observer (Observer const &)=delete
	Observer (Observer &&)=delete
Observer &	operator= (Observer const &)=delete
Observer &	operator= (Observer &&)=delete
void	registerProducts (ProductDescriptions &, ModuleDescription const &)
void	fillDescriptions (ModuleDescription const &)
fhicl::ParameterSetID	selectorConfig () const
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	analyze (const art::Event &evt) override
void	endJob () override

Private Attributes
std::string	fSeedModuleLabel
std::string	fMCTrackModuleLabel
int	fDump
double	fMinMCKE
double	fMinMCLen
double	fMatchColinearity
double	fMatchDisp
bool	fIgnoreSign
std::map< int, MCHists >	fMCHistMap
std::map< int, RecoHists >	fRecoHistMap
int	fNumEvent

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >
using	ModuleType = EDAnalyzer
Protected Member Functions inherited from art::Observer
std::string const &	processName () const
bool	wantAllEvents () const noexcept
bool	wantEvent (ScheduleID id, Event const &e) const
Handle< TriggerResults >	getTriggerResults (Event const &e) const
	Observer (fhicl::ParameterSet const &config)
	Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
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 230 of file SeedAna_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 548 of file SeedAna_module.cc.

    : EDAnalyzer(pset)
    , fSeedModuleLabel(pset.get<std::string>("SeedModuleLabel"))
    , fMCTrackModuleLabel(pset.get<std::string>("MCTrackModuleLabel"))
    , 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"))
    , fIgnoreSign(pset.get<bool>("IgnoreSign"))
    , fNumEvent(0)
  {
    mf::LogInfo("SeedAna") << "SeedAna configured with the following parameters:\n"
                           << "  SeedModuleLabel = " << fSeedModuleLabel << "\n"
                           << "  MCTrackModuleLabel = " << fMCTrackModuleLabel << "\n"
                           << "  Dump = " << fDump << "\n"
                           << "  MinMCKE = " << fMinMCKE << "\n"
                           << "  MinMCLen = " << fMinMCLen;
  }

Member Function Documentation

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

overrideprivate

Definition at line 574 of file SeedAna_module.cc.

  {
    ++fNumEvent;

    // Optional dump stream.

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

    // Make sure histograms are booked.

    bool mc = !evt.isRealData();

    // Get seed handle.

    art::Handle<std::vector<recob::Seed>> seedh;
    evt.getByLabel(fSeedModuleLabel, seedh);

    // Seed->mc track matching map.

    std::map<const recob::Seed*, int> seedmap;

    if (mc) {

      // Get MCTracks.

      art::Handle<std::vector<sim::MCTrack>> mctrackh;
      evt.getByLabel(fMCTrackModuleLabel, mctrackh);

      // 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 particles.
      auto const detProp =
        art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);

      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

          // 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(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) {

              // 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 << "\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());

              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(mcstartmom.Mag());
              mchists.fHmclen->Fill(plen);

              // Loop over seeds and do matching.

              int nmatch = 0;
              if (seedh.isValid()) {

                // Loop over seeds.

                int nseed = seedh->size();
                for (int i = 0; i < nseed; ++i) {
                  art::Ptr<recob::Seed> pseed(seedh, i);
                  const recob::Seed& seed = *pseed;
                  if (seedmap.count(&seed) == 0) seedmap[&seed] = -1;

                  // Get parameters of this seed.

                  TVector3 pos;
                  TVector3 dir;
                  double err[3];
                  seed.GetPoint(&pos[0], err);
                  seed.GetDirection(&dir[0], err);

                  // Calculate the global-to-local rotation matrix.
                  // Copied from Track.cxx.

                  TMatrixD rot(3, 3);
                  double dirmag = dir.Mag();
                  double diryz = std::sqrt(dir.Y() * dir.Y() + dir.Z() * dir.Z());

                  double sinth = dir.X() / dirmag;
                  double costh = diryz / dirmag;
                  double sinphi = 0.;
                  double cosphi = 1.;
                  if (diryz != 0) {
                    sinphi = -dir.Y() / diryz;
                    cosphi = dir.Z() / diryz;
                  }
                  rot(0, 0) = costh;
                  rot(1, 0) = 0.;
                  rot(2, 0) = sinth;
                  rot(0, 1) = sinth * sinphi;
                  rot(1, 1) = cosphi;
                  rot(2, 1) = -costh * sinphi;
                  rot(0, 2) = -sinth * cosphi;
                  rot(1, 2) = sinphi;
                  rot(2, 2) = costh * cosphi;

                  // Get best matching mc trajectory point.

                  int itraj = mcmatch(detProp, mctrk, seed);
                  if (itraj >= 0) {

                    // Get mc relative position and direction at matching trajectory point.

                    TVector3 mcpos = mctrk[itraj].Position().Vect() - pos;
                    TVector3 mcmom = mctrk[itraj].Momentum().Vect();
                    mcpos[0] += mcdx;

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

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

                    if (mcmoml.Z() < 0.) mcmoml = -mcmoml;
                    double costh = 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::sqrt(u0 * u0 + v0 * v0);

                    // Fill matching histograms.

                    mchists.fHduvcosth->Fill(costh, uv0);
                    if (std::abs(uv0) < fMatchDisp) {

                      // Fill slope matching histograms.

                      mchists.fHmcdudw->Fill(dudw);
                      mchists.fHmcdvdw->Fill(dvdw);
                    }
                    mchists.fHcosth->Fill(costh);
                    if (costh > fMatchColinearity) {

                      // Fill displacement matching histograms.

                      mchists.fHmcu->Fill(u0);
                      mchists.fHmcv->Fill(v0);
                      mchists.fHmcw->Fill(w);

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

                        // Now we have passed all matching cuts and we have a matching
                        // mc particle + seed pair.

                        ++nmatch;
                        seedmap[&seed] = mctrk.TrackID();

                        // Fill matched seed histograms (seed multiplicity).

                        mchists.fHmstartx->Fill(mcstart.X());
                        mchists.fHmstarty->Fill(mcstart.Y());
                        mchists.fHmstartz->Fill(mcstart.Z());
                        mchists.fHmendx->Fill(mcend.X());
                        mchists.fHmendy->Fill(mcend.Y());
                        mchists.fHmendz->Fill(mcend.Z());
                        mchists.fHmtheta->Fill(mcstartmom.Theta());
                        mchists.fHmphi->Fill(mcstartmom.Phi());
                        mchists.fHmtheta_xz->Fill(mctheta_xz);
                        mchists.fHmtheta_yz->Fill(mctheta_yz);
                        mchists.fHmmom->Fill(mcstartmom.Mag());
                        mchists.fHmlen->Fill(plen);
                      }
                    }
                  }
                }

                // If we found at least one matched seed, fill good
                // particle histograms.

                if (nmatch > 0) {
                  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(mcstartmom.Mag());
                  mchists.fHglen->Fill(plen);
                }
              }
            }
          }
        }
      }
    }

    // Loop over seeds and fill reco-only seed histograms.

    if (seedh.isValid()) {

      // Loop over seeds.

      int nseed = seedh->size();

      if (nseed > 0 && pdump != 0) {
        *pdump << "\nReconstructed Seeds\n"
               << "           MCid           x         y         z          dx        dy        dz "
                  "          p\n"
               << "--------------------------------------------------------------------------------"
                  "-----------\n";
      }

      for (int i = 0; i < nseed; ++i) {
        art::Ptr<recob::Seed> pseed(seedh, i);
        const recob::Seed& seed = *pseed;

        // Fill histograms involving reco seeds only.

        TVector3 pos;
        TVector3 dir;
        double err[3];
        seed.GetPoint(&pos[0], err);
        seed.GetDirection(&dir[0], err);
        double mdir = dir.Mag();
        if (mdir != 0.) { dir *= (1. / mdir); }

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

        // Dump seed information here.

        if (pdump) {
          int mcid = seedmap[&seed];
          *pdump << std::setw(15) << mcid << "  " << std::fixed << std::setprecision(2)
                 << std::setw(10) << pos[0] << std::setw(10) << pos[1] << std::setw(10) << pos[2]
                 << "  " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
                 << std::setw(10) << dir[1] << std::setw(10) << dir[2] << "\n";
        }

        // Fill histograms.

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

        rhists.fHx->Fill(pos.X());
        rhists.fHy->Fill(pos.Y());
        rhists.fHz->Fill(pos.Z());
        rhists.fHdist->Fill(dpos);
        rhists.fHtheta->Fill(dir.Theta());
        rhists.fHphi->Fill(dir.Phi());
        rhists.fHtheta_xz->Fill(theta_xz);
        rhists.fHtheta_yz->Fill(theta_yz);
      }
    }
  }

void trkf::SeedAna::endJob ( )

overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 934 of file SeedAna_module.cc.

  {
    // Print summary.

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

    // Fill multiplicity histograms.

    for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
         ++i) {
      const MCHists& mchists = i->second;
      mulcalc(mchists.fHmstartx, mchists.fHmcstartx, mchists.fHmulstartx);
      mulcalc(mchists.fHmstarty, mchists.fHmcstarty, mchists.fHmulstarty);
      mulcalc(mchists.fHmstartz, mchists.fHmcstartz, mchists.fHmulstartz);
      mulcalc(mchists.fHmendx, mchists.fHmcendx, mchists.fHmulendx);
      mulcalc(mchists.fHmendy, mchists.fHmcendy, mchists.fHmulendy);
      mulcalc(mchists.fHmendz, mchists.fHmcendz, mchists.fHmulendz);
      mulcalc(mchists.fHmtheta, mchists.fHmctheta, mchists.fHmultheta);
      mulcalc(mchists.fHmphi, mchists.fHmcphi, mchists.fHmulphi);
      mulcalc(mchists.fHmtheta_xz, mchists.fHmctheta_xz, mchists.fHmultheta_xz);
      mulcalc(mchists.fHmtheta_yz, mchists.fHmctheta_yz, mchists.fHmultheta_yz);
      mulcalc(mchists.fHmmom, mchists.fHmcmom, mchists.fHmulmom);
      mulcalc(mchists.fHmlen, mchists.fHmclen, mchists.fHmullen);
    }
    // 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.fHglen, mchists.fHmclen, mchists.fHelen);
    }
  }

Member Data Documentation

int trkf::SeedAna::fDump

private

Definition at line 352 of file SeedAna_module.cc.

bool trkf::SeedAna::fIgnoreSign

private

Definition at line 357 of file SeedAna_module.cc.

double trkf::SeedAna::fMatchColinearity

private

Definition at line 355 of file SeedAna_module.cc.

double trkf::SeedAna::fMatchDisp

private

Definition at line 356 of file SeedAna_module.cc.

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

private

Definition at line 361 of file SeedAna_module.cc.

std::string trkf::SeedAna::fMCTrackModuleLabel

private

Definition at line 351 of file SeedAna_module.cc.

double trkf::SeedAna::fMinMCKE

private

Definition at line 353 of file SeedAna_module.cc.

double trkf::SeedAna::fMinMCLen

private

Definition at line 354 of file SeedAna_module.cc.

int trkf::SeedAna::fNumEvent

private

Definition at line 366 of file SeedAna_module.cc.

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

private

Definition at line 362 of file SeedAna_module.cc.

std::string trkf::SeedAna::fSeedModuleLabel

private

Definition at line 350 of file SeedAna_module.cc.

---

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

• larreco/larreco/TrackFinder/SeedAna_module.cc