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