717 std::unique_ptr<mf::LogInfo> pdump;
720 pdump = std::make_unique<mf::LogInfo>(
"TrackAna");
732 std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;
738 if (mc && mctrackh.
isValid()) {
740 selected_mctracks.reserve(mctrackh->size());
745 *pdump <<
"MC Tracks\n" 746 <<
" Id pdg x y z dx dy dz " 748 <<
"--------------------------------------------------------------------------------" 757 imctrk != mctrackh->end();
780 double mctime = mctrk.
Start().
T();
781 double mcdx = mctime * 1.e-3 * detProp.DriftVelocity();
789 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
798 selected_mctracks.push_back(std::make_pair(&mctrk, -1));
803 double pstart = mcstartmom.Mag();
804 double pend = mcendmom.Mag();
805 *pdump <<
"\nOffset" << std::setw(3) << mctrk.
TrackID() << std::setw(6)
806 << mctrk.
PdgCode() <<
" " << std::fixed << std::setprecision(2)
807 << std::setw(10) << mcdx <<
"\nStart " << std::setw(3) << mctrk.
TrackID()
808 << std::setw(6) << mctrk.
PdgCode() <<
" " << std::fixed
809 << std::setprecision(2) << std::setw(10) << mcstart.X() << std::setw(10)
810 << mcstart.Y() << std::setw(10) << mcstart.Z();
812 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
813 << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
814 << std::setw(10) << mcstartmom[2] / pstart;
817 *pdump << std::setw(32) <<
" ";
818 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
819 *pdump <<
"\nEnd " << std::setw(3) << mctrk.
TrackID() << std::setw(6)
820 << mctrk.
PdgCode() <<
" " << std::fixed << std::setprecision(2)
821 << std::setw(10) << mcend.X() << std::setw(10) << mcend.Y() << std::setw(10)
824 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
825 << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend
826 << std::setw(10) << mcendmom[2] / pend;
829 *pdump << std::setw(32) <<
" ";
830 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
831 <<
"\nLength: " << plen <<
"\n";
837 std::ostringstream ostr;
843 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
844 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
845 double mcmom = mcstartmom.Mag();
847 double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
849 mchists.fHmcstartx->Fill(mcstart.X());
850 mchists.fHmcstarty->Fill(mcstart.Y());
851 mchists.fHmcstartz->Fill(mcstart.Z());
852 mchists.fHmcendx->Fill(mcend.X());
853 mchists.fHmcendy->Fill(mcend.Y());
854 mchists.fHmcendz->Fill(mcend.Z());
855 mchists.fHmctheta->Fill(mcstartmom.Theta());
856 mchists.fHmcphi->Fill(mcstartmom.Phi());
857 mchists.fHmctheta_xz->Fill(mctheta_xz);
858 mchists.fHmctheta_yz->Fill(mctheta_yz);
859 mchists.fHmcmom->Fill(mcmom);
860 mchists.fHmcmoml->Fill(mcmom);
861 mchists.fHmcke->Fill(mcke);
862 mchists.fHmckel->Fill(mcke);
863 mchists.fHmclen->Fill(plen);
864 mchists.fHmclens->Fill(plen);
882 std::vector<art::Ptr<recob::Hit>> allhits;
884 allhits.reserve(hith->size());
885 for (
unsigned int i = 0; i < hith->size(); ++i) {
886 allhits.emplace_back(hith, i);
899 mf::LogDebug(
"TrackAna") <<
"TrackAna read " << trackvh->size()
900 <<
" vectors of Stitched PtrVectorsof tracks.";
907 *pdump <<
"\nReconstructed Tracks\n" 908 <<
" Id MCid x y z dx dy dz " 910 <<
"--------------------------------------------------------------------------------" 916 int ntrack = trackh->size();
917 for (
int i = 0; i < ntrack; ++i) {
923 std::vector<art::Ptr<recob::Hit>> trackhits;
924 tkhit_find.get(i, trackhits);
934 TVector3 pos = track.
Vertex<TVector3>();
936 TVector3
end = track.
End<TVector3>();
940 double dpos = bdist(pos);
941 double dend = bdist(end);
942 double tlen = length(track);
943 double theta_xz = std::atan2(dir.X(), dir.Z());
944 double theta_yz = std::atan2(dir.Y(), dir.Z());
949 rhists.fHstartx->Fill(pos.X());
950 rhists.fHstarty->Fill(pos.Y());
951 rhists.fHstartz->Fill(pos.Z());
952 rhists.fHstartd->Fill(dpos);
953 rhists.fHendx->Fill(end.X());
954 rhists.fHendy->Fill(end.Y());
955 rhists.fHendz->Fill(end.Z());
956 rhists.fHendd->Fill(dend);
957 rhists.fHtheta->Fill(dir.Theta());
958 rhists.fHphi->Fill(dir.Phi());
959 rhists.fHtheta_xz->Fill(theta_xz);
960 rhists.fHtheta_yz->Fill(theta_yz);
964 rhists.fHmom->Fill(mom);
965 rhists.fHmoml->Fill(mom);
966 rhists.fHlen->Fill(tlen);
967 rhists.fHlens->Fill(tlen);
975 for (
int swap = 0;
swap < 2; ++
swap) {
985 int start_point = (
swap == 0 ? 0 : ntraj - 1);
991 rot(1, 0) = -rot(1, 0);
992 rot(2, 0) = -rot(2, 0);
993 rot(1, 1) = -rot(1, 1);
994 rot(2, 1) = -rot(2, 1);
995 rot(1, 2) = -rot(1, 2);
996 rot(2, 2) = -rot(2, 2);
998 pos = track.
End<TVector3>();
1000 end = track.
Vertex<TVector3>();
1006 theta_xz = std::atan2(dir.X(), dir.Z());
1007 theta_yz = std::atan2(dir.Y(), dir.Z());
1018 for (
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1019 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1024 throw cet::exception(
"TrackAna") <<
"no particle with ID=" << pdg <<
"\n";
1025 const MCHists& mchists = iMCHistMap->second;
1029 double mctime = mctrk.
Start().
T();
1030 double mcdx = mctime * 1.e-3 * detProp.DriftVelocity();
1037 TVector3 mcstartmom;
1039 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1043 TVector3 mcpos{mcstart.X() - pos[0], mcstart.Y() - pos[1], mcstart.Z() - pos[2]};
1048 TVector3 mcmoml = rot * mcstartmom;
1049 TVector3 mcposl = rot * mcpos;
1051 double colinearity = mcmoml.Z() / mcmoml.Mag();
1053 double u = mcposl.X();
1054 double v = mcposl.Y();
1055 double w = mcposl.Z();
1057 double pu = mcmoml.X();
1058 double pv = mcmoml.Y();
1059 double pw = mcmoml.Z();
1061 double dudw = pu / pw;
1062 double dvdw = pv / pw;
1064 double u0 = u - w * dudw;
1065 double v0 = v - w * dvdw;
1066 double uv0 = std::hypot(u0, v0);
1068 mchists.fHduvcosth->Fill(colinearity, uv0);
1074 mchists.fHmcdudw->Fill(dudw);
1075 mchists.fHmcdvdw->Fill(dvdw);
1076 mchists.fHdudwpull->Fill(dudw / std::sqrt(cov(2, 2)));
1077 mchists.fHdvdwpull->Fill(dvdw / std::sqrt(cov(3, 3)));
1079 mchists.fHcosth->Fill(colinearity);
1084 mchists.fHmcu->Fill(u0);
1085 mchists.fHmcv->Fill(v0);
1086 mchists.fHmcw->Fill(w);
1087 mchists.fHupull->Fill(u0 / std::sqrt(cov(0, 0)));
1088 mchists.fHvpull->Fill(v0 / std::sqrt(cov(1, 1)));
1094 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
1095 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
1096 double mcmom = mcstartmom.Mag();
1098 double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
1100 mchists.fHstartdx->Fill(pos.X() - mcstart.X());
1101 mchists.fHstartdy->Fill(pos.Y() - mcstart.Y());
1102 mchists.fHstartdz->Fill(pos.Z() - mcstart.Z());
1103 mchists.fHenddx->Fill(end.X() - mcend.X());
1104 mchists.fHenddy->Fill(end.Y() - mcend.Y());
1105 mchists.fHenddz->Fill(end.Z() - mcend.Z());
1106 mchists.fHlvsl->Fill(plen, tlen);
1107 mchists.fHdl->Fill(tlen - plen);
1108 mchists.fHpvsp->Fill(mcmom, mom);
1109 double dp = mom - mcmom;
1110 mchists.fHdp->Fill(dp);
1111 mchists.fHppull->Fill(dp / std::sqrt(cov(4, 4)));
1113 mchists.fHpvspc->Fill(mcmom, mom);
1114 mchists.fHdpc->Fill(dp);
1115 mchists.fHppullc->Fill(dp / std::sqrt(cov(4, 4)));
1129 std::set<int> tkidset;
1130 tkidset.insert(mcid);
1132 clockData, tkidset, trackhits, allhits,
geo::k3D);
1134 mchists.fHHitEff->Fill(hiteff);
1135 mchists.fHHitPurity->Fill(hitpurity);
1139 mchists.fHgstartx->Fill(mcstart.X());
1140 mchists.fHgstarty->Fill(mcstart.Y());
1141 mchists.fHgstartz->Fill(mcstart.Z());
1142 mchists.fHgendx->Fill(mcend.X());
1143 mchists.fHgendy->Fill(mcend.Y());
1144 mchists.fHgendz->Fill(mcend.Z());
1145 mchists.fHgtheta->Fill(mcstartmom.Theta());
1146 mchists.fHgphi->Fill(mcstartmom.Phi());
1147 mchists.fHgtheta_xz->Fill(mctheta_xz);
1148 mchists.fHgtheta_yz->Fill(mctheta_yz);
1149 mchists.fHgmom->Fill(mcmom);
1150 mchists.fHgmoml->Fill(mcmom);
1151 mchists.fHgke->Fill(mcke);
1152 mchists.fHgkel->Fill(mcke);
1153 mchists.fHglen->Fill(plen);
1154 mchists.fHglens->Fill(plen);
1157 selected_mctracks[imc].second = i;
1161 float KE = ptkl->
E() - ptkl->
Mass();
1162 std::string KEUnits =
" GeV";
1169 << evt.
run() <<
"." << evt.
event() <<
" Match MCTkID " << std::setw(6)
1170 << mctrk.
TrackID() <<
" Origin " << mctrk.
Origin() <<
" PDG" << std::setw(5)
1171 << mctrk.
PdgCode() <<
" KE" << std::setw(4) << (int)KE << KEUnits
1172 <<
" RecoTrkID " << track.
ID() <<
" hitEff " << std::setprecision(2)
1173 << hiteff <<
" hitPur " << hitpurity;
1175 for (
unsigned short ipl = 0; ipl < geom->
Nplanes(); ++ipl) {
1178 auto const sTick = detProp.ConvertXToTicks(mcstart.X(), planeID);
1180 auto const eTick = detProp.ConvertXToTicks(mcend.X(), planeID);
1182 <<
" Wire:Tick in Pln " << ipl <<
" W:T " << sWire <<
":" << sTick
1183 <<
" - " << eWire <<
":" << eTick;
1195 TVector3 pos = track.
Vertex<TVector3>();
1197 TVector3 end = track.
End<TVector3>();
1203 *pdump <<
"\nOffset" << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " 1204 << std::fixed << std::setprecision(2) << std::setw(10) << trackdx <<
"\nStart " 1205 << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " << std::fixed
1206 << std::setprecision(2) << std::setw(10) << pos[0] << std::setw(10) << pos[1]
1207 << std::setw(10) << pos[2];
1209 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
1210 << std::setw(10) << dir[1] << std::setw(10) << dir[2];
1213 *pdump << std::setw(32) <<
" ";
1214 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
1215 *pdump <<
"\nEnd " << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " 1216 << std::fixed << std::setprecision(2) << std::setw(10) << end[0] << std::setw(10)
1217 << end[1] << std::setw(10) << end[2];
1219 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10) << enddir[0]
1220 << std::setw(10) << enddir[1] << std::setw(10) << enddir[2];
1223 *pdump << std::setw(32) <<
" ";
1224 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
1225 <<
"\nLength: " << tlen <<
"\n";
1233 for (
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1234 if (selected_mctracks[imc].
second >= 0)
continue;
1235 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1237 float KE = ptkl->
E() - ptkl->
Mass();
1238 std::string KEUnits =
" GeV";
1246 TVector3 mcstartmom, mcendmom;
1247 double mcdx = mctrk.
Start().
T() * 1.e-3 * detProp.DriftVelocity();
1248 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1250 << evt.
run() <<
"." << evt.
event() <<
" NoMat MCTkID " << std::setw(6) << mctrk.
TrackID()
1251 <<
" Origin " << mctrk.
Origin() <<
" PDG" << std::setw(5) << mctrk.
PdgCode() <<
" KE" 1252 << std::setw(4) << (int)KE << KEUnits <<
" Length " << std::fixed << std::setprecision(1)
1256 for (
unsigned short ipl = 0; ipl < geom->
Nplanes(); ++ipl) {
1259 auto const sTick = detProp.ConvertXToTicks(mcstart.X(), planeID);
1261 auto const eTick = detProp.ConvertXToTicks(mcend.X(), planeID);
1262 mf::LogVerbatim(
"TrackAna") <<
" Wire:Tick in Pln " << ipl <<
" W:T " << sWire <<
":" 1263 << sTick <<
" - " << eWire <<
":" << eTick;
double E(const int i=0) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
simb::Origin_t Origin() const
double EndMomentum() const
double VertexMomentum() const
const simb::MCParticle * TrackIdToParticle_P(int id) const
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
The data type to uniquely identify a Plane.
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
constexpr auto abs(T v)
Returns the absolute value of the argument.
Vector_t VertexDirection() const
Access to track direction at different points.
const SMatrixSym55 & VertexCovariance() const
Access to covariance matrices.
std::string fStitchModuleLabel
WireID_t Wire
Index of the wire within its plane.
std::map< int, MCHists > fMCHistMap
3-dimensional objects, potentially hits, clusters, prongs, etc.
Rotation_t GlobalToLocalRotationAtPoint(size_t p) const
bool isValid() const noexcept
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
double HitCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
void anaStitch(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp)
Point_t const & Vertex() const
Access to track position at different points.
const SMatrixSym55 & EndCovariance() const
Access to covariance matrices.
std::string fMCTrackModuleLabel
EventNumber_t event() const
const TLorentzVector & Momentum() const
simb::Origin_t fOriginValue
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Vector_t EndDirection() const
Access to track direction at different points.
std::map< int, RecoHists > fRecoHistMap
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
Point_t const & End() const
Access to track position at different points.
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
const MCStep & Start() const
std::string fTrackModuleLabel
unsigned int TrackID() const
second_as<> second
Type of time stored in seconds, in double precision.
std::string fHitModuleLabel
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
cet::coded_exception< error, detail::translate > exception