LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
CCHitFinderAlg.cxx
Go to the documentation of this file.
1 
13 // class header
15 
16 // C/C++ standard libraries
17 #include <algorithm> // std::sort(), std::copy()
18 #include <array>
19 #include <cmath> // std::sqrt(), std::abs()
20 #include <iomanip>
21 #include <iostream>
22 #include <sstream>
23 #include <utility> // std::pair<>, std::make_pair()
24 
25 // framework libraries
27 
28 // LArSoft Includes
30 #include "lardata/Utilities/SimpleFits.h" // lar::util::GaussianFit<>
31 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
32 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
33 
34 // ROOT Includes
35 #include "TF1.h"
36 #include "TGraph.h"
37 
38 namespace hit {
39 
40  constexpr unsigned int CCHitFinderAlg::MaxGaussians; // definition
41 
42  //------------------------------------------------------------------------------
44  : FitCache{"GausFitCache_CCHitFinderAlg"}
45  {
46  if (pset.has_key("MinSigInd"))
48  << "CCHitFinderAlg: Using no-longer-valid fcl input: MinSigInd, MinSigCol, etc";
49 
50  fMinPeak = pset.get<std::vector<float>>("MinPeak");
51  fMinRMS = pset.get<std::vector<float>>("MinRMS");
52  fMaxBumps = pset.get<unsigned short>("MaxBumps");
53  fMaxXtraHits = pset.get<unsigned short>("MaxXtraHits");
54  fChiSplit = pset.get<float>("ChiSplit");
55  fChiNorms = pset.get<std::vector<float>>("ChiNorms");
56  fUseFastFit = pset.get<bool>("UseFastFit", false);
57  fUseChannelFilter = pset.get<bool>("UseChannelFilter", true);
58  fStudyHits = pset.get<bool>("StudyHits", false);
59  // The following variables are only used in StudyHits mode
60  fUWireRange = pset.get<std::vector<short>>("UWireRange");
61  fUTickRange = pset.get<std::vector<short>>("UTickRange");
62  fVWireRange = pset.get<std::vector<short>>("VWireRange");
63  fVTickRange = pset.get<std::vector<short>>("VTickRange");
64  fWWireRange = pset.get<std::vector<short>>("WWireRange");
65  fWTickRange = pset.get<std::vector<short>>("WTickRange");
66 
67  if (fMinPeak.size() != fMinRMS.size()) {
68  mf::LogError("CCTF") << "MinPeak size != MinRMS size";
69  return;
70  }
71 
72  if (fMaxBumps > MaxGaussians) {
73  // MF_LOG_WARNING will point the user to this line of code.
74  // Any value of MaxGaussians can be used.
75  // That value is defined in the header file.
76  MF_LOG_WARNING("CCHitFinderAlg")
77  << "CCHitFinder algorithm is currently hard-coded to support at most " << MaxGaussians
78  << " bumps per region of interest, but " << fMaxBumps << " have been requested.\n"
79  << "We are forcing the parameter to " << MaxGaussians
80  << ". If this is not acceptable, increase CCHitFinderAlg::MaxGaussians"
81  << " value and recompile.";
82  fMaxBumps = MaxGaussians;
83  } // if too many gaussians
84 
87 
88  // sanity check for StudyHits mode
89  if (fStudyHits) {
90  if (fUWireRange.size() != 2 || fUTickRange.size() != 2 || fVWireRange.size() != 2 ||
91  fVTickRange.size() != 2 || fWWireRange.size() != 2 || fWTickRange.size() != 2) {
92  mf::LogError("CCHF") << "Invalid vector size for StudyHits. Must be 2";
93  }
94  } // fStudyHits
95  }
96 
97  //------------------------------------------------------------------------------
98  void CCHitFinderAlg::RunCCHitFinder(std::vector<recob::Wire> const& Wires)
99  {
100  allhits.clear();
101 
102  constexpr unsigned short maxticks = 1000;
103  float* ticks = new float[maxticks];
104  // define the ticks array used for fitting
105  for (unsigned short ii = 0; ii < maxticks; ++ii) {
106  ticks[ii] = ii;
107  }
108  float* signl = new float[maxticks];
109  float adcsum = 0;
110  // initialize the vectors for the hit study
111  if (fStudyHits) StudyHits(0);
112  bool first;
113 
114  lariov::ChannelStatusProvider const& channelStatus =
116 
117  for (size_t wireIter = 0; wireIter < Wires.size(); wireIter++) {
118  recob::Wire const& theWire = Wires[wireIter];
119  theChannel = theWire.Channel();
120  // ignore bad channels
121  if (channelStatus.IsBad(theChannel)) continue;
122 
123  std::vector<geo::WireID> wids = wireReadoutGeom->ChannelToWire(theChannel);
124  thePlane = wids[0].Plane;
125  if (thePlane > fMinPeak.size() - 1) {
126  mf::LogError("CCHF") << "MinPeak vector too small for plane " << thePlane;
127  return;
128  }
129  theWireNum = wids[0].Wire;
130  HitChannelInfo_t WireInfo{&theWire, wids[0], wireReadoutGeom->SignalType(theWire.Channel())};
131 
132  // minimum number of time samples
133  unsigned short minSamples = 2 * fMinRMS[thePlane];
134 
135  // factor used to normalize the chi/dof fits for each plane
137 
138  std::vector<float> signal(theWire.Signal());
139 
140  unsigned short nabove = 0;
141  unsigned short tstart = 0;
142  unsigned short maxtime = signal.size() - 2;
143  // find the min time when the signal is below threshold
144  unsigned short mintime = 3;
145  for (unsigned short time = 3; time < maxtime; ++time) {
146  if (signal[time] < fMinPeak[thePlane]) {
147  mintime = time;
148  break;
149  }
150  }
151  for (unsigned short time = mintime; time < maxtime; ++time) {
152  if (signal[time] > fMinPeak[thePlane]) {
153  if (nabove == 0) tstart = time;
154  ++nabove;
155  }
156  else {
157  // check for a wide enough signal above threshold
158  if (nabove > minSamples) {
159  // skip this wire if the RAT is too long
160  if (nabove > maxticks)
161  mf::LogError("CCHitFinder")
162  << "Long RAT " << nabove << " " << maxticks << " No signal on wire " << theWireNum
163  << " after time " << time;
164  if (nabove > maxticks) break;
165  unsigned short npt = 0;
166  // look for bumps to inform the fit
167  bumps.clear();
168  adcsum = 0;
169  for (unsigned short ii = tstart; ii < time; ++ii) {
170  signl[npt] = signal[ii];
171  adcsum += signl[npt];
172  if (signal[ii] > signal[ii - 1] && signal[ii - 1] > signal[ii - 2] &&
173  signal[ii] > signal[ii + 1] && signal[ii + 1] > signal[ii + 2])
174  bumps.push_back(npt);
175  ++npt;
176  }
177  // decide if this RAT should be studied
178  if (fStudyHits) StudyHits(1, npt, signl, tstart);
179  // just make a crude hit if too many bumps
180  if (bumps.size() > fMaxBumps) {
181  MakeCrudeHit(npt, ticks, signl);
182  StoreHits(tstart, npt, WireInfo, adcsum);
183  nabove = 0;
184  continue;
185  }
186  // start looking for hits with the found bumps
187  unsigned short nHitsFit = bumps.size();
188  unsigned short nfit = 0;
189  chidof = 0.;
190  dof = -1;
191  bool HitStored = false;
192  unsigned short nMaxFit = bumps.size() + fMaxXtraHits;
193  // only used in StudyHits mode
194  first = true;
195  while (nHitsFit <= nMaxFit) {
196 
197  FitNG(nHitsFit, npt, ticks, signl);
198  if (fStudyHits && first && SelRAT) {
199  first = false;
200  StudyHits(2, npt, signl, tstart);
201  }
202  // good chisq so store it
203  if (chidof < fChiSplit) {
204  StoreHits(tstart, npt, WireInfo, adcsum);
205  HitStored = true;
206  break;
207  }
208  // the previous fit was better, so revert to it and store it
209  ++nHitsFit;
210  ++nfit;
211  } // nHitsFit < fMaxXtraHits
212  if (!HitStored && npt < maxticks) {
213  // failed all fitting. Make a crude hit
214  MakeCrudeHit(npt, ticks, signl);
215  StoreHits(tstart, npt, WireInfo, adcsum);
216  }
217  else if (nHitsFit > 0)
218  FinalFitStats.AddMultiGaus(nHitsFit);
219  } // nabove > minSamples
220  nabove = 0;
221  } // signal < fMinPeak
222  } // time
223  } // wireIter
224 
225  // print out
226  if (fStudyHits) StudyHits(4);
227 
228  delete[] ticks;
229  delete[] signl;
230 
231  } //RunCCHitFinder
232 
234  bool CCHitFinderAlg::FastGaussianFit(unsigned short npt,
235  float const* ticks,
236  float const* signl,
237  std::array<double, 3>& params,
238  std::array<double, 3>& paramerrors,
239  float& chidof)
240  {
241  // parameters: amplitude, mean, sigma
242 
243  lar::util::GaussianFit<double> fitter; // probably "double" is overdoing
244 
245  // apply a time shift so that the center of the time interval is 0
246  const float time_shift = (ticks[npt - 1] + ticks[0]) / 2.F;
247 
248  // fill the input data, no uncertainty
249  for (size_t i = 0; i < npt; ++i) {
250  if (signl[i] <= 0) {
251  MF_LOG_DEBUG("CCHitFinderAlg")
252  << "Non-positive charge encountered. Backing up to ROOT fit.";
253  return false;
254  }
255  // we could freely add a Poisson uncertainty (as third parameter)
256  fitter.add(ticks[i] - time_shift, signl[i]);
257  } // for
258 
259  // we might have found that we don't want the fast fit after all...
260  if (!fitter.FillResults(params, paramerrors)) {
261  // something went wrong...
262  MF_LOG_DEBUG("CCHitFinderAlg") << "Fast Gaussian fit failed.";
263  return false;
264  }
265 
266  // note that this is not the full chi^2, but it is the chi^2 of the
267  // parabolic fit underlying the Gaussian one
268  const double chi2 = fitter.ChiSquare();
269  chidof = chi2 / fitter.NDF();
270 
271  // remove the time shift
272  params[1] += time_shift; // mean
273 
274  // GP: inflate the uncertainties on the fit parameters according to chi2/NDF
275  // (not sure if this is in any way justified)
276  if (chidof > 1.)
277  for (double& par : paramerrors)
278  par *= std::sqrt(chidof);
279 
280  return true;
281  } // FastGaussianFit()
282 
284  void CCHitFinderAlg::FitNG(unsigned short nGaus, unsigned short npt, float* ticks, float* signl)
285  {
286  // Fit the signal to n Gaussians
287 
288  dof = npt - 3 * nGaus;
289 
290  chidof = 9999.;
291 
292  if (dof < 3) return;
293  if (bumps.size() == 0) return;
294 
295  // load the fit into a temp vector
296  std::vector<double> partmp;
297  std::vector<double> partmperr;
298 
299  //
300  // if it is possible, we try first with the quick single Gaussian fit
301  //
303 
304  bool bNeedROOTfit = (nGaus > 1) || !fUseFastFit;
305  if (!bNeedROOTfit) {
306  // so, we need only one puny Gaussian;
307  std::array<double, 3> params, paramerrors;
308 
310 
311  if (FastGaussianFit(npt, ticks, signl, params, paramerrors, chidof)) {
312  // success? copy the results in the proper structures
313  partmp.resize(3);
314  std::copy(params.begin(), params.end(), partmp.begin());
315  partmperr.resize(3);
316  std::copy(paramerrors.begin(), paramerrors.end(), partmperr.begin());
317  }
318  else
319  bNeedROOTfit = true; // if we fail, let's schedule ROOT to back us up
320 
321  if (!bNeedROOTfit) FinalFitStats.AddFast();
322 
323  } // if we don't need ROOT to fit
324 
325  if (bNeedROOTfit) {
326  // we may land here either because the simple Gaussian fit did not work (either
327  // failed, or we chose not to trust it) or because the fit is multi-Gaussian
328 
329  // define the fit string to pass to TF1
330 
331  std::stringstream numConv;
332  std::string eqn = "gaus";
333  if (nGaus > 1) eqn = "gaus(0)";
334  for (unsigned short ii = 3; ii < nGaus * 3; ii += 3) {
335  eqn.append(" + gaus(");
336  numConv.str("");
337  numConv << ii;
338  eqn.append(numConv.str());
339  eqn.append(")");
340  }
341 
342  auto Gn = std::make_unique<TF1>("gn", eqn.c_str());
343  TGraph* fitn = new TGraph(npt, ticks, signl);
344 
345  // put in the bump parameters. Assume that nGaus >= bumps.size()
346  for (unsigned short ii = 0; ii < bumps.size(); ++ii) {
347  unsigned short index = ii * 3;
348  unsigned short bumptime = bumps[ii];
349  double amp = signl[bumptime];
350  Gn->SetParameter(index, amp);
351  Gn->SetParLimits(index, 0., 9999.);
352  Gn->SetParameter(index + 1, (double)bumptime);
353  Gn->SetParLimits(index + 1, 0, (double)npt);
354  Gn->SetParameter(index + 2, (double)fMinRMS[thePlane]);
355  Gn->SetParLimits(index + 2, 1., 3 * (double)fMinRMS[thePlane]);
356  } // ii bumps
357 
358  // search for other bumps that may be hidden by the already found ones
359  for (unsigned short ii = bumps.size(); ii < nGaus; ++ii) {
360  // bump height must exceed fMinPeak
361  float big = fMinPeak[thePlane];
362  unsigned short imbig = 0;
363  for (unsigned short jj = 0; jj < npt; ++jj) {
364  float diff = signl[jj] - Gn->Eval((Double_t)jj, 0, 0, 0);
365  if (diff > big) {
366  big = diff;
367  imbig = jj;
368  }
369  } // jj
370  if (imbig > 0) {
371  // set the parameters for the bump
372  unsigned short index = ii * 3;
373  Gn->SetParameter(index, (double)big);
374  Gn->SetParLimits(index, 0., 9999.);
375  Gn->SetParameter(index + 1, (double)imbig);
376  Gn->SetParLimits(index + 1, 0, (double)npt);
377  Gn->SetParameter(index + 2, (double)fMinRMS[thePlane]);
378  Gn->SetParLimits(index + 2, 1., 5 * (double)fMinRMS[thePlane]);
379  } // imbig > 0
380  } // ii
381 
382  // W = set weights to 1, N = no drawing or storing, Q = quiet
383  // B = bounded parameters
384  fitn->Fit(&*Gn, "WNQB");
385 
386  for (unsigned short ipar = 0; ipar < 3 * nGaus; ++ipar) {
387  partmp.push_back(Gn->GetParameter(ipar));
388  partmperr.push_back(Gn->GetParError(ipar));
389  }
390  chidof = Gn->GetChisquare() / (dof * chinorm);
391 
392  delete fitn;
393  // delete Gn;
394 
395  } // if ROOT fit
396 
397  // Sort by increasing time if necessary
398  if (nGaus > 1) {
399  std::vector<std::pair<unsigned short, unsigned short>> times;
400  // fill the sort vector
401  for (unsigned short ii = 0; ii < nGaus; ++ii) {
402  unsigned short index = ii * 3;
403  times.push_back(std::make_pair(partmp[index + 1], ii));
404  } // ii
405  std::sort(times.begin(), times.end());
406  // see if re-arranging is necessary
407  bool sortem = false;
408  for (unsigned short ii = 0; ii < nGaus; ++ii) {
409  if (times[ii].second != ii) {
410  sortem = true;
411  break;
412  }
413  } // ii
414  if (sortem) {
415  // temp temp vectors for putting things in the right time order
416  std::vector<double> partmpt;
417  std::vector<double> partmperrt;
418  for (unsigned short ii = 0; ii < nGaus; ++ii) {
419  unsigned short index = times[ii].second * 3;
420  partmpt.push_back(partmp[index]);
421  partmpt.push_back(partmp[index + 1]);
422  partmpt.push_back(partmp[index + 2]);
423  partmperrt.push_back(partmperr[index]);
424  partmperrt.push_back(partmperr[index + 1]);
425  partmperrt.push_back(partmperr[index + 2]);
426  } // ii
427  partmp = partmpt;
428  partmperr = partmperrt;
429  } // sortem
430  } // nGaus > 1
431 
432  // ensure that the fit is reasonable
433  bool fitok = true;
434  for (unsigned short ii = 0; ii < nGaus; ++ii) {
435  unsigned short index = ii * 3;
436  // ensure that the fitted time is within the signal bounds
437  short fittime = partmp[index + 1];
438  if (fittime < 0 || fittime > npt - 1) {
439  fitok = false;
440  break;
441  }
442  // ensure that the signal peak is large enough
443  if (partmp[index] < fMinPeak[thePlane]) {
444  fitok = false;
445  break;
446  }
447  // ensure that the RMS is large enough but not too large
448  float rms = partmp[index + 2];
449  if (rms < 0.5 * fMinRMS[thePlane] || rms > 5 * fMinRMS[thePlane]) {
450  fitok = false;
451  break;
452  }
453  // ensure that the hits are not too similar in time (< 2 ticks)
454  for (unsigned short jj = 0; jj < nGaus; ++jj) {
455  if (jj == ii) continue;
456  unsigned short jndex = jj * 3;
457  float timediff = std::abs(partmp[jndex + 1] - partmp[index + 1]);
458  if (timediff < 2.) {
459  fitok = false;
460  break;
461  }
462  }
463  if (!fitok) break;
464  }
465 
466  if (fitok) {
467  par = partmp;
468  parerr = partmperr;
469  }
470  else {
471  chidof = 9999.;
472  dof = -1;
473  }
474  } // FitNG
475 
477  void CCHitFinderAlg::MakeCrudeHit(unsigned short npt, float* ticks, float* signl)
478  {
479  // make a single crude hit if fitting failed
480  float sumS = 0.;
481  float sumST = 0.;
482  for (unsigned short ii = 0; ii < npt; ++ii) {
483  sumS += signl[ii];
484  sumST += signl[ii] * ticks[ii];
485  }
486  float mean = sumST / sumS;
487  float rms = 0.;
488  for (unsigned short ii = 0; ii < npt; ++ii) {
489  float arg = ticks[ii] - mean;
490  rms += signl[ii] * arg * arg;
491  }
492  rms = std::sqrt(rms / sumS);
493  float amp = sumS / (Sqrt2Pi * rms);
494  par.clear();
495  par.push_back(amp);
496  par.push_back(mean);
497  par.push_back(rms);
498  // need to do the errors better
499  parerr.clear();
500  float amperr = npt;
501  float meanerr = std::sqrt(1 / sumS);
502  float rmserr = 0.2 * rms;
503  parerr.push_back(amperr);
504  parerr.push_back(meanerr);
505  parerr.push_back(rmserr);
506  chidof = 9999.;
507  dof = -1;
508  } // MakeCrudeHit
509 
511  void CCHitFinderAlg::StoreHits(unsigned short TStart,
512  unsigned short npt,
513  HitChannelInfo_t info,
514  float adcsum)
515  {
516  // store the hits in the struct
517  size_t nhits = par.size() / 3;
518 
519  if (allhits.max_size() - allhits.size() < nhits) {
520  mf::LogError("CCHitFinder") << "Too many hits: existing " << allhits.size() << " plus new "
521  << nhits << " beyond the maximum " << allhits.max_size();
522  return;
523  }
524 
525  if (nhits == 0) return;
526 
527  // fill RMS for single hits
528  if (fStudyHits) StudyHits(3);
529 
530  const float loTime = TStart;
531  const float hiTime = TStart + npt;
532 
533  // Find sum of the areas of all Gaussians
534  float gsum = 0.;
535  for (size_t hit = 0; hit < nhits; ++hit) {
536  const unsigned short index = 3 * hit;
537  gsum += Sqrt2Pi * par[index] * par[index + 2];
538  }
539  for (size_t hit = 0; hit < nhits; ++hit) {
540  const size_t index = 3 * hit;
541  const float charge = Sqrt2Pi * par[index] * par[index + 2];
542  const float charge_err =
543  SqrtPi * (parerr[index] * par[index + 2] + par[index] * parerr[index + 2]);
544 
545  allhits.emplace_back(info.wire->Channel(), // channel
546  loTime, // start_tick
547  hiTime, // end_tick
548  par[index + 1] + TStart, // peak_time
549  parerr[index + 1], // sigma_peak_time
550  par[index + 2], // rms
551  par[index], // peak_amplitude
552  parerr[index], // sigma_peak_amplitude
553  adcsum * charge / gsum, // ROIsummedADC
554  adcsum * charge / gsum, // HitsummedADC NOT CORRECTLY FILLED
555  charge, // hit_integral
556  charge_err, // hit_sigma_integral
557  nhits, // multiplicity
558  hit, // local_index
559  chidof, // goodness_of_fit
560  dof, // dof
561  info.wire->View(), // view
562  info.sigType, // signal_type
563  info.wireID // wireID
564  );
565  } // hit
566  } // StoreHits
567 
569  void CCHitFinderAlg::StudyHits(unsigned short flag,
570  unsigned short npt,
571  float* signl,
572  unsigned short tstart)
573  {
574  // study hits in user-selected ranges of wires and ticks in each plane. The user should identify
575  // a shallow-angle isolated track, e.g. using the event display, to determine the wire/tick ranges.
576  // One hit should be reconstructed on each wire when the hit finding fcl parameters are set correctly.
577  // The intent of this study is to determine the correct fcl parameters. The flag variable determines
578  // the operation performed.
579  // flag = 0: Initialize study vectors
580  // flag = 1: Set SelRat true if the Region Above Threshold resides within a wire/hit range
581  // flag = 2: Find the maximum signal and calculate the RMS. Also find the low and high ticks of signals
582  // in the wire range to allow a later calculation of the track angle. This isn't strictly
583  // necessary for the study and presumes that the user has selected compatible regions in each plane.
584  // flag = 3: Accumulate the RMS from the first Gaussian fit
585  // flag = 4: Calculate recommended fcl parameters and print the results to the screen
586 
587  // init
588  if (flag == 0) {
589  for (unsigned short ipl = 0; ipl < 3; ++ipl) {
590  // Average chisq of the first fit on a single bump in each plane
591  bumpChi.push_back(0.);
592  // Average RMS of the dump
593  bumpRMS.push_back(0.);
594  // The number of single bumps in each plane
595  bumpCnt.push_back(0.);
596  // number of RATs
597  RATCnt.push_back(0);
598  // The number of single hits found in each plane
599  hitCnt.push_back(0.);
600  // Average reconstructed hit RMS
601  hitRMS.push_back(0.);
602  // lo/hi wire/time
603  loWire.push_back(9999.);
604  loTime.push_back(0.);
605  hiWire.push_back(-1.);
606  hiTime.push_back(0.);
607  } // ii
608  return;
609  } // flag == 0
610 
611  if (flag == 1) {
612  SelRAT = false;
613  if (thePlane == 0) {
614  if (theWireNum > fUWireRange[0] && theWireNum < fUWireRange[1] && tstart > fUTickRange[0] &&
615  tstart < fUTickRange[1]) {
616  SelRAT = true;
617  RATCnt[thePlane] += 1;
618  }
619  return;
620  } // thePlane == 0
621  if (thePlane == 1) {
622  if (theWireNum > fVWireRange[0] && theWireNum < fVWireRange[1] && tstart > fVTickRange[0] &&
623  tstart < fVTickRange[1]) {
624  SelRAT = true;
625  RATCnt[thePlane] += 1;
626  }
627  return;
628  } // thePlane == 1
629  if (thePlane == 2) {
630  if (theWireNum > fWWireRange[0] && theWireNum < fWWireRange[1] && tstart > fWTickRange[0] &&
631  tstart < fWTickRange[1]) {
632  SelRAT = true;
633  RATCnt[thePlane] += 1;
634  }
635  return;
636  } // thePlane == 2
637  } // flag == 1
638 
639  if (flag == 2) {
640  if (!SelRAT) return;
641  // in this section we find the low/hi wire/time for a signal. This can be used to calculate
642  // the slope dT/dW to study hit width, fraction of crude hits, etc vs dT/dW
643  float big = 0.;
644  float imbig = 0.;
645  for (unsigned short ii = 0; ii < npt; ++ii) {
646  if (signl[ii] > big) {
647  big = signl[ii];
648  imbig = ii;
649  }
650  } // ii
651  // require a significant PH
652  if (big > fMinPeak[0]) {
653  // get the Lo info
654  if (theWireNum < loWire[thePlane]) {
656  loTime[thePlane] = tstart + imbig;
657  }
658  // get the Hi info
659  if (theWireNum > hiWire[thePlane]) {
661  hiTime[thePlane] = tstart + imbig;
662  }
663  } // big > fMinPeak[0]
664  if (bumps.size() == 1 && chidof < 9999.) {
665  bumpCnt[thePlane] += bumps.size();
666  bumpChi[thePlane] += chidof;
667  // calculate the average bin
668  float sumt = 0.;
669  float sum = 0.;
670  for (unsigned short ii = 0; ii < npt; ++ii) {
671  sum += signl[ii];
672  sumt += signl[ii] * ii;
673  } // ii
674  float aveb = sumt / sum;
675  // now calculate the RMS
676  sumt = 0.;
677  for (unsigned short ii = 0; ii < npt; ++ii) {
678  float dbin = (float)ii - aveb;
679  sumt += signl[ii] * dbin * dbin;
680  } // ii
681  bumpRMS[thePlane] += std::sqrt(sumt / sum);
682  } // bumps.size() == 1 && chidof < 9999.
683  return;
684  } // flag == 2
685 
686  // fill info for single hits
687  if (flag == 3) {
688  if (!SelRAT) return;
689  if (par.size() == 3) {
690  hitCnt[thePlane] += 1;
691  hitRMS[thePlane] += par[2];
692  }
693  return;
694  }
695 
696  if (flag == 4) {
697  // The goal is to adjust the fcl inputs so that the number of single
698  // hits found is ~equal to the number of single bumps found for shallow
699  // angle tracks. The ChiNorm inputs should be adjusted so the average
700  // chisq/DOF is ~1 in each plane.
701  std::cout << "Check lo and hi W/T for each plane" << std::endl;
702  for (unsigned short ipl = 0; ipl < 3; ++ipl) {
703  std::cout << ipl << " lo " << loWire[ipl] << " " << loTime[ipl] << " hi " << hiWire[ipl]
704  << " " << hiTime[ipl] << std::endl;
705  }
706  std::cout << " ipl nRAT bCnt bChi bRMS hCnt hRMS dT/dW New_ChiNorm" << std::endl;
707  for (unsigned short ipl = 0; ipl < 3; ++ipl) {
708  if (bumpCnt[ipl] > 0) {
709  bumpChi[ipl] = bumpChi[ipl] / (float)bumpCnt[ipl];
710  bumpRMS[ipl] = bumpRMS[ipl] / (float)bumpCnt[ipl];
711  hitRMS[ipl] = hitRMS[ipl] / (float)hitCnt[ipl];
712  // calculate the slope
713  float dTdW = std::abs((hiTime[ipl] - loTime[ipl]) / (hiWire[ipl] - loWire[ipl]));
714  std::cout << ipl << std::right << std::setw(5) << RATCnt[ipl] << std::setw(5)
715  << bumpCnt[ipl] << std::setw(7) << std::fixed << std::setprecision(2)
716  << bumpChi[ipl] << std::setw(7) << bumpRMS[ipl] << std::setw(7) << hitCnt[ipl]
717  << std::setw(7) << std::setprecision(1) << hitRMS[ipl] << std::setw(7) << dTdW
718  << std::setw(7) << std::setprecision(2) << bumpChi[ipl] * fChiNorms[ipl]
719  << std::endl;
720  } //
721  } // ipl
722  std::cout << "nRAT is the number of Regions Above Threshold (RAT) used in the study.\n";
723  std::cout << "bCnt is the number of single bumps that were successfully fitted \n";
724  std::cout << "bChi is the average chisq/DOF of the first fit\n";
725  std::cout << "bRMS is the average calculated RMS of the bumps\n";
726  std::cout << "hCnt is the number of RATs that have a single hit\n";
727  std::cout << "hRMS is the average RMS from the Gaussian fit -> use this value for "
728  "fMinRMS[plane] in the fcl file\n";
729  std::cout << "dTdW is the slope of the track\n";
730  std::cout
731  << "New_ChiNorm is the recommended values of ChiNorm that should be used in the fcl file\n";
732  bumpChi.clear();
733  bumpRMS.clear();
734  bumpCnt.clear();
735  RATCnt.clear();
736  hitRMS.clear();
737  hitCnt.clear();
738  loWire.clear();
739  loTime.clear();
740  hiWire.clear();
741  hiTime.clear();
742  }
743  } // StudyHits
744 
746  void CCHitFinderAlg::FitStats_t::Reset(unsigned int nGaus)
747  {
748  if (nGaus == 0) return;
749  MultiGausFits.resize(nGaus);
750  std::fill(MultiGausFits.begin(), MultiGausFits.end(), 0);
751  FastFits = 0;
752  }
753 
755  {
756  ++MultiGausFits[std::min(nGaus, (unsigned int)MultiGausFits.size()) - 1];
757  }
758 
759 } // namespace hit
unsigned short fMaxBumps
void RunCCHitFinder(std::vector< recob::Wire > const &Wires)
std::vector< short > fWTickRange
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
bool add(Data_t x, Data_t y, Data_t sy=Data_t(1.0))
Clears all the input statistics.
Definition: SimpleFits.h:1719
void AddMultiGaus(unsigned int nGaus)
float fChiSplit
Estimated noise error on the Signal.
static constexpr float SqrtPi
void FitNG(unsigned short nGaus, unsigned short npt, float *ticks, float *signl)
raw::ChannelID_t theChannel
std::vector< float > loWire
void StoreHits(unsigned short TStart, unsigned short npt, HitChannelInfo_t info, float adcsum)
std::vector< float > hiWire
std::vector< unsigned short > bumps
std::vector< float > fChiNorms
std::vector< float > bumpRMS
bool fUseFastFit
whether to attempt using a fast fit on single gauss.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CCHitFinderAlg(fhicl::ParameterSet const &pset)
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
FitStats_t TriedFitStats
counts of the tried fits
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
tick ticks
Alias for common language habits.
Definition: electronics.h:76
Classes performing simple fits.
geo::View_t View() const
Returns the view the channel belongs to.
Definition: Wire.h:219
virtual bool FillResults(FitParameters_t &params, FitMatrix_t &Xmat, Data_t &det, FitMatrix_t &Smat) const override
Fills the specified parameters.
Definition: SimpleFits.h:1803
std::vector< short > fUTickRange
std::vector< recob::Hit > allhits
unsigned short thePlane
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:223
geo::WireReadoutGeom const * wireReadoutGeom
void Reset(unsigned int nGaus)
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
"Fast" Gaussian fit
Definition: SimpleFits.h:993
std::vector< float > hiTime
std::vector< short > fUWireRange
virtual int NDF() const override
Returns the degrees of freedom in the determination of the fit.
Definition: SimpleFits.h:1144
std::vector< int > bumpCnt
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
std::vector< float > loTime
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
static bool FastGaussianFit(unsigned short npt, float const *ticks, float const *signl, std::array< double, 3 > &params, std::array< double, 3 > &paramerrors, float &chidof)
Performs a "fast" fit.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
void StudyHits(unsigned short flag, unsigned short npt=0, float *signl=nullptr, unsigned short tstart=0)
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
std::vector< short > fWWireRange
unsigned short theWireNum
std::vector< float > fMinPeak
#define MF_LOG_DEBUG(id)
Class holding the regions of interest of signal from a channel.
Definition: Wire.h:116
std::vector< int > hitCnt
static constexpr unsigned int MaxGaussians
static constexpr float Sqrt2Pi
std::vector< short > fVWireRange
std::vector< short > fVTickRange
virtual Data_t ChiSquare() const override
Returns the of the original fit.
Definition: SimpleFits.h:1135
unsigned short fMaxXtraHits
void MakeCrudeHit(unsigned short npt, float *ticks, float *signl)
std::vector< float > hitRMS
#define MF_LOG_WARNING(category)
Double_t sum
Definition: plot.C:31
std::vector< int > RATCnt
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
FitStats_t FinalFitStats
counts of the good fits
std::vector< double > parerr
std::vector< float > fMinRMS
std::vector< float > bumpChi
exchange data about the originating wire
std::vector< double > par
art framework interface to geometry description