LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
FFTHitFinder_module.cc
Go to the documentation of this file.
1 //
3 // FFTHitFinder class
4 //
5 // pagebri3@msu.edu
6 //
7 // This algorithm is designed to find hits on wires after deconvolution
8 // with an average shape used as the input response.
10 
11 // C/C++ standard library
12 #include <numeric> // std::accumulate
13 #include <string>
14 
15 // Framework includes
22 
23 // LArSoft Includes
25 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
29 
30 // ROOT Includes
31 #include "TDecompSVD.h"
32 #include "TF1.h"
33 #include "TH1D.h"
34 #include "TMath.h"
35 
36 namespace hit {
37 
38  class FFTHitFinder : public art::EDProducer {
39 
40  public:
41  explicit FFTHitFinder(fhicl::ParameterSet const& pset);
42 
43  private:
44  void produce(art::Event& evt) override;
45 
46  std::string fCalDataModuleLabel;
47  double fMinSigInd;
48  double fMinSigCol;
49  double fIndWidth;
50  double fColWidth;
51  double fIndMinWidth;
52  double fColMinWidth;
55  std::vector<double> fAreaNorms;
56 
57  }; // class FFTHitFinder
58 
59  //-------------------------------------------------
61  {
62  fCalDataModuleLabel = pset.get<std::string>("CalDataModuleLabel");
63  fMinSigInd = pset.get<double>("MinSigInd");
64  fMinSigCol = pset.get<double>("MinSigCol");
65  fIndWidth = pset.get<double>("IndWidth");
66  fColWidth = pset.get<double>("ColWidth");
67  fIndMinWidth = pset.get<double>("IndMinWidth");
68  fColMinWidth = pset.get<double>("ColMinWidth");
69  fMaxMultiHit = pset.get<int>("MaxMultiHit");
70  fAreaMethod = pset.get<int>("AreaMethod");
71  fAreaNorms = pset.get<std::vector<double>>("AreaNorms");
72 
73  // let HitCollectionCreator declare that we are going to produce
74  // hits and associations with wires and raw digits
75  // (with no particular product label)
77  }
78 
79  // This algorithm uses the fact that deconvolved signals are very smooth
80  // and looks for hits as areas between local minima that have signal above
81  // threshold.
82  //-------------------------------------------------
84  {
85  // this object contains the hit collection
86  // and its associations to wires and raw digits:
88 
89  // Read in the wire List object(s).
91  evt.getByLabel(fCalDataModuleLabel, wireVecHandle);
92  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
93 
94  // also get the raw digits associated with wires
95  art::FindOneP<raw::RawDigit> WireToRawDigits(wireVecHandle, evt, fCalDataModuleLabel);
96 
97  std::vector<int> startTimes; // stores time of 1st local minimum
98  std::vector<int> maxTimes; // stores time of local maximum
99  std::vector<int> endTimes; // stores time of 2nd local minimum
100  int time = 0; // current time bin
101  int minTimeHolder = 0; // current start time
102  raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
103  bool maxFound = false; // Flag for whether a peak > threshold has been found
104  double threshold = 0.; // minimum signal size for id'ing a hit
105  double fitWidth = 0.; // hit fit width initial value
106  double minWidth = 0.; // minimum hit width
107  std::string eqn = "gaus(0)"; // string for equation to fit
108  geo::SigType_t sigType = geo::kInduction; // type of plane we are looking at
109  std::stringstream numConv;
110 
111  //loop over wires
112  for (unsigned int wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
113 
114  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
115  startTimes.clear();
116  maxTimes.clear();
117  endTimes.clear();
118  std::vector<float> signal(wire->Signal());
119  std::vector<float>::iterator timeIter; // iterator for time bins
120  time = 0;
121  minTimeHolder = 0;
122  maxFound = false;
123  channel = wire->Channel();
124  sigType = wireReadoutGeom.SignalType(channel);
125 
126  //Set the appropriate signal widths and thresholds
127  if (sigType == geo::kInduction) {
128  threshold = fMinSigInd;
129  fitWidth = fIndWidth;
130  minWidth = fIndMinWidth;
131  }
132  else if (sigType == geo::kCollection) {
133  threshold = fMinSigCol;
134  fitWidth = fColWidth;
135  minWidth = fColMinWidth;
136  }
137  // loop over signal
138  for (timeIter = signal.begin(); timeIter + 2 < signal.end(); timeIter++) {
139  //test if timeIter+1 is a local minimum
140  if (*timeIter > *(timeIter + 1) && *(timeIter + 1) < *(timeIter + 2)) {
141  //only add points if already found a local max above threshold.
142  if (maxFound) {
143  endTimes.push_back(time + 1);
144  maxFound = false;
145  //keep these in case new hit starts right away
146  minTimeHolder = time + 2;
147  }
148  else
149  minTimeHolder = time + 1;
150  }
151  //if not a minimum, test if we are at a local maximum
152  //if so, and the max value is above threshold, add it and proceed.
153  else if (*timeIter < *(timeIter + 1) && *(timeIter + 1) > *(timeIter + 2) &&
154  *(timeIter + 1) > threshold) {
155  maxFound = true;
156  maxTimes.push_back(time + 1);
157  startTimes.push_back(minTimeHolder);
158  }
159  time++;
160  } //end loop over signal vec
161 
162  //if no inflection found before end, but peak found add end point
163  while (maxTimes.size() > endTimes.size())
164  endTimes.push_back(signal.size() - 1);
165  if (startTimes.size() == 0) continue;
166 
167  //All code below does the fitting, adding of hits
168  //to the hit vector and when all wires are complete
169  //saving them
170  double totSig(0); // stores the total hit signal
171  double startT(0); // stores the start time
172  double endT(0); // stores the end time
173  int numHits(0); // number of consecutive hits being fitted
174  int size(0); // size of data vector for fit
175  int hitIndex(0); // index of current hit in sequence
176  double amplitude(0), position(0), width(0); //fit parameters
177  double amplitudeErr(0), positionErr(0), widthErr(0); //fit errors
178  double goodnessOfFit(0), chargeErr(0); //Chi2/NDF and error on charge
179  double minPeakHeight(0); //lowest peak height in multi-hit fit
180 
181  //stores gaussian paramters first index is the hit number
182  //the second refers to height, position, and width respectively
183  std::vector<double> hitSig;
184 
185  //add found hits to hit vector
186  while (hitIndex < (signed)startTimes.size()) {
187 
188  startT = endT = 0;
189  numHits = 1;
190  minPeakHeight = signal[maxTimes[hitIndex]];
191 
192  //consider adding pulse to group of consecutive hits if:
193  //1 less than max consecutive hits
194  //2 we are not at the last point in the signal vector
195  //3 the height of the dip between the two is greater than threshold/2
196  //4 and there is no gap between them
197  while (numHits < fMaxMultiHit && numHits + hitIndex < (signed)endTimes.size() &&
198  signal[endTimes[hitIndex + numHits - 1]] > threshold / 2.0 &&
199  startTimes[hitIndex + numHits] - endTimes[hitIndex + numHits - 1] < 2) {
200 
201  if (signal[maxTimes[hitIndex + numHits]] < minPeakHeight)
202  minPeakHeight = signal[maxTimes[hitIndex + numHits]];
203 
204  ++numHits;
205  }
206 
207  //finds the first point > 1/2 the smallest peak
208  startT = startTimes[hitIndex];
209 
210  while (signal[(int)startT] < minPeakHeight / 2.0)
211  ++startT;
212 
213  //finds the first point from the end > 1/2 the smallest peak
214  endT = endTimes[hitIndex + numHits - 1];
215 
216  while (signal[(int)endT] < minPeakHeight / 2.0)
217  --endT;
218  size = (int)(endT - startT);
219  TH1D hitSignal("hitSignal", "", size, startT, endT);
220  for (int i = (int)startT; i < (int)endT; ++i)
221  hitSignal.Fill(i, signal[i]);
222 
223  //build the TFormula
224  eqn = "gaus(0)";
225 
226  for (int i = 3; i < numHits * 3; i += 3) {
227  eqn.append("+gaus(");
228  numConv.str("");
229  numConv << i;
230  eqn.append(numConv.str());
231  eqn.append(")");
232  }
233 
234  TF1 gSum("gSum", eqn.c_str(), 0, size);
235 
236  if (numHits > 1) {
237  TArrayD data(numHits * numHits);
238  TVectorD amps(numHits);
239  for (int i = 0; i < numHits; ++i) {
240  amps[i] = signal[maxTimes[hitIndex + i]];
241  for (int j = 0; j < numHits; j++)
242  data[i + numHits * j] =
243  TMath::Gaus(maxTimes[hitIndex + j], maxTimes[hitIndex + i], fitWidth);
244  } //end loop over hits
245 
246  //This section uses a linear approximation in order to get an
247  //initial value of the individual hit amplitudes
248  try {
249  TMatrixD h(numHits, numHits);
250  h.Use(numHits, numHits, data.GetArray());
251  TDecompSVD a(h);
252  a.Solve(amps);
253  }
254  catch (...) {
255  mf::LogInfo("FFTHitFinder") << "TDcompSVD failed";
256  hitIndex += numHits;
257  continue;
258  }
259 
260  for (int i = 0; i < numHits; ++i) {
261  //if the approximation makes a peak vanish
262  //set initial height as average of threshold and
263  //raw peak height
264  if (amps[i] > 0)
265  amplitude = amps[i];
266  else
267  amplitude = 0.5 * (threshold + signal[maxTimes[hitIndex + i]]);
268  gSum.SetParameter(3 * i, amplitude);
269  gSum.SetParameter(1 + 3 * i, maxTimes[hitIndex + i]);
270  gSum.SetParameter(2 + 3 * i, fitWidth);
271  gSum.SetParLimits(3 * i, 0.0, 3.0 * amplitude);
272  gSum.SetParLimits(1 + 3 * i, startT, endT);
273  gSum.SetParLimits(2 + 3 * i, 0.0, 10.0 * fitWidth);
274  } //end loop over hits
275  } //end if numHits > 1
276  else {
277  gSum.SetParameters(signal[maxTimes[hitIndex]], maxTimes[hitIndex], fitWidth);
278  gSum.SetParLimits(0, 0.0, 1.5 * signal[maxTimes[hitIndex]]);
279  gSum.SetParLimits(1, startT, endT);
280  gSum.SetParLimits(2, 0.0, 10.0 * fitWidth);
281  }
282 
284  hitSignal.Fit(&gSum, "QNRW", "", startT, endT);
285  for (int hitNumber = 0; hitNumber < numHits; ++hitNumber) {
286  totSig = 0;
287  if (gSum.GetParameter(3 * hitNumber) > threshold / 2.0 &&
288  gSum.GetParameter(3 * hitNumber + 2) > minWidth) {
289  amplitude = gSum.GetParameter(3 * hitNumber);
290  position = gSum.GetParameter(3 * hitNumber + 1);
291  width = gSum.GetParameter(3 * hitNumber + 2);
292  amplitudeErr = gSum.GetParError(3 * hitNumber);
293  positionErr = gSum.GetParError(3 * hitNumber + 1);
294  widthErr = gSum.GetParError(3 * hitNumber + 2);
295  goodnessOfFit = gSum.GetChisquare() / (double)gSum.GetNDF();
296  int DoF = gSum.GetNDF();
297 
298  //estimate error from area of Gaussian
299  chargeErr = std::sqrt(TMath::Pi()) * (amplitudeErr * width + widthErr * amplitude);
300 
301  hitSig.resize(size);
302 
303  for (int sigPos = 0; sigPos < size; ++sigPos) {
304  hitSig[sigPos] = amplitude * TMath::Gaus(sigPos + startT, position, width);
305  totSig += hitSig[(int)sigPos];
306  }
307 
308  if (fAreaMethod)
309  totSig = std::sqrt(2 * TMath::Pi()) * amplitude * width / fAreaNorms[(size_t)sigType];
310 
311  // get the WireID for this hit
312  std::vector<geo::WireID> wids = wireReadoutGeom.ChannelToWire(channel);
314  // for now, just take the first option returned from ChannelToWire
315  geo::WireID wid = wids[0];
316 
317  // make the hit
318  recob::HitCreator hit(*wire, // wire
319  wid, // wireID
320  (int)startT, // start_tick
321  (int)endT, // end_tick
322  width, // rms
323  position, // peak_time
324  positionErr, // sigma_peak_time
325  amplitude, // peak_amplitude
326  amplitudeErr, // sigma_peak_amplitude
327  totSig, // hit_integral
328  chargeErr, // hit_sigma_integral
329  std::accumulate // summedADC
330  (signal.begin() + (int)startT, signal.begin() + (int)endT, 0.),
331  1, // multiplicity
332  -1, // local_index
334  goodnessOfFit, // goodness_of_fit
335  DoF // dof
336  );
337 
338  // get the object associated with the original hit
339  art::Ptr<raw::RawDigit> rawdigits = WireToRawDigits.at(wireIter);
340 
341  hcol.emplace_back(hit.move(), wire, rawdigits);
342  } //end if over threshold
343  } //end loop over hits
344  hitIndex += numHits;
345  } // end while on hitIndex<(signed)startTimes.size()
346 
347  } // while on Wires
348 
349  // put the hit collection and associations into the event
350  hcol.put_into(evt);
351 
352  } // End of produce()
353 
355 
356 } // end of hit namespace
double fIndMinWidth
Minimum induction hit width.
intermediate_table::iterator iterator
FFTHitFinder(fhicl::ParameterSet const &pset)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
double fMinSigInd
Induction signal height threshold.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
double fIndWidth
Initial width for induction fit.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
int fMaxMultiHit
maximum hits for multi fit
double fColWidth
Initial width for collection fit.
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:261
std::vector< double > fAreaNorms
factors for converting area to same units as peak height
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
Helper functions to create a hit.
std::string fCalDataModuleLabel
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:31
A class handling a collection of hits and its associations.
Definition: HitCreator.h:489
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Signal from induction planes.
Definition: geo_types.h:147
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:223
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
double fColMinWidth
Minimum collection hit width.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:299
void produce(art::Event &evt) override
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:622
Detector simulation of raw signals on wires.
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
Definition: Wire.cxx:30
ProducesCollector & producesCollector() noexcept
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
double fMinSigCol
Collection signal height threshold.
Declaration of basic channel signal object.
int fAreaMethod
Type of area calculation.
TCEvent evt
Definition: DataStructs.cxx:8
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
Definition: fwd.h:26
Signal from collection planes.
Definition: geo_types.h:148