SimulationDrawer.cxx
Go to the documentation of this file.
1 ///
2 /// \file SimulationDrawer.cxx
3 /// \brief Render the objects from the Simulation package
4 /// \author messier@indiana.edu
5 ///
6 
7 #include <algorithm>
8 #include <iomanip>
9 
10 #include "TDatabasePDG.h"
11 #include "TLatex.h"
12 #include "TLine.h"
13 #include "TPolyLine.h"
14 #include "TPolyLine3D.h"
15 #include "TPolyMarker.h"
16 #include "TPolyMarker3D.h"
17 
33 #include "nuevdb/EventDisplayBase/View2D.h"
34 #include "nuevdb/EventDisplayBase/View3D.h"
37 
42 
43 namespace {
44  // Utility function to make uniform error messages.
45  void
46  writeErrMsg(const char* fcn, cet::exception const& e)
47  {
48  mf::LogWarning("SimulationDrawer") << "SimulationDrawer::" << fcn << " failed with message:\n"
49  << e;
50  }
51 }
52 
53 namespace evd {
54 
56  {
57  // For now only draw cryostat=0.
59  minx = 1e9;
60  maxx = -1e9;
61  miny = 1e9;
62  maxy = -1e9;
63  minz = 1e9;
64  maxz = -1e9;
65 
66  // This is looking to find the range of the complete active volume... this may not be the
67  // best way to do this...
68  for (size_t cryoIdx = 0; cryoIdx < geom->Ncryostats(); cryoIdx++) {
69  const geo::CryostatGeo& cryoGeo = geom->Cryostat(cryoIdx);
70 
71  for (size_t tpcIdx = 0; tpcIdx < cryoGeo.NTPC(); tpcIdx++) {
72  const geo::TPCGeo& tpc = cryoGeo.TPC(tpcIdx);
73 
74  mf::LogDebug("SimulationDrawer")
75  << "Cryo/TPC idx: " << cryoIdx << "/" << tpcIdx << ", TPC center: " << tpc.GetCenter()[0]
76  << ", " << tpc.GetCenter()[1] << ", " << tpc.GetCenter()[2] << std::endl;
77  mf::LogDebug("SimulationDrawer")
78  << " TPC Active center: " << tpc.GetActiveVolumeCenter()[0] << ", "
79  << tpc.GetActiveVolumeCenter()[1] << ", " << tpc.GetActiveVolumeCenter()[2]
80  << ", H/W/L: " << tpc.ActiveHalfHeight() << "/" << tpc.ActiveHalfWidth() << "/"
81  << tpc.ActiveLength() << std::endl;
82 
83  if (minx > tpc.GetCenter()[0] - tpc.HalfWidth())
84  minx = tpc.GetCenter()[0] - tpc.HalfWidth();
85  if (maxx < tpc.GetCenter()[0] + tpc.HalfWidth())
86  maxx = tpc.GetCenter()[0] + tpc.HalfWidth();
87  if (miny > tpc.GetCenter()[1] - tpc.HalfHeight())
88  miny = tpc.GetCenter()[1] - tpc.HalfHeight();
89  if (maxy < tpc.GetCenter()[1] + tpc.HalfHeight())
90  maxy = tpc.GetCenter()[1] + tpc.HalfHeight();
91  if (minz > tpc.GetCenter()[2] - tpc.Length() / 2.)
92  minz = tpc.GetCenter()[2] - tpc.Length() / 2.;
93  if (maxz < tpc.GetCenter()[2] + tpc.Length() / 2.)
94  maxz = tpc.GetCenter()[2] + tpc.Length() / 2.;
95 
96  mf::LogDebug("SimulationDrawer")
97  << " minx/maxx: " << minx << "/" << maxx << ", miny/maxy: " << miny << "/" << maxy
98  << ", minz/miny: " << minz << "/" << maxz << std::endl;
99  }
100  }
101  }
102 
103  //......................................................................
104 
106 
107  //......................................................................
108 
109  void
111  {
112 
113  if (evt.isRealData()) return;
114 
116  // Skip drawing if option is turned off
117  if (!drawopt->fShowMCTruthText) return;
118 
119  std::vector<const simb::MCTruth*> mctruth;
120  this->GetMCTruth(evt, mctruth);
121 
122  for (unsigned int i = 0; i < mctruth.size(); ++i) {
123  std::string mctext;
124  bool firstin = true;
125  bool firstout = true;
127  std::string incoming;
128  std::string outgoing;
129  // Label cosmic rays -- others are pretty obvious
130  if (mctruth[i]->Origin() == simb::kCosmicRay) origin = "c-ray: ";
131  int jmax = TMath::Min(20, mctruth[i]->NParticles());
132  for (int j = 0; j < jmax; ++j) {
133  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
134  char buff[1024];
135  if (p.P() > 0.05) {
136  sprintf(buff,
137  "#color[%d]{%s #scale[0.75]{[%.1f GeV/c]}}",
140  p.P());
141  }
142  else {
143  sprintf(buff,
144  "#color[%d]{%s}",
147  }
148  if (p.StatusCode() == 0) {
149  if (firstin == false) incoming += " + ";
150  incoming += buff;
151  firstin = false;
152  }
153  if (p.StatusCode() == 1) {
154  if (firstout == false) outgoing += " + ";
155  outgoing += buff;
156  firstout = false;
157  }
158  } // loop on j particles
159  if (origin == "" && incoming == "") { mctext = outgoing; }
160  else {
161  mctext = origin + incoming + " #rightarrow " + outgoing;
162  }
163  TLatex& latex = view->AddLatex(0.03, 0.2, mctext.c_str());
164  latex.SetTextSize(0.6);
165 
166  } // loop on i mctruth objects
167  }
168 
169  //......................................................................
170 
171  void
172  SimulationDrawer::MCTruthLongText(const art::Event& evt, evdb::View2D* /*view*/)
173  {
174  if (evt.isRealData()) return;
175 
177  // Skip drawing if option is turned off
178  if (!drawopt->fShowMCTruthText) return;
179 
180  std::vector<const simb::MCTruth*> mctruth;
181  this->GetMCTruth(evt, mctruth);
182  std::cout << "\nMCTruth Ptcl trackID PDG P T Moth Process\n";
183  for (unsigned int i = 0; i < mctruth.size(); ++i) {
184  for (int j = 0; j < mctruth[i]->NParticles(); ++j) {
185  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
186  if (p.StatusCode() == 0 || p.StatusCode() == 1) {
187  int KE = 1000 * (p.E() - p.Mass());
188  std::cout << std::right << std::setw(7) << i << std::setw(5) << j << std::setw(8)
189  << p.TrackId() << " " << std::setw(14) << Style::LatexName(p.PdgCode())
190  << std::setw(7) << int(1000 * p.P()) << std::setw(7) << KE << std::setw(7)
191  << p.Mother() << " " << p.Process() << "\n";
192  }
193  /*
194  std::cout << Style::LatexName(p.PdgCode())
195  << "\t\t" << p.E() << " GeV"
196  << "\t" << "(" << p.P() << " GeV/c)"
197  << std::endl;
198 */
199  } // loop on j particles in list
200  }
201  std::cout << "Note: Momentum, P, and kinetic energy, T, in MeV/c\n";
202  } // MCTruthLongText
203 
204  //......................................................................
205  //this is the method you would use to color code hits by the MC truth pdg value
206  void
207  SimulationDrawer::MCTruthVectors2D(const art::Event& evt, evdb::View2D* view, unsigned int plane)
208  {
209  if (evt.isRealData()) return;
210 
211  const spacecharge::SpaceCharge* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
212 
214  // If the option is turned off, there's nothing to do
215  if (drawopt->fShowMCTruthVectors > 3) {
216  std::cout << "Unsupported ShowMCTruthVectors option (> 2)\n";
217  return;
218  }
219 
222  // get the x position of the plane in question
223  double xyz1[3] = {0.};
224  double xyz2[3] = {0.};
225 
226  // shift the truth by a fixed amount so it doesn't overlay the reco
227  double xShift = -2;
228 
229  static bool first = true;
230  if (first) {
231  std::cout
232  << "******** Show MCTruth (Genie) particles when ShowMCTruthVectors = 1 or 3 ******** \n";
233  std::cout << " MCTruth vectors corrected for space charge? " << sce->EnableCorrSCE()
234  << " and shifted by " << xShift << " cm in X\n";
235  std::cout << " Neutrons and photons drawn with dotted lines. \n";
236  std::cout << " Red = e+/-, nue, nuebar. Blue = mu+/-, numu, numubar. Green = tau+/-, nutau, "
237  "nutaubar.\n";
238  std::cout << " Yellow = photons. Magenta = pions, protons and nuetrons.\n";
239  std::cout << "******** Show MCParticle (Geant) decay photons (e.g. from pizeros) when "
240  "ShowMCTruthVectors = 2 or 3 ******** \n";
241  std::cout << " Photons > 50 MeV are drawn as dotted lines corrected for space charge and "
242  "are not shifted.\n";
243  std::cout << " Decay photon end points are drawn at 2 interaction lengths (44 cm) from the "
244  "start position.\n";
245  std::cout << " Color: Green = (50 < E_photon < 100 MeV), Blue = (100 MeV < E_photon < 200 "
246  "MeV), Red = (E_photon > 300 MeV).\n";
247  }
248 
249  bool showTruth = (drawopt->fShowMCTruthVectors == 1 || drawopt->fShowMCTruthVectors == 3);
250  bool showPhotons = (drawopt->fShowMCTruthVectors > 1);
251 
252  auto const detProp =
254 
255  if (showTruth) {
256  // Unpack and draw the MC vectors
257  std::vector<const simb::MCTruth*> mctruth;
258  this->GetMCTruth(evt, mctruth);
259 
260  for (size_t i = 0; i < mctruth.size(); ++i) {
261  if (mctruth[i]->Origin() == simb::kCosmicRay) continue;
262  for (int j = 0; j < mctruth[i]->NParticles(); ++j) {
263  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
264 
265  // Skip all but incoming and out-going particles
266  if (!(p.StatusCode() == 0 || p.StatusCode() == 1)) continue;
267 
268  double r = p.P() * 10.0; // Scale length so 10 cm = 1 GeV/c
269 
270  if (p.StatusCode() == 0) r = -r; // Flip for incoming particles
271 
272  auto gptStart = geo::Point_t(p.Vx(), p.Vy(), p.Vz());
273  geo::Point_t sceOffset{0, 0, 0};
274  if (sce->EnableCorrSCE()) sceOffset = sce->GetPosOffsets(gptStart);
275  xyz1[0] = p.Vx() - sceOffset.X();
276  xyz1[1] = p.Vy() + sceOffset.Y();
277  xyz1[2] = p.Vz() + sceOffset.Z();
278  xyz2[0] = xyz1[0] + r * p.Px() / p.P();
279  xyz2[1] = xyz1[1] + r * p.Py() / p.P();
280  xyz2[2] = xyz1[2] + r * p.Pz() / p.P();
281 
282  double w1 =
283  geo->WireCoordinate(xyz1[1], xyz1[2], (int)plane, rawopt->fTPC, rawopt->fCryostat);
284  double w2 =
285  geo->WireCoordinate(xyz2[1], xyz2[2], (int)plane, rawopt->fTPC, rawopt->fCryostat);
286 
287  double time =
288  detProp.ConvertXToTicks(xyz1[0] + xShift, (int)plane, rawopt->fTPC, rawopt->fCryostat);
289  double time2 =
290  detProp.ConvertXToTicks(xyz2[0] + xShift, (int)plane, rawopt->fTPC, rawopt->fCryostat);
291 
292  if (rawopt->fAxisOrientation < 1) {
293  TLine& l = view->AddLine(w1, time, w2, time2);
295  }
296  else {
297  TLine& l = view->AddLine(time, w1, time2, w2);
299  }
300  //
301 
302  } // loop on j particles in list
303  } // loop on truths
304  } // showTruth
305 
306  if (showPhotons) {
307  // draw pizero photons with T > 30 MeV
309  sim::ParticleList const& plist = pi_serv->ParticleList();
310  if (plist.empty()) return;
311  // photon interaction length approximate
312  double r = 44;
313  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
314  simb::MCParticle* p = (*ipart).second;
315  int trackID = p->TrackId();
316  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
317  if (theTruth->Origin() == simb::kCosmicRay) continue;
318  if (p->PdgCode() != 22) continue;
319  if (p->Process() != "Decay") continue;
320  int TMeV = 1000 * (p->E() - p->Mass());
321  if (TMeV < 30) continue;
322  auto gptStart = geo::Point_t(p->Vx(), p->Vy(), p->Vz());
323  geo::Point_t sceOffset{0, 0, 0};
324  if (sce->EnableCorrSCE()) sceOffset = sce->GetPosOffsets(gptStart);
325  xyz1[0] = p->Vx() - sceOffset.X();
326  xyz1[1] = p->Vy() + sceOffset.Y();
327  xyz1[2] = p->Vz() + sceOffset.Z();
328  xyz2[0] = xyz1[0] + r * p->Px() / p->P();
329  xyz2[1] = xyz1[1] + r * p->Py() / p->P();
330  xyz2[2] = xyz1[2] + r * p->Pz() / p->P();
331  double w1 =
332  geo->WireCoordinate(xyz1[1], xyz1[2], (int)plane, rawopt->fTPC, rawopt->fCryostat);
333  double t1 = detProp.ConvertXToTicks(xyz1[0], (int)plane, rawopt->fTPC, rawopt->fCryostat);
334  double w2 =
335  geo->WireCoordinate(xyz2[1], xyz2[2], (int)plane, rawopt->fTPC, rawopt->fCryostat);
336  double t2 = detProp.ConvertXToTicks(xyz2[0], (int)plane, rawopt->fTPC, rawopt->fCryostat);
337  TLine& l = view->AddLine(w1, t1, w2, t2);
338  l.SetLineWidth(2);
339  l.SetLineStyle(kDotted);
340  if (TMeV < 100) { l.SetLineColor(kGreen); }
341  else if (TMeV < 200) {
342  l.SetLineColor(kBlue);
343  }
344  else {
345  l.SetLineColor(kRed);
346  }
347  } // ipart
348  } // showPhotons
349 
350  first = false;
351 
352  } // MCTruthVectors2D
353 
354  //......................................................................
355  //this method draws the true particle trajectories in 3D
356  void
357  SimulationDrawer::MCTruth3D(const art::Event& evt, evdb::View3D* view)
358  {
359  if (evt.isRealData()) return;
360 
362  // If the option is turned off, there's nothing to do
363  if (!drawopt->fShowMCTruthTrajectories) return;
364 
365  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
366  auto const detProp =
369 
370  // get the particles from the Geant4 step
371  std::vector<const simb::MCParticle*> plist;
372  this->GetParticle(evt, plist);
373 
374  // Define a couple of colors for neutrals and if we gray it out...
375  int neutralColor(12);
376  int grayedColor(15);
377  int neutrinoColor(38);
378 
379  // Use the LArVoxelList to get the true energy deposition locations as opposed to using MCTrajectories
380  const sim::LArVoxelList voxels =
382 
383  mf::LogDebug("SimulationDrawer")
384  << "Starting loop over " << plist.size() << " McParticles, voxel list size is "
385  << voxels.size() << std::endl;
386 
387  // Using the voxel information can be slow (see previous implementation of this code).
388  // In order to speed things up we have modified the strategy:
389  // 1) Make one pass through the list of voxels
390  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
391  // which is done by keeping a map between the MCParticle and a vector of positions
392  // 3) Then loop through the map to draw the particle trajectories.
393  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
394  // more loop to make a map of track id's and MCParticles.
395 
396  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
397  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
398 
399  // Should we display the trajectories too?
400  double minPartEnergy(0.01);
401 
402  for (size_t p = 0; p < plist.size(); ++p) {
403  trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
404 
405  // Quick loop through to draw trajectories...
406  if (drawopt->fShowMCTruthTrajectories) {
407  // Is there an associated McTrajectory?
408  const simb::MCParticle* mcPart = plist[p];
409  const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
410 
411  int pdgCode(mcPart->PdgCode());
412  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
413  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
414  double partCharge = partPDG ? partPDG->Charge() : 0.;
415  double partEnergy = mcPart->E();
416 
417  if (!drawopt->fShowMCTruthColors) colorIdx = grayedColor;
418 
419  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000) {
420  double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
421  double xOffset = 0.;
422  double xPosMinTick = 0.;
423  double xPosMaxTick = std::numeric_limits<double>::max();
424 
425  // collect the points from this particle
426  int numTrajPoints = mcTraj.size();
427 
428  std::unique_ptr<double[]> hitPositions(new double[3 * numTrajPoints]);
429  int hitCount(0);
430 
431  for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
432  double xPos = mcTraj.X(hitIdx);
433  double yPos = mcTraj.Y(hitIdx);
434  double zPos = mcTraj.Z(hitIdx);
435 
436  // If we have cosmic rays then we need to get the offset which allows translating from
437  // when they were generated vs when they were tracked.
438  // Note that this also explicitly checks that they are in a TPC volume
439  geo::Point_t hitLocation(xPos, yPos, zPos);
440 
441  try {
442  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
443  geo::PlaneID planeID(tpcID, 0);
444 
445  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
446  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
447  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
448 
449  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
450  }
451  catch (...) {
452  continue;
453  }
454 
455  // Now move the hit position to correspond to the timing
456  xPos += xOffset;
457 
458  // Check fiducial limits
459  if (xPos > xPosMinTick && xPos < xPosMaxTick) {
460  // Check for space charge offsets
461  // if (spaceCharge->EnableSimEfieldSCE())
462  // {
463  // std::vector<double> offsetVec = spaceCharge->GetPosOffsets(xPos,yPos,zPos);
464  // xPos += offsetVec[0] - 0.7;
465  // yPos -= offsetVec[1];
466  // zPos -= offsetVec[2];
467  // }
468 
469  hitPositions[3 * hitCount] = xPos;
470  hitPositions[3 * hitCount + 1] = yPos;
471  hitPositions[3 * hitCount + 2] = zPos;
472  hitCount++;
473  }
474  }
475 
476  TPolyLine3D& pl(view->AddPolyLine3D(1, colorIdx, 1, 1));
477 
478  // Draw neutrals as a gray dotted line to help fade into background a bit...
479  if (partCharge == 0.) {
480  pl.SetLineColor(neutralColor);
481  pl.SetLineStyle(3);
482  pl.SetLineWidth(1);
483  }
484  pl.SetPolyLine(hitCount, hitPositions.get(), "");
485  }
486  }
487  }
488 
489  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
490  std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
491 
493  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
494  const sim::LArVoxelData& vxd = (*vxitr).second;
495 
496  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
497  if (vxd.Energy(partIdx) > drawopt->fMinEnergyDeposition) {
498  int trackId = vxd.TrackID(partIdx);
499 
500  // It can be in some instances that mcPart here could be zero.
501  const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
502 
503  partToPosMap[mcPart].push_back(std::vector<double>(3));
504 
505  partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
506  partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
507  partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
508  }
509  } // end if this track id is in the current voxel
510  } // end loop over voxels
511 
512  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
513  // draw the trajectories
514  std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
515 
516  for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
517  partToPosMapItr++) {
518  // Recover the McParticle, we'll need to access several data members so may as well dereference it
519  const simb::MCParticle* mcPart = partToPosMapItr->first;
520 
521  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
522  if (!mcPart || partToPosMapItr->second.empty()) continue;
523 
524  double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
525  double xOffset = 0.;
526  double xPosMinTick = 0.;
527  double xPosMaxTick = std::numeric_limits<double>::max();
528 
529  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
530  int markerIdx(kFullDotSmall);
531  int markerSize(2);
532 
533  if (!drawopt->fShowMCTruthFullSize) {
534  colorIdx = grayedColor;
535  markerIdx = kDot;
536  markerSize = 1;
537  }
538 
539  std::unique_ptr<double[]> hitPositions(new double[3 * partToPosMapItr->second.size()]);
540  int hitCount(0);
541 
542  // Now loop over points and add to trajectory
543  for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
544  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
545 
546  // Check xOffset state and set if necessary
547  geo::Point_t hitLocation(posVec[0], posVec[1], posVec[2]);
548 
549  try {
550  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
551  geo::PlaneID planeID(tpcID, 0);
552 
553  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
554  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
555  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
556 
557  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
558  }
559  catch (...) {
560  continue;
561  }
562 
563  double xCoord = posVec[0] + xOffset;
564 
565  // If a voxel records an energy deposit then must have been in the TPC
566  // But because things get shifted still need to cut off if outside drift
567  if (xCoord > xPosMinTick && xCoord < xPosMaxTick) {
568  hitPositions[3 * hitCount] = xCoord;
569  hitPositions[3 * hitCount + 1] = posVec[1];
570  hitPositions[3 * hitCount + 2] = posVec[2];
571  hitCount++;
572  }
573  }
574 
575  TPolyMarker3D& pm = view->AddPolyMarker3D(1, colorIdx, markerIdx, markerSize);
576  pm.SetPolyMarker(hitCount, hitPositions.get(), markerIdx);
577  }
578 
579  // Finally, let's see if we can draw the incoming particle from the MCTruth information
580  std::vector<const simb::MCTruth*> mctruth;
581  this->GetMCTruth(evt, mctruth);
582 
583  // Loop through the MCTruth vector
584  for (unsigned int idx = 0; idx < mctruth.size(); idx++) {
585  // Go through each MCTruth object in the list
586  for (int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++) {
587  const simb::MCParticle& mcPart = mctruth[idx]->GetParticle(particleIdx);
588 
589  // A negative mother id indicates the "primary" particle
590  if (mcPart.Mother() == -1 && mcPart.StatusCode() == 0) {
591  mf::LogDebug("SimulationDrawer") << mcPart << std::endl;
592 
593  // Get position vector
594  TVector3 particlePosition(mcPart.Vx(), mcPart.Vy(), mcPart.Vz());
595 
596  // Get direction vector (in opposite direction)
597  TVector3 oppPartDir(-mcPart.Px(), -mcPart.Py(), -mcPart.Pz());
598 
599  if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
600 
601  double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
602 
603  // No point in drawing if arc length is zero (e.g. Ar nucleus)
604  if (arcLenToDraw > 0.) {
605  // Draw the line, use an off color to be unique
606  TPolyLine3D& pl(view->AddPolyLine3D(2, neutrinoColor, 1, 2));
607 
608  pl.SetPoint(0, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
609 
610  particlePosition += std::min(arcLenToDraw + 10., 1000.) * oppPartDir;
611 
612  pl.SetPoint(1, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
613  }
614  }
615  // The particles we want to draw will be early in the list so break out if we didn't find them
616  else
617  break;
618  } // loop on particles in list
619  }
620 
621  return;
622  }
623 
624  //......................................................................
625  //this method draws the true particle trajectories in 3D Ortho view.
626  void
628  evd::OrthoProj_t proj,
629  double msize,
630  evdb::View2D* view)
631  {
632  if (evt.isRealData()) return;
633 
635 
636  // If the option is turned off, there's nothing to do
637  if (!drawopt->fShowMCTruthTrajectories) return;
638 
639  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
640  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
641  auto const detProp =
643 
644  // get the particles from the Geant4 step
645  std::vector<const simb::MCParticle*> plist;
646  this->GetParticle(evt, plist);
647 
648  // Useful variables
649 
650  double xMinimum(-1. * (maxx - minx));
651  double xMaximum(2. * (maxx - minx));
652 
653  // Use the LArVoxelList to get the true energy deposition locations as
654  // opposed to using MCTrajectories
655  // GetLarVoxelList should be adjusted to take an InputTag instead of strange.
656  const sim::LArVoxelList voxels =
658 
659  mf::LogDebug("SimulationDrawer")
660  << "Starting loop over " << plist.size() << " McParticles, voxel list size is "
661  << voxels.size() << std::endl;
662 
663  // Using the voxel information can be slow (see previous implementation of this code).
664  // In order to speed things up we have modified the strategy:
665  // 1) Make one pass through the list of voxels
666  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
667  // which is done by keeping a map between the MCParticle and a vector of positions
668  // 3) Then loop through the map to draw the particle trajectories.
669  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
670  // more loop to make a map of track id's and MCParticles.
671 
672  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
673  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
674 
675  // Should we display the trajectories too?
676  bool displayMcTrajectories(true);
677  double minPartEnergy(0.025);
678 
679  double tpcminx = 1.0;
680  double tpcmaxx = -1.0;
681  double xOffset = 0.0;
682  double g4Ticks = 0.0;
683  double coeff = 0.0;
684  double readoutwindowsize = 0.0;
685  double vtx[3] = {0.0, 0.0, 0.0};
686  for (size_t p = 0; p < plist.size(); ++p) {
687  trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
688 
689  // Quick loop through to drawn trajectories...
690  if (displayMcTrajectories) {
691  // Is there an associated McTrajectory?
692  const simb::MCParticle* mcPart = plist[p];
693  const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
694 
695  int pdgCode(mcPart->PdgCode());
696  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
697  double partCharge = partPDG ? partPDG->Charge() : 0.;
698  double partEnergy = mcPart->E();
699 
700  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000) {
701  // collect the points from this particle
702  int numTrajPoints = mcTraj.size();
703 
704  std::unique_ptr<double[]> hitPosX(new double[numTrajPoints]);
705  std::unique_ptr<double[]> hitPosY(new double[numTrajPoints]);
706  std::unique_ptr<double[]> hitPosZ(new double[numTrajPoints]);
707  int hitCount(0);
708 
709  double xPos = mcTraj.X(0);
710  double yPos = mcTraj.Y(0);
711  double zPos = mcTraj.Z(0);
712 
713  tpcminx = 1.0;
714  tpcmaxx = -1.0;
715  xOffset = 0.0;
716  g4Ticks = 0.0;
717  vtx[0] = 0.0;
718  vtx[1] = 0.0;
719  vtx[2] = 0.0;
720  coeff = 0.0;
721  readoutwindowsize = 0.0;
722  for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
723  xPos = mcTraj.X(hitIdx);
724  yPos = mcTraj.Y(hitIdx);
725  zPos = mcTraj.Z(hitIdx);
726 
727  // If the original simulated hit did not occur in the TPC volume then don't draw it
728  if (xPos < minx || xPos > maxx || yPos < miny || yPos > maxy || zPos < minz ||
729  zPos > maxz)
730  continue;
731 
732  if ((xPos < tpcminx) || (xPos > tpcmaxx)) {
733  vtx[0] = xPos;
734  vtx[1] = yPos;
735  vtx[2] = zPos;
736  geo::TPCID tpcid = geom->FindTPCAtPosition(vtx);
737  unsigned int cryo = geom->FindCryostatAtPosition(vtx);
738 
739  if (tpcid.isValid) {
740  unsigned int tpc = tpcid.TPC;
741  const geo::TPCGeo& tpcgeo = geom->GetElement(tpcid);
742  tpcminx = tpcgeo.MinX();
743  tpcmaxx = tpcgeo.MaxX();
744 
745  coeff = detProp.GetXTicksCoefficient(tpc, cryo);
746  readoutwindowsize =
747  detProp.ConvertTicksToX(detProp.ReadOutWindowSize(), 0, tpc, cryo);
748 
749  // The following is meant to get the correct offset for drawing
750  // the particle trajectory In particular, the cosmic rays will
751  // not be correctly placed without this
752  g4Ticks = clockData.TPCG4Time2Tick(mcPart->T()) +
753  detProp.GetXTicksOffset(0, tpc, cryo) - trigger_offset(clockData);
754 
755  xOffset = detProp.ConvertTicksToX(g4Ticks, 0, tpc, cryo);
756  }
757  else {
758  xOffset = 0;
759  tpcminx = 1.0;
760  tpcmaxx = -1.0;
761  coeff = 0.0;
762  readoutwindowsize = 0.0;
763  }
764  }
765 
766  // Now move the hit position to correspond to the timing
767  xPos += xOffset;
768 
769  bool inreadoutwindow = false;
770  if (coeff < 0) {
771  if ((xPos > readoutwindowsize) && (xPos < tpcmaxx)) inreadoutwindow = true;
772  }
773  else if (coeff > 0) {
774  if ((xPos > tpcminx) && (xPos < readoutwindowsize)) inreadoutwindow = true;
775  }
776 
777  if (!inreadoutwindow) continue;
778 
779  // Check fiducial limits
780  if (xPos > xMinimum && xPos < xMaximum) {
781  hitPosX[hitCount] = xPos;
782  hitPosY[hitCount] = yPos;
783  hitPosZ[hitCount] = zPos;
784  hitCount++;
785  }
786  }
787 
788  TPolyLine& pl = view->AddPolyLine(
789  1, evd::Style::ColorFromPDG(mcPart->PdgCode()), 1, 1); //kFullCircle, msize);
790 
791  // Draw neutrals as a gray dotted line to help fade into background a bit...
792  if (partCharge == 0.) {
793  pl.SetLineColor(13);
794  pl.SetLineStyle(3);
795  pl.SetLineWidth(1);
796  }
797 
798  if (proj == evd::kXY)
799  pl.SetPolyLine(hitCount, hitPosX.get(), hitPosY.get(), "");
800  else if (proj == evd::kXZ)
801  pl.SetPolyLine(hitCount, hitPosZ.get(), hitPosX.get(), "");
802  else if (proj == evd::kYZ)
803  pl.SetPolyLine(hitCount, hitPosZ.get(), hitPosY.get(), "");
804  }
805  }
806  }
807 
808  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
809  std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
810 
812  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
813  const sim::LArVoxelData& vxd = (*vxitr).second;
814 
815  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
816  if (vxd.Energy(partIdx) > drawopt->fMinEnergyDeposition) {
817  int trackId = vxd.TrackID(partIdx);
818 
819  // It can be in some instances that mcPart here could be zero.
820  const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
821 
822  partToPosMap[mcPart].push_back(std::vector<double>(3));
823 
824  partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
825  partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
826  partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
827  }
828  } // end if this track id is in the current voxel
829  } // end loop over voxels
830 
831  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
832  // draw the trajectories
833  std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
834 
835  for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
836  partToPosMapItr++) {
837  // Recover the McParticle, we'll need to access several data members so may as well dereference it
838  const simb::MCParticle* mcPart = partToPosMapItr->first;
839 
840  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
841  if (!mcPart || partToPosMapItr->second.empty()) continue;
842 
843  tpcminx = 1.0;
844  tpcmaxx = -1.0;
845  xOffset = 0.0;
846  g4Ticks = 0.0;
847  std::vector<std::array<double, 3>> posVecCorr;
848  posVecCorr.reserve(partToPosMapItr->second.size());
849  coeff = 0.0;
850  readoutwindowsize = 0.0;
851 
852  // Now loop over points and add to trajectory
853  for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
854  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
855 
856  if ((posVec[0] < tpcminx) || (posVec[0] > tpcmaxx)) {
857  vtx[0] = posVec[0];
858  vtx[1] = posVec[1];
859  vtx[2] = posVec[2];
860  geo::TPCID tpcid = geom->FindTPCAtPosition(vtx);
861  unsigned int cryo = geom->FindCryostatAtPosition(vtx);
862 
863  if (tpcid.isValid) {
864  unsigned int tpc = tpcid.TPC;
865 
866  const geo::TPCGeo& tpcgeo = geom->GetElement(tpcid);
867  tpcminx = tpcgeo.MinX();
868  tpcmaxx = tpcgeo.MaxX();
869 
870  coeff = detProp.GetXTicksCoefficient(tpc, cryo);
871  readoutwindowsize = detProp.ConvertTicksToX(detProp.ReadOutWindowSize(), 0, tpc, cryo);
872  // The following is meant to get the correct offset for drawing the
873  // particle trajectory In particular, the cosmic rays will not be
874  // correctly placed without this
875  g4Ticks = clockData.TPCG4Time2Tick(mcPart->T()) +
876  detProp.GetXTicksOffset(0, tpc, cryo) - trigger_offset(clockData);
877 
878  xOffset = detProp.ConvertTicksToX(g4Ticks, 0, tpc, cryo);
879  }
880  else {
881  xOffset = 0;
882  tpcminx = 1.0;
883  tpcmaxx = -1.0;
884  coeff = 0.0;
885  readoutwindowsize = 0.0;
886  }
887  }
888 
889  double xCoord = posVec[0] + xOffset;
890 
891  bool inreadoutwindow = false;
892  if (coeff < 0) {
893  if ((xCoord > readoutwindowsize) && (xCoord < tpcmaxx)) inreadoutwindow = true;
894  }
895  else if (coeff > 0) {
896  if ((xCoord > tpcminx) && (xCoord < readoutwindowsize)) inreadoutwindow = true;
897  }
898 
899  if (inreadoutwindow && (xCoord > xMinimum && xCoord < xMaximum)) {
900  posVecCorr.push_back({{xCoord, posVec[1], posVec[2]}});
901  }
902  }
903 
904  TPolyMarker& pm = view->AddPolyMarker(posVecCorr.size(),
906  kFullDotMedium,
907  2); //kFullCircle, msize);
908 
909  for (size_t p = 0; p < posVecCorr.size(); ++p) {
910  if (proj == evd::kXY)
911  pm.SetPoint(p, posVecCorr[p][0], posVecCorr[p][1]);
912  else if (proj == evd::kXZ)
913  pm.SetPoint(p, posVecCorr[p][2], posVecCorr[p][0]);
914  else if (proj == evd::kYZ)
915  pm.SetPoint(p, posVecCorr[p][2], posVecCorr[p][1]);
916  }
917  }
918 
919  return;
920  }
921 
922  //......................................................................
923  int
924  SimulationDrawer::GetParticle(const art::Event& evt, std::vector<const simb::MCParticle*>& plist)
925  {
926  plist.clear();
927 
928  if (evt.isRealData()) return 0;
929 
931 
932  std::vector<const simb::MCParticle*> temp;
933 
935  // use get by Type because there should only be one collection of these in the event
936  try {
937  evt.getView(drawopt->fG4ModuleLabel, plcol);
938  for (unsigned int i = 0; i < plcol.vals().size(); ++i) {
939  temp.push_back(plcol.vals().at(i));
940  }
941  temp.swap(plist);
942  }
943  catch (cet::exception& e) {
944  writeErrMsg("GetRawDigits", e);
945  }
946 
947  return plist.size();
948  }
949 
950  //......................................................................
951 
952  int
953  SimulationDrawer::GetMCTruth(const art::Event& evt, std::vector<const simb::MCTruth*>& mcvec)
954  {
955  mcvec.clear();
956 
957  if (evt.isRealData()) return 0;
958 
959  std::vector<const simb::MCTruth*> temp;
960 
961  std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
962  // use get by Type because there should only be one collection of these in the event
963  try {
964  //evt.getManyByType(mctcol);
965  mctcol = evt.getMany<std::vector<simb::MCTruth>>();
966  for (size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
967  art::Handle<std::vector<simb::MCTruth>> mclistHandle = mctcol[mctc];
968 
969  for (size_t i = 0; i < mclistHandle->size(); ++i) {
970  temp.push_back(&(mclistHandle->at(i)));
971  }
972  }
973  temp.swap(mcvec);
974  }
975  catch (cet::exception& e) {
976  writeErrMsg("GetMCTruth", e);
977  }
978 
979  return mcvec.size();
980  }
981 
982  //......................................................................
983 
984  void
985  SimulationDrawer::HiLite(int trkId, bool dohilite)
986  {
987  fHighlite[trkId] = dohilite;
988  }
989 
990 } // namespace
991 ////////////////////////////////////////////////////////////////////////
std::map< int, bool > fHighlite
double E(const int i=0) const
Definition: MCParticle.h:233
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
double Z(const size_type i) const
Definition: MCTrajectory.h:151
double X(const size_type i) const
Definition: MCTrajectory.h:149
int PdgCode() const
Definition: MCParticle.h:212
Point GetActiveVolumeCenter() const
Returns the center of the TPC active volume in world coordinates [cm].
Definition: TPCGeo.h:288
double Py(const int i=0) const
Definition: MCParticle.h:231
void MCTruthLongText(const art::Event &evt, evdb::View2D *view)
CryostatGeo const & GetElement(geo::CryostatID const &cryoid) const
unsigned int fTPC
TPC number to draw, typically set by TWQProjectionView.
void HiLite(int trkId, bool hlt=true)
size_type size() const
Definition: LArVoxelList.h:140
Container of LAr voxel information.
std::string string
Definition: nybbler.cc:12
double ActiveHalfHeight() const
Half height (associated with y coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:99
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
int Mother() const
Definition: MCParticle.h:213
simb::Origin_t Origin() const
Definition: MCTruth.h:74
double Mass() const
Definition: MCParticle.h:239
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double Px(const int i=0) const
Definition: MCParticle.h:230
static sim::LArVoxelList GetLArVoxelList(const art::Event &evt, std::string moduleLabel)
int GetMCTruth(const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
art::InputTag fG4ModuleLabel
module label producing sim::SimChannel objects
int StatusCode() const
Definition: MCParticle.h:211
geo::TPCID PositionToTPCID(geo::Point_t const &point) const
Returns the ID of the TPC at specified location.
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
intermediate_table::const_iterator const_iterator
static void FromPDG(TLine &line, int pdgcode)
Definition: Style.cxx:139
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
OrthoProj_t
Definition: OrthoProj.h:12
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
std::string Process() const
Definition: MCParticle.h:215
void MCTruth3D(const art::Event &evt, evdb::View3D *view)
Particle class.
static QStrList * l
Definition: config.cpp:1044
std::string encode() const
Definition: InputTag.cc:97
bool empty() const
Definition: MCTrajectory.h:167
art framework interface to geometry description
iterator begin()
Definition: LArVoxelList.h:131
int TrackId() const
Definition: MCParticle.h:210
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
list_type::const_iterator const_iterator
Definition: LArVoxelList.h:79
virtual bool EnableCorrSCE() const =0
auto & vals() noexcept
Definition: View.h:68
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
std::string const & label() const noexcept
Definition: InputTag.cc:79
bool isRealData() const
virtual geo::Vector_t GetPosOffsets(geo::Point_t const &point) const =0
const double e
double Y(const size_type i) const
Definition: MCTrajectory.h:150
Encapsulates the information we want store for a voxel.
LArSoft includes.
Definition: InfoTransfer.h:33
void swap(Handle< T > &a, Handle< T > &b)
unsigned int fCryostat
Cryostat number to draw, typically set by TWQProjectionView.
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
double Y() const
Definition: LArVoxelID.cxx:79
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
double ActiveHalfWidth() const
Half width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:95
int GetParticle(const art::Event &evt, std::vector< const simb::MCParticle * > &plist)
double P(const int i=0) const
Definition: MCParticle.h:234
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:65
double T(const int i=0) const
Definition: MCParticle.h:224
p
Definition: test.py:223
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
size_type NumberParticles() const
Definition: LArVoxelData.h:196
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: fwd.h:46
static int max(int a, int b)
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
Description of geometry of one entire detector.
double ActiveLength() const
Length (associated with z coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:103
double Z() const
Definition: LArVoxelID.cxx:86
void MCTruthShortText(const art::Event &evt, evdb::View2D *view)
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
const sim::ParticleList & ParticleList() const
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
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:184
double X() const
Definition: LArVoxelID.cxx:72
double Vx(const int i=0) const
Definition: MCParticle.h:221
mapped_type Energy() const
Definition: LArVoxelData.h:203
art::InputTag fSimChannelLabel
SimChannels may be independent of MC stuff.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
size_type size() const
Definition: MCTrajectory.h:166
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
static const char * LatexName(int pdgcode)
Convert PDG code to a latex string (root-style)
Definition: Style.cxx:12
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Pz(const int i=0) const
Definition: MCParticle.h:232
int fAxisOrientation
0 = TDC values on y-axis, wire number on x-axis, 1 = swapped
double Vz(const int i=0) const
Definition: MCParticle.h:223
int trigger_offset(DetectorClocksData const &data)
void MCTruthOrtho(const art::Event &evt, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
void MCTruthVectors2D(const art::Event &evt, evdb::View2D *view, unsigned int plane)
sim::LArVoxelID VoxelID() const
Definition: LArVoxelData.h:195
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:500
TCEvent evt
Definition: DataStructs.cxx:7
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
LArSoft geometry interface.
Definition: ChannelGeo.h:16
geo::CryostatID::CryostatID_t FindCryostatAtPosition(geo::Point_t const &worldLoc) const
Returns the index of the cryostat at specified location.
double Vy(const int i=0) const
Definition: MCParticle.h:222
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:107
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Cosmic rays.
Definition: MCTruth.h:24
QTextStream & endl(QTextStream &s)
Encapsulate the construction of a single detector plane.
const key_type & TrackID(const size_type) const
Point GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.h:779