Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::TrackAna Class Reference
Inheritance diagram for trkf::TrackAna:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Classes

struct  MCHists
 
struct  RecoHists
 

Public Member Functions

 TrackAna (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
 
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
 

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

660  : EDAnalyzer(pset)
661  , fTrackModuleLabel(pset.get<std::string>("TrackModuleLabel"))
662  , fMCTrackModuleLabel(pset.get<std::string>("MCTrackModuleLabel"))
663  , fSpacepointModuleLabel(pset.get<std::string>("SpacepointModuleLabel"))
664  , fStitchModuleLabel(pset.get<std::string>("StitchModuleLabel"))
665  , fTrkSpptAssocModuleLabel(pset.get<std::string>("TrkSpptAssocModuleLabel"))
666  , fHitSpptAssocModuleLabel(pset.get<std::string>("HitSpptAssocModuleLabel"))
667  , fHitModuleLabel(pset.get<std::string>("HitModuleLabel"))
668  , fDump(pset.get<int>("Dump"))
669  , fMinMCKE(pset.get<double>("MinMCKE"))
670  , fMinMCLen(pset.get<double>("MinMCLen"))
671  , fMatchColinearity(pset.get<double>("MatchColinearity"))
672  , fMatchDisp(pset.get<double>("MatchDisp"))
673  , fWMatchDisp(pset.get<double>("WMatchDisp"))
674  , fMatchLength(pset.get<double>("MatchLength"))
675  , fIgnoreSign(pset.get<bool>("IgnoreSign"))
676  , fStitchedAnalysis(pset.get<bool>("StitchedAnalysis", false))
677  , fOrigin(pset.get<std::string>("MCTrackOrigin", "Any"))
678  , fPrintLevel(pset.get<int>("PrintLevel", 0))
679  , fNumEvent(0)
680  {
681 
682  // Decide whether to check MCTrack origin
683  fCheckOrigin = false;
684  fOriginValue = simb::kUnknown;
685  if (fOrigin.find("Beam") != std::string::npos) {
686  fCheckOrigin = true;
687  fOriginValue = simb::kBeamNeutrino;
688  }
689  else if (fOrigin.find("Cosmic") != std::string::npos) {
690  fCheckOrigin = true;
691  fOriginValue = simb::kCosmicRay;
692  }
693  else if (fOrigin.find("Super") != std::string::npos) {
694  fCheckOrigin = true;
695  fOriginValue = simb::kSuperNovaNeutrino;
696  }
697  else if (fOrigin.find("Single") != std::string::npos) {
698  fCheckOrigin = true;
699  fOriginValue = simb::kSingleParticle;
700  }
701 
702  // Report.
703 
704  mf::LogInfo("TrackAna") << "TrackAna configured with the following parameters:\n"
705  << " TrackModuleLabel = " << fTrackModuleLabel << "\n"
706  << " MCTrackModuleLabel = " << fMCTrackModuleLabel << "\n"
707  << " StitchModuleLabel = " << fStitchModuleLabel << "\n"
708  << " TrkSpptAssocModuleLabel = " << fTrkSpptAssocModuleLabel << "\n"
709  << " HitSpptAssocModuleLabel = " << fHitSpptAssocModuleLabel << "\n"
710  << " HitModuleLabel = " << fHitModuleLabel << "\n"
711  << " Dump = " << fDump << "\n"
712  << " MinMCKE = " << fMinMCKE << "\n"
713  << " MinMCLen = " << fMinMCLen << " Origin = " << fOrigin
714  << " Origin value " << fOriginValue;
715  }
trkf::TrackAna::fHitSpptAssocModuleLabel
std::string fHitSpptAssocModuleLabel
Definition: TrackAna_module.cc:417
trkf::TrackAna::fTrkSpptAssocModuleLabel
std::string fTrkSpptAssocModuleLabel
Definition: TrackAna_module.cc:416
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
trkf::TrackAna::fSpacepointModuleLabel
std::string fSpacepointModuleLabel
Definition: TrackAna_module.cc:414
trkf::TrackAna::fStitchModuleLabel
std::string fStitchModuleLabel
Definition: TrackAna_module.cc:415
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
simb::kSingleParticle
single particles thrown at the detector
Definition: MCTruth.h:26
trkf::TrackAna::fMCTrackModuleLabel
std::string fMCTrackModuleLabel
Definition: TrackAna_module.cc:413
trkf::TrackAna::fOriginValue
simb::Origin_t fOriginValue
Definition: TrackAna_module.cc:432
trkf::TrackAna::fOrigin
std::string fOrigin
Definition: TrackAna_module.cc:430
simb::kSuperNovaNeutrino
Supernova neutrinos.
Definition: MCTruth.h:25
trkf::TrackAna::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: TrackAna_module.cc:412
trkf::TrackAna::fHitModuleLabel
std::string fHitModuleLabel
Definition: TrackAna_module.cc:418
simb::kCosmicRay
Cosmic rays.
Definition: MCTruth.h:24
simb::kBeamNeutrino
Beam neutrinos.
Definition: MCTruth.h:23

Member Function Documentation

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

Definition at line 718 of file TrackAna_module.cc.

719  {
720  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
721  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
722  art::ServiceHandle<geo::Geometry const> geom;
723 
724  ++fNumEvent;
725 
726  // Optional dump stream.
727 
728  std::unique_ptr<mf::LogInfo> pdump;
729  if (fDump > 0) {
730  --fDump;
731  pdump = std::make_unique<mf::LogInfo>("TrackAna");
732  }
733 
734  // Make sure histograms are booked.
735 
736  bool mc = !evt.isRealData();
737 
738  // Get MCTracks.
739 
740  art::Handle<std::vector<sim::MCTrack>> mctrackh;
741  evt.getByLabel(fMCTrackModuleLabel, mctrackh);
742  // pair of MCTrack and index of matched reco track
743  std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;
744 
745  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
746  auto const detProp =
747  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
748 
749  if (mc && mctrackh.isValid()) {
750 
751  selected_mctracks.reserve(mctrackh->size());
752 
753  // Dump MCTracks.
754 
755  if (pdump) {
756  *pdump << "MC Tracks\n"
757  << " Id pdg x y z dx dy dz "
758  " p\n"
759  << "--------------------------------------------------------------------------------"
760  "-----------\n";
761  }
762 
763  // Loop over mc tracks, and fill histograms that depend only
764  // on mc information. Also, fill a secondary list of mc tracks
765  // that pass various selection criteria.
766 
767  for (std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();
768  imctrk != mctrackh->end();
769  ++imctrk) {
770  const sim::MCTrack& mctrk = *imctrk;
771  int pdg = mctrk.PdgCode();
772  if (fIgnoreSign) pdg = std::abs(pdg);
773 
774  // Ignore everything except stable charged nonshowering particles.
775 
776  int apdg = std::abs(pdg);
777  if (apdg == 13 || // Muon
778  apdg == 211 || // Charged pion
779  apdg == 321 || // Charged kaon
780  apdg == 2212) { // (Anti)proton
781 
782  // check MC track origin?
783  if (fCheckOrigin && mctrk.Origin() != fOriginValue) continue;
784 
785  // Apply minimum energy cut.
786 
787  if (mctrk.Start().E() >= mctrk.Start().Momentum().Mag() + 1000. * fMinMCKE) {
788 
789  // Calculate the x offset due to nonzero mc particle time.
790 
791  double mctime = mctrk.Start().T(); // nsec
792  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
793 
794  // Calculate the length of this mc particle inside the fiducial volume.
795 
796  TVector3 mcstart;
797  TVector3 mcend;
798  TVector3 mcstartmom;
799  TVector3 mcendmom;
800  double plen =
801  length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
802 
803  // Apply minimum fiducial length cut. Always reject particles that have
804  // zero length in the tpc regardless of the configured cut.
805 
806  if (plen > 0. && plen > fMinMCLen) {
807 
808  // This is a good mc particle (capable of making a track).
809 
810  selected_mctracks.push_back(std::make_pair(&mctrk, -1));
811 
812  // Dump MC particle information here.
813 
814  if (pdump) {
815  double pstart = mcstartmom.Mag();
816  double pend = mcendmom.Mag();
817  *pdump << "\nOffset" << std::setw(3) << mctrk.TrackID() << std::setw(6)
818  << mctrk.PdgCode() << " " << std::fixed << std::setprecision(2)
819  << std::setw(10) << mcdx << "\nStart " << std::setw(3) << mctrk.TrackID()
820  << std::setw(6) << mctrk.PdgCode() << " " << std::fixed
821  << std::setprecision(2) << std::setw(10) << mcstart[0] << std::setw(10)
822  << mcstart[1] << std::setw(10) << mcstart[2];
823  if (pstart > 0.) {
824  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10)
825  << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
826  << std::setw(10) << mcstartmom[2] / pstart;
827  }
828  else
829  *pdump << std::setw(32) << " ";
830  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
831  *pdump << "\nEnd " << std::setw(3) << mctrk.TrackID() << std::setw(6)
832  << mctrk.PdgCode() << " " << std::fixed << std::setprecision(2)
833  << std::setw(10) << mcend[0] << std::setw(10) << mcend[1] << std::setw(10)
834  << mcend[2];
835  if (pend > 0.01) {
836  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10)
837  << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend
838  << std::setw(10) << mcendmom[2] / pend;
839  }
840  else
841  *pdump << std::setw(32) << " ";
842  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
843  << "\nLength: " << plen << "\n";
844  }
845 
846  // Fill histograms.
847 
848  if (fMCHistMap.count(pdg) == 0) {
849  std::ostringstream ostr;
850  ostr << "MC" << (fIgnoreSign ? "All" : (pdg > 0 ? "Pos" : "Neg")) << std::abs(pdg);
851  fMCHistMap.emplace(pdg, MCHists{ostr.str()});
852  }
853  const MCHists& mchists = fMCHistMap.at(pdg);
854 
855  double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
856  double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
857  double mcmom = mcstartmom.Mag();
858  double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
859  double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
860 
861  mchists.fHmcstartx->Fill(mcstart.X());
862  mchists.fHmcstarty->Fill(mcstart.Y());
863  mchists.fHmcstartz->Fill(mcstart.Z());
864  mchists.fHmcendx->Fill(mcend.X());
865  mchists.fHmcendy->Fill(mcend.Y());
866  mchists.fHmcendz->Fill(mcend.Z());
867  mchists.fHmctheta->Fill(mcstartmom.Theta());
868  mchists.fHmcphi->Fill(mcstartmom.Phi());
869  mchists.fHmctheta_xz->Fill(mctheta_xz);
870  mchists.fHmctheta_yz->Fill(mctheta_yz);
871  mchists.fHmcmom->Fill(mcmom);
872  mchists.fHmcmoml->Fill(mcmom);
873  mchists.fHmcke->Fill(mcke);
874  mchists.fHmckel->Fill(mcke);
875  mchists.fHmclen->Fill(plen);
876  mchists.fHmclens->Fill(plen);
877  }
878  }
879  }
880  }
881  } //mc
882 
883  // Get tracks and spacepoints and hits
884  art::Handle<std::vector<recob::Track>> trackh;
885  art::Handle<std::vector<art::PtrVector<recob::Track>>> trackvh;
886  art::Handle<std::vector<recob::Hit>> hith;
887 
888  evt.getByLabel(fTrackModuleLabel, trackh);
889  evt.getByLabel(fStitchModuleLabel, trackvh);
890  evt.getByLabel(fHitModuleLabel, hith);
891 
892  // Extract all hits into a vector of art::Ptrs (the format used by back tracker).
893 
894  std::vector<art::Ptr<recob::Hit>> allhits;
895  if (hith.isValid()) {
896  allhits.reserve(hith->size());
897  for (unsigned int i = 0; i < hith->size(); ++i) {
898  allhits.emplace_back(hith, i);
899  }
900  }
901 
902  // Construct FindManyP object to be used for finding track-hit associations.
903 
904  art::FindManyP<recob::Hit> tkhit_find(trackh, evt, fTrackModuleLabel);
905 
906  // This new top part of TrackAna between two long lines of ************s
907  // is particular to analyzing Stitched Tracks.
908  // *******************************************************************//
909 
910  if (trackvh.isValid() && fStitchedAnalysis) {
911  mf::LogDebug("TrackAna") << "TrackAna read " << trackvh->size()
912  << " vectors of Stitched PtrVectorsof tracks.";
913  anaStitch(evt, clockData, detProp);
914  }
915 
916  if (trackh.isValid()) {
917 
918  if (pdump) {
919  *pdump << "\nReconstructed Tracks\n"
920  << " Id MCid x y z dx dy dz "
921  " p\n"
922  << "--------------------------------------------------------------------------------"
923  "-----------\n";
924  }
925 
926  // Loop over tracks.
927 
928  int ntrack = trackh->size();
929  for (int i = 0; i < ntrack; ++i) {
930  art::Ptr<recob::Track> ptrack(trackh, i);
931  const recob::Track& track = *ptrack;
932 
933  // Extract hits associated with this track.
934 
935  std::vector<art::Ptr<recob::Hit>> trackhits;
936  tkhit_find.get(i, trackhits);
937 
938  // Calculate the x offset due to nonzero reconstructed time.
939 
940  double trackdx = 0;
941 
942  // Fill histograms involving reco tracks only.
943 
944  int ntraj = track.NumberTrajectoryPoints();
945  if (ntraj > 0) {
946  TVector3 pos = track.Vertex<TVector3>();
947  TVector3 dir = track.VertexDirection<TVector3>();
948  TVector3 end = track.End<TVector3>();
949  pos[0] += trackdx;
950  end[0] += trackdx;
951 
952  double dpos = bdist(pos);
953  double dend = bdist(end);
954  double tlen = length(track);
955  double theta_xz = std::atan2(dir.X(), dir.Z());
956  double theta_yz = std::atan2(dir.Y(), dir.Z());
957 
958  if (fRecoHistMap.count(0) == 0) fRecoHistMap.emplace(0, RecoHists{"Reco"});
959  const RecoHists& rhists = fRecoHistMap.at(0);
960 
961  rhists.fHstartx->Fill(pos.X());
962  rhists.fHstarty->Fill(pos.Y());
963  rhists.fHstartz->Fill(pos.Z());
964  rhists.fHstartd->Fill(dpos);
965  rhists.fHendx->Fill(end.X());
966  rhists.fHendy->Fill(end.Y());
967  rhists.fHendz->Fill(end.Z());
968  rhists.fHendd->Fill(dend);
969  rhists.fHtheta->Fill(dir.Theta());
970  rhists.fHphi->Fill(dir.Phi());
971  rhists.fHtheta_xz->Fill(theta_xz);
972  rhists.fHtheta_yz->Fill(theta_yz);
973 
974  double mom = 0.;
975  if (track.HasMomentum()) mom = track.VertexMomentum();
976  rhists.fHmom->Fill(mom);
977  rhists.fHmoml->Fill(mom);
978  rhists.fHlen->Fill(tlen);
979  rhists.fHlens->Fill(tlen);
980 
981  // Id of matching mc particle.
982 
983  int mcid = -1;
984 
985  // Loop over direction.
986 
987  for (int swap = 0; swap < 2; ++swap) {
988 
989  // Analyze reversed tracks only if start momentum = end momentum.
990 
991  if (swap != 0 && track.HasMomentum() &&
992  std::abs(track.VertexMomentum() - track.EndMomentum()) > 1.e-3)
993  continue;
994 
995  // Calculate the global-to-local rotation matrix.
996 
997  int start_point = (swap == 0 ? 0 : ntraj - 1);
998  TMatrixD rot = track.GlobalToLocalRotationAtPoint<TMatrixD>(start_point);
999 
1000  // Update track data for reversed case.
1001 
1002  if (swap != 0) {
1003  rot(1, 0) = -rot(1, 0);
1004  rot(2, 0) = -rot(2, 0);
1005  rot(1, 1) = -rot(1, 1);
1006  rot(2, 1) = -rot(2, 1);
1007  rot(1, 2) = -rot(1, 2);
1008  rot(2, 2) = -rot(2, 2);
1009 
1010  pos = track.End<TVector3>();
1011  dir = -track.EndDirection<TVector3>();
1012  end = track.Vertex<TVector3>();
1013  pos[0] += trackdx;
1014  end[0] += trackdx;
1015 
1016  dpos = bdist(pos);
1017  dend = bdist(end);
1018  theta_xz = std::atan2(dir.X(), dir.Z());
1019  theta_yz = std::atan2(dir.Y(), dir.Z());
1020 
1021  if (track.HasMomentum()) mom = track.EndMomentum();
1022  }
1023 
1024  // Get covariance matrix.
1025 
1026  const auto& cov = (swap == 0 ? track.VertexCovariance() : track.EndCovariance());
1027 
1028  // Loop over track-like mc particles.
1029 
1030  for (unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1031  const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1032  int pdg = mctrk.PdgCode();
1033  if (fIgnoreSign) pdg = std::abs(pdg);
1034  auto iMCHistMap = fMCHistMap.find(pdg);
1035  if (iMCHistMap == fMCHistMap.end())
1036  throw cet::exception("TrackAna") << "no particle with ID=" << pdg << "\n";
1037  const MCHists& mchists = iMCHistMap->second;
1038 
1039  // Calculate the x offset due to nonzero mc particle time.
1040 
1041  double mctime = mctrk.Start().T(); // nsec
1042  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
1043 
1044  // Calculate the points where this mc particle enters and leaves the
1045  // fiducial volume, and the length in the fiducial volume.
1046 
1047  TVector3 mcstart;
1048  TVector3 mcend;
1049  TVector3 mcstartmom;
1050  TVector3 mcendmom;
1051  double plen =
1052  length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1053 
1054  // Get the displacement of this mc particle in the global coordinate system.
1055 
1056  TVector3 mcpos = mcstart - pos;
1057 
1058  // Rotate the momentum and position to the
1059  // track-local coordinate system.
1060 
1061  TVector3 mcmoml = rot * mcstartmom;
1062  TVector3 mcposl = rot * mcpos;
1063 
1064  double colinearity = mcmoml.Z() / mcmoml.Mag();
1065 
1066  double u = mcposl.X();
1067  double v = mcposl.Y();
1068  double w = mcposl.Z();
1069 
1070  double pu = mcmoml.X();
1071  double pv = mcmoml.Y();
1072  double pw = mcmoml.Z();
1073 
1074  double dudw = pu / pw;
1075  double dvdw = pv / pw;
1076 
1077  double u0 = u - w * dudw;
1078  double v0 = v - w * dvdw;
1079  double uv0 = std::hypot(u0, v0);
1080 
1081  mchists.fHduvcosth->Fill(colinearity, uv0);
1082 
1083  if (std::abs(uv0) < fMatchDisp) {
1084 
1085  // Fill slope matching histograms.
1086 
1087  mchists.fHmcdudw->Fill(dudw);
1088  mchists.fHmcdvdw->Fill(dvdw);
1089  mchists.fHdudwpull->Fill(dudw / std::sqrt(cov(2, 2)));
1090  mchists.fHdvdwpull->Fill(dvdw / std::sqrt(cov(3, 3)));
1091  }
1092  mchists.fHcosth->Fill(colinearity);
1093  if (colinearity > fMatchColinearity) {
1094 
1095  // Fill displacement matching histograms.
1096 
1097  mchists.fHmcu->Fill(u0);
1098  mchists.fHmcv->Fill(v0);
1099  mchists.fHmcw->Fill(w);
1100  mchists.fHupull->Fill(u0 / std::sqrt(cov(0, 0)));
1101  mchists.fHvpull->Fill(v0 / std::sqrt(cov(1, 1)));
1102 
1103  if (std::abs(uv0) < fMatchDisp) {
1104 
1105  // Fill matching histograms.
1106 
1107  double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
1108  double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
1109  double mcmom = mcstartmom.Mag();
1110  double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
1111  double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
1112 
1113  mchists.fHstartdx->Fill(pos.X() - mcstart.X());
1114  mchists.fHstartdy->Fill(pos.Y() - mcstart.Y());
1115  mchists.fHstartdz->Fill(pos.Z() - mcstart.Z());
1116  mchists.fHenddx->Fill(end.X() - mcend.X());
1117  mchists.fHenddy->Fill(end.Y() - mcend.Y());
1118  mchists.fHenddz->Fill(end.Z() - mcend.Z());
1119  mchists.fHlvsl->Fill(plen, tlen);
1120  mchists.fHdl->Fill(tlen - plen);
1121  mchists.fHpvsp->Fill(mcmom, mom);
1122  double dp = mom - mcmom;
1123  mchists.fHdp->Fill(dp);
1124  mchists.fHppull->Fill(dp / std::sqrt(cov(4, 4)));
1125  if (std::abs(dpos) >= 5. && std::abs(dend) >= 5.) {
1126  mchists.fHpvspc->Fill(mcmom, mom);
1127  mchists.fHdpc->Fill(dp);
1128  mchists.fHppullc->Fill(dp / std::sqrt(cov(4, 4)));
1129  }
1130 
1131  // Count this track as well-reconstructed if it is matched to an
1132  // mc particle (which it is if get here), and if in addition the
1133  // starting position (w) matches and the reconstructed track length
1134  // is more than 0.5 of the mc particle trajectory length.
1135 
1136  bool good = std::abs(w) <= fWMatchDisp && tlen > fMatchLength * plen;
1137  if (good) {
1138  mcid = mctrk.TrackID();
1139 
1140  // Calculate and fill hit efficiency and purity.
1141 
1142  std::set<int> tkidset;
1143  tkidset.insert(mcid);
1144  double hiteff = bt_serv->HitCollectionEfficiency(
1145  clockData, tkidset, trackhits, allhits, geo::k3D);
1146  double hitpurity = bt_serv->HitCollectionPurity(clockData, tkidset, trackhits);
1147  mchists.fHHitEff->Fill(hiteff);
1148  mchists.fHHitPurity->Fill(hitpurity);
1149 
1150  // Fill efficiency numerator histograms.
1151 
1152  mchists.fHgstartx->Fill(mcstart.X());
1153  mchists.fHgstarty->Fill(mcstart.Y());
1154  mchists.fHgstartz->Fill(mcstart.Z());
1155  mchists.fHgendx->Fill(mcend.X());
1156  mchists.fHgendy->Fill(mcend.Y());
1157  mchists.fHgendz->Fill(mcend.Z());
1158  mchists.fHgtheta->Fill(mcstartmom.Theta());
1159  mchists.fHgphi->Fill(mcstartmom.Phi());
1160  mchists.fHgtheta_xz->Fill(mctheta_xz);
1161  mchists.fHgtheta_yz->Fill(mctheta_yz);
1162  mchists.fHgmom->Fill(mcmom);
1163  mchists.fHgmoml->Fill(mcmom);
1164  mchists.fHgke->Fill(mcke);
1165  mchists.fHgkel->Fill(mcke);
1166  mchists.fHglen->Fill(plen);
1167  mchists.fHglens->Fill(plen);
1168 
1169  // set the match flag
1170  selected_mctracks[imc].second = i;
1171 
1172  if (fPrintLevel > 0) {
1173  const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mcid);
1174  float KE = ptkl->E() - ptkl->Mass();
1175  std::string KEUnits = " GeV";
1176  if (mctrk.Origin() != simb::kCosmicRay) {
1177  // MeV for low energy particles
1178  KE *= 1000;
1179  KEUnits = " MeV";
1180  }
1181  mf::LogVerbatim("TrackAna")
1182  << evt.run() << "." << evt.event() << " Match MCTkID " << std::setw(6)
1183  << mctrk.TrackID() << " Origin " << mctrk.Origin() << " PDG" << std::setw(5)
1184  << mctrk.PdgCode() << " KE" << std::setw(4) << (int)KE << KEUnits
1185  << " RecoTrkID " << track.ID() << " hitEff " << std::setprecision(2)
1186  << hiteff << " hitPur " << hitpurity;
1187  int sWire, sTick, eWire, eTick;
1188  // this won't work for DUNE
1189  for (unsigned short ipl = 0; ipl < geom->Nplanes(); ++ipl) {
1190  sWire = geom->NearestWire(mcstart, ipl, 0, 0);
1191  sTick = detProp.ConvertXToTicks(mcstart[0], ipl, 0, 0);
1192  eWire = geom->NearestWire(mcend, ipl, 0, 0);
1193  eTick = detProp.ConvertXToTicks(mcend[0], ipl, 0, 0);
1194  mf::LogVerbatim("TrackAna")
1195  << " Wire:Tick in Pln " << ipl << " W:T " << sWire << ":" << sTick
1196  << " - " << eWire << ":" << eTick;
1197  } // ipl
1198  } // fPrintLevel == 2
1199  } // good
1200  }
1201  }
1202  }
1203  }
1204 
1205  // Dump track information here.
1206 
1207  if (pdump) {
1208  TVector3 pos = track.Vertex<TVector3>();
1209  TVector3 dir = track.VertexDirection<TVector3>();
1210  TVector3 end = track.End<TVector3>();
1211  pos[0] += trackdx;
1212  end[0] += trackdx;
1213  TVector3 enddir = track.EndDirection<TVector3>();
1214  double pstart = track.VertexMomentum();
1215  double pend = track.EndMomentum();
1216  *pdump << "\nOffset" << std::setw(3) << track.ID() << std::setw(6) << mcid << " "
1217  << std::fixed << std::setprecision(2) << std::setw(10) << trackdx << "\nStart "
1218  << std::setw(3) << track.ID() << std::setw(6) << mcid << " " << std::fixed
1219  << std::setprecision(2) << std::setw(10) << pos[0] << std::setw(10) << pos[1]
1220  << std::setw(10) << pos[2];
1221  if (pstart > 0.) {
1222  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
1223  << std::setw(10) << dir[1] << std::setw(10) << dir[2];
1224  }
1225  else
1226  *pdump << std::setw(32) << " ";
1227  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
1228  *pdump << "\nEnd " << std::setw(3) << track.ID() << std::setw(6) << mcid << " "
1229  << std::fixed << std::setprecision(2) << std::setw(10) << end[0] << std::setw(10)
1230  << end[1] << std::setw(10) << end[2];
1231  if (pend > 0.01) {
1232  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10) << enddir[0]
1233  << std::setw(10) << enddir[1] << std::setw(10) << enddir[2];
1234  }
1235  else
1236  *pdump << std::setw(32) << " ";
1237  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
1238  << "\nLength: " << tlen << "\n";
1239  }
1240  }
1241  }
1242  } // i
1243 
1244  // print out un-matched MC tracks
1245  if (fPrintLevel > 0) {
1246  for (unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1247  if (selected_mctracks[imc].second >= 0) continue;
1248  const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1249  const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mctrk.TrackID());
1250  float KE = ptkl->E() - ptkl->Mass();
1251  std::string KEUnits = " GeV";
1252  if (mctrk.Origin() != simb::kCosmicRay) {
1253  // MeV for low energy particles
1254  KE *= 1000;
1255  KEUnits = " MeV";
1256  }
1257  // find the start/end wire:time in each plane
1258  TVector3 mcstart, mcend, mcstartmom, mcendmom;
1259  double mcdx = mctrk.Start().T() * 1.e-3 * detProp.DriftVelocity(); // cm
1260  double plen = length(clockData, detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1261  mf::LogVerbatim("TrackAna")
1262  << evt.run() << "." << evt.event() << " NoMat MCTkID " << std::setw(6) << mctrk.TrackID()
1263  << " Origin " << mctrk.Origin() << " PDG" << std::setw(5) << mctrk.PdgCode() << " KE"
1264  << std::setw(4) << (int)KE << KEUnits << " Length " << std::fixed << std::setprecision(1)
1265  << plen << " cm";
1266  if (fPrintLevel > 1) {
1267  int sWire, sTick, eWire, eTick;
1268  // this won't work for DUNE
1269  for (unsigned short ipl = 0; ipl < geom->Nplanes(); ++ipl) {
1270  sWire = geom->NearestWire(mcstart, ipl, 0, 0);
1271  sTick = detProp.ConvertXToTicks(mcstart[0], ipl, 0, 0);
1272  eWire = geom->NearestWire(mcend, ipl, 0, 0);
1273  eTick = detProp.ConvertXToTicks(mcend[0], ipl, 0, 0);
1274  mf::LogVerbatim("TrackAna") << " Wire:Tick in Pln " << ipl << " W:T " << sWire << ":"
1275  << sTick << " - " << eWire << ":" << eTick;
1276  } // ipl
1277  } // fPrintLevel > 1
1278  } // imc
1279  } // fPrintLevel > 0
1280  }
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
sim::MCTrack::Origin
simb::Origin_t Origin() const
Definition: MCTrack.h:40
recob::Track::EndMomentum
double EndMomentum() const
Definition: Track.h:144
art::DataViewImpl::event
EventNumber_t event() const
Definition: DataViewImpl.cc:85
recob::Track::VertexMomentum
double VertexMomentum() const
Definition: Track.h:142
string
std::string string
Definition: nybbler.cc:12
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService_service.cc:128
cheat::BackTrackerService::HitCollectionEfficiency
double HitCollectionEfficiency(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits, std::vector< art::Ptr< recob::Hit >> const &allhits, geo::View_t const &view) const
Definition: BackTrackerService_service.cc:354
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
recob::Track::HasMomentum
bool HasMomentum() const
Definition: Track.h:140
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
recob::Track::VertexDirection
Vector_t VertexDirection() const
Definition: Track.h:132
recob::Track::VertexCovariance
const SMatrixSym55 & VertexCovariance() const
Definition: Track.h:154
trkf::TrackAna::fStitchModuleLabel
std::string fStitchModuleLabel
Definition: TrackAna_module.cc:415
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
sim::MCStep::T
double T() const
Definition: MCStep.h:45
dir
string dir
Definition: INukeNucleonCorr.cxx:61
setprecision
Q_EXPORT QTSManip setprecision(int p)
Definition: qtextstream.h:343
trkf::TrackAna::fMCHistMap
std::map< int, MCHists > fMCHistMap
Definition: TrackAna_module.cc:437
cet::hypot
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
geo::k3D
3-dimensional objects, potentially hits, clusters, prongs, etc.
Definition: geo_types.h:135
recob::Track::GlobalToLocalRotationAtPoint
Rotation_t GlobalToLocalRotationAtPoint(size_t p) const
Definition: Track.h:190
art::Handle::isValid
bool isValid() const noexcept
Definition: Handle.h:191
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2848
e
const double e
Definition: gUpMuFluxGen.cxx:165
art::swap
void swap(Handle< T > &a, Handle< T > &b)
cheat::BackTrackerService::HitCollectionPurity
double HitCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
Definition: BackTrackerService_service.cc:332
geo::GeometryCore::NearestWire
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
Definition: GeometryCore.cxx:1141
trkf::TrackAna::anaStitch
void anaStitch(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp)
Definition: TrackAna_module.cc:1283
recob::Track::Vertex
Point_t const & Vertex() const
Definition: Track.h:124
recob::Track::EndCovariance
const SMatrixSym55 & EndCovariance() const
Definition: Track.h:155
trkf::TrackAna::fMCTrackModuleLabel
std::string fMCTrackModuleLabel
Definition: TrackAna_module.cc:413
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
sim::MCTrack::PdgCode
int PdgCode() const
Definition: MCTrack.h:41
sim::MCStep::Momentum
const TLorentzVector & Momentum() const
Definition: MCStep.h:38
setw
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
trkf::TrackAna::fOriginValue
simb::Origin_t fOriginValue
Definition: TrackAna_module.cc:432
recob::Track::ID
int ID() const
Definition: Track.h:198
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:208
recob::Track::EndDirection
Vector_t EndDirection() const
Definition: Track.h:133
trkf::TrackAna::fRecoHistMap
std::map< int, RecoHists > fRecoHistMap
Definition: TrackAna_module.cc:438
sim::MCStep::E
double E() const
Definition: MCStep.h:49
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
recob::Track::End
Point_t const & End() const
Definition: Track.h:125
sim::MCTrack::Start
const MCStep & Start() const
Definition: MCTrack.h:44
trkf::TrackAna::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: TrackAna_module.cc:412
sim::MCTrack::TrackID
unsigned int TrackID() const
Definition: MCTrack.h:42
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
trkf::TrackAna::fHitModuleLabel
std::string fHitModuleLabel
Definition: TrackAna_module.cc:418
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
simb::kCosmicRay
Cosmic rays.
Definition: MCTruth.h:24
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.

1286  {
1287  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
1288  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1289  art::ServiceHandle<geo::Geometry const> geom;
1290 
1291  std::map<int, std::map<int, art::PtrVector<recob::Hit>>> hitmap; // trkID, otrk, hitvec
1292  std::map<int, int> KEmap; // length traveled in det [cm]?, trkID want to sort by KE
1293  bool mc = !evt.isRealData();
1294  art::Handle<std::vector<recob::Track>> trackh;
1295  art::Handle<std::vector<recob::SpacePoint>> sppth;
1296  art::Handle<std::vector<art::PtrVector<recob::Track>>> trackvh;
1297 
1298  evt.getByLabel(fTrackModuleLabel, trackh);
1299  evt.getByLabel(fSpacepointModuleLabel, sppth);
1300  evt.getByLabel(fStitchModuleLabel, trackvh);
1301  int ntv(trackvh->size());
1302 
1303  std::vector<art::PtrVector<recob::Track>>::const_iterator cti = trackvh->begin();
1304 
1305  if (trackh.isValid()) {
1306  art::FindManyP<recob::SpacePoint> fswhole(trackh, evt, fTrkSpptAssocModuleLabel);
1307  int nsppts_assnwhole = fswhole.size();
1308  std::cout << "TrackAna: Number of clumps of Spacepoints from Assn for all Tracks: "
1309  << nsppts_assnwhole << std::endl;
1310  }
1311 
1312  if (fRecoHistMap.count(0) == 0) {
1313  fRecoHistMap.emplace(0, RecoHists{"Reco"});
1314  std::cout << "\n"
1315  << "\t\t TrkAna: Fresh fRecoHistMap[0] ******* \n"
1316  << std::endl;
1317  }
1318  const RecoHists& rhistsStitched = fRecoHistMap.at(0);
1319 
1320  std::vector<std::vector<unsigned int>> NtrkIdsAll;
1321  std::vector<double> ntvsorted;
1322  hitmap.clear();
1323  KEmap.clear();
1324 
1325  // Look at the components of the stitched tracks. Grab their sppts/hits from Assns.
1326  for (int o = 0; o < ntv; ++o) // o for outer
1327  {
1328  const art::PtrVector<recob::Track> pvtrack(*(cti++));
1329  int ntrack = pvtrack.size();
1330  std::vector<std::vector<unsigned int>> NtrkId_Hit; // hit IDs in inner tracks
1331  std::vector<unsigned int> vecMode;
1332  art::FindManyP<recob::SpacePoint> fs(pvtrack, evt, fTrkSpptAssocModuleLabel);
1333 
1334  for (int i = 0; i < ntrack; ++i) {
1335 
1336  // From gdb> ptype fs, the vector of Ptr<SpacePoint>s it appears is
1337  // grabbed after fs.at(0)
1338  bool assns(true);
1339  try {
1340  // Got Spacepoints from this Track; now get Hits from those
1341  // Spacepoints.
1342 
1343  int nsppts_assn = fs.at(i).size();
1344 
1345  const auto& sppt = fs.at(i); //.at(0);
1346  // since we're in a try and worried about failure, we won't pull the following
1347  // FindManyP out of the loop.
1348  art::FindManyP<recob::Hit> fh(sppt, evt, fHitSpptAssocModuleLabel);
1349 
1350  // Importantly, loop on all sppts, though they don't all contribute to the track.
1351  // As opposed to looping on the trajectory pts, which is a lower number.
1352  // Also, important, in job in whch this runs I set TrackKal3DSPS parameter MaxPass=1,
1353  // cuz I don't want merely the sparse set of sppts as follows from the uncontained
1354  // MS-measurement in 2nd pass.
1355  std::vector<unsigned int> vecNtrkIds;
1356  for (int is = 0; is < nsppts_assn; ++is) {
1357  int nhits = fh.at(is).size(); // should be 2 or 3: number of planes.
1358  for (int ih = 0; ih < nhits; ++ih) {
1359  const auto& hit = fh.at(is).at(ih); // Our vector is after the .at(is) this time.
1360  if (hit->SignalType() != geo::kCollection) continue;
1361  rhistsStitched.fHHitChg->Fill(hit->Integral());
1362  rhistsStitched.fHHitWidth->Fill(2. * hit->RMS());
1363  if (mc) {
1364  std::vector<sim::TrackIDE> tids = bt_serv->HitToTrackIDEs(clockData, hit);
1365  // more here.
1366  // Loop over track ids.
1367  bool justOne(true); // Only take first trk that contributed to this hit
1368  for (std::vector<sim::TrackIDE>::const_iterator itid = tids.begin();
1369  itid != tids.end();
1370  ++itid) {
1371  int trackID = std::abs(itid->trackID);
1372  hitmap[trackID][o].push_back(hit);
1373 
1374  if (justOne) {
1375  vecNtrkIds.push_back(trackID);
1376  justOne = false;
1377  }
1378  // Add hit to PtrVector corresponding to this track id.
1379  rhistsStitched.fHHitTrkId->Fill(trackID);
1380  const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(trackID);
1381  if (!part) break;
1382 
1383  rhistsStitched.fHHitPdg->Fill(part->PdgCode());
1384  // This really needs to be indexed as KE deposited in volTPC, not just KE. EC, 24-July-2014.
1385 
1386  TVector3 mcstart;
1387  TVector3 mcend;
1388  TVector3 mcstartmom;
1389  TVector3 mcendmom;
1390  double mctime = part->T(); // nsec
1391  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
1392 
1393  double plen =
1394  length(clockData, detProp, *part, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1395 
1396  KEmap[(int)(1e6 * plen)] = trackID; // multiple assignment but always the same, so
1397  // fine.
1398  }
1399 
1400  } // mc
1401  } // hits
1402 
1403  } // spacepoints
1404 
1405  if (mc) {
1406  NtrkId_Hit.push_back(vecNtrkIds);
1407  // Find the trkID mode for this i^th track
1408  unsigned int ii(1);
1409  int max(-12), n(1), ind(0);
1410  std::sort(vecNtrkIds.begin(), vecNtrkIds.end());
1411  std::vector<unsigned int> strkIds(vecNtrkIds);
1412  while (ii < vecNtrkIds.size()) {
1413  if (strkIds.at(ii) != strkIds.at(ii - 1)) { n = 1; }
1414  else {
1415  n++;
1416  }
1417  if (n > max) {
1418  max = n;
1419  ind = ii;
1420  }
1421  ii++;
1422  }
1423  unsigned int mode(sim::NoParticleId);
1424  if (strkIds.begin() != strkIds.end()) mode = strkIds.at(ind);
1425  vecMode.push_back(mode);
1426 
1427  if (strkIds.size() != 0)
1428  rhistsStitched.fModeFrac->Fill((double)max / (double)strkIds.size());
1429  else
1430  rhistsStitched.fModeFrac->Fill(-1.0);
1431  } // mc
1432 
1433  } // end try
1434  catch (cet::exception& x) {
1435  assns = false;
1436  }
1437  if (!assns) throw cet::exception("TrackAna") << "Bad Associations. \n";
1438 
1439  } // i
1440 
1441  if (mc) {
1442  // one vector per o trk, for all modes of stitched i trks
1443  NtrkIdsAll.push_back(vecMode);
1444 
1445  std::unique(NtrkIdsAll.back().begin(), NtrkIdsAll.back().end());
1446  double sum(0.0);
1447  for (auto const val : NtrkIdsAll.back()) {
1448  sum += hitmap[val][o].size();
1449  }
1450  ntvsorted.push_back(sum);
1451  }
1452 
1453  //
1454  } // o
1455 
1456  int vtmp(0);
1457  // get KEmap indices by most energetic first, least last.
1458  for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1459  vtmp++;
1460  }
1461 
1462  // get o trk indices by their hits. Most populous track first, least last.
1463  for (auto const o : fsort_indexes(ntvsorted)) {
1464  int v(0);
1465  // get KEmap indices by longest trajectory first, least last.
1466  for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1467  int val = it->second; // grab trkIDs in order, since they're sorted by KE
1468  rhistsStitched.fNTrkIdTrks3->Fill(o, v, hitmap[val][o].size());
1469  v++;
1470  }
1471  }
1472 
1473  // In how many o tracks did each trkId appear? Histo it. Would like it to be
1474  // precisely 1. Histo it vs. particle KE.
1475  flattener flat(NtrkIdsAll);
1476  std::vector<unsigned int>& v = flat;
1477  // auto const it ( std::unique(v.begin(),v.end()) ); // never use this it,
1478  // perhaps.
1479  for (auto const val : v) {
1480  if (val != (unsigned int)sim::NoParticleId) {
1481  const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(val);
1482  double T(part->E() - 0.001 * part->Mass());
1483  rhistsStitched.fNTrkIdTrks->Fill(std::count(v.begin(), v.end(), val));
1484  rhistsStitched.fNTrkIdTrks2->Fill(std::count(v.begin(), v.end(), val), T);
1485  }
1486  else {
1487  rhistsStitched.fNTrkIdTrks2->Fill(-1.0, 0.0);
1488  }
1489  }
1490  }
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
cheat::BackTrackerService::HitToTrackIDEs
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTrackerService_service.cc:173
trkf::TrackAna::fHitSpptAssocModuleLabel
std::string fHitSpptAssocModuleLabel
Definition: TrackAna_module.cc:417
trkf::TrackAna::fTrkSpptAssocModuleLabel
std::string fTrkSpptAssocModuleLabel
Definition: TrackAna_module.cc:416
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService_service.cc:128
trkf::TrackAna::fSpacepointModuleLabel
std::string fSpacepointModuleLabel
Definition: TrackAna_module.cc:414
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
genie::units::fs
static constexpr double fs
Definition: Units.h:100
trkf::TrackAna::fStitchModuleLabel
std::string fStitchModuleLabel
Definition: TrackAna_module.cc:415
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
art::Handle::isValid
bool isValid() const noexcept
Definition: Handle.h:191
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
trkf::TrackAna::fsort_indexes
std::vector< size_t > fsort_indexes(const std::vector< double > &v)
Definition: TrackAna_module.cc:1526
sim::NoParticleId
static const int NoParticleId
Definition: sim.h:28
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
trkf::TrackAna::fRecoHistMap
std::map< int, RecoHists > fRecoHistMap
Definition: TrackAna_module.cc:438
reco_momentum_tuples.count
int count
Definition: reco_momentum_tuples.py:202
train.x
list x
Definition: train.py:276
trkf::TrackAna::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: TrackAna_module.cc:412
pandora_metadata.ih
ih
Definition: pandora_metadata.py:81
make_pion_angle.flat
def flat(h_in)
Definition: make_pion_angle.py:38
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
geo::kCollection
Signal from collection planes.
Definition: geo_types.h:146
void trkf::TrackAna::endJob ( )
overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 1493 of file TrackAna_module.cc.

1494  {
1495  // Print summary.
1496 
1497  mf::LogInfo("TrackAna") << "TrackAna statistics:\n"
1498  << " Number of events = " << fNumEvent;
1499 
1500  // Fill efficiency histograms.
1501 
1502  for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
1503  ++i) {
1504  const MCHists& mchists = i->second;
1505  effcalc(mchists.fHgstartx, mchists.fHmcstartx, mchists.fHestartx);
1506  effcalc(mchists.fHgstarty, mchists.fHmcstarty, mchists.fHestarty);
1507  effcalc(mchists.fHgstartz, mchists.fHmcstartz, mchists.fHestartz);
1508  effcalc(mchists.fHgendx, mchists.fHmcendx, mchists.fHeendx);
1509  effcalc(mchists.fHgendy, mchists.fHmcendy, mchists.fHeendy);
1510  effcalc(mchists.fHgendz, mchists.fHmcendz, mchists.fHeendz);
1511  effcalc(mchists.fHgtheta, mchists.fHmctheta, mchists.fHetheta);
1512  effcalc(mchists.fHgphi, mchists.fHmcphi, mchists.fHephi);
1513  effcalc(mchists.fHgtheta_xz, mchists.fHmctheta_xz, mchists.fHetheta_xz);
1514  effcalc(mchists.fHgtheta_yz, mchists.fHmctheta_yz, mchists.fHetheta_yz);
1515  effcalc(mchists.fHgmom, mchists.fHmcmom, mchists.fHemom);
1516  effcalc(mchists.fHgmoml, mchists.fHmcmoml, mchists.fHemoml);
1517  effcalc(mchists.fHgke, mchists.fHmcke, mchists.fHeke);
1518  effcalc(mchists.fHgkel, mchists.fHmckel, mchists.fHekel);
1519  effcalc(mchists.fHglen, mchists.fHmclen, mchists.fHelen);
1520  effcalc(mchists.fHglens, mchists.fHmclens, mchists.fHelens);
1521  }
1522  }
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
trkf::TrackAna::fMCHistMap
std::map< int, MCHists > fMCHistMap
Definition: TrackAna_module.cc:437
std::vector< size_t > trkf::TrackAna::fsort_indexes ( const std::vector< double > &  v)
private

Definition at line 1526 of file TrackAna_module.cc.

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

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:
Generated by   doxygen 1.8.11