LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
TrackAna_module.cc
Go to the documentation of this file.
1 //
2 // Name: TrackAna_module.cc
3 //
4 // Purpose: Module TrackAna.
5 //
6 // Configuration parameters.
7 //
8 // TrackModuleLabel: Track module label.
9 // MCTrackModuleLabel: MCTrack module label.
10 // MinMCKE: Minimum MC particle kinetic energy.
11 // MatchColinearity: Minimum colinearity for mc-track matching.
12 // MatchDisp: Maximum uv displacement for mc-track matching.
13 // WMatchDisp: maximum w displacement for mc-track matching.
14 // MatchLength: Minimum length fraction for mc-track matching.
15 //
16 // Created: 2-Aug-2011 H. Greenlee
17 //
18 
19 #include <cmath>
20 #include <iomanip>
21 #include <iostream>
22 #include <map>
23 #include <memory>
24 #include <sstream>
25 
30 #include "art_root_io/TFileService.h"
33 #include "cetlib_except/exception.h"
35 
48 
49 #include "TH1F.h"
50 #include "TH2F.h"
51 #include "TMatrixD.h"
52 
53 namespace {
54 
55  // Calculate distance to boundary.
56  //----------------------------------------------------------------------------
57  double bdist(const TVector3& pos)
58  {
59  // Get geometry.
60 
61  auto const& tpc = art::ServiceHandle<geo::Geometry const>()->TPC({0, 0});
62 
63  double d1 = pos.X(); // Distance to right side (wires).
64  double d2 = 2. * tpc.HalfWidth() - pos.X(); // Distance to left side (cathode).
65  double d3 = pos.Y() + tpc.HalfHeight(); // Distance to bottom.
66  double d4 = tpc.HalfHeight() - pos.Y(); // Distance to top.
67  double d5 = pos.Z(); // Distance to front.
68  double d6 = tpc.Length() - pos.Z(); // Distance to back.
69 
70  return std::min({d1, d2, d3, d4, d5, d6});
71  }
72 
73  // Length of reconstructed track.
74  //----------------------------------------------------------------------------
75  double length(const recob::Track& track)
76  {
77  return track.Length();
78  }
79 
80  // Length of MC particle.
81  //----------------------------------------------------------------------------
82  double length(detinfo::DetectorPropertiesData const& detProp,
83  const simb::MCParticle& part,
84  double dx,
85  geo::Point_t& start,
87  TVector3& startmom,
88  TVector3& endmom)
89  {
90  auto const& tpc = art::ServiceHandle<geo::Geometry const>()->TPC({0, 0});
91 
92  // Get fiducial volume boundary.
93  double xmin = 0.;
94  double xmax = 2. * tpc.HalfWidth();
95  double ymin = -tpc.HalfHeight();
96  double ymax = tpc.HalfHeight();
97  double zmin = 0.;
98  double zmax = tpc.Length();
99  double result = 0.;
100  TVector3 disp;
101  int n = part.NumberTrajectoryPoints();
102  bool first = true;
103 
104  for (int i = 0; i < n; ++i) {
105  TVector3 pos = part.Position(i).Vect();
106 
107  // Make fiducial cuts. Require the particle to be within the physical volume of
108  // the tpc, and also require the apparent x position to be within the expanded
109  // readout frame.
110 
111  if (pos.X() >= xmin && pos.X() <= xmax && pos.Y() >= ymin && pos.Y() <= ymax &&
112  pos.Z() >= zmin && pos.Z() <= zmax) {
113  pos[0] += dx;
114  double ticks = detProp.ConvertXToTicks(pos[0], 0, 0, 0);
115  if (ticks >= 0. && ticks < detProp.ReadOutWindowSize()) {
116  if (first) {
117  start = pos;
118  startmom = part.Momentum(i).Vect();
119  }
120  else {
121  disp -= pos;
122  result += disp.Mag();
123  }
124  first = false;
125  disp = pos;
126  end = pos;
127  endmom = part.Momentum(i).Vect();
128  }
129  }
130  }
131 
132  return result;
133  }
134 
135  // Length of MCTrack.
136  // In this function, the extracted start and end momenta are converted to GeV
137  // (MCTrack stores momenta in Mev).
138  //----------------------------------------------------------------------------
139  double length(detinfo::DetectorPropertiesData const& detProp,
140  const sim::MCTrack& mctrk,
141  double dx,
142  geo::Point_t& start,
143  geo::Point_t& end,
144  TVector3& startmom,
145  TVector3& endmom)
146  {
147  auto const& tpc = art::ServiceHandle<geo::Geometry const>()->TPC({0, 0});
148 
149  // Get fiducial volume boundary.
150 
151  double xmin = 0.;
152  double xmax = 2. * tpc.HalfWidth();
153  double ymin = -tpc.HalfHeight();
154  double ymax = tpc.HalfHeight();
155  double zmin = 0.;
156  double zmax = tpc.Length();
157  double result = 0.;
158  TVector3 disp;
159  int n = mctrk.size();
160  bool first = true;
161 
162  for (int i = 0; i < n; ++i) {
163  TVector3 pos = mctrk[i].Position().Vect();
164 
165  // Make fiducial cuts. Require the particle to be within the physical volume of
166  // the tpc, and also require the apparent x position to be within the expanded
167  // readout frame.
168 
169  if (pos.X() >= xmin && pos.X() <= xmax && pos.Y() >= ymin && pos.Y() <= ymax &&
170  pos.Z() >= zmin && pos.Z() <= zmax) {
171  pos[0] += dx;
172  double ticks = detProp.ConvertXToTicks(pos[0], 0, 0, 0);
173  if (ticks >= 0. && ticks < detProp.ReadOutWindowSize()) {
174  if (first) {
175  start = pos;
176  startmom = 0.001 * mctrk[i].Momentum().Vect();
177  }
178  else {
179  disp -= pos;
180  result += disp.Mag();
181  }
182  first = false;
183  disp = pos;
184  end = pos;
185  endmom = 0.001 * mctrk[i].Momentum().Vect();
186  }
187  }
188  }
189 
190  return result;
191  }
192 
193  // Fill efficiency histogram assuming binomial errors.
194 
195  void effcalc(const TH1* hnum, const TH1* hden, TH1* heff)
196  {
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";
202 
203  // Loop over bins, including underflow and overflow.
204 
205  for (int ibin = 0; ibin <= nbins + 1; ++ibin) {
206  double num = hnum->GetBinContent(ibin);
207  double den = hden->GetBinContent(ibin);
208  if (den == 0.) {
209  heff->SetBinContent(ibin, 0.);
210  heff->SetBinError(ibin, 0.);
211  }
212  else {
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);
219  }
220  }
221 
222  heff->SetMinimum(0.);
223  heff->SetMaximum(1.05);
224  heff->SetMarkerStyle(20);
225  }
226 
227  class flattener : public std::vector<unsigned int> {
228  public:
229  flattener() = default;
230 
231  flattener(const std::vector<std::vector<unsigned int>>& input) { Convert(input); }
232 
233  void Convert(const std::vector<std::vector<unsigned int>>& input)
234  {
235  clear();
236  size_t length = 0;
237  for (auto const& vec : input)
238  length += vec.size();
239  reserve(length);
240 
241  for (auto const& vec : input)
242  for (auto const& value : vec)
243  push_back(value);
244  }
245  }; // end class flattener
246 
247 } // end namespace
248 
249 namespace trkf {
250 
251  class TrackAna : public art::EDAnalyzer {
252  public:
253  // Embedded structs.
254 
255  // Struct for histograms that depend on reco track only.
256 
257  struct RecoHists {
258  explicit RecoHists(const std::string& subdir);
259 
260  // Pure reco track histograms.
261 
262  TH1F* fHstartx{nullptr}; // Starting x position.
263  TH1F* fHstarty{nullptr}; // Starting y position.
264  TH1F* fHstartz{nullptr}; // Starting z position.
265  TH1F* fHstartd{nullptr}; // Starting distance to boundary.
266  TH1F* fHendx{nullptr}; // Ending x position.
267  TH1F* fHendy{nullptr}; // Ending y position.
268  TH1F* fHendz{nullptr}; // Ending z position.
269  TH1F* fHendd{nullptr}; // Ending distance to boundary.
270  TH1F* fHtheta{nullptr}; // Theta.
271  TH1F* fHphi{nullptr}; // Phi.
272  TH1F* fHtheta_xz{nullptr}; // Theta_xz.
273  TH1F* fHtheta_yz{nullptr}; // Theta_yz.
274  TH1F* fHmom{nullptr}; // Momentum.
275  TH1F* fHmoml{nullptr}; // Momentum (low momentum).
276  TH1F* fHlen{nullptr}; // Length.
277  TH1F* fHlens{nullptr}; // Length (short tracks).
278 
279  // Histograms for the consituent Hits
280 
281  TH1F* fHHitChg{nullptr}; // hit charge
282  TH1F* fHHitWidth{nullptr}; // hit width
283  TH1F* fHHitPdg{nullptr}; // Pdg primarily responsible.
284  TH1F* fHHitTrkId{nullptr}; // TrkId
285  TH1F* fModeFrac{nullptr}; // mode fraction
286  TH1F* fNTrkIdTrks{nullptr}; // # of stitched tracks in which unique TrkId appears
287  TH2F* fNTrkIdTrks2{nullptr};
288  TH2F* fNTrkIdTrks3{nullptr};
289  };
290 
291  // Struct for mc particles and mc-matched tracks.
292 
293  struct MCHists {
294  explicit MCHists(const std::string& subdir);
295 
296  // Reco-MC matching.
297 
298  TH2F* fHduvcosth{nullptr}; // 2D mc vs. data matching, duv vs. cos(theta).
299  TH1F* fHcosth{nullptr}; // 1D direction matching, cos(theta).
300  TH1F* fHmcu{nullptr}; // 1D endpoint truth u.
301  TH1F* fHmcv{nullptr}; // 1D endpoint truth v.
302  TH1F* fHmcw{nullptr}; // 1D endpoint truth w.
303  TH1F* fHupull{nullptr}; // 1D endpoint u pull.
304  TH1F* fHvpull{nullptr}; // 1D endpoint v pull.
305  TH1F* fHmcdudw{nullptr}; // Truth du/dw.
306  TH1F* fHmcdvdw{nullptr}; // Truth dv/dw.
307  TH1F* fHdudwpull{nullptr}; // du/dw pull.
308  TH1F* fHdvdwpull{nullptr}; // dv/dw pull.
309  TH1F* fHHitEff{nullptr}; // Hit efficiency.
310  TH1F* fHHitPurity{nullptr}; // Hit purity.
311 
312  // Histograms for matched tracks.
313 
314  TH1F* fHstartdx{nullptr}; // Start dx.
315  TH1F* fHstartdy{nullptr}; // Start dy.
316  TH1F* fHstartdz{nullptr}; // Start dz.
317  TH1F* fHenddx{nullptr}; // End dx.
318  TH1F* fHenddy{nullptr}; // End dy.
319  TH1F* fHenddz{nullptr}; // End dz.
320  TH2F* fHlvsl{nullptr}; // MC vs. reco length.
321  TH1F* fHdl{nullptr}; // Delta(length).
322  TH2F* fHpvsp{nullptr}; // MC vs. reco momentum.
323  TH2F* fHpvspc{nullptr}; // MC vs. reco momentum (contained tracks).
324  TH1F* fHdp{nullptr}; // Momentum difference.
325  TH1F* fHdpc{nullptr}; // Momentum difference (contained tracks).
326  TH1F* fHppull{nullptr}; // Momentum pull.
327  TH1F* fHppullc{nullptr}; // Momentum pull (contained tracks).
328 
329  // Pure MC particle histograms (efficiency denominator).
330 
331  TH1F* fHmcstartx{nullptr}; // Starting x position.
332  TH1F* fHmcstarty{nullptr}; // Starting y position.
333  TH1F* fHmcstartz{nullptr}; // Starting z position.
334  TH1F* fHmcendx{nullptr}; // Ending x position.
335  TH1F* fHmcendy{nullptr}; // Ending y position.
336  TH1F* fHmcendz{nullptr}; // Ending z position.
337  TH1F* fHmctheta{nullptr}; // Theta.
338  TH1F* fHmcphi{nullptr}; // Phi.
339  TH1F* fHmctheta_xz{nullptr}; // Theta_xz.
340  TH1F* fHmctheta_yz{nullptr}; // Theta_yz.
341  TH1F* fHmcmom{nullptr}; // Momentum.
342  TH1F* fHmcmoml{nullptr}; // Momentum (low momentum).
343  TH1F* fHmcke{nullptr}; // Kinetic energy.
344  TH1F* fHmckel{nullptr}; // Kinetic energy (low energy).
345  TH1F* fHmclen{nullptr}; // Length.
346  TH1F* fHmclens{nullptr}; // Length (short tracks).
347 
348  // Histograms for well-reconstructed matched tracks (efficiency numerator).
349 
350  TH1F* fHgstartx{nullptr}; // Starting x position.
351  TH1F* fHgstarty{nullptr}; // Starting y position.
352  TH1F* fHgstartz{nullptr}; // Starting z position.
353  TH1F* fHgendx{nullptr}; // Ending x position.
354  TH1F* fHgendy{nullptr}; // Ending y position.
355  TH1F* fHgendz{nullptr}; // Ending z position.
356  TH1F* fHgtheta{nullptr}; // Theta.
357  TH1F* fHgphi{nullptr}; // Phi.
358  TH1F* fHgtheta_xz{nullptr}; // Theta_xz.
359  TH1F* fHgtheta_yz{nullptr}; // Theta_yz.
360  TH1F* fHgmom{nullptr}; // Momentum.
361  TH1F* fHgmoml{nullptr}; // Momentum (low momentum).
362  TH1F* fHgke{nullptr}; // Kinetic energy.
363  TH1F* fHgkel{nullptr}; // Kinetic energy (low momentum).
364  TH1F* fHglen{nullptr}; // Length.
365  TH1F* fHglens{nullptr}; // Length (short tracks).
366 
367  // Efficiency histograms.
368 
369  TH1F* fHestartx{nullptr}; // Starting x position.
370  TH1F* fHestarty{nullptr}; // Starting y position.
371  TH1F* fHestartz{nullptr}; // Starting z position.
372  TH1F* fHeendx{nullptr}; // Ending x position.
373  TH1F* fHeendy{nullptr}; // Ending y position.
374  TH1F* fHeendz{nullptr}; // Ending z position.
375  TH1F* fHetheta{nullptr}; // Theta.
376  TH1F* fHephi{nullptr}; // Phi.
377  TH1F* fHetheta_xz{nullptr}; // Theta_xz.
378  TH1F* fHetheta_yz{nullptr}; // Theta_yz.
379  TH1F* fHemom{nullptr}; // Momentum.
380  TH1F* fHemoml{nullptr}; // Momentum (low momentum).
381  TH1F* fHeke{nullptr}; // Kinetic energy.
382  TH1F* fHekel{nullptr}; // Kinetic energy (low momentum).
383  TH1F* fHelen{nullptr}; // Length.
384  TH1F* fHelens{nullptr}; // Length (short tracks).
385  };
386 
387  // Constructors, destructor
388 
389  explicit TrackAna(fhicl::ParameterSet const& pset);
390 
391  private:
392  void analyze(const art::Event& evt) override;
393  void endJob() override;
394 
395  void anaStitch(const art::Event& evt,
396  detinfo::DetectorClocksData const& clockData,
397  detinfo::DetectorPropertiesData const& detProp);
398 
399  std::vector<size_t> fsort_indexes(const std::vector<double>& v);
400 
401  // Fcl Attributes.
402 
403  std::string fTrackModuleLabel;
404  std::string fMCTrackModuleLabel;
406  std::string fStitchModuleLabel;
409  std::string fHitModuleLabel;
410 
411  int fDump; // Number of events to dump to debug message facility.
412  double fMinMCKE; // Minimum MC particle kinetic energy (GeV).
413  double fMinMCLen; // Minimum MC particle length in tpc (cm).
414  double fMatchColinearity; // Minimum matching colinearity.
415  double fMatchDisp; // Maximum matching displacement.
416  double fWMatchDisp; // Maximum matching displacement in the w direction.
417  double fMatchLength; // Minimum length fraction.
418  bool fIgnoreSign; // Ignore sign of mc particle if true.
419  bool fStitchedAnalysis; // if true, do the whole drill-down from stitched track to assd hits
420 
421  std::string fOrigin;
424  int fPrintLevel; // 0 = none, 1 = event summary, 2 = track detail
425 
426  // Histograms.
427 
428  std::map<int, MCHists> fMCHistMap; // Indexed by pdg id.
429  std::map<int, RecoHists> fRecoHistMap; // Indexed by pdg id.
430 
431  // Statistics.
432 
434  };
435 
437 
438  // RecoHists methods.
439 
440  TrackAna::RecoHists::RecoHists(const std::string& subdir)
441  {
442  // Get services.
443 
444  auto const& tpc = art::ServiceHandle<geo::Geometry const>()->TPC({0, 0});
446 
447  // Make histogram directory.
448 
449  art::TFileDirectory topdir = tfs->mkdir("trkana", "TrackAna histograms");
450  art::TFileDirectory dir = topdir.mkdir(subdir);
451 
452  // Book histograms.
453 
454  fHstartx = dir.make<TH1F>(
455  "xstart", "X Start Position", 100, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
456  fHstarty =
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());
459  fHstartd =
460  dir.make<TH1F>("dstart", "Start Position Distance to Boundary", 100, -10., tpc.HalfWidth());
461  fHendx =
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());
465  fHendd =
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);
479  fModeFrac =
480  dir.make<TH1F>("hmodefrac",
481  "quasi-Purity: Fraction of component tracks with the Track mode value",
482  20,
483  0.01,
484  1.01);
485  fNTrkIdTrks =
486  dir.make<TH1F>("hntrkids",
487  "quasi-Efficiency: Number of stitched tracks in which TrkId appears",
488  20,
489  0.,
490  +10.0);
491  fNTrkIdTrks2 = dir.make<TH2F>("hntrkids2",
492  "Number of stitched tracks in which TrkId appears vs KE [GeV]",
493  20,
494  0.,
495  +10.0,
496  20,
497  0.0,
498  1.5);
499  fNTrkIdTrks3 = dir.make<TH2F>("hntrkids3",
500  "MC Track vs Reco Track, wtd by nhits on Collection Plane",
501  10,
502  -0.5,
503  9.5,
504  10,
505  -0.5,
506  9.5);
507  fNTrkIdTrks3->GetXaxis()->SetTitle("Sorted-by-Descending-CPlane-Hits outer Track Number");
508  fNTrkIdTrks3->GetYaxis()->SetTitle("Sorted-by-Descending-True-Length G4Track");
509  }
510 
511  // MCHists methods.
512 
513  TrackAna::MCHists::MCHists(const std::string& subdir)
514  {
515  // Get services.
516 
517  auto const& tpc = art::ServiceHandle<geo::Geometry const>()->TPC({0, 0});
519 
520  // Make histogram directory.
521 
522  art::TFileDirectory topdir = tfs->mkdir("trkana", "TrackAna histograms");
523  art::TFileDirectory dir = topdir.mkdir(subdir);
524 
525  // Book histograms.
526 
527  fHduvcosth =
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",
549  100,
550  0.,
551  1.1 * tpc.Length(),
552  100,
553  0.,
554  1.1 * tpc.Length());
555  fHdl = dir.make<TH1F>("dl", "Track Length Minus MC Particle Length", 100, -50., 50.);
556  fHpvsp =
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.);
564 
565  fHmcstartx = dir.make<TH1F>(
566  "mcxstart", "MC X Start Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
567  fHmcstarty =
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());
572  fHmcendy =
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());
585 
586  fHgstartx = dir.make<TH1F>(
587  "gxstart", "Good X Start Position", 10, -2. * tpc.HalfWidth(), 4. * tpc.HalfWidth());
588  fHgstarty =
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());
593  fHgendy =
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());
606 
607  fHestartx = dir.make<TH1F>("exstart",
608  "Efficiency vs. X Start Position",
609  10,
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());
630  }
631 
632  TrackAna::TrackAna(const fhicl::ParameterSet& pset)
633  : EDAnalyzer(pset)
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))
652  , fNumEvent(0)
653  {
654 
655  // Decide whether to check MCTrack origin
656  fCheckOrigin = false;
658  if (fOrigin.find("Beam") != std::string::npos) {
659  fCheckOrigin = true;
661  }
662  else if (fOrigin.find("Cosmic") != std::string::npos) {
663  fCheckOrigin = true;
665  }
666  else if (fOrigin.find("Super") != std::string::npos) {
667  fCheckOrigin = true;
669  }
670  else if (fOrigin.find("Single") != std::string::npos) {
671  fCheckOrigin = true;
673  }
674 
675  // Report.
676 
677  mf::LogInfo("TrackAna") << "TrackAna configured with the following parameters:\n"
678  << " TrackModuleLabel = " << fTrackModuleLabel << "\n"
679  << " MCTrackModuleLabel = " << fMCTrackModuleLabel << "\n"
680  << " StitchModuleLabel = " << fStitchModuleLabel << "\n"
681  << " TrkSpptAssocModuleLabel = " << fTrkSpptAssocModuleLabel << "\n"
682  << " HitSpptAssocModuleLabel = " << fHitSpptAssocModuleLabel << "\n"
683  << " HitModuleLabel = " << fHitModuleLabel << "\n"
684  << " Dump = " << fDump << "\n"
685  << " MinMCKE = " << fMinMCKE << "\n"
686  << " MinMCLen = " << fMinMCLen << " Origin = " << fOrigin
687  << " Origin value " << fOriginValue;
688  }
689 
691  {
694  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
695 
696  ++fNumEvent;
697 
698  // Optional dump stream.
699 
700  std::unique_ptr<mf::LogInfo> pdump;
701  if (fDump > 0) {
702  --fDump;
703  pdump = std::make_unique<mf::LogInfo>("TrackAna");
704  }
705 
706  // Make sure histograms are booked.
707 
708  bool mc = !evt.isRealData();
709 
710  // Get MCTracks.
711 
713  evt.getByLabel(fMCTrackModuleLabel, mctrackh);
714  // pair of MCTrack and index of matched reco track
715  std::vector<std::pair<const sim::MCTrack*, int>> selected_mctracks;
716 
717  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
718  auto const detProp =
720 
721  if (mc && mctrackh.isValid()) {
722 
723  selected_mctracks.reserve(mctrackh->size());
724 
725  // Dump MCTracks.
726 
727  if (pdump) {
728  *pdump << "MC Tracks\n"
729  << " Id pdg x y z dx dy dz "
730  " p\n"
731  << "--------------------------------------------------------------------------------"
732  "-----------\n";
733  }
734 
735  // Loop over mc tracks, and fill histograms that depend only
736  // on mc information. Also, fill a secondary list of mc tracks
737  // that pass various selection criteria.
738 
739  for (std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();
740  imctrk != mctrackh->end();
741  ++imctrk) {
742  const sim::MCTrack& mctrk = *imctrk;
743  int pdg = mctrk.PdgCode();
744  if (fIgnoreSign) pdg = std::abs(pdg);
745 
746  // Ignore everything except stable charged nonshowering particles.
747 
748  int apdg = std::abs(pdg);
749  if (apdg == 13 || // Muon
750  apdg == 211 || // Charged pion
751  apdg == 321 || // Charged kaon
752  apdg == 2212) { // (Anti)proton
753 
754  // check MC track origin?
755  if (fCheckOrigin && mctrk.Origin() != fOriginValue) continue;
756 
757  // Apply minimum energy cut.
758 
759  if (mctrk.Start().E() >= mctrk.Start().Momentum().Mag() + 1000. * fMinMCKE) {
760 
761  // Calculate the x offset due to nonzero mc particle time.
762 
763  double mctime = mctrk.Start().T(); // nsec
764  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
765 
766  // Calculate the length of this mc particle inside the fiducial volume.
767 
768  geo::Point_t mcstart;
769  geo::Point_t mcend;
770  TVector3 mcstartmom;
771  TVector3 mcendmom;
772  double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
773 
774  // Apply minimum fiducial length cut. Always reject particles that have
775  // zero length in the tpc regardless of the configured cut.
776 
777  if (plen > 0. && plen > fMinMCLen) {
778 
779  // This is a good mc particle (capable of making a track).
780 
781  selected_mctracks.push_back(std::make_pair(&mctrk, -1));
782 
783  // Dump MC particle information here.
784 
785  if (pdump) {
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();
794  if (pstart > 0.) {
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;
798  }
799  else
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)
805  << mcend.Z();
806  if (pend > 0.01) {
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;
810  }
811  else
812  *pdump << std::setw(32) << " ";
813  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
814  << "\nLength: " << plen << "\n";
815  }
816 
817  // Fill histograms.
818 
819  if (fMCHistMap.count(pdg) == 0) {
820  std::ostringstream ostr;
821  ostr << "MC" << (fIgnoreSign ? "All" : (pdg > 0 ? "Pos" : "Neg")) << std::abs(pdg);
822  fMCHistMap.emplace(pdg, MCHists{ostr.str()});
823  }
824  const MCHists& mchists = fMCHistMap.at(pdg);
825 
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();
829  double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
830  double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
831 
832  mchists.fHmcstartx->Fill(mcstart.X());
833  mchists.fHmcstarty->Fill(mcstart.Y());
834  mchists.fHmcstartz->Fill(mcstart.Z());
835  mchists.fHmcendx->Fill(mcend.X());
836  mchists.fHmcendy->Fill(mcend.Y());
837  mchists.fHmcendz->Fill(mcend.Z());
838  mchists.fHmctheta->Fill(mcstartmom.Theta());
839  mchists.fHmcphi->Fill(mcstartmom.Phi());
840  mchists.fHmctheta_xz->Fill(mctheta_xz);
841  mchists.fHmctheta_yz->Fill(mctheta_yz);
842  mchists.fHmcmom->Fill(mcmom);
843  mchists.fHmcmoml->Fill(mcmom);
844  mchists.fHmcke->Fill(mcke);
845  mchists.fHmckel->Fill(mcke);
846  mchists.fHmclen->Fill(plen);
847  mchists.fHmclens->Fill(plen);
848  }
849  }
850  }
851  }
852  } //mc
853 
854  // Get tracks and spacepoints and hits
858 
859  evt.getByLabel(fTrackModuleLabel, trackh);
860  evt.getByLabel(fStitchModuleLabel, trackvh);
861  evt.getByLabel(fHitModuleLabel, hith);
862 
863  // Extract all hits into a vector of art::Ptrs (the format used by back tracker).
864 
865  std::vector<art::Ptr<recob::Hit>> allhits;
866  if (hith.isValid()) {
867  allhits.reserve(hith->size());
868  for (unsigned int i = 0; i < hith->size(); ++i) {
869  allhits.emplace_back(hith, i);
870  }
871  }
872 
873  // Construct FindManyP object to be used for finding track-hit associations.
874 
875  art::FindManyP<recob::Hit> tkhit_find(trackh, evt, fTrackModuleLabel);
876 
877  // This new top part of TrackAna between two long lines of ************s
878  // is particular to analyzing Stitched Tracks.
879  // *******************************************************************//
880 
881  if (trackvh.isValid() && fStitchedAnalysis) {
882  mf::LogDebug("TrackAna") << "TrackAna read " << trackvh->size()
883  << " vectors of Stitched PtrVectorsof tracks.";
884  anaStitch(evt, clockData, detProp);
885  }
886 
887  if (trackh.isValid()) {
888 
889  if (pdump) {
890  *pdump << "\nReconstructed Tracks\n"
891  << " Id MCid x y z dx dy dz "
892  " p\n"
893  << "--------------------------------------------------------------------------------"
894  "-----------\n";
895  }
896 
897  // Loop over tracks.
898 
899  int ntrack = trackh->size();
900  for (int i = 0; i < ntrack; ++i) {
901  art::Ptr<recob::Track> ptrack(trackh, i);
902  const recob::Track& track = *ptrack;
903 
904  // Extract hits associated with this track.
905 
906  std::vector<art::Ptr<recob::Hit>> trackhits;
907  tkhit_find.get(i, trackhits);
908 
909  // Calculate the x offset due to nonzero reconstructed time.
910 
911  double trackdx = 0;
912 
913  // Fill histograms involving reco tracks only.
914 
915  int ntraj = track.NumberTrajectoryPoints();
916  if (ntraj > 0) {
917  TVector3 pos = track.Vertex<TVector3>();
918  TVector3 dir = track.VertexDirection<TVector3>();
919  TVector3 end = track.End<TVector3>();
920  pos[0] += trackdx;
921  end[0] += trackdx;
922 
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());
928 
929  if (fRecoHistMap.count(0) == 0) fRecoHistMap.emplace(0, RecoHists{"Reco"});
930  const RecoHists& rhists = fRecoHistMap.at(0);
931 
932  rhists.fHstartx->Fill(pos.X());
933  rhists.fHstarty->Fill(pos.Y());
934  rhists.fHstartz->Fill(pos.Z());
935  rhists.fHstartd->Fill(dpos);
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());
942  rhists.fHtheta_xz->Fill(theta_xz);
943  rhists.fHtheta_yz->Fill(theta_yz);
944 
945  double mom = 0.;
946  if (track.HasMomentum()) mom = track.VertexMomentum();
947  rhists.fHmom->Fill(mom);
948  rhists.fHmoml->Fill(mom);
949  rhists.fHlen->Fill(tlen);
950  rhists.fHlens->Fill(tlen);
951 
952  // Id of matching mc particle.
953 
954  int mcid = -1;
955 
956  // Loop over direction.
957 
958  for (int swap = 0; swap < 2; ++swap) {
959 
960  // Analyze reversed tracks only if start momentum = end momentum.
961 
962  if (swap != 0 && track.HasMomentum() &&
963  std::abs(track.VertexMomentum() - track.EndMomentum()) > 1.e-3)
964  continue;
965 
966  // Calculate the global-to-local rotation matrix.
967 
968  int start_point = (swap == 0 ? 0 : ntraj - 1);
969  TMatrixD rot = track.GlobalToLocalRotationAtPoint<TMatrixD>(start_point);
970 
971  // Update track data for reversed case.
972 
973  if (swap != 0) {
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);
980 
981  pos = track.End<TVector3>();
982  dir = -track.EndDirection<TVector3>();
983  end = track.Vertex<TVector3>();
984  pos[0] += trackdx;
985  end[0] += trackdx;
986 
987  dpos = bdist(pos);
988  dend = bdist(end);
989  theta_xz = std::atan2(dir.X(), dir.Z());
990  theta_yz = std::atan2(dir.Y(), dir.Z());
991 
992  if (track.HasMomentum()) mom = track.EndMomentum();
993  }
994 
995  // Get covariance matrix.
996 
997  const auto& cov = (swap == 0 ? track.VertexCovariance() : track.EndCovariance());
998 
999  // Loop over track-like mc particles.
1000 
1001  for (unsigned int imc = 0; imc < selected_mctracks.size(); ++imc) {
1002  const sim::MCTrack& mctrk = *selected_mctracks[imc].first;
1003  int pdg = mctrk.PdgCode();
1004  if (fIgnoreSign) pdg = std::abs(pdg);
1005  auto iMCHistMap = fMCHistMap.find(pdg);
1006  if (iMCHistMap == fMCHistMap.end())
1007  throw cet::exception("TrackAna") << "no particle with ID=" << pdg << "\n";
1008  const MCHists& mchists = iMCHistMap->second;
1009 
1010  // Calculate the x offset due to nonzero mc particle time.
1011 
1012  double mctime = mctrk.Start().T(); // nsec
1013  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
1014 
1015  // Calculate the points where this mc particle enters and leaves the
1016  // fiducial volume, and the length in the fiducial volume.
1017 
1018  geo::Point_t mcstart;
1019  geo::Point_t mcend;
1020  TVector3 mcstartmom;
1021  TVector3 mcendmom;
1022  double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1023 
1024  // Get the displacement of this mc particle in the global coordinate system.
1025 
1026  TVector3 mcpos{mcstart.X() - pos[0], mcstart.Y() - pos[1], mcstart.Z() - pos[2]};
1027 
1028  // Rotate the momentum and position to the
1029  // track-local coordinate system.
1030 
1031  TVector3 mcmoml = rot * mcstartmom;
1032  TVector3 mcposl = rot * mcpos;
1033 
1034  double colinearity = mcmoml.Z() / mcmoml.Mag();
1035 
1036  double u = mcposl.X();
1037  double v = mcposl.Y();
1038  double w = mcposl.Z();
1039 
1040  double pu = mcmoml.X();
1041  double pv = mcmoml.Y();
1042  double pw = mcmoml.Z();
1043 
1044  double dudw = pu / pw;
1045  double dvdw = pv / pw;
1046 
1047  double u0 = u - w * dudw;
1048  double v0 = v - w * dvdw;
1049  double uv0 = std::hypot(u0, v0);
1050 
1051  mchists.fHduvcosth->Fill(colinearity, uv0);
1052 
1053  if (std::abs(uv0) < fMatchDisp) {
1054 
1055  // Fill slope matching histograms.
1056 
1057  mchists.fHmcdudw->Fill(dudw);
1058  mchists.fHmcdvdw->Fill(dvdw);
1059  mchists.fHdudwpull->Fill(dudw / std::sqrt(cov(2, 2)));
1060  mchists.fHdvdwpull->Fill(dvdw / std::sqrt(cov(3, 3)));
1061  }
1062  mchists.fHcosth->Fill(colinearity);
1063  if (colinearity > fMatchColinearity) {
1064 
1065  // Fill displacement matching histograms.
1066 
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)));
1072 
1073  if (std::abs(uv0) < fMatchDisp) {
1074 
1075  // Fill matching histograms.
1076 
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();
1080  double mcmass = 0.001 * mctrk.Start().Momentum().Mag();
1081  double mcke = mcmom * mcmom / (std::hypot(mcmom, mcmass) + mcmass);
1082 
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)));
1095  if (std::abs(dpos) >= 5. && std::abs(dend) >= 5.) {
1096  mchists.fHpvspc->Fill(mcmom, mom);
1097  mchists.fHdpc->Fill(dp);
1098  mchists.fHppullc->Fill(dp / std::sqrt(cov(4, 4)));
1099  }
1100 
1101  // Count this track as well-reconstructed if it is matched to an
1102  // mc particle (which it is if get here), and if in addition the
1103  // starting position (w) matches and the reconstructed track length
1104  // is more than 0.5 of the mc particle trajectory length.
1105 
1106  bool good = std::abs(w) <= fWMatchDisp && tlen > fMatchLength * plen;
1107  if (good) {
1108  mcid = mctrk.TrackID();
1109 
1110  // Calculate and fill hit efficiency and purity.
1111 
1112  std::set<int> tkidset;
1113  tkidset.insert(mcid);
1114  double hiteff = bt_serv->HitCollectionEfficiency(
1115  clockData, tkidset, trackhits, allhits, geo::k3D);
1116  double hitpurity = bt_serv->HitCollectionPurity(clockData, tkidset, trackhits);
1117  mchists.fHHitEff->Fill(hiteff);
1118  mchists.fHHitPurity->Fill(hitpurity);
1119 
1120  // Fill efficiency numerator histograms.
1121 
1122  mchists.fHgstartx->Fill(mcstart.X());
1123  mchists.fHgstarty->Fill(mcstart.Y());
1124  mchists.fHgstartz->Fill(mcstart.Z());
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());
1130  mchists.fHgtheta_xz->Fill(mctheta_xz);
1131  mchists.fHgtheta_yz->Fill(mctheta_yz);
1132  mchists.fHgmom->Fill(mcmom);
1133  mchists.fHgmoml->Fill(mcmom);
1134  mchists.fHgke->Fill(mcke);
1135  mchists.fHgkel->Fill(mcke);
1136  mchists.fHglen->Fill(plen);
1137  mchists.fHglens->Fill(plen);
1138 
1139  // set the match flag
1140  selected_mctracks[imc].second = i;
1141 
1142  if (fPrintLevel > 0) {
1143  const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mcid);
1144  float KE = ptkl->E() - ptkl->Mass();
1145  std::string KEUnits = " GeV";
1146  if (mctrk.Origin() != simb::kCosmicRay) {
1147  // MeV for low energy particles
1148  KE *= 1000;
1149  KEUnits = " MeV";
1150  }
1151  mf::LogVerbatim("TrackAna")
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;
1157  // this won't work for DUNE
1158  for (unsigned short ipl = 0; ipl < wireReadoutGeom.Nplanes(); ++ipl) {
1159  geo::PlaneID const planeID{0, 0, ipl};
1160  auto const& plane = wireReadoutGeom.Plane(planeID);
1161  auto const sWire = plane.NearestWireID(mcstart).Wire;
1162  auto const sTick = detProp.ConvertXToTicks(mcstart.X(), planeID);
1163  auto const eWire = plane.NearestWireID(mcend).Wire;
1164  auto const eTick = detProp.ConvertXToTicks(mcend.X(), planeID);
1165  mf::LogVerbatim("TrackAna")
1166  << " Wire:Tick in Pln " << ipl << " W:T " << sWire << ":" << sTick
1167  << " - " << eWire << ":" << eTick;
1168  } // ipl
1169  } // fPrintLevel == 2
1170  } // good
1171  }
1172  }
1173  }
1174  }
1175 
1176  // Dump track information here.
1177 
1178  if (pdump) {
1179  TVector3 pos = track.Vertex<TVector3>();
1180  TVector3 dir = track.VertexDirection<TVector3>();
1181  TVector3 end = track.End<TVector3>();
1182  pos[0] += trackdx;
1183  end[0] += trackdx;
1184  TVector3 enddir = track.EndDirection<TVector3>();
1185  double pstart = track.VertexMomentum();
1186  double pend = track.EndMomentum();
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];
1192  if (pstart > 0.) {
1193  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10) << dir[0]
1194  << std::setw(10) << dir[1] << std::setw(10) << dir[2];
1195  }
1196  else
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];
1202  if (pend > 0.01) {
1203  *pdump << " " << std::fixed << std::setprecision(3) << std::setw(10) << enddir[0]
1204  << std::setw(10) << enddir[1] << std::setw(10) << enddir[2];
1205  }
1206  else
1207  *pdump << std::setw(32) << " ";
1208  *pdump << std::setw(12) << std::fixed << std::setprecision(3) << pend
1209  << "\nLength: " << tlen << "\n";
1210  }
1211  }
1212  }
1213  } // i
1214 
1215  // print out un-matched MC tracks
1216  if (fPrintLevel > 0) {
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;
1220  const simb::MCParticle* ptkl = pi_serv->TrackIdToParticle_P(mctrk.TrackID());
1221  float KE = ptkl->E() - ptkl->Mass();
1222  std::string KEUnits = " GeV";
1223  if (mctrk.Origin() != simb::kCosmicRay) {
1224  // MeV for low energy particles
1225  KE *= 1000;
1226  KEUnits = " MeV";
1227  }
1228  // find the start/end wire:time in each plane
1229  geo::Point_t mcstart, mcend;
1230  TVector3 mcstartmom, mcendmom;
1231  double mcdx = mctrk.Start().T() * 1.e-3 * detProp.DriftVelocity(); // cm
1232  double plen = length(detProp, mctrk, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1233  mf::LogVerbatim("TrackAna")
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)
1237  << plen << " cm";
1238  if (fPrintLevel > 1) {
1239  // this won't work for DUNE
1240  for (unsigned short ipl = 0; ipl < wireReadoutGeom.Nplanes(); ++ipl) {
1241  geo::PlaneID const& planeID{0, 0, ipl};
1242  auto const& plane = wireReadoutGeom.Plane(planeID);
1243  auto const sWire = plane.NearestWireID(mcstart).Wire;
1244  auto const sTick = detProp.ConvertXToTicks(mcstart.X(), planeID);
1245  auto const eWire = plane.NearestWireID(mcend).Wire;
1246  auto const eTick = detProp.ConvertXToTicks(mcend.X(), planeID);
1247  mf::LogVerbatim("TrackAna") << " Wire:Tick in Pln " << ipl << " W:T " << sWire << ":"
1248  << sTick << " - " << eWire << ":" << eTick;
1249  } // ipl
1250  } // fPrintLevel > 1
1251  } // imc
1252  } // fPrintLevel > 0
1253  }
1254 
1256  detinfo::DetectorClocksData const& clockData,
1257  detinfo::DetectorPropertiesData const& detProp)
1258  {
1261 
1262  std::map<int, std::map<int, art::PtrVector<recob::Hit>>> hitmap; // trkID, otrk, hitvec
1263  std::map<int, int> KEmap; // length traveled in det [cm]?, trkID want to sort by KE
1264  bool mc = !evt.isRealData();
1268 
1269  evt.getByLabel(fTrackModuleLabel, trackh);
1270  evt.getByLabel(fSpacepointModuleLabel, sppth);
1271  evt.getByLabel(fStitchModuleLabel, trackvh);
1272  int ntv(trackvh->size());
1273 
1274  std::vector<art::PtrVector<recob::Track>>::const_iterator cti = trackvh->begin();
1275 
1276  if (trackh.isValid()) {
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;
1281  }
1282 
1283  if (fRecoHistMap.count(0) == 0) {
1284  fRecoHistMap.emplace(0, RecoHists{"Reco"});
1285  std::cout << "\n"
1286  << "\t\t TrkAna: Fresh fRecoHistMap[0] ******* \n"
1287  << std::endl;
1288  }
1289  const RecoHists& rhistsStitched = fRecoHistMap.at(0);
1290 
1291  std::vector<std::vector<unsigned int>> NtrkIdsAll;
1292  std::vector<double> ntvsorted;
1293  hitmap.clear();
1294  KEmap.clear();
1295 
1296  // Look at the components of the stitched tracks. Grab their sppts/hits from Assns.
1297  for (int o = 0; o < ntv; ++o) // o for outer
1298  {
1299  const art::PtrVector<recob::Track> pvtrack(*(cti++));
1300  int ntrack = pvtrack.size();
1301  std::vector<std::vector<unsigned int>> NtrkId_Hit; // hit IDs in inner tracks
1302  std::vector<unsigned int> vecMode;
1304 
1305  for (int i = 0; i < ntrack; ++i) {
1306 
1307  // From gdb> ptype fs, the vector of Ptr<SpacePoint>s it appears is
1308  // grabbed after fs.at(0)
1309  bool assns(true);
1310  try {
1311  // Got Spacepoints from this Track; now get Hits from those
1312  // Spacepoints.
1313 
1314  int nsppts_assn = fs.at(i).size();
1315 
1316  const auto& sppt = fs.at(i); //.at(0);
1317  // since we're in a try and worried about failure, we won't pull the following
1318  // FindManyP out of the loop.
1320 
1321  // Importantly, loop on all sppts, though they don't all contribute to the track.
1322  // As opposed to looping on the trajectory pts, which is a lower number.
1323  // Also, important, in job in whch this runs I set TrackKal3DSPS parameter MaxPass=1,
1324  // cuz I don't want merely the sparse set of sppts as follows from the uncontained
1325  // MS-measurement in 2nd pass.
1326  std::vector<unsigned int> vecNtrkIds;
1327  for (int is = 0; is < nsppts_assn; ++is) {
1328  int nhits = fh.at(is).size(); // should be 2 or 3: number of planes.
1329  for (int ih = 0; ih < nhits; ++ih) {
1330  const auto& hit = fh.at(is).at(ih); // Our vector is after the .at(is) this time.
1331  if (hit->SignalType() != geo::kCollection) continue;
1332  rhistsStitched.fHHitChg->Fill(hit->Integral());
1333  rhistsStitched.fHHitWidth->Fill(2. * hit->RMS());
1334  if (mc) {
1335  std::vector<sim::TrackIDE> tids = bt_serv->HitToTrackIDEs(clockData, hit);
1336  // more here.
1337  // Loop over track ids.
1338  bool justOne(true); // Only take first trk that contributed to this hit
1339  for (std::vector<sim::TrackIDE>::const_iterator itid = tids.begin();
1340  itid != tids.end();
1341  ++itid) {
1342  int trackID = std::abs(itid->trackID);
1343  hitmap[trackID][o].push_back(hit);
1344 
1345  if (justOne) {
1346  vecNtrkIds.push_back(trackID);
1347  justOne = false;
1348  }
1349  // Add hit to PtrVector corresponding to this track id.
1350  rhistsStitched.fHHitTrkId->Fill(trackID);
1351  const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(trackID);
1352  if (!part) break;
1353 
1354  rhistsStitched.fHHitPdg->Fill(part->PdgCode());
1355  // This really needs to be indexed as KE deposited in volTPC, not just KE. EC, 24-July-2014.
1356 
1357  geo::Point_t mcstart;
1358  geo::Point_t mcend;
1359  TVector3 mcstartmom;
1360  TVector3 mcendmom;
1361  double mctime = part->T(); // nsec
1362  double mcdx = mctime * 1.e-3 * detProp.DriftVelocity(); // cm
1363 
1364  double plen = length(detProp, *part, mcdx, mcstart, mcend, mcstartmom, mcendmom);
1365 
1366  KEmap[(int)(1e6 * plen)] = trackID; // multiple assignment but always the same, so
1367  // fine.
1368  }
1369 
1370  } // mc
1371  } // hits
1372 
1373  } // spacepoints
1374 
1375  if (mc) {
1376  NtrkId_Hit.push_back(vecNtrkIds);
1377  // Find the trkID mode for this i^th track
1378  unsigned int ii(1);
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; }
1384  else {
1385  n++;
1386  }
1387  if (n > max) {
1388  max = n;
1389  ind = ii;
1390  }
1391  ii++;
1392  }
1393  unsigned int mode(sim::NoParticleId);
1394  if (strkIds.begin() != strkIds.end()) mode = strkIds.at(ind);
1395  vecMode.push_back(mode);
1396 
1397  if (strkIds.size() != 0)
1398  rhistsStitched.fModeFrac->Fill((double)max / (double)strkIds.size());
1399  else
1400  rhistsStitched.fModeFrac->Fill(-1.0);
1401  } // mc
1402 
1403  } // end try
1404  catch (cet::exception& x) {
1405  assns = false;
1406  }
1407  if (!assns) throw cet::exception("TrackAna") << "Bad Associations. \n";
1408 
1409  } // i
1410 
1411  if (mc) {
1412  // one vector per o trk, for all modes of stitched i trks
1413  auto& back = NtrkIdsAll.emplace_back(vecMode);
1414 
1415  auto const last = std::unique(back.begin(), back.end());
1416  back.erase(last, back.end());
1417  double sum(0.0);
1418  for (auto const val : back) {
1419  sum += hitmap[val][o].size();
1420  }
1421  ntvsorted.push_back(sum);
1422  }
1423 
1424  //
1425  } // o
1426 
1427  int vtmp(0);
1428  // get KEmap indices by most energetic first, least last.
1429  for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1430  vtmp++;
1431  }
1432 
1433  // get o trk indices by their hits. Most populous track first, least last.
1434  for (auto const o : fsort_indexes(ntvsorted)) {
1435  int v(0);
1436  // get KEmap indices by longest trajectory first, least last.
1437  for (auto it = KEmap.rbegin(); it != KEmap.rend(); ++it) {
1438  int val = it->second; // grab trkIDs in order, since they're sorted by KE
1439  rhistsStitched.fNTrkIdTrks3->Fill(o, v, hitmap[val][o].size());
1440  v++;
1441  }
1442  }
1443 
1444  // In how many o tracks did each trkId appear? Histo it. Would like it to be
1445  // precisely 1. Histo it vs. particle KE.
1446  flattener flat(NtrkIdsAll);
1447  std::vector<unsigned int>& v = flat;
1448  // auto const it ( std::unique(v.begin(),v.end()) ); // never use this it,
1449  // perhaps.
1450  for (auto const val : v) {
1451  if (val != (unsigned int)sim::NoParticleId) {
1452  const simb::MCParticle* part = pi_serv->TrackIdToParticle_P(val);
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);
1456  }
1457  else {
1458  rhistsStitched.fNTrkIdTrks2->Fill(-1.0, 0.0);
1459  }
1460  }
1461  }
1462 
1464  {
1465  // Print summary.
1466 
1467  mf::LogInfo("TrackAna") << "TrackAna statistics:\n"
1468  << " Number of events = " << fNumEvent;
1469 
1470  // Fill efficiency histograms.
1471 
1472  for (std::map<int, MCHists>::const_iterator i = fMCHistMap.begin(); i != fMCHistMap.end();
1473  ++i) {
1474  const MCHists& mchists = i->second;
1475  effcalc(mchists.fHgstartx, mchists.fHmcstartx, mchists.fHestartx);
1476  effcalc(mchists.fHgstarty, mchists.fHmcstarty, mchists.fHestarty);
1477  effcalc(mchists.fHgstartz, mchists.fHmcstartz, mchists.fHestartz);
1478  effcalc(mchists.fHgendx, mchists.fHmcendx, mchists.fHeendx);
1479  effcalc(mchists.fHgendy, mchists.fHmcendy, mchists.fHeendy);
1480  effcalc(mchists.fHgendz, mchists.fHmcendz, mchists.fHeendz);
1481  effcalc(mchists.fHgtheta, mchists.fHmctheta, mchists.fHetheta);
1482  effcalc(mchists.fHgphi, mchists.fHmcphi, mchists.fHephi);
1483  effcalc(mchists.fHgtheta_xz, mchists.fHmctheta_xz, mchists.fHetheta_xz);
1484  effcalc(mchists.fHgtheta_yz, mchists.fHmctheta_yz, mchists.fHetheta_yz);
1485  effcalc(mchists.fHgmom, mchists.fHmcmom, mchists.fHemom);
1486  effcalc(mchists.fHgmoml, mchists.fHmcmoml, mchists.fHemoml);
1487  effcalc(mchists.fHgke, mchists.fHmcke, mchists.fHeke);
1488  effcalc(mchists.fHgkel, mchists.fHmckel, mchists.fHekel);
1489  effcalc(mchists.fHglen, mchists.fHmclen, mchists.fHelen);
1490  effcalc(mchists.fHglens, mchists.fHmclens, mchists.fHelens);
1491  }
1492  }
1493 
1494  // Stole this from online. Returns indices sorted by corresponding vector values.
1495  std::vector<size_t> TrackAna::fsort_indexes(const std::vector<double>& v)
1496  {
1497  // initialize original index locations
1498  std::vector<size_t> idx(v.size());
1499  for (size_t i = 0; i != idx.size(); ++i)
1500  idx[i] = i;
1501  // sort indexes based on comparing values in v
1502  std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) {
1503  return v[i1] > v[i2];
1504  }); // Want most occupied trks first. was <, EC, 21-July.
1505  return idx;
1506  }
1507 
1508 }
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:234
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
simb::Origin_t Origin() const
Definition: MCTrack.h:38
double EndMomentum() const
Definition: Track.h:178
int PdgCode() const
Definition: MCParticle.h:213
double VertexMomentum() const
Definition: Track.h:176
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
Float_t den
Definition: plot.C:35
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
double Mass() const
Definition: MCParticle.h:240
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
enum simb::_ev_origin Origin_t
event origin types
bool HasMomentum() const
Definition: Track.h:174
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
constexpr auto abs(T v)
Returns the absolute value of the argument.
Vector_t VertexDirection() const
Access to track direction at different points.
Definition: Track.h:166
const SMatrixSym55 & VertexCovariance() const
Access to covariance matrices.
Definition: Track.h:191
STL namespace.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::string fStitchModuleLabel
intermediate_table::const_iterator const_iterator
double T() const
Definition: MCStep.h:41
bool isRealData() const
Definition: Event.cc:53
Particle class.
tick ticks
Alias for common language habits.
Definition: electronics.h:76
std::map< int, MCHists > fMCHistMap
3-dimensional objects, potentially hits, clusters, prongs, etc.
Definition: geo_types.h:137
Rotation_t GlobalToLocalRotationAtPoint(size_t p) const
Definition: Track.h:230
bool isValid() const noexcept
Definition: Handle.h:203
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
TString part[npart]
Definition: Style.C:32
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
std::vector< size_t > fsort_indexes(const std::vector< double > &v)
static const int NoParticleId
Definition: sim.h:21
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
Definition: MCTruth.h:26
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.
Definition: Track.h:158
double T(const int i=0) const
Definition: MCParticle.h:225
double ConvertXToTicks(double X, int p, int t, int c) const
const SMatrixSym55 & EndCovariance() const
Access to covariance matrices.
Definition: Track.h:192
std::string fMCTrackModuleLabel
EventNumber_t event() const
Definition: Event.cc:41
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
Class def header for mctrack data container.
double value
Definition: spectrum.C:18
size_type size() const
Definition: PtrVector.h:302
int PdgCode() const
Definition: MCTrack.h:39
void analyze(const art::Event &evt) override
Detector simulation of raw signals on wires.
const TLorentzVector & Momentum() const
Definition: MCStep.h:34
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.
Definition: geo_vectors.h:180
int ID() const
Definition: Track.h:244
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
std::string fOrigin
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.
TDirectory * dir
Definition: macro.C:5
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Vector_t EndDirection() const
Access to track direction at different points.
Definition: Track.h:167
std::map< int, RecoHists > fRecoHistMap
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
void endJob() override
Supernova neutrinos.
Definition: MCTruth.h:25
double E() const
Definition: MCStep.h:45
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:159
const MCStep & Start() const
Definition: MCTrack.h:42
Char_t n[5]
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:453
std::string fTrackModuleLabel
unsigned int TrackID() const
Definition: MCTrack.h:40
TCEvent evt
Definition: DataStructs.cxx:8
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
Float_t track
Definition: plot.C:35
Float_t w
Definition: plot.C:20
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:
Definition: Track.h:49
vec_iX clear()
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Cosmic rays.
Definition: MCTruth.h:24
Signal from collection planes.
Definition: geo_types.h:148
Beam neutrinos.
Definition: MCTruth.h:23