539 std::unique_ptr<mf::LogInfo> pdump;
542 pdump = std::make_unique<mf::LogInfo>(
"TrackAna");
556 std::map<const recob::Seed*, int> seedmap;
568 *pdump <<
"MC Tracks\n" 569 <<
" Id pdg x y z dx dy dz " 571 <<
"--------------------------------------------------------------------------------" 580 for (
auto const& mctrk : *mctrackh) {
581 int pdg = mctrk.PdgCode();
595 if (mctrk.Start().E() < mctrk.Start().Momentum().Mag() + 1000. *
fMinMCKE) {
continue; }
599 double mctime = mctrk.Start().T();
600 double mcdx = mctime * 1.e-3 * detProp.DriftVelocity();
608 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
613 if (plen <= std::max(
fMinMCLen, 0.)) {
continue; }
618 double pstart = mcstartmom.Mag();
619 double pend = mcendmom.Mag();
620 *pdump <<
"\nOffset" << std::setw(3) << mctrk.TrackID() << std::setw(6) << mctrk.PdgCode()
621 <<
" " << std::fixed << std::setprecision(2) << std::setw(10) << mcdx
622 <<
"\nStart " << std::setw(3) << mctrk.TrackID() << std::setw(6) << mctrk.PdgCode()
623 <<
" " << std::fixed << std::setprecision(2) << std::setw(10) << mcstart[0]
624 << std::setw(10) << mcstart[1] << std::setw(10) << mcstart[2];
626 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
627 << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
628 << std::setw(10) << mcstartmom[2] / pstart;
631 *pdump << std::setw(32) <<
" ";
632 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
633 *pdump <<
"\nEnd " << std::setw(3) << mctrk.TrackID() << std::setw(6) << mctrk.PdgCode()
634 <<
" " << std::fixed << std::setprecision(2) << std::setw(10) << mcend[0]
635 << std::setw(10) << mcend[1] << std::setw(10) << mcend[2];
637 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
638 << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend << std::setw(10)
639 << mcendmom[2] / pend;
642 *pdump << std::setw(32) <<
" ";
643 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend <<
"\n";
649 std::ostringstream ostr;
655 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
656 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
658 mchists.fHmcstartx->Fill(mcstart.X());
659 mchists.fHmcstarty->Fill(mcstart.Y());
660 mchists.fHmcstartz->Fill(mcstart.Z());
661 mchists.fHmcendx->Fill(mcend.X());
662 mchists.fHmcendy->Fill(mcend.Y());
663 mchists.fHmcendz->Fill(mcend.Z());
664 mchists.fHmctheta->Fill(mcstartmom.Theta());
665 mchists.fHmcphi->Fill(mcstartmom.Phi());
666 mchists.fHmctheta_xz->Fill(mctheta_xz);
667 mchists.fHmctheta_yz->Fill(mctheta_yz);
668 mchists.fHmcmom->Fill(mcstartmom.Mag());
669 mchists.fHmclen->Fill(plen);
674 if (!seedh.
isValid()) {
continue; }
678 int nseed = seedh->size();
679 for (
int i = 0; i < nseed; ++i) {
682 if (seedmap.count(&seed) == 0) seedmap[&seed] = -1;
696 double dirmag = dir.Mag();
697 double diryz = std::sqrt(dir.Y() * dir.Y() + dir.Z() * dir.Z());
699 double sinth = dir.X() / dirmag;
700 double costh = diryz / dirmag;
704 sinphi = -dir.Y() / diryz;
705 cosphi = dir.Z() / diryz;
710 rot(0, 1) = sinth * sinphi;
712 rot(2, 1) = -costh * sinphi;
713 rot(0, 2) = -sinth * cosphi;
715 rot(2, 2) = costh * cosphi;
719 int itraj = mcmatch(detProp, mctrk, seed);
720 if (itraj < 0) {
continue; }
724 TVector3 mcpos = mctrk[itraj].Position().Vect() - pos;
725 TVector3 mcmom = mctrk[itraj].Momentum().Vect();
731 TVector3 mcmoml = rot * mcmom;
732 TVector3 mcposl = rot * mcpos;
734 if (mcmoml.Z() < 0.) mcmoml = -mcmoml;
735 costh = mcmoml.Z() / mcmoml.Mag();
737 double u = mcposl.X();
738 double v = mcposl.Y();
739 double w = mcposl.Z();
741 double pu = mcmoml.X();
742 double pv = mcmoml.Y();
743 double pw = mcmoml.Z();
745 double dudw = pu / pw;
746 double dvdw = pv / pw;
748 double u0 = u - w * dudw;
749 double v0 = v - w * dvdw;
750 double uv0 = std::sqrt(u0 * u0 + v0 * v0);
754 mchists.fHduvcosth->Fill(costh, uv0);
759 mchists.fHmcdudw->Fill(dudw);
760 mchists.fHmcdvdw->Fill(dvdw);
762 mchists.fHcosth->Fill(costh);
767 mchists.fHmcu->Fill(u0);
768 mchists.fHmcv->Fill(v0);
769 mchists.fHmcw->Fill(w);
777 seedmap[&
seed] = mctrk.TrackID();
781 mchists.fHmstartx->Fill(mcstart.X());
782 mchists.fHmstarty->Fill(mcstart.Y());
783 mchists.fHmstartz->Fill(mcstart.Z());
784 mchists.fHmendx->Fill(mcend.X());
785 mchists.fHmendy->Fill(mcend.Y());
786 mchists.fHmendz->Fill(mcend.Z());
787 mchists.fHmtheta->Fill(mcstartmom.Theta());
788 mchists.fHmphi->Fill(mcstartmom.Phi());
789 mchists.fHmtheta_xz->Fill(mctheta_xz);
790 mchists.fHmtheta_yz->Fill(mctheta_yz);
791 mchists.fHmmom->Fill(mcstartmom.Mag());
792 mchists.fHmlen->Fill(plen);
798 if (nmatch == 0) {
continue; }
800 mchists.fHgstartx->Fill(mcstart.X());
801 mchists.fHgstarty->Fill(mcstart.Y());
802 mchists.fHgstartz->Fill(mcstart.Z());
803 mchists.fHgendx->Fill(mcend.X());
804 mchists.fHgendy->Fill(mcend.Y());
805 mchists.fHgendz->Fill(mcend.Z());
806 mchists.fHgtheta->Fill(mcstartmom.Theta());
807 mchists.fHgphi->Fill(mcstartmom.Phi());
808 mchists.fHgtheta_xz->Fill(mctheta_xz);
809 mchists.fHgtheta_yz->Fill(mctheta_yz);
810 mchists.fHgmom->Fill(mcstartmom.Mag());
811 mchists.fHglen->Fill(plen);
821 int nseed = seedh->size();
823 if (nseed > 0 && pdump != 0) {
824 *pdump <<
"\nReconstructed Seeds\n" 825 <<
" MCid x y z dx dy dz " 827 <<
"--------------------------------------------------------------------------------" 831 for (
int i = 0; i < nseed; ++i) {
842 double mdir = dir.Mag();
843 if (mdir != 0.) { dir *= (1. / mdir); }
845 double dpos = bdist(pos);
846 double theta_xz = std::atan2(dir.X(), dir.Z());
847 double theta_yz = std::atan2(dir.Y(), dir.Z());
852 int mcid = seedmap[&
seed];
853 *pdump << std::setw(15) << mcid <<
" " << std::fixed << std::setprecision(2)
854 << std::setw(10) << pos[0] << std::setw(10) << pos[1] << std::setw(10) << pos[2]
855 <<
" " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
856 << std::setw(10) << dir[1] << std::setw(10) << dir[2] <<
"\n";
864 rhists.fHx->Fill(pos.X());
865 rhists.fHy->Fill(pos.Y());
866 rhists.fHz->Fill(pos.Z());
867 rhists.fHdist->Fill(dpos);
868 rhists.fHtheta->Fill(dir.Theta());
869 rhists.fHphi->Fill(dir.Phi());
870 rhists.fHtheta_xz->Fill(theta_xz);
871 rhists.fHtheta_yz->Fill(theta_yz);
std::string fMCTrackModuleLabel
void GetPoint(double *Pt, double *Err) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::string fSeedModuleLabel
bool isValid() const noexcept
std::map< int, RecoHists > fRecoHistMap
std::map< int, MCHists > fMCHistMap
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
void GetDirection(double *Dir, double *Err) const