All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RecoBaseDrawer.cxx
Go to the documentation of this file.
1 /// \file RecoBaseDrawer.cxx
2 /// \brief Class to aid in the rendering of RecoBase objects
3 /// \author brebel@fnal.gov
4 
5 #include <cmath>
6 #include <limits>
7 #include <map>
8 #include <stdint.h>
9 
10 #include "TBox.h"
11 #include "TH1.h"
12 #include "TLine.h"
13 #include "TMarker.h"
14 #include "TPolyLine.h"
15 #include "TPolyLine3D.h"
16 #include "TPolyMarker.h"
17 #include "TPolyMarker3D.h"
18 #include "TRotation.h"
19 #include "TText.h"
20 #include "TVector3.h"
21 
22 #include "larcore/Geometry/Geometry.h"
23 #include "larcorealg/Geometry/CryostatGeo.h"
24 #include "larcorealg/Geometry/Exceptions.h"
25 #include "larcorealg/Geometry/PlaneGeo.h"
26 #include "larcorealg/Geometry/TPCGeo.h"
27 #include "larcorealg/Geometry/WireGeo.h"
28 #include "lardata/ArtDataHelper/MVAReader.h"
29 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
30 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
31 #include "lardata/RecoBaseProxy/Track.h"
32 #include "lardataobj/AnalysisBase/CosmicTag.h"
33 #include "lardataobj/RecoBase/Cluster.h"
34 #include "lardataobj/RecoBase/Edge.h"
35 #include "lardataobj/RecoBase/EndPoint2D.h"
36 #include "lardataobj/RecoBase/Event.h"
37 #include "lardataobj/RecoBase/Hit.h"
38 #include "lardataobj/RecoBase/OpFlash.h"
39 #include "lardataobj/RecoBase/PCAxis.h"
40 #include "lardataobj/RecoBase/PFParticle.h"
41 #include "lardataobj/RecoBase/Seed.h"
42 #include "lardataobj/RecoBase/Shower.h"
43 #include "lardataobj/RecoBase/SpacePoint.h"
44 #include "lardataobj/RecoBase/Track.h"
45 #include "lardataobj/RecoBase/Vertex.h"
46 #include "lardataobj/RecoBase/Wire.h"
47 #include "lareventdisplay/EventDisplay/3DDrawers/ISpacePoints3D.h"
48 #include "lareventdisplay/EventDisplay/ColorDrawingOptions.h"
49 #include "lareventdisplay/EventDisplay/RawDrawingOptions.h"
50 #include "lareventdisplay/EventDisplay/RecoBaseDrawer.h"
51 #include "lareventdisplay/EventDisplay/RecoDrawingOptions.h"
52 #include "lareventdisplay/EventDisplay/eventdisplay.h"
53 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
54 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
55 #include "nuevdb/EventDisplayBase/EventHolder.h"
56 #include "nuevdb/EventDisplayBase/View2D.h"
57 #include "nuevdb/EventDisplayBase/View3D.h"
58 
59 #include "art/Framework/Principal/Event.h"
60 #include "art/Framework/Principal/Handle.h"
61 #include "art/Framework/Services/Registry/ServiceHandle.h"
62 #include "art/Utilities/make_tool.h"
63 #include "canvas/Persistency/Common/FindMany.h"
64 #include "canvas/Persistency/Common/Ptr.h"
65 #include "canvas/Persistency/Common/PtrVector.h"
66 #include "cetlib_except/exception.h"
67 #include "messagefacility/MessageLogger/MessageLogger.h"
68 
69 namespace {
70  // Utility function to make uniform error messages.
71  void
72  writeErrMsg(const char* fcn, cet::exception const& e)
73  {
74  mf::LogWarning("RecoBaseDrawer") << "RecoBaseDrawer::" << fcn << " failed with message:\n" << e;
75  }
76 } // namespace
77 
78 namespace evd {
79 
80  //......................................................................
81  RecoBaseDrawer::RecoBaseDrawer()
82  {
83  art::ServiceHandle<geo::Geometry const> geo;
84  art::ServiceHandle<evd::RawDrawingOptions const> rawOptions;
85  art::ServiceHandle<evd::RecoDrawingOptions const> recoOptions;
86 
87  fWireMin.resize(0);
88  fWireMax.resize(0);
89  fTimeMin.resize(0);
90  fTimeMax.resize(0);
91  fRawCharge.resize(0);
92  fConvertedCharge.resize(0);
93 
94  // set the list of channels in this detector
95  for (size_t t = 0; t < geo->NTPC(); ++t) {
96  unsigned int nplanes = geo->Nplanes(t);
97  fWireMin.resize(nplanes, -1);
98  fWireMax.resize(nplanes, -1);
99  fTimeMin.resize(nplanes, -1);
100  fTimeMax.resize(nplanes, -1);
101  fRawCharge.resize(nplanes, 0);
102  fConvertedCharge.resize(nplanes, 0);
103  for (size_t p = 0; p < geo->Nplanes(t); ++p) {
104  fWireMin[p] = 0;
105  fWireMax[p] = geo->TPC(t).Plane(p).Nwires();
106  fTimeMin[p] = 0;
107  fTimeMax[p] = rawOptions->fTicks;
108  } // end loop over planes
109  } // end loop over TPCs
110 
111  fAllSpacePointDrawer =
112  art::make_tool<evdb_tool::ISpacePoints3D>(recoOptions->fAllSpacePointDrawerParams);
113  fSpacePointDrawer =
114  art::make_tool<evdb_tool::ISpacePoints3D>(recoOptions->fSpacePointDrawerParams);
115  }
116 
117  //......................................................................
118  RecoBaseDrawer::~RecoBaseDrawer() {}
119 
120  //......................................................................
121  void
122  RecoBaseDrawer::Wire2D(const art::Event& evt, evdb::View2D* view, unsigned int plane)
123  {
124  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
125  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
126  art::ServiceHandle<geo::Geometry const> geo;
127  art::ServiceHandle<evd::ColorDrawingOptions const> cst;
128 
129  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
130 
131  lariov::ChannelStatusProvider const& channelStatus =
132  art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();
133 
134  int ticksPerPoint = rawOpt->fTicksPerPoint;
135 
136  // to make det independent later:
137  double mint = 5000;
138  double maxt = 0;
139  double minw = 5000;
140  double maxw = 0;
141 
142  geo::PlaneID pid(rawOpt->fCryostat, rawOpt->fTPC, plane);
143 
144  for (size_t imod = 0; imod < recoOpt->fWireLabels.size(); ++imod) {
145  art::InputTag const which = recoOpt->fWireLabels[imod];
146 
147  art::PtrVector<recob::Wire> wires;
148  this->GetWires(evt, which, wires);
149 
150  if (wires.size() < 1) continue;
151 
152  for (size_t i = 0; i < wires.size(); ++i) {
153 
154  uint32_t channel = wires[i]->Channel();
155 
156  if (!rawOpt->fSeeBadChannels && channelStatus.IsBad(channel)) continue;
157 
158  std::vector<geo::WireID> wireids = geo->ChannelToWire(channel);
159 
160  geo::SigType_t sigType = geo->SignalType(channel);
161 
162  for (auto const& wid : wireids) {
163  if (wid.planeID() != pid) continue;
164 
165  double wire = 1. * wid.Wire;
166  double tick = 0;
167  // get the unpacked ROIs
168  std::vector<float> wirSig = wires[i]->Signal();
169  if (wirSig.size() == 0) continue;
170  // get an iterator over the adc values
171  std::vector<float>::const_iterator itr = wirSig.begin();
172  while (itr != wirSig.end()) {
173  int ticksUsed = 0;
174  double tdcsum = 0.;
175  double adcsum = 0.;
176  while (ticksUsed < ticksPerPoint && itr != wirSig.end()) {
177  tdcsum += tick;
178  adcsum += (1. * (*itr));
179  ++ticksUsed;
180  tick += 1.;
181  itr++; // this advance of the iterator is sufficient for the external loop too
182  }
183  double adc = adcsum / ticksPerPoint;
184  double tdc = tdcsum / ticksPerPoint;
185 
186  if (TMath::Abs(adc) < rawOpt->fMinSignal) continue;
187  if (tdc > rawOpt->fTicks) continue;
188 
189  int co = 0;
190  double sf = 1.;
191  double q0 = 1000.0;
192 
193  co = cst->CalQ(sigType).GetColor(adc);
194  if (rawOpt->fScaleDigitsByCharge) {
195  sf = sqrt(adc / q0);
196  if (sf > 1.0) sf = 1.0;
197  }
198 
199  if (wire < minw) minw = wire;
200  if (wire > maxw) maxw = wire;
201  if (tdc < mint) mint = tdc;
202  if (tdc > maxt) maxt = tdc;
203 
204  if (rawOpt->fAxisOrientation < 1) {
205  TBox& b1 = view->AddBox(wire - sf * 0.5,
206  tdc - sf * 0.5 * ticksPerPoint,
207  wire + sf * 0.5,
208  tdc + sf * 0.5 * ticksPerPoint);
209  b1.SetFillStyle(1001);
210  b1.SetFillColor(co);
211  b1.SetBit(kCannotPick);
212  }
213  else {
214  TBox& b1 = view->AddBox(tdc - sf * 0.5 * ticksPerPoint,
215  wire - sf * 0.5,
216  tdc + sf * 0.5 * ticksPerPoint,
217  wire + sf * 0.5);
218  b1.SetFillStyle(1001);
219  b1.SetFillColor(co);
220  b1.SetBit(kCannotPick);
221  }
222  } // end loop over samples
223  } //end loop over wire segments
224  } //end loop over wires
225  } // end loop over wire module labels
226 
227  fWireMin[plane] = minw;
228  fWireMax[plane] = maxw;
229  fTimeMin[plane] = mint;
230  fTimeMax[plane] = maxt;
231 
232  // Add a loop to draw dead wires in 2D display
233  double startTick(50.);
234  double endTick((rawOpt->fTicks - 50.) * ticksPerPoint);
235 
236  for (size_t wireNo = 0; wireNo < geo->Nwires(pid); wireNo++) {
237  raw::ChannelID_t channel =
238  geo->PlaneWireToChannel(geo::WireID(rawOpt->fCryostat, rawOpt->fTPC, plane, wireNo));
239 
240  if (!rawOpt->fSeeBadChannels && channelStatus.IsBad(channel)) {
241  double wire = 1. * wireNo;
242  TLine& line = view->AddLine(wire, startTick, wire, endTick);
243  line.SetLineColor(kGray);
244  line.SetLineWidth(1.0);
245  line.SetBit(kCannotPick);
246  }
247  }
248  }
249 
250  //......................................................................
251  ///
252  /// Render Hit objects on a 2D viewing canvas
253  ///
254  /// @param evt : Event handle to get data objects from
255  /// @param view : Pointer to view to draw on
256  /// @param plane : plane number of view
257  ///
258  int
259  RecoBaseDrawer::Hit2D(const art::Event& evt,
260  detinfo::DetectorPropertiesData const& detProp,
261  evdb::View2D* view,
262  unsigned int plane)
263  {
264  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
265  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
266  art::ServiceHandle<geo::Geometry const> geo;
267 
268  int nHitsDrawn(0);
269 
270  if (recoOpt->fDrawHits == 0) return nHitsDrawn;
271  if (rawOpt->fDrawRawDataOrCalibWires < 1) return nHitsDrawn;
272 
273  fRawCharge[plane] = 0;
274  fConvertedCharge[plane] = 0;
275 
276  for (size_t imod = 0; imod < recoOpt->fHitLabels.size(); ++imod) {
277  art::InputTag const which = recoOpt->fHitLabels[imod];
278 
279  std::vector<const recob::Hit*> hits;
280  this->GetHits(evt, which, hits, plane);
281 
282  // Display all hits on the two 2D views provided
283  for (auto itr : hits) {
284 
285  if (itr->WireID().TPC != rawOpt->fTPC || itr->WireID().Cryostat != rawOpt->fCryostat)
286  continue;
287 
288  // Try to get the "best" charge measurement, ie. the one last in
289  // the calibration chain
290  fRawCharge[itr->WireID().Plane] += itr->PeakAmplitude();
291  double dQdX = itr->PeakAmplitude() / geo->WirePitch() / detProp.ElectronsToADC();
292  fConvertedCharge[itr->WireID().Plane] += detProp.BirksCorrection(dQdX);
293  } // loop on hits
294 
295  nHitsDrawn = this->Hit2D(hits, kBlack, view, recoOpt->fDrawAllWireIDs);
296 
297  } // loop on imod folders
298 
299  return nHitsDrawn;
300  }
301 
302  //......................................................................
303  ///
304  /// Render Hit objects on a 2D viewing canvas
305  ///
306  /// @param hits : vector of hits for the veiw
307  /// @param color : color of associated cluster/prong
308  /// @param view : Pointer to view to draw on
309  ///
310  /// assumes the hits are all from the correct plane for the given view
311  int
312  RecoBaseDrawer::Hit2D(std::vector<const recob::Hit*> hits,
313  int color,
314  evdb::View2D* view,
315  bool allWireIDs,
316  bool drawConnectingLines,
317  int lineWidth)
318  {
319  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
320  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
321  art::ServiceHandle<geo::Geometry const> geo;
322 
323  unsigned int w = 0;
324  unsigned int wold = 0;
325  float timeold = 0.;
326 
327  if (color == -1) color = recoOpt->fSelectedHitColor;
328 
329  int nHitsDrawn(0);
330 
331  for (const auto& hit : hits) {
332  // Note that the WireID in the hit object is useless for those detectors where a channel can correspond to
333  // more than one plane/wire. So our plan is to recover the list of wire IDs from the channel number and
334  // loop over those (if there are any)
335  // However, we need to preserve the option for drawing hits only associated to the wireID it contains
336  std::vector<geo::WireID> wireIDs;
337 
338  if (allWireIDs)
339  wireIDs = geo->ChannelToWire(hit->Channel());
340  else
341  wireIDs.push_back(hit->WireID());
342 
343  // Loop to find match
344  for (const auto& wireID : wireIDs) {
345  if (wireID.TPC != rawOpt->fTPC || wireID.Cryostat != rawOpt->fCryostat) continue;
346 
347  if (std::isnan(hit->PeakTime()) || std::isnan(hit->Integral())) {
348  std::cout << "====>> Found hit with a NAN, channel: " << hit->Channel()
349  << ", start/end: " << hit->StartTick() << "/" << hit->EndTick()
350  << ", chisquare: " << hit->GoodnessOfFit() << std::endl;
351  }
352 
353  if (hit->PeakTime() > rawOpt->fTicks) continue;
354 
355  w = wireID.Wire;
356 
357  // Try to get the "best" charge measurement, ie. the one last in
358  // the calibration chain
359  float time = hit->PeakTime();
360  float rms = 0.5 * hit->RMS();
361 
362  if (rawOpt->fAxisOrientation < 1) {
363  TBox& b1 = view->AddBox(w - 0.5, time - rms, w + 0.5, time + rms);
364  if (drawConnectingLines && nHitsDrawn > 0) {
365  TLine& l = view->AddLine(w, time, wold, timeold);
366  l.SetLineColor(color);
367  l.SetBit(kCannotPick);
368  }
369  b1.SetFillStyle(0);
370  b1.SetBit(kCannotPick);
371  b1.SetLineColor(color);
372  b1.SetLineWidth(lineWidth);
373  }
374  else {
375  TBox& b1 = view->AddBox(time - rms, w - 0.5, time + rms, w + 0.5);
376  if (drawConnectingLines && nHitsDrawn > 0) {
377  TLine& l = view->AddLine(time, w, timeold, wold);
378  l.SetLineColor(color);
379  l.SetBit(kCannotPick);
380  }
381  b1.SetFillStyle(0);
382  b1.SetBit(kCannotPick);
383  b1.SetLineColor(color);
384  b1.SetLineWidth(lineWidth);
385  }
386  wold = w;
387  timeold = time;
388  nHitsDrawn++;
389  }
390  } // loop on hits
391 
392  return nHitsDrawn;
393  }
394 
395  //........................................................................
396  int
397  RecoBaseDrawer::Hit2D(std::vector<const recob::Hit*> hits, evdb::View2D* view, float cosmicscore)
398  {
399  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
400  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
401  art::ServiceHandle<geo::Geometry const> geo;
402 
403  unsigned int w(0);
404  unsigned int wold(0);
405  float timeold(0.);
406  int nHitsDrawn(0);
407 
408  for (const auto& hit : hits) {
409  // check that we are in the correct TPC
410  // the view should tell use we are in the correct plane
411  if (hit->WireID().TPC != rawOpt->fTPC || hit->WireID().Cryostat != rawOpt->fCryostat)
412  continue;
413 
414  w = hit->WireID().Wire;
415 
416  // Try to get the "best" charge measurement, ie. the one last in
417  // the calibration chain
418  float time = hit->PeakTime();
419 
420  if (rawOpt->fAxisOrientation < 1) {
421  if (nHitsDrawn > 0) {
422  TLine& l = view->AddLine(w, time + 100, wold, timeold + 100);
423  l.SetLineWidth(3);
424  l.SetLineColor(1);
425  if (cosmicscore > 0.5) l.SetLineColor(kMagenta);
426  l.SetBit(kCannotPick);
427  }
428  }
429  else {
430  if (nHitsDrawn > 0) {
431  TLine& l = view->AddLine(time + 20, w, timeold + 20, wold);
432  l.SetLineColor(1);
433  if (cosmicscore > 0.5) l.SetLineStyle(2);
434  l.SetBit(kCannotPick);
435  }
436  }
437 
438  wold = w;
439  timeold = time;
440  nHitsDrawn++;
441  } // loop on hits
442 
443  return nHitsDrawn;
444  }
445 
446  //........................................................................
447  int
448  RecoBaseDrawer::GetRegionOfInterest(int plane, int& minw, int& maxw, int& mint, int& maxt)
449  {
450  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
451  art::ServiceHandle<geo::Geometry const> geo;
452 
453  if ((unsigned int)plane > fWireMin.size()) {
454  mf::LogWarning("RecoBaseDrawer")
455  << " Requested plane " << plane << " is larger than those available ";
456  return -1;
457  }
458 
459  minw = fWireMin[plane];
460  maxw = fWireMax[plane];
461  mint = fTimeMin[plane];
462  maxt = fTimeMax[plane];
463 
464  //make values a bit larger, but make sure they don't go out of bounds
465  minw = (minw - 30 < 0) ? 0 : minw - 30;
466  mint = (mint - 10 < 0) ? 0 : mint - 10;
467 
468  int fTicks = rawOpt->fTicks;
469 
470  maxw = (maxw + 10 > (int)geo->Nwires(plane)) ? geo->Nwires(plane) : maxw + 10;
471  maxt = (maxt + 10 > fTicks) ? fTicks : maxt + 10;
472 
473  return 0;
474  }
475 
476  //......................................................................
477  void
478  RecoBaseDrawer::GetChargeSum(int plane, double& charge, double& convcharge)
479  {
480  charge = fRawCharge[plane];
481  convcharge = fConvertedCharge[plane];
482 
483  return;
484  }
485 
486  //......................................................................
487  void
488  RecoBaseDrawer::EndPoint2D(const art::Event& evt, evdb::View2D* view, unsigned int plane)
489  {
490  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
491  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
492  art::ServiceHandle<geo::Geometry const> geo;
493 
494  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
495  if (recoOpt->fDraw2DEndPoints == 0) return;
496 
497  geo::View_t gview = geo->TPC(rawOpt->fTPC).Plane(plane).View();
498 
499  for (size_t imod = 0; imod < recoOpt->fEndPoint2DLabels.size(); ++imod) {
500  art::InputTag const which = recoOpt->fEndPoint2DLabels[imod];
501 
502  art::PtrVector<recob::EndPoint2D> ep2d;
503  this->GetEndPoint2D(evt, which, ep2d);
504 
505  for (size_t iep = 0; iep < ep2d.size(); ++iep) {
506  // only worry about end points with the correct view
507  if (ep2d[iep]->View() != gview) continue;
508 
509  ///\todo - have to verify that we are in the right TPC, but to do that we
510  // need to be sure that all EndPoint2D objects have filled the required information
511 
512  // draw cluster with unique marker
513  // Place this cluster's unique marker at the hit's location
514  int color = evd::kColor[ep2d[iep]->ID() % evd::kNCOLS];
515 
516  double x = ep2d[iep]->WireID().Wire;
517  double y = ep2d[iep]->DriftTime();
518 
519  if (rawOpt->fAxisOrientation > 0) {
520  x = ep2d[iep]->DriftTime();
521  y = ep2d[iep]->WireID().Wire;
522  }
523 
524  TMarker& strt = view->AddMarker(x, y, color, 30, 2.0);
525  strt.SetMarkerColor(color);
526  // BB: draw the ID
527  if (recoOpt->fDraw2DEndPoints > 1) {
528  std::string s = "2V" + std::to_string(ep2d[iep]->ID());
529  char const* txt = s.c_str();
530  TText& vtxID = view->AddText(x, y + 20, txt);
531  vtxID.SetTextColor(color);
532  vtxID.SetTextSize(0.05);
533  }
534 
535  } // loop on iep end points
536  } // loop on imod folders
537 
538  return;
539  }
540 
541  //......................................................................
542  void
543  RecoBaseDrawer::OpFlash2D(const art::Event& evt,
544  detinfo::DetectorClocksData const& clockData,
545  detinfo::DetectorPropertiesData const& detProp,
546  evdb::View2D* view,
547  unsigned int plane)
548  {
549  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
550  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
551 
552  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
553  if (recoOpt->fDrawOpFlashes == 0) return;
554 
555  art::ServiceHandle<geo::Geometry const> geo;
556  geo::PlaneID pid(rawOpt->fCryostat, rawOpt->fTPC, plane);
557 
558  for (size_t imod = 0; imod < recoOpt->fOpFlashLabels.size(); ++imod) {
559  const art::InputTag which = recoOpt->fOpFlashLabels[imod];
560 
561  art::PtrVector<recob::OpFlash> opflashes;
562  this->GetOpFlashes(evt, which, opflashes);
563 
564  if (opflashes.size() < 1) continue;
565 
566  int NFlashes = opflashes.size();
567  //double TopCoord = 1000;
568 
569  MF_LOG_VERBATIM("RecoBaseDrawer") << "Total " << NFlashes << " flashes.";
570 
571  // project each seed into this view
572  for (size_t iof = 0; iof < opflashes.size(); ++iof) {
573  if (opflashes[iof]->TotalPE() < recoOpt->fFlashMinPE) continue;
574  if (opflashes[iof]->Time() < recoOpt->fFlashTMin) continue;
575  if (opflashes[iof]->Time() > recoOpt->fFlashTMax) continue;
576  int Color = evd::kColor[(iof) % evd::kNCOLS];
577  MF_LOG_VERBATIM("RecoBaseDrawer")
578  << "Flash t: " << opflashes[iof]->Time() << "\t y,z : " << opflashes[iof]->YCenter()
579  << ", " << opflashes[iof]->ZCenter() << " \t PE :" << opflashes[iof]->TotalPE();
580 
581  float flashtick =
582  opflashes[iof]->Time() / sampling_rate(clockData) * 1e3 + detProp.GetXTicksOffset(pid);
583  float wire0 = FLT_MAX;
584  float wire1 = FLT_MIN;
585 
586  //Find the 4 corners and convert them to wire numbers
587  std::vector<TVector3> points;
588  points.push_back(TVector3(0,
589  opflashes[iof]->YCenter() - opflashes[iof]->YWidth(),
590  opflashes[iof]->ZCenter() - opflashes[iof]->ZWidth()));
591  points.push_back(TVector3(0,
592  opflashes[iof]->YCenter() - opflashes[iof]->YWidth(),
593  opflashes[iof]->ZCenter() + opflashes[iof]->ZWidth()));
594  points.push_back(TVector3(0,
595  opflashes[iof]->YCenter() + opflashes[iof]->YWidth(),
596  opflashes[iof]->ZCenter() - opflashes[iof]->ZWidth()));
597  points.push_back(TVector3(0,
598  opflashes[iof]->YCenter() + opflashes[iof]->YWidth(),
599  opflashes[iof]->ZCenter() + opflashes[iof]->ZWidth()));
600 
601  for (size_t i = 0; i < points.size(); ++i) {
602  geo::WireID wireID;
603  try {
604  wireID = geo->NearestWireID(points[i], pid);
605  }
606  catch (geo::InvalidWireError const& e) {
607  wireID = e.suggestedWireID(); // pick the closest valid wire
608  }
609  if (wireID.Wire < wire0) wire0 = wireID.Wire;
610  if (wireID.Wire > wire1) wire1 = wireID.Wire;
611  }
612  if (rawOpt->fAxisOrientation > 0) {
613  TLine& line = view->AddLine(flashtick, wire0, flashtick, wire1);
614  line.SetLineWidth(2);
615  line.SetLineStyle(2);
616  line.SetLineColor(Color);
617  }
618  else {
619  TLine& line = view->AddLine(wire0, flashtick, wire1, flashtick);
620  line.SetLineWidth(2);
621  line.SetLineStyle(2);
622  line.SetLineColor(Color);
623  }
624  } // loop on opflashes
625  } // loop on imod folders
626 
627  return;
628  }
629 
630  //......................................................................
631  void
632  RecoBaseDrawer::Seed2D(const art::Event& evt,
633  detinfo::DetectorPropertiesData const& detProp,
634  evdb::View2D* view,
635  unsigned int plane)
636  {
637  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
638  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
639  art::ServiceHandle<geo::Geometry const> geo;
640 
641  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
642  if (recoOpt->fDrawSeeds == 0) return;
643 
644  for (size_t imod = 0; imod < recoOpt->fSeedLabels.size(); ++imod) {
645  art::InputTag const which = recoOpt->fSeedLabels[imod];
646 
647  art::PtrVector<recob::Seed> seeds;
648  this->GetSeeds(evt, which, seeds);
649 
650  if (seeds.size() < 1) continue;
651 
652  // project each seed into this view
653  for (size_t isd = 0; isd < seeds.size(); ++isd) {
654  double SeedPoint[3];
655  double SeedDir[3];
656  double SeedPointErr[3];
657  double SeedDirErr[3];
658  double SeedEnd1[3];
659  double SeedEnd2[3];
660 
661  seeds[isd]->GetPoint(SeedPoint, SeedPointErr);
662  seeds[isd]->GetDirection(SeedDir, SeedDirErr);
663 
664  SeedEnd1[0] = SeedPoint[0] + SeedDir[0];
665  SeedEnd1[1] = SeedPoint[1] + SeedDir[1];
666  SeedEnd1[2] = SeedPoint[2] + SeedDir[2];
667 
668  SeedEnd2[0] = SeedPoint[0] - SeedDir[0];
669  SeedEnd2[1] = SeedPoint[1] - SeedDir[1];
670  SeedEnd2[2] = SeedPoint[2] - SeedDir[2];
671 
672  // Draw seed on evd
673  // int color = kColor[seeds[isd]->ID()%kNCOLS];
674  int color = evd::kColor[0];
675  unsigned int wirepoint = 0;
676  unsigned int wireend1 = 0;
677  unsigned int wireend2 = 0;
678  try {
679  wirepoint = geo->NearestWire(SeedPoint, plane, rawOpt->fTPC, rawOpt->fCryostat);
680  }
681  catch (cet::exception& e) {
682  wirepoint = atoi(e.explain_self().substr(e.explain_self().find("#") + 1, 5).c_str());
683  }
684  try {
685  wireend1 = geo->NearestWire(SeedEnd1, plane, rawOpt->fTPC, rawOpt->fCryostat);
686  }
687  catch (cet::exception& e) {
688  wireend1 = atoi(e.explain_self().substr(e.explain_self().find("#") + 1, 5).c_str());
689  }
690  try {
691  wireend2 = geo->NearestWire(SeedEnd2, plane, rawOpt->fTPC, rawOpt->fCryostat);
692  }
693  catch (cet::exception& e) {
694  wireend2 = atoi(e.explain_self().substr(e.explain_self().find("#") + 1, 5).c_str());
695  }
696 
697  double x = wirepoint;
698  double y = detProp.ConvertXToTicks(SeedPoint[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
699  double x1 = wireend1;
700  double y1 = detProp.ConvertXToTicks(SeedEnd1[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
701  double x2 = wireend2;
702  double y2 = detProp.ConvertXToTicks(SeedEnd2[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
703 
704  if (rawOpt->fAxisOrientation > 0) {
705  x = detProp.ConvertXToTicks(SeedPoint[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
706  y = wirepoint;
707  x1 = detProp.ConvertXToTicks(SeedEnd1[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
708  y1 = wireend1;
709  x2 = detProp.ConvertXToTicks(SeedEnd2[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
710  y2 = wireend2;
711  }
712 
713  TMarker& strt = view->AddMarker(x, y, color, 4, 1.5);
714  TLine& line = view->AddLine(x1, y1, x2, y2);
715  strt.SetMarkerColor(color);
716  line.SetLineColor(color);
717  line.SetLineWidth(2.0);
718  } // loop on seeds
719  } // loop on imod folders
720 
721  return;
722  }
723 
724  //......................................................................
725  void
726  RecoBaseDrawer::Slice2D(const art::Event& evt,
727  detinfo::DetectorPropertiesData const& detProp,
728  evdb::View2D* view,
729  unsigned int plane)
730  {
731  // Color code hits associated with Slices
732  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
733  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
734  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
735  if (recoOpt->fDrawSlices == 0) return;
736 
737  art::ServiceHandle<geo::Geometry const> geo;
738 
739  static bool first = true;
740  if (first) {
741  std::cout
742  << "******** DrawSlices: 0 = none, 1 = color coded, 2 = color coded + ID at slice center\n";
743  std::cout << " 3 = open circle at slice center with size proportional to the AspectRatio. "
744  "Closed circles";
745  std::cout << " at the slice ends with connecting dotted lines\n";
746  first = false;
747  }
748  unsigned int c = rawOpt->fCryostat;
749  unsigned int t = rawOpt->fTPC;
750  geo::PlaneID planeID(c, t, plane);
751 
752  for (size_t imod = 0; imod < recoOpt->fSliceLabels.size(); ++imod) {
753  art::InputTag const which = recoOpt->fSliceLabels[imod];
754  art::PtrVector<recob::Slice> slices;
755  this->GetSlices(evt, which, slices);
756  if (slices.size() < 1) continue;
757  art::FindMany<recob::Hit> fmh(slices, evt, which);
758  for (size_t isl = 0; isl < slices.size(); ++isl) {
759  int slcID(std::abs(slices[isl]->ID()));
760  int color(evd::kColor[slcID % evd::kNCOLS]);
761  if (recoOpt->fDrawSlices < 3) {
762  // draw color-coded hits
763  std::vector<const recob::Hit*> hits = fmh.at(isl);
764  std::vector<const recob::Hit*> hits_on_plane;
765  for (auto hit : hits) {
766  if (hit->WireID().Plane == plane) { hits_on_plane.push_back(hit); }
767  }
768  if (this->Hit2D(hits_on_plane, color, view, false, false) < 1) continue;
769  if (recoOpt->fDrawSlices == 2) {
770  geo::Point_t slicePos(
771  slices[isl]->Center().X(), slices[isl]->Center().Y(), slices[isl]->Center().Z());
772  double tick = detProp.ConvertXToTicks(slices[isl]->Center().X(), planeID);
773  double wire = geo->WireCoordinate(slicePos, planeID);
774  std::string s = std::to_string(slcID);
775  char const* txt = s.c_str();
776  TText& slcID = view->AddText(wire, tick, txt);
777  slcID.SetTextSize(0.05);
778  slcID.SetTextColor(color);
779  } // draw ID
780  }
781  else {
782  // draw the center, end points and direction vector
783  geo::Point_t slicePos(
784  slices[isl]->Center().X(), slices[isl]->Center().Y(), slices[isl]->Center().Z());
785  double tick = detProp.ConvertXToTicks(slices[isl]->Center().X(), planeID);
786  double wire = geo->WireCoordinate(slicePos, planeID);
787  float markerSize = 1;
788  if (slices[isl]->AspectRatio() > 0) {
789  markerSize = 1 / slices[isl]->AspectRatio();
790  if (markerSize > 3) markerSize = 3;
791  }
792  TMarker& ctr = view->AddMarker(wire, tick, color, 24, markerSize);
793  ctr.SetMarkerColor(color);
794  // npts, color, width, style
795  TPolyLine& pline = view->AddPolyLine(2, color, 2, 3);
796  geo::Point_t slicePos0(
797  slices[isl]->End0Pos().X(), slices[isl]->End0Pos().Y(), slices[isl]->End0Pos().Z());
798  tick = detProp.ConvertXToTicks(slices[isl]->End0Pos().X(), planeID);
799  wire = geo->WireCoordinate(slicePos0, planeID);
800  TMarker& end0 = view->AddMarker(wire, tick, color, 20, 1.0);
801  end0.SetMarkerColor(color);
802  pline.SetPoint(0, wire, tick);
803  geo::Point_t slicePos1(
804  slices[isl]->End1Pos().X(), slices[isl]->End1Pos().Y(), slices[isl]->End1Pos().Z());
805  tick = detProp.ConvertXToTicks(slices[isl]->End1Pos().X(), plane, t, c);
806  wire = geo->WireCoordinate(slicePos1, planeID);
807  TMarker& end1 = view->AddMarker(wire, tick, color, 20, 1.0);
808  end1.SetMarkerColor(color);
809  pline.SetPoint(1, wire, tick);
810  }
811  } // isl
812 
813  } // imod
814 
815  } // Slice2D
816  //......................................................................
817  void
818  RecoBaseDrawer::Cluster2D(const art::Event& evt,
819  detinfo::DetectorClocksData const& clockData,
820  detinfo::DetectorPropertiesData const& detProp,
821  evdb::View2D* view,
822  unsigned int plane)
823  {
824  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
825  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
826  art::ServiceHandle<geo::Geometry const> geo;
827 
828  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
829  if (recoOpt->fDrawClusters == 0) return;
830 
831  geo::View_t gview = geo->TPC(rawOpt->fTPC).Plane(plane).View();
832 
833  // if user sets "DrawClusters" to 2, draw the clusters differently:
834  // bool drawAsMarkers = (recoOpt->fDrawClusters == 1 ||
835  // recoOpt->fDrawClusters == 3);
836  bool drawAsMarkers = recoOpt->fDrawClusters != 2;
837 
838  // draw connecting lines between cluster hits?
839  bool drawConnectingLines = (recoOpt->fDrawClusters >= 3);
840 
841  static bool first = true;
842  if (first) {
843  std::cout << "******** DrawClusters: 0 = none, 1 = cluster hits, 2 = unique marker, 3 = "
844  "cluster hits with connecting lines.\n";
845  std::cout << " 4 = with T<cluster or trajectory ID> P<PFParticle ID> color-matched. "
846  "Unmatched cluster IDs shown in black.\n";
847  std::cout << " Color scheme: By cluster ID in each plane or by PFParticle ID (Self) if a "
848  "PFParticle - Cluster association exists.\n";
849  first = false;
850  }
851 
852  for (size_t imod = 0; imod < recoOpt->fClusterLabels.size(); ++imod) {
853  art::InputTag const which = recoOpt->fClusterLabels[imod];
854 
855  art::PtrVector<recob::Cluster> clust;
856  this->GetClusters(evt, which, clust);
857 
858  if (clust.size() < 1) continue;
859 
860  // We want to draw the hits that are associated to "free" space points (non clustered)
861  // This is done here, before drawing the hits on clusters so they will be "under" the cluster
862  // hits (since spacepoints could be made from a used 2D hit but then not used themselves)
863  // Get the space points created by the PFParticle producer
864  std::vector<art::Ptr<recob::SpacePoint>> spacePointVec;
865  this->GetSpacePoints(evt, which, spacePointVec);
866 
867  // No space points no continue
868  if (spacePointVec.size() > 0) {
869  // Add the relations to recover associations cluster hits
870  art::FindManyP<recob::Hit> spHitAssnVec(spacePointVec, evt, which);
871 
872  if (spHitAssnVec.isValid()) {
873  // Create a local hit vector...
874  std::vector<const recob::Hit*> freeHitVec;
875 
876  // loop through space points looking for those that are free
877  for (const auto& spacePointPtr : spacePointVec) {
878  if (spacePointPtr->Chisq() < -99.) {
879  // Recover associated hits
880  const std::vector<art::Ptr<recob::Hit>>& hitVec =
881  spHitAssnVec.at(spacePointPtr.key());
882 
883  for (const auto& hitPtr : hitVec) {
884  if (hitPtr.get()->WireID().Plane != plane) continue;
885 
886  freeHitVec.push_back(hitPtr.get());
887  }
888  }
889  }
890 
891  // Draw the free hits in gray
892  this->Hit2D(freeHitVec, kGray, view, false, false, false);
893  }
894  }
895 
896  // Ok, now proceed with our normal processing of hits on clusters
897  art::FindMany<recob::Hit> fmh(clust, evt, which);
898  art::FindManyP<recob::PFParticle> fmc(clust, evt, which);
899 
900  for (size_t ic = 0; ic < clust.size(); ++ic) {
901  // only worry about clusters with the correct view
902  // if(clust[ic]->View() != gview) continue;
903  if (clust[ic]->Plane().Plane != plane) continue;
904 
905  // see if we can set the color index in a sensible fashion
906  int clusterIdx(std::abs(clust[ic]->ID()));
907  int colorIdx(clusterIdx % evd::kNCOLS);
908  bool pfpAssociation = false;
909  int pfpIndex = INT_MAX;
910  float cosmicscore = FLT_MIN;
911 
912  if (fmc.isValid()) {
913  std::vector<art::Ptr<recob::PFParticle>> pfplist = fmc.at(ic);
914  // Use the first one
915  if (!pfplist.empty()) {
916  clusterIdx = pfplist[0]->Self();
917  colorIdx = clusterIdx % evd::kNCOLS;
918  pfpAssociation = true;
919  pfpIndex = pfplist[0]->Self();
920  //Get cosmic score
921  if (recoOpt->fDrawCosmicTags) {
922  art::FindManyP<anab::CosmicTag> fmct(pfplist, evt, which);
923  if (fmct.isValid()) {
924  std::vector<art::Ptr<anab::CosmicTag>> ctlist = fmct.at(0);
925  if (!ctlist.empty()) {
926  //std::cout<<"cosmic tag "<<ctlist[0]->CosmicScore()<<std::endl;
927  cosmicscore = ctlist[0]->CosmicScore();
928  }
929  }
930  }
931  } // pfplist is not empty
932  }
933 
934  std::vector<const recob::Hit*> hits = fmh.at(ic);
935 
936  if (drawAsMarkers) {
937  // draw cluster with unique marker
938  // Place this cluster's unique marker at the hit's location
939  int color = evd::kColor[colorIdx];
940 
941  // If there are no hits in this cryostat/TPC then we skip the rest
942  // That no hits were drawn is the sign for this
943  if (this->Hit2D(hits, color, view, false, drawConnectingLines) < 1) continue;
944 
945  if (recoOpt->fDrawCosmicTags && cosmicscore != FLT_MIN)
946  this->Hit2D(hits, view, cosmicscore);
947 
948  if (recoOpt->fDrawClusters > 3) {
949  // BB: draw the cluster ID
950  //std::string s = std::to_string(clusterIdx);
951  // TY: change to draw cluster id instead of index
952  // std::string s = std::to_string(clusterIdx) + "," + std::to_string(clust[ic]->ID());
953  // BB: Put a T in front to denote a trajectory ID
954  std::string s = "T" + std::to_string(clust[ic]->ID());
955  // append the PFP index + 1 (sort of the ID)
956  if (pfpIndex != INT_MAX) s = s + " P" + std::to_string(pfpIndex + 1);
957  char const* txt = s.c_str();
958  double wire = 0.5 * (clust[ic]->StartWire() + clust[ic]->EndWire());
959  double tick = 20 + 0.5 * (clust[ic]->StartTick() + clust[ic]->EndTick());
960  TText& clID = view->AddText(wire, tick, txt);
961  clID.SetTextSize(0.05);
962  if (pfpAssociation) { clID.SetTextColor(color); }
963  else {
964  clID.SetTextColor(kBlack);
965  }
966  } // recoOpt->fDrawClusters > 3
967  }
968  else {
969 
970  // default "outline" method:
971  std::vector<double> tpts, wpts;
972 
973  this->GetClusterOutlines(hits, tpts, wpts, plane);
974 
975  int lcolor = 9; // line color
976  int fcolor = 9; // fill color
977  int width = 2; // line width
978  int style = 1; // 1=solid line style
979  if (view != 0) {
980  TPolyLine& p1 = view->AddPolyLine(wpts.size(), lcolor, width, style);
981  TPolyLine& p2 = view->AddPolyLine(wpts.size(), lcolor, width, style);
982  p1.SetOption("f");
983  p1.SetFillStyle(3003);
984  p1.SetFillColor(fcolor);
985  for (size_t i = 0; i < wpts.size(); ++i) {
986  if (rawOpt->fAxisOrientation < 1) {
987  p1.SetPoint(i, wpts[i], tpts[i]);
988  p2.SetPoint(i, wpts[i], tpts[i]);
989  }
990  else {
991  p1.SetPoint(i, tpts[i], wpts[i]);
992  p2.SetPoint(i, tpts[i], wpts[i]);
993  }
994  } // loop on i points in ZX view
995  } // if we have a cluster in the ZX view
996  } // end if outline mode
997 
998  // draw the direction cosine of the cluster as well as it's starting point
999  // (average of the start and end angle -- by default they are the same value)
1000  // thetawire is the angle measured CW from +z axis to wire
1001  //double thetawire = geo->TPC(t).Plane(plane).Wire(0).ThetaZ();
1002  double wirePitch = geo->WirePitch(gview);
1003  double driftvelocity = detProp.DriftVelocity(); // cm/us
1004  double timetick = sampling_rate(clockData) * 1e-3; // time sample in us
1005  // rotate coord system CCW around x-axis by pi-thetawire
1006  // new yprime direction is perpendicular to the wire direction
1007  // in the same plane as the wires and in the direction of
1008  // increasing wire number
1009  //use yprime-component of dir cos in rotated coord sys to get
1010  // dTdW (number of time ticks per unit of wire pitch)
1011  //double rotang = 3.1416-thetawire;
1012  this->Draw2DSlopeEndPoints(
1013  clust[ic]->StartWire(),
1014  clust[ic]->StartTick(),
1015  clust[ic]->EndWire(),
1016  clust[ic]->EndTick(),
1017  std::tan((clust[ic]->StartAngle() + clust[ic]->EndAngle()) / 2.) * wirePitch /
1018  driftvelocity / timetick,
1019  evd::kColor[colorIdx],
1020  view);
1021 
1022  } // loop on ic clusters
1023  } // loop on imod folders
1024 
1025  return;
1026  }
1027 
1028  //......................................................................
1029  void
1030  RecoBaseDrawer::Draw2DSlopeEndPoints(double xStart,
1031  double yStart,
1032  double xEnd,
1033  double yEnd,
1034  double slope,
1035  int color,
1036  evdb::View2D* view)
1037  {
1038  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1039  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1040 
1041  if (recoOpt->fDraw2DSlopeEndPoints < 1) return;
1042 
1043  double x1 = xStart;
1044  double y1 = yStart;
1045  double x2 = xEnd;
1046  double slope1 = slope;
1047 
1048  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1049  if (rawOpt->fAxisOrientation > 0) {
1050  x1 = yStart;
1051  y1 = xStart;
1052  x2 = yEnd;
1053  if (std::abs(slope) > 0.)
1054  slope1 = 1. / slope;
1055  else
1056  slope1 = 1.e6;
1057  }
1058 
1059  double deltaX = 0.5 * (x2 - x1);
1060  double xm = x1 + deltaX;
1061  double ym = y1 + deltaX * slope;
1062 
1063  TMarker& strt = view->AddMarker(xm, ym, color, kFullCircle, 1.0);
1064  strt.SetMarkerColor(color); // stupid line to shut up compiler warning
1065 
1066  // double stublen = 50.0 ;
1067  double stublen = 2. * deltaX;
1068  TLine& l = view->AddLine(x1, y1, x1 + stublen, y1 + slope1 * stublen);
1069  l.SetLineColor(color);
1070  l.SetLineWidth(1); //2);
1071 
1072  return;
1073  }
1074 
1075  //......................................................................
1076  void
1077  RecoBaseDrawer::Draw2DSlopeEndPoints(double x,
1078  double y,
1079  double slope,
1080  int color,
1081  evdb::View2D* view)
1082  {
1083  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1084  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1085 
1086  if (recoOpt->fDraw2DSlopeEndPoints < 1) return;
1087 
1088  double x1 = x;
1089  double y1 = y;
1090  double slope1 = slope;
1091 
1092  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1093  if (rawOpt->fAxisOrientation > 0) {
1094  x1 = y;
1095  y1 = x;
1096  if (std::abs(slope) > 0.)
1097  slope1 = 1. / slope;
1098  else
1099  slope1 = 1.e6;
1100  }
1101 
1102  TMarker& strt = view->AddMarker(x1, y1, color, kFullStar, 2.0);
1103  strt.SetMarkerColor(color); // stupid line to shut up compiler warning
1104 
1105  // double stublen = 50.0 ;
1106  double stublen = 300.0;
1107  TLine& l = view->AddLine(x1, y1, x1 + stublen, y1 + slope1 * stublen);
1108  l.SetLineColor(color);
1109  l.SetLineWidth(2);
1110  l.SetLineStyle(2);
1111 
1112  return;
1113  }
1114 
1115  //......................................................................
1116  void
1117  RecoBaseDrawer::Draw2DSlopeEndPoints(double x,
1118  double y,
1119  double cosx,
1120  double cosy,
1121  int color,
1122  evdb::View2D* view)
1123  {
1124  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1125  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1126 
1127  if (recoOpt->fDraw2DSlopeEndPoints < 1) return;
1128 
1129  double x1 = x;
1130  double y1 = y;
1131  double cosx1 = cosx;
1132  double cosy1 = cosy;
1133 
1134  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1135  if (rawOpt->fAxisOrientation > 0) {
1136  x1 = y;
1137  y1 = x;
1138  cosx1 = cosy;
1139  cosy1 = cosx;
1140  }
1141 
1142  TMarker& strt = view->AddMarker(x1, y1, color, kFullStar, 2.0);
1143  strt.SetMarkerColor(color); // stupid line to shut up compiler warning
1144 
1145  // double stublen = 50.0 ;
1146  double stublen = 300.0;
1147  TLine& l = view->AddLine(x1, y1, x1 + stublen * cosx1, y1 + stublen * cosy1);
1148  l.SetLineColor(color);
1149  l.SetLineWidth(2);
1150  l.SetLineStyle(2);
1151 
1152  return;
1153  }
1154 
1155  //......................................................................
1156  ///
1157  /// Make a set of points which outline a cluster
1158  ///
1159  /// @param c : Reco base cluster to outline
1160  /// @param wpts : wire values of the outlines
1161  /// @param tpts : tdc values of the outlines
1162  /// @param plane : plane number
1163  ///
1164  void
1165  RecoBaseDrawer::GetClusterOutlines(std::vector<const recob::Hit*>& hits,
1166  std::vector<double>& wpts,
1167  std::vector<double>& tpts,
1168  unsigned int plane)
1169  {
1170  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1171 
1172  // Map wire numbers to highest and lowest in the plane
1173  std::map<unsigned int, double> wlo, whi;
1174  // On first pass, initialize
1175  for (size_t j = 0; j < hits.size(); ++j) {
1176  // check that we are on the correct plane and TPC
1177  if (hits[j]->WireID().Plane != plane || hits[j]->WireID().TPC != rawOpt->fTPC ||
1178  hits[j]->WireID().Cryostat != rawOpt->fCryostat)
1179  continue;
1180 
1181  wlo[hits[j]->WireID().Wire] = hits[j]->PeakTime();
1182  whi[hits[j]->WireID().Wire] = hits[j]->PeakTime();
1183  }
1184 
1185  double t = 0.;
1186 
1187  // Finalize on second pass
1188  for (size_t j = 0; j < hits.size(); ++j) {
1189  t = hits[j]->PeakTime();
1190 
1191  if (t < wlo[hits[j]->WireID().Wire]) wlo[hits[j]->WireID().Wire] = t;
1192  if (t > whi[hits[j]->WireID().Wire]) whi[hits[j]->WireID().Wire] = t;
1193  }
1194 
1195  // Loop over wires and low times to make lines along bottom
1196  // edge. Work from upstream edge to downstream edge
1197  std::map<unsigned int, double>::iterator itr(wlo.begin());
1198  std::map<unsigned int, double>::iterator itrEnd(wlo.end());
1199  for (; itr != itrEnd; ++itr) {
1200  unsigned int w = itr->first;
1201  t = itr->second;
1202 
1203  wpts.push_back(1. * w - 0.1);
1204  tpts.push_back(t - 0.1);
1205  wpts.push_back(1. * w + 0.1);
1206  tpts.push_back(t - 0.1);
1207  }
1208 
1209  // Loop over planes and high cells to make lines along top
1210  // edge. Work from downstream edge toward upstream edge
1211  std::map<unsigned int, double>::reverse_iterator ritr(whi.rbegin());
1212  std::map<unsigned int, double>::reverse_iterator ritrEnd(whi.rend());
1213  for (; ritr != ritrEnd; ++ritr) {
1214  unsigned int w = ritr->first;
1215  t = ritr->second;
1216 
1217  wpts.push_back(1. * w + 0.1);
1218  tpts.push_back(t + 0.1);
1219  wpts.push_back(1. * w - 0.1);
1220  tpts.push_back(t + 0.1);
1221  }
1222 
1223  // Add link to starting point to close the box
1224  wpts.push_back(wpts[0]);
1225  tpts.push_back(tpts[0]);
1226 
1227  return;
1228  }
1229 
1230  //......................................................................
1231  void
1232  RecoBaseDrawer::DrawProng2D(detinfo::DetectorPropertiesData const& detProp,
1233  std::vector<const recob::Hit*>& hits,
1234  evdb::View2D* view,
1235  unsigned int plane,
1236  TVector3 const& startPos,
1237  TVector3 const& startDir,
1238  int id,
1239  float cscore)
1240  {
1241  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1242  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1243  art::ServiceHandle<geo::Geometry const> geo;
1244 
1245  unsigned int c = rawOpt->fCryostat;
1246  unsigned int t = rawOpt->fTPC;
1247  geo::PlaneID planeID(c, t, plane);
1248  geo::Point_t localPos(startPos.X(), startPos.Y(), startPos.Z());
1249 
1250  int color(evd::kColor2[id % evd::kNCOLS]);
1251  int lineWidth(1);
1252 
1253  if (cscore > 0.1 && recoOpt->fDrawCosmicTags) {
1254  color = kRed;
1255  if (cscore < 0.6) color = kMagenta;
1256  lineWidth = 3;
1257  }
1258  else if (cscore < -10000) { //shower hits
1259  lineWidth = 3;
1260  }
1261 
1262  // first draw the hits
1263  if (cscore < -1000) { //shower
1264  this->Hit2D(hits, color, view, false, false, lineWidth);
1265  if (recoOpt->fDrawShowers >= 1) {
1266  //draw the shower ID at the beginning of shower
1267  std::string s = std::to_string(id);
1268  char const* txt = s.c_str();
1269  double tick = 30 + detProp.ConvertXToTicks(startPos.X(), planeID);
1270  double wire = geo->WireCoordinate(localPos, planeID);
1271  TText& shwID = view->AddText(wire, tick, txt);
1272  shwID.SetTextColor(evd::kColor2[id % evd::kNCOLS]);
1273  shwID.SetTextSize(0.1);
1274  }
1275  }
1276  else
1277  this->Hit2D(hits, color, view, false, false, lineWidth);
1278 
1279  double tick0 = detProp.ConvertXToTicks(startPos.X(), planeID);
1280  double wire0 = geo->WireCoordinate(localPos, planeID);
1281 
1282  localPos = geo::Point_t(startPos + startDir); // Huh? what is this?
1283 
1284  double tick1 = detProp.ConvertXToTicks((startPos + startDir).X(), planeID);
1285  double wire1 = geo->WireCoordinate(localPos, planeID);
1286  double cost = 0;
1287  double cosw = 0;
1288  double ds = sqrt(pow(tick0 - tick1, 2) + pow(wire0 - wire1, 2));
1289 
1290  if (ds) {
1291  cost = (tick1 - tick0) / ds;
1292  cosw = (wire1 - wire0) / ds;
1293  }
1294 
1295  this->Draw2DSlopeEndPoints(wire0, tick0, cosw, cost, evd::kColor[id % evd::kNCOLS], view);
1296 
1297  return;
1298  }
1299 
1300  //......................................................................
1301  void
1302  RecoBaseDrawer::DrawTrack2D(detinfo::DetectorClocksData const& clockData,
1303  detinfo::DetectorPropertiesData const& detProp,
1304  std::vector<const recob::Hit*>& hits,
1305  evdb::View2D* view,
1306  unsigned int plane,
1307  const recob::Track* track,
1308  int color,
1309  int lineWidth)
1310  {
1311  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1312  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1313  art::ServiceHandle<geo::Geometry const> geo;
1314  unsigned int c = rawOpt->fCryostat;
1315  unsigned int t = rawOpt->fTPC;
1316 
1317  // first draw the hits
1318  this->Hit2D(hits, color, view, false, true, lineWidth);
1319 
1320  const auto& startPos = track->Vertex();
1321  const auto& startDir = track->VertexDirection();
1322 
1323  // prepare to draw prongs
1324  double local[3] = {0.};
1325  double world[3] = {0.};
1326  geo->Cryostat(c).TPC(t).Plane(plane).LocalToWorld(local, world);
1327  world[1] = startPos.Y();
1328  world[2] = startPos.Z();
1329 
1330  // convert the starting position and direction from 3D to 2D coordinates
1331  double tick = detProp.ConvertXToTicks(startPos.X(), plane, t, c);
1332  double wire = 0.;
1333  try {
1334  wire = 1. * geo->NearestWire(world, plane, t, c);
1335  }
1336  catch (cet::exception& e) {
1337  wire = 1. * atoi(e.explain_self().substr(e.explain_self().find("#") + 1, 5).c_str());
1338  }
1339 
1340  // thetawire is the angle measured CW from +z axis to wire
1341  double thetawire = geo->TPC(t).Plane(plane).Wire(0).ThetaZ();
1342  double wirePitch = geo->WirePitch(hits[0]->View());
1343  double driftvelocity = detProp.DriftVelocity(); // cm/us
1344  double timetick = sampling_rate(clockData) * 1e-3; // time sample in us
1345  // rotate coord system CCW around x-axis by pi-thetawire
1346  // new yprime direction is perpendicular to the wire direction
1347  // in the same plane as the wires and in the direction of
1348  // increasing wire number
1349  //use yprime-component of dir cos in rotated coord sys to get
1350  // dTdW (number of time ticks per unit of wire pitch)
1351  double rotang = 3.1416 - thetawire;
1352  double yprime = std::cos(rotang) * startDir.Y() + std::sin(rotang) * startDir.Z();
1353  double dTdW = startDir.X() * wirePitch / driftvelocity / timetick / yprime;
1354 
1355  this->Draw2DSlopeEndPoints(wire, tick, dTdW, color, view);
1356 
1357  // Draw a line to the hit positions, starting from the vertex
1358  size_t nTrackHits = track->NumberTrajectoryPoints();
1359  //TPolyLine& pl = view->AddPolyLine(track->CountValidPoints(),1,1,0); //kColor[id%evd::kNCOLS],1,0);
1360  TPolyLine& pl = view->AddPolyLine(0, 1, 1, 0); //kColor[id%evd::kNCOLS],1,0);
1361 
1362  size_t vidx = 0;
1363  for (size_t idx = 0; idx < nTrackHits; idx++) {
1364  if (track->HasValidPoint(idx) == 0) continue;
1365  const auto& hitPos = track->LocationAtPoint(idx);
1366 
1367  // Use "world" from above
1368  world[1] = hitPos.Y();
1369  world[2] = hitPos.Z();
1370 
1371  // convert the starting position and direction from 3D to 2D coordinates
1372  double tickHit = detProp.ConvertXToTicks(hitPos.X(), plane, t, c);
1373  double wireHit = 0.;
1374  try {
1375  wireHit = 1. * geo->NearestWire(world, plane, t, c);
1376  }
1377  catch (cet::exception& e) {
1378  wireHit = 1. * atoi(e.explain_self().substr(e.explain_self().find("#") + 1, 5).c_str());
1379  }
1380  const size_t tpc = geo->FindTPCAtPosition(hitPos).TPC;
1381  const size_t cryo = geo->FindCryostatAtPosition(hitPos);
1382  if (tpc == t && cryo == c) { pl.SetPoint(vidx++, wireHit, tickHit); }
1383  }
1384  //pl.SetPolyLine(vidx);
1385 
1386  return;
1387  }
1388 
1389  //......................................................................
1390  void
1391  RecoBaseDrawer::Prong2D(const art::Event& evt,
1392  detinfo::DetectorClocksData const& clockData,
1393  detinfo::DetectorPropertiesData const& detProp,
1394  evdb::View2D* view,
1395  unsigned int plane)
1396  {
1397  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1398  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1399  art::ServiceHandle<geo::Geometry const> geo;
1400 
1401  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1402 
1403  geo::View_t gview = geo->TPC(rawOpt->fTPC).Plane(plane).View();
1404 
1405  // annoying for now, but have to have multiple copies of basically the
1406  // same code to draw prongs, showers and tracks so that we can use
1407  // the art::Assns to get the hits and clusters.
1408 
1409  unsigned int cstat = rawOpt->fCryostat;
1410  unsigned int tpc = rawOpt->fTPC;
1411  geo::PlaneID planeID(cstat, tpc, plane);
1412  int tid = 0;
1413 
1414  if (recoOpt->fDrawTracks != 0) {
1415  for (size_t imod = 0; imod < recoOpt->fTrackLabels.size(); ++imod) {
1416  art::InputTag const which = recoOpt->fTrackLabels[imod];
1417 
1418  art::View<recob::Track> track;
1419  this->GetTracks(evt, which, track);
1420 
1421  if (track.vals().size() < 1) continue;
1422 
1423  art::FindMany<recob::Hit> fmh(track, evt, which);
1424 
1425  art::InputTag const whichTag(
1426  recoOpt->fCosmicTagLabels.size() > imod ? recoOpt->fCosmicTagLabels[imod] : "");
1427  art::FindManyP<anab::CosmicTag> cosmicTrackTags(track, evt, whichTag);
1428 
1429  auto tracksProxy = proxy::getCollection<proxy::Tracks>(evt, which);
1430 
1431  // loop over the prongs and get the clusters and hits associated with
1432  // them. only keep those that are in this view
1433  for (size_t t = 0; t < track.vals().size(); ++t) {
1434  // Check for possible issue
1435  if (track.vals().at(t)->NumberTrajectoryPoints() == 0) {
1436  std::cout << "***** Track with no trajectory points ********" << std::endl;
1437  continue;
1438  }
1439 
1440  if (recoOpt->fDrawTracks > 1) {
1441  // BB: draw the track ID at the end of the track
1442  geo::Point_t trackPos(track.vals().at(t)->End().X(),
1443  track.vals().at(t)->End().Y(),
1444  track.vals().at(t)->End().Z());
1445  double tick = 30 + detProp.ConvertXToTicks(trackPos.X(), plane, tpc, cstat);
1446  double wire = geo->WireCoordinate(trackPos, geo::PlaneID(cstat, tpc, plane));
1447  tid =
1448  track.vals().at(t)->ID() &
1449  65535; //this is a hack for PMA track id which uses the 16th bit to identify shower-like track.;
1450  std::string s = std::to_string(tid);
1451  char const* txt = s.c_str();
1452  TText& trkID = view->AddText(wire, tick, txt);
1453  trkID.SetTextColor(evd::kColor[tid % evd::kNCOLS]);
1454  trkID.SetTextSize(0.1);
1455  }
1456 
1457  float Score = -999;
1458  if (cosmicTrackTags.isValid()) {
1459  if (cosmicTrackTags.at(t).size() > 0) {
1460  art::Ptr<anab::CosmicTag> currentTag = cosmicTrackTags.at(t).at(0);
1461  Score = currentTag->CosmicScore();
1462  }
1463  }
1464 
1465  std::vector<const recob::Hit*> hits;
1466  if (track.vals().at(t)->NumberTrajectoryPoints() == fmh.at(t).size()) {
1467  auto tp = tracksProxy[t];
1468  for (auto point : tp.points()) {
1469  if (!point.isPointValid()) continue;
1470  hits.push_back(point.hit());
1471  }
1472  }
1473  else {
1474  hits = fmh.at(t);
1475  }
1476  // only get the hits for the current view
1477  std::vector<const recob::Hit*>::iterator itr = hits.begin();
1478  while (itr < hits.end()) {
1479  if ((*itr)->View() != gview)
1480  hits.erase(itr);
1481  else
1482  itr++;
1483  }
1484 
1485  const recob::Track* aTrack(track.vals().at(t));
1486  int color(evd::kColor[(aTrack->ID() & 65535) % evd::kNCOLS]);
1487  int lineWidth(1);
1488 
1489  if (Score > 0.1 && recoOpt->fDrawCosmicTags) {
1490  color = kRed;
1491  if (Score < 0.6) color = kMagenta;
1492  lineWidth = 3;
1493  }
1494  else if (Score < -10000) { //shower hits
1495  lineWidth = 3;
1496  }
1497 
1498  this->DrawTrack2D(clockData, detProp, hits, view, plane, aTrack, color, lineWidth);
1499  } // end loop over prongs
1500  } // end loop over labels
1501  } // end draw tracks
1502 
1503  if (recoOpt->fDrawShowers != 0) {
1504  static bool first = true;
1505 
1506  if (first) {
1507  std::cout << "DrawShower options: \n";
1508  std::cout << " 1 = Hits in shower color-coded by the shower ID\n";
1509  std::cout << " 2 = Same as 1 + shower axis and circle representing the shower cone\n";
1510  std::cout << " Black cone = shower start dE/dx < 1 MeV/cm (< 1/2 MIP)\n";
1511  std::cout << " Blue cone = shower start dE/dx < 3 MeV/cm (~1 MIP)\n";
1512  std::cout << " Green cone = shower start 3 MeV/cm < dE/dx < 5 MeV/cm (~2 MIP)\n";
1513  std::cout << " Red cone = shower start 5 MeV/cm < dE/dx (>2 MIP)\n";
1514  first = false;
1515  }
1516  for (size_t imod = 0; imod < recoOpt->fShowerLabels.size(); ++imod) {
1517  art::InputTag const which = recoOpt->fShowerLabels[imod];
1518 
1519  art::View<recob::Shower> shower;
1520  this->GetShowers(evt, which, shower);
1521  if (shower.vals().size() < 1) continue;
1522 
1523  art::FindMany<recob::Hit> fmh(shower, evt, which);
1524 
1525  // loop over the prongs and get the clusters and hits associated with
1526  // them. only keep those that are in this view
1527  for (size_t s = 0; s < shower.vals().size(); ++s) {
1528 
1529  std::vector<const recob::Hit*> hits = fmh.at(s);
1530  // only get the hits for the current view
1531  std::vector<const recob::Hit*>::iterator itr = hits.begin();
1532  while (itr < hits.end()) {
1533  if ((*itr)->View() != gview)
1534  hits.erase(itr);
1535  else
1536  itr++;
1537  }
1538  if (recoOpt->fDrawShowers > 1) {
1539  // BB draw a line between the start and end points and a "circle" that represents
1540  // the shower cone angle at the end point
1541  if (!shower.vals().at(s)->has_length()) continue;
1542  if (!shower.vals().at(s)->has_open_angle()) continue;
1543 
1544  TVector3 startPos = shower.vals().at(s)->ShowerStart();
1545  TVector3 dir = shower.vals().at(s)->Direction();
1546  double length = shower.vals().at(s)->Length();
1547  double openAngle = shower.vals().at(s)->OpenAngle();
1548 
1549  // Find the center of the cone base
1550  TVector3 endPos = startPos + length * dir;
1551 
1552  geo::Point_t localStart(startPos);
1553  geo::Point_t localEnd(endPos);
1554 
1555  double swire = geo->WireCoordinate(localStart, planeID);
1556  double stick = detProp.ConvertXToTicks(startPos.X(), planeID);
1557  double ewire = geo->WireCoordinate(localEnd, planeID);
1558  double etick = detProp.ConvertXToTicks(endPos.X(), planeID);
1559  TLine& coneLine = view->AddLine(swire, stick, ewire, etick);
1560  // color coding by dE/dx
1561  std::vector<double> dedxVec = shower.vals().at(s)->dEdx();
1562  // float dEdx = shower.vals().at(s)->dEdx()[plane];
1563  // use black for too-low dE/dx
1564  int color = kBlack;
1565  if (plane < dedxVec.size()) {
1566  if (dedxVec[plane] > 1 && dedxVec[plane] < 3) {
1567  // use blue for ~1 MIP
1568  color = kBlue;
1569  }
1570  else if (dedxVec[plane] < 5) {
1571  // use green for ~2 MIP
1572  color = kGreen;
1573  }
1574  else {
1575  // use red for >~ 2 MIP
1576  color = kRed;
1577  }
1578  }
1579  coneLine.SetLineColor(color);
1580 
1581  // Now find the 3D circle that represents the base of the cone
1582  double radius = length * openAngle;
1583  auto coneRim = Circle3D(endPos, dir, radius);
1584  TPolyLine& pline = view->AddPolyLine(coneRim.size(), color, 2, 0);
1585  // project these points into the plane
1586  for (unsigned short ipt = 0; ipt < coneRim.size(); ++ipt) {
1587  geo::Point_t localPos(coneRim[ipt][0], coneRim[ipt][1], coneRim[ipt][2]);
1588 
1589  double wire = geo->WireCoordinate(localPos, planeID);
1590  double tick = detProp.ConvertXToTicks(coneRim[ipt][0], planeID);
1591  pline.SetPoint(ipt, wire, tick);
1592  } // ipt
1593  }
1594  this->DrawProng2D(detProp,
1595  hits,
1596  view,
1597  plane,
1598  shower.vals().at(s)->ShowerStart(),
1599  shower.vals().at(s)->Direction(),
1600  s,
1601  -10001); //use -10001 to increase shower hit size
1602 
1603  } // end loop over prongs
1604  } // end loop over labels
1605  } // end draw showers
1606 
1607  return;
1608  }
1609 
1610  //......................................................................
1611  void
1612  RecoBaseDrawer::DrawTrackVertexAssns2D(const art::Event& evt,
1613  detinfo::DetectorClocksData const& clockData,
1614  detinfo::DetectorPropertiesData const& detProp,
1615  evdb::View2D* view,
1616  unsigned int plane)
1617  {
1618  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1619  if (!recoOpt->fDrawTrackVertexAssns) return;
1620 
1621  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1622  art::ServiceHandle<geo::Geometry const> geo;
1623 
1624  geo::View_t gview = geo->TPC(rawOpt->fTPC).Plane(plane).View();
1625 
1626  // annoying for now, but have to have multiple copies of basically the
1627  // same code to draw prongs, showers and tracks so that we can use
1628  // the art::Assns to get the hits and clusters.
1629 
1630  unsigned int cstat = rawOpt->fCryostat;
1631  unsigned int tpc = rawOpt->fTPC;
1632  geo::PlaneID planeID(cstat, tpc, plane);
1633  int tid = 0;
1634 
1635  for (size_t imod = 0; imod < recoOpt->fTrkVtxTrackLabels.size(); ++imod) {
1636  art::InputTag const which = recoOpt->fTrkVtxTrackLabels[imod];
1637 
1638  art::View<recob::Track> trackCol;
1639  this->GetTracks(evt, which, trackCol);
1640 
1641  if (trackCol.vals().size() < 1) continue;
1642 
1643  // Recover associations output from the filter
1644  std::unique_ptr<art::Assns<recob::Vertex, recob::Track>> vertexTrackAssociations(
1645  new art::Assns<recob::Vertex, recob::Track>);
1646 
1647  // Recover a handle to the collection of associations between vertices and tracks
1648  // This is a bit non-standard way to do this but trying to avoid complications
1649  art::Handle<art::Assns<recob::Vertex, recob::Track>> vertexTrackAssnsHandle;
1650 
1651  evt.getByLabel(recoOpt->fTrkVtxFilterLabels[imod], vertexTrackAssnsHandle);
1652 
1653  if (vertexTrackAssnsHandle->size() < 1) continue;
1654 
1655  // Get the rest of the associations in the standard way
1656  art::FindMany<recob::Hit> fmh(trackCol, evt, which);
1657 
1658  art::FindManyP<anab::CosmicTag> cosmicTrackTags(
1659  trackCol, evt, recoOpt->fTrkVtxCosmicLabels[imod]);
1660 
1661  auto tracksProxy = proxy::getCollection<proxy::Tracks>(evt, which);
1662 
1663  // Need to keep track of vertices unfortunately
1664  int lastVtxIdx(-1);
1665  int color(kRed);
1666 
1667  std::cout << "==> Neutrino Candidate drawing for tagger "
1668  << recoOpt->fTrkVtxFilterLabels[imod] << std::endl;
1669 
1670  // Now we can iterate over the vertex/track associations and do some drawing
1671  for (const auto& vertexTrackAssn : *vertexTrackAssnsHandle) {
1672  // Start by drawing the vertex
1673  art::Ptr<recob::Vertex> vertex = vertexTrackAssn.first;
1674 
1675  if (vertex->ID() != lastVtxIdx) {
1676  // BB: draw polymarker at the vertex position in this plane
1677  double xyz[3];
1678 
1679  vertex->XYZ(xyz);
1680 
1681  geo::Point_t localXYZ(xyz[0], xyz[1], xyz[2]);
1682 
1683  double wire = geo->WireCoordinate(localXYZ, planeID);
1684  double time = detProp.ConvertXToTicks(xyz[0], planeID);
1685 
1686  TMarker& strt = view->AddMarker(wire, time, color, 24, 3.0);
1687  strt.SetMarkerColor(color);
1688 
1689  std::cout << " --> Drawing vertex id: " << vertex->ID() << std::endl;
1690  }
1691 
1692  lastVtxIdx = vertex->ID();
1693 
1694  const art::Ptr<recob::Track>& track = vertexTrackAssn.second;
1695 
1696  // BB: draw the track ID at the end of the track
1697  double x = track->End().X();
1698  geo::Point_t trackEnd(track->End());
1699  double tick = 30 + detProp.ConvertXToTicks(x, planeID);
1700  double wire = geo->WireCoordinate(trackEnd, planeID);
1701 
1702  tid = track->ID() & 65535;
1703 
1704  std::cout << " --> Drawing Track id: " << tid << std::endl;
1705 
1706  std::string s = std::to_string(tid);
1707  char const* txt = s.c_str();
1708 
1709  TText& trkID = view->AddText(wire, tick, txt);
1710  trkID.SetTextColor(color);
1711  trkID.SetTextSize(0.1);
1712 
1713  float cosmicScore = -999;
1714  if (cosmicTrackTags.isValid()) {
1715  if (cosmicTrackTags.at(track.key()).size() > 0) {
1716  art::Ptr<anab::CosmicTag> currentTag = cosmicTrackTags.at(track.key()).at(0);
1717  cosmicScore = currentTag->CosmicScore();
1718  }
1719  }
1720 
1721  std::vector<const recob::Hit*> hits;
1722  if (track->NumberTrajectoryPoints() == fmh.at(track.key()).size()) {
1723  auto tp = tracksProxy[track.key()];
1724  for (auto point : tp.points()) {
1725  if (!point.isPointValid()) continue;
1726  hits.push_back(point.hit());
1727  }
1728  }
1729  else {
1730  hits = fmh.at(track.key());
1731  }
1732  // only get the hits for the current view
1733  std::vector<const recob::Hit*>::iterator itr = hits.begin();
1734  while (itr < hits.end()) {
1735  if ((*itr)->View() != gview)
1736  hits.erase(itr);
1737  else
1738  itr++;
1739  }
1740 
1741  int lineWidth(1);
1742 
1743  if (cosmicScore > 0.1) {
1744  color = kRed;
1745  if (cosmicScore < 0.6) color = kMagenta;
1746  lineWidth = 3;
1747  }
1748  else if (cosmicScore < -10000) { //shower hits
1749  lineWidth = 3;
1750  }
1751 
1752  this->DrawTrack2D(clockData, detProp, hits, view, plane, track.get(), color, lineWidth);
1753 
1754  } // end loop over vertex/track associations
1755 
1756  } // end loop over labels
1757  }
1758 
1759  //......................................................................
1760  void
1761  RecoBaseDrawer::Vertex2D(const art::Event& evt,
1762  detinfo::DetectorPropertiesData const& detProp,
1763  evdb::View2D* view,
1764  unsigned int plane)
1765  {
1766  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1767  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1768 
1769  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1770  if (recoOpt->fDrawVertices == 0) return;
1771 
1772  art::ServiceHandle<geo::Geometry const> geo;
1773  static bool first = true;
1774 
1775  if (first) {
1776  std::cout << "******** DrawVertices: Open circles color coded across all planes. Set "
1777  "DrawVertices > 1 to display the vertex ID\n";
1778  first = false;
1779  }
1780 
1781  for (size_t imod = 0; imod < recoOpt->fVertexLabels.size(); ++imod) {
1782  art::InputTag const which = recoOpt->fVertexLabels[imod];
1783 
1784  art::PtrVector<recob::Vertex> vertex;
1785  this->GetVertices(evt, which, vertex);
1786 
1787  if (vertex.size() < 1) continue;
1788 
1789  double local[3] = {0., 0., 0.};
1790  double world[3] = {0., 0., 0.};
1791  const geo::TPCGeo& tpc = geo->TPC(rawOpt->fTPC);
1792  tpc.LocalToWorld(local, world);
1793  double minxyz[3], maxxyz[3];
1794  minxyz[0] = world[0] - geo->DetHalfWidth(rawOpt->fTPC, rawOpt->fCryostat);
1795  maxxyz[0] = world[0] + geo->DetHalfWidth(rawOpt->fTPC, rawOpt->fCryostat);
1796  minxyz[1] = world[1] - geo->DetHalfWidth(rawOpt->fTPC, rawOpt->fCryostat);
1797  maxxyz[1] = world[1] + geo->DetHalfWidth(rawOpt->fTPC, rawOpt->fCryostat);
1798  minxyz[2] = world[2] - geo->DetLength(rawOpt->fTPC, rawOpt->fCryostat) / 2;
1799  maxxyz[2] = world[2] + geo->DetLength(rawOpt->fTPC, rawOpt->fCryostat) / 2;
1800 
1801  for (size_t v = 0; v < vertex.size(); ++v) {
1802  // ensure the vertex is inside the current tpc
1803  double xyz[3];
1804  vertex[v]->XYZ(xyz);
1805  if (xyz[0] < minxyz[0] || xyz[0] > maxxyz[0]) continue;
1806  if (xyz[1] < minxyz[1] || xyz[1] > maxxyz[1]) continue;
1807  if (xyz[2] < minxyz[2] || xyz[2] > maxxyz[2]) continue;
1808 
1809  geo::Point_t localPos(xyz[0], xyz[1], xyz[2]);
1810 
1811  // BB: draw polymarker at the vertex position in this plane
1812  double wire =
1813  geo->WireCoordinate(localPos, geo::PlaneID(rawOpt->fCryostat, rawOpt->fTPC, plane));
1814  double time = detProp.ConvertXToTicks(xyz[0], plane, rawOpt->fTPC, rawOpt->fCryostat);
1815  int color = evd::kColor[vertex[v]->ID() % evd::kNCOLS];
1816  TMarker& strt = view->AddMarker(wire, time, color, 24, 1.0);
1817  strt.SetMarkerColor(color);
1818 
1819  // BB: draw the vertex ID
1820  if (recoOpt->fDrawVertices > 1) {
1821  std::string s = "3V" + std::to_string(vertex[v]->ID());
1822  char const* txt = s.c_str();
1823  TText& vtxID = view->AddText(wire, time + 30, txt);
1824  vtxID.SetTextColor(color);
1825  vtxID.SetTextSize(0.05);
1826  }
1827  } // end loop over vertices to draw from this label
1828  } // end loop over vertex module lables
1829 
1830  return;
1831  }
1832 
1833  //......................................................................
1834  void
1835  RecoBaseDrawer::Event2D(const art::Event& evt, evdb::View2D* view, unsigned int plane)
1836  {
1837  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1838  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1839  art::ServiceHandle<geo::Geometry const> geo;
1840 
1841  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1842 
1843  if (recoOpt->fDrawEvents != 0) {
1844  geo::View_t gview = geo->TPC(rawOpt->fTPC).Plane(plane).View();
1845 
1846  for (unsigned int imod = 0; imod < recoOpt->fEventLabels.size(); ++imod) {
1847  art::InputTag const which = recoOpt->fEventLabels[imod];
1848 
1849  art::PtrVector<recob::Event> event;
1850  this->GetEvents(evt, which, event);
1851 
1852  if (event.size() < 1) continue;
1853 
1854  art::FindMany<recob::Hit> fmh(event, evt, which);
1855 
1856  for (size_t e = 0; e < event.size(); ++e) {
1857  std::vector<const recob::Hit*> hits;
1858 
1859  hits = fmh.at(e);
1860 
1861  // only get the hits for the current view
1862  std::vector<const recob::Hit*>::iterator itr = hits.begin();
1863  while (itr < hits.end()) {
1864  if ((*itr)->View() != gview)
1865  hits.erase(itr);
1866  else
1867  itr++;
1868  }
1869 
1870  this->Hit2D(hits, evd::kColor[event[e]->ID() % evd::kNCOLS], view, false, true);
1871  } // end loop over events
1872  } // end loop over event module lables
1873  } // end if we are drawing events
1874 
1875  return;
1876  }
1877 
1878  //......................................................................
1879  void
1880  RecoBaseDrawer::Seed3D(const art::Event& evt, evdb::View3D* view)
1881  {
1882  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1883  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1884 
1885  std::vector<art::InputTag> labels;
1886  if (recoOpt->fDrawSeeds != 0)
1887  for (size_t imod = 0; imod < recoOpt->fSeedLabels.size(); ++imod)
1888  labels.push_back(recoOpt->fSeedLabels[imod]);
1889 
1890  for (size_t imod = 0; imod < labels.size(); ++imod) {
1891  art::InputTag const which = labels[imod];
1892 
1893  art::PtrVector<recob::Seed> seeds;
1894  this->GetSeeds(evt, which, seeds);
1895 
1896  int color = 0;
1897 
1898  if (seeds.size() < 1) continue;
1899 
1900  TPolyMarker3D& pmrk = view->AddPolyMarker3D(seeds.size(), color, 4, 1);
1901 
1902  for (size_t iseed = 0; iseed != seeds.size(); ++iseed) {
1903  double pt[3], pterr[3], dir[3], direrr[3];
1904  seeds.at(iseed)->GetPoint(pt, pterr);
1905  seeds.at(iseed)->GetDirection(dir, direrr);
1906 
1907  double end1[3], end2[3];
1908  for (int i = 0; i != 3; ++i) {
1909  end1[i] = pt[i] + dir[i];
1910  end2[i] = pt[i] - dir[i];
1911  }
1912 
1913  TPolyLine3D& pline = view->AddPolyLine3D(2, color, 2, 0);
1914 
1915  pmrk.SetPoint(iseed, pt[0], pt[1], pt[2]);
1916  pline.SetPoint(0, end1[0], end1[1], end1[2]);
1917  pline.SetPoint(1, end2[0], end2[1], end2[2]);
1918  } // end loop over seeds
1919  } // end loop over module labels
1920 
1921  return;
1922  }
1923 
1924  //......................................................................
1925  void
1926  RecoBaseDrawer::SeedOrtho(const art::Event& evt, evd::OrthoProj_t proj, evdb::View2D* view)
1927  {
1928  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1929  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1930 
1931  std::vector<art::InputTag> labels;
1932  if (recoOpt->fDrawSeeds != 0)
1933  for (size_t imod = 0; imod < recoOpt->fSeedLabels.size(); ++imod)
1934  labels.push_back(recoOpt->fSeedLabels[imod]);
1935 
1936  for (size_t imod = 0; imod < labels.size(); ++imod) {
1937  art::InputTag const which = labels[imod];
1938 
1939  art::PtrVector<recob::Seed> seeds;
1940  this->GetSeeds(evt, which, seeds);
1941 
1942  int color = 0;
1943 
1944  for (size_t iseed = 0; iseed != seeds.size(); ++iseed) {
1945  double pt[3], pterr[3], dir[3], direrr[3];
1946  seeds.at(iseed)->GetPoint(pt, pterr);
1947  seeds.at(iseed)->GetDirection(dir, direrr);
1948 
1949  double end1[3], end2[3];
1950  for (int i = 0; i != 3; ++i) {
1951  end1[i] = pt[i] + dir[i];
1952  end2[i] = pt[i] - dir[i];
1953  }
1954 
1955  if (proj == evd::kXY) {
1956  TMarker& strt = view->AddMarker(pt[1], pt[0], color, 4, 1.5);
1957  TLine& line = view->AddLine(end1[1], end1[0], end2[1], end2[0]);
1958  strt.SetMarkerColor(evd::kColor[color]);
1959  line.SetLineColor(evd::kColor[color]);
1960  line.SetLineWidth(2.0);
1961  }
1962  else if (proj == evd::kXZ) {
1963  TMarker& strt = view->AddMarker(pt[2], pt[0], color, 4, 1.5);
1964  TLine& line = view->AddLine(end1[2], end1[0], end2[2], end2[0]);
1965  strt.SetMarkerColor(evd::kColor[color]);
1966  line.SetLineColor(evd::kColor[color]);
1967  line.SetLineWidth(2.0);
1968  }
1969  else {
1970  if (proj != evd::kYZ)
1971  throw cet::exception("RecoBaseDrawer:SeedOrtho")
1972  << "projection is not YZ as expected\n";
1973 
1974  TMarker& strt = view->AddMarker(pt[2], pt[1], color, 4, 1.5);
1975  TLine& line = view->AddLine(end1[2], end1[1], end2[2], end2[1]);
1976  strt.SetMarkerColor(evd::kColor[color]);
1977  line.SetLineColor(evd::kColor[color]);
1978  line.SetLineWidth(2.0);
1979  }
1980  }
1981  }
1982  }
1983 
1984  //......................................................................
1985  void
1986  RecoBaseDrawer::SpacePoint3D(const art::Event& evt, evdb::View3D* view)
1987  {
1988  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
1989  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
1990 
1991  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
1992 
1993  std::vector<art::InputTag> labels;
1994  if (recoOpt->fDrawSpacePoints != 0) {
1995  for (size_t imod = 0; imod < recoOpt->fSpacePointLabels.size(); ++imod)
1996  labels.push_back(recoOpt->fSpacePointLabels[imod]);
1997  }
1998 
1999  for (size_t imod = 0; imod < labels.size(); ++imod) {
2000  art::InputTag const which = labels[imod];
2001 
2002  std::vector<art::Ptr<recob::SpacePoint>> spts;
2003  this->GetSpacePoints(evt, which, spts);
2004  int color = 10 * imod + 11;
2005 
2006  color = 0;
2007 
2008  // std::vector<const recob::SpacePoint*> sptsVec;
2009  //
2010  // sptsVec.resize(spts.size());
2011  // for(const auto& spt : spts){
2012  // std::cout<<spt<<" "<<*spt<<" "<<&*spt<<std::endl;
2013  // sptsVec.push_back(&*spt);
2014  // std::cout<<sptsVec.back()<<std::endl;
2015  // }
2016  fAllSpacePointDrawer->Draw(spts, view, color, kFullDotMedium, 1);
2017  }
2018 
2019  return;
2020  }
2021 
2022  //......................................................................
2023  void
2024  RecoBaseDrawer::PFParticle3D(const art::Event& evt, evdb::View3D* view)
2025  {
2026  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2027  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2028 
2029  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
2030  if (recoOpt->fDrawPFParticles < 1) return;
2031 
2032  // The plan is to loop over the list of possible particles
2033  for (size_t imod = 0; imod < recoOpt->fPFParticleLabels.size(); ++imod) {
2034  art::InputTag const which = recoOpt->fPFParticleLabels[imod];
2035  art::InputTag const assns = recoOpt->fSpacePointLabels[imod];
2036 
2037  // Start off by recovering our 3D Clusters for this label
2038  art::PtrVector<recob::PFParticle> pfParticleVec;
2039  this->GetPFParticles(evt, which, pfParticleVec);
2040 
2041  mf::LogDebug("RecoBaseDrawer")
2042  << "RecoBaseDrawer: number PFParticles to draw: " << pfParticleVec.size() << std::endl;
2043 
2044  // Make sure we have some clusters
2045  if (pfParticleVec.size() < 1) continue;
2046 
2047  // Get the space points created by the PFParticle producer
2048  std::vector<art::Ptr<recob::SpacePoint>> spacePointVec;
2049  this->GetSpacePoints(evt, assns, spacePointVec);
2050 
2051  // Recover the edges
2052  std::vector<art::Ptr<recob::Edge>> edgeVec;
2053  if (recoOpt->fDrawEdges) this->GetEdges(evt, assns, edgeVec);
2054 
2055  // No space points no continue
2056  if (spacePointVec.empty()) continue;
2057 
2058  // Add the relations to recover associations cluster hits
2059  art::FindManyP<recob::SpacePoint> edgeSpacePointAssnsVec(edgeVec, evt, assns);
2060  art::FindManyP<recob::SpacePoint> spacePointAssnVec(pfParticleVec, evt, assns);
2061  art::FindManyP<recob::Hit> spHitAssnVec(spacePointVec, evt, assns);
2062  art::FindManyP<recob::Edge> edgeAssnsVec(pfParticleVec, evt, assns);
2063 
2064  // If no valid space point associations then nothing to do
2065  if (!spacePointAssnVec.isValid()) continue;
2066 
2067  // Need the PCA info as well
2068  art::FindMany<recob::PCAxis> pcAxisAssnVec(pfParticleVec, evt, which);
2069 
2070  // Want CR tagging info
2071  // Note the cosmic tags come from a different producer - we assume that the producers are
2072  // matched in the fcl label vectors!
2073  art::InputTag cosmicTagLabel =
2074  imod < recoOpt->fCosmicTagLabels.size() ? recoOpt->fCosmicTagLabels[imod] : "";
2075  art::FindMany<anab::CosmicTag> pfCosmicAssns(pfParticleVec, evt, cosmicTagLabel);
2076 
2077  // We also want to drive display of tracks but have the same issue with production... so follow the
2078  // same prescription.
2079  art::InputTag trackTagLabel =
2080  imod < recoOpt->fTrackLabels.size() ? recoOpt->fTrackLabels[imod] : "";
2081  art::FindMany<recob::Track> pfTrackAssns(pfParticleVec, evt, trackTagLabel);
2082 
2083  // Commence looping over possible clusters
2084  for (size_t idx = 0; idx < pfParticleVec.size(); idx++) {
2085  // Recover cluster
2086  const art::Ptr<recob::PFParticle> pfParticle(pfParticleVec.at(idx));
2087 
2088  // Drawing will be done recursively down the chain of hieirarchy... so we want to begin
2089  // with only "primary" particles, if we find one that isn't then we skip
2090  if (!pfParticle->IsPrimary()) continue;
2091 
2092  // Call the recursive drawing routine
2093  DrawPFParticle3D(pfParticle,
2094  pfParticleVec,
2095  spacePointVec,
2096  edgeAssnsVec,
2097  spacePointAssnVec,
2098  edgeSpacePointAssnsVec,
2099  spHitAssnVec,
2100  pfTrackAssns,
2101  pcAxisAssnVec,
2102  pfCosmicAssns,
2103  0,
2104  view);
2105  }
2106  }
2107 
2108  return;
2109  }
2110 
2111  float
2112  RecoBaseDrawer::SpacePointChiSq(const std::vector<art::Ptr<recob::Hit>>& hitVec) const
2113  {
2114  float hitChiSq(0.);
2115 
2116  bool usePlane[] = {false, false, false};
2117  float peakTimeVec[] = {0., 0., 0.};
2118  float peakSigmaVec[] = {0., 0., 0.};
2119  float aveSum(0.);
2120  float weightSum(0.);
2121 
2122  // Temp ad hoc correction to investigate...
2123  std::map<size_t, double> planeOffsetMap;
2124 
2125  planeOffsetMap[0] = 0.;
2126  planeOffsetMap[1] = 4.;
2127  planeOffsetMap[2] = 8.;
2128 
2129  for (const auto& hit : hitVec) {
2130  if (!hit) continue;
2131 
2132  float peakTime = hit->PeakTime() - planeOffsetMap[hit->WireID().Plane];
2133  float peakRMS = hit->RMS();
2134 
2135  aveSum += peakTime / (peakRMS * peakRMS);
2136  weightSum += 1. / (peakRMS * peakRMS);
2137 
2138  peakTimeVec[hit->WireID().Plane] = peakTime;
2139  peakSigmaVec[hit->WireID().Plane] = peakRMS;
2140  usePlane[hit->WireID().Plane] = true;
2141  }
2142 
2143  aveSum /= weightSum;
2144 
2145  for (int idx = 0; idx < 3; idx++) {
2146  if (usePlane[idx]) {
2147  float deltaTime = peakTimeVec[idx] - aveSum;
2148  float sigmaPeakTimeSq = peakSigmaVec[idx] * peakSigmaVec[idx];
2149 
2150  hitChiSq += deltaTime * deltaTime / sigmaPeakTimeSq;
2151  }
2152  }
2153 
2154  return hitChiSq;
2155  }
2156 
2157  void
2158  RecoBaseDrawer::DrawPFParticle3D(const art::Ptr<recob::PFParticle>& pfPart,
2159  const art::PtrVector<recob::PFParticle>& pfParticleVec,
2160  const std::vector<art::Ptr<recob::SpacePoint>>& spacePointVec,
2161  const art::FindManyP<recob::Edge>& edgeAssnsVec,
2162  const art::FindManyP<recob::SpacePoint>& spacePointAssnVec,
2163  const art::FindManyP<recob::SpacePoint>& edgeSPAssnVec,
2164  const art::FindManyP<recob::Hit>& spHitAssnVec,
2165  const art::FindMany<recob::Track>& trackAssnVec,
2166  const art::FindMany<recob::PCAxis>& pcAxisAssnVec,
2167  const art::FindMany<anab::CosmicTag>& cosmicTagAssnVec,
2168  int depth,
2169  evdb::View3D* view)
2170  {
2171  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2172  art::ServiceHandle<evd::ColorDrawingOptions const> cst;
2173 
2174  // First let's draw the hits associated to this cluster
2175  const std::vector<art::Ptr<recob::SpacePoint>>& hitsVec(spacePointAssnVec.at(pfPart.key()));
2176 
2177  // Use the particle ID to determine the color to draw the points
2178  // Ok, this is what we would like to do eventually but currently all particles are the same...
2179  bool isCosmic(false);
2180  int colorIdx(evd::kColor[pfPart->Self() % evd::kNCOLS]);
2181 
2182  // Recover cosmic tag info if any
2183  if (cosmicTagAssnVec.isValid() && recoOpt->fDrawPFParticles > 3) {
2184  std::vector<const anab::CosmicTag*> pfCosmicTagVec = cosmicTagAssnVec.at(pfPart.key());
2185 
2186  if (!pfCosmicTagVec.empty()) {
2187  const anab::CosmicTag* cosmicTag = pfCosmicTagVec.front();
2188 
2189  if (cosmicTag->CosmicScore() > 0.6) isCosmic = true;
2190  }
2191  }
2192 
2193  // Reset color index if a cosmic
2194  if (isCosmic) colorIdx = 12;
2195 
2196  if (!hitsVec.empty() && recoOpt->fDraw3DSpacePoints)
2197  fSpacePointDrawer->Draw(hitsVec, view, 1, kFullDotLarge, 0.25, &spHitAssnVec);
2198  /*
2199  {
2200  using HitPosition = std::array<double,6>;
2201  std::map<int,std::vector<HitPosition>> colorToHitMap;
2202 
2203  float minHitChiSquare(0.);
2204  float maxHitChiSquare(2.);
2205  float hitChiSqScale((cst->fRecoQHigh[geo::kCollection] - cst->fRecoQLow[geo::kCollection]) / (maxHitChiSquare - minHitChiSquare));
2206 
2207  for(const auto& spacePoint : hitsVec)
2208  {
2209  const double* pos = spacePoint->XYZ();
2210  const double* err = spacePoint->ErrXYZ();
2211 
2212  bool storeHit(false);
2213  int chargeColorIdx(0);
2214  float spacePointChiSq(spacePoint->Chisq());
2215 
2216  if (recoOpt->fDraw3DSpacePointHeatMap)
2217  {
2218  storeHit = true;
2219 
2220  float hitChiSq = std::max(minHitChiSquare, std::min(maxHitChiSquare, spacePointChiSq));
2221 
2222  float chgFactor = cst->fRecoQHigh[geo::kCollection] - hitChiSqScale * hitChiSq;
2223  //float chgFactor = delTScaleFctr * hitChiSq + cst->fRecoQLow[geo::kCollection];
2224 
2225  chargeColorIdx = cst->CalQ(geo::kCollection).GetColor(chgFactor);
2226  }
2227  else
2228  {
2229  if (spacePointChiSq > 0. && !recoOpt->fSkeletonOnly) // All cluster hits which are unmarked
2230  {
2231  storeHit = true;
2232  }
2233  else if (spacePointChiSq > -2.) // Skeleton hits
2234  {
2235  chargeColorIdx = 5;
2236  storeHit = true;
2237  }
2238  else if (spacePointChiSq > -3.) // Pure edge hits
2239  {
2240  if (chargeColorIdx < 0) chargeColorIdx = !isCosmic ? 3 : colorIdx + 3;
2241  storeHit = true;
2242  }
2243  else if (spacePointChiSq > -4.) // Skeleton hits which are also edge hits
2244  {
2245  if (chargeColorIdx < 0) chargeColorIdx = !isCosmic ? 0 : colorIdx + 3;
2246  storeHit = true;
2247  }
2248  else if (spacePoint->Chisq() > -5.) // Hits which form seeds for tracks
2249  {
2250  if (chargeColorIdx < 0) chargeColorIdx = !isCosmic ? 5 : colorIdx + 3;
2251  storeHit = true;
2252  }
2253  else if (!recoOpt->fSkeletonOnly) // hits made from pairs
2254  {
2255  chargeColorIdx = 15;
2256  storeHit = true;
2257  }
2258 
2259  if (chargeColorIdx < 0) chargeColorIdx = 0;
2260  }
2261 
2262  if (storeHit) colorToHitMap[chargeColorIdx].push_back(HitPosition()={{pos[0],pos[1],pos[2],err[3],err[3],err[5]}});
2263  }
2264 
2265  size_t nHitsDrawn(0);
2266 
2267  for(auto& hitPair : colorToHitMap)
2268  {
2269  //TPolyMarker3D& pm = view->AddPolyMarker3D(hitPair.second.size(), hitPair.first, kFullDotMedium, 3);
2270  TPolyMarker3D& pm = view->AddPolyMarker3D(hitPair.second.size(), hitPair.first, kFullDotLarge, 0.25); //kFullDotLarge, 0.3);
2271  for (const auto& hit : hitPair.second) pm.SetNextPoint(hit[0],hit[1],hit[2]);
2272 
2273  nHitsDrawn += hitPair.second.size();
2274  }
2275  }
2276 */
2277 
2278  // Now try to draw any associated edges
2279  if (edgeAssnsVec.isValid() && recoOpt->fDraw3DEdges) {
2280  const std::vector<art::Ptr<recob::Edge>>& edgeVec(edgeAssnsVec.at(pfPart.key()));
2281 
2282  if (!edgeVec.empty()) {
2283  TPolyMarker3D& pm = view->AddPolyMarker3D(
2284  2 * edgeVec.size(), colorIdx, kFullDotMedium, 1.25); //kFullDotLarge, 0.5);
2285 
2286  for (const auto& edge : edgeVec) {
2287  try {
2288  const std::vector<art::Ptr<recob::SpacePoint>>& spacePointVec(
2289  edgeSPAssnVec.at(edge.key()));
2290 
2291  if (spacePointVec.size() != 2) {
2292  std::cout << "Space Point vector associated to edge is not of size 2: "
2293  << spacePointVec.size() << std::endl;
2294  continue;
2295  }
2296 
2297  const recob::SpacePoint* firstSP = spacePointVec[0].get();
2298  const recob::SpacePoint* secondSP = spacePointVec[1].get();
2299 
2300  TVector3 startPoint(firstSP->XYZ()[0], firstSP->XYZ()[1], firstSP->XYZ()[2]);
2301  TVector3 endPoint(secondSP->XYZ()[0], secondSP->XYZ()[1], secondSP->XYZ()[2]);
2302  TVector3 lineVec(endPoint - startPoint);
2303 
2304  pm.SetNextPoint(startPoint[0], startPoint[1], startPoint[2]);
2305  pm.SetNextPoint(endPoint[0], endPoint[1], endPoint[2]);
2306 
2307  double length = lineVec.Mag();
2308 
2309  if (length == 0.) {
2310  std::cout << "Edge length is zero, index 1: " << edge->FirstPointID()
2311  << ", index 2: " << edge->SecondPointID() << std::endl;
2312  continue;
2313  }
2314 
2315  double minLen = std::max(2.01, length);
2316 
2317  if (minLen > length) {
2318  lineVec.SetMag(1.);
2319 
2320  startPoint += -0.5 * (minLen - length) * lineVec;
2321  endPoint += 0.5 * (minLen - length) * lineVec;
2322  }
2323 
2324  // Get a polyline object to draw from the first to the second space point
2325  TPolyLine3D& pl = view->AddPolyLine3D(2, colorIdx, 4, 1);
2326 
2327  pl.SetPoint(0, startPoint[0], startPoint[1], startPoint[2]);
2328  pl.SetPoint(1, endPoint[0], endPoint[1], endPoint[2]);
2329  }
2330  catch (...) {
2331  continue;
2332  }
2333  }
2334  }
2335  }
2336 
2337  // Draw associated tracks
2338  if (trackAssnVec.isValid()) {
2339  std::vector<const recob::Track*> trackVec(trackAssnVec.at(pfPart.key()));
2340 
2341  if (!trackVec.empty()) {
2342  for (const auto& track : trackVec)
2343  DrawTrack3D(*track, view, colorIdx, kFullDotLarge, 0.5);
2344  }
2345  }
2346 
2347  // Look up the PCA info
2348  if (pcAxisAssnVec.isValid() && recoOpt->fDraw3DPCAAxes) {
2349  std::vector<const recob::PCAxis*> pcaVec(pcAxisAssnVec.at(pfPart.key()));
2350 
2351  if (!pcaVec.empty()) {
2352  // For each axis we are going to draw a solid line between two points
2353  int numPoints(2);
2354  //int lineWidth[2] = { 3, 1};
2355  int lineWidth[2] = {2, 1};
2356  int lineStyle[2] = {1, 13};
2357  int lineColor[2] = {colorIdx, 18};
2358  //int markStyle[2] = { 4, 4};
2359  int markStyle[2] = {kFullDotLarge, kFullDotLarge};
2360  double markSize[2] = {0.5, 0.2};
2361  int pcaIdx(0);
2362 
2363  if (!isCosmic) lineColor[1] = colorIdx;
2364 
2365  // The order of axes in the returned association vector is arbitrary... the "first" axis is
2366  // better and we can divine that by looking at the axis id's (the best will have been made first)
2367  if (pcaVec.size() > 1 && pcaVec.front()->getID() > pcaVec.back()->getID())
2368  std::reverse(pcaVec.begin(), pcaVec.end());
2369 
2370  for (const auto& pca : pcaVec) {
2371  // We need the mean position
2372  const double* avePosition = pca->getAvePosition();
2373 
2374  // Let's draw a marker at the interesting points
2375  int pmrkIdx(0);
2376  TPolyMarker3D& pmrk =
2377  view->AddPolyMarker3D(7, lineColor[pcaIdx], markStyle[pcaIdx], markSize[pcaIdx]);
2378 
2379  pmrk.SetPoint(pmrkIdx++, avePosition[0], avePosition[1], avePosition[2]);
2380 
2381  // Loop over pca dimensions
2382  for (int dimIdx = 0; dimIdx < 3; dimIdx++) {
2383  // Oh please oh please give me an instance of a poly line...
2384  TPolyLine3D& pl = view->AddPolyLine3D(
2385  numPoints, lineColor[pcaIdx], lineWidth[pcaIdx], lineStyle[pcaIdx]);
2386 
2387  // We will use the eigen value to give the length of the line we're going to plot
2388  double eigenValue = pca->getEigenValues()[dimIdx];
2389 
2390  // Make sure a valid eigenvalue
2391  if (eigenValue > 0) {
2392  // Really want the root of the eigen value
2393  eigenValue = 3. * sqrt(eigenValue);
2394 
2395  // Recover the eigenvector
2396  const std::vector<double>& eigenVector = pca->getEigenVectors()[dimIdx];
2397 
2398  // Set the first point
2399  double xl = avePosition[0] - 0.5 * eigenValue * eigenVector[0];
2400  double yl = avePosition[1] - 0.5 * eigenValue * eigenVector[1];
2401  double zl = avePosition[2] - 0.5 * eigenValue * eigenVector[2];
2402 
2403  pl.SetPoint(0, xl, yl, zl);
2404  pmrk.SetPoint(pmrkIdx++, xl, yl, zl);
2405 
2406  // Set the second point
2407  double xu = avePosition[0] + 0.5 * eigenValue * eigenVector[0];
2408  double yu = avePosition[1] + 0.5 * eigenValue * eigenVector[1];
2409  double zu = avePosition[2] + 0.5 * eigenValue * eigenVector[2];
2410 
2411  pl.SetPoint(1, xu, yu, zu);
2412  pmrk.SetPoint(pmrkIdx++, xu, yu, zu);
2413  }
2414  }
2415 
2416  // By convention we will have drawn the "best" axis first
2417  if (recoOpt->fBestPCAAxisOnly) break;
2418 
2419  pcaIdx++;
2420  }
2421  }
2422  }
2423 
2424  // Now let's loop over daughters and call drawing routine for them
2425  if (pfPart->NumDaughters() > 0) {
2426  depth++;
2427 
2428  // std::string indent(depth, ' ');
2429 
2430  // std::cout << indent << "-drawPFParticle3D--> pfPart idx: " << pfPart->Self() << ", daughter list size: " << pfPart->Daughters().size() << std::endl;
2431 
2432  for (const auto& daughterIdx : pfPart->Daughters()) {
2433  DrawPFParticle3D(pfParticleVec.at(daughterIdx),
2434  pfParticleVec,
2435  spacePointVec,
2436  edgeAssnsVec,
2437  spacePointAssnVec,
2438  edgeSPAssnVec,
2439  spHitAssnVec,
2440  trackAssnVec,
2441  pcAxisAssnVec,
2442  cosmicTagAssnVec,
2443  depth,
2444  view);
2445  }
2446  }
2447 
2448  return;
2449  }
2450 
2451  //......................................................................
2452  void
2453  RecoBaseDrawer::Edge3D(const art::Event& evt, evdb::View3D* view)
2454  {
2455  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2456 
2457  if (recoOpt->fDrawEdges < 1) return;
2458 
2459  // The plan is to loop over the list of possible particles
2460  for (size_t imod = 0; imod < recoOpt->fEdgeLabels.size(); ++imod) {
2461  art::InputTag const which = recoOpt->fEdgeLabels[imod];
2462 
2463  // Start off by recovering our 3D Clusters for this label
2464  std::vector<art::Ptr<recob::Edge>> edgeVec;
2465  this->GetEdges(evt, which, edgeVec);
2466 
2467  mf::LogDebug("RecoBaseDrawer")
2468  << "RecoBaseDrawer: number Edges to draw: " << edgeVec.size() << std::endl;
2469 
2470  if (!edgeVec.empty()) {
2471  // Get the space points created by the PFParticle producer
2472  std::vector<art::Ptr<recob::SpacePoint>> spacePointVec;
2473  this->GetSpacePoints(evt, which, spacePointVec);
2474 
2475  // First draw the space points (all of them), then circle back on the edges...
2476  int colorIdx(41); //2);
2477 
2478  TPolyMarker3D& pm = view->AddPolyMarker3D(
2479  spacePointVec.size(), colorIdx, kFullDotMedium, 0.5); //kFullDotLarge, 0.5);
2480 
2481  for (const auto& spacePoint : spacePointVec) {
2482  TVector3 spPosition(spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]);
2483 
2484  pm.SetNextPoint(spPosition[0], spPosition[1], spPosition[2]);
2485  }
2486 
2487  // Now draw the edges
2488  for (const auto& edge : edgeVec) {
2489  art::Ptr<recob::SpacePoint> firstSP = spacePointVec.at(edge->FirstPointID());
2490  art::Ptr<recob::SpacePoint> secondSP = spacePointVec.at(edge->SecondPointID());
2491 
2492  if (firstSP->ID() != edge->FirstPointID() || secondSP->ID() != edge->SecondPointID()) {
2493  mf::LogDebug("RecoBaseDrawer")
2494  << "Edge: Space point index mismatch, first: " << firstSP->ID() << ", "
2495  << edge->FirstPointID() << ", second: " << secondSP->ID() << ", "
2496  << edge->SecondPointID() << std::endl;
2497  continue;
2498  }
2499 
2500  TVector3 startPoint(firstSP->XYZ()[0], firstSP->XYZ()[1], firstSP->XYZ()[2]);
2501  TVector3 endPoint(secondSP->XYZ()[0], secondSP->XYZ()[1], secondSP->XYZ()[2]);
2502  TVector3 lineVec(endPoint - startPoint);
2503 
2504  double length = lineVec.Mag();
2505 
2506  if (length == 0.) {
2507  // std::cout << "Edge length is zero, index 1: " << edge->FirstPointID() << ", index 2: " << edge->SecondPointID() << std::endl;
2508  continue;
2509  }
2510 
2511  // Get a polyline object to draw from the first to the second space point
2512  // TPolyLine3D& pl = view->AddPolyLine3D(2, colorIdx, 1, 1); //4, 1);
2513  //
2514  // pl.SetPoint(0, startPoint[0], startPoint[1], startPoint[2]);
2515  // pl.SetPoint(1, endPoint[0], endPoint[1], endPoint[2]);
2516  TPolyMarker3D& fakeLine = view->AddPolyMarker3D(10, 5, kFullDotMedium, 1.0);
2517 
2518  lineVec.SetMag(1.);
2519 
2520  for (int idx = 1; idx <= 10; idx++) {
2521  TVector3 plotPoint = startPoint + 0.1 * double(idx) * length * lineVec;
2522 
2523  fakeLine.SetNextPoint(plotPoint[0], plotPoint[1], plotPoint[2]);
2524  }
2525  }
2526  }
2527  }
2528 
2529  // Draw any associated Extreme Points
2530  for (size_t imod = 0; imod < recoOpt->fExtremePointLabels.size(); ++imod) {
2531  art::InputTag const which = recoOpt->fExtremePointLabels[imod];
2532 
2533  // Start off by recovering our 3D Clusters for this label
2534  std::vector<art::Ptr<recob::SpacePoint>> spacePointVec;
2535  this->GetSpacePoints(evt, which, spacePointVec);
2536 
2537  mf::LogDebug("RecoBaseDrawer")
2538  << "RecoBaseDrawer: number Extreme points to draw: " << spacePointVec.size() << std::endl;
2539 
2540  if (!spacePointVec.empty()) {
2541  // First draw the space points (all of them), then circle back on the edges...
2542  int colorIdx(kYellow);
2543 
2544  TPolyMarker3D& pm = view->AddPolyMarker3D(
2545  spacePointVec.size(), colorIdx, kFullDotLarge, 1.0); //kFullDotLarge, 0.5);
2546 
2547  for (const auto& spacePoint : spacePointVec) {
2548  TVector3 spPosition(spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]);
2549 
2550  pm.SetNextPoint(spPosition[0], spPosition[1], spPosition[2]);
2551  }
2552  }
2553  }
2554 
2555  return;
2556  }
2557 
2558  //......................................................................
2559  void
2560  RecoBaseDrawer::Prong3D(const art::Event& evt, evdb::View3D* view)
2561  {
2562  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2563  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2564 
2565  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
2566 
2567  // annoying for now, but have to have multiple copies of basically the
2568  // same code to draw prongs, showers and tracks so that we can use
2569  // the art::Assns to get the hits and clusters.
2570 
2571  // Tracks.
2572 
2573  if (recoOpt->fDrawTracks > 2) {
2574  for (size_t imod = 0; imod < recoOpt->fTrackLabels.size(); ++imod) {
2575  art::InputTag which = recoOpt->fTrackLabels[imod];
2576  art::View<recob::Track> trackView;
2577  this->GetTracks(evt, which, trackView);
2578  if (!trackView.isValid())
2579  continue; //Prevent potential segmentation fault if no tracks found. aoliv23@lsu.edu
2580 
2581  art::PtrVector<recob::Track> trackVec;
2582 
2583  trackView.fill(trackVec);
2584 
2585  art::InputTag const cosmicTagLabel(
2586  recoOpt->fCosmicTagLabels.size() > imod ? recoOpt->fCosmicTagLabels[imod] : "");
2587  art::FindMany<anab::CosmicTag> cosmicTagAssnVec(trackVec, evt, cosmicTagLabel);
2588 
2589  for (const auto& track : trackVec) {
2590  int color = evd::kColor[track.key() % evd::kNCOLS];
2591  int marker = kFullDotLarge;
2592  float size = 2.0;
2593 
2594  // Check if a CosmicTag object is available
2595 
2596  // Recover cosmic tag info if any
2597  if (cosmicTagAssnVec.isValid()) {
2598  std::vector<const anab::CosmicTag*> tkCosmicTagVec = cosmicTagAssnVec.at(track.key());
2599 
2600  if (!tkCosmicTagVec.empty()) {
2601  const anab::CosmicTag* cosmicTag = tkCosmicTagVec.front();
2602 
2603  // If tagged as Cosmic then neutralize the color
2604  if (cosmicTag->CosmicScore() > 0.6) {
2605  color = 14;
2606  size = 0.5;
2607  }
2608  }
2609  }
2610 
2611  // Draw track using only embedded information.
2612 
2613  DrawTrack3D(*track, view, color, marker, size);
2614  }
2615  }
2616  }
2617 
2618  // Showers.
2619 
2620  if (recoOpt->fDrawShowers != 0) {
2621  for (size_t imod = 0; imod < recoOpt->fShowerLabels.size(); ++imod) {
2622  art::InputTag which = recoOpt->fShowerLabels[imod];
2623  art::View<recob::Shower> shower;
2624  this->GetShowers(evt, which, shower);
2625 
2626  for (size_t s = 0; s < shower.vals().size(); ++s) {
2627  const recob::Shower* pshower = shower.vals().at(s);
2628  int color = pshower->ID();
2629  DrawShower3D(*pshower, color, view);
2630  }
2631  }
2632  }
2633 
2634  return;
2635  }
2636 
2637  //......................................................................
2638  void
2639  RecoBaseDrawer::DrawTrack3D(const recob::Track& track,
2640  evdb::View3D* view,
2641  int color,
2642  int marker,
2643  float size)
2644  {
2645  // Get options.
2646  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2647 
2648  if (recoOpt->fDrawTrackSpacePoints) {
2649  // Use brute force to find the module label and index of this
2650  // track, so that we can find associated space points and draw
2651  // them.
2652  const art::Event* evt = evdb::EventHolder::Instance()->GetEvent();
2653  //std::vector<art::Handle<std::vector<recob::Track>>> handles;
2654  //evt->getManyByType(handles);
2655  auto handles = evt->getMany<std::vector<recob::Track>>();
2656 
2657  for (auto ih : handles) {
2658  const art::Handle<std::vector<recob::Track>> handle = ih;
2659 
2660  if (handle.isValid()) {
2661  const std::string& which = handle.provenance()->moduleLabel();
2662  art::FindManyP<recob::SpacePoint> fmsp(handle, *evt, which);
2663 
2664  if (fmsp.isValid() && fmsp.size() > 0) {
2665  int n = handle->size();
2666  float spSize = 0.5 * size;
2667 
2668  for (int i = 0; i < n; ++i) {
2669  art::Ptr<recob::Track> p(handle, i);
2670  if (&*p == &track) {
2671  std::vector<art::Ptr<recob::SpacePoint>> spts = fmsp.at(i);
2672  fSpacePointDrawer->Draw(spts, view, color, marker, spSize);
2673  }
2674  }
2675  }
2676  }
2677  }
2678  }
2679 
2680  if (recoOpt->fDrawTrackTrajectoryPoints) {
2681  // Draw trajectory points.
2682  int np = track.NumberTrajectoryPoints();
2683 
2684  int lineSize = size;
2685 
2686  if (lineSize < 1) lineSize = 1;
2687 
2688  // Make and fill a special polymarker for the head of the track
2689  //TPolyMarker3D& pmStart = view->AddPolyMarker3D(1, color, 4, 3);
2690  TPolyMarker3D& pmStart = view->AddPolyMarker3D(1, 0, marker, 2. * size);
2691 
2692  const auto& firstPos = track.LocationAtPoint(0);
2693  pmStart.SetPoint(0, firstPos.X(), firstPos.Y(), firstPos.Z());
2694 
2695  // Make and fill a polymarker.
2696  TPolyMarker3D& pm = view->AddPolyMarker3D(track.CountValidPoints(), color, 1, 3);
2697 
2698  for (int p = 0; p < np; ++p) {
2699  if (!track.HasValidPoint(p)) continue;
2700 
2701  const auto& pos = track.LocationAtPoint(p);
2702  pm.SetPoint(p, pos.X(), pos.Y(), pos.Z());
2703  }
2704 
2705  // As we are a track, should we not be drawing a line here?
2706  TPolyLine3D& pl = view->AddPolyLine3D(track.CountValidPoints(), color, lineSize, 7);
2707 
2708  for (int p = 0; p < np; ++p) {
2709  if (!track.HasValidPoint(p)) continue;
2710 
2711  const auto pos = track.LocationAtPoint(p);
2712 
2713  pl.SetPoint(p, pos.X(), pos.Y(), pos.Z());
2714  }
2715 
2716  if (recoOpt->fDrawTrackTrajectoryPoints > 4) {
2717  // Can we add the track direction at each point?
2718  // This won't work for the last point... but let's try
2719  auto startPos = track.LocationAtPoint(0);
2720  auto startDir = track.DirectionAtPoint(0);
2721 
2722  for (int p = 1; p < np; ++p) {
2723  if (!track.HasValidPoint(p)) continue;
2724 
2725  TPolyLine3D& pl = view->AddPolyLine3D(2, (color + 1) % evd::kNCOLS, size, 7); //1, 3);
2726 
2727  auto nextPos = track.LocationAtPoint(p);
2728  auto deltaPos = nextPos - startPos;
2729  double arcLen = deltaPos.Dot(
2730  startDir); // arc len to plane containing next point perpendicular to current point dir
2731 
2732  if (arcLen < 0.) arcLen = 3.;
2733 
2734  auto endPoint = startPos + arcLen * startDir;
2735 
2736  pl.SetPoint(0, startPos.X(), startPos.Y(), startPos.Z());
2737  pl.SetPoint(1, endPoint.X(), endPoint.Y(), endPoint.Z());
2738 
2739  startPos = nextPos;
2740  startDir = track.DirectionAtPoint(p);
2741  }
2742  }
2743  }
2744 
2745  return;
2746  }
2747 
2748  //......................................................................
2749  void
2750  RecoBaseDrawer::DrawShower3D(const recob::Shower& shower, int color, evdb::View3D* view)
2751  {
2752  // Use brute force to find the module label and index of this
2753  // shower, so that we can find associated space points and draw
2754  // them.
2755  // B. Baller: Catch an exception if there are no space points and draw a cone instead.
2756 
2757  const art::Event* evt = evdb::EventHolder::Instance()->GetEvent();
2758  //std::vector<art::Handle<std::vector<recob::Shower>>> handles;
2759  //evt->getManyByType(handles);
2760  auto handles = evt->getMany<std::vector<recob::Shower>>();
2761 
2762  bool noSpts = false;
2763 
2764  for (auto ih : handles) {
2765  const art::Handle<std::vector<recob::Shower>> handle = ih;
2766 
2767  if (handle.isValid()) {
2768 
2769  const std::string& which = handle.provenance()->moduleLabel();
2770  art::FindManyP<recob::SpacePoint> fmsp(handle, *evt, which);
2771 
2772  int n = handle->size();
2773  for (int i = 0; i < n; ++i) {
2774  art::Ptr<recob::Shower> p(handle, i);
2775  if (&*p == &shower) {
2776  // BB catch if no space points
2777  std::vector<art::Ptr<recob::SpacePoint>> spts;
2778  try {
2779  spts = fmsp.at(i);
2780  fSpacePointDrawer->Draw(spts, view, color);
2781  }
2782  catch (...) {
2783  noSpts = true;
2784  continue;
2785  } // catch
2786  } // shower
2787  } // i
2788  } // ih
2789  }
2790 
2791  if (noSpts && shower.has_length() && shower.has_open_angle()) {
2792  std::cout << "No space points associated with the shower. Drawing a cone instead\n";
2793  color = kRed;
2794  auto& dedx = shower.dEdx();
2795  if (!dedx.empty()) {
2796  double dedxAve = 0;
2797  for (auto& dedxInPln : dedx)
2798  dedxAve += dedxInPln;
2799  dedxAve /= (double)dedx.size();
2800  // Use blue for ~1 MIP
2801  color = kBlue;
2802  // use green for ~2 MIP
2803  if (dedxAve > 3 && dedxAve < 5) color = kGreen;
2804  }
2805  double radius = shower.Length() * shower.OpenAngle();
2806  TVector3 startPos = shower.ShowerStart();
2807  TVector3 endPos = startPos + shower.Length() * shower.Direction();
2808  auto coneRim = Circle3D(endPos, shower.Direction(), radius);
2809  TPolyLine3D& pl = view->AddPolyLine3D(coneRim.size(), color, 2, 0);
2810  for (unsigned short ipt = 0; ipt < coneRim.size(); ++ipt) {
2811  auto& pt = coneRim[ipt];
2812  pl.SetPoint(ipt, pt[0], pt[1], pt[2]);
2813  }
2814  // Draw a line from the start position to each point on the rim
2815  for (unsigned short ipt = 0; ipt < coneRim.size(); ++ipt) {
2816  TPolyLine3D& panel = view->AddPolyLine3D(2, color, 2, 0);
2817  panel.SetPoint(0, startPos.X(), startPos.Y(), startPos.Z());
2818  panel.SetPoint(1, coneRim[ipt][0], coneRim[ipt][1], coneRim[ipt][2]);
2819  } // ipt
2820 
2821  } // no space points
2822 
2823  return;
2824  }
2825 
2826  //......................................................................
2827  std::vector<std::array<double, 3>>
2828  RecoBaseDrawer::Circle3D(const TVector3& centerPos, const TVector3& axisDir, const double& radius)
2829  {
2830  // B. Baller Create a polyline circle in 3D
2831 
2832  // Make the rotation matrix to transform into the circle coordinate system
2833  TRotation r;
2834  r.RotateX(axisDir.X());
2835  r.RotateY(axisDir.Y());
2836  r.RotateZ(axisDir.Z());
2837  constexpr unsigned short nRimPts = 16;
2838  std::vector<std::array<double, 3>> rimPts(nRimPts + 1);
2839  for (unsigned short iang = 0; iang < nRimPts; ++iang) {
2840  double rimAngle = iang * 2 * M_PI / (float)nRimPts;
2841  TVector3 rim = {0, 0, 1};
2842  rim.SetX(radius * cos(rimAngle));
2843  rim.SetY(radius * sin(rimAngle));
2844  rim.SetZ(0);
2845  rim.Transform(r);
2846  rim += centerPos;
2847  for (unsigned short ixyz = 0; ixyz < 3; ++ixyz)
2848  rimPts[iang][ixyz] = rim[ixyz];
2849  } // iang
2850  // close the circle
2851  rimPts[nRimPts] = rimPts[0];
2852  return rimPts;
2853  } // PolyLineCircle
2854 
2855  //......................................................................
2856  void
2857  RecoBaseDrawer::Vertex3D(const art::Event& evt, evdb::View3D* view)
2858  {
2859  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2860  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2861 
2862  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
2863 
2864  if (recoOpt->fDrawVertices != 0) {
2865 
2866  for (size_t imod = 0; imod < recoOpt->fVertexLabels.size(); ++imod) {
2867  art::InputTag const which = recoOpt->fVertexLabels[imod];
2868 
2869  art::PtrVector<recob::Vertex> vertex;
2870  this->GetVertices(evt, which, vertex);
2871 
2872  art::FindManyP<recob::Track> fmt(vertex, evt, which);
2873  art::FindManyP<recob::Shower> fms(vertex, evt, which);
2874 
2875  for (size_t v = 0; v < vertex.size(); ++v) {
2876 
2877  if (fmt.isValid()) {
2878  std::vector<art::Ptr<recob::Track>> tracks = fmt.at(v);
2879 
2880  // grab the Prongs from the vertex and draw those
2881  for (size_t t = 0; t < tracks.size(); ++t)
2882  this->DrawTrack3D(*(tracks[t]), view, vertex[v]->ID());
2883  }
2884 
2885  if (fms.isValid()) {
2886  std::vector<art::Ptr<recob::Shower>> showers = fms.at(v);
2887 
2888  for (size_t s = 0; s < showers.size(); ++s)
2889  this->DrawShower3D(*(showers[s]), vertex[v]->ID(), view);
2890  }
2891 
2892  double xyz[3] = {0.};
2893  vertex[v]->XYZ(xyz);
2894  TPolyMarker3D& pm =
2895  view->AddPolyMarker3D(1, evd::kColor[vertex[v]->ID() % evd::kNCOLS], 29, 6);
2896  pm.SetPoint(0, xyz[0], xyz[1], xyz[2]);
2897 
2898  } // end loop over vertices to draw from this label
2899  } // end loop over vertex module lables
2900  } // end if we are drawing vertices
2901 
2902  return;
2903  }
2904 
2905  //......................................................................
2906  void
2907  RecoBaseDrawer::Event3D(const art::Event& evt, evdb::View3D* view)
2908  {
2909  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2910  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2911 
2912  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
2913  if (recoOpt->fDrawEvents != 0) {
2914 
2915  for (size_t imod = 0; imod < recoOpt->fEventLabels.size(); ++imod) {
2916  art::InputTag const which = recoOpt->fEventLabels[imod];
2917 
2918  art::PtrVector<recob::Event> event;
2919  this->GetEvents(evt, which, event);
2920 
2921  if (event.size() < 1) continue;
2922 
2923  art::FindManyP<recob::Vertex> fmvp(event, evt, which);
2924  art::FindMany<recob::Vertex> fmv(event, evt, which);
2925 
2926  for (size_t e = 0; e < event.size(); ++e) {
2927 
2928  // grab the vertices for this event
2929  std::vector<art::Ptr<recob::Vertex>> vertex = fmvp.at(e);
2930 
2931  if (vertex.size() < 1) continue;
2932 
2933  art::FindManyP<recob::Track> fmt(vertex, evt, recoOpt->fVertexLabels[0]);
2934  art::FindManyP<recob::Shower> fms(vertex, evt, recoOpt->fVertexLabels[0]);
2935 
2936  for (size_t v = 0; v < vertex.size(); ++v) {
2937 
2938  /// \todo need a better way to grab the vertex module labels,
2939  // right now assume there is only 1 in the list
2940  std::vector<art::Ptr<recob::Track>> tracks = fmt.at(v);
2941  std::vector<art::Ptr<recob::Shower>> showers = fms.at(v);
2942 
2943  // grab the Prongs from the vertex and draw those
2944  for (size_t t = 0; t < tracks.size(); ++t)
2945  this->DrawTrack3D(*(tracks[t]), view, event[e]->ID());
2946 
2947  for (size_t s = 0; s < showers.size(); ++s)
2948  this->DrawShower3D(*(showers[s]), event[e]->ID(), view);
2949 
2950  } // end loop over vertices from this event
2951 
2952  double xyz[3] = {0.};
2953  std::vector<const recob::Vertex*> vts = fmv.at(e);
2954 
2955  event[e]->PrimaryVertex(vts)->XYZ(xyz);
2956  TPolyMarker3D& pm =
2957  view->AddPolyMarker3D(1, evd::kColor[event[e]->ID() % evd::kNCOLS], 29, 6);
2958  pm.SetPoint(0, xyz[0], xyz[1], xyz[2]);
2959 
2960  } // end loop over events
2961  } // end loop over event module lables
2962  } // end if we are drawing events
2963 
2964  return;
2965  }
2966  //......................................................................
2967  void
2968  RecoBaseDrawer::Slice3D(const art::Event& evt, evdb::View3D* view)
2969  {
2970  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2971  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
2972 
2973  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
2974  if (recoOpt->fDrawSlices < 1) return;
2975  if (recoOpt->fDrawSliceSpacePoints < 1) return;
2976  for (size_t imod = 0; imod < recoOpt->fSliceLabels.size(); ++imod) {
2977  art::InputTag const which = recoOpt->fSliceLabels[imod];
2978  art::PtrVector<recob::Slice> slices;
2979  this->GetSlices(evt, which, slices);
2980  if (slices.size() < 1) continue;
2981  art::FindManyP<recob::SpacePoint> fmsp(slices, evt, which);
2982  for (size_t isl = 0; isl < slices.size(); ++isl) {
2983  int slcID = std::abs(slices[isl]->ID());
2984  int color = evd::kColor[slcID % evd::kNCOLS];
2985  std::vector<art::Ptr<recob::SpacePoint>> spts = fmsp.at(isl);
2986  fSpacePointDrawer->Draw(spts, view, color, kFullDotLarge, 2);
2987  }
2988  }
2989  }
2990  //......................................................................
2991  void
2992  RecoBaseDrawer::OpFlashOrtho(const art::Event& evt,
2993  detinfo::DetectorClocksData const& clockData,
2994  detinfo::DetectorPropertiesData const& detProp,
2995  evd::OrthoProj_t proj,
2996  evdb::View2D* view)
2997  {
2998  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
2999  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3000  art::ServiceHandle<geo::Geometry const> geo;
3001 
3002  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
3003  if (recoOpt->fDrawOpFlashes == 0) return;
3004 
3005  geo::PlaneID pid(rawOpt->fCryostat, rawOpt->fTPC, 0);
3006 
3007  double minx = 1e9;
3008  double maxx = -1e9;
3009  for (size_t i = 0; i < geo->NTPC(); ++i) {
3010  double local[3] = {0., 0., 0.};
3011  double world[3] = {0., 0., 0.};
3012  const geo::TPCGeo& tpc = geo->TPC(i);
3013  tpc.LocalToWorld(local, world);
3014  if (minx > world[0] - geo->DetHalfWidth(i)) minx = world[0] - geo->DetHalfWidth(i);
3015  if (maxx < world[0] + geo->DetHalfWidth(i)) maxx = world[0] + geo->DetHalfWidth(i);
3016  }
3017 
3018  for (size_t imod = 0; imod < recoOpt->fOpFlashLabels.size(); ++imod) {
3019  const art::InputTag which = recoOpt->fOpFlashLabels[imod];
3020 
3021  art::PtrVector<recob::OpFlash> opflashes;
3022  this->GetOpFlashes(evt, which, opflashes);
3023 
3024  if (opflashes.size() < 1) continue;
3025 
3026  int NFlashes = opflashes.size();
3027 
3028  // project each seed into this view
3029  for (int iof = 0; iof < NFlashes; ++iof) {
3030 
3031  if (opflashes[iof]->TotalPE() < recoOpt->fFlashMinPE) continue;
3032  if (opflashes[iof]->Time() < recoOpt->fFlashTMin) continue;
3033  if (opflashes[iof]->Time() > recoOpt->fFlashTMax) continue;
3034 
3035  double YCentre = opflashes[iof]->YCenter();
3036  double YHalfWidth = opflashes[iof]->YWidth();
3037  double ZCentre = opflashes[iof]->ZCenter();
3038  double ZHalfWidth = opflashes[iof]->ZWidth();
3039 
3040  int Colour = evd::kColor[(iof) % evd::kNCOLS];
3041 
3042  if (proj == evd::kXY) {
3043  TBox& b1 = view->AddBox(YCentre - YHalfWidth, minx, YCentre + YHalfWidth, maxx);
3044  b1.SetFillStyle(3004 + (iof % 3));
3045  b1.SetFillColor(Colour);
3046  }
3047  else if (proj == evd::kXZ) {
3048  float xflash = detProp.ConvertTicksToX(
3049  opflashes[iof]->Time() / sampling_rate(clockData) * 1e3 + detProp.GetXTicksOffset(pid),
3050  pid);
3051  TLine& line = view->AddLine(ZCentre - ZHalfWidth, xflash, ZCentre + ZHalfWidth, xflash);
3052  line.SetLineWidth(2);
3053  line.SetLineStyle(2);
3054  line.SetLineColor(Colour);
3055  }
3056  else if (proj == evd::kYZ) {
3057  TBox& b1 = view->AddBox(
3058  ZCentre - ZHalfWidth, YCentre - YHalfWidth, ZCentre + ZHalfWidth, YCentre + YHalfWidth);
3059  b1.SetFillStyle(3004 + (iof % 3));
3060  b1.SetFillColor(Colour);
3061  view->AddMarker(ZCentre, YCentre, Colour, 4, 1.5);
3062  }
3063 
3064  } // Flashes with this label
3065  } // Vector of OpFlash labels
3066  }
3067  //......................................................................
3068  void
3069  RecoBaseDrawer::VertexOrtho(const art::PtrVector<recob::Vertex>& vertex,
3070  evd::OrthoProj_t proj,
3071  evdb::View2D* view,
3072  int marker)
3073  {
3074  for (size_t v = 0; v < vertex.size(); ++v) {
3075 
3076  double xyz[3] = {0.};
3077  vertex[v]->XYZ(xyz);
3078 
3079  int color = evd::kColor[vertex[v]->ID() % evd::kNCOLS];
3080 
3081  if (proj == evd::kXY) {
3082  TMarker& strt = view->AddMarker(xyz[1], xyz[0], color, marker, 1.0);
3083  strt.SetMarkerColor(color);
3084  }
3085  else if (proj == evd::kXZ) {
3086  TMarker& strt = view->AddMarker(xyz[2], xyz[0], color, marker, 1.0);
3087  strt.SetMarkerColor(color);
3088  }
3089  else if (proj == evd::kYZ) {
3090  TMarker& strt = view->AddMarker(xyz[2], xyz[1], color, marker, 1.0);
3091  strt.SetMarkerColor(color);
3092  }
3093  }
3094  }
3095  void
3096  RecoBaseDrawer::VertexOrtho(const art::Event& evt, evd::OrthoProj_t proj, evdb::View2D* view)
3097  {
3098  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
3099  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3100  art::ServiceHandle<geo::Geometry const> geo;
3101 
3102  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
3103  if (recoOpt->fDrawVertices == 0) return;
3104 
3105  for (size_t imod = 0; imod < recoOpt->fVertexLabels.size(); ++imod) {
3106  art::InputTag const which = recoOpt->fVertexLabels[imod];
3107 
3108  art::PtrVector<recob::Vertex> vertex;
3109  this->GetVertices(evt, which, vertex);
3110  this->VertexOrtho(vertex, proj, view, 24);
3111 
3112  //this->GetVertices(evt, art::InputTag(which.label(), "kink", which.process()), vertex);
3113  //this->VertexOrtho(vertex, proj, view, 27);
3114 
3115  //this->GetVertices(evt, art::InputTag(which.label(), "node", which.process()), vertex);
3116  //this->VertexOrtho(vertex, proj, view, 22);
3117  }
3118  return;
3119  }
3120 
3121  //......................................................................
3122  void
3123  RecoBaseDrawer::SpacePointOrtho(const art::Event& evt,
3124  evd::OrthoProj_t proj,
3125  double msize,
3126  evdb::View2D* view)
3127  {
3128  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
3129  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3130 
3131  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
3132 
3133  std::vector<art::InputTag> labels;
3134  if (recoOpt->fDrawSpacePoints != 0) {
3135  for (size_t imod = 0; imod < recoOpt->fSpacePointLabels.size(); ++imod)
3136  labels.push_back(recoOpt->fSpacePointLabels[imod]);
3137  }
3138 
3139  for (size_t imod = 0; imod < labels.size(); ++imod) {
3140  art::InputTag const which = labels[imod];
3141 
3142  std::vector<art::Ptr<recob::SpacePoint>> spts;
3143  this->GetSpacePoints(evt, which, spts);
3144  int color = imod;
3145 
3146  this->DrawSpacePointOrtho(spts, color, proj, msize, view);
3147  }
3148 
3149  return;
3150  }
3151 
3152  //......................................................................
3153  void
3154  RecoBaseDrawer::PFParticleOrtho(const art::Event& evt,
3155  evd::OrthoProj_t proj,
3156  double msize,
3157  evdb::View2D* view)
3158  {
3159  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
3160  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3161 
3162  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
3163  if (recoOpt->fDrawPFParticles < 1) return;
3164 
3165  // The plan is to loop over the list of possible particles
3166  for (size_t imod = 0; imod < recoOpt->fPFParticleLabels.size(); ++imod) {
3167  art::InputTag const which = recoOpt->fPFParticleLabels[imod];
3168 
3169  // Start off by recovering our 3D Clusters for this label
3170  art::PtrVector<recob::PFParticle> pfParticleVec;
3171  this->GetPFParticles(evt, which, pfParticleVec);
3172 
3173  // Make sure we have some clusters
3174  if (pfParticleVec.size() < 1) continue;
3175 
3176  // Add the relations to recover associations cluster hits
3177  art::FindMany<recob::SpacePoint> spacePointAssnVec(pfParticleVec, evt, which);
3178 
3179  // If no valid space point associations then nothing to do
3180  if (!spacePointAssnVec.isValid()) continue;
3181 
3182  // Need the PCA info as well
3183  art::FindMany<recob::PCAxis> pcAxisAssnVec(pfParticleVec, evt, which);
3184 
3185  if (!pcAxisAssnVec.isValid()) continue;
3186 
3187  // Commence looping over possible clusters
3188  for (size_t idx = 0; idx < pfParticleVec.size(); idx++) {
3189  // Recover cluster
3190  const art::Ptr<recob::PFParticle> pfParticle(pfParticleVec.at(idx));
3191 
3192  // Drawing will be done recursively down the chain of hieirarchy... so we want to begin
3193  // with only "primary" particles, if we find one that isn't then we skip
3194  if (!pfParticle->IsPrimary()) continue;
3195 
3196  // Call the recursive drawing routine
3197  DrawPFParticleOrtho(
3198  pfParticle, pfParticleVec, spacePointAssnVec, pcAxisAssnVec, 0, proj, view);
3199  }
3200  }
3201 
3202  return;
3203  }
3204 
3205  void
3206  RecoBaseDrawer::DrawPFParticleOrtho(const art::Ptr<recob::PFParticle>& pfPart,
3207  const art::PtrVector<recob::PFParticle>& pfParticleVec,
3208  const art::FindMany<recob::SpacePoint>& spacePointAssnVec,
3209  const art::FindMany<recob::PCAxis>& pcAxisAssnVec,
3210  int depth,
3211  evd::OrthoProj_t proj,
3212  evdb::View2D* view)
3213  {
3214  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3215 
3216  // First let's draw the hits associated to this cluster
3217  const std::vector<const recob::SpacePoint*>& hitsVec(spacePointAssnVec.at(pfPart->Self()));
3218 
3219  // Use the particle ID to determine the color to draw the points
3220  // Ok, this is what we would like to do eventually but currently all particles are the same...
3221  // int colorIdx = evd::Style::ColorFromPDG(pfPart->PdgCode());
3222  int colorIdx = evd::kColor[pfPart->Self() % evd::kNCOLS];
3223 
3224  if (!hitsVec.empty()) {
3225  std::vector<const recob::SpacePoint*> hitPosVec;
3226  std::vector<const recob::SpacePoint*> skeletonPosVec;
3227  std::vector<const recob::SpacePoint*> skelEdgePosVec;
3228  std::vector<const recob::SpacePoint*> edgePosVec;
3229  std::vector<const recob::SpacePoint*> seedPosVec;
3230  std::vector<const recob::SpacePoint*> pairPosVec;
3231 
3232  for (const auto& spacePoint : hitsVec) {
3233  if (spacePoint->Chisq() > 0.)
3234  hitPosVec.push_back(spacePoint);
3235  else if (spacePoint->Chisq() == -1.)
3236  skeletonPosVec.push_back(spacePoint);
3237  else if (spacePoint->Chisq() == -3.)
3238  skelEdgePosVec.push_back(spacePoint);
3239  else if (spacePoint->Chisq() == -4.)
3240  seedPosVec.push_back(spacePoint);
3241  else if (spacePoint->Chisq() > -10.)
3242  edgePosVec.push_back(spacePoint);
3243  else
3244  pairPosVec.push_back(spacePoint);
3245  }
3246 
3247  int hitIdx(0);
3248 
3249  if (!recoOpt->fSkeletonOnly) {
3250  TPolyMarker& pm1 = view->AddPolyMarker(hitPosVec.size(), colorIdx, kFullDotMedium, 1);
3251  for (const auto* spacePoint : hitPosVec) {
3252  const double* pos = spacePoint->XYZ();
3253 
3254  if (proj == evd::kXY)
3255  pm1.SetPoint(hitIdx++, pos[0], pos[1]);
3256  else if (proj == evd::kXZ)
3257  pm1.SetPoint(hitIdx++, pos[2], pos[0]);
3258  else if (proj == evd::kYZ)
3259  pm1.SetPoint(hitIdx++, pos[2], pos[1]);
3260  }
3261 
3262  hitIdx = 0;
3263 
3264  TPolyMarker& pm2 = view->AddPolyMarker(edgePosVec.size(), 28, kFullDotMedium, 1);
3265  for (const auto* spacePoint : edgePosVec) {
3266  const double* pos = spacePoint->XYZ();
3267 
3268  if (proj == evd::kXY)
3269  pm2.SetPoint(hitIdx++, pos[0], pos[1]);
3270  else if (proj == evd::kXZ)
3271  pm2.SetPoint(hitIdx++, pos[2], pos[0]);
3272  else if (proj == evd::kYZ)
3273  pm2.SetPoint(hitIdx++, pos[2], pos[1]);
3274  }
3275 
3276  hitIdx = 0;
3277 
3278  TPolyMarker& pm3 = view->AddPolyMarker(pairPosVec.size(), 2, kFullDotMedium, 1);
3279  for (const auto* spacePoint : pairPosVec) {
3280  const double* pos = spacePoint->XYZ();
3281 
3282  if (proj == evd::kXY)
3283  pm3.SetPoint(hitIdx++, pos[0], pos[1]);
3284  else if (proj == evd::kXZ)
3285  pm3.SetPoint(hitIdx++, pos[2], pos[0]);
3286  else if (proj == evd::kYZ)
3287  pm3.SetPoint(hitIdx++, pos[2], pos[1]);
3288  }
3289  }
3290 
3291  hitIdx = 0;
3292 
3293  TPolyMarker& pm4 = view->AddPolyMarker(skeletonPosVec.size(), 1, kFullDotMedium, 1);
3294  for (const auto* spacePoint : skeletonPosVec) {
3295  const double* pos = spacePoint->XYZ();
3296 
3297  if (proj == evd::kXY)
3298  pm4.SetPoint(hitIdx++, pos[0], pos[1]);
3299  else if (proj == evd::kXZ)
3300  pm4.SetPoint(hitIdx++, pos[2], pos[0]);
3301  else if (proj == evd::kYZ)
3302  pm4.SetPoint(hitIdx++, pos[2], pos[1]);
3303  }
3304 
3305  hitIdx = 0;
3306 
3307  TPolyMarker& pm5 = view->AddPolyMarker(skelEdgePosVec.size(), 3, kFullDotMedium, 1);
3308  for (const auto* spacePoint : skelEdgePosVec) {
3309  const double* pos = spacePoint->XYZ();
3310 
3311  if (proj == evd::kXY)
3312  pm5.SetPoint(hitIdx++, pos[0], pos[1]);
3313  else if (proj == evd::kXZ)
3314  pm5.SetPoint(hitIdx++, pos[2], pos[0]);
3315  else if (proj == evd::kYZ)
3316  pm5.SetPoint(hitIdx++, pos[2], pos[1]);
3317  }
3318 
3319  hitIdx = 0;
3320 
3321  TPolyMarker& pm6 = view->AddPolyMarker(seedPosVec.size(), 6, kFullDotMedium, 1);
3322  for (const auto* spacePoint : seedPosVec) {
3323  const double* pos = spacePoint->XYZ();
3324 
3325  if (proj == evd::kXY)
3326  pm6.SetPoint(hitIdx++, pos[0], pos[1]);
3327  else if (proj == evd::kXZ)
3328  pm6.SetPoint(hitIdx++, pos[2], pos[0]);
3329  else if (proj == evd::kYZ)
3330  pm6.SetPoint(hitIdx++, pos[2], pos[1]);
3331  }
3332  }
3333 
3334  // Look up the PCA info
3335  if (pcAxisAssnVec.isValid()) {
3336  std::vector<const recob::PCAxis*> pcaVec(pcAxisAssnVec.at(pfPart->Self()));
3337 
3338  if (!pcaVec.empty()) {
3339  // For each axis we are going to draw a solid line between two points
3340  int numPoints(2);
3341  int lineWidth[2] = {3, 1};
3342  int lineStyle[2] = {1, 13};
3343  int lineColor[2] = {colorIdx, 18};
3344  int markStyle[2] = {4, 4};
3345  int pcaIdx(0);
3346 
3347  // The order of axes in the returned association vector is arbitrary... the "first" axis is
3348  // better and we can divine that by looking at the axis id's (the best will have been made first)
3349  if (pcaVec.size() > 1 && pcaVec.front()->getID() > pcaVec.back()->getID())
3350  std::reverse(pcaVec.begin(), pcaVec.end());
3351 
3352  for (const auto& pca : pcaVec) {
3353  // We need the mean position
3354  const double* avePosition = pca->getAvePosition();
3355 
3356  // Let's draw a marker at the interesting points
3357  int pmrkIdx(0);
3358  TPolyMarker& pmrk = view->AddPolyMarker(7, lineColor[pcaIdx], markStyle[pcaIdx], 1);
3359 
3360  if (proj == evd::kXY)
3361  pmrk.SetPoint(pmrkIdx++, avePosition[0], avePosition[1]);
3362  else if (proj == evd::kXZ)
3363  pmrk.SetPoint(pmrkIdx++, avePosition[2], avePosition[0]);
3364  else if (proj == evd::kYZ)
3365  pmrk.SetPoint(pmrkIdx++, avePosition[2], avePosition[1]);
3366 
3367  // Loop over pca dimensions
3368  for (int dimIdx = 0; dimIdx < 3; dimIdx++) {
3369  // Oh please oh please give me an instance of a poly line...
3370  TPolyLine& pl =
3371  view->AddPolyLine(numPoints, lineColor[pcaIdx], lineWidth[pcaIdx], lineStyle[pcaIdx]);
3372 
3373  // We will use the eigen value to give the length of the line we're going to plot
3374  double eigenValue = pca->getEigenValues()[dimIdx];
3375 
3376  // Make sure a valid eigenvalue
3377  if (eigenValue > 0) {
3378  // Really want the root of the eigen value
3379  eigenValue = 3. * sqrt(eigenValue);
3380 
3381  // Recover the eigenvector
3382  const std::vector<double>& eigenVector = pca->getEigenVectors()[dimIdx];
3383 
3384  // Set the first point
3385  double xl = avePosition[0] - 0.5 * eigenValue * eigenVector[0];
3386  double yl = avePosition[1] - 0.5 * eigenValue * eigenVector[1];
3387  double zl = avePosition[2] - 0.5 * eigenValue * eigenVector[2];
3388 
3389  if (proj == evd::kXY) {
3390  pl.SetPoint(0, xl, yl);
3391  pmrk.SetPoint(pmrkIdx++, xl, yl);
3392  }
3393  else if (proj == evd::kXZ) {
3394  pl.SetPoint(0, zl, xl);
3395  pmrk.SetPoint(pmrkIdx++, zl, xl);
3396  }
3397  else if (proj == evd::kYZ) {
3398  pl.SetPoint(0, zl, yl);
3399  pmrk.SetPoint(pmrkIdx++, zl, yl);
3400  }
3401 
3402  // Set the second point
3403  double xu = avePosition[0] + 0.5 * eigenValue * eigenVector[0];
3404  double yu = avePosition[1] + 0.5 * eigenValue * eigenVector[1];
3405  double zu = avePosition[2] + 0.5 * eigenValue * eigenVector[2];
3406 
3407  if (proj == evd::kXY) {
3408  pl.SetPoint(1, xu, yu);
3409  pmrk.SetPoint(pmrkIdx++, xu, yu);
3410  }
3411  else if (proj == evd::kXZ) {
3412  pl.SetPoint(1, zu, xu);
3413  pmrk.SetPoint(pmrkIdx++, zu, xu);
3414  }
3415  else if (proj == evd::kYZ) {
3416  pl.SetPoint(1, zu, yu);
3417  pmrk.SetPoint(pmrkIdx++, zu, yu);
3418  }
3419  }
3420  }
3421 
3422  // By convention we will have drawn the "best" axis first
3423  if (recoOpt->fBestPCAAxisOnly) break;
3424 
3425  pcaIdx++;
3426  }
3427  }
3428  }
3429 
3430  // Now let's loop over daughters and call drawing routine for them
3431  if (pfPart->NumDaughters() > 0) {
3432  depth++;
3433 
3434  for (const auto& daughterIdx : pfPart->Daughters()) {
3435  DrawPFParticleOrtho(pfParticleVec.at(daughterIdx),
3436  pfParticleVec,
3437  spacePointAssnVec,
3438  pcAxisAssnVec,
3439  depth,
3440  proj,
3441  view);
3442  }
3443  }
3444 
3445  return;
3446  }
3447 
3448  //......................................................................
3449  void
3450  RecoBaseDrawer::ProngOrtho(const art::Event& evt,
3451  evd::OrthoProj_t proj,
3452  double msize,
3453  evdb::View2D* view)
3454  {
3455  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
3456  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3457 
3458  if (rawOpt->fDrawRawDataOrCalibWires < 1) return;
3459 
3460  // annoying for now, but have to have multiple copies of basically the
3461  // same code to draw prongs, showers and tracks so that we can use
3462  // the art::Assns to get the hits and clusters.
3463 
3464  // Tracks.
3465 
3466  if (recoOpt->fDrawTracks != 0) {
3467  for (size_t imod = 0; imod < recoOpt->fTrackLabels.size(); ++imod) {
3468  art::InputTag which = recoOpt->fTrackLabels[imod];
3469  art::View<recob::Track> track;
3470  this->GetTracks(evt, which, track);
3471 
3472  for (size_t t = 0; t < track.vals().size(); ++t) {
3473  const recob::Track* ptrack = track.vals().at(t);
3474  int color = ptrack->ID() & 65535;
3475 
3476  // Draw track using only embedded information.
3477 
3478  DrawTrackOrtho(*ptrack, color, proj, msize, view);
3479  }
3480  }
3481  }
3482 
3483  // Showers.
3484 
3485  if (recoOpt->fDrawShowers != 0) {
3486  for (size_t imod = 0; imod < recoOpt->fShowerLabels.size(); ++imod) {
3487  art::InputTag which = recoOpt->fShowerLabels[imod];
3488  art::View<recob::Shower> shower;
3489  this->GetShowers(evt, which, shower);
3490 
3491  for (size_t s = 0; s < shower.vals().size(); ++s) {
3492  const recob::Shower* pshower = shower.vals().at(s);
3493  int color = pshower->ID();
3494  DrawShowerOrtho(*pshower, color, proj, msize, view);
3495  }
3496  }
3497  }
3498 
3499  return;
3500  }
3501 
3502  //......................................................................
3503  void
3504  RecoBaseDrawer::DrawSpacePointOrtho(std::vector<art::Ptr<recob::SpacePoint>>& spts,
3505  int color,
3506  evd::OrthoProj_t proj,
3507  double msize,
3508  evdb::View2D* view,
3509  int mode)
3510  {
3511  // Get services.
3512 
3513  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3514 
3515  // Organize space points into separate collections according to the color
3516  // we want them to be.
3517  // If option If option fColorSpacePointsByChisq is false, this means
3518  // having a single collection with color inherited from the prong
3519  // (specified by the argument color).
3520 
3521  std::map<int, std::vector<art::Ptr<recob::SpacePoint>>> spmap; // Indexed by color.
3522 
3523  for (auto& pspt : spts) {
3524 
3525  // By default use event display palette.
3526 
3527  int spcolor = evd::kColor[color % evd::kNCOLS];
3528  if (mode == 1) { //shower hits
3529  spcolor = evd::kColor2[color % evd::kNCOLS];
3530  }
3531  // For rainbow effect, choose root colors in range [51,100].
3532  // We are using 100=best (red), 51=worst (blue).
3533 
3534  if (recoOpt->fColorSpacePointsByChisq) {
3535  spcolor = 100 - 2.5 * pspt->Chisq();
3536  if (spcolor < 51) spcolor = 51;
3537  if (spcolor > 100) spcolor = 100;
3538  }
3539  spmap[spcolor].push_back(pspt);
3540  }
3541 
3542  // Loop over colors.
3543  // Note that larger (=better) space points are plotted on
3544  // top for optimal visibility.
3545 
3546  for (auto icolor : spmap) {
3547  int spcolor = icolor.first;
3548  std::vector<art::Ptr<recob::SpacePoint>>& psps = icolor.second;
3549 
3550  // Make and fill a polymarker.
3551 
3552  TPolyMarker& pm = view->AddPolyMarker(psps.size(), spcolor, kFullCircle, msize);
3553  for (size_t s = 0; s < psps.size(); ++s) {
3554  const recob::SpacePoint& spt = *psps[s];
3555  const double* xyz = spt.XYZ();
3556  switch (proj) {
3557  case evd::kXY: pm.SetPoint(s, xyz[0], xyz[1]); break;
3558  case evd::kXZ: pm.SetPoint(s, xyz[2], xyz[0]); break;
3559  case evd::kYZ: pm.SetPoint(s, xyz[2], xyz[1]); break;
3560  default:
3561  throw cet::exception("RecoBaseDrawer")
3562  << __func__ << ": unknown projection #" << ((int)proj) << "\n";
3563  } // switch
3564  }
3565  }
3566 
3567  return;
3568  }
3569 
3570  //......................................................................
3571  void
3572  RecoBaseDrawer::DrawTrackOrtho(const recob::Track& track,
3573  int color,
3574  evd::OrthoProj_t proj,
3575  double msize,
3576  evdb::View2D* view)
3577  {
3578  // Get options.
3579 
3580  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
3581 
3582  if (recoOpt->fDrawTrackSpacePoints) {
3583 
3584  // Use brute force to find the module label and index of this
3585  // track, so that we can find associated space points and draw
3586  // them.
3587 
3588  const art::Event* evt = evdb::EventHolder::Instance()->GetEvent();
3589  //std::vector<art::Handle<std::vector<recob::Track>>> handles;
3590  //evt->getManyByType(handles);
3591  auto handles = evt->getMany<std::vector<recob::Track>>();
3592  for (auto ih : handles) {
3593  const art::Handle<std::vector<recob::Track>> handle = ih;
3594  if (handle.isValid()) {
3595  const std::string& which = handle.provenance()->moduleLabel();
3596  art::FindManyP<recob::SpacePoint> fmsp(handle, *evt, which);
3597 
3598  int n = handle->size();
3599  for (int i = 0; i < n; ++i) {
3600  art::Ptr<recob::Track> p(handle, i);
3601  if (&*p == &track) {
3602  std::vector<art::Ptr<recob::SpacePoint>> spts = fmsp.at(i);
3603  DrawSpacePointOrtho(spts, color, proj, msize, view);
3604  }
3605  }
3606  }
3607  }
3608  }
3609  if (recoOpt->fDrawTrackTrajectoryPoints) {
3610 
3611  // Draw trajectory points.
3612 
3613  int np = track.NumberTrajectoryPoints();
3614  int vp = track.CountValidPoints();
3615 
3616  // Make and fill a polymarker.
3617 
3618  TPolyMarker& pm =
3619  view->AddPolyMarker(vp, evd::kColor[color % evd::kNCOLS], kFullCircle, msize);
3620  TPolyLine& pl = view->AddPolyLine(vp, evd::kColor[color % evd::kNCOLS], 2, 0);
3621  for (int p = 0; p < np; ++p) {
3622  if (track.HasValidPoint(p) == 0) continue;
3623  const auto& pos = track.LocationAtPoint(p);
3624  switch (proj) {
3625  case evd::kXY:
3626  pm.SetPoint(p, pos.X(), pos.Y());
3627  pl.SetPoint(p, pos.X(), pos.Y());
3628  break;
3629  case evd::kXZ:
3630  pm.SetPoint(p, pos.Z(), pos.X());
3631  pl.SetPoint(p, pos.Z(), pos.X());
3632  break;
3633  case evd::kYZ:
3634  pm.SetPoint(p, pos.Z(), pos.Y());
3635  pl.SetPoint(p, pos.Z(), pos.Y());
3636  break;
3637  default:
3638  throw cet::exception("RecoBaseDrawer")
3639  << __func__ << ": unknown projection #" << ((int)proj) << "\n";
3640  } // switch
3641  } // p
3642  // BB: draw the track ID at the end of the track
3643  if (recoOpt->fDrawTracks > 1) {
3644  int tid =
3645  track.ID() &
3646  65535; //this is a hack for PMA track id which uses the 16th bit to identify shower-like track.
3647  std::string s = std::to_string(tid);
3648  char const* txt = s.c_str();
3649  double x = track.End().X();
3650  double y = track.End().Y();
3651  double z = track.End().Z();
3652  if (proj == evd::kXY) {
3653  TText& trkID = view->AddText(x, y, txt);
3654  trkID.SetTextColor(evd::kColor[tid % evd::kNCOLS]);
3655  trkID.SetTextSize(0.03);
3656  }
3657  else if (proj == evd::kXZ) {
3658  TText& trkID = view->AddText(z, x, txt);
3659  trkID.SetTextColor(evd::kColor[tid % evd::kNCOLS]);
3660  trkID.SetTextSize(0.03);
3661  }
3662  else if (proj == evd::kYZ) {
3663  TText& trkID = view->AddText(z, y, txt);
3664  trkID.SetTextColor(evd::kColor[tid % evd::kNCOLS]);
3665  trkID.SetTextSize(0.03);
3666  } // proj
3667  } // recoOpt->fDrawTracks > 1
3668  }
3669 
3670  return;
3671  }
3672 
3673  //......................................................................
3674  void
3675  RecoBaseDrawer::DrawShowerOrtho(const recob::Shower& shower,
3676  int color,
3677  evd::OrthoProj_t proj,
3678  double msize,
3679  evdb::View2D* view)
3680  {
3681  // Use brute force to find the module label and index of this
3682  // shower, so that we can find associated space points and draw
3683  // them.
3684 
3685  const art::Event* evt = evdb::EventHolder::Instance()->GetEvent();
3686  //std::vector<art::Handle<std::vector<recob::Shower>>> handles;
3687  //evt->getManyByType(handles);
3688  auto handles = evt->getMany<std::vector<recob::Shower>>();
3689  for (auto ih : handles) {
3690  const art::Handle<std::vector<recob::Shower>> handle = ih;
3691  if (handle.isValid()) {
3692  const std::string& which = handle.provenance()->moduleLabel();
3693  art::FindManyP<recob::SpacePoint> fmsp(handle, *evt, which);
3694  if (!fmsp.isValid()) continue;
3695  int n = handle->size();
3696  for (int i = 0; i < n; ++i) {
3697  art::Ptr<recob::Shower> p(handle, i);
3698  if (&*p == &shower) {
3699  switch (proj) {
3700  case evd::kXY:
3701  view->AddMarker(p->ShowerStart().X(),
3702  p->ShowerStart().Y(),
3703  evd::kColor2[color % evd::kNCOLS],
3704  5,
3705  2.0);
3706  break;
3707  case evd::kXZ:
3708  view->AddMarker(p->ShowerStart().Z(),
3709  p->ShowerStart().X(),
3710  evd::kColor2[color % evd::kNCOLS],
3711  5,
3712  2.0);
3713  break;
3714  case evd::kYZ:
3715  view->AddMarker(p->ShowerStart().Z(),
3716  p->ShowerStart().Y(),
3717  evd::kColor2[color % evd::kNCOLS],
3718  5,
3719  2.0);
3720  break;
3721  default:
3722  throw cet::exception("RecoBaseDrawer")
3723  << __func__ << ": unknown projection #" << ((int)proj) << "\n";
3724  } // switch
3725 
3726  if (fmsp.isValid()) {
3727  std::vector<art::Ptr<recob::SpacePoint>> spts = fmsp.at(i);
3728  DrawSpacePointOrtho(spts, color, proj, msize, view, 1);
3729  }
3730  }
3731  }
3732  }
3733  }
3734 
3735  return;
3736  }
3737 
3738  //......................................................................
3739  int
3740  RecoBaseDrawer::GetWires(const art::Event& evt,
3741  const art::InputTag& which,
3742  art::PtrVector<recob::Wire>& wires)
3743  {
3744  wires.clear();
3745 
3746  art::Handle<std::vector<recob::Wire>> wcol;
3747  art::PtrVector<recob::Wire> temp;
3748 
3749  try {
3750  evt.getByLabel(which, wcol);
3751 
3752  for (unsigned int i = 0; i < wcol->size(); ++i) {
3753  art::Ptr<recob::Wire> w(wcol, i);
3754  temp.push_back(w);
3755  }
3756  temp.swap(wires);
3757  }
3758  catch (cet::exception& e) {
3759  writeErrMsg("GetWires", e);
3760  }
3761 
3762  return wires.size();
3763  }
3764 
3765  //......................................................................
3766  int
3767  RecoBaseDrawer::GetHits(const art::Event& evt,
3768  const art::InputTag& which,
3769  std::vector<const recob::Hit*>& hits,
3770  unsigned int plane)
3771  {
3772  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
3773  art::ServiceHandle<geo::Geometry const> geo;
3774 
3775  hits.clear();
3776 
3777  std::vector<const recob::Hit*> temp;
3778 
3779  try {
3780  evt.getView(which, temp);
3781  for (const auto& hit : temp) {
3782  // Note that the WireID in the hit object is useless for those detectors where a channel can correspond to
3783  // more than one plane/wire. So our plan is to recover the list of wire IDs from the channel number and
3784  // loop over those (if there are any)
3785  const std::vector<geo::WireID>& wireIDs = geo->ChannelToWire(hit->Channel());
3786 
3787  // Loop to find match
3788  for (const auto& wireID : wireIDs) {
3789  if (wireID.Plane == plane && wireID.TPC == rawOpt->fTPC &&
3790  wireID.Cryostat == rawOpt->fCryostat)
3791  hits.push_back(hit);
3792  }
3793  }
3794  }
3795  catch (cet::exception& e) {
3796  writeErrMsg("GetHits", e);
3797  }
3798 
3799  return hits.size();
3800  }
3801 
3802  //......................................................................
3803  int
3804  RecoBaseDrawer::GetSlices(const art::Event& evt,
3805  const art::InputTag& which,
3806  art::PtrVector<recob::Slice>& slices)
3807  {
3808  slices.clear();
3809  art::PtrVector<recob::Slice> temp;
3810 
3811  art::Handle<std::vector<recob::Slice>> slcCol;
3812 
3813  try {
3814  evt.getByLabel(which, slcCol);
3815  temp.reserve(slcCol->size());
3816  for (unsigned int i = 0; i < slcCol->size(); ++i) {
3817  art::Ptr<recob::Slice> slc(slcCol, i);
3818  temp.push_back(slc);
3819  }
3820  temp.swap(slices);
3821  }
3822  catch (cet::exception& e) {
3823  writeErrMsg("GetSlices", e);
3824  }
3825 
3826  return slices.size();
3827  }
3828 
3829  //......................................................................
3830  int
3831  RecoBaseDrawer::GetClusters(const art::Event& evt,
3832  const art::InputTag& which,
3833  art::PtrVector<recob::Cluster>& clust)
3834  {
3835  clust.clear();
3836  art::PtrVector<recob::Cluster> temp;
3837 
3838  art::Handle<std::vector<recob::Cluster>> clcol;
3839 
3840  try {
3841  evt.getByLabel(which, clcol);
3842  temp.reserve(clcol->size());
3843  for (unsigned int i = 0; i < clcol->size(); ++i) {
3844  art::Ptr<recob::Cluster> cl(clcol, i);
3845  temp.push_back(cl);
3846  }
3847  temp.swap(clust);
3848  }
3849  catch (cet::exception& e) {
3850  writeErrMsg("GetClusters", e);
3851  }
3852 
3853  return clust.size();
3854  }
3855 
3856  //......................................................................
3857  int
3858  RecoBaseDrawer::GetPFParticles(const art::Event& evt,
3859  const art::InputTag& which,
3860  art::PtrVector<recob::PFParticle>& clust)
3861  {
3862  clust.clear();
3863  art::PtrVector<recob::PFParticle> temp;
3864 
3865  art::Handle<std::vector<recob::PFParticle>> clcol;
3866 
3867  try {
3868  evt.getByLabel(which, clcol);
3869  for (unsigned int i = 0; i < clcol->size(); ++i) {
3870  art::Ptr<recob::PFParticle> cl(clcol, i);
3871  temp.push_back(cl);
3872  }
3873  temp.swap(clust);
3874  }
3875  catch (cet::exception& e) {
3876  writeErrMsg("GetPFParticles", e);
3877  }
3878 
3879  return clust.size();
3880  }
3881 
3882  //......................................................................
3883  int
3884  RecoBaseDrawer::GetEndPoint2D(const art::Event& evt,
3885  const art::InputTag& which,
3886  art::PtrVector<recob::EndPoint2D>& ep2d)
3887  {
3888  ep2d.clear();
3889  art::PtrVector<recob::EndPoint2D> temp;
3890 
3891  art::Handle<std::vector<recob::EndPoint2D>> epcol;
3892 
3893  try {
3894  evt.getByLabel(which, epcol);
3895  for (unsigned int i = 0; i < epcol->size(); ++i) {
3896  art::Ptr<recob::EndPoint2D> ep(epcol, i);
3897  temp.push_back(ep);
3898  }
3899  temp.swap(ep2d);
3900  }
3901  catch (cet::exception& e) {
3902  writeErrMsg("GetEndPoint2D", e);
3903  }
3904 
3905  return ep2d.size();
3906  }
3907 
3908  //......................................................................
3909 
3910  int
3911  RecoBaseDrawer::GetOpFlashes(const art::Event& evt,
3912  const art::InputTag& which,
3913  art::PtrVector<recob::OpFlash>& opflashes)
3914  {
3915  opflashes.clear();
3916  art::PtrVector<recob::OpFlash> temp;
3917 
3918  art::Handle<std::vector<recob::OpFlash>> opflashcol;
3919 
3920  try {
3921  evt.getByLabel(which, opflashcol);
3922  for (unsigned int i = 0; i < opflashcol->size(); ++i) {
3923  art::Ptr<recob::OpFlash> opf(opflashcol, i);
3924  temp.push_back(opf);
3925  }
3926  temp.swap(opflashes);
3927  }
3928  catch (cet::exception& e) {
3929  writeErrMsg("GetOpFlashes", e);
3930  }
3931 
3932  return opflashes.size();
3933  }
3934 
3935  //......................................................................
3936 
3937  int
3938  RecoBaseDrawer::GetSeeds(const art::Event& evt,
3939  const art::InputTag& which,
3940  art::PtrVector<recob::Seed>& seeds)
3941  {
3942  seeds.clear();
3943  art::PtrVector<recob::Seed> temp;
3944 
3945  art::Handle<std::vector<recob::Seed>> seedcol;
3946 
3947  try {
3948  evt.getByLabel(which, seedcol);
3949  for (unsigned int i = 0; i < seedcol->size(); ++i) {
3950  art::Ptr<recob::Seed> sd(seedcol, i);
3951  temp.push_back(sd);
3952  }
3953  temp.swap(seeds);
3954  }
3955  catch (cet::exception& e) {
3956  writeErrMsg("GetSeeds", e);
3957  }
3958 
3959  return seeds.size();
3960  }
3961 
3962  //......................................................................
3963  int
3964  RecoBaseDrawer::GetSpacePoints(const art::Event& evt,
3965  const art::InputTag& which,
3966  std::vector<art::Ptr<recob::SpacePoint>>& spts)
3967  {
3968  spts.clear();
3969  art::Handle<std::vector<recob::SpacePoint>> spcol;
3970  if (evt.getByLabel(which, spcol)) art::fill_ptr_vector(spts, spcol);
3971 
3972  return spts.size();
3973  }
3974 
3975  //......................................................................
3976  int
3977  RecoBaseDrawer::GetEdges(const art::Event& evt,
3978  const art::InputTag& which,
3979  std::vector<art::Ptr<recob::Edge>>& edges)
3980  {
3981  edges.clear();
3982 
3983  art::Handle<std::vector<recob::Edge>> edgeCol;
3984 
3985  evt.getByLabel(which, edgeCol);
3986 
3987  for (unsigned int i = 0; i < edgeCol->size(); ++i)
3988  edges.emplace_back(edgeCol, i);
3989 
3990  return edges.size();
3991  }
3992 
3993  //......................................................................
3994  int
3995  RecoBaseDrawer::GetTracks(const art::Event& evt,
3996  const art::InputTag& which,
3997  art::View<recob::Track>& track)
3998  {
3999  try {
4000  evt.getView(which, track);
4001  }
4002  catch (cet::exception& e) {
4003  writeErrMsg("GetTracks", e);
4004  }
4005 
4006  return track.vals().size();
4007  }
4008 
4009  //......................................................................
4010  int
4011  RecoBaseDrawer::GetShowers(const art::Event& evt,
4012  const art::InputTag& which,
4013  art::View<recob::Shower>& shower)
4014  {
4015  try {
4016  evt.getView(which, shower);
4017  }
4018  catch (cet::exception& e) {
4019  writeErrMsg("GetShowers", e);
4020  }
4021 
4022  return shower.vals().size();
4023  }
4024 
4025  //......................................................................
4026  int
4027  RecoBaseDrawer::GetVertices(const art::Event& evt,
4028  const art::InputTag& which,
4029  art::PtrVector<recob::Vertex>& vertex)
4030  {
4031  vertex.clear();
4032  art::PtrVector<recob::Vertex> temp;
4033 
4034  art::Handle<std::vector<recob::Vertex>> vcol;
4035 
4036  try {
4037  evt.getByLabel(which, vcol);
4038  for (size_t i = 0; i < vcol->size(); ++i) {
4039  art::Ptr<recob::Vertex> v(vcol, i);
4040  temp.push_back(v);
4041  }
4042  temp.swap(vertex);
4043  }
4044  catch (cet::exception& e) {
4045  writeErrMsg("GetVertices", e);
4046  }
4047 
4048  return vertex.size();
4049  }
4050 
4051  //......................................................................
4052  int
4053  RecoBaseDrawer::GetEvents(const art::Event& evt,
4054  const art::InputTag& which,
4055  art::PtrVector<recob::Event>& event)
4056  {
4057  event.clear();
4058  art::PtrVector<recob::Event> temp;
4059 
4060  art::Handle<std::vector<recob::Event>> ecol;
4061 
4062  try {
4063  evt.getByLabel(which, ecol);
4064  for (size_t i = 0; i < ecol->size(); ++i) {
4065  art::Ptr<recob::Event> e(ecol, i);
4066  temp.push_back(e);
4067  }
4068  temp.swap(event);
4069  }
4070  catch (cet::exception& e) {
4071  writeErrMsg("GetEvents", e);
4072  }
4073 
4074  return event.size();
4075  }
4076 
4077  //......................................................................
4078  int
4079  RecoBaseDrawer::CountHits(const art::Event& evt,
4080  const art::InputTag& which,
4081  unsigned int cryostat,
4082  unsigned int tpc,
4083  unsigned int plane)
4084  {
4085  std::vector<const recob::Hit*> temp;
4086  int NumberOfHitsBeforeThisPlane = 0;
4087  evt.getView(
4088  which,
4089  temp); //temp.size() = total number of hits for this event (number of all hits in all Cryostats, TPC's, planes and wires)
4090  for (size_t t = 0; t < temp.size(); ++t) {
4091  if (temp[t]->WireID().Cryostat == cryostat && temp[t]->WireID().TPC == tpc &&
4092  temp[t]->WireID().Plane == plane)
4093  break;
4094  NumberOfHitsBeforeThisPlane++;
4095  }
4096  return NumberOfHitsBeforeThisPlane;
4097  }
4098 
4099  //......................................................................
4100  void
4101  RecoBaseDrawer::FillTQHisto(const art::Event& evt,
4102  unsigned int plane,
4103  unsigned int wire,
4104  TH1F* histo)
4105  {
4106  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
4107  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
4108  art::ServiceHandle<geo::Geometry const> geo;
4109 
4110  float minSig(std::numeric_limits<float>::max());
4111  float maxSig(std::numeric_limits<float>::lowest());
4112  bool setLimits(false);
4113 
4114  // Check if we're supposed to draw raw hits at all
4115  if (rawOpt->fDrawRawDataOrCalibWires == 0) return;
4116 
4117  for (size_t imod = 0; imod < recoOpt->fWireLabels.size(); ++imod) {
4118  art::InputTag const which = recoOpt->fWireLabels[imod];
4119 
4120  art::PtrVector<recob::Wire> wires;
4121  this->GetWires(evt, which, wires);
4122 
4123  for (size_t i = 0; i < wires.size(); ++i) {
4124 
4125  std::vector<geo::WireID> wireids = geo->ChannelToWire(wires[i]->Channel());
4126 
4127  bool goodWID = false;
4128  for (auto const& wid : wireids) {
4129  // check for correct plane, wire and tpc
4130  if (wid.Plane == plane && wid.Wire == wire && wid.TPC == rawOpt->fTPC &&
4131  wid.Cryostat == rawOpt->fCryostat)
4132  goodWID = true;
4133  }
4134  if (!goodWID) continue;
4135 
4136  std::vector<float> wirSig = wires[i]->Signal();
4137  for (unsigned int ii = 0; ii < wirSig.size(); ++ii) {
4138  // histo->SetLineColor(imod+4);
4139  // histo->Fill(1.*ii, wirSig[ii]);
4140  minSig = std::min(minSig, wirSig[ii]);
4141  maxSig = std::max(maxSig, wirSig[ii]);
4142  }
4143 
4144  setLimits = true;
4145  } //end loop over wires
4146  } //end loop over wire modules
4147 
4148  if (setLimits) {
4149  histo->SetMaximum(1.2 * maxSig);
4150  histo->SetMinimum(1.2 * minSig);
4151  }
4152 
4153  return;
4154  }
4155 
4156  //......................................................................
4157  void
4158  RecoBaseDrawer::FillQHisto(const art::Event& evt, unsigned int plane, TH1F* histo)
4159  {
4160  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
4161  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
4162  art::ServiceHandle<geo::Geometry const> geo;
4163 
4164  // Check if we're supposed to draw raw hits at all
4165  if (rawOpt->fDrawRawDataOrCalibWires == 0) return;
4166 
4167  for (size_t imod = 0; imod < recoOpt->fWireLabels.size(); ++imod) {
4168  art::InputTag const which = recoOpt->fWireLabels[imod];
4169 
4170  art::PtrVector<recob::Wire> wires;
4171  this->GetWires(evt, which, wires);
4172 
4173  for (unsigned int i = 0; i < wires.size(); ++i) {
4174 
4175  std::vector<geo::WireID> wireids = geo->ChannelToWire(wires[i]->Channel());
4176 
4177  bool goodWID = false;
4178  for (auto const& wid : wireids) {
4179  // check for correct plane, wire and tpc
4180  if (wid.Plane == plane && wid.TPC == rawOpt->fTPC && wid.Cryostat == rawOpt->fCryostat)
4181  goodWID = true;
4182  }
4183  if (!goodWID) continue;
4184  std::vector<float> wirSig = wires[i]->Signal();
4185  for (unsigned int ii = 0; ii < wirSig.size(); ++ii)
4186  histo->Fill(wirSig[ii]);
4187  /*
4188  for(size_t s = 0; s < wires[i]->NSignal(); ++s)
4189  histo->Fill(wires[i]->Signal()[s]);
4190 */
4191 
4192  } //end loop over raw hits
4193  } //end loop over Wire modules
4194 
4195  return;
4196  }
4197 
4198  //......................................................................
4199  void
4200  RecoBaseDrawer::FillTQHistoDP(const art::Event& evt,
4201  unsigned int plane,
4202  unsigned int wire,
4203  TH1F* histo,
4204  std::vector<double>& htau1,
4205  std::vector<double>& htau2,
4206  std::vector<double>& hitamplitudes,
4207  std::vector<double>& hpeaktimes,
4208  std::vector<int>& hstartT,
4209  std::vector<int>& hendT,
4210  std::vector<int>& hNMultiHit,
4211  std::vector<int>& hLocalHitIndex)
4212  {
4213  art::ServiceHandle<evd::RawDrawingOptions const> rawOpt;
4214  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
4215  art::ServiceHandle<geo::Geometry const> geo;
4216 
4217  // Check if we're supposed to draw raw hits at all
4218  if (rawOpt->fDrawRawDataOrCalibWires == 0) return;
4219 
4220  for (size_t imod = 0; imod < recoOpt->fWireLabels.size(); ++imod) {
4221  art::InputTag const which = recoOpt->fWireLabels[imod];
4222 
4223  art::PtrVector<recob::Wire> wires;
4224  this->GetWires(evt, which, wires);
4225 
4226  for (size_t i = 0; i < wires.size(); ++i) {
4227 
4228  std::vector<geo::WireID> wireids = geo->ChannelToWire(wires[i]->Channel());
4229 
4230  bool goodWID = false;
4231  for (auto const& wid : wireids) {
4232  if (wid.Plane == plane && wid.Wire == wire && wid.TPC == rawOpt->fTPC &&
4233  wid.Cryostat == rawOpt->fCryostat)
4234  goodWID = true;
4235  }
4236 
4237  if (!goodWID) continue;
4238 
4239  std::vector<float> wirSig = wires[i]->Signal();
4240  for (unsigned int ii = 0; ii < wirSig.size(); ++ii)
4241  histo->Fill(1. * ii, wirSig[ii]);
4242  break;
4243  } //end loop over wires
4244  } //end loop over wire modules
4245 
4246  for (size_t imod = 0; imod < recoOpt->fHitLabels.size(); ++imod) {
4247  art::InputTag const which = recoOpt->fHitLabels[imod];
4248 
4249  std::vector<const recob::Hit*> hits;
4250  this->GetHits(evt, which, hits, plane);
4251 
4252  auto hitResults = anab::FVectorReader<recob::Hit, 4>::create(evt, "dprawhit");
4253  const auto& fitParams = hitResults->vectors();
4254 
4255  int FitParamsOffset = CountHits(evt, which, rawOpt->fCryostat, rawOpt->fTPC, plane);
4256 
4257  for (size_t i = 0; i < hits.size(); ++i) {
4258  // check for correct wire. Plane, cryostat and tpc were checked in GetHits
4259  if (hits[i]->WireID().Wire != wire) continue;
4260 
4261  hpeaktimes.push_back(fitParams[FitParamsOffset + i][0]);
4262  htau1.push_back(fitParams[FitParamsOffset + i][1]);
4263  htau2.push_back(fitParams[FitParamsOffset + i][2]);
4264  hitamplitudes.push_back(fitParams[FitParamsOffset + i][3]);
4265  hstartT.push_back(hits[i]->StartTick());
4266  hendT.push_back(hits[i]->EndTick());
4267  hNMultiHit.push_back(hits[i]->Multiplicity());
4268  hLocalHitIndex.push_back(hits[i]->LocalIndex());
4269  } //end loop over reco hits
4270  } //end loop over HitFinding modules
4271 
4272  return;
4273  }
4274 
4275  //......................................................................
4276  //double RecoBaseDrawer::EvalExpoFit(double x,
4277  // double tau1,
4278  // double tau2,
4279  // double amplitude,
4280  // double peaktime)
4281  //{
4282  //return (amplitude * exp(0.4*(x-peaktime)/tau1) / ( 1 + exp(0.4*(x-peaktime)/tau2) ) );
4283  //}
4284 
4285  //......................................................................
4286  //double RecoBaseDrawer::EvalMultiExpoFit(double x,
4287  // int HitNumber,
4288  // int NHits,
4289  // std::vector<double> tau1,
4290  // std::vector<double> tau2,
4291  // std::vector<double> amplitude,
4292  // std::vector<double> peaktime)
4293  //{
4294  // double x_sum = 0.;
4295  //
4296  // for(int i = HitNumber; i < HitNumber+NHits; i++)
4297  // {
4298  // x_sum += (amplitude[i] * exp(0.4*(x-peaktime[i])/tau1[i]) / ( 1 + exp(0.4*(x-peaktime[i])/tau2[i]) ) );
4299  // }
4300  //
4301  //return x_sum;
4302  //}
4303 
4304 } // namespace evd
4305 ////////////////////////////////////////////////////////////////////////
evd::RecoBaseDrawer::fSpacePointDrawer
ISpacePointDrawerPtr fSpacePointDrawer
Definition: RecoBaseDrawer.h:332
evd::RecoBaseDrawer::SpacePointOrtho
void SpacePointOrtho(const art::Event &evt, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3123
vertex
process_name vertex
Definition: cheaterreco.fcl:51
evd::RecoBaseDrawer::DrawShowerOrtho
void DrawShowerOrtho(const recob::Shower &shower, int color, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3675
evd::RecoBaseDrawer::SpacePoint3D
void SpacePoint3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:1986
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
evd::RecoBaseDrawer::GetTracks
int GetTracks(const art::Event &evt, const art::InputTag &which, art::View< recob::Track > &track)
Definition: RecoBaseDrawer.cxx:3995
evd::RecoBaseDrawer::fRawCharge
std::vector< double > fRawCharge
Sum of Raw Charge.
Definition: RecoBaseDrawer.h:339
evd::RecoBaseDrawer::FillTQHisto
void FillTQHisto(const art::Event &evt, unsigned int plane, unsigned int wire, TH1F *histo)
Definition: RecoBaseDrawer.cxx:4101
evd::RecoBaseDrawer::PFParticleOrtho
void PFParticleOrtho(const art::Event &evt, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3154
labels
* labels
Definition: build-larbase-py2.sh:59
Edge.h
An object to define a &quot;edge&quot; which is used to connect space points in a triangulation algorithm...
evd::RecoBaseDrawer::DrawProng2D
void DrawProng2D(detinfo::DetectorPropertiesData const &detProp, std::vector< const recob::Hit * > &hits, evdb::View2D *view, unsigned int plane, TVector3 const &startPos, TVector3 const &startDir, int id, float cscore=-5)
Definition: RecoBaseDrawer.cxx:1232
Plane
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
Definition: rawdigitfilter_icarus.fcl:8
recob::Shower::has_length
bool has_length() const
Returns whether the shower has a valid length.
Definition: Shower.h:211
evd::RecoBaseDrawer::fAllSpacePointDrawer
ISpacePointDrawerPtr fAllSpacePointDrawer
Definition: RecoBaseDrawer.h:331
tracks
ClusterModuleLabel join with tracks
Definition: trackfinderalgorithms.fcl:56
CryostatGeo.h
Encapsulate the construction of a single cyostat.
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
lariov::ChannelStatusProvider::IsBad
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
evd::RecoBaseDrawer::Prong3D
void Prong3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2560
evd::RecoBaseDrawer::SpacePointChiSq
float SpacePointChiSq(const std::vector< art::Ptr< recob::Hit >> &) const
Definition: RecoBaseDrawer.cxx:2112
evd::RecoBaseDrawer::DrawPFParticle3D
void DrawPFParticle3D(const art::Ptr< recob::PFParticle > &pfPart, const art::PtrVector< recob::PFParticle > &pfParticleVec, const std::vector< art::Ptr< recob::SpacePoint >> &spacePointVec, const art::FindManyP< recob::Edge > &edgeAssnsVec, const art::FindManyP< recob::SpacePoint > &spacePointAssnsVec, const art::FindManyP< recob::SpacePoint > &edgeSPAssnVec, const art::FindManyP< recob::Hit > &spHitAssnVec, const art::FindMany< recob::Track > &trackAssnVec, const art::FindMany< recob::PCAxis > &pcAxisAssnVec, const art::FindMany< anab::CosmicTag > &cosmicTagAssnVec, int depth, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2158
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
evd::RecoBaseDrawer::GetSpacePoints
int GetSpacePoints(const art::Event &evt, const art::InputTag &which, std::vector< art::Ptr< recob::SpacePoint >> &spts)
Definition: RecoBaseDrawer.cxx:3964
detinfo::DetectorPropertiesData::GetXTicksOffset
double GetXTicksOffset(int p, int t, int c) const
Definition: DetectorPropertiesData.cc:148
evd::RecoBaseDrawer::SeedOrtho
void SeedOrtho(const art::Event &evt, evd::OrthoProj_t proj, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:1926
recob::Shower::Direction
const TVector3 & Direction() const
Definition: Shower.h:189
Hit.h
Declaration of signal hit object.
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:126
p
pdgs p
Definition: selectors.fcl:22
recob::Shower::has_open_angle
bool has_open_angle() const
Returns whether the shower has a valid opening angle.
Definition: Shower.h:210
recob::Shower::dEdx
const std::vector< double > & dEdx() const
Definition: Shower.h:203
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
evd::RecoBaseDrawer::Vertex3D
void Vertex3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2857
recob::Track::HasValidPoint
bool HasValidPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:111
evd::RecoBaseDrawer::Event3D
void Event3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2907
anab::FVectorReader::create
static std::unique_ptr< FVectorReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:29
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: bj/lardataobj/RecoBase/Track.h:102
StartTick
BEGIN_PROLOG StartTick
Definition: decoderTools_icarus.fcl:7
recob::Track::VertexDirection
Vector_t VertexDirection() const
Definition: bj/lardataobj/RecoBase/Track.h:132
evd::RecoBaseDrawer::ProngOrtho
void ProngOrtho(const art::Event &evt, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3450
evd::RecoBaseDrawer::GetClusters
int GetClusters(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Cluster > &clust)
Definition: RecoBaseDrawer.cxx:3831
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
evd::RecoBaseDrawer::DrawShower3D
void DrawShower3D(const recob::Shower &shower, int color, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2750
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
evd::kNCOLS
static const int kNCOLS
Definition: eventdisplay.h:10
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
util.style
def style
Definition: util.py:237
evd::RecoBaseDrawer::DrawSpacePointOrtho
void DrawSpacePointOrtho(std::vector< art::Ptr< recob::SpacePoint >> &spts, int color, evd::OrthoProj_t proj, double msize, evdb::View2D *view, int mode=0)
0: track, 1: shower
Definition: RecoBaseDrawer.cxx:3504
evd::RecoBaseDrawer::Edge3D
void Edge3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2453
hit
process_name hit
Definition: cheaterreco.fcl:51
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
evd::OrthoProj_t
OrthoProj_t
Definition: OrthoProj.h:12
evd::RecoBaseDrawer::GetEvents
int GetEvents(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Event > &event)
Definition: RecoBaseDrawer.cxx:4053
evd::RecoBaseDrawer::Prong2D
void Prong2D(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:1391
shower
process_name shower
Definition: cheaterreco.fcl:51
anab::CosmicTag::CosmicScore
float & CosmicScore()
Definition: CosmicTag.h:61
Track.h
Offers proxy::Tracks and proxy::Track class for recob::Track access.
evd::RecoBaseDrawer::fConvertedCharge
std::vector< double > fConvertedCharge
Sum of Charge Converted using Birks&#39; formula.
Definition: RecoBaseDrawer.h:340
recob::Shower::Length
double Length() const
Definition: Shower.h:201
evd::RecoBaseDrawer::Slice2D
void Slice2D(const art::Event &evt, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:726
evd::RecoBaseDrawer::Cluster2D
void Cluster2D(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:818
evd::RecoBaseDrawer::DrawTrackVertexAssns2D
void DrawTrackVertexAssns2D(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:1612
ColorDrawingOptions.h
The color scales used by the event display.
local
then local
Definition: make_xml_premcc1.0.sh:85
recob::Track::CountValidPoints
unsigned int CountValidPoints() const
Definition: bj/lardataobj/RecoBase/Track.h:112
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
evd::RecoBaseDrawer::Seed3D
void Seed3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:1880
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
evd::RecoBaseDrawer::CountHits
int CountHits(const art::Event &evt, const art::InputTag &which, unsigned int cryostat, unsigned int tpc, unsigned int plane)
Definition: RecoBaseDrawer.cxx:4079
Exceptions.h
Collection of exceptions for Geometry system.
evd::RecoBaseDrawer::Seed2D
void Seed2D(const art::Event &evt, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:632
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
evd::RecoBaseDrawer::GetEndPoint2D
int GetEndPoint2D(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::EndPoint2D > &ep2d)
Definition: RecoBaseDrawer.cxx:3884
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
evd::RecoBaseDrawer::FillTQHistoDP
void FillTQHistoDP(const art::Event &evt, unsigned int plane, unsigned int wire, TH1F *histo, std::vector< double > &htau1, std::vector< double > &htau2, std::vector< double > &hitamplitudes, std::vector< double > &hpeaktimes, std::vector< int > &hstartT, std::vector< int > &hendT, std::vector< int > &hNMultiHit, std::vector< int > &hLocalHitIndex)
Definition: RecoBaseDrawer.cxx:4200
recob::Shower::OpenAngle
double OpenAngle() const
Definition: Shower.h:202
geo::SigType_t
enum geo::_plane_sigtype SigType_t
evd::RecoBaseDrawer::PFParticle3D
void PFParticle3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2024
mode
const char mode
Definition: noise_ana.cxx:20
evd::RecoBaseDrawer::Draw2DSlopeEndPoints
void Draw2DSlopeEndPoints(double xStart, double yStart, double xEnd, double yEnd, double slope, int color, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:1030
recob::Track::Vertex
Point_t const & Vertex() const
Definition: bj/lardataobj/RecoBase/Track.h:124
evd::RecoBaseDrawer::Hit2D
int Hit2D(const art::Event &evt, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:259
util::quantities::tick
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
detinfo::DetectorPropertiesData::ConvertXToTicks
double ConvertXToTicks(double X, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:99
evd::RecoBaseDrawer::GetSlices
int GetSlices(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Slice > &slices)
Definition: RecoBaseDrawer.cxx:3804
evd::RecoBaseDrawer::Slice3D
void Slice3D(const art::Event &evt, evdb::View3D *view)
Definition: RecoBaseDrawer.cxx:2968
evd::RecoBaseDrawer::DrawTrack2D
void DrawTrack2D(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< const recob::Hit * > &hits, evdb::View2D *view, unsigned int plane, const recob::Track *track, int color, int lineWidth)
Definition: RecoBaseDrawer.cxx:1302
evd::RecoBaseDrawer::fWireMax
std::vector< int > fWireMax
highest wire in interesting region for each plane
Definition: RecoBaseDrawer.h:335
evd::RecoBaseDrawer::GetHits
int GetHits(const art::Event &evt, const art::InputTag &which, std::vector< const recob::Hit * > &hits, unsigned int plane)
Definition: RecoBaseDrawer.cxx:3767
evd::RecoBaseDrawer::Circle3D
std::vector< std::array< double, 3 > > Circle3D(const TVector3 &pos, const TVector3 &axisDir, const double &radius)
Definition: RecoBaseDrawer.cxx:2828
evd::RecoBaseDrawer::GetWires
int GetWires(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Wire > &wires)
Definition: RecoBaseDrawer.cxx:3740
recob::SpacePoint::XYZ
const Double32_t * XYZ() const
Definition: SpacePoint.h:76
evd::RecoBaseDrawer::GetOpFlashes
int GetOpFlashes(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::OpFlash > &opflash)
Definition: RecoBaseDrawer.cxx:3911
evd::RecoBaseDrawer::GetSeeds
int GetSeeds(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Seed > &seed)
Definition: RecoBaseDrawer.cxx:3938
evd::RecoBaseDrawer::GetShowers
int GetShowers(const art::Event &evt, const art::InputTag &which, art::View< recob::Shower > &shower)
Definition: RecoBaseDrawer.cxx:4011
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
detinfo::DetectorPropertiesData::DriftVelocity
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Definition: DetectorPropertiesData.cc:23
evd::RecoBaseDrawer::fWireMin
std::vector< int > fWireMin
lowest wire in interesting region for each plane
Definition: RecoBaseDrawer.h:334
Cluster.h
Declaration of cluster object.
RecoBaseDrawer.h
Class to aid in the rendering of RecoBase objects.
evd::kColor
static const int kColor[kNCOLS]
Definition: eventdisplay.h:11
Track.h
Provides recob::Track data product.
evd::RecoBaseDrawer::GetChargeSum
void GetChargeSum(int plane, double &charge, double &convcharge)
Definition: RecoBaseDrawer.cxx:478
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
dropbox.dir
tuple dir
Definition: dropbox.py:28
WireGeo.h
Encapsulate the geometry of a wire.
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
tca::seeds
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:14
evd::RecoBaseDrawer::DrawTrack3D
void DrawTrack3D(const recob::Track &track, evdb::View3D *view, int color, int marker=1, float size=2.)
Definition: RecoBaseDrawer.cxx:2639
eventdisplay.h
Place to keep constants for event display.
recob::Shower::ID
int ID() const
Definition: Shower.h:187
PlaneGeo.h
Encapsulate the construction of a single detector plane.
Event.h
Represents the Event implemented as a self balancing binary search tree.
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
evd::RecoBaseDrawer::FillQHisto
void FillQHisto(const art::Event &evt, unsigned int plane, TH1F *histo)
Definition: RecoBaseDrawer.cxx:4158
detinfo::DetectorPropertiesData::BirksCorrection
double BirksCorrection(double dQdX) const
dQ/dX in electrons/cm, returns dE/dX in MeV/cm.
Definition: DetectorPropertiesData.cc:29
evd::RecoBaseDrawer::GetPFParticles
int GetPFParticles(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::PFParticle > &pfpart)
Definition: RecoBaseDrawer.cxx:3858
evd::RecoBaseDrawer::EndPoint2D
void EndPoint2D(const art::Event &evt, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:488
evd::RecoBaseDrawer::Wire2D
void Wire2D(const art::Event &evt, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:122
ChannelStatusProvider.h
Interface for experiment-specific channel quality info provider.
e
do i e
Definition: updateVersionList.sh:68
evd::RecoBaseDrawer::DrawPFParticleOrtho
void DrawPFParticleOrtho(const art::Ptr< recob::PFParticle > &pfPart, const art::PtrVector< recob::PFParticle > &pfParticleVec, const art::FindMany< recob::SpacePoint > &spacePointAssnsVec, const art::FindMany< recob::PCAxis > &pcAxisAssnVec, int depth, evd::OrthoProj_t proj, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3206
geo::InvalidWireError
Exception thrown on invalid wire number.
Definition: Exceptions.h:42
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
recob::Track::End
Point_t const & End() const
Definition: bj/lardataobj/RecoBase/Track.h:125
Wire.h
Declaration of basic channel signal object.
evd::RecoBaseDrawer::OpFlash2D
void OpFlash2D(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:543
evd::RecoBaseDrawer::GetRegionOfInterest
int GetRegionOfInterest(int plane, int &minw, int &maxw, int &mint, int &maxt)
Definition: RecoBaseDrawer.cxx:448
recob::Track::DirectionAtPoint
Vector_t DirectionAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:134
evd::RecoBaseDrawer::Event2D
void Event2D(const art::Event &evt, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:1835
evd::RecoBaseDrawer::GetClusterOutlines
void GetClusterOutlines(std::vector< const recob::Hit * > &hits, std::vector< double > &tpts, std::vector< double > &wpts, unsigned int plane)
Definition: RecoBaseDrawer.cxx:1165
evd::RecoBaseDrawer::OpFlashOrtho
void OpFlashOrtho(const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, evd::OrthoProj_t proj, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:2992
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
evd::RecoBaseDrawer::Vertex2D
void Vertex2D(const art::Event &evt, detinfo::DetectorPropertiesData const &detProp, evdb::View2D *view, unsigned int plane)
Definition: RecoBaseDrawer.cxx:1761
ChannelStatusService.h
Interface for experiment-specific service for channel quality info.
evd::RecoBaseDrawer::DrawTrackOrtho
void DrawTrackOrtho(const recob::Track &track, int color, evd::OrthoProj_t proj, double msize, evdb::View2D *view)
Definition: RecoBaseDrawer.cxx:3572
evd::kColor2
static const int kColor2[kNCOLS]
Definition: eventdisplay.h:13
geo::InvalidWireError::suggestedWireID
geo::WireID suggestedWireID() const
Returns a better wire ID.
Definition: Exceptions.h:120
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
geo::Point_t
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
evd::RecoBaseDrawer::fTimeMax
std::vector< int > fTimeMax
highest time in interesting region for each plane
Definition: RecoBaseDrawer.h:337
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
geo::TPCGeo::LocalToWorld
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
evd::RecoBaseDrawer::GetVertices
int GetVertices(const art::Event &evt, const art::InputTag &which, art::PtrVector< recob::Vertex > &vertex)
Definition: RecoBaseDrawer.cxx:4027
p1
physics associatedGroupsWithLeft p1
Definition: test_associatedgroupswithleft_full.fcl:24
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
evd::RecoBaseDrawer::VertexOrtho
void VertexOrtho(const art::PtrVector< recob::Vertex > &vertex, evd::OrthoProj_t proj, evdb::View2D *view, int marker)
Definition: RecoBaseDrawer.cxx:3069
evd::RecoBaseDrawer::GetEdges
int GetEdges(const art::Event &evt, const art::InputTag &which, std::vector< art::Ptr< recob::Edge >> &edges)
Definition: RecoBaseDrawer.cxx:3977
TPCGeo.h
Encapsulate the construction of a single detector plane.
evd::RecoBaseDrawer::fTimeMin
std::vector< int > fTimeMin
lowest time in interesting region for each plane
Definition: RecoBaseDrawer.h:336
Generated by   doxygen 1.8.5