All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::SeedAna Class Reference
Inheritance diagram for trkf::SeedAna:

Classes

struct  MCHists
 
struct  RecoHists
 

Public Member Functions

 SeedAna (fhicl::ParameterSet const &pset)
 

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
 

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.

554  : EDAnalyzer(pset)
555  , fSeedModuleLabel(pset.get<std::string>("SeedModuleLabel"))
556  , fMCTrackModuleLabel(pset.get<std::string>("MCTrackModuleLabel"))
557  , fDump(pset.get<int>("Dump"))
558  , fMinMCKE(pset.get<double>("MinMCKE"))
559  , fMinMCLen(pset.get<double>("MinMCLen"))
560  , fMatchColinearity(pset.get<double>("MatchColinearity"))
561  , fMatchDisp(pset.get<double>("MatchDisp"))
562  , fIgnoreSign(pset.get<bool>("IgnoreSign"))
563  , fNumEvent(0)
564  {
565  mf::LogInfo("SeedAna") << "SeedAna configured with the following parameters:\n"
566  << " SeedModuleLabel = " << fSeedModuleLabel << "\n"
567  << " MCTrackModuleLabel = " << fMCTrackModuleLabel << "\n"
568  << " Dump = " << fDump << "\n"
569  << " MinMCKE = " << fMinMCKE << "\n"
570  << " MinMCLen = " << fMinMCLen;
571  }
trkf::SeedAna::fMCTrackModuleLabel
std::string fMCTrackModuleLabel
Definition: SeedAna_module.cc:351
trkf::SeedAna::fSeedModuleLabel
std::string fSeedModuleLabel
Definition: SeedAna_module.cc:350
trkf::SeedAna::fMatchColinearity
double fMatchColinearity
Definition: SeedAna_module.cc:355

Member Function Documentation

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

Definition at line 574 of file SeedAna_module.cc.

580  {
581  ++fNumEvent;
582 
583  // Optional dump stream.
584 
585  std::unique_ptr<mf::LogInfo> pdump;
586  if (fDump > 0) {
587  --fDump;
588  pdump = std::unique_ptr<mf::LogInfo>(new mf::LogInfo("TrackAna"));
589  }
590 
591  // Make sure histograms are booked.
592 
593  bool mc = !evt.isRealData();
594 
595  // Get seed handle.
596 
597  art::Handle<std::vector<recob::Seed>> seedh;
598  evt.getByLabel(fSeedModuleLabel, seedh);
599 
600  // Seed->mc track matching map.
601 
602  std::map<const recob::Seed*, int> seedmap;
603 
604  if (mc) {
605 
606  // Get MCTracks.
607 
608  art::Handle<std::vector<sim::MCTrack>> mctrackh;
609  evt.getByLabel(fMCTrackModuleLabel, mctrackh);
610 
611  // Dump MCTracks.
612 
613  if (pdump) {
614  *pdump << "MC Tracks\n"
615  << " Id pdg x y z dx dy dz "
616  " p\n"
617  << "--------------------------------------------------------------------------------"
618  "-----------\n";
619  }
620 
621  // Loop over mc tracks, and fill histograms that depend only
622  // on mc particles.
623  auto const detProp =
624  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
625 
626  for (std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();
627  imctrk != mctrackh->end();
628  ++imctrk) {
629  const sim::MCTrack& mctrk = *imctrk;
630  int pdg = mctrk.PdgCode();
631  if (fIgnoreSign) pdg = std::abs(pdg);
632 
633  // Ignore everything except stable charged nonshowering particles.
634 
635  int apdg = std::abs(pdg);
636  if (apdg == 13 || // Muon
637  apdg == 211 || // Charged pion
638  apdg == 321 || // Charged kaon
639  apdg == 2212) { // (Anti)proton
640 
641  // Apply minimum energy cut.
642 
643  if (mctrk.Start().E() >= mctrk.Start().Momentum().Mag() + 1000. * fMinMCKE) {
644 
645  // Calculate the x offset due to nonzero mc particle time.
646 
647  double mctime = mctrk.Start().T(); // nsec
648  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
649 
650  // Calculate the length of this mc particle inside the fiducial volume.
651 
652  TVector3 mcstart;
653  TVector3 mcend;
654  TVector3 mcstartmom;
655  TVector3 mcendmom;
656  double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
657 
658  // Apply minimum fiducial length cut. Always reject particles that have
659  // zero length in the tpc regardless of the configured cut.
660 
661  if (plen > 0. && plen > fMinMCLen) {
662 
663  // Dump MC particle information here.
664 
665  if (pdump) {
666  double pstart = mcstartmom.Mag();
667  double pend = mcendmom.Mag();
668  *pdump << "\nOffset" << std::setw(3) << mctrk.TrackID() << std::setw(6)
669  << mctrk.PdgCode() << " " << std::fixed << std::setprecision(2)
670  << std::setw(10) << mcdx << "\nStart " << std::setw(3) << mctrk.TrackID()
671  << std::setw(6) << mctrk.PdgCode() << " " << std::fixed
672  << std::setprecision(2) << std::setw(10) << mcstart[0] << std::setw(10)
673  << mcstart[1] << std::setw(10) << mcstart[2];
674  if (pstart > 0.) {
675  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10)
676  << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
677  << std::setw(10) << mcstartmom[2] / pstart;
678  }
679  else
680  *pdump << std::setw(32) << " ";
681  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
682  *pdump << "\nEnd " << std::setw(3) << mctrk.TrackID() << std::setw(6)
683  << mctrk.PdgCode() << " " << std::fixed << std::setprecision(2)
684  << std::setw(10) << mcend[0] << std::setw(10) << mcend[1] << std::setw(10)
685  << mcend[2];
686  if (pend > 0.01) {
687  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10)
688  << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend
689  << std::setw(10) << mcendmom[2] / pend;
690  }
691  else
692  *pdump << std::setw(32) << " ";
693  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend << "\n";
694  }
695 
696  // Fill histograms.
697 
698  if (fMCHistMap.count(pdg) == 0) {
699  std::ostringstream ostr;
700  ostr << "MC" << (fIgnoreSign ? "All" : (pdg > 0 ? "Pos" : "Neg")) << std::abs(pdg);
701  fMCHistMap.emplace(pdg, MCHists{ostr.str()});
702  }
703  const MCHists& mchists = fMCHistMap.at(pdg);
704 
705  double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
706  double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
707 
708  mchists.fHmcstartx->Fill(mcstart.X());
709  mchists.fHmcstarty->Fill(mcstart.Y());
710  mchists.fHmcstartz->Fill(mcstart.Z());
711  mchists.fHmcendx->Fill(mcend.X());
712  mchists.fHmcendy->Fill(mcend.Y());
713  mchists.fHmcendz->Fill(mcend.Z());
714  mchists.fHmctheta->Fill(mcstartmom.Theta());
715  mchists.fHmcphi->Fill(mcstartmom.Phi());
716  mchists.fHmctheta_xz->Fill(mctheta_xz);
717  mchists.fHmctheta_yz->Fill(mctheta_yz);
718  mchists.fHmcmom->Fill(mcstartmom.Mag());
719  mchists.fHmclen->Fill(plen);
720 
721  // Loop over seeds and do matching.
722 
723  int nmatch = 0;
724  if (seedh.isValid()) {
725 
726  // Loop over seeds.
727 
728  int nseed = seedh->size();
729  for (int i = 0; i < nseed; ++i) {
730  art::Ptr<recob::Seed> pseed(seedh, i);
731  const recob::Seed& seed = *pseed;
732  if (seedmap.count(&seed) == 0) seedmap[&seed] = -1;
733 
734  // Get parameters of this seed.
735 
736  TVector3 pos;
737  TVector3 dir;
738  double err[3];
739  seed.GetPoint(&pos[0], err);
740  seed.GetDirection(&dir[0], err);
741 
742  // Calculate the global-to-local rotation matrix.
743  // Copied from Track.cxx.
744 
745  TMatrixD rot(3, 3);
746  double dirmag = dir.Mag();
747  double diryz = std::sqrt(dir.Y() * dir.Y() + dir.Z() * dir.Z());
748 
749  double sinth = dir.X() / dirmag;
750  double costh = diryz / dirmag;
751  double sinphi = 0.;
752  double cosphi = 1.;
753  if (diryz != 0) {
754  sinphi = -dir.Y() / diryz;
755  cosphi = dir.Z() / diryz;
756  }
757  rot(0, 0) = costh;
758  rot(1, 0) = 0.;
759  rot(2, 0) = sinth;
760  rot(0, 1) = sinth * sinphi;
761  rot(1, 1) = cosphi;
762  rot(2, 1) = -costh * sinphi;
763  rot(0, 2) = -sinth * cosphi;
764  rot(1, 2) = sinphi;
765  rot(2, 2) = costh * cosphi;
766 
767  // Get best matching mc trajectory point.
768 
769  int itraj = mcmatch(detProp, mctrk, seed);
770  if (itraj >= 0) {
771 
772  // Get mc relative position and direction at matching trajectory point.
773 
774  TVector3 mcpos = mctrk[itraj].Position().Vect() - pos;
775  TVector3 mcmom = mctrk[itraj].Momentum().Vect();
776  mcpos[0] += mcdx;
777 
778  // Rotate the momentum and position to the
779  // seed-local coordinate system.
780 
781  TVector3 mcmoml = rot * mcmom;
782  TVector3 mcposl = rot * mcpos;
783 
784  if (mcmoml.Z() < 0.) mcmoml = -mcmoml;
785  double costh = mcmoml.Z() / mcmoml.Mag();
786 
787  double u = mcposl.X();
788  double v = mcposl.Y();
789  double w = mcposl.Z();
790 
791  double pu = mcmoml.X();
792  double pv = mcmoml.Y();
793  double pw = mcmoml.Z();
794 
795  double dudw = pu / pw;
796  double dvdw = pv / pw;
797 
798  double u0 = u - w * dudw;
799  double v0 = v - w * dvdw;
800  double uv0 = std::sqrt(u0 * u0 + v0 * v0);
801 
802  // Fill matching histograms.
803 
804  mchists.fHduvcosth->Fill(costh, uv0);
805  if (std::abs(uv0) < fMatchDisp) {
806 
807  // Fill slope matching histograms.
808 
809  mchists.fHmcdudw->Fill(dudw);
810  mchists.fHmcdvdw->Fill(dvdw);
811  }
812  mchists.fHcosth->Fill(costh);
813  if (costh > fMatchColinearity) {
814 
815  // Fill displacement matching histograms.
816 
817  mchists.fHmcu->Fill(u0);
818  mchists.fHmcv->Fill(v0);
819  mchists.fHmcw->Fill(w);
820 
821  if (std::abs(uv0) < fMatchDisp) {
822 
823  // Now we have passed all matching cuts and we have a matching
824  // mc particle + seed pair.
825 
826  ++nmatch;
827  seedmap[&seed] = mctrk.TrackID();
828 
829  // Fill matched seed histograms (seed multiplicity).
830 
831  mchists.fHmstartx->Fill(mcstart.X());
832  mchists.fHmstarty->Fill(mcstart.Y());
833  mchists.fHmstartz->Fill(mcstart.Z());
834  mchists.fHmendx->Fill(mcend.X());
835  mchists.fHmendy->Fill(mcend.Y());
836  mchists.fHmendz->Fill(mcend.Z());
837  mchists.fHmtheta->Fill(mcstartmom.Theta());
838  mchists.fHmphi->Fill(mcstartmom.Phi());
839  mchists.fHmtheta_xz->Fill(mctheta_xz);
840  mchists.fHmtheta_yz->Fill(mctheta_yz);
841  mchists.fHmmom->Fill(mcstartmom.Mag());
842  mchists.fHmlen->Fill(plen);
843  }
844  }
845  }
846  }
847 
848  // If we found at least one matched seed, fill good
849  // particle histograms.
850 
851  if (nmatch > 0) {
852  mchists.fHgstartx->Fill(mcstart.X());
853  mchists.fHgstarty->Fill(mcstart.Y());
854  mchists.fHgstartz->Fill(mcstart.Z());
855  mchists.fHgendx->Fill(mcend.X());
856  mchists.fHgendy->Fill(mcend.Y());
857  mchists.fHgendz->Fill(mcend.Z());
858  mchists.fHgtheta->Fill(mcstartmom.Theta());
859  mchists.fHgphi->Fill(mcstartmom.Phi());
860  mchists.fHgtheta_xz->Fill(mctheta_xz);
861  mchists.fHgtheta_yz->Fill(mctheta_yz);
862  mchists.fHgmom->Fill(mcstartmom.Mag());
863  mchists.fHglen->Fill(plen);
864  }
865  }
866  }
867  }
868  }
869  }
870  }
871 
872  // Loop over seeds and fill reco-only seed histograms.
873 
874  if (seedh.isValid()) {
875 
876  // Loop over seeds.
877 
878  int nseed = seedh->size();
879 
880  if (nseed > 0 && pdump != 0) {
881  *pdump << "\nReconstructed Seeds\n"
882  << " MCid x y z dx dy dz "
883  " p\n"
884  << "--------------------------------------------------------------------------------"
885  "-----------\n";
886  }
887 
888  for (int i = 0; i < nseed; ++i) {
889  art::Ptr<recob::Seed> pseed(seedh, i);
890  const recob::Seed& seed = *pseed;
891 
892  // Fill histograms involving reco seeds only.
893 
894  TVector3 pos;
895  TVector3 dir;
896  double err[3];
897  seed.GetPoint(&pos[0], err);
898  seed.GetDirection(&dir[0], err);
899  double mdir = dir.Mag();
900  if (mdir != 0.) { dir *= (1. / mdir); }
901 
902  double dpos = bdist(pos);
903  double theta_xz = std::atan2(dir.X(), dir.Z());
904  double theta_yz = std::atan2(dir.Y(), dir.Z());
905 
906  // Dump seed information here.
907 
908  if (pdump) {
909  int mcid = seedmap[&seed];
910  *pdump << std::setw(15) << mcid << " " << std::fixed << std::setprecision(2)
911  << std::setw(10) << pos[0] << std::setw(10) << pos[1] << std::setw(10) << pos[2]
912  << " " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
913  << std::setw(10) << dir[1] << std::setw(10) << dir[2] << "\n";
914  }
915 
916  // Fill histograms.
917 
918  if (fRecoHistMap.count(0) == 0) fRecoHistMap.emplace(0, RecoHists{"Reco"});
919  const RecoHists& rhists = fRecoHistMap.at(0);
920 
921  rhists.fHx->Fill(pos.X());
922  rhists.fHy->Fill(pos.Y());
923  rhists.fHz->Fill(pos.Z());
924  rhists.fHdist->Fill(dpos);
925  rhists.fHtheta->Fill(dir.Theta());
926  rhists.fHphi->Fill(dir.Phi());
927  rhists.fHtheta_xz->Fill(theta_xz);
928  rhists.fHtheta_yz->Fill(theta_yz);
929  }
930  }
931  }
pdg
var pdg
Definition: selectors.fcl:14
trkf::SeedAna::fMCTrackModuleLabel
std::string fMCTrackModuleLabel
Definition: SeedAna_module.cc:351
evd::err
EResult err(const char *call)
Definition: WebEVDServer.cxx:140
recob::Seed::GetPoint
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:108
trkf::SeedAna::fSeedModuleLabel
std::string fSeedModuleLabel
Definition: SeedAna_module.cc:350
sim::MCStep::T
double T() const
Definition: MCStep.h:45
trkf::SeedAna::fMatchColinearity
double fMatchColinearity
Definition: SeedAna_module.cc:355
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
gen_crt_frags.pv
tuple pv
Definition: gen_crt_frags.py:1419
seed
unsigned int seed
Definition: NeutronOsc_module.cc:129
trkf::SeedAna::fRecoHistMap
std::map< int, RecoHists > fRecoHistMap
Definition: SeedAna_module.cc:362
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
sim::MCTrack::PdgCode
int PdgCode() const
Definition: MCTrack.h:41
sim::MCStep::Momentum
const TLorentzVector & Momentum() const
Definition: MCStep.h:38
dropbox.dir
tuple dir
Definition: dropbox.py:28
trkf::SeedAna::fMCHistMap
std::map< int, MCHists > fMCHistMap
Definition: SeedAna_module.cc:361
sim::MCStep::E
double E() const
Definition: MCStep.h:49
sim::MCTrack::Start
const MCStep & Start() const
Definition: MCTrack.h:44
sim::MCTrack::TrackID
unsigned int TrackID() const
Definition: MCTrack.h:42
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
recob::Seed::GetDirection
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:98
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
void trkf::SeedAna::endJob ( )
overrideprivate

Definition at line 934 of file SeedAna_module.cc.

938  {
939  // Print summary.
940 
941  mf::LogInfo("SeedAna") << "SeedAna statistics:\n"
942  << " Number of events = " << fNumEvent;
943 
944  // Fill multiplicity histograms.
945 
946  for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
947  ++i) {
948  const MCHists& mchists = i->second;
949  mulcalc(mchists.fHmstartx, mchists.fHmcstartx, mchists.fHmulstartx);
950  mulcalc(mchists.fHmstarty, mchists.fHmcstarty, mchists.fHmulstarty);
951  mulcalc(mchists.fHmstartz, mchists.fHmcstartz, mchists.fHmulstartz);
952  mulcalc(mchists.fHmendx, mchists.fHmcendx, mchists.fHmulendx);
953  mulcalc(mchists.fHmendy, mchists.fHmcendy, mchists.fHmulendy);
954  mulcalc(mchists.fHmendz, mchists.fHmcendz, mchists.fHmulendz);
955  mulcalc(mchists.fHmtheta, mchists.fHmctheta, mchists.fHmultheta);
956  mulcalc(mchists.fHmphi, mchists.fHmcphi, mchists.fHmulphi);
957  mulcalc(mchists.fHmtheta_xz, mchists.fHmctheta_xz, mchists.fHmultheta_xz);
958  mulcalc(mchists.fHmtheta_yz, mchists.fHmctheta_yz, mchists.fHmultheta_yz);
959  mulcalc(mchists.fHmmom, mchists.fHmcmom, mchists.fHmulmom);
960  mulcalc(mchists.fHmlen, mchists.fHmclen, mchists.fHmullen);
961  }
962  // Fill efficiency histograms.
963 
964  for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
965  ++i) {
966  const MCHists& mchists = i->second;
967  effcalc(mchists.fHgstartx, mchists.fHmcstartx, mchists.fHestartx);
968  effcalc(mchists.fHgstarty, mchists.fHmcstarty, mchists.fHestarty);
969  effcalc(mchists.fHgstartz, mchists.fHmcstartz, mchists.fHestartz);
970  effcalc(mchists.fHgendx, mchists.fHmcendx, mchists.fHeendx);
971  effcalc(mchists.fHgendy, mchists.fHmcendy, mchists.fHeendy);
972  effcalc(mchists.fHgendz, mchists.fHmcendz, mchists.fHeendz);
973  effcalc(mchists.fHgtheta, mchists.fHmctheta, mchists.fHetheta);
974  effcalc(mchists.fHgphi, mchists.fHmcphi, mchists.fHephi);
975  effcalc(mchists.fHgtheta_xz, mchists.fHmctheta_xz, mchists.fHetheta_xz);
976  effcalc(mchists.fHgtheta_yz, mchists.fHmctheta_yz, mchists.fHetheta_yz);
977  effcalc(mchists.fHgmom, mchists.fHmcmom, mchists.fHemom);
978  effcalc(mchists.fHglen, mchists.fHmclen, mchists.fHelen);
979  }
980  }
trkf::SeedAna::fMCHistMap
std::map< int, MCHists > fMCHistMap
Definition: SeedAna_module.cc:361

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