LArSoft  v10_04_05
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 27 of file SimulationDrawer.h.

Constructor & Destructor Documentation

evd::SimulationDrawer::SimulationDrawer ( )

Definition at line 56 of file SimulationDrawer.cxx.

References maxx, maxy, maxz, minx, miny, and minz.

57  {
58  // For now only draw cryostat=0.
60  minx = 1e9;
61  maxx = -1e9;
62  miny = 1e9;
63  maxy = -1e9;
64  minz = 1e9;
65  maxz = -1e9;
66 
67  // This is looking to find the range of the complete active volume... this may not be the
68  // best way to do this...
69  for (auto const& tpc : geom->Iterate<geo::TPCGeo>()) {
70  auto const tpc_center = tpc.Center();
71  auto const active_volume_center = tpc.GetActiveVolumeCenter();
72  mf::LogDebug("SimulationDrawer") << tpc.ID() << ", TPC center: " << tpc_center.X() << ", "
73  << tpc_center.Y() << ", " << tpc_center.Z() << std::endl;
74  mf::LogDebug("SimulationDrawer")
75  << " TPC Active center: " << active_volume_center.X() << ", "
76  << active_volume_center.Y() << ", " << active_volume_center.Z()
77  << ", H/W/L: " << tpc.ActiveHalfHeight() << "/" << tpc.ActiveHalfWidth() << "/"
78  << tpc.ActiveLength() << std::endl;
79 
80  if (minx > tpc_center.X() - tpc.HalfWidth()) minx = tpc_center.X() - tpc.HalfWidth();
81  if (maxx < tpc_center.X() + tpc.HalfWidth()) maxx = tpc_center.X() + tpc.HalfWidth();
82  if (miny > tpc_center.Y() - tpc.HalfHeight()) miny = tpc_center.Y() - tpc.HalfHeight();
83  if (maxy < tpc_center.Y() + tpc.HalfHeight()) maxy = tpc_center.Y() + tpc.HalfHeight();
84  if (minz > tpc_center.Z() - tpc.Length() / 2.) minz = tpc_center.Z() - tpc.Length() / 2.;
85  if (maxz < tpc_center.Z() + tpc.Length() / 2.) maxz = tpc_center.Z() + tpc.Length() / 2.;
86 
87  mf::LogDebug("SimulationDrawer")
88  << " minx/maxx: " << minx << "/" << maxx << ", miny/maxy: " << miny << "/" << maxy
89  << ", minz/miny: " << minz << "/" << maxz << std::endl;
90  }
91  }
Geometry information for a single TPC.
Definition: TPCGeo.h:33
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
evd::SimulationDrawer::~SimulationDrawer ( )

Definition at line 95 of file SimulationDrawer.cxx.

95 {}

Member Function Documentation

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

Definition at line 914 of file SimulationDrawer.cxx.

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

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

915  {
916  mcvec.clear();
917 
918  if (evt.isRealData()) return 0;
919 
920  std::vector<const simb::MCTruth*> temp;
921 
922  std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
923  // use get by Type because there should only be one collection of these in the event
924  try {
925  //evt.getManyByType(mctcol);
926  mctcol = evt.getMany<std::vector<simb::MCTruth>>();
927  for (size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
928  art::Handle<std::vector<simb::MCTruth>> mclistHandle = mctcol[mctc];
929 
930  for (size_t i = 0; i < mclistHandle->size(); ++i) {
931  temp.push_back(&(mclistHandle->at(i)));
932  }
933  }
934  temp.swap(mcvec);
935  }
936  catch (cet::exception& e) {
937  writeErrMsg("GetMCTruth", e);
938  }
939 
940  return mcvec.size();
941  }
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 885 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().

887  {
888  plist.clear();
889 
890  if (evt.isRealData()) return 0;
891 
893 
894  std::vector<const simb::MCParticle*> temp;
895 
897  // use get by Type because there should only be one collection of these in the event
898  try {
899  evt.getView(drawopt->fG4ModuleLabel, plcol);
900  for (unsigned int i = 0; i < plcol.vals().size(); ++i) {
901  temp.push_back(plcol.vals().at(i));
902  }
903  temp.swap(plist);
904  }
905  catch (cet::exception& e) {
906  writeErrMsg("GetRawDigits", e);
907  }
908 
909  return plist.size();
910  }
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 945 of file SimulationDrawer.cxx.

References fHighlite.

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

Definition at line 328 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().

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

161  {
162  if (evt.isRealData()) return;
163 
165  // Skip drawing if option is turned off
166  if (!drawopt->fShowMCTruthText) return;
167 
168  std::vector<const simb::MCTruth*> mctruth;
169  this->GetMCTruth(evt, mctruth);
170  std::cout << "\nMCTruth Ptcl trackID PDG P T Moth Process\n";
171  for (unsigned int i = 0; i < mctruth.size(); ++i) {
172  for (int j = 0; j < mctruth[i]->NParticles(); ++j) {
173  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
174  if (p.StatusCode() == 0 || p.StatusCode() == 1) {
175  int KE = 1000 * (p.E() - p.Mass());
176  std::cout << std::right << std::setw(7) << i << std::setw(5) << j << std::setw(8)
177  << p.TrackId() << " " << std::setw(14) << Style::LatexName(p.PdgCode())
178  << std::setw(7) << int(1000 * p.P()) << std::setw(7) << KE << std::setw(7)
179  << p.Mother() << " " << p.Process() << "\n";
180  }
181  } // loop on j particles in list
182  }
183  std::cout << "Note: Momentum, P, and kinetic energy, T, in MeV/c\n";
184  } // 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 597 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().

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

Definition at line 99 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().

100  {
101 
102  if (evt.isRealData()) return;
103 
105  // Skip drawing if option is turned off
106  if (!drawopt->fShowMCTruthText) return;
107 
108  std::vector<const simb::MCTruth*> mctruth;
109  this->GetMCTruth(evt, mctruth);
110 
111  for (unsigned int i = 0; i < mctruth.size(); ++i) {
112  std::string mctext;
113  bool firstin = true;
114  bool firstout = true;
115  std::string origin;
116  std::string incoming;
117  std::string outgoing;
118  // Label cosmic rays -- others are pretty obvious
119  if (mctruth[i]->Origin() == simb::kCosmicRay) origin = "c-ray: ";
120  int jmax = TMath::Min(20, mctruth[i]->NParticles());
121  for (int j = 0; j < jmax; ++j) {
122  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
123  char buff[1024];
124  if (p.P() > 0.05) {
125  sprintf(buff,
126  "#color[%d]{%s #scale[0.75]{[%.1f GeV/c]}}",
129  p.P());
130  }
131  else {
132  sprintf(buff,
133  "#color[%d]{%s}",
136  }
137  if (p.StatusCode() == 0) {
138  if (firstin == false) incoming += " + ";
139  incoming += buff;
140  firstin = false;
141  }
142  if (p.StatusCode() == 1) {
143  if (firstout == false) outgoing += " + ";
144  outgoing += buff;
145  firstout = false;
146  }
147  } // loop on j particles
148  if (origin == "" && incoming == "") { mctext = outgoing; }
149  else {
150  mctext = origin + incoming + " #rightarrow " + outgoing;
151  }
152  TLatex& latex = view->AddLatex(0.03, 0.2, mctext.c_str());
153  latex.SetTextSize(0.6);
154 
155  } // loop on i mctruth objects
156  }
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 188 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, Get, 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(), and simb::MCParticle::Vz().

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

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

Referenced by HiLite().

double evd::SimulationDrawer::maxx

Definition at line 46 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::maxy

Definition at line 48 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::maxz

Definition at line 50 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::minx

Definition at line 45 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().

double evd::SimulationDrawer::miny

Definition at line 47 of file SimulationDrawer.h.

Referenced by SimulationDrawer().

double evd::SimulationDrawer::minz

Definition at line 49 of file SimulationDrawer.h.

Referenced by MCTruthOrtho(), and SimulationDrawer().


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