LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DrawGausHits_tool.cc
Go to the documentation of this file.
1 
10 
13 
19 #include "cetlib_except/exception.h"
21 
22 #include "TF1.h"
23 #include "TPolyLine.h"
24 
25 #include <fstream>
26 
27 namespace evdb_tool {
28 
29  class DrawGausHits : public IWFHitDrawer {
30  public:
31  explicit DrawGausHits(const fhicl::ParameterSet& pset);
32 
33  ~DrawGausHits();
34 
35  void configure(const fhicl::ParameterSet& pset) override;
36  void Draw(evdb::View2D&, raw::ChannelID_t&) const override;
37 
38  private:
39  struct HitParams_t {
40  float hitCenter;
41  float hitSigma;
42  float hitHeight;
43  float hitStart;
44  float hitEnd;
45  };
46 
47  using ROIHitParamsVec = std::vector<HitParams_t>;
48  using HitParamsVec = std::vector<ROIHitParamsVec>;
49 
52  std::vector<int> fColorVec;
53  mutable std::vector<std::unique_ptr<TF1>> fHitFunctionVec;
54  };
55 
56  //----------------------------------------------------------------------
57  // Constructor.
59  {
60  configure(pset);
61  }
62 
64 
66  {
67  fNumPoints = pset.get<int>("NumPoints", 1000);
68  fFloatBaseline = pset.get<bool>("FloatBaseline", false);
69 
70  fColorVec.clear();
71  fColorVec.push_back(kRed);
72  fColorVec.push_back(kGreen);
73  fColorVec.push_back(kMagenta);
74  fColorVec.push_back(kCyan);
75 
76  fHitFunctionVec.clear();
77 
78  return;
79  }
80 
81  void DrawGausHits::Draw(evdb::View2D& view2D, raw::ChannelID_t& channel) const
82  {
84 
85  //grab the singleton with the event
87  if (!event) return;
88 
89  fHitFunctionVec.clear();
90 
91  for (size_t imod = 0; imod < recoOpt->fHitLabels.size(); ++imod) {
92  // Step one is to recover the hits for this label that match the input channel
93  art::InputTag const which = recoOpt->fHitLabels[imod];
94 
96  event->getByLabel(which, hitVecHandle);
97 
98  // Get a container for the subset of hits we are drawing
100 
101  try {
102  for (size_t hitIdx = 0; hitIdx < hitVecHandle->size(); hitIdx++) {
103  art::Ptr<recob::Hit> hit(hitVecHandle, hitIdx);
104 
105  if (hit->Channel() == channel) hitPtrVec.push_back(hit);
106  }
107  }
108  catch (cet::exception& e) {
109  mf::LogWarning("DrawGausHits") << "DrawGausHits"
110  << " failed with message:\n"
111  << e;
112  }
113 
114  if (hitPtrVec.empty()) continue;
115 
116  // Apparently you cannot trust some hit producers to put the hits in the correct order!
117  std::sort(hitPtrVec.begin(), hitPtrVec.end(), [](const auto& left, const auto& right) {
118  return left->PeakTime() < right->PeakTime();
119  });
120 
121  // Get associations to wires
122  art::FindManyP<recob::Wire> wireAssnsVec(hitPtrVec, *event, which);
123  std::vector<float> wireDataVec;
124 
125  // Recover the full (zero-padded outside ROI's) deconvolved waveform for this wire
126  if (wireAssnsVec.isValid() && wireAssnsVec.size() > 0) {
127  auto hwafp = wireAssnsVec.at(0).front();
128  if (!hwafp.isNull() && hwafp.isAvailable()) { wireDataVec = hwafp->Signal(); }
129  }
130 
131  // Now go through and process the hits back into the hit parameters
132  using ROIHitParamsVec = std::vector<HitParams_t>;
133  using HitParamsVec = std::vector<ROIHitParamsVec>;
134 
135  // Get an initial container for common hits on ROI
136  HitParamsVec hitParamsVec;
137  ROIHitParamsVec roiHitParamsVec;
138  raw::TDCtick_t lastEndTick(10000);
139 
140  for (const auto& hit : hitPtrVec) {
141  // check roi end condition
142  if (hit->PeakTime() - 3. * hit->RMS() > lastEndTick) {
143  if (!roiHitParamsVec.empty()) hitParamsVec.push_back(roiHitParamsVec);
144  roiHitParamsVec.clear();
145  }
146 
147  HitParams_t hitParams;
148 
149  hitParams.hitCenter = hit->PeakTime();
150  hitParams.hitSigma = hit->RMS();
151  hitParams.hitHeight = hit->PeakAmplitude();
152  hitParams.hitStart = hit->StartTick(); //PeakTime() - 3. * hit->RMS();
153  hitParams.hitEnd = hit->EndTick(); //PeakTime() + 3. * hit->RMS();
154 
155  lastEndTick = hitParams.hitEnd;
156 
157  roiHitParamsVec.emplace_back(hitParams);
158 
159  } //end loop over reco hits
160 
161  // Just in case (probably never called...)
162  if (!roiHitParamsVec.empty()) hitParamsVec.push_back(roiHitParamsVec);
163 
164  size_t roiCount(0);
165 
166  for (const auto& roiHitParamsVec : hitParamsVec) {
167  // Create a histogram here...
168  double roiStart = roiHitParamsVec.front().hitStart;
169  double roiStop = roiHitParamsVec.back().hitEnd;
170 
171  std::string funcString = "gaus(0)";
172  std::string funcName = Form("hitshape_%05zu_c%02zu", size_t(channel), roiCount++);
173 
174  for (size_t idx = 1; idx < roiHitParamsVec.size(); idx++)
175  funcString += "+gaus(" + std::to_string(3 * idx) + ")";
176 
177  // Include a baseline
178  float baseline(0.);
179 
180  if (fFloatBaseline && !wireDataVec.empty()) baseline = wireDataVec.at(roiStart);
181 
182  funcString += "+" + std::to_string(baseline);
183 
184  fHitFunctionVec.emplace_back(
185  std::make_unique<TF1>(funcName.c_str(), funcString.c_str(), roiStart, roiStop));
186  TF1& hitFunc = *fHitFunctionVec.back().get();
187 
188  hitFunc.SetLineColor(fColorVec[imod % fColorVec.size()]);
189 
190  size_t idx(0);
191  for (const auto& hitParams : roiHitParamsVec) {
192  hitFunc.SetParameter(idx + 0, hitParams.hitHeight);
193  hitFunc.SetParameter(idx + 1, hitParams.hitCenter);
194  hitFunc.SetParameter(idx + 2, hitParams.hitSigma);
195 
196  TPolyLine& hitHeight = view2D.AddPolyLine(2, kBlack, 1, 1);
197 
198  hitHeight.SetPoint(0, hitParams.hitCenter, baseline);
199  hitHeight.SetPoint(1, hitParams.hitCenter, hitParams.hitHeight + baseline);
200 
201  hitHeight.Draw("same");
202 
203  TPolyLine& hitSigma = view2D.AddPolyLine(2, kGray, 1, 1);
204 
205  hitSigma.SetPoint(
206  0, hitParams.hitCenter - hitParams.hitSigma, 0.6 * hitParams.hitHeight + baseline);
207  hitSigma.SetPoint(
208  1, hitParams.hitCenter + hitParams.hitSigma, 0.6 * hitParams.hitHeight + baseline);
209 
210  hitSigma.Draw("same");
211 
212  idx += 3;
213  }
214 
215  hitFunc.Draw("same");
216  }
217  } //end loop over HitFinding modules
218 
219  return;
220  }
221 
223 }
std::vector< HitParams_t > ROIHitParamsVec
void configure(const fhicl::ParameterSet &pset) override
const art::Event * GetEvent() const
Definition: EventHolder.cxx:45
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
TPolyLine & AddPolyLine(int n, int c, int w, int s)
Definition: View2D.cxx:210
iterator begin()
Definition: PtrVector.h:217
Declaration of signal hit object.
std::vector< ROIHitParamsVec > HitParamsVec
A collection of drawable 2-D objects.
Singleton to hold the current art::Event for the event display.
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
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.
static EventHolder * Instance()
Definition: EventHolder.cxx:15
T get(std::string const &key) const
Definition: ParameterSet.h:314
iterator end()
Definition: PtrVector.h:231
bool empty() const
Definition: PtrVector.h:330
Detector simulation of raw signals on wires.
std::vector< int > fColorVec
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Declaration of basic channel signal object.
void Draw(evdb::View2D &, raw::ChannelID_t &) const override
std::vector< art::InputTag > fHitLabels
module labels that produced hits
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
Float_t e
Definition: plot.C:35
std::vector< std::unique_ptr< TF1 > > fHitFunctionVec
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:268
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.
DrawGausHits(const fhicl::ParameterSet &pset)