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