SimulationDrawer.cxx
Go to the documentation of this file.
1 ///
2 /// \file SimulationDrawer.cxx
3 /// \brief Render the objects from the Simulation package
4 /// \author messier@indiana.edu
5 /// \version $Id: SimulationDrawer.cxx,v 1.2 2010/11/11 18:11:22 p-novaart Exp $
6 ///
7 
8 #include <iomanip>
9 #include <algorithm>
10 
11 #include "TParticle.h"
12 #include "TLatex.h"
13 #include "TPolyMarker3D.h"
14 #include "TPolyMarker.h"
15 #include "TPolyLine.h"
16 #include "TPolyLine3D.h"
17 #include "TDatabasePDG.h"
18 
20 #include "nuevdb/EventDisplayBase/View2D.h"
21 #include "nuevdb/EventDisplayBase/View3D.h"
22 #include "CoreUtils/ServiceUtil.h"
23 #include "Geometry/GeometryGAr.h"
26 #include "EventDisplay/EVD/Style.h"
30 #include "DetectorInfo/DetectorProperties.h"
31 
36 
37 namespace {
38  // Utility function to make uniform error messages.
39  void writeErrMsg(const char* fcn,
40  cet::exception const& e)
41  {
42  mf::LogWarning("SimulationDrawer") << "SimulationDrawer::" << fcn
43  << " failed with message:\n"
44  << e;
45  }
46 }
47 
48 namespace gar {
49  namespace evd{
50 
52  {
54  minx = 1e9;
55  maxx = -1e9;
56  miny = 1e9;
57  maxy = -1e9;
58  minz = 1e9;
59  maxz = -1e9;
60  double world[3] = {geom->TPCXCent(),geom->TPCYCent(),geom->TPCZCent()};
61 
62  if (minx>world[0]-geom->TPCRadius())
63  minx = world[0]-geom->TPCRadius();
64  if (maxx<world[0]+geom->TPCRadius())
65  maxx = world[0]+geom->TPCRadius();
66  if (miny>world[1]-geom->TPCRadius())
67  miny = world[1]-geom->TPCRadius();
68  if (maxy<world[1]+geom->TPCRadius())
69  maxy = world[1]+geom->TPCRadius();
70  if (minz>world[2]-geom->TPCLength()/2.)
71  minz = world[2]-geom->TPCLength()/2.;
72  if (maxz<world[2]+geom->TPCLength()/2.)
73  maxz = world[2]+geom->TPCLength()/2.;
74  }
75 
76  //......................................................................
77 
79  {
80  }
81 
82  //......................................................................
84  evdb::View2D* view)
85  {
86 
87  if( evt.isRealData() ) return;
88 
90  // Skip drawing if option is turned off
91  if (!drawopt->fShowMCTruthText) return;
92 
93  std::vector<const simb::MCTruth*> mctruth;
94  this->GetMCTruth(evt, mctruth);
95 
96  for (unsigned int i = 0; i<mctruth.size(); ++i) {
97  std::string mctext;
98  bool firstin = true;
99  bool firstout = true;
101  std::string incoming;
102  std::string outgoing;
103  // Label cosmic rays -- others are pretty obvious
104  if (mctruth[i]->Origin()==simb::kCosmicRay) origin = "c-ray: ";
105  int jmax = TMath::Min(20,mctruth[i]->NParticles());
106  for (int j=0; j<jmax; ++j) {
107  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
108  char buff[1024];
109  if (p.P()>0.05) {
110  sprintf(buff,"#color[%d]{%s #scale[0.75]{[%.1f GeV/c]}}",
113  p.P());
114  }
115  else {
116  sprintf(buff,"#color[%d]{%s}",
119  }
120  if (p.StatusCode()==0) {
121  if (firstin==false) incoming += " + ";
122  incoming += buff;
123  firstin = false;
124  }
125  if (p.StatusCode()==1) {
126  if (firstout==false) outgoing += " + ";
127  outgoing += buff;
128  firstout = false;
129  }
130  } // loop on j particles
131  if (origin=="" && incoming=="") {
132  mctext = outgoing;
133  }
134  else {
135  mctext = origin+incoming+" #rightarrow "+outgoing;
136  }
137  TLatex& latex = view->AddLatex(0.03, 0.2, mctext.c_str());
138  latex.SetTextSize(0.6);
139 
140  } // loop on i mctruth objects
141  }
142 
143  //......................................................................
145  evdb::View2D* /*view*/)
146  {
147  if( evt.isRealData() ) return;
148 
150  // Skip drawing if option is turned off
151  if (!drawopt->fShowMCTruthText) return;
152 
153  std::vector<const simb::MCTruth*> mctruth;
154  this->GetMCTruth(evt, mctruth);
155  std::cout<<"\nMCTruth Ptcl trackID PDG P T Moth Process\n";
156  for (unsigned int i=0; i<mctruth.size(); ++i) {
157  for (int j=0; j<mctruth[i]->NParticles(); ++j) {
158  const simb::MCParticle& p = mctruth[i]->GetParticle(j);
159  if(p.StatusCode() == 0 || p.StatusCode() == 1) {
160  int KE = 1000 * (p.E() - p.Mass());
161  std::cout<<std::right<<std::setw(7)<<i<<std::setw(5)<<j
162  <<std::setw(8)<<p.TrackId()
163  <<" "<<std::setw(14)<<Style::LatexName(p.PdgCode())
164  <<std::setw(7)<<int(1000 * p.P())
165  <<std::setw(7)<<KE<<std::setw(7)<<p.Mother()
166  <<" "<<p.Process()
167  <<"\n";
168  }
169  } // loop on j particles in list
170  }
171  std::cout<<"Note: Momentum, P, and kinetic energy, T, in MeV/c\n";
172  } // MCTruthLongText
173 
174 
175  //......................................................................
176  //this method draws the true particle trajectories in 3D
178  evdb::View3D* view)
179  {
180  if( evt.isRealData() ) return;
181 
183 
185  double xcent = geo->TPCXCent();
186  double ycent = geo->TPCYCent();
187  double zcent = geo->TPCZCent();
188  xcent = 0;
189  ycent = 0;
190  zcent = 0;
191 
192  // get the particles from the Geant4 step
193  std::vector<const simb::MCParticle*> plist;
194  this->GetParticle(evt, plist);
195 
196  // Useful variables
197  double xMinimum(-1.*(maxx-minx));
198  double xMaximum( 2.*(maxx-minx));
199 
200  // Define a couple of colors for neutrals and if we gray it out...
201  int neutralColor(12);
202  int grayedColor(15);
203  int neutrinoColor(38);
204 
205  // Using the voxel information can be slow (see previous implementation of this code).
206  // In order to speed things up we have modified the strategy:
207  // 1) Make one pass through the list of voxels
208  // 2) For each voxel, keep track of the MCParticle contributing energy to it and it's position
209  // which is done by keeping a map between the MCParticle and a vector of positions
210  // 3) Then loop through the map to draw the particle trajectories.
211  // One caveat is the need for MCParticles... and the voxels contain the track ids. So we'll need one
212  // more loop to make a map of track id's and MCParticles.
213 
214  // First up is to build the map between track id's and associated MCParticles so we can recover when looping over voxels
215  std::map<int, const simb::MCParticle*> trackToMcParticleMap;
216 
217  // Should we display the trajectories too?
218  double minPartEnergy(0.01);
219 
220  std::cout << "MC particle listing (limit 20)" << std::endl;
221  std::cout << " ID PDG Px Py Pz E (GeV) origin (cm) " << std::endl;
222 
223  for(size_t p = 0; p < plist.size(); ++p){
224 
225  trackToMcParticleMap[plist[p]->TrackId()] = plist[p];
226 
227  const simb::MCParticle* mcPart = plist[p];
228  const simb::MCTrajectory& mcTraj = mcPart->Trajectory();
229 
230  int pdgCode(mcPart->PdgCode());
231  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
232  TParticlePDG* partPDG(TDatabasePDG::Instance()->GetParticle(pdgCode));
233  double partCharge = partPDG ? partPDG->Charge() : 0.;
234  double partEnergy = mcPart->E();
235 
236  if (p<20)
237  {
238  printf("%6d %12d %8.3f %8.3f %8.3f %8.3f ",(int) p,pdgCode,mcPart->Px(),mcPart->Py(),mcPart->Pz(),partEnergy);
239  if (!mcTraj.empty())
240  {
241  double xPos = mcTraj.X(0) + xcent;
242  double yPos = mcTraj.Y(0) + ycent;
243  double zPos = mcTraj.Z(0) + zcent;
244  printf("%8.3f %8.3f %8.3f\n",xPos, yPos, zPos);
245  }
246  else
247  {
248  printf("\n");
249  }
250  }
251 
252  // If the option is turned off, stop here after printing it out
253  if (!drawopt->fShowMCTruthTrajectories) continue;
254 
255  // cut out all neutrals if asked to do so
256  if ( !drawopt->fShowNeutrals && partCharge > -0.1 && partCharge < 0.1) continue;
257 
258  // energy cut on photons
259  if (pdgCode == 22 && partEnergy < drawopt->fPhotonEnergyCut) continue;
260 
261  // energy cut on neutrons
262  if (pdgCode == 2112 && partEnergy < drawopt->fNeutronEnergyCut) continue;
263 
264  // energy cut on all other neutral particles, not photons or neutrons
265  if ( partEnergy < drawopt->fOtherNeutralEnergyCut && partCharge > -0.1 && partCharge < 0.1
266  && pdgCode != 22 && pdgCode != 2112) continue;
267 
268  if (!drawopt->fShowMCTruthColors) colorIdx = grayedColor;
269 
270  if (pdgCode >= 1000000000) continue;
271  if (pdgCode == 111) continue; // don't draw pizeros
272 
273  //if ( partCharge > -0.1 && partCharge < 0.1 )
274  //{
275  // std::cout << "Drawing neutral particle: " << pdgCode << " " << partEnergy << std::endl;
276  //}
277 
278 
279  if (!mcTraj.empty() && partEnergy > minPartEnergy && mcPart->TrackId() < 100000000){
280  // The following is meant to get the correct offset for drawing the particle trajectory
281  // In particular, the cosmic rays will not be correctly placed without this
282 
283 
284  // collect the points from this particle
285  int numTrajPoints = mcTraj.size();
286 
287  std::unique_ptr<double[]> hitPositions(new double[3*numTrajPoints]);
288  int hitCount(0);
289 
290  for(int hitIdx = 0; hitIdx < numTrajPoints; hitIdx++){
291  double xPos = mcTraj.X(hitIdx) + xcent;
292  double yPos = mcTraj.Y(hitIdx) + ycent;
293  double zPos = mcTraj.Z(hitIdx) + zcent;
294 
295  // If the original simulated hit did not occur in the TPC volume then don't draw it
296  if (xPos < minx || xPos > maxx || yPos < miny || yPos > maxy|| zPos < minz || zPos > maxz) continue;
297 
298  // Check fiducial limits
299  // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y) to the polyline coords
300  if (xPos > xMinimum && xPos < xMaximum){
301  hitPositions[3*hitCount ] = zPos;
302  hitPositions[3*hitCount + 1] = xPos;
303  hitPositions[3*hitCount + 2] = yPos;
304  hitCount++;
305  }
306  }
307 
308  TPolyLine3D& pl(view->AddPolyLine3D(1, colorIdx, 1, 1));
309 
310  // Draw neutrals as a gray dotted line to help fade into background a bit...
311  if (partCharge == 0.){
312  pl.SetLineColor(neutralColor);
313  pl.SetLineStyle(3);
314  pl.SetLineWidth(1);
315  }
316  pl.SetPolyLine(hitCount, hitPositions.get(), "");
317  }
318  }
319 
320  // Now we set up and build the map between MCParticles and a vector of positions obtained from the voxels
321  std::map<const simb::MCParticle*, std::vector<std::vector<double> > > partToPosMap;
322 
323  // Finally ready for the main event! Simply loop through the map between MCParticle and positions to
324  // draw the trajectories
325  std::map<const simb::MCParticle*, std::vector<std::vector<double> > >::iterator partToPosMapItr;
326 
327  for(partToPosMapItr = partToPosMap.begin(); partToPosMapItr != partToPosMap.end(); ++partToPosMapItr)
328  {
329  // Recover the McParticle, we'll need to access several data members so may as well dereference it
330  const simb::MCParticle* mcPart = partToPosMapItr->first;
331 
332  // Apparently, it can happen that we get a null pointer here or maybe no points to plot
333  if (!mcPart || partToPosMapItr->second.empty()) continue;
334 
335  // The following is meant to get the correct offset for drawing the particle trajectory
336  // In particular, the cosmic rays will not be correctly placed without this
337 
338  int colorIdx(evd::Style::ColorFromPDG(mcPart->PdgCode()));
339  int markerIdx(kFullDotSmall);
340  int markerSize(2);
341 
342  if (!drawopt->fShowMCTruthFullSize){
343  colorIdx = grayedColor;
344  markerIdx = kDot;
345  markerSize = 1;
346  }
347 
348  TPolyMarker3D& pm = view->AddPolyMarker3D(partToPosMapItr->second.size(), colorIdx, markerIdx, markerSize);
349 
350  // Now loop over points and add to trajectory
351  for(size_t posIdx = 0; posIdx < partToPosMapItr->second.size(); posIdx++){
352  const std::vector<double>& posVec = partToPosMapItr->second[posIdx];
353 
354  double xCoord = posVec[0];
355 
356  // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y)
357  if (xCoord > xMinimum && xCoord < xMaximum)
358  pm.SetPoint(posIdx, posVec[2], xCoord, posVec[1]);
359  }
360  }
361 
362  // Finally, let's see if we can draw the incoming particle from the MCTruth information
363  std::vector<const simb::MCTruth*> mctruth;
364  this->GetMCTruth(evt, mctruth);
365 
366  // Loop through the MCTruth vector
367  for (unsigned int idx = 0; idx < mctruth.size(); idx++){
368  // Go through each MCTruth object in the list
369  for (int particleIdx = 0; particleIdx < mctruth[idx]->NParticles(); particleIdx++){
370  const simb::MCParticle& mcPart = mctruth[idx]->GetParticle(particleIdx);
371 
372  // A negative mother id indicates the "primary" particle
373  if(mcPart.Mother() == -1 && mcPart.StatusCode() == 0){
374  mf::LogDebug("SimulationDrawer") << mcPart << std::endl;
375 
376  // Get position vector
377  TVector3 particlePosition(mcPart.Vx(),mcPart.Vy(),mcPart.Vz());
378 
379  // Get direction vector (in opposite direction)
380  TVector3 oppPartDir(-mcPart.Px(),-mcPart.Py(),-mcPart.Pz());
381 
382  if (oppPartDir.Mag2() > 0.) oppPartDir.SetMag(1.);
383 
384  //double arcLenToDraw = -particlePosition.Z() / oppPartDir.CosTheta();
385  double arcLenToDraw = 50.0; // always draw 50 cm of neutrinos
386 
387  // No point in drawing if arc length is zero (e.g. Ar nucleus)
388  if (arcLenToDraw > 0.){
389  // Draw the line, use an off color to be unique
390  TPolyLine3D& pl(view->AddPolyLine3D(2, neutrinoColor, 1, 2));
391 
392  // nb. coordinates are so Y is up, but ROOT thinks Z is up, so report (z,x,y)
393 
394  pl.SetPoint(0,particlePosition.Z()+zcent,particlePosition.X()+xcent,particlePosition.Y()+ycent);
395 
396  particlePosition += std::min(arcLenToDraw + 10.,1000.) * oppPartDir;
397 
398  pl.SetPoint(1,particlePosition.Z()+zcent,particlePosition.X()+xcent,particlePosition.Y()+ycent);
399  }
400  }
401  // The particles we want to draw will be early in the list so break out if we didn't find them
402  else break;
403  } // loop on particles in list
404  }
405 
406  return;
407  }
408 
409  //......................................................................
411  std::vector<const simb::MCParticle*>& plist)
412  {
413  plist.clear();
414 
415  if( evt.isRealData() ) return 0;
416 
418 
419  std::vector<const simb::MCParticle*> temp;
420 
422  // use get by Type because there should only be one collection of these in the event
423  try{
424  evt.getView(drawopt->fG4ModuleLabel, plcol);
425  for(unsigned int i = 0; i < plcol.vals().size(); ++i){
426  temp.push_back(plcol.vals().at(i));
427  }
428  temp.swap(plist);
429  }
430  catch(cet::exception& e){
431  writeErrMsg("GetRawDigits", e);
432  }
433 
434  return plist.size();
435 
436  }
437 
438  //......................................................................
440  std::vector<const simb::MCTruth*>& mcvec)
441  {
442  mcvec.clear();
443 
444  if( evt.isRealData() ) return 0;
445 
446  std::vector<const simb::MCTruth*> temp;
447 
448  // use get by Type because there should only be one collection of these in the event
449  try{
450  auto mctcol = evt.getMany<std::vector<simb::MCTruth> >();
451  for(size_t mctc = 0; mctc < mctcol.size(); ++mctc){
452  art::Handle< std::vector<simb::MCTruth> > mclistHandle = mctcol[mctc];
453 
454  for(size_t i = 0; i < mclistHandle->size(); ++i){
455  temp.push_back(&(mclistHandle->at(i)));
456  }
457  }
458  temp.swap(mcvec);
459  }
460  catch(cet::exception& e){
461  writeErrMsg("GetMCTruth", e);
462  }
463 
464  return mcvec.size();
465  }
466 
467  //......................................................................
468  void SimulationDrawer::HiLite(int trkId, bool dohilite)
469  {
470  fHighlite[trkId] = dohilite;
471  }
472  }
473 }// namespace
474 ////////////////////////////////////////////////////////////////////////
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
std::map< int, bool > fHighlite
int GetMCTruth(const art::Event &evt, std::vector< const simb::MCTruth * > &mctruth)
double Py(const int i=0) const
Definition: MCParticle.h:231
void MCTruth3D(const art::Event &evt, evdb::View3D *view)
static const char * LatexName(int pdgcode)
Convert PDG code to a latex string (root-style)
Definition: Style.cxx:14
std::string string
Definition: nybbler.cc:12
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
int Mother() const
Definition: MCParticle.h:213
double Mass() const
Definition: MCParticle.h:239
double Px(const int i=0) const
Definition: MCParticle.h:230
art::InputTag fG4ModuleLabel
module label producing sim::SimChannel objects
int StatusCode() const
Definition: MCParticle.h:211
void MCTruthLongText(const art::Event &evt, evdb::View2D *view)
std::string Process() const
Definition: MCParticle.h:215
Particle class.
void HiLite(int trkId, bool hlt=true)
bool empty() const
Definition: MCTrajectory.h:167
int TrackId() const
Definition: MCParticle.h:210
auto & vals() noexcept
Definition: View.h:68
bool isRealData() const
void MCTruthShortText(const art::Event &evt, evdb::View2D *view)
const double e
double Y(const size_type i) const
Definition: MCTrajectory.h:150
LArSoft includes.
Definition: InfoTransfer.h:33
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
double P(const int i=0) const
Definition: MCParticle.h:234
static int ColorFromPDG(int pdgcode)
Definition: Style.cxx:67
p
Definition: test.py:223
Definition: fwd.h:46
int GetParticle(const art::Event &evt, std::vector< const simb::MCParticle * > &plist)
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
General GArSoft Utilities.
double Vx(const int i=0) const
Definition: MCParticle.h:221
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
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
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:500
TCEvent evt
Definition: DataStructs.cxx:7
LArSoft geometry interface.
Definition: ChannelGeo.h:16
art framework interface to geometry description
double Vy(const int i=0) const
Definition: MCParticle.h:222
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Cosmic rays.
Definition: MCTruth.h:24
QTextStream & endl(QTextStream &s)