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