DrawLArVoxel3D_tool.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file DrawLArVoxel3D_tool.cc
3 /// \author T. Usher
4 ////////////////////////////////////////////////////////////////////////
5 
8 
12 #include "cetlib_except/exception.h"
14 
21 
26 
27 #include "TDatabasePDG.h"
28 #include "TPolyLine3D.h"
29 #include "TPolyMarker3D.h"
30 
31 namespace evdb_tool {
32 
33  class DrawLArVoxel3D : public ISim3DDrawer {
34  public:
35  explicit DrawLArVoxel3D(const fhicl::ParameterSet&);
36 
37  void Draw(const art::Event&, evdb::View3D*) const override;
38 
39  private:
40  int GetMCTruth(const art::Event&, std::vector<const simb::MCTruth*>&) const;
41  };
42 
43  //----------------------------------------------------------------------
44  // Constructor.
46 
47  void
48  DrawLArVoxel3D::Draw(const art::Event& evt, evdb::View3D* view) const
49  {
51 
52  // If the option is turned off, there's nothing to do
53  if (!drawOpt->fShowSimChannelInfo) return;
54 
55  // If the option is turned off, there's nothing to do
56  if (!drawOpt->fShowMCTruthTrajectories) return;
57 
58  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
59  auto const detProp =
62 
63  // Recover a handle to the collection of MCParticles
65 
66  evt.getByLabel(drawOpt->fG4ModuleLabel, mcParticleHandle);
67 
68  if (!mcParticleHandle.isValid()) return;
69 
70  // Define a couple of colors for neutrals and if we gray it out...
71  int neutralColor(12);
72  int grayedColor(15);
73  int neutrinoColor(38);
74 
75  // Use the LArVoxelList to get the true energy deposition locations as opposed to using MCTrajectories
76  const sim::LArVoxelList voxels =
78 
79  mf::LogDebug("SimulationDrawer")
80  << "Starting loop over " << mcParticleHandle->size() << " McParticles, voxel list size is "
81  << voxels.size() << std::endl;
82 
83  // Using the voxel information can be slow (see previous implementation of this code).
84  // In order to speed things up we have modified the strategy:
85  // 1) Make one pass through the list of voxels
86  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
87  // which is done by keeping a map between the MCParticle and a vector of positions
88  // 3) Then loop through the map to draw the particle trajectories.
89  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
90  // more loop to make a map of track id's and MCParticles.
91 
92  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
93  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
94 
95  // Should we display the trajectories too?
96  double minPartEnergy(0.01);
97 
98  for (size_t p = 0; p < mcParticleHandle->size(); ++p) {
99  art::Ptr<simb::MCParticle> mcParticle(mcParticleHandle, p);
100 
101  trackToMcParticleMap[mcParticle->TrackId()] = mcParticle.get();
102 
103  // Quick loop through to draw trajectories...
104  if (drawOpt->fShowMCTruthTrajectories) {
105  // Is there an associated McTrajectory?
106  const simb::MCTrajectory& mcTraj = mcParticle->Trajectory();
107 
108  int pdgCode(mcParticle->PdgCode());
109  int colorIdx(evd::Style::ColorFromPDG(mcParticle->PdgCode()));
110  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
111  double partCharge = partPDG ? partPDG->Charge() : 0.;
112  double partEnergy = mcParticle->E();
113 
114  if (!drawOpt->fShowMCTruthColors) colorIdx = grayedColor;
115 
116  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcParticle->TrackId() < 100000000) {
117  double g4Ticks(clockData.TPCG4Time2Tick(mcParticle->T()) - trigger_offset(clockData));
118  double xOffset(0.);
119  double xPosMinTick = 0.;
120  double xPosMaxTick = std::numeric_limits<double>::max();
121 
122  // collect the points from this particle
123  int numTrajPoints = mcTraj.size();
124 
125  std::unique_ptr<double[]> hitPositions(new double[3 * numTrajPoints]);
126  int hitCount(0);
127 
128  for (int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++) {
129  double xPos = mcTraj.X(hitIdx);
130  double yPos = mcTraj.Y(hitIdx);
131  double zPos = mcTraj.Z(hitIdx);
132 
133  // If we have cosmic rays then we need to get the offset which allows translating from
134  // when they were generated vs when they were tracked.
135  // Note that this also explicitly checks that they are in a TPC volume
136  geo::Point_t hitLocation(xPos, yPos, zPos);
137 
138  try {
139  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
140  geo::PlaneID planeID(tpcID, 0);
141 
142  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
143  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
144  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
145 
146  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
147  }
148  catch (...) {
149  continue;
150  }
151 
152  // Now move the hit position to correspond to the timing
153  xPos += xOffset;
154 
155  // Check fiducial limits
156  if (xPos > xPosMinTick && xPos < xPosMaxTick) {
157  hitPositions[3 * hitCount] = xPos;
158  hitPositions[3 * hitCount + 1] = yPos;
159  hitPositions[3 * hitCount + 2] = zPos;
160  hitCount++;
161  }
162  }
163 
164  TPolyLine3D& pl(view->AddPolyLine3D(1, colorIdx, 1, 1));
165 
166  // Draw neutrals as a gray dotted line to help fade into background a bit...
167  if (partCharge == 0.) {
168  pl.SetLineColor(neutralColor);
169  pl.SetLineStyle(3);
170  pl.SetLineWidth(1);
171  }
172  pl.SetPolyLine(hitCount, hitPositions.get(), "");
173  }
174  }
175  }
176 
177  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
178  std::map<const simb::MCParticle*, std::vector<std::vector<double>>> partToPosMap;
179 
181  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++) {
182  const sim::LArVoxelData& vxd = (*vxitr).second;
183 
184  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++) {
185  if (vxd.Energy(partIdx) > drawOpt->fMinEnergyDeposition) {
186  int trackId = vxd.TrackID(partIdx);
187 
188  // It can be in some instances that mcPart here could be zero.
189  const simb::MCParticle* mcPart = trackToMcParticleMap[trackId];
190 
191  partToPosMap[mcPart].push_back(std::vector<double>(3));
192 
193  partToPosMap[mcPart].back()[0] = vxd.VoxelID().X();
194  partToPosMap[mcPart].back()[1] = vxd.VoxelID().Y();
195  partToPosMap[mcPart].back()[2] = vxd.VoxelID().Z();
196  }
197  } // end if this track id is in the current voxel
198  } // end loop over voxels
199 
200  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
201  // draw the trajectories
202  std::map<const simb::MCParticle*, std::vector<std::vector<double>>>::iterator partToPosMapItr;
203 
204  for (partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end();
205  partToPosMapItr++) {
206  // Recover the McParticle, we'll need to access several data members so may as well dereference it
207  const simb::MCParticle* mcPart = partToPosMapItr->first;
208 
209  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
210  if (!mcPart || partToPosMapItr->second.empty()) continue;
211 
212  double g4Ticks(clockData.TPCG4Time2Tick(mcPart->T()) - trigger_offset(clockData));
213  double xOffset = 0.;
214  double xPosMinTick = 0.;
215  double xPosMaxTick = std::numeric_limits<double>::max();
216 
217  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
218  int markerIdx(kFullDotSmall);
219  int markerSize(2);
220 
221  if (!drawOpt->fShowMCTruthFullSize) {
222  colorIdx = grayedColor;
223  markerIdx = kDot;
224  markerSize = 1;
225  }
226 
227  std::unique_ptr<double[]> hitPositions(new double[3 * partToPosMapItr->second.size()]);
228  int hitCount(0);
229 
230  // Now loop over points and add to trajectory
231  for (size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++) {
232  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
233 
234  // Check xOffset state and set if necessary
235  geo::Point_t hitLocation(posVec[0], posVec[1], posVec[2]);
236 
237  try {
238  geo::TPCID tpcID = geom->PositionToTPCID(hitLocation);
239  geo::PlaneID planeID(tpcID, 0);
240 
241  xPosMinTick = detProp.ConvertTicksToX(0, planeID);
242  xPosMaxTick = detProp.ConvertTicksToX(detProp.NumberTimeSamples(), planeID);
243  xOffset = detProp.ConvertTicksToX(g4Ticks, planeID) - xPosMinTick;
244 
245  if (xPosMaxTick < xPosMinTick) std::swap(xPosMinTick, xPosMaxTick);
246  }
247  catch (...) {
248  continue;
249  }
250 
251  double xCoord = posVec[0] + xOffset;
252 
253  // If a voxel records an energy deposit then must have been in the TPC
254  // But because things get shifted still need to cut off if outside drift
255  if (xCoord > xPosMinTick && xCoord < xPosMaxTick) {
256  hitPositions[3 * hitCount] = xCoord;
257  hitPositions[3 * hitCount + 1] = posVec[1];
258  hitPositions[3 * hitCount + 2] = posVec[2];
259  hitCount++;
260  }
261  }
262 
263  TPolyMarker3D& pm = view->AddPolyMarker3D(1, colorIdx, markerIdx, markerSize);
264  pm.SetPolyMarker(hitCount, hitPositions.get(), markerIdx);
265  }
266 
267  // Finally, let's see if we can draw the incoming particle from the MCTruth information
268  std::vector<const simb::MCTruth*> mctruth;
269  this->GetMCTruth(evt, mctruth);
270 
271  // Loop through the MCTruth vector
272  for (unsigned int idx = 0; idx < mctruth.size(); idx++) {
273  // Go through each MCTruth object in the list
274  for (int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++) {
275  const simb::MCParticle& mcPart = mctruth[idx]->GetParticle(particleIdx);
276 
277  // A negative mother id indicates the "primary" particle
278  if (mcPart.Mother() == -1 && mcPart.StatusCode() == 0) {
279  mf::LogDebug("SimulationDrawer") << mcPart << std::endl;
280 
281  // Get position vector
282  TVector3 particlePosition(mcPart.Vx(), mcPart.Vy(), mcPart.Vz());
283 
284  // Get direction vector (in opposite direction)
285  TVector3 oppPartDir(-mcPart.Px(), -mcPart.Py(), -mcPart.Pz());
286 
287  if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
288 
289  double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
290 
291  // No point in drawing if arc length is zero (e.g. Ar nucleus)
292  if (arcLenToDraw > 0.) {
293  // Draw the line, use an off color to be unique
294  TPolyLine3D& pl(view->AddPolyLine3D(2, neutrinoColor, 1, 2));
295 
296  pl.SetPoint(0, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
297 
298  particlePosition += std::min(arcLenToDraw + 10., 1000.) * oppPartDir;
299 
300  pl.SetPoint(1, particlePosition.X(), particlePosition.Y(), particlePosition.Z());
301  }
302  }
303  // The particles we want to draw will be early in the list so break out if we didn't find them
304  else
305  break;
306  } // loop on particles in list
307  }
308 
309  return;
310  }
311 
312  //......................................................................
313  int
314  DrawLArVoxel3D::GetMCTruth(const art::Event& evt, std::vector<const simb::MCTruth*>& mcvec) const
315  {
316  mcvec.clear();
317 
318  if (evt.isRealData()) return 0;
319 
320  std::vector<const simb::MCTruth*> temp;
321 
322  std::vector<art::Handle<std::vector<simb::MCTruth>>> mctcol;
323 
324  // use get by Type because there should only be one collection of these in the event
325  try {
326  //evt.getManyByType(mctcol);
327  mctcol = evt.getMany<std::vector<simb::MCTruth>>();
328 
329  for (size_t mctc = 0; mctc < mctcol.size(); ++mctc) {
330  art::Handle<std::vector<simb::MCTruth>> mclistHandle = mctcol[mctc];
331 
332  for (size_t i = 0; i < mclistHandle->size(); ++i) {
333  temp.push_back(&(mclistHandle->at(i)));
334  }
335  }
336  temp.swap(mcvec);
337  }
338  catch (cet::exception& e) {
339  mf::LogWarning("DrawLArVoxel3D") << "GetMCTruth:"
340  << " failed with message:\n"
341  << e;
342  }
343 
344  return mcvec.size();
345  }
346 
348 }
double E(const int i=0) const
Definition: MCParticle.h:233
double Z(const size_type i) const
Definition: MCTrajectory.h:151
double X(const size_type i) const
Definition: MCTrajectory.h:149
int PdgCode() const
Definition: MCParticle.h:212
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
double Py(const int i=0) const
Definition: MCParticle.h:231
size_type size() const
Definition: LArVoxelList.h:140
Container of LAr voxel information.
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
int Mother() const
Definition: MCParticle.h:213
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
double Px(const int i=0) const
Definition: MCParticle.h:230
static sim::LArVoxelList GetLArVoxelList(const art::Event &evt, std::string moduleLabel)
art::InputTag fG4ModuleLabel
module label producing sim::SimChannel objects
int StatusCode() const
Definition: MCParticle.h:211
geo::TPCID PositionToTPCID(geo::Point_t const &point) const
Returns the ID of the TPC at specified location.
Particle class.
std::string encode() const
Definition: InputTag.cc:97
bool empty() const
Definition: MCTrajectory.h:167
art framework interface to geometry description
iterator begin()
Definition: LArVoxelList.h:131
int TrackId() const
Definition: MCParticle.h:210
list_type::const_iterator const_iterator
Definition: LArVoxelList.h:79
bool isValid() const noexcept
Definition: Handle.h:191
bool isRealData() const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
const double e
double Y(const size_type i) const
Definition: MCTrajectory.h:150
Encapsulates the information we want store for a voxel.
void swap(Handle< T > &a, Handle< T > &b)
double Y() const
Definition: LArVoxelID.cxx:79
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:65
double T(const int i=0) const
Definition: MCParticle.h:224
p
Definition: test.py:223
size_type NumberParticles() const
Definition: LArVoxelData.h:196
static int max(int a, int b)
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double Z() const
Definition: LArVoxelID.cxx:86
DrawLArVoxel3D(const fhicl::ParameterSet &)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
double X() const
Definition: LArVoxelID.cxx:72
double Vx(const int i=0) const
Definition: MCParticle.h:221
mapped_type Energy() const
Definition: LArVoxelData.h:203
art::InputTag fSimChannelLabel
SimChannels may be independent of MC stuff.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
size_type size() const
Definition: MCTrajectory.h:166
void Draw(const art::Event &, evdb::View3D *) const override
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Pz(const int i=0) const
Definition: MCParticle.h:232
double Vz(const int i=0) const
Definition: MCParticle.h:223
int trigger_offset(DetectorClocksData const &data)
sim::LArVoxelID VoxelID() const
Definition: LArVoxelData.h:195
TCEvent evt
Definition: DataStructs.cxx:7
T const * get() const
Definition: Ptr.h:149
int GetMCTruth(const art::Event &, std::vector< const simb::MCTruth * > &) const
double Vy(const int i=0) const
Definition: MCParticle.h:222
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
const key_type & TrackID(const size_type) const