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