30 #include "art_root_io/TFileService.h" 33 #include "cetlib_except/exception.h" 57 double bdist(
const TVector3& pos)
64 double d2 = 2. * tpc.HalfWidth() - pos.X();
65 double d3 = pos.Y() + tpc.HalfHeight();
66 double d4 = tpc.HalfHeight() - pos.Y();
68 double d6 = tpc.Length() - pos.Z();
70 return std::min({d1, d2, d3, d4, d5, d6});
94 double xmax = 2. * tpc.HalfWidth();
95 double ymin = -tpc.HalfHeight();
96 double ymax = tpc.HalfHeight();
98 double zmax = tpc.Length();
104 for (
int i = 0; i <
n; ++i) {
105 TVector3 pos = part.
Position(i).Vect();
111 if (pos.X() >= xmin && pos.X() <= xmax && pos.Y() >= ymin && pos.Y() <= ymax &&
112 pos.Z() >= zmin && pos.Z() <= zmax) {
122 result += disp.Mag();
152 double xmax = 2. * tpc.HalfWidth();
153 double ymin = -tpc.HalfHeight();
154 double ymax = tpc.HalfHeight();
156 double zmax = tpc.Length();
159 int n = mctrk.size();
162 for (
int i = 0; i <
n; ++i) {
163 TVector3 pos = mctrk[i].Position().Vect();
169 if (pos.X() >= xmin && pos.X() <= xmax && pos.Y() >= ymin && pos.Y() <= ymax &&
170 pos.Z() >= zmin && pos.Z() <= zmax) {
176 startmom = 0.001 * mctrk[i].Momentum().Vect();
180 result += disp.Mag();
185 endmom = 0.001 * mctrk[i].Momentum().Vect();
195 void effcalc(
const TH1* hnum,
const TH1* hden, TH1* heff)
197 int nbins = hnum->GetNbinsX();
198 if (nbins != hden->GetNbinsX())
199 throw cet::exception(
"TrackAna") <<
"effcalc[" __FILE__
"]: incompatible histograms (I)\n";
200 if (nbins != heff->GetNbinsX())
201 throw cet::exception(
"TrackAna") <<
"effcalc[" __FILE__
"]: incompatible histograms (II)\n";
205 for (
int ibin = 0; ibin <= nbins + 1; ++ibin) {
206 double num = hnum->GetBinContent(ibin);
207 double den = hden->GetBinContent(ibin);
209 heff->SetBinContent(ibin, 0.);
210 heff->SetBinError(ibin, 0.);
213 double eff = num /
den;
214 if (eff < 0.) eff = 0.;
215 if (eff > 1.) eff = 1.;
216 double err = std::sqrt(eff * (1. - eff) / den);
217 heff->SetBinContent(ibin, eff);
218 heff->SetBinError(ibin, err);
222 heff->SetMinimum(0.);
223 heff->SetMaximum(1.05);
224 heff->SetMarkerStyle(20);
227 class flattener :
public std::vector<unsigned int> {
229 flattener() =
default;
231 flattener(
const std::vector<std::vector<unsigned int>>& input) { Convert(input); }
233 void Convert(
const std::vector<std::vector<unsigned int>>& input)
237 for (
auto const& vec : input)
238 length += vec.size();
241 for (
auto const& vec : input)
242 for (
auto const&
value : vec)
258 explicit RecoHists(
const std::string& subdir);
262 TH1F* fHstartx{
nullptr};
263 TH1F* fHstarty{
nullptr};
264 TH1F* fHstartz{
nullptr};
265 TH1F* fHstartd{
nullptr};
266 TH1F* fHendx{
nullptr};
267 TH1F* fHendy{
nullptr};
268 TH1F* fHendz{
nullptr};
269 TH1F* fHendd{
nullptr};
270 TH1F* fHtheta{
nullptr};
271 TH1F* fHphi{
nullptr};
272 TH1F* fHtheta_xz{
nullptr};
273 TH1F* fHtheta_yz{
nullptr};
274 TH1F* fHmom{
nullptr};
275 TH1F* fHmoml{
nullptr};
276 TH1F* fHlen{
nullptr};
277 TH1F* fHlens{
nullptr};
281 TH1F* fHHitChg{
nullptr};
282 TH1F* fHHitWidth{
nullptr};
283 TH1F* fHHitPdg{
nullptr};
284 TH1F* fHHitTrkId{
nullptr};
285 TH1F* fModeFrac{
nullptr};
286 TH1F* fNTrkIdTrks{
nullptr};
287 TH2F* fNTrkIdTrks2{
nullptr};
288 TH2F* fNTrkIdTrks3{
nullptr};
294 explicit MCHists(
const std::string& subdir);
298 TH2F* fHduvcosth{
nullptr};
299 TH1F* fHcosth{
nullptr};
300 TH1F* fHmcu{
nullptr};
301 TH1F* fHmcv{
nullptr};
302 TH1F* fHmcw{
nullptr};
303 TH1F* fHupull{
nullptr};
304 TH1F* fHvpull{
nullptr};
305 TH1F* fHmcdudw{
nullptr};
306 TH1F* fHmcdvdw{
nullptr};
307 TH1F* fHdudwpull{
nullptr};
308 TH1F* fHdvdwpull{
nullptr};
309 TH1F* fHHitEff{
nullptr};
310 TH1F* fHHitPurity{
nullptr};
314 TH1F* fHstartdx{
nullptr};
315 TH1F* fHstartdy{
nullptr};
316 TH1F* fHstartdz{
nullptr};
317 TH1F* fHenddx{
nullptr};
318 TH1F* fHenddy{
nullptr};
319 TH1F* fHenddz{
nullptr};
320 TH2F* fHlvsl{
nullptr};
322 TH2F* fHpvsp{
nullptr};
323 TH2F* fHpvspc{
nullptr};
325 TH1F* fHdpc{
nullptr};
326 TH1F* fHppull{
nullptr};
327 TH1F* fHppullc{
nullptr};
331 TH1F* fHmcstartx{
nullptr};
332 TH1F* fHmcstarty{
nullptr};
333 TH1F* fHmcstartz{
nullptr};
334 TH1F* fHmcendx{
nullptr};
335 TH1F* fHmcendy{
nullptr};
336 TH1F* fHmcendz{
nullptr};
337 TH1F* fHmctheta{
nullptr};
338 TH1F* fHmcphi{
nullptr};
339 TH1F* fHmctheta_xz{
nullptr};
340 TH1F* fHmctheta_yz{
nullptr};
341 TH1F* fHmcmom{
nullptr};
342 TH1F* fHmcmoml{
nullptr};
343 TH1F* fHmcke{
nullptr};
344 TH1F* fHmckel{
nullptr};
345 TH1F* fHmclen{
nullptr};
346 TH1F* fHmclens{
nullptr};
350 TH1F* fHgstartx{
nullptr};
351 TH1F* fHgstarty{
nullptr};
352 TH1F* fHgstartz{
nullptr};
353 TH1F* fHgendx{
nullptr};
354 TH1F* fHgendy{
nullptr};
355 TH1F* fHgendz{
nullptr};
356 TH1F* fHgtheta{
nullptr};
357 TH1F* fHgphi{
nullptr};
358 TH1F* fHgtheta_xz{
nullptr};
359 TH1F* fHgtheta_yz{
nullptr};
360 TH1F* fHgmom{
nullptr};
361 TH1F* fHgmoml{
nullptr};
362 TH1F* fHgke{
nullptr};
363 TH1F* fHgkel{
nullptr};
364 TH1F* fHglen{
nullptr};
365 TH1F* fHglens{
nullptr};
369 TH1F* fHestartx{
nullptr};
370 TH1F* fHestarty{
nullptr};
371 TH1F* fHestartz{
nullptr};
372 TH1F* fHeendx{
nullptr};
373 TH1F* fHeendy{
nullptr};
374 TH1F* fHeendz{
nullptr};
375 TH1F* fHetheta{
nullptr};
376 TH1F* fHephi{
nullptr};
377 TH1F* fHetheta_xz{
nullptr};
378 TH1F* fHetheta_yz{
nullptr};
379 TH1F* fHemom{
nullptr};
380 TH1F* fHemoml{
nullptr};
381 TH1F* fHeke{
nullptr};
382 TH1F* fHekel{
nullptr};
383 TH1F* fHelen{
nullptr};
384 TH1F* fHelens{
nullptr};
393 void endJob()
override;
399 std::vector<size_t> fsort_indexes(
const std::vector<double>& v);
449 art::TFileDirectory topdir = tfs->mkdir(
"trkana",
"TrackAna histograms");
450 art::TFileDirectory
dir = topdir.mkdir(subdir);
454 fHstartx = dir.make<TH1F>(
455 "xstart",
"X Start Position", 100, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
457 dir.make<TH1F>(
"ystart",
"Y Start Position", 100, -tpc.HalfHeight(), tpc.HalfHeight());
458 fHstartz = dir.make<TH1F>(
"zstart",
"Z Start Position", 100, 0., tpc.Length());
460 dir.make<TH1F>(
"dstart",
"Start Position Distance to Boundary", 100, -10., tpc.HalfWidth());
462 dir.make<TH1F>(
"xend",
"X End Position", 100, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
463 fHendy = dir.make<TH1F>(
"yend",
"Y End Position", 100, -tpc.HalfHeight(), tpc.HalfHeight());
464 fHendz = dir.make<TH1F>(
"zend",
"Z End Position", 100, 0., tpc.Length());
466 dir.make<TH1F>(
"dend",
"End Position Distance to Boundary", 100, -10., tpc.HalfWidth());
467 fHtheta = dir.make<TH1F>(
"theta",
"Theta", 100, 0., 3.142);
468 fHphi = dir.make<TH1F>(
"phi",
"Phi", 100, -3.142, 3.142);
469 fHtheta_xz = dir.make<TH1F>(
"theta_xz",
"Theta_xz", 100, -3.142, 3.142);
470 fHtheta_yz = dir.make<TH1F>(
"theta_yz",
"Theta_yz", 100, -3.142, 3.142);
471 fHmom = dir.make<TH1F>(
"mom",
"Momentum", 100, 0., 10.);
472 fHmoml = dir.make<TH1F>(
"moml",
"Momentum", 100, 0., 1.);
473 fHlen = dir.make<TH1F>(
"len",
"Track Length", 100, 0., 1.1 * tpc.Length());
474 fHlens = dir.make<TH1F>(
"lens",
"Track Length", 100, 0., 0.1 * tpc.Length());
475 fHHitChg = dir.make<TH1F>(
"hchg",
"Hit Charge (ADC counts)", 100, 0., 4000.);
476 fHHitWidth = dir.make<TH1F>(
"hwid",
"Hit Width (ticks)", 40, 0., 20.);
477 fHHitPdg = dir.make<TH1F>(
"hpdg",
"Hit Pdg code", 5001, -2500.5, +2500.5);
478 fHHitTrkId = dir.make<TH1F>(
"htrkid",
"Hit Track ID", 401, -200.5, +200.5);
480 dir.make<TH1F>(
"hmodefrac",
481 "quasi-Purity: Fraction of component tracks with the Track mode value",
486 dir.make<TH1F>(
"hntrkids",
487 "quasi-Efficiency: Number of stitched tracks in which TrkId appears",
491 fNTrkIdTrks2 = dir.make<TH2F>(
"hntrkids2",
492 "Number of stitched tracks in which TrkId appears vs KE [GeV]",
499 fNTrkIdTrks3 = dir.make<TH2F>(
"hntrkids3",
500 "MC Track vs Reco Track, wtd by nhits on Collection Plane",
507 fNTrkIdTrks3->GetXaxis()->SetTitle(
"Sorted-by-Descending-CPlane-Hits outer Track Number");
508 fNTrkIdTrks3->GetYaxis()->SetTitle(
"Sorted-by-Descending-True-Length G4Track");
513 TrackAna::MCHists::MCHists(
const std::string& subdir)
522 art::TFileDirectory topdir = tfs->mkdir(
"trkana",
"TrackAna histograms");
523 art::TFileDirectory
dir = topdir.mkdir(subdir);
528 dir.make<TH2F>(
"duvcosth",
"Delta(uv) vs. Colinearity", 100, 0.95, 1., 100, 0., 1.);
529 fHcosth = dir.make<TH1F>(
"colin",
"Colinearity", 100, 0.95, 1.);
530 fHmcu = dir.make<TH1F>(
"mcu",
"MC Truth U", 100, -5., 5.);
531 fHmcv = dir.make<TH1F>(
"mcv",
"MC Truth V", 100, -5., 5.);
532 fHmcw = dir.make<TH1F>(
"mcw",
"MC Truth W", 100, -20., 20.);
533 fHupull = dir.make<TH1F>(
"dupull",
"U Pull", 100, -20., 20.);
534 fHvpull = dir.make<TH1F>(
"dvpull",
"V Pull", 100, -20., 20.);
535 fHmcdudw = dir.make<TH1F>(
"mcdudw",
"MC Truth U Slope", 100, -0.2, 0.2);
536 fHmcdvdw = dir.make<TH1F>(
"mcdvdw",
"MV Truth V Slope", 100, -0.2, 0.2);
537 fHdudwpull = dir.make<TH1F>(
"dudwpull",
"U Slope Pull", 100, -10., 10.);
538 fHdvdwpull = dir.make<TH1F>(
"dvdwpull",
"V Slope Pull", 100, -10., 10.);
539 fHHitEff = dir.make<TH1F>(
"hiteff",
"MC Hit Efficiency", 100, 0., 1.0001);
540 fHHitPurity = dir.make<TH1F>(
"hitpurity",
"MC Hit Purity", 100, 0., 1.0001);
541 fHstartdx = dir.make<TH1F>(
"startdx",
"Start Delta x", 100, -10., 10.);
542 fHstartdy = dir.make<TH1F>(
"startdy",
"Start Delta y", 100, -10., 10.);
543 fHstartdz = dir.make<TH1F>(
"startdz",
"Start Delta z", 100, -10., 10.);
544 fHenddx = dir.make<TH1F>(
"enddx",
"End Delta x", 100, -10., 10.);
545 fHenddy = dir.make<TH1F>(
"enddy",
"End Delta y", 100, -10., 10.);
546 fHenddz = dir.make<TH1F>(
"enddz",
"End Delta z", 100, -10., 10.);
547 fHlvsl = dir.make<TH2F>(
"lvsl",
548 "Reco Length vs. MC Truth Length",
555 fHdl = dir.make<TH1F>(
"dl",
"Track Length Minus MC Particle Length", 100, -50., 50.);
557 dir.make<TH2F>(
"pvsp",
"Reco Momentum vs. MC Truth Momentum", 100, 0., 5., 100, 0., 5.);
558 fHpvspc = dir.make<TH2F>(
559 "pvspc",
"Reco Momentum vs. MC Truth Momentum (Contained Tracks)", 100, 0., 5., 100, 0., 5.);
560 fHdp = dir.make<TH1F>(
"dp",
"Reco-MC Momentum Difference", 100, -5., 5.);
561 fHdpc = dir.make<TH1F>(
"dpc",
"Reco-MC Momentum Difference (Contained Tracks)", 100, -5., 5.);
562 fHppull = dir.make<TH1F>(
"ppull",
"Momentum Pull", 100, -10., 10.);
563 fHppullc = dir.make<TH1F>(
"ppullc",
"Momentum Pull (Contained Tracks)", 100, -10., 10.);
565 fHmcstartx = dir.make<TH1F>(
566 "mcxstart",
"MC X Start Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
568 dir.make<TH1F>(
"mcystart",
"MC Y Start Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
569 fHmcstartz = dir.make<TH1F>(
"mczstart",
"MC Z Start Position", 10, 0., tpc.Length());
570 fHmcendx = dir.make<TH1F>(
571 "mcxend",
"MC X End Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
573 dir.make<TH1F>(
"mcyend",
"MC Y End Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
574 fHmcendz = dir.make<TH1F>(
"mczend",
"MC Z End Position", 10, 0., tpc.Length());
575 fHmctheta = dir.make<TH1F>(
"mctheta",
"MC Theta", 20, 0., 3.142);
576 fHmcphi = dir.make<TH1F>(
"mcphi",
"MC Phi", 10, -3.142, 3.142);
577 fHmctheta_xz = dir.make<TH1F>(
"mctheta_xz",
"MC Theta_xz", 40, -3.142, 3.142);
578 fHmctheta_yz = dir.make<TH1F>(
"mctheta_yz",
"MC Theta_yz", 40, -3.142, 3.142);
579 fHmcmom = dir.make<TH1F>(
"mcmom",
"MC Momentum", 10, 0., 10.);
580 fHmcmoml = dir.make<TH1F>(
"mcmoml",
"MC Momentum", 10, 0., 1.);
581 fHmcke = dir.make<TH1F>(
"mcke",
"MC Kinetic Energy", 10, 0., 10.);
582 fHmckel = dir.make<TH1F>(
"mckel",
"MC Kinetic Energy", 10, 0., 1.);
583 fHmclen = dir.make<TH1F>(
"mclen",
"MC Particle Length", 10, 0., 1.1 * tpc.Length());
584 fHmclens = dir.make<TH1F>(
"mclens",
"MC Particle Length", 10, 0., 0.1 * tpc.Length());
586 fHgstartx = dir.make<TH1F>(
587 "gxstart",
"Good X Start Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
589 dir.make<TH1F>(
"gystart",
"Good Y Start Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
590 fHgstartz = dir.make<TH1F>(
"gzstart",
"Good Z Start Position", 10, 0., tpc.Length());
591 fHgendx = dir.make<TH1F>(
592 "gxend",
"Good X End Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
594 dir.make<TH1F>(
"gyend",
"Good Y End Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
595 fHgendz = dir.make<TH1F>(
"gzend",
"Good Z End Position", 10, 0., tpc.Length());
596 fHgtheta = dir.make<TH1F>(
"gtheta",
"Good Theta", 20, 0., 3.142);
597 fHgphi = dir.make<TH1F>(
"gphi",
"Good Phi", 10, -3.142, 3.142);
598 fHgtheta_xz = dir.make<TH1F>(
"gtheta_xz",
"Good Theta_xz", 40, -3.142, 3.142);
599 fHgtheta_yz = dir.make<TH1F>(
"gtheta_yz",
"Good Theta_yz", 40, -3.142, 3.142);
600 fHgmom = dir.make<TH1F>(
"gmom",
"Good Momentum", 10, 0., 10.);
601 fHgmoml = dir.make<TH1F>(
"gmoml",
"Good Momentum", 10, 0., 1.);
602 fHgke = dir.make<TH1F>(
"gke",
"Good Kinetic Energy", 10, 0., 10.);
603 fHgkel = dir.make<TH1F>(
"gkel",
"Good Kinetic Energy", 10, 0., 1.);
604 fHglen = dir.make<TH1F>(
"glen",
"Good Particle Length", 10, 0., 1.1 * tpc.Length());
605 fHglens = dir.make<TH1F>(
"glens",
"Good Particle Length", 10, 0., 0.1 * tpc.Length());
607 fHestartx = dir.make<TH1F>(
"exstart",
608 "Efficiency vs. X Start Position",
610 -2. * tpc.HalfWidth(),
611 4. * tpc.HalfWidth());
612 fHestarty = dir.make<TH1F>(
613 "eystart",
"Efficiency vs. Y Start Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
614 fHestartz = dir.make<TH1F>(
"ezstart",
"Efficiency vs. Z Start Position", 10, 0., tpc.Length());
615 fHeendx = dir.make<TH1F>(
616 "exend",
"Efficiency vs. X End Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
617 fHeendy = dir.make<TH1F>(
618 "eyend",
"Efficiency vs. Y End Position", 10, -tpc.HalfHeight(), tpc.HalfHeight());
619 fHeendz = dir.make<TH1F>(
"ezend",
"Efficiency vs. Z End Position", 10, 0., tpc.Length());
620 fHetheta = dir.make<TH1F>(
"etheta",
"Efficiency vs. Theta", 20, 0., 3.142);
621 fHephi = dir.make<TH1F>(
"ephi",
"Efficiency vs. Phi", 10, -3.142, 3.142);
622 fHetheta_xz = dir.make<TH1F>(
"etheta_xz",
"Efficiency vs. Theta_xz", 40, -3.142, 3.142);
623 fHetheta_yz = dir.make<TH1F>(
"etheta_yz",
"Efficiency vs. Theta_yz", 40, -3.142, 3.142);
624 fHemom = dir.make<TH1F>(
"emom",
"Efficiency vs. Momentum", 10, 0., 10.);
625 fHemoml = dir.make<TH1F>(
"emoml",
"Efficiency vs. Momentum", 10, 0., 1.);
626 fHeke = dir.make<TH1F>(
"eke",
"Efficiency vs. Kinetic Energy", 10, 0., 10.);
627 fHekel = dir.make<TH1F>(
"ekel",
"Efficiency vs. Kinetic Energy", 10, 0., 1.);
628 fHelen = dir.make<TH1F>(
"elen",
"Efficiency vs. Particle Length", 10, 0., 1.1 * tpc.Length());
629 fHelens = dir.make<TH1F>(
"elens",
"Efficiency vs. Particle Length", 10, 0., 0.1 * tpc.Length());
634 , fTrackModuleLabel(pset.
get<
std::string>(
"TrackModuleLabel"))
635 , fMCTrackModuleLabel(pset.
get<
std::string>(
"MCTrackModuleLabel"))
636 , fSpacepointModuleLabel(pset.
get<
std::string>(
"SpacepointModuleLabel"))
637 , fStitchModuleLabel(pset.
get<
std::string>(
"StitchModuleLabel"))
638 , fTrkSpptAssocModuleLabel(pset.
get<
std::string>(
"TrkSpptAssocModuleLabel"))
639 , fHitSpptAssocModuleLabel(pset.
get<
std::string>(
"HitSpptAssocModuleLabel"))
640 , fHitModuleLabel(pset.
get<
std::string>(
"HitModuleLabel"))
641 , fDump(pset.
get<int>(
"Dump"))
642 , fMinMCKE(pset.
get<double>(
"MinMCKE"))
643 , fMinMCLen(pset.
get<double>(
"MinMCLen"))
644 , fMatchColinearity(pset.
get<double>(
"MatchColinearity"))
645 , fMatchDisp(pset.
get<double>(
"MatchDisp"))
646 , fWMatchDisp(pset.
get<double>(
"WMatchDisp"))
647 , fMatchLength(pset.
get<double>(
"MatchLength"))
648 , fIgnoreSign(pset.
get<bool>(
"IgnoreSign"))
649 , fStitchedAnalysis(pset.
get<bool>(
"StitchedAnalysis", false))
650 , fOrigin(pset.
get<
std::string>(
"MCTrackOrigin",
"Any"))
651 , fPrintLevel(pset.
get<int>(
"PrintLevel", 0))
658 if (
fOrigin.find(
"Beam") != std::string::npos) {
662 else if (
fOrigin.find(
"Cosmic") != std::string::npos) {
666 else if (
fOrigin.find(
"Super") != std::string::npos) {
670 else if (
fOrigin.find(
"Single") != std::string::npos) {
677 mf::LogInfo(
"TrackAna") <<
"TrackAna configured with the following parameters:\n" 684 <<
" Dump = " <<
fDump <<
"\n" 685 <<
" MinMCKE = " <<
fMinMCKE <<
"\n" 700 std::unique_ptr<mf::LogInfo> pdump;
703 pdump = std::make_unique<mf::LogInfo>(
"TrackAna");
715 std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;
721 if (mc && mctrackh.
isValid()) {
723 selected_mctracks.reserve(mctrackh->size());
728 *pdump <<
"MC Tracks\n" 729 <<
" Id pdg x y z dx dy dz " 731 <<
"--------------------------------------------------------------------------------" 740 imctrk != mctrackh->end();
763 double mctime = mctrk.
Start().
T();
772 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
781 selected_mctracks.push_back(std::make_pair(&mctrk, -1));
786 double pstart = mcstartmom.Mag();
787 double pend = mcendmom.Mag();
788 *pdump <<
"\nOffset" << std::setw(3) << mctrk.
TrackID() << std::setw(6)
789 << mctrk.
PdgCode() <<
" " << std::fixed << std::setprecision(2)
790 << std::setw(10) << mcdx <<
"\nStart " << std::setw(3) << mctrk.
TrackID()
791 << std::setw(6) << mctrk.
PdgCode() <<
" " << std::fixed
792 << std::setprecision(2) << std::setw(10) << mcstart.X() << std::setw(10)
793 << mcstart.Y() << std::setw(10) << mcstart.Z();
795 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
796 << mcstartmom[0] / pstart << std::setw(10) << mcstartmom[1] / pstart
797 << std::setw(10) << mcstartmom[2] / pstart;
800 *pdump << std::setw(32) <<
" ";
801 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
802 *pdump <<
"\nEnd " << std::setw(3) << mctrk.
TrackID() << std::setw(6)
803 << mctrk.
PdgCode() <<
" " << std::fixed << std::setprecision(2)
804 << std::setw(10) << mcend.X() << std::setw(10) << mcend.Y() << std::setw(10)
807 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10)
808 << mcendmom[0] / pend << std::setw(10) << mcendmom[1] / pend
809 << std::setw(10) << mcendmom[2] / pend;
812 *pdump << std::setw(32) <<
" ";
813 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
814 <<
"\nLength: " << plen <<
"\n";
820 std::ostringstream ostr;
826 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
827 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
828 double mcmom = mcstartmom.Mag();
830 double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
838 mchists.
fHmctheta->Fill(mcstartmom.Theta());
839 mchists.
fHmcphi->Fill(mcstartmom.Phi());
844 mchists.
fHmcke->Fill(mcke);
865 std::vector<art::Ptr<recob::Hit>> allhits;
867 allhits.reserve(hith->size());
868 for (
unsigned int i = 0; i < hith->size(); ++i) {
869 allhits.emplace_back(hith, i);
882 mf::LogDebug(
"TrackAna") <<
"TrackAna read " << trackvh->size()
883 <<
" vectors of Stitched PtrVectorsof tracks.";
890 *pdump <<
"\nReconstructed Tracks\n" 891 <<
" Id MCid x y z dx dy dz " 893 <<
"--------------------------------------------------------------------------------" 899 int ntrack = trackh->size();
900 for (
int i = 0; i < ntrack; ++i) {
906 std::vector<art::Ptr<recob::Hit>> trackhits;
907 tkhit_find.
get(i, trackhits);
917 TVector3 pos = track.
Vertex<TVector3>();
919 TVector3 end = track.
End<TVector3>();
923 double dpos = bdist(pos);
924 double dend = bdist(end);
925 double tlen = length(track);
926 double theta_xz = std::atan2(dir.X(), dir.Z());
927 double theta_yz = std::atan2(dir.Y(), dir.Z());
936 rhists.
fHendx->Fill(end.X());
937 rhists.
fHendy->Fill(end.Y());
938 rhists.
fHendz->Fill(end.Z());
939 rhists.
fHendd->Fill(dend);
940 rhists.
fHtheta->Fill(dir.Theta());
941 rhists.
fHphi->Fill(dir.Phi());
947 rhists.
fHmom->Fill(mom);
949 rhists.
fHlen->Fill(tlen);
950 rhists.
fHlens->Fill(tlen);
958 for (
int swap = 0; swap < 2; ++swap) {
968 int start_point = (swap == 0 ? 0 : ntraj - 1);
974 rot(1, 0) = -rot(1, 0);
975 rot(2, 0) = -rot(2, 0);
976 rot(1, 1) = -rot(1, 1);
977 rot(2, 1) = -rot(2, 1);
978 rot(1, 2) = -rot(1, 2);
979 rot(2, 2) = -rot(2, 2);
981 pos = track.
End<TVector3>();
983 end = track.
Vertex<TVector3>();
989 theta_xz = std::atan2(dir.X(), dir.Z());
990 theta_yz = std::atan2(dir.Y(), dir.Z());
1001 for (
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1002 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1007 throw cet::exception(
"TrackAna") <<
"no particle with ID=" << pdg <<
"\n";
1008 const MCHists& mchists = iMCHistMap->second;
1012 double mctime = mctrk.
Start().
T();
1020 TVector3 mcstartmom;
1022 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1026 TVector3 mcpos{mcstart.X() - pos[0], mcstart.Y() - pos[1], mcstart.Z() - pos[2]};
1031 TVector3 mcmoml = rot * mcstartmom;
1032 TVector3 mcposl = rot * mcpos;
1034 double colinearity = mcmoml.Z() / mcmoml.Mag();
1036 double u = mcposl.X();
1037 double v = mcposl.Y();
1038 double w = mcposl.Z();
1040 double pu = mcmoml.X();
1041 double pv = mcmoml.Y();
1042 double pw = mcmoml.Z();
1044 double dudw = pu / pw;
1045 double dvdw = pv / pw;
1047 double u0 = u - w * dudw;
1048 double v0 = v - w * dvdw;
1049 double uv0 = std::hypot(u0, v0);
1059 mchists.
fHdudwpull->Fill(dudw / std::sqrt(cov(2, 2)));
1060 mchists.
fHdvdwpull->Fill(dvdw / std::sqrt(cov(3, 3)));
1062 mchists.
fHcosth->Fill(colinearity);
1067 mchists.
fHmcu->Fill(u0);
1068 mchists.
fHmcv->Fill(v0);
1069 mchists.
fHmcw->Fill(w);
1070 mchists.
fHupull->Fill(u0 / std::sqrt(cov(0, 0)));
1071 mchists.
fHvpull->Fill(v0 / std::sqrt(cov(1, 1)));
1077 double mctheta_xz = std::atan2(mcstartmom.X(), mcstartmom.Z());
1078 double mctheta_yz = std::atan2(mcstartmom.Y(), mcstartmom.Z());
1079 double mcmom = mcstartmom.Mag();
1081 double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
1083 mchists.
fHstartdx->Fill(pos.X() - mcstart.X());
1084 mchists.
fHstartdy->Fill(pos.Y() - mcstart.Y());
1085 mchists.
fHstartdz->Fill(pos.Z() - mcstart.Z());
1086 mchists.
fHenddx->Fill(end.X() - mcend.X());
1087 mchists.
fHenddy->Fill(end.Y() - mcend.Y());
1088 mchists.
fHenddz->Fill(end.Z() - mcend.Z());
1089 mchists.
fHlvsl->Fill(plen, tlen);
1090 mchists.
fHdl->Fill(tlen - plen);
1091 mchists.
fHpvsp->Fill(mcmom, mom);
1092 double dp = mom - mcmom;
1093 mchists.
fHdp->Fill(dp);
1094 mchists.
fHppull->Fill(dp / std::sqrt(cov(4, 4)));
1096 mchists.
fHpvspc->Fill(mcmom, mom);
1097 mchists.
fHdpc->Fill(dp);
1098 mchists.
fHppullc->Fill(dp / std::sqrt(cov(4, 4)));
1112 std::set<int> tkidset;
1113 tkidset.insert(mcid);
1115 clockData, tkidset, trackhits, allhits,
geo::k3D);
1125 mchists.
fHgendx->Fill(mcend.X());
1126 mchists.
fHgendy->Fill(mcend.Y());
1127 mchists.
fHgendz->Fill(mcend.Z());
1128 mchists.
fHgtheta->Fill(mcstartmom.Theta());
1129 mchists.
fHgphi->Fill(mcstartmom.Phi());
1132 mchists.
fHgmom->Fill(mcmom);
1134 mchists.
fHgke->Fill(mcke);
1135 mchists.
fHgkel->Fill(mcke);
1136 mchists.
fHglen->Fill(plen);
1140 selected_mctracks[imc].second = i;
1144 float KE = ptkl->
E() - ptkl->
Mass();
1145 std::string KEUnits =
" GeV";
1152 << evt.
run() <<
"." << evt.
event() <<
" Match MCTkID " << std::setw(6)
1153 << mctrk.
TrackID() <<
" Origin " << mctrk.
Origin() <<
" PDG" << std::setw(5)
1154 << mctrk.
PdgCode() <<
" KE" << std::setw(4) << (int)KE << KEUnits
1155 <<
" RecoTrkID " << track.
ID() <<
" hitEff " << std::setprecision(2)
1156 << hiteff <<
" hitPur " << hitpurity;
1158 for (
unsigned short ipl = 0; ipl < wireReadoutGeom.Nplanes(); ++ipl) {
1160 auto const& plane = wireReadoutGeom.
Plane(planeID);
1161 auto const sWire = plane.NearestWireID(mcstart).Wire;
1163 auto const eWire = plane.NearestWireID(mcend).Wire;
1166 <<
" Wire:Tick in Pln " << ipl <<
" W:T " << sWire <<
":" << sTick
1167 <<
" - " << eWire <<
":" << eTick;
1179 TVector3 pos = track.
Vertex<TVector3>();
1181 TVector3 end = track.
End<TVector3>();
1187 *pdump <<
"\nOffset" << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " 1188 << std::fixed << std::setprecision(2) << std::setw(10) << trackdx <<
"\nStart " 1189 << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " << std::fixed
1190 << std::setprecision(2) << std::setw(10) << pos[0] << std::setw(10) << pos[1]
1191 << std::setw(10) << pos[2];
1193 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
1194 << std::setw(10) << dir[1] << std::setw(10) << dir[2];
1197 *pdump << std::setw(32) <<
" ";
1198 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pstart;
1199 *pdump <<
"\nEnd " << std::setw(3) << track.
ID() << std::setw(6) << mcid <<
" " 1200 << std::fixed << std::setprecision(2) << std::setw(10) << end[0] << std::setw(10)
1201 << end[1] << std::setw(10) << end[2];
1203 *pdump <<
" " << std::fixed << std::setprecision(3) << std::setw(10) << enddir[0]
1204 << std::setw(10) << enddir[1] << std::setw(10) << enddir[2];
1207 *pdump << std::setw(32) <<
" ";
1208 *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
1209 <<
"\nLength: " << tlen <<
"\n";
1217 for (
unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1218 if (selected_mctracks[imc].
second >= 0)
continue;
1219 const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1221 float KE = ptkl->
E() - ptkl->
Mass();
1222 std::string KEUnits =
" GeV";
1230 TVector3 mcstartmom, mcendmom;
1232 double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1234 << evt.
run() <<
"." << evt.
event() <<
" NoMat MCTkID " << std::setw(6) << mctrk.
TrackID()
1235 <<
" Origin " << mctrk.
Origin() <<
" PDG" << std::setw(5) << mctrk.
PdgCode() <<
" KE" 1236 << std::setw(4) << (int)KE << KEUnits <<
" Length " << std::fixed << std::setprecision(1)
1240 for (
unsigned short ipl = 0; ipl < wireReadoutGeom.Nplanes(); ++ipl) {
1242 auto const& plane = wireReadoutGeom.
Plane(planeID);
1243 auto const sWire = plane.NearestWireID(mcstart).Wire;
1245 auto const eWire = plane.NearestWireID(mcend).Wire;
1247 mf::LogVerbatim(
"TrackAna") <<
" Wire:Tick in Pln " << ipl <<
" W:T " << sWire <<
":" 1248 << sTick <<
" - " << eWire <<
":" << eTick;
1262 std::map<int, std::map<int, art::PtrVector<recob::Hit>>> hitmap;
1263 std::map<int, int> KEmap;
1272 int ntv(trackvh->size());
1274 std::vector<art::PtrVector<recob::Track>>
::const_iterator cti = trackvh->begin();
1278 int nsppts_assnwhole = fswhole.size();
1279 std::cout <<
"TrackAna: Number of clumps of Spacepoints from Assn for all Tracks: " 1280 << nsppts_assnwhole << std::endl;
1286 <<
"\t\t TrkAna: Fresh fRecoHistMap[0] ******* \n" 1291 std::vector<std::vector<unsigned int>> NtrkIdsAll;
1292 std::vector<double> ntvsorted;
1297 for (
int o = 0; o < ntv; ++o)
1300 int ntrack = pvtrack.
size();
1301 std::vector<std::vector<unsigned int>> NtrkId_Hit;
1302 std::vector<unsigned int> vecMode;
1305 for (
int i = 0; i < ntrack; ++i) {
1314 int nsppts_assn = fs.at(i).size();
1316 const auto&
sppt = fs.at(i);
1326 std::vector<unsigned int> vecNtrkIds;
1327 for (
int is = 0; is < nsppts_assn; ++is) {
1328 int nhits = fh.at(is).size();
1329 for (
int ih = 0; ih < nhits; ++ih) {
1330 const auto&
hit = fh.at(is).at(ih);
1342 int trackID =
std::abs(itid->trackID);
1343 hitmap[trackID][o].push_back(
hit);
1346 vecNtrkIds.push_back(trackID);
1359 TVector3 mcstartmom;
1361 double mctime = part->
T();
1364 double plen = length(detProp, *part, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1366 KEmap[(int)(1e6 * plen)] = trackID;
1376 NtrkId_Hit.push_back(vecNtrkIds);
1379 int max(-12),
n(1), ind(0);
1380 std::sort(vecNtrkIds.begin(), vecNtrkIds.end());
1381 std::vector<unsigned int> strkIds(vecNtrkIds);
1382 while (ii < vecNtrkIds.size()) {
1383 if (strkIds.at(ii) != strkIds.at(ii - 1)) { n = 1; }
1394 if (strkIds.begin() != strkIds.end()) mode = strkIds.at(ind);
1395 vecMode.push_back(mode);
1397 if (strkIds.size() != 0)
1398 rhistsStitched.
fModeFrac->Fill((
double)max / (
double)strkIds.size());
1407 if (!assns)
throw cet::exception(
"TrackAna") <<
"Bad Associations. \n";
1413 auto& back = NtrkIdsAll.emplace_back(vecMode);
1415 auto const last = std::unique(back.begin(), back.end());
1416 back.erase(last, back.end());
1418 for (
auto const val : back) {
1419 sum += hitmap[val][o].size();
1421 ntvsorted.push_back(sum);
1429 for (
auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1437 for (
auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1438 int val = it->second;
1446 flattener flat(NtrkIdsAll);
1447 std::vector<unsigned int>& v = flat;
1450 for (
auto const val : v) {
1453 double T(part->
E() - 0.001 * part->
Mass());
1454 rhistsStitched.
fNTrkIdTrks->Fill(std::count(v.begin(), v.end(), val));
1455 rhistsStitched.
fNTrkIdTrks2->Fill(std::count(v.begin(), v.end(), val), T);
1467 mf::LogInfo(
"TrackAna") <<
"TrackAna statistics:\n" 1474 const MCHists& mchists = i->second;
1498 std::vector<size_t> idx(v.size());
1499 for (
size_t i = 0; i != idx.size(); ++i)
1502 std::sort(idx.begin(), idx.end(), [&v](
size_t i1,
size_t i2) {
1503 return v[i1] > v[i2];
double E(const int i=0) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
unsigned int NumberTrajectoryPoints() const
const TLorentzVector & Position(const int i=0) const
simb::Origin_t Origin() const
double EndMomentum() const
double VertexMomentum() const
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::string fHitSpptAssocModuleLabel
std::string fTrkSpptAssocModuleLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const simb::MCParticle * TrackIdToParticle_P(int id) const
std::string fSpacepointModuleLabel
Declaration of signal hit object.
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.
enum simb::_ev_origin Origin_t
event origin types
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.
unsigned int ReadOutWindowSize() const
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
std::string fStitchModuleLabel
tick ticks
Alias for common language habits.
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.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
double Length(size_t p=0) const
Access to various track properties.
#define DEFINE_ART_MODULE(klass)
std::vector< size_t > fsort_indexes(const std::vector< double > &v)
static const int NoParticleId
double HitCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
single particles thrown at the detector
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.
double T(const int i=0) const
double ConvertXToTicks(double X, int p, int t, int c) const
const SMatrixSym55 & EndCovariance() const
Access to covariance matrices.
std::string fMCTrackModuleLabel
EventNumber_t event() const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
PlaneID_t Plane
Index of the plane within its TPC.
Class def header for mctrack data container.
void analyze(const art::Event &evt) override
Detector simulation of raw signals on wires.
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.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Contains all timing reference information for the detector.
Provides recob::Track data product.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Vector_t EndDirection() const
Access to track direction at different points.
std::map< int, RecoHists > fRecoHistMap
const TLorentzVector & Momentum(const int i=0) const
Point_t const & End() const
Access to track position at different points.
const MCStep & Start() const
size_type get(size_type i, reference item, data_reference data) const
std::string fTrackModuleLabel
unsigned int TrackID() const
second_as<> second
Type of time stored in seconds, in double precision.
std::string fHitModuleLabel
Tools and modules for checking out the basics of the Monte Carlo.
art framework interface to geometry description
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
Signal from collection planes.