853 std::unique_ptr<mf::LogInfo> pdump;
856 pdump = std::unique_ptr<mf::LogInfo>(
new mf::LogInfo(
"TrackAna"));
868 std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;
872 selected_mctracks.reserve(mctrackh->size());
877 *pdump <<
"MC Tracks\n" 878 <<
" Id pdg x y z dx dy dz p\n" 879 <<
"-------------------------------------------------------------------------------------------\n";
887 imctrk != mctrackh->end(); ++imctrk) {
895 int apdg = std::abs(pdg);
910 double mctime = mctrk.
Start().
T();
919 double plen = length(mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
928 selected_mctracks.push_back(std::make_pair(&mctrk, -1));
933 double pstart = mcstartmom.Mag();
934 double pend = mcendmom.Mag();
936 << std::setw(3) << mctrk.
TrackID()
937 << std::setw(6) << mctrk.
PdgCode()
939 << std::fixed << std::setprecision(2)
940 << std::setw(10) << mcdx
942 << std::setw(3) << mctrk.
TrackID()
943 << std::setw(6) << mctrk.
PdgCode()
945 << std::fixed << std::setprecision(2)
946 << std::setw(10) << mcstart[0]
947 << std::setw(10) << mcstart[1]
948 << std::setw(10) << mcstart[2];
951 << std::fixed << std::setprecision(3)
952 << std::setw(10) << mcstartmom[0]/pstart
953 << std::setw(10) << mcstartmom[1]/pstart
954 << std::setw(10) << mcstartmom[2]/pstart;
957 *pdump << std::setw(32) <<
" ";
958 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
960 << std::setw(3) << mctrk.
TrackID()
961 << std::setw(6) << mctrk.
PdgCode()
963 << std::fixed << std::setprecision(2)
964 << std::setw(10) << mcend[0]
965 << std::setw(10) << mcend[1]
966 << std::setw(10) << mcend[2];
969 << std::fixed << std::setprecision(3)
970 << std::setw(10) << mcendmom[0]/pend
971 << std::setw(10) << mcendmom[1]/pend
972 << std::setw(10) << mcendmom[2]/pend;
975 *pdump << std::setw(32)<<
" ";
976 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
977 <<
"\nLength: " << plen <<
"\n";
983 std::ostringstream ostr;
984 ostr <<
"MC" << (
fIgnoreSign ?
"All" : (pdg > 0 ?
"Pos" :
"Neg")) << std::abs(pdg);
989 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
990 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
991 double mcmom = mcstartmom.Mag();
993 double mcke = mcmom*mcmom / (std::sqrt(mcmom*mcmom + mcmass*mcmass) + mcmass);
995 mchists.fHmcstartx->Fill(mcstart.X());
996 mchists.fHmcstarty->Fill(mcstart.Y());
997 mchists.fHmcstartz->Fill(mcstart.Z());
998 mchists.fHmcendx->Fill(mcend.X());
999 mchists.fHmcendy->Fill(mcend.Y());
1000 mchists.fHmcendz->Fill(mcend.Z());
1001 mchists.fHmctheta->Fill(mcstartmom.Theta());
1002 mchists.fHmcphi->Fill(mcstartmom.Phi());
1003 mchists.fHmctheta_xz->Fill(mctheta_xz);
1004 mchists.fHmctheta_yz->Fill(mctheta_yz);
1005 mchists.fHmcmom->Fill(mcmom);
1006 mchists.fHmcmoml->Fill(mcmom);
1007 mchists.fHmcke->Fill(mcke);
1008 mchists.fHmckel->Fill(mcke);
1009 mchists.fHmclen->Fill(plen);
1010 mchists.fHmclens->Fill(plen);
1029 std::vector<art::Ptr<recob::Hit> > allhits;
1031 allhits.reserve(hith->size());
1032 for(
unsigned int i=0; i<hith->size(); ++i) {
1033 allhits.emplace_back(hith, i);
1048 <<
"TrackAna read " << trackvh->size()
1049 <<
" vectors of Stitched PtrVectorsof tracks.";
1056 *pdump <<
"\nReconstructed Tracks\n" 1057 <<
" Id MCid x y z dx dy dz p\n" 1058 <<
"-------------------------------------------------------------------------------------------\n";
1063 int ntrack = trackh->size();
1064 for(
int i = 0; i < ntrack; ++i) {
1070 std::vector<art::Ptr<recob::Hit> > trackhits;
1071 tkhit_find.get(i, trackhits);
1084 TVector3 pos = track.
Vertex<TVector3>();
1086 TVector3
end = track.
End<TVector3>();
1090 double dpos = bdist(pos);
1091 double dend = bdist(end);
1092 double tlen = length(track);
1093 double theta_xz = std::atan2(dir.X(), dir.Z());
1094 double theta_yz = std::atan2(dir.Y(), dir.Z());
1100 rhists.fHstartx->Fill(pos.X());
1101 rhists.fHstarty->Fill(pos.Y());
1102 rhists.fHstartz->Fill(pos.Z());
1103 rhists.fHstartd->Fill(dpos);
1104 rhists.fHendx->Fill(end.X());
1105 rhists.fHendy->Fill(end.Y());
1106 rhists.fHendz->Fill(end.Z());
1107 rhists.fHendd->Fill(dend);
1108 rhists.fHtheta->Fill(dir.Theta());
1109 rhists.fHphi->Fill(dir.Phi());
1110 rhists.fHtheta_xz->Fill(theta_xz);
1111 rhists.fHtheta_yz->Fill(theta_yz);
1116 rhists.fHmom->Fill(mom);
1117 rhists.fHmoml->Fill(mom);
1118 rhists.fHlen->Fill(tlen);
1119 rhists.fHlens->Fill(tlen);
1137 int start_point = (
swap == 0 ? 0 : ntraj-1);
1143 rot(1, 0) = -rot(1, 0);
1144 rot(2, 0) = -rot(2, 0);
1145 rot(1, 1) = -rot(1, 1);
1146 rot(2, 1) = -rot(2, 1);
1147 rot(1, 2) = -rot(1, 2);
1148 rot(2, 2) = -rot(2, 2);
1150 pos = track.
End<TVector3>();
1152 end = track.
Vertex<TVector3>();
1158 theta_xz = std::atan2(dir.X(), dir.Z());
1159 theta_yz = std::atan2(dir.Y(), dir.Z());
1170 for(
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1171 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1176 throw cet::exception(
"TrackAna") <<
"no particle with ID=" << pdg <<
"\n";
1177 const MCHists& mchists = iMCHistMap->second;
1181 double mctime = mctrk.
Start().
T();
1189 TVector3 mcstartmom;
1191 double plen = length(mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1195 TVector3 mcpos = mcstart - pos;
1200 TVector3 mcmoml = rot * mcstartmom;
1201 TVector3 mcposl = rot * mcpos;
1203 double colinearity = mcmoml.Z() / mcmoml.Mag();
1205 double u = mcposl.X();
1206 double v = mcposl.Y();
1207 double w = mcposl.Z();
1209 double pu = mcmoml.X();
1210 double pv = mcmoml.Y();
1211 double pw = mcmoml.Z();
1213 double dudw = pu / pw;
1214 double dvdw = pv / pw;
1216 double u0 = u - w * dudw;
1217 double v0 = v - w * dvdw;
1218 double uv0 = std::sqrt(u0*u0 + v0*v0);
1220 mchists.fHduvcosth->Fill(colinearity, uv0);
1226 mchists.fHmcdudw->Fill(dudw);
1227 mchists.fHmcdvdw->Fill(dvdw);
1228 mchists.fHdudwpull->Fill(dudw / std::sqrt(cov(2,2)));
1229 mchists.fHdvdwpull->Fill(dvdw / std::sqrt(cov(3,3)));
1231 mchists.fHcosth->Fill(colinearity);
1236 mchists.fHmcu->Fill(u0);
1237 mchists.fHmcv->Fill(v0);
1238 mchists.fHmcw->Fill(w);
1239 mchists.fHupull->Fill(u0 / std::sqrt(cov(0,0)));
1240 mchists.fHvpull->Fill(v0 / std::sqrt(cov(1,1)));
1246 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
1247 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
1248 double mcmom = mcstartmom.Mag();
1250 double mcke = mcmom*mcmom / (std::sqrt(mcmom*mcmom + mcmass*mcmass) + mcmass);
1252 mchists.fHstartdx->Fill(pos.X() - mcstart.X());
1253 mchists.fHstartdy->Fill(pos.Y() - mcstart.Y());
1254 mchists.fHstartdz->Fill(pos.Z() - mcstart.Z());
1255 mchists.fHenddx->Fill(end.X() - mcend.X());
1256 mchists.fHenddy->Fill(end.Y() - mcend.Y());
1257 mchists.fHenddz->Fill(end.Z() - mcend.Z());
1258 mchists.fHlvsl->Fill(plen, tlen);
1259 mchists.fHdl->Fill(tlen - plen);
1260 mchists.fHpvsp->Fill(mcmom, mom);
1261 double dp = mom - mcmom;
1262 mchists.fHdp->Fill(dp);
1263 mchists.fHppull->Fill(dp / std::sqrt(cov(4,4)));
1264 if(std::abs(dpos) >= 5. && std::abs(dend) >= 5.) {
1265 mchists.fHpvspc->Fill(mcmom, mom);
1266 mchists.fHdpc->Fill(dp);
1267 mchists.fHppullc->Fill(dp / std::sqrt(cov(4,4)));
1282 std::set<int> tkidset;
1283 tkidset.insert(mcid);
1285 bt_serv->HitCollectionEfficiency(tkidset, trackhits, allhits,
geo::k3D);
1286 double hitpurity = bt_serv->HitCollectionPurity(tkidset, trackhits);
1287 mchists.fHHitEff->Fill(hiteff);
1288 mchists.fHHitPurity->Fill(hitpurity);
1292 mchists.fHgstartx->Fill(mcstart.X());
1293 mchists.fHgstarty->Fill(mcstart.Y());
1294 mchists.fHgstartz->Fill(mcstart.Z());
1295 mchists.fHgendx->Fill(mcend.X());
1296 mchists.fHgendy->Fill(mcend.Y());
1297 mchists.fHgendz->Fill(mcend.Z());
1298 mchists.fHgtheta->Fill(mcstartmom.Theta());
1299 mchists.fHgphi->Fill(mcstartmom.Phi());
1300 mchists.fHgtheta_xz->Fill(mctheta_xz);
1301 mchists.fHgtheta_yz->Fill(mctheta_yz);
1302 mchists.fHgmom->Fill(mcmom);
1303 mchists.fHgmoml->Fill(mcmom);
1304 mchists.fHgke->Fill(mcke);
1305 mchists.fHgkel->Fill(mcke);
1306 mchists.fHglen->Fill(plen);
1307 mchists.fHglens->Fill(plen);
1310 selected_mctracks[imc].second = i;
1314 float KE = ptkl->
E() - ptkl->
Mass();
1315 std::string KEUnits =
" GeV";
1323 <<
" Match MCTkID "<<std::setw(6)<<mctrk.
TrackID()
1324 <<
" Origin "<<mctrk.
Origin()
1325 <<
" PDG"<<std::setw(5)<<mctrk.
PdgCode()
1326 <<
" KE"<<std::setw(4)<<(int)KE<<KEUnits
1327 <<
" RecoTrkID "<<track.
ID()
1328 <<
" hitEff "<<std::setprecision(2)<<hiteff<<
" hitPur "<<hitpurity;
1329 int sWire, sTick, eWire, eTick;
1331 for(
unsigned short ipl = 0; ipl < geom->
Nplanes(); ++ipl) {
1336 mf::LogVerbatim(
"TrackAna")<<
" Wire:Tick in Pln "<<ipl<<
" W:T "<<sWire<<
":"<<sTick<<
" - "<<eWire<<
":"<<eTick;
1348 TVector3 pos = track.
Vertex<TVector3>();
1350 TVector3 end = track.
End<TVector3>();
1356 *pdump <<
"\nOffset" 1357 << std::setw(3) << track.
ID()
1358 << std::setw(6) << mcid
1360 << std::fixed << std::setprecision(2)
1361 << std::setw(10) << trackdx
1363 << std::setw(3) << track.
ID()
1364 << std::setw(6) << mcid
1366 << std::fixed << std::setprecision(2)
1367 << std::setw(10) << pos[0]
1368 << std::setw(10) << pos[1]
1369 << std::setw(10) << pos[2];
1372 << std::fixed << std::setprecision(3)
1373 << std::setw(10) << dir[0]
1374 << std::setw(10) << dir[1]
1375 << std::setw(10) << dir[2];
1378 *pdump << std::setw(32) <<
" ";
1379 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
1381 << std::setw(3) << track.
ID()
1382 << std::setw(6) << mcid
1384 << std::fixed << std::setprecision(2)
1385 << std::setw(10) << end[0]
1386 << std::setw(10) << end[1]
1387 << std::setw(10) << end[2];
1390 << std::fixed << std::setprecision(3)
1391 << std::setw(10) << enddir[0]
1392 << std::setw(10) << enddir[1]
1393 << std::setw(10) << enddir[2];
1396 *pdump << std::setw(32)<<
" ";
1397 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
1398 <<
"\nLength: " << tlen <<
"\n";
1406 for(
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1407 if(selected_mctracks[imc].second >= 0)
continue;
1408 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1410 float KE = ptkl->
E() - ptkl->
Mass();
1411 std::string KEUnits =
" GeV";
1418 TVector3 mcstart, mcend, mcstartmom, mcendmom;
1420 double plen = length(mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1422 <<
" NoMat MCTkID "<<std::setw(6)<<mctrk.
TrackID()
1423 <<
" Origin "<<mctrk.
Origin()
1424 <<
" PDG"<<std::setw(5)<<mctrk.
PdgCode()
1425 <<
" KE"<<std::setw(4)<<(int)KE<<KEUnits
1426 <<
" Length "<<std::fixed<<std::setprecision(1)<<plen<<
" cm";
1428 int sWire, sTick, eWire, eTick;
1430 for(
unsigned short ipl = 0; ipl < geom->
Nplanes(); ++ipl) {
1436 <<
" W:T "<<sWire<<
":"<<sTick<<
" - "<<eWire<<
":"<<eTick;
double E(const int i=0) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
const simb::MCParticle * TrackIdToParticle_P(int const &id)
simb::Origin_t Origin() const
double EndMomentum() const
double VertexMomentum() const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Vector_t VertexDirection() const
Access to track direction at different points.
const SMatrixSym55 & VertexCovariance() const
Access to covariance matrices.
std::string fStitchModuleLabel
void anaStitch(const art::Event &evt)
std::map< int, MCHists > fMCHistMap
virtual double ConvertXToTicks(double X, int p, int t, int c) const =0
3-dimensional objects, potentially hits, clusters, prongs, etc.
Rotation_t GlobalToLocalRotationAtPoint(size_t p) const
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
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.
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
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.
const MCStep & Start() const
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
std::string fTrackModuleLabel
unsigned int TrackID() const
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