LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
evd::SimulationDrawer Class Reference

#include "SimulationDrawer.h"

Public Member Functions

 SimulationDrawer ()
 
 ~SimulationDrawer ()
 
void MCTruthShortText (const art::Event &evt, evdb::View2D *view)
 
void MCTruthLongText (const art::Event &evt, evdb::View2D *view)
 
void MCTruthVectors2D (const art::Event &evt, evdb::View2D *view, unsigned int plane)
 
void MCTruth3D (const art::Event &evt, evdb::View3D *view)
 
void MCTruthOrtho (const art::Event &evt, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
 
void HiLite (int trkId, bool hlt=true)
 

Public Attributes

double minx
 
double maxx
 
double miny
 
double maxy
 
double minz
 
double maxz
 

Private Member Functions

int GetMCTruth (const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
 
int GetParticle (const art::Event &evt, std::vector< const simb::MCParticle * > &plist)
 

Private Attributes

std::map< int, bool > fHighlite
 

Detailed Description

Definition at line 29 of file SimulationDrawer.h.

Constructor & Destructor Documentation

evd::SimulationDrawer::SimulationDrawer ( )

Definition at line 55 of file SimulationDrawer.cxx.

References geo::GeometryCore::Iterate(), maxx, maxy, maxz, minx, miny, and minz.

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 (auto const& tpc : geom->Iterate<geo::TPCGeo>()) {
69  auto const tpc_center = tpc.Center();
70  auto const active_volume_center = tpc.GetActiveVolumeCenter();
71  mf::LogDebug("SimulationDrawer") << tpc.ID() << ", TPC center: " << tpc_center.X() << ", "
72  << tpc_center.Y() << ", " << tpc_center.Z() << std::endl;
73  mf::LogDebug("SimulationDrawer")
74  << " TPC Active center: " << active_volume_center.X() << ", "
75  << active_volume_center.Y() << ", " << active_volume_center.Z()
76  << ", H/W/L: " << tpc.ActiveHalfHeight() << "/" << tpc.ActiveHalfWidth() << "/"
77  << tpc.ActiveLength() << std::endl;
78 
79  if (minx > tpc_center.X() - tpc.HalfWidth()) minx = tpc_center.X() - tpc.HalfWidth();
80  if (maxx < tpc_center.X() + tpc.HalfWidth()) maxx = tpc_center.X() + tpc.HalfWidth();
81  if (miny > tpc_center.Y() - tpc.HalfHeight()) miny = tpc_center.Y() - tpc.HalfHeight();
82  if (maxy < tpc_center.Y() + tpc.HalfHeight()) maxy = tpc_center.Y() + tpc.HalfHeight();
83  if (minz > tpc_center.Z() - tpc.Length() / 2.) minz = tpc_center.Z() - tpc.Length() / 2.;
84  if (maxz < tpc_center.Z() + tpc.Length() / 2.) maxz = tpc_center.Z() + tpc.Length() / 2.;
85 
86  mf::LogDebug("SimulationDrawer")
87  << " minx/maxx: " << minx << "/" << maxx << ", miny/maxy: " << miny << "/" << maxy
88  << ", minz/miny: " << minz << "/" << maxz << std::endl;
89  }
90  }
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
Geometry information for a single TPC.
Definition: TPCGeo.h:36
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
evd::SimulationDrawer::~SimulationDrawer ( )

Definition at line 94 of file SimulationDrawer.cxx.

94 {}

Member Function Documentation

int evd::SimulationDrawer::GetMCTruth ( const art::Event evt,
std::vector< const simb::MCTruth * > &  mctruth 
)
private

Definition at line 913 of file SimulationDrawer.cxx.

References e, art::ProductRetriever::getMany(), and art::Event::isRealData().

Referenced by MCTruth3D(), MCTruthLongText(), MCTruthShortText(), and MCTruthVectors2D().

914  {
915  mcvec.clear();
916 
917  if (evt.isRealData()) return 0;
918 
919  std::vector<const simb::MCTruth*> temp;
920 
921  std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
922  // use get by Type because there should only be one collection of these in the event
923  try {
924  //evt.getManyByType(mctcol);
925  mctcol = evt.getMany<std::vector<simb::MCTruth>>();
926  for (size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
927  art::Handle<std::vector<simb::MCTruth>> mclistHandle = mctcol[mctc];
928 
929  for (size_t i = 0; i < mclistHandle->size(); ++i) {
930  temp.push_back(&(mclistHandle->at(i)));
931  }
932  }
933  temp.swap(mcvec);
934  }
935  catch (cet::exception& e) {
936  writeErrMsg("GetMCTruth", e);
937  }
938 
939  return mcvec.size();
940  }
bool isRealData() const
Definition: Event.cc:53
Float_t e
Definition: plot.C:35
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
int evd::SimulationDrawer::GetParticle ( const art::Event evt,
std::vector< const simb::MCParticle * > &  plist 
)
private

Definition at line 884 of file SimulationDrawer.cxx.

References e, evd::SimulationDrawingOptions::fG4ModuleLabel, art::ProductRetriever::getView(), art::Event::isRealData(), and art::View< T >::vals().

Referenced by MCTruth3D(), and MCTruthOrtho().

886  {
887  plist.clear();
888 
889  if (evt.isRealData()) return 0;
890 
892 
893  std::vector<const simb::MCParticle*> temp;
894 
896  // use get by Type because there should only be one collection of these in the event
897  try {
898  evt.getView(drawopt->fG4ModuleLabel, plcol);
899  for (unsigned int i = 0; i < plcol.vals().size(); ++i) {
900  temp.push_back(plcol.vals().at(i));
901  }
902  temp.swap(plist);
903  }
904  catch (cet::exception& e) {
905  writeErrMsg("GetRawDigits", e);
906  }
907 
908  return plist.size();
909  }
art::InputTag fG4ModuleLabel
module label producing sim::SimChannel objects
bool isRealData() const
Definition: Event.cc:53
auto & vals() noexcept
Definition: View.h:68
Definition: fwd.h:46
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Float_t e
Definition: plot.C:35
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evd::SimulationDrawer::HiLite ( int  trkId,
bool  hlt = true 
)

Definition at line 944 of file SimulationDrawer.cxx.

References fHighlite.

945  {
946  fHighlite[trkId] = dohilite;
947  }
std::map< int, bool > fHighlite
void evd::SimulationDrawer::MCTruth3D ( const art::Event evt,
evdb::View3D view 
)

Definition at line 327 of file SimulationDrawer.cxx.

References evdb::View3D::AddPolyLine3D(), evdb::View3D::AddPolyMarker3D(), sim::LArVoxelList::begin(), evd::Style::ColorFromPDG(), simb::MCParticle::E(), simb::MCTrajectory::empty(), sim::LArVoxelList::end(), sim::LArVoxelData::Energy(), evd::SimulationDrawingOptions::fG4ModuleLabel, evd::SimulationDrawingOptions::fMinEnergyDeposition, evd::SimulationDrawingOptions::fShowMCTruthColors, evd::SimulationDrawingOptions::fShowMCTruthFullSize, evd::SimulationDrawingOptions::fShowMCTruthTrajectories, sim::SimListUtils::GetLArVoxelList(), GetMCTruth(), GetParticle(), art::Event::isRealData(), art::InputTag::label(), simb::MCParticle::Mother(), sim::LArVoxelData::NumberParticles(), simb::MCParticle::PdgCode(), geo::GeometryCore::PositionToTPCID(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), simb::MCTrajectory::size(), sim::LArVoxelList::size(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::TrackId(), sim::LArVoxelData::TrackID(), simb::MCParticle::Trajectory(), detinfo::trigger_offset(), sim::LArVoxelData::VoxelID(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), simb::MCTrajectory::X(), sim::LArVoxelID::X(), simb::MCTrajectory::Y(), sim::LArVoxelID::Y(), simb::MCTrajectory::Z(), and sim::LArVoxelID::Z().

328  {
329  if (evt.isRealData()) return;
330 
332  // If the option is turned off, there's nothing to do
333  if (!drawopt->fShowMCTruthTrajectories) return;
334 
335  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
336  auto const detProp =
339 
340  // get the particles from the Geant4 step
341  std::vector<const simb::MCParticle*> plist;
342  this->GetParticle(evt, plist);
343 
344  // Define a couple of colors for neutrals and if we gray it out...
345  int neutralColor(12);
346  int grayedColor(15);
347  int neutrinoColor(38);
348 
349  // Use the LArVoxelList to get the true energy deposition locations as opposed to using MCTrajectories
350  const sim::LArVoxelList voxels =
352 
353  mf::LogDebug("SimulationDrawer")
354  << "Starting loop over " << plist.size() << " McParticles, voxel list size is "
355  << voxels.size() << std::endl;
356 
357  // Using the voxel information can be slow (see previous implementation of this code).
358  // In order to speed things up we have modified the strategy:
359  // 1) Make one pass through the list of voxels
360  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
361  // which is done by keeping a map between the MCParticle and a vector of positions
362  // 3) Then loop through the map to draw the particle trajectories.
363  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
364  // more loop to make a map of track id's and MCParticles.
365 
366  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
367  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
368 
369  // Should we display the trajectories too?
370  double minPartEnergy(0.01);
371 
372  for (size_t p = 0; p < plist.size(); ++p) {
373  trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
374 
375  // Quick loop through to draw trajectories...
376  if (drawopt->fShowMCTruthTrajectories) {
377  // Is there an associated McTrajectory?
378  const simb::MCParticle* mcPart = plist[p];
379  const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
380 
381  int pdgCode(mcPart->PdgCode());
382  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
383  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
384  double partCharge = partPDG ? partPDG->Charge() : 0.;
385  double partEnergy = mcPart->E();
386 
387  if (!drawopt->fShowMCTruthColors) colorIdx = grayedColor;
388 
389  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000) {
390  double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
391  double xOffset = 0.;
392  double xPosMinTick = 0.;
393  double xPosMaxTick = std::numeric_limits<double>::max();
394 
395  // collect the points from this particle
396  int numTrajPoints = mcTraj.size();
397 
398  std::unique_ptr<double[]> hitPositions(new double[3 * numTrajPoints]);
399  int hitCount(0);
400 
401  for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
402  double xPos = mcTraj.X(hitIdx);
403  double yPos = mcTraj.Y(hitIdx);
404  double zPos = mcTraj.Z(hitIdx);
405 
406  // If we have cosmic rays then we need to get the offset which allows translating from
407  // when they were generated vs when they were tracked.
408  // Note that this also explicitly checks that they are in a TPC volume
409  geo::Point_t hitLocation(xPos, yPos, zPos);
410 
411  try {
412  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
413  geo::PlaneID planeID(tpcID, 0);
414 
415  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
416  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
417  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
418 
419  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
420  }
421  catch (...) {
422  continue;
423  }
424 
425  // Now move the hit position to correspond to the timing
426  xPos += xOffset;
427 
428  // Check fiducial limits
429  if (xPos > xPosMinTick && xPos < xPosMaxTick) {
430  // Check for space charge offsets
431  // if (spaceCharge->EnableSimEfieldSCE())
432  // {
433  // std::vector<double> offsetVec = spaceCharge->GetPosOffsets(xPos,yPos,zPos);
434  // xPos += offsetVec[0] - 0.7;
435  // yPos -= offsetVec[1];
436  // zPos -= offsetVec[2];
437  // }
438 
439  hitPositions[3 * hitCount] = xPos;
440  hitPositions[3 * hitCount + 1] = yPos;
441  hitPositions[3 * hitCount + 2] = zPos;
442  hitCount++;
443  }
444  }
445 
446  TPolyLine3D& pl(view->AddPolyLine3D(1, colorIdx, 1, 1));
447 
448  // Draw neutrals as a gray dotted line to help fade into background a bit...
449  if (partCharge == 0.) {
450  pl.SetLineColor(neutralColor);
451  pl.SetLineStyle(3);
452  pl.SetLineWidth(1);
453  }
454  pl.SetPolyLine(hitCount, hitPositions.get(), "");
455  }
456  }
457  }
458 
459  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
460  std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
461 
463  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
464  const sim::LArVoxelData& vxd = (*vxitr).second;
465 
466  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
467  if (vxd.Energy(partIdx) > drawopt->fMinEnergyDeposition) {
468  int trackId = vxd.TrackID(partIdx);
469 
470  // It can be in some instances that mcPart here could be zero.
471  const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
472 
473  partToPosMap[mcPart].push_back(std::vector<double>(3));
474 
475  partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
476  partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
477  partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
478  }
479  } // end if this track id is in the current voxel
480  } // end loop over voxels
481 
482  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
483  // draw the trajectories
484  std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
485 
486  for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
487  partToPosMapItr++) {
488  // Recover the McParticle, we'll need to access several data members so may as well dereference it
489  const simb::MCParticle* mcPart = partToPosMapItr->first;
490 
491  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
492  if (!mcPart || partToPosMapItr->second.empty()) continue;
493 
494  double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
495  double xOffset = 0.;
496  double xPosMinTick = 0.;
497  double xPosMaxTick = std::numeric_limits<double>::max();
498 
499  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
500  int markerIdx(kFullDotSmall);
501  int markerSize(2);
502 
503  if (!drawopt->fShowMCTruthFullSize) {
504  colorIdx = grayedColor;
505  markerIdx = kDot;
506  markerSize = 1;
507  }
508 
509  std::unique_ptr<double[]> hitPositions(new double[3 * partToPosMapItr->second.size()]);
510  int hitCount(0);
511 
512  // Now loop over points and add to trajectory
513  for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
514  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
515 
516  // Check xOffset state and set if necessary
517  geo::Point_t hitLocation(posVec[0], posVec[1], posVec[2]);
518 
519  try {
520  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
521  geo::PlaneID planeID(tpcID, 0);
522 
523  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
524  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
525  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
526 
527  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
528  }
529  catch (...) {
530  continue;
531  }
532 
533  double xCoord = posVec[0] + xOffset;
534 
535  // If a voxel records an energy deposit then must have been in the TPC
536  // But because things get shifted still need to cut off if outside drift
537  if (xCoord > xPosMinTick && xCoord < xPosMaxTick) {
538  hitPositions[3 * hitCount] = xCoord;
539  hitPositions[3 * hitCount + 1] = posVec[1];
540  hitPositions[3 * hitCount + 2] = posVec[2];
541  hitCount++;
542  }
543  }
544 
545  TPolyMarker3D& pm = view->AddPolyMarker3D(1, colorIdx, markerIdx, markerSize);
546  pm.SetPolyMarker(hitCount, hitPositions.get(), markerIdx);
547  }
548 
549  // Finally, let's see if we can draw the incoming particle from the MCTruth information
550  std::vector<const simb::MCTruth*> mctruth;
551  this->GetMCTruth(evt, mctruth);
552 
553  // Loop through the MCTruth vector
554  for (unsigned int idx = 0; idx < mctruth.size(); idx++) {
555  // Go through each MCTruth object in the list
556  for (int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++) {
557  const simb::MCParticle& mcPart = mctruth[idx]->GetParticle(particleIdx);
558 
559  // A negative mother id indicates the "primary" particle
560  if (mcPart.Mother() == -1 && mcPart.StatusCode() == 0) {
561  mf::LogDebug("SimulationDrawer") << mcPart << std::endl;
562 
563  // Get position vector
564  TVector3 particlePosition(mcPart.Vx(), mcPart.Vy(), mcPart.Vz());
565 
566  // Get direction vector (in opposite direction)
567  TVector3 oppPartDir(-mcPart.Px(), -mcPart.Py(), -mcPart.Pz());
568 
569  if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
570 
571  double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
572 
573  // No point in drawing if arc length is zero (e.g. Ar nucleus)
574  if (arcLenToDraw > 0.) {
575  // Draw the line, use an off color to be unique
576  TPolyLine3D& pl(view->AddPolyLine3D(2, neutrinoColor, 1, 2));
577 
578  pl.SetPoint(0, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
579 
580  particlePosition += std::min(arcLenToDraw + 10., 1000.) * oppPartDir;
581 
582  pl.SetPoint(1, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
583  }
584  }
585  // The particles we want to draw will be early in the list so break out if we didn't find them
586  else
587  break;
588  } // loop on particles in list
589  }
590 
591  return;
592  }
double E(const int i=0) const
Definition: MCParticle.h:234
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:213
double Py(const int i=0) const
Definition: MCParticle.h:232
size_type size() const
Definition: LArVoxelList.h:139
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:254
int Mother() const
Definition: MCParticle.h:214
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
double Px(const int i=0) const
Definition: MCParticle.h:231
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:212
bool isRealData() const
Definition: Event.cc:53
bool empty() const
Definition: MCTrajectory.h:167
iterator begin()
Definition: LArVoxelList.h:130
int TrackId() const
Definition: MCParticle.h:211
list_type::const_iterator const_iterator
Definition: LArVoxelList.h:78
TPolyLine3D & AddPolyLine3D(int n, int c, int w, int s)
Definition: View3D.cxx:105
std::string const & label() const noexcept
Definition: InputTag.cc:79
double Y(const size_type i) const
Definition: MCTrajectory.h:150
double Y() const
Definition: LArVoxelID.cxx:75
int GetParticle(const art::Event &evt, std::vector< const simb::MCParticle * > &plist)
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:65
double T(const int i=0) const
Definition: MCParticle.h:225
size_type NumberParticles() const
Definition: LArVoxelData.h:200
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
double Z() const
Definition: LArVoxelID.cxx:82
double X() const
Definition: LArVoxelID.cxx:68
double Vx(const int i=0) const
Definition: MCParticle.h:222
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
mapped_type Energy() const
Definition: LArVoxelData.h:215
size_type size() const
Definition: MCTrajectory.h:166
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
double Pz(const int i=0) const
Definition: MCParticle.h:233
double Vz(const int i=0) const
Definition: MCParticle.h:224
int trigger_offset(DetectorClocksData const &data)
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
sim::LArVoxelID VoxelID() const
Definition: LArVoxelData.h:196
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
double Vy(const int i=0) const
Definition: MCParticle.h:223
TPolyMarker3D & AddPolyMarker3D(int n, int c, int st, double sz)
Definition: View3D.cxx:75
const key_type & TrackID(const size_type) const
void evd::SimulationDrawer::MCTruthLongText ( const art::Event evt,
evdb::View2D view 
)

Definition at line 159 of file SimulationDrawer.cxx.

References simb::MCParticle::E(), evd::SimulationDrawingOptions::fShowMCTruthText, GetMCTruth(), art::Event::isRealData(), evd::Style::LatexName(), simb::MCParticle::Mass(), simb::MCParticle::Mother(), simb::MCParticle::P(), simb::MCParticle::PdgCode(), simb::MCParticle::Process(), art::right(), simb::MCParticle::StatusCode(), and simb::MCParticle::TrackId().

Referenced by evd::MCBriefPad::Draw().

160  {
161  if (evt.isRealData()) return;
162 
164  // Skip drawing if option is turned off
165  if (!drawopt->fShowMCTruthText) return;
166 
167  std::vector<const simb::MCTruth*> mctruth;
168  this->GetMCTruth(evt, mctruth);
169  std::cout << "\nMCTruth Ptcl trackID PDG P T Moth Process\n";
170  for (unsigned int i = 0; i < mctruth.size(); ++i) {
171  for (int j = 0; j < mctruth[i]->NParticles(); ++j) {
172  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
173  if (p.StatusCode() == 0 || p.StatusCode() == 1) {
174  int KE = 1000 * (p.E() - p.Mass());
175  std::cout << std::right << std::setw(7) << i << std::setw(5) << j << std::setw(8)
176  << p.TrackId() << " " << std::setw(14) << Style::LatexName(p.PdgCode())
177  << std::setw(7) << int(1000 * p.P()) << std::setw(7) << KE << std::setw(7)
178  << p.Mother() << " " << p.Process() << "\n";
179  }
180  } // loop on j particles in list
181  }
182  std::cout << "Note: Momentum, P, and kinetic energy, T, in MeV/c\n";
183  } // MCTruthLongText
double E(const int i=0) const
Definition: MCParticle.h:234
int PdgCode() const
Definition: MCParticle.h:213
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
int Mother() const
Definition: MCParticle.h:214
double Mass() const
Definition: MCParticle.h:240
int GetMCTruth(const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
int StatusCode() const
Definition: MCParticle.h:212
bool isRealData() const
Definition: Event.cc:53
std::string Process() const
Definition: MCParticle.h:216
int TrackId() const
Definition: MCParticle.h:211
double P(const int i=0) const
Definition: MCParticle.h:235
static const char * LatexName(int pdgcode)
Convert PDG code to a latex string (root-style)
Definition: Style.cxx:12
void evd::SimulationDrawer::MCTruthOrtho ( const art::Event evt,
evd::OrthoProj_t  proj,
double  msize,
evdb::View2D view 
)

Definition at line 596 of file SimulationDrawer.cxx.

References evdb::View2D::AddPolyLine(), evdb::View2D::AddPolyMarker(), sim::LArVoxelList::begin(), evd::Style::ColorFromPDG(), geo::CryostatID::Cryostat, simb::MCParticle::E(), simb::MCTrajectory::empty(), art::InputTag::encode(), sim::LArVoxelList::end(), sim::LArVoxelData::Energy(), geo::GeometryCore::FindTPCAtPosition(), evd::SimulationDrawingOptions::fMinEnergyDeposition, evd::SimulationDrawingOptions::fShowMCTruthTrajectories, evd::SimulationDrawingOptions::fSimChannelLabel, geo::GeometryCore::GetElement(), sim::SimListUtils::GetLArVoxelList(), GetParticle(), art::Event::isRealData(), geo::CryostatID::isValid, evd::kXY, evd::kXZ, evd::kYZ, maxx, geo::BoxBoundedGeo::MaxX(), maxy, maxz, minx, geo::BoxBoundedGeo::MinX(), minz, sim::LArVoxelData::NumberParticles(), simb::MCParticle::PdgCode(), geo::GeometryCore::PositionToCryostatID(), simb::MCTrajectory::size(), sim::LArVoxelList::size(), simb::MCParticle::T(), geo::TPCID::TPC, simb::MCParticle::TrackId(), sim::LArVoxelData::TrackID(), simb::MCParticle::Trajectory(), detinfo::trigger_offset(), sim::LArVoxelData::VoxelID(), simb::MCTrajectory::X(), sim::LArVoxelID::X(), simb::MCTrajectory::Y(), sim::LArVoxelID::Y(), simb::MCTrajectory::Z(), and sim::LArVoxelID::Z().

Referenced by evd::Ortho3DPad::Draw().

600  {
601  if (evt.isRealData()) return;
602 
604 
605  // If the option is turned off, there's nothing to do
606  if (!drawopt->fShowMCTruthTrajectories) return;
607 
608  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
609  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
610  auto const detProp =
612 
613  // get the particles from the Geant4 step
614  std::vector<const simb::MCParticle*> plist;
615  this->GetParticle(evt, plist);
616 
617  // Useful variables
618 
619  double xMinimum(-1. * (maxx - minx));
620  double xMaximum(2. * (maxx - minx));
621 
622  // Use the LArVoxelList to get the true energy deposition locations as
623  // opposed to using MCTrajectories
624  // GetLarVoxelList should be adjusted to take an InputTag instead of strange.
625  const sim::LArVoxelList voxels =
627 
628  mf::LogDebug("SimulationDrawer")
629  << "Starting loop over " << plist.size() << " McParticles, voxel list size is "
630  << voxels.size() << std::endl;
631 
632  // Using the voxel information can be slow (see previous implementation of this code).
633  // In order to speed things up we have modified the strategy:
634  // 1) Make one pass through the list of voxels
635  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
636  // which is done by keeping a map between the MCParticle and a vector of positions
637  // 3) Then loop through the map to draw the particle trajectories.
638  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
639  // more loop to make a map of track id's and MCParticles.
640 
641  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
642  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
643 
644  // Should we display the trajectories too?
645  bool displayMcTrajectories(true);
646  double minPartEnergy(0.025);
647 
648  double tpcminx = 1.0;
649  double tpcmaxx = -1.0;
650  double xOffset = 0.0;
651  double g4Ticks = 0.0;
652  double coeff = 0.0;
653  double readoutwindowsize = 0.0;
654  for (size_t p = 0; p < plist.size(); ++p) {
655  trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
656 
657  // Quick loop through to drawn trajectories...
658  if (displayMcTrajectories) {
659  // Is there an associated McTrajectory?
660  const simb::MCParticle* mcPart = plist[p];
661  const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
662 
663  int pdgCode(mcPart->PdgCode());
664  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
665  double partCharge = partPDG ? partPDG->Charge() : 0.;
666  double partEnergy = mcPart->E();
667 
668  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000) {
669  // collect the points from this particle
670  int numTrajPoints = mcTraj.size();
671 
672  std::unique_ptr<double[]> hitPosX(new double[numTrajPoints]);
673  std::unique_ptr<double[]> hitPosY(new double[numTrajPoints]);
674  std::unique_ptr<double[]> hitPosZ(new double[numTrajPoints]);
675  int hitCount(0);
676 
677  double xPos = mcTraj.X(0);
678  double yPos = mcTraj.Y(0);
679  double zPos = mcTraj.Z(0);
680 
681  tpcminx = 1.0;
682  tpcmaxx = -1.0;
683  xOffset = 0.0;
684  g4Ticks = 0.0;
685  coeff = 0.0;
686  readoutwindowsize = 0.0;
687  for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
688  xPos = mcTraj.X(hitIdx);
689  yPos = mcTraj.Y(hitIdx);
690  zPos = mcTraj.Z(hitIdx);
691 
692  // If the original simulated hit did not occur in the TPC volume then don't draw it
693  if (xPos < minx || xPos > maxx || yPos < miny || yPos > maxy || zPos < minz ||
694  zPos > maxz)
695  continue;
696 
697  if ((xPos < tpcminx) || (xPos > tpcmaxx)) {
698  geo::Point_t const vtx{xPos, yPos, zPos};
699  geo::TPCID tpcid = geom->FindTPCAtPosition(vtx);
700  unsigned int cryo = geom->PositionToCryostatID(vtx).Cryostat;
701 
702  if (tpcid.isValid) {
703  unsigned int tpc = tpcid.TPC;
704  const geo::TPCGeo& tpcgeo = geom->GetElement(tpcid);
705  tpcminx = tpcgeo.MinX();
706  tpcmaxx = tpcgeo.MaxX();
707 
708  coeff = detProp.GetXTicksCoefficient(tpc, cryo);
709  readoutwindowsize =
710  detProp.ConvertTicksToX(detProp.ReadOutWindowSize(), 0, tpc, cryo);
711 
712  // The following is meant to get the correct offset for drawing
713  // the particle trajectory In particular, the cosmic rays will
714  // not be correctly placed without this
715  g4Ticks = clockData.TPCG4Time2Tick(mcPart->T()) +
716  detProp.GetXTicksOffset(0, tpc, cryo) - trigger_offset(clockData);
717 
718  xOffset = detProp.ConvertTicksToX(g4Ticks, 0, tpc, cryo);
719  }
720  else {
721  xOffset = 0;
722  tpcminx = 1.0;
723  tpcmaxx = -1.0;
724  coeff = 0.0;
725  readoutwindowsize = 0.0;
726  }
727  }
728 
729  // Now move the hit position to correspond to the timing
730  xPos += xOffset;
731 
732  bool inreadoutwindow = false;
733  if (coeff < 0) {
734  if ((xPos > readoutwindowsize) && (xPos < tpcmaxx)) inreadoutwindow = true;
735  }
736  else if (coeff > 0) {
737  if ((xPos > tpcminx) && (xPos < readoutwindowsize)) inreadoutwindow = true;
738  }
739 
740  if (!inreadoutwindow) continue;
741 
742  // Check fiducial limits
743  if (xPos > xMinimum && xPos < xMaximum) {
744  hitPosX[hitCount] = xPos;
745  hitPosY[hitCount] = yPos;
746  hitPosZ[hitCount] = zPos;
747  hitCount++;
748  }
749  }
750 
751  TPolyLine& pl = view->AddPolyLine(
752  1, evd::Style::ColorFromPDG(mcPart->PdgCode()), 1, 1); //kFullCircle, msize);
753 
754  // Draw neutrals as a gray dotted line to help fade into background a bit...
755  if (partCharge == 0.) {
756  pl.SetLineColor(13);
757  pl.SetLineStyle(3);
758  pl.SetLineWidth(1);
759  }
760 
761  if (proj == evd::kXY)
762  pl.SetPolyLine(hitCount, hitPosX.get(), hitPosY.get(), "");
763  else if (proj == evd::kXZ)
764  pl.SetPolyLine(hitCount, hitPosZ.get(), hitPosX.get(), "");
765  else if (proj == evd::kYZ)
766  pl.SetPolyLine(hitCount, hitPosZ.get(), hitPosY.get(), "");
767  }
768  }
769  }
770 
771  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
772  std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
773 
775  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
776  const sim::LArVoxelData& vxd = (*vxitr).second;
777 
778  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
779  if (vxd.Energy(partIdx) > drawopt->fMinEnergyDeposition) {
780  int trackId = vxd.TrackID(partIdx);
781 
782  // It can be in some instances that mcPart here could be zero.
783  const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
784 
785  partToPosMap[mcPart].push_back(std::vector<double>(3));
786 
787  partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
788  partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
789  partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
790  }
791  } // end if this track id is in the current voxel
792  } // end loop over voxels
793 
794  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
795  // draw the trajectories
796  std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
797 
798  for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
799  partToPosMapItr++) {
800  // Recover the McParticle, we'll need to access several data members so may as well dereference it
801  const simb::MCParticle* mcPart = partToPosMapItr->first;
802 
803  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
804  if (!mcPart || partToPosMapItr->second.empty()) continue;
805 
806  tpcminx = 1.0;
807  tpcmaxx = -1.0;
808  xOffset = 0.0;
809  g4Ticks = 0.0;
810  std::vector<std::array<double, 3>> posVecCorr;
811  posVecCorr.reserve(partToPosMapItr->second.size());
812  coeff = 0.0;
813  readoutwindowsize = 0.0;
814 
815  // Now loop over points and add to trajectory
816  for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
817  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
818 
819  if ((posVec[0] < tpcminx) || (posVec[0] > tpcmaxx)) {
820  geo::Point_t const vtx{posVec[0], posVec[1], posVec[2]};
821  geo::TPCID tpcid = geom->FindTPCAtPosition(vtx);
822  unsigned int cryo = geom->PositionToCryostatID(vtx).Cryostat;
823 
824  if (tpcid.isValid) {
825  unsigned int tpc = tpcid.TPC;
826 
827  const geo::TPCGeo& tpcgeo = geom->GetElement(tpcid);
828  tpcminx = tpcgeo.MinX();
829  tpcmaxx = tpcgeo.MaxX();
830 
831  coeff = detProp.GetXTicksCoefficient(tpc, cryo);
832  readoutwindowsize = detProp.ConvertTicksToX(detProp.ReadOutWindowSize(), 0, tpc, cryo);
833  // The following is meant to get the correct offset for drawing the
834  // particle trajectory In particular, the cosmic rays will not be
835  // correctly placed without this
836  g4Ticks = clockData.TPCG4Time2Tick(mcPart->T()) +
837  detProp.GetXTicksOffset(0, tpc, cryo) - trigger_offset(clockData);
838 
839  xOffset = detProp.ConvertTicksToX(g4Ticks, 0, tpc, cryo);
840  }
841  else {
842  xOffset = 0;
843  tpcminx = 1.0;
844  tpcmaxx = -1.0;
845  coeff = 0.0;
846  readoutwindowsize = 0.0;
847  }
848  }
849 
850  double xCoord = posVec[0] + xOffset;
851 
852  bool inreadoutwindow = false;
853  if (coeff < 0) {
854  if ((xCoord > readoutwindowsize) && (xCoord < tpcmaxx)) inreadoutwindow = true;
855  }
856  else if (coeff > 0) {
857  if ((xCoord > tpcminx) && (xCoord < readoutwindowsize)) inreadoutwindow = true;
858  }
859 
860  if (inreadoutwindow && (xCoord > xMinimum && xCoord < xMaximum)) {
861  posVecCorr.push_back({{xCoord, posVec[1], posVec[2]}});
862  }
863  }
864 
865  TPolyMarker& pm = view->AddPolyMarker(posVecCorr.size(),
867  kFullDotMedium,
868  2); //kFullCircle, msize);
869 
870  for (size_t p = 0; p < posVecCorr.size(); ++p) {
871  if (proj == evd::kXY)
872  pm.SetPoint(p, posVecCorr[p][0], posVecCorr[p][1]);
873  else if (proj == evd::kXZ)
874  pm.SetPoint(p, posVecCorr[p][2], posVecCorr[p][0]);
875  else if (proj == evd::kYZ)
876  pm.SetPoint(p, posVecCorr[p][2], posVecCorr[p][1]);
877  }
878  }
879 
880  return;
881  }
double E(const int i=0) const
Definition: MCParticle.h:234
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:213
TPolyLine & AddPolyLine(int n, int c, int w, int s)
Definition: View2D.cxx:210
size_type size() const
Definition: LArVoxelList.h:139
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:254
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:90
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
Geometry information for a single TPC.
Definition: TPCGeo.h:36
static sim::LArVoxelList GetLArVoxelList(const art::Event &evt, std::string moduleLabel)
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:93
bool isRealData() const
Definition: Event.cc:53
CryostatGeo const & GetElement(CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.h:465
std::string encode() const
Definition: InputTag.cc:97
bool empty() const
Definition: MCTrajectory.h:167
iterator begin()
Definition: LArVoxelList.h:130
int TrackId() const
Definition: MCParticle.h:211
list_type::const_iterator const_iterator
Definition: LArVoxelList.h:78
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
double Y(const size_type i) const
Definition: MCTrajectory.h:150
double Y() const
Definition: LArVoxelID.cxx:75
int GetParticle(const art::Event &evt, std::vector< const simb::MCParticle * > &plist)
TPolyMarker & AddPolyMarker(int n, int c, int st, double sz)
Definition: View2D.cxx:157
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:65
double T(const int i=0) const
Definition: MCParticle.h:225
size_type NumberParticles() const
Definition: LArVoxelData.h:200
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
double Z() const
Definition: LArVoxelID.cxx:82
double X() const
Definition: LArVoxelID.cxx:68
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
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
mapped_type Energy() const
Definition: LArVoxelData.h:215
art::InputTag fSimChannelLabel
SimChannels may be independent of MC stuff.
size_type size() const
Definition: MCTrajectory.h:166
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Float_t proj
Definition: plot.C:35
int trigger_offset(DetectorClocksData const &data)
sim::LArVoxelID VoxelID() const
Definition: LArVoxelData.h:196
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
const key_type & TrackID(const size_type) const
void evd::SimulationDrawer::MCTruthShortText ( const art::Event evt,
evdb::View2D view 
)

Definition at line 98 of file SimulationDrawer.cxx.

References evdb::View2D::AddLatex(), evd::Style::ColorFromPDG(), evd::SimulationDrawingOptions::fShowMCTruthText, GetMCTruth(), art::Event::isRealData(), simb::kCosmicRay, evd::Style::LatexName(), geo::origin(), simb::MCParticle::P(), simb::MCParticle::PdgCode(), and simb::MCParticle::StatusCode().

Referenced by evd::MCBriefPad::Draw().

99  {
100 
101  if (evt.isRealData()) return;
102 
104  // Skip drawing if option is turned off
105  if (!drawopt->fShowMCTruthText) return;
106 
107  std::vector<const simb::MCTruth*> mctruth;
108  this->GetMCTruth(evt, mctruth);
109 
110  for (unsigned int i = 0; i < mctruth.size(); ++i) {
111  std::string mctext;
112  bool firstin = true;
113  bool firstout = true;
114  std::string origin;
115  std::string incoming;
116  std::string outgoing;
117  // Label cosmic rays -- others are pretty obvious
118  if (mctruth[i]->Origin() == simb::kCosmicRay) origin = "c-ray: ";
119  int jmax = TMath::Min(20, mctruth[i]->NParticles());
120  for (int j = 0; j < jmax; ++j) {
121  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
122  char buff[1024];
123  if (p.P() > 0.05) {
124  sprintf(buff,
125  "#color[%d]{%s #scale[0.75]{[%.1f GeV/c]}}",
128  p.P());
129  }
130  else {
131  sprintf(buff,
132  "#color[%d]{%s}",
135  }
136  if (p.StatusCode() == 0) {
137  if (firstin == false) incoming += " + ";
138  incoming += buff;
139  firstin = false;
140  }
141  if (p.StatusCode() == 1) {
142  if (firstout == false) outgoing += " + ";
143  outgoing += buff;
144  firstout = false;
145  }
146  } // loop on j particles
147  if (origin == "" && incoming == "") { mctext = outgoing; }
148  else {
149  mctext = origin + incoming + " #rightarrow " + outgoing;
150  }
151  TLatex& latex = view->AddLatex(0.03, 0.2, mctext.c_str());
152  latex.SetTextSize(0.6);
153 
154  } // loop on i mctruth objects
155  }
int PdgCode() const
Definition: MCParticle.h:213
int GetMCTruth(const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
int StatusCode() const
Definition: MCParticle.h:212
bool isRealData() const
Definition: Event.cc:53
double P(const int i=0) const
Definition: MCParticle.h:235
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:65
TLatex & AddLatex(double x, double y, const char *text)
Definition: View2D.cxx:308
static const char * LatexName(int pdgcode)
Convert PDG code to a latex string (root-style)
Definition: Style.cxx:12
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:229
Cosmic rays.
Definition: MCTruth.h:24
void evd::SimulationDrawer::MCTruthVectors2D ( const art::Event evt,
evdb::View2D view,
unsigned int  plane 
)

Definition at line 187 of file SimulationDrawer.cxx.

References evdb::View2D::AddLine(), sim::ParticleList::begin(), simb::MCParticle::E(), sim::ParticleList::empty(), spacecharge::SpaceCharge::EnableCorrSCE(), sim::ParticleList::end(), evd::RawDrawingOptions::fAxisOrientation, evd::RawDrawingOptions::fCryostat, evd::Style::FromPDG(), evd::RawDrawingOptions::fTPC, GetMCTruth(), spacecharge::SpaceCharge::GetPosOffsets(), ipart, art::Event::isRealData(), simb::kCosmicRay, simb::MCParticle::Mass(), simb::MCTruth::Origin(), simb::MCParticle::P(), cheat::ParticleInventoryService::ParticleList(), simb::MCParticle::PdgCode(), simb::MCParticle::Process(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), r, simb::MCParticle::StatusCode(), t1, t2, simb::MCParticle::TrackId(), cheat::ParticleInventoryService::TrackIdToMCTruth_P(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), and geo::GeometryCore::WireCoordinate().

Referenced by evd::TWireProjPad::Draw().

190  {
191  if (evt.isRealData()) return;
192 
193  const spacecharge::SpaceCharge* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
194 
196  // If the option is turned off, there's nothing to do
197  if (drawopt->fShowMCTruthVectors > 3) {
198  std::cout << "Unsupported ShowMCTruthVectors option (> 2)\n";
199  return;
200  }
201 
204  geo::PlaneID const planeID{rawopt->fCryostat, rawopt->fTPC, plane};
205 
206  // shift the truth by a fixed amount so it doesn't overlay the reco
207  double const xShift = -2;
208 
209  static bool first = true;
210  if (first) {
211  std::cout
212  << "******** Show MCTruth (Genie) particles when ShowMCTruthVectors = 1 or 3 ******** \n";
213  std::cout << " MCTruth vectors corrected for space charge? " << sce->EnableCorrSCE()
214  << " and shifted by " << xShift << " cm in X\n";
215  std::cout << " Neutrons and photons drawn with dotted lines. \n";
216  std::cout << " Red = e+/-, nue, nuebar. Blue = mu+/-, numu, numubar. Green = tau+/-, nutau, "
217  "nutaubar.\n";
218  std::cout << " Yellow = photons. Magenta = pions, protons and nuetrons.\n";
219  std::cout << "******** Show MCParticle (Geant) decay photons (e.g. from pizeros) when "
220  "ShowMCTruthVectors = 2 or 3 ******** \n";
221  std::cout << " Photons > 50 MeV are drawn as dotted lines corrected for space charge and "
222  "are not shifted.\n";
223  std::cout << " Decay photon end points are drawn at 2 interaction lengths (44 cm) from the "
224  "start position.\n";
225  std::cout << " Color: Green = (50 < E_photon < 100 MeV), Blue = (100 MeV < E_photon < 200 "
226  "MeV), Red = (E_photon > 300 MeV).\n";
227  }
228 
229  bool showTruth = (drawopt->fShowMCTruthVectors == 1 || drawopt->fShowMCTruthVectors == 3);
230  bool showPhotons = (drawopt->fShowMCTruthVectors > 1);
231 
232  auto const detProp =
234 
235  if (showTruth) {
236  // Unpack and draw the MC vectors
237  std::vector<const simb::MCTruth*> mctruth;
238  this->GetMCTruth(evt, mctruth);
239 
240  for (size_t i = 0; i < mctruth.size(); ++i) {
241  if (mctruth[i]->Origin() == simb::kCosmicRay) continue;
242  for (int j = 0; j < mctruth[i]->NParticles(); ++j) {
243  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
244 
245  // Skip all but incoming and out-going particles
246  if (!(p.StatusCode() == 0 || p.StatusCode() == 1)) continue;
247 
248  double r = p.P() * 10.0; // Scale length so 10 cm = 1 GeV/c
249 
250  if (p.StatusCode() == 0) r = -r; // Flip for incoming particles
251 
252  auto gptStart = geo::Point_t(p.Vx(), p.Vy(), p.Vz());
253  geo::Point_t sceOffset{0, 0, 0};
254  if (sce->EnableCorrSCE()) sceOffset = sce->GetPosOffsets(gptStart);
255  geo::Point_t const xyz1{
256  p.Vx() - sceOffset.X(), p.Vy() + sceOffset.Y(), p.Vz() + sceOffset.Z()};
257  geo::Point_t const xyz2{xyz1.X() + r * p.Px() / p.P(),
258  xyz1.Y() + r * p.Py() / p.P(),
259  xyz1.Z() + r * p.Pz() / p.P()};
260  double w1 = geo->WireCoordinate(xyz1, planeID);
261  double w2 = geo->WireCoordinate(xyz2, planeID);
262 
263  double time = detProp.ConvertXToTicks(xyz1.X() + xShift, planeID);
264  double time2 = detProp.ConvertXToTicks(xyz2.X() + xShift, planeID);
265 
266  if (rawopt->fAxisOrientation < 1) {
267  TLine& l = view->AddLine(w1, time, w2, time2);
269  }
270  else {
271  TLine& l = view->AddLine(time, w1, time2, w2);
273  }
274  //
275 
276  } // loop on j particles in list
277  } // loop on truths
278  } // showTruth
279 
280  if (showPhotons) {
281  // draw pizero photons with T > 30 MeV
283  sim::ParticleList const& plist = pi_serv->ParticleList();
284  if (plist.empty()) return;
285  // photon interaction length approximate
286  double r = 44;
287  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
288  simb::MCParticle* p = (*ipart).second;
289  int trackID = p->TrackId();
290  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
291  if (theTruth->Origin() == simb::kCosmicRay) continue;
292  if (p->PdgCode() != 22) continue;
293  if (p->Process() != "Decay") continue;
294  int TMeV = 1000 * (p->E() - p->Mass());
295  if (TMeV < 30) continue;
296  auto gptStart = geo::Point_t(p->Vx(), p->Vy(), p->Vz());
297  geo::Point_t sceOffset{0, 0, 0};
298  if (sce->EnableCorrSCE()) sceOffset = sce->GetPosOffsets(gptStart);
299  geo::Point_t const xyz1{
300  p->Vx() - sceOffset.X(), p->Vy() + sceOffset.Y(), p->Vz() + sceOffset.Z()};
301  geo::Point_t const xyz2{xyz1.X() + r * p->Px() / p->P(),
302  xyz1.Y() + r * p->Py() / p->P(),
303  xyz1.Z() + r * p->Pz() / p->P()};
304  double w1 = geo->WireCoordinate(xyz1, planeID);
305  double t1 = detProp.ConvertXToTicks(xyz1.X(), planeID);
306  double w2 = geo->WireCoordinate(xyz2, planeID);
307  double t2 = detProp.ConvertXToTicks(xyz2.X(), planeID);
308  TLine& l = view->AddLine(w1, t1, w2, t2);
309  l.SetLineWidth(2);
310  l.SetLineStyle(kDotted);
311  if (TMeV < 100) { l.SetLineColor(kGreen); }
312  else if (TMeV < 200) {
313  l.SetLineColor(kBlue);
314  }
315  else {
316  l.SetLineColor(kRed);
317  }
318  } // ipart
319  } // showPhotons
320 
321  first = false;
322 
323  } // MCTruthVectors2D
double E(const int i=0) const
Definition: MCParticle.h:234
TRandom r
Definition: spectrum.C:23
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
int PdgCode() const
Definition: MCParticle.h:213
bool empty() const
Definition: ParticleList.h:314
double Py(const int i=0) const
Definition: MCParticle.h:232
unsigned int fTPC
TPC number to draw, typically set by TWQProjectionView.
TTree * t1
Definition: plottest35.C:26
Int_t ipart
Definition: Style.C:10
TLine & AddLine(double x1, double y1, double x2, double y2)
Definition: View2D.cxx:187
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
simb::Origin_t Origin() const
Definition: MCTruth.h:74
double Mass() const
Definition: MCParticle.h:240
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
double Px(const int i=0) const
Definition: MCParticle.h:231
int GetMCTruth(const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
int StatusCode() const
Definition: MCParticle.h:212
static void FromPDG(TLine &line, int pdgcode)
Definition: Style.cxx:131
bool isRealData() const
Definition: Event.cc:53
std::string Process() const
Definition: MCParticle.h:216
int TrackId() const
Definition: MCParticle.h:211
virtual bool EnableCorrSCE() const =0
virtual geo::Vector_t GetPosOffsets(geo::Point_t const &point) const =0
unsigned int fCryostat
Cryostat number to draw, typically set by TWQProjectionView.
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
double P(const int i=0) const
Definition: MCParticle.h:235
iterator begin()
Definition: ParticleList.h:305
TTree * t2
Definition: plottest35.C:36
const sim::ParticleList & ParticleList() const
double Vx(const int i=0) const
Definition: MCParticle.h:222
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
double Pz(const int i=0) const
Definition: MCParticle.h:233
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:224
Namespace collecting geometry-related classes utilities.
double Vy(const int i=0) const
Definition: MCParticle.h:223
Cosmic rays.
Definition: MCTruth.h:24

Member Data Documentation

std::map<int, bool> evd::SimulationDrawer::fHighlite
private

Definition at line 59 of file SimulationDrawer.h.

Referenced by HiLite().

double evd::SimulationDrawer::maxx

Definition at line 48 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::maxy

Definition at line 50 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::maxz

Definition at line 52 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::minx

Definition at line 47 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::miny

Definition at line 49 of file SimulationDrawer.h.

Referenced by SimulationDrawer().

double evd::SimulationDrawer::minz

Definition at line 51 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().


The documentation for this class was generated from the following files: