LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
DPRawHitFinder_module.cc
Go to the documentation of this file.
1 //
3 // DPRawHitFinder class
4 //
5 // christoph.alt@cern.ch
6 //
7 // This algorithm is designed to find hits on raw waveforms in collection planes (dual phase/single phase)
8 // It is based on GausHitFinder.
9 // -----------------------------------
10 //
11 // 1. The algorithm walks along the wire and looks for pulses above threshold.
12 // 2. Depending on distance and minimum ADC count between peaks inside the same ROI,
13 // peaks can be grouped. Grouped peaks are fitted together (see later).
14 // 3. Inside each group of peaks, look for pairs of peaks that are very close, with
15 // one peak being much smaller than the other one (in integral and amplitude).
16 // If such a pair of peaks is found, merge the two peaks to one peak.
17 //
18 // For pulse trains with #peaks <= MaxMultiHit and width < MaxGroupLength:
19 // 4. Fit n double exponentials to each group of peaks, where n is the number
20 // of peaks inside this group.
21 // 5. If the Chi2/NDF returned > Chi2NDFRetry, attempt to fit n+1 double exponentials
22 // to the group of peaks by adding a peak close to the maximum deviation between
23 // fit and waveform. If this is a better fit it then uses the parameters of this
24 // fit to characterize the "hit" object. If not, try n+2 exponentials and so on.
25 // Stop when Chi2/NDF is good or 3 times the number of inital exponentials is reached.
26 //
27 // If Chi2/NDF is still bad or if #peaks > MaxMultiHit or width > MaxGroupLength:
28 // 6. Split pulse into equally long hits.
29 //
30 // The parameters of the fit are saved in a feature vector by using MVAWriter to
31 // draw the fitted function in the event display.
32 //
34 
35 // C/C++ standard library
36 #include <algorithm> // std::accumulate()
37 #include <cmath>
38 #include <memory> // std::unique_ptr()
39 #include <string>
40 #include <utility> // std::move()
41 
42 // Framework includes
43 #include "art/Framework/Core/EDProducer.h"
44 #include "art/Framework/Core/ModuleMacros.h"
45 #include "art/Framework/Principal/Event.h"
46 #include "art/Framework/Services/Registry/ServiceHandle.h"
47 #include "art/Framework/Services/System/TriggerNamesService.h"
48 #include "art_root_io/TFileService.h"
49 #include "canvas/Persistency/Common/FindOneP.h"
50 #include "canvas/Utilities/InputTag.h"
51 #include "fhiclcpp/ParameterSet.h"
52 #include "messagefacility/MessageLogger/MessageLogger.h"
53 
54 // LArSoft Includes
55 #include "larcore/Geometry/WireReadout.h"
56 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
57 #include "lardata/ArtDataHelper/HitCreator.h"
58 #include "lardata/ArtDataHelper/MVAWriter.h"
59 #include "lardataobj/RecoBase/Hit.h"
60 #include "lardataobj/RecoBase/Wire.h"
61 
62 // ROOT Includes
63 #include "TF1.h"
64 #include "TH1F.h"
65 #include "TMath.h"
66 
67 namespace hit {
68  class DPRawHitFinder : public art::EDProducer {
69 
70  public:
71  explicit DPRawHitFinder(fhicl::ParameterSet const& pset);
72 
73  private:
74  void produce(art::Event& evt) override;
75  void beginJob() override;
76 
77  using TimeValsVec = std::vector<std::tuple<int, int, int>>; // start, max, end of a peak
78  using PeakTimeWidVec = std::vector<
79  std::tuple<int, int, int, int>>; // max, width, start, end of a peak within a group
80  using MergedTimeWidVec =
81  std::vector<std::tuple<int,
82  int,
83  PeakTimeWidVec,
84  int>>; // start, end of group of peaks, PeakTimeWidVec, NFluctuations
85  using PeakDevVec = std::vector<std::tuple<double, int, int, int>>;
86  using ParameterVec = std::vector<std::pair<double, double>>; //< parameter/error vec
87 
88  void findCandidatePeaks(std::vector<float>::const_iterator startItr,
89  std::vector<float>::const_iterator stopItr,
90  TimeValsVec& timeValsVec,
91  float& PeakMin,
92  int firstTick) const;
93 
94  int EstimateFluctuations(const std::vector<float> fsignalVec,
95  int peakStart,
96  int peakMean,
97  int peakEnd);
98 
99  void mergeCandidatePeaks(const std::vector<float> signalVec, TimeValsVec, MergedTimeWidVec&);
100 
101  // ### This function will fit N-Exponentials to a TH1D where N is set ###
102  // ### by the number of peaks found in the pulse ###
103 
104  void FitExponentials(const std::vector<float> fSignalVector,
105  const PeakTimeWidVec fPeakVals,
106  int fStartTime,
107  int fEndTime,
108  ParameterVec& fparamVec,
109  double& fchi2PerNDF,
110  int& fNDF,
111  bool fSameShape);
112 
113  void FindPeakWithMaxDeviation(const std::vector<float> fSignalVector,
114  int fNPeaks,
115  int fStartTime,
116  int fEndTime,
117  bool fSameShape,
118  ParameterVec fparamVec,
119  PeakTimeWidVec fpeakVals,
120  PeakDevVec& fPeakDev);
121 
122  std::string CreateFitFunction(int fNPeaks, bool fSameShape);
123 
124  void AddPeak(std::tuple<double, int, int, int> fPeakDevCand, PeakTimeWidVec& fpeakValsTemp);
125 
126  void SplitPeak(std::tuple<double, int, int, int> fPeakDevCand, PeakTimeWidVec& fpeakValsTemp);
127 
128  double WidthFunc(double fPeakMean,
129  double fPeakAmp,
130  double fPeakTau1,
131  double fPeakTau2,
132  double fStartTime,
133  double fEndTime,
134  double fPeakMeanTrue);
135 
136  double ChargeFunc(double fPeakMean,
137  double fPeakAmp,
138  double fPeakTau1,
139  double fPeakTau2,
140  double fChargeNormFactor,
141  double fPeakMeanTrue);
142 
143  void FillOutHitParameterVector(const std::vector<double>& input, std::vector<double>& output);
144 
145  void doBinAverage(const std::vector<float>& inputVec,
146  std::vector<float>& outputVec,
147  size_t binsToAverage) const;
148 
149  void reBin(const std::vector<float>& inputVec,
150  std::vector<float>& outputVec,
151  size_t nBinsToCombine) const;
152 
153  struct Comp {
154  // important: we need two overloads, because the comparison
155  // needs to be done both ways to check for equality
156 
157  bool operator()(std::tuple<int, int, int, int> p, int s) const { return std::get<0>(p) < s; }
158  bool operator()(int s, std::tuple<int, int, int, int> p) const { return s < std::get<0>(p); }
159  };
160 
161  std::string fCalDataModuleLabel;
162 
163  //FHiCL parameter (see "hitfindermodules.fcl" for details)
164  int fLogLevel;
165  float fMinSig;
166  int fTicksToStopPeakFinder;
167  int fMinWidth;
168  double fMinADCSum;
169  double fMinADCSumOverWidth;
170  int fMaxMultiHit;
171  int fMaxGroupLength;
172  double fChargeNorm;
173  bool fTryNplus1Fits;
174  double fChi2NDFRetry;
175  double fChi2NDFRetryFactorMultiHits;
176  double fChi2NDFMax;
177  double fChi2NDFMaxFactorMultiHits;
178  size_t fNumBinsToAverage;
179  double fMinTau;
180  double fMaxTau;
181  double fFitPeakMeanRange;
182  int fGroupMaxDistance;
183  double fGroupMinADC;
184  bool fSameShape;
185  bool fDoMergePeaks;
186  double fMergeADCSumThreshold;
187  double fMergeMaxADCThreshold;
188  double fMinRelativePeakHeightLeft;
189  double fMinRelativePeakHeightRight;
190  double fMergeMaxADCLimit;
191  double fWidthNormalization;
192  int fLongMaxHits;
193  int fLongPulseWidth;
194  int fMaxFluctuations;
195 
196  art::InputTag
197  fNewHitsTag; // tag of hits produced by this module, need to have it for fit parameter data products
198  anab::FVectorWriter<4> fHitParamWriter; // helper for saving hit fit parameters in data products
199 
200  TH1F* fFirstChi2;
201  TH1F* fChi2;
202 
203  }; // class DPRawHitFinder
204 
205  //-------------------------------------------------
206  //-------------------------------------------------
207  DPRawHitFinder::DPRawHitFinder(fhicl::ParameterSet const& pset)
208  : EDProducer{pset}
209  , fNewHitsTag(pset.get<std::string>("module_label"),
210  "",
211  art::ServiceHandle<art::TriggerNamesService const>()->getProcessName())
212  , fHitParamWriter(producesCollector())
213  {
214  fLogLevel = pset.get<int>("LogLevel");
215  fCalDataModuleLabel = pset.get<std::string>("CalDataModuleLabel");
216  fMaxMultiHit = pset.get<int>("MaxMultiHit");
217  fMaxGroupLength = pset.get<int>("MaxGroupLength");
218  fTryNplus1Fits = pset.get<bool>("TryNplus1Fits");
219  fChi2NDFRetry = pset.get<double>("Chi2NDFRetry");
220  fChi2NDFRetryFactorMultiHits = pset.get<double>("Chi2NDFRetryFactorMultiHits");
221  fChi2NDFMax = pset.get<double>("Chi2NDFMax");
222  fChi2NDFMaxFactorMultiHits = pset.get<double>("Chi2NDFMaxFactorMultiHits");
223  fNumBinsToAverage = pset.get<size_t>("NumBinsToAverage");
224  fMinSig = pset.get<float>("MinSig");
225  fMinWidth = pset.get<int>("MinWidth");
226  fMinADCSum = pset.get<double>("MinADCSum");
227  fMinADCSumOverWidth = pset.get<double>("MinADCSumOverWidth");
228  fChargeNorm = pset.get<double>("ChargeNorm");
229  fTicksToStopPeakFinder = pset.get<double>("TicksToStopPeakFinder");
230  fMinTau = pset.get<double>("MinTau");
231  fMaxTau = pset.get<double>("MaxTau");
232  fFitPeakMeanRange = pset.get<double>("FitPeakMeanRange");
233  fGroupMaxDistance = pset.get<int>("GroupMaxDistance");
234  fGroupMinADC = pset.get<int>("GroupMinADC");
235  fSameShape = pset.get<bool>("SameShape");
236  fDoMergePeaks = pset.get<bool>("DoMergePeaks");
237  fMergeADCSumThreshold = pset.get<double>("MergeADCSumThreshold");
238  fMergeMaxADCThreshold = pset.get<double>("MergeMaxADCThreshold");
239  fMinRelativePeakHeightLeft = pset.get<double>("MinRelativePeakHeightLeft");
240  fMinRelativePeakHeightRight = pset.get<double>("MinRelativePeakHeightRight");
241  fMergeMaxADCLimit = pset.get<double>("MergeMaxADCLimit");
242  fWidthNormalization = pset.get<double>("WidthNormalization");
243  fLongMaxHits = pset.get<double>("LongMaxHits");
244  fLongPulseWidth = pset.get<double>("LongPulseWidth");
245  fMaxFluctuations = pset.get<double>("MaxFluctuations");
246 
247  // let HitCollectionCreator declare that we are going to produce
248  // hits and associations with wires and raw digits
249  // (with no particular product label)
250  recob::HitCollectionCreator::declare_products(producesCollector());
251 
252  // declare that data products with feature vectors describing
253  // hits is going to be produced
254  fHitParamWriter.produces_using<recob::Hit>();
255 
256  } // DPRawHitFinder::DPRawHitFinder()
257 
258  //-------------------------------------------------
259  //-------------------------------------------------
260  void DPRawHitFinder::FillOutHitParameterVector(const std::vector<double>& input,
261  std::vector<double>& output)
262  {
263  if (input.size() == 0)
264  throw std::runtime_error(
265  "DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
266 
267  const unsigned int N_PLANES = art::ServiceHandle<geo::WireReadout>()->Get().Nplanes();
268 
269  if (input.size() == 1)
270  output.resize(N_PLANES, input[0]);
271  else if (input.size() == N_PLANES)
272  output = input;
273  else
274  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config "
275  "vector size !=1 and !=N_PLANES.");
276  }
277 
278  //-------------------------------------------------
279  //-------------------------------------------------
280  void DPRawHitFinder::beginJob()
281  {
282  // get access to the TFile service
283  art::ServiceHandle<art::TFileService const> tfs;
284 
285  // ======================================
286  // === Hit Information for Histograms ===
287  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 5000);
288  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 5000);
289  }
290 
291  //-------------------------------------------------
292  void DPRawHitFinder::produce(art::Event& evt)
293  {
294  //==================================================================================================
295  TH1::AddDirectory(kFALSE);
296 
297  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
298 
299  // ###############################################
300  // ### Making a ptr vector to put on the event ###
301  // ###############################################
302  // this contains the hit collection
303  // and its associations to wires and raw digits
304  recob::HitCollectionCreator hcol(evt);
305 
306  // start collection of fit parameters, initialize metadata describing it
307  auto hitID =
308  fHitParamWriter.initOutputs<recob::Hit>(fNewHitsTag, {"t0", "tau1", "tau2", "ampl"});
309 
310  // ##########################################
311  // ### Reading in the Wire List object(s) ###
312  // ##########################################
313  art::Handle<std::vector<recob::Wire>> wireVecHandle;
314  evt.getByLabel(fCalDataModuleLabel, wireVecHandle);
315 
316  // #################################################################
317  // ### Reading in the RawDigit associated with these wires, too ###
318  // #################################################################
319  art::FindOneP<raw::RawDigit> RawDigits(wireVecHandle, evt, fCalDataModuleLabel);
320  // Channel Number
321  raw::ChannelID_t channel = raw::InvalidChannelID;
322 
323  //##############################
324  //### Looping over the wires ###
325  //##############################
326  for (size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
327  // ####################################
328  // ### Getting this particular wire ###
329  // ####################################
330  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
331  art::Ptr<raw::RawDigit> rawdigits = RawDigits.at(wireIter);
332  // --- Setting Channel Number and Signal type ---
333  channel = wire->Channel();
334  // get the WireID for this hit
335  std::vector<geo::WireID> wids = wireReadoutGeom.ChannelToWire(channel);
336  // for now, just take the first option returned from ChannelToWire
337  geo::WireID wid = wids[0];
338 
339  if (fLogLevel >= 1) {
340  std::cout << std::endl;
341  std::cout << std::endl;
342  std::cout << std::endl;
343  std::cout << "-----------------------------------------------------------------------------"
344  "------------------------------"
345  << std::endl;
346  std::cout << "Channel: " << channel << std::endl;
347  std::cout << "Cryostat: " << wid.Cryostat << std::endl;
348  std::cout << "TPC: " << wid.TPC << std::endl;
349  std::cout << "Plane: " << wid.Plane << std::endl;
350  std::cout << "Wire: " << wid.Wire << std::endl;
351  }
352 
353  // #################################################
354  // ### Set up to loop over ROI's for this wire ###
355  // #################################################
356  const recob::Wire::RegionsOfInterest_t& signalROI = wire->SignalROI();
357 
358  int CountROI = 0;
359 
360  for (const auto& range : signalROI.get_ranges()) {
361  // #################################################
362  // ### Getting a vector of signals for this wire ###
363  // #################################################
364  const std::vector<float>& signal = range.data();
365 
366  // ROI start time
367  raw::TDCtick_t roiFirstBinTick = range.begin_index();
368  MergedTimeWidVec mergedVec;
369 
370  // ###########################################################
371  // ### If option set do bin averaging before finding peaks ###
372  // ###########################################################
373 
374  if (fNumBinsToAverage > 1) {
375  std::vector<float> timeAve;
376  doBinAverage(signal, timeAve, fNumBinsToAverage);
377 
378  // ###################################################################
379  // ### Search current averaged ROI for candidate peaks and widths ###
380  // ###################################################################
381  TimeValsVec timeValsVec;
382  findCandidatePeaks(timeAve.begin(), timeAve.end(), timeValsVec, fMinSig, 0);
383 
384  // ####################################################
385  // ### If no startTime hit was found skip this wire ###
386  // ####################################################
387  if (timeValsVec.empty()) continue;
388 
389  // #############################################################
390  // ### Merge potentially overlapping peaks and do multi fit ###
391  // #############################################################
392  mergeCandidatePeaks(timeAve, timeValsVec, mergedVec);
393  }
394 
395  // ###########################################################
396  // ### Otherwise, operate directonly on signal vector ###
397  // ###########################################################
398  else {
399  // ##########################################################
400  // ### Search current ROI for candidate peaks and widths ###
401  // ##########################################################
402  TimeValsVec timeValsVec;
403  findCandidatePeaks(signal.begin(), signal.end(), timeValsVec, fMinSig, 0);
404 
405  if (fLogLevel >= 1) {
406  std::cout << std::endl;
407  std::cout << std::endl;
408  std::cout << "-------------------- ROI #" << CountROI << " -------------------- "
409  << std::endl;
410  if (timeValsVec.size() == 1)
411  std::cout << "ROI #" << CountROI << " (" << timeValsVec.size()
412  << " peak): ROIStartTick: " << range.offset
413  << " ROIEndTick:" << range.offset + range.size() << std::endl;
414  else
415  std::cout << "ROI #" << CountROI << " (" << timeValsVec.size()
416  << " peaks): ROIStartTick: " << range.offset
417  << " ROIEndTick:" << range.offset + range.size() << std::endl;
418  CountROI++;
419  }
420 
421  if (fLogLevel >= 2) {
422  int CountPeak = 0;
423  for (auto const& timeValsTmp : timeValsVec) {
424  std::cout << "Peak #" << CountPeak
425  << ": PeakStartTick: " << range.offset + std::get<0>(timeValsTmp)
426  << " PeakMaxTick: " << range.offset + std::get<1>(timeValsTmp)
427  << " PeakEndTick: " << range.offset + std::get<2>(timeValsTmp)
428  << std::endl;
429  CountPeak++;
430  }
431  }
432  // ####################################################
433  // ### If no startTime hit was found skip this wire ###
434  // ####################################################
435  if (timeValsVec.empty()) continue;
436 
437  // #############################################################
438  // ### Merge potentially overlapping peaks and do multi fit ###
439  // #############################################################
440  mergeCandidatePeaks(signal, timeValsVec, mergedVec);
441  }
442 
443  // #######################################################
444  // ### Creating the parameter vector for the new pulse ###
445  // #######################################################
446  ParameterVec paramVec;
447 
448  // === Number of Exponentials to try ===
449  int NumberOfPeaksBeforeFit = 0;
450  unsigned int nExponentialsForFit = 0;
451  double chi2PerNDF = 0.;
452  int NDF = 0;
453 
454  unsigned int NumberOfMergedVecs = mergedVec.size();
455 
456  // ################################################################
457  // ### Lets loop over the groups of peaks we found on this wire ###
458  // ################################################################
459 
460  for (unsigned int j = 0; j < NumberOfMergedVecs; j++) {
461  int startT = std::get<0>(mergedVec.at(j));
462  int endT = std::get<1>(mergedVec.at(j));
463  int width = endT + 1 - startT;
464  PeakTimeWidVec& peakVals = std::get<2>(mergedVec.at(j));
465 
466  int NFluctuations = std::get<3>(mergedVec.at(j));
467 
468  if (fLogLevel >= 3) {
469  std::cout << std::endl;
470  if (peakVals.size() == 1)
471  std::cout << "- Group #" << j << " (" << peakVals.size()
472  << " peak): GroupStartTick: " << range.offset + startT
473  << " GroupEndTick: " << range.offset + endT << std::endl;
474  else
475  std::cout << "- Group #" << j << " (" << peakVals.size()
476  << " peaks): GroupStartTick: " << range.offset + startT
477  << " GroupEndTick: " << range.offset + endT << std::endl;
478  std::cout << "Fluctuations in this group: " << NFluctuations << std::endl;
479  int CountPeakInGroup = 0;
480  for (auto const& peakValsTmp : peakVals) {
481  std::cout << "Peak #" << CountPeakInGroup << " in group #" << j
482  << ": PeakInGroupStartTick: " << range.offset + std::get<2>(peakValsTmp)
483  << " PeakInGroupMaxTick: " << range.offset + std::get<0>(peakValsTmp)
484  << " PeakInGroupEndTick: " << range.offset + std::get<3>(peakValsTmp)
485  << std::endl;
486  CountPeakInGroup++;
487  }
488  }
489 
490  // ### Getting rid of noise hits ###
491  if (width < fMinWidth ||
492  (double)std::accumulate(signal.begin() + startT, signal.begin() + endT + 1, 0) <
493  fMinADCSum ||
494  (double)std::accumulate(signal.begin() + startT, signal.begin() + endT + 1, 0) /
495  width <
496  fMinADCSumOverWidth) {
497  if (fLogLevel >= 3) {
498  std::cout << "Delete this group of peaks because width, integral or width/intergral "
499  "is too small."
500  << std::endl;
501  }
502  continue;
503  }
504 
505  // #####################################################################################################
506  // ### Only attempt to fit if number of peaks <= fMaxMultiHit and if group length <= fMaxGroupLength ###
507  // #####################################################################################################
508  NumberOfPeaksBeforeFit = peakVals.size();
509  nExponentialsForFit = peakVals.size();
510  chi2PerNDF = 0.;
511  NDF = 0;
512  if (NumberOfPeaksBeforeFit <= fMaxMultiHit && width <= fMaxGroupLength &&
513  NFluctuations <= fMaxFluctuations) {
514  // #####################################################
515  // ### Calling the function for fitting Exponentials ###
516  // #####################################################
517  paramVec.clear();
518  FitExponentials(signal, peakVals, startT, endT, paramVec, chi2PerNDF, NDF, fSameShape);
519 
520  if (fLogLevel >= 4) {
521  std::cout << std::endl;
522  std::cout << "--- First fit ---" << std::endl;
523  if (nExponentialsForFit == 1)
524  std::cout << "- Fitted " << nExponentialsForFit << " peak in group #" << j << ":"
525  << std::endl;
526  else
527  std::cout << "- Fitted " << nExponentialsForFit << " peaks in group #" << j << ":"
528  << std::endl;
529  std::cout << "chi2/ndf = " << chi2PerNDF << std::endl;
530 
531  if (fSameShape) {
532  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
533  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
534 
535  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
536  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2 * (i + 1)].first
537  << std::endl;
538  std::cout << "Peak #" << i
539  << ": t0 [ticks] = " << range.offset + paramVec[2 * (i + 1) + 1].first
540  << std::endl;
541  }
542  }
543  else {
544  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
545  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[4 * i + 2].first
546  << std::endl;
547  std::cout << "Peak #" << i
548  << ": t0 [ticks] = " << range.offset + paramVec[4 * i + 3].first
549  << std::endl;
550  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4 * i].first
551  << std::endl;
552  std::cout << "Peak #" << i << ": tau2 [mus] = " << paramVec[4 * i + 1].first
553  << std::endl;
554  }
555  }
556  }
557 
558  // If the chi2 is infinite then there is a real problem so we bail
559  if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) continue;
560 
561  fFirstChi2->Fill(chi2PerNDF);
562 
563  // ########################################################
564  // ### Trying extra Exponentials for an initial bad fit ###
565  // ########################################################
566 
567  if ((fTryNplus1Fits && nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFRetry) ||
568  (fTryNplus1Fits && nExponentialsForFit > 1 &&
569  chi2PerNDF > fChi2NDFRetryFactorMultiHits * fChi2NDFRetry)) {
570  unsigned int nExponentialsBeforeRefit = nExponentialsForFit;
571  unsigned int nExponentialsAfterRefit = nExponentialsForFit;
572  double oldChi2PerNDF = chi2PerNDF;
573  double chi2PerNDF2;
574  int NDF2;
575  bool RefitSuccess;
576  PeakTimeWidVec peakValsTemp;
577  while ((nExponentialsForFit == 1 &&
578  nExponentialsAfterRefit < 3 * nExponentialsBeforeRefit &&
579  chi2PerNDF > fChi2NDFRetry) ||
580  (nExponentialsForFit > 1 &&
581  nExponentialsAfterRefit < 3 * nExponentialsBeforeRefit &&
582  chi2PerNDF > fChi2NDFRetryFactorMultiHits * fChi2NDFRetry)) {
583  RefitSuccess = false;
584  PeakDevVec PeakDev;
585  FindPeakWithMaxDeviation(signal,
586  nExponentialsForFit,
587  startT,
588  endT,
589  fSameShape,
590  paramVec,
591  peakVals,
592  PeakDev);
593 
594  //Add peak and re-fit
595  for (auto& PeakDevCand : PeakDev) {
596  chi2PerNDF2 = 0.;
597  NDF2 = 0.;
598  ParameterVec paramVecRefit;
599  peakValsTemp = peakVals;
600 
601  AddPeak(PeakDevCand, peakValsTemp);
602  FitExponentials(signal,
603  peakValsTemp,
604  startT,
605  endT,
606  paramVecRefit,
607  chi2PerNDF2,
608  NDF2,
609  fSameShape);
610 
611  if (chi2PerNDF2 < chi2PerNDF) {
612  paramVec = paramVecRefit;
613  peakVals = peakValsTemp;
614  nExponentialsForFit = peakVals.size();
615  chi2PerNDF = chi2PerNDF2;
616  NDF = NDF2;
617  nExponentialsAfterRefit++;
618  RefitSuccess = true;
619  break;
620  }
621  }
622 
623  //Split peak and re-fit
624  if (RefitSuccess == false) {
625  for (auto& PeakDevCand : PeakDev) {
626  chi2PerNDF2 = 0.;
627  NDF2 = 0.;
628  ParameterVec paramVecRefit;
629  peakValsTemp = peakVals;
630 
631  SplitPeak(PeakDevCand, peakValsTemp);
632  FitExponentials(signal,
633  peakValsTemp,
634  startT,
635  endT,
636  paramVecRefit,
637  chi2PerNDF2,
638  NDF2,
639  fSameShape);
640 
641  if (chi2PerNDF2 < chi2PerNDF) {
642  paramVec = paramVecRefit;
643  peakVals = peakValsTemp;
644  nExponentialsForFit = peakVals.size();
645  chi2PerNDF = chi2PerNDF2;
646  NDF = NDF2;
647  nExponentialsAfterRefit++;
648  RefitSuccess = true;
649  break;
650  }
651  }
652  }
653 
654  if (RefitSuccess == false) { break; }
655  }
656 
657  if (fLogLevel >= 5) {
658  std::cout << std::endl;
659  std::cout << "--- Refit ---" << std::endl;
660  if (chi2PerNDF == oldChi2PerNDF)
661  std::cout << "chi2/ndf didn't improve. Keep first fit." << std::endl;
662  else {
663  std::cout << "- Added peaks to group #" << j << ". This group now has "
664  << nExponentialsForFit << " peaks:" << std::endl;
665  std::cout << "- Group #" << j << " (" << peakVals.size()
666  << " peaks): GroupStartTick: " << range.offset + startT
667  << " GroupEndTick: " << range.offset + endT << std::endl;
668 
669  int CountPeakInGroup = 0;
670  for (auto const& peakValsTmp : peakVals) {
671  std::cout << "Peak #" << CountPeakInGroup << " in group #" << j
672  << ": PeakInGroupStartTick: "
673  << range.offset + std::get<2>(peakValsTmp)
674  << " PeakInGroupMaxTick: "
675  << range.offset + std::get<0>(peakValsTmp)
676  << " PeakInGroupEndTick: "
677  << range.offset + std::get<3>(peakValsTmp) << std::endl;
678  CountPeakInGroup++;
679  }
680 
681  std::cout << "chi2/ndf = " << chi2PerNDF << std::endl;
682 
683  if (fSameShape) {
684  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
685  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
686 
687  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
688  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2 * (i + 1)].first
689  << std::endl;
690  std::cout << "Peak #" << i << ": t0 [ticks] = "
691  << range.offset + paramVec[2 * (i + 1) + 1].first << std::endl;
692  }
693  }
694  else {
695  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
696  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[4 * i + 2].first
697  << std::endl;
698  std::cout << "Peak #" << i
699  << ": t0 [ticks] = " << range.offset + paramVec[4 * i + 3].first
700  << std::endl;
701  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4 * i].first
702  << std::endl;
703  std::cout << "Peak #" << i << ": tau2 [mus] = " << paramVec[4 * i + 1].first
704  << std::endl;
705  }
706  }
707  }
708  }
709  }
710 
711  // #######################################################
712  // ### Loop through returned peaks and make recob hits ###
713  // #######################################################
714 
715  int numHits(0);
716  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
717  //Extract fit parameters for this hit
718  double peakTau1;
719  double peakTau2;
720  double peakAmp;
721  double peakMean;
722 
723  if (fSameShape) {
724  peakTau1 = paramVec[0].first;
725  peakTau2 = paramVec[1].first;
726  peakAmp = paramVec[2 * (i + 1)].first;
727  peakMean = paramVec[2 * (i + 1) + 1].first;
728  }
729  else {
730  peakTau1 = paramVec[4 * i].first;
731  peakTau2 = paramVec[4 * i + 1].first;
732  peakAmp = paramVec[4 * i + 2].first;
733  peakMean = paramVec[4 * i + 3].first;
734  }
735 
736  //Highest ADC count in peak = peakAmpTrue
737  double peakAmpTrue = signal[std::get<0>(peakVals.at(i))];
738  double peakAmpErr = 1.;
739 
740  //Determine peak position of fitted function (= peakMeanTrue)
741  TF1 Exponentials("Exponentials",
742  "( [0] * exp(0.4*(x-[1])/[2]) / ( 1 + exp(0.4*(x-[1])/[3]) ) )",
743  startT,
744  endT);
745  Exponentials.SetParameter(0, peakAmp);
746  Exponentials.SetParameter(1, peakMean);
747  Exponentials.SetParameter(2, peakTau1);
748  Exponentials.SetParameter(3, peakTau2);
749  double peakMeanTrue = Exponentials.GetMaximumX(startT, endT);
750  Exponentials.Delete();
751 
752  //Calculate width (=FWHM)
753  double peakWidth =
754  WidthFunc(peakMean, peakAmp, peakTau1, peakTau2, startT, endT, peakMeanTrue);
755  peakWidth /=
756  fWidthNormalization; //from FWHM to "standard deviation": standard deviation = FWHM/(2*sqrt(2*ln(2)))
757 
758  // Extract fit parameter errors
759  double peakMeanErr;
760 
761  if (fSameShape) { peakMeanErr = paramVec[2 * (i + 1) + 1].second; }
762  else {
763  peakMeanErr = paramVec[4 * i + 3].second;
764  }
765  double peakWidthErr = 0.1 * peakWidth;
766 
767  // ### Charge ###
768  double charge =
769  ChargeFunc(peakMean, peakAmp, peakTau1, peakTau2, fChargeNorm, peakMeanTrue);
770  double chargeErr =
771  std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
772 
773  // ### limits for getting sum of ADC counts
774  int startTthisHit = std::get<2>(peakVals.at(i));
775  int endTthisHit = std::get<3>(peakVals.at(i));
776  std::vector<float>::const_iterator sumStartItr = signal.begin() + startTthisHit;
777  std::vector<float>::const_iterator sumEndItr = signal.begin() + endTthisHit;
778 
779  // ### Sum of ADC counts
780  double sumADC = std::accumulate(sumStartItr, sumEndItr + 1, 0.);
781 
782  //Check if fit returns reasonable values and ich chi2 is below threshold
783  if (peakWidth <= 0 || charge <= 0. || charge != charge ||
784  (nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax) ||
785  (nExponentialsForFit >= 2 &&
786  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)) {
787  if (fLogLevel >= 1) {
788  std::cout << std::endl;
789  std::cout << "WARNING: For peak #" << i << " in this group:" << std::endl;
790  if (peakWidth <= 0 || charge <= 0. || charge != charge)
791  std::cout << "Fit function returned width < 0 or charge < 0 or charge = nan."
792  << std::endl;
793  if ((nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax) ||
794  (nExponentialsForFit >= 2 &&
795  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)) {
796  std::cout << std::endl;
797  std::cout << "WARNING: For fit of this group (" << NumberOfPeaksBeforeFit
798  << " peaks before refit, " << nExponentialsForFit
799  << " peaks after refit): " << std::endl;
800  if (nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax)
801  std::cout << "chi2/ndf of this fit (" << chi2PerNDF
802  << ") is higher than threshold (" << fChi2NDFMax << ")."
803  << std::endl;
804  if (nExponentialsForFit >= 2 &&
805  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)
806  std::cout << "chi2/ndf of this fit (" << chi2PerNDF
807  << ") is higher than threshold ("
808  << fChi2NDFMaxFactorMultiHits * fChi2NDFMax << ")." << std::endl;
809  }
810  std::cout << "---> DO NOT create hit object from fit parameters but use peak "
811  "values instead."
812  << std::endl;
813  std::cout << "---> Set fit parameter so that a sharp peak with a width of 1 tick "
814  "is shown in the event display. This indicates that the fit failed."
815  << std::endl;
816  }
817  peakWidth =
818  (((double)endTthisHit - (double)startTthisHit) / 4.) /
819  fWidthNormalization; //~4 is the factor between FWHM and full width of the hit (last bin - first bin). no drift: 4.4, 6m drift: 3.7
820  peakMeanErr = peakWidth / 2;
821  charge = sumADC;
822  peakMeanTrue = std::get<0>(peakVals.at(i));
823 
824  //set the fit values to make it visible in the event display that this fit failed
825  peakMean = peakMeanTrue;
826  peakTau1 = 0.008;
827  peakTau2 = 0.0065;
828  peakAmp = 20.;
829  }
830 
831  // Create the hit
832  recob::HitCreator hitcreator(
833  *wire, // wire reference
834  wid, // wire ID
835  startTthisHit + roiFirstBinTick, // start_tick TODO check
836  endTthisHit + roiFirstBinTick, // end_tick TODO check
837  peakWidth, // rms
838  peakMeanTrue + roiFirstBinTick, // peak_time
839  peakMeanErr, // sigma_peak_time
840  peakAmpTrue, // peak_amplitude
841  peakAmpErr, // sigma_peak_amplitude
842  charge, // hit_integral
843  chargeErr, // hit_sigma_integral
844  sumADC, // summedADC FIXME
845  nExponentialsForFit, // multiplicity
846  numHits, // local_index TODO check that the order is correct
847  chi2PerNDF, // goodness_of_fit
848  NDF // dof
849  );
850 
851  if (fLogLevel >= 6) {
852  std::cout << std::endl;
853  std::cout << "- Created hit object for peak #" << i
854  << " in this group with the following parameters (obtained from fit):"
855  << std::endl;
856  std::cout << "HitStartTick: " << startTthisHit + roiFirstBinTick << std::endl;
857  std::cout << "HitEndTick: " << endTthisHit + roiFirstBinTick << std::endl;
858  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
859  std::cout << "HitMeanTick: " << peakMeanTrue + roiFirstBinTick << " +- "
860  << peakMeanErr << std::endl;
861  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr
862  << std::endl;
863  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr
864  << std::endl;
865  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
866  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
867  std::cout << "HitIndex in group: " << numHits << std::endl;
868  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
869  std::cout << "HitNDF: " << NDF << std::endl;
870  }
871 
872  const recob::Hit hit(hitcreator.move());
873 
874  hcol.emplace_back(std::move(hit), wire, rawdigits);
875  // add fit parameters associated to the hit just pushed to the collection
876  std::array<float, 4> fitParams;
877  fitParams[0] = peakMean + roiFirstBinTick;
878  fitParams[1] = peakTau1;
879  fitParams[2] = peakTau2;
880  fitParams[3] = peakAmp;
881  fHitParamWriter.addVector(hitID, fitParams);
882  numHits++;
883  } // <---End loop over Exponentials
884  } // <---End if(NumberOfPeaksBeforeFit <= fMaxMultiHit && width <= fMaxGroupLength), then fit
885 
886  // #######################################################
887  // ### If too large then force alternate solution ###
888  // ### - Make n hits from pulse train where n will ###
889  // ### depend on the fhicl parameter fLongPulseWidth ###
890  // ### Also do this if chi^2 is too large ###
891  // #######################################################
892  if (NumberOfPeaksBeforeFit > fMaxMultiHit || (width > fMaxGroupLength) ||
893  NFluctuations > fMaxFluctuations) {
894 
895  int nHitsInThisGroup = (endT - startT + 1) / fLongPulseWidth;
896 
897  if (nHitsInThisGroup > fLongMaxHits) {
898  nHitsInThisGroup = fLongMaxHits;
899  fLongPulseWidth = (endT - startT + 1) / nHitsInThisGroup;
900  }
901 
902  if (nHitsInThisGroup * fLongPulseWidth < (endT - startT + 1)) nHitsInThisGroup++;
903 
904  int firstTick = startT;
905  int lastTick = std::min(endT, firstTick + fLongPulseWidth - 1);
906 
907  if (fLogLevel >= 1) {
908  if (NumberOfPeaksBeforeFit > fMaxMultiHit) {
909  std::cout << std::endl;
910  std::cout << "WARNING: Number of peaks in this group (" << NumberOfPeaksBeforeFit
911  << ") is higher than threshold (" << fMaxMultiHit << ")." << std::endl;
912  std::cout
913  << "---> DO NOT fit. Split group of peaks into hits with equal length instead."
914  << std::endl;
915  }
916  if (width > fMaxGroupLength) {
917  std::cout << std::endl;
918  std::cout << "WARNING: group of peak is longer (" << width
919  << " ticks) than threshold (" << fMaxGroupLength << " ticks)."
920  << std::endl;
921  std::cout
922  << "---> DO NOT fit. Split group of peaks into hits with equal length instead."
923  << std::endl;
924  }
925  if (NFluctuations > fMaxFluctuations) {
926  std::cout << std::endl;
927  std::cout << "WARNING: fluctuations (" << NFluctuations
928  << ") higher than threshold (" << fMaxFluctuations << ")." << std::endl;
929  std::cout
930  << "---> DO NOT fit. Split group of peaks into hits with equal length instead."
931  << std::endl;
932  }
933  std::cout << "---> Group goes from tick " << roiFirstBinTick + startT << " to "
934  << roiFirstBinTick + endT << ". Split group into ("
935  << roiFirstBinTick + endT << " - " << roiFirstBinTick + startT << ")/"
936  << fLongPulseWidth << " = " << (endT - startT) << "/" << fLongPulseWidth
937  << " = " << nHitsInThisGroup << " peaks (" << fLongPulseWidth
938  << " = LongPulseWidth), or maximum LongMaxHits = " << fLongMaxHits
939  << " peaks." << std::endl;
940  }
941 
942  for (int hitIdx = 0; hitIdx < nHitsInThisGroup; hitIdx++) {
943  // This hit parameters
944  double peakWidth =
945  ((lastTick - firstTick) / 4.) /
946  fWidthNormalization; //~4 is the factor between FWHM and full width of the hit (last bin - first bin). no drift: 4.4, 6m drift: 3.7
947  double peakMeanTrue = (firstTick + lastTick) / 2.;
948  if (NumberOfPeaksBeforeFit == 1 && nHitsInThisGroup == 1)
949  peakMeanTrue = std::get<0>(peakVals.at(
950  0)); //if only one peak was found, we want the mean of this peak to be the tick with the max. ADC count
951  double peakMeanErr = (lastTick - firstTick) / 2.;
952  double sumADC =
953  std::accumulate(signal.begin() + firstTick, signal.begin() + lastTick + 1, 0.);
954  double charge = sumADC;
955  double chargeErr = 0.1 * sumADC;
956  double peakAmpTrue = 0;
957 
958  for (int tick = firstTick; tick <= lastTick; tick++) {
959  if (signal[tick] > peakAmpTrue) peakAmpTrue = signal[tick];
960  }
961 
962  double peakAmpErr = 1.;
963  nExponentialsForFit = nHitsInThisGroup;
964  NDF = -1;
965  chi2PerNDF = -1.;
966  //set the fit values to make it visible in the event display that this fit failed
967  double peakMean = peakMeanTrue - 2;
968  double peakTau1 = 0.008;
969  double peakTau2 = 0.0065;
970  double peakAmp = 20.;
971 
972  recob::HitCreator hitcreator(
973  *wire, // wire reference
974  wid, // wire ID
975  firstTick + roiFirstBinTick, // start_tick TODO check
976  lastTick + roiFirstBinTick, // end_tick TODO check
977  peakWidth, // rms
978  peakMeanTrue + roiFirstBinTick, // peak_time
979  peakMeanErr, // sigma_peak_time
980  peakAmpTrue, // peak_amplitude
981  peakAmpErr, // sigma_peak_amplitude
982  charge, // hit_integral
983  chargeErr, // hit_sigma_integral
984  sumADC, // summedADC FIXME
985  nExponentialsForFit, // multiplicity
986  hitIdx, // local_index TODO check that the order is correct
987  chi2PerNDF, // goodness_of_fit
988  NDF // dof
989  );
990 
991  if (fLogLevel >= 6) {
992  std::cout << std::endl;
993  std::cout
994  << "- Created hit object for peak #" << hitIdx
995  << " in this group with the following parameters (obtained from waveform):"
996  << std::endl;
997  std::cout << "HitStartTick: " << firstTick + roiFirstBinTick << std::endl;
998  std::cout << "HitEndTick: " << lastTick + roiFirstBinTick << std::endl;
999  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
1000  std::cout << "HitMeanTick: " << peakMeanTrue + roiFirstBinTick << " +- "
1001  << peakMeanErr << std::endl;
1002  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr
1003  << std::endl;
1004  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr
1005  << std::endl;
1006  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
1007  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
1008  std::cout << "HitIndex in group: " << hitIdx << std::endl;
1009  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
1010  std::cout << "HitNDF: " << NDF << std::endl;
1011  }
1012  const recob::Hit hit(hitcreator.move());
1013  hcol.emplace_back(std::move(hit), wire, rawdigits);
1014 
1015  std::array<float, 4> fitParams;
1016  fitParams[0] = peakMean + roiFirstBinTick;
1017  fitParams[1] = peakTau1;
1018  fitParams[2] = peakTau2;
1019  fitParams[3] = peakAmp;
1020  fHitParamWriter.addVector(hitID, fitParams);
1021 
1022  // set for next loop
1023  firstTick = lastTick + 1;
1024  lastTick = std::min(firstTick + fLongPulseWidth - 1, endT);
1025 
1026  } //<---Hits in this group
1027  } //<---End if #peaks > MaxMultiHit
1028  fChi2->Fill(chi2PerNDF);
1029  } //<---End loop over merged candidate hits
1030  } //<---End looping over ROI's
1031  } //<---End looping over all the wires
1032 
1033  //==================================================================================================
1034  // End of the event
1035 
1036  // move the hit collection and the associations into the event
1037  hcol.put_into(evt);
1038 
1039  // and put hit fit parameters together with metadata into the event
1040  fHitParamWriter.saveOutputs(evt);
1041 
1042  } // End of produce()
1043 
1044  // --------------------------------------------------------------------------------------------
1045  // Initial finding of candidate peaks
1046  // --------------------------------------------------------------------------------------------
1047  void hit::DPRawHitFinder::findCandidatePeaks(std::vector<float>::const_iterator startItr,
1048  std::vector<float>::const_iterator stopItr,
1049  std::vector<std::tuple<int, int, int>>& timeValsVec,
1050  float& PeakMin,
1051  int firstTick) const
1052  {
1053  // Need a minimum number of ticks to do any work here
1054  if (std::distance(startItr, stopItr) > 4) {
1055  // Find the highest peak in the range given
1056  auto maxItr = std::max_element(startItr, stopItr);
1057 
1058  float maxValue = *maxItr;
1059  int maxTime = std::distance(startItr, maxItr);
1060 
1061  if (maxValue >= PeakMin) {
1062  // backwards to find first bin for this candidate hit
1063  auto firstItr = std::distance(startItr, maxItr) > 2 ? maxItr - 1 : startItr;
1064  bool PeakStartIsHere = true;
1065 
1066  while (firstItr != startItr) {
1067  //Check for inflection point & ADC count <= 0
1068  PeakStartIsHere = true;
1069  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1070  if (*firstItr >= *(firstItr - i)) {
1071  PeakStartIsHere = false;
1072  break;
1073  }
1074  }
1075  if (*firstItr <= 0 || PeakStartIsHere) break;
1076  firstItr--;
1077  }
1078 
1079  int firstTime = std::distance(startItr, firstItr);
1080 
1081  // Recursive call to find all candidate hits earlier than this peak
1082  findCandidatePeaks(startItr, firstItr - 1, timeValsVec, PeakMin, firstTick);
1083 
1084  // forwards to find last bin for this candidate hit
1085  auto lastItr = std::distance(maxItr, stopItr) > 2 ? maxItr + 1 : stopItr - 1;
1086  bool PeakEndIsHere = true;
1087 
1088  while (lastItr != stopItr) {
1089  //Check for inflection point & ADC count <= 0
1090  PeakEndIsHere = true;
1091  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1092  if (*lastItr >= *(lastItr + i)) {
1093  PeakEndIsHere = false;
1094  break;
1095  }
1096  }
1097  if (*lastItr <= 0 || PeakEndIsHere) break;
1098  lastItr++;
1099  }
1100 
1101  int lastTime = std::distance(startItr, lastItr);
1102 
1103  // Now save this candidate's start and max time info
1104  timeValsVec.push_back(
1105  std::make_tuple(firstTick + firstTime, firstTick + maxTime, firstTick + lastTime));
1106 
1107  // Recursive call to find all candidate hits later than this peak
1108  findCandidatePeaks(lastItr + 1,
1109  stopItr,
1110  timeValsVec,
1111  PeakMin,
1112  firstTick + std::distance(startItr, lastItr + 1));
1113  }
1114  }
1115 
1116  return;
1117  }
1118 
1119  // --------------------------------------------------------------------------------------------
1120  // Merging of nearby candidate peaks
1121  // --------------------------------------------------------------------------------------------
1122 
1123  void hit::DPRawHitFinder::mergeCandidatePeaks(const std::vector<float> signalVec,
1124  TimeValsVec timeValsVec,
1125  MergedTimeWidVec& mergedVec)
1126  {
1127  // ################################################################
1128  // ### Lets loop over the candidate pulses we found in this ROI ###
1129  // ################################################################
1130  auto timeValsVecItr = timeValsVec.begin();
1131  unsigned int PeaksInThisMergedPeak = 0;
1132  //int EndTickOfPreviousMergedPeak=0;
1133  while (timeValsVecItr != timeValsVec.end()) {
1134  PeakTimeWidVec peakVals;
1135 
1136  // Setting the start, peak, and end time of the pulse
1137  auto& timeVal = *timeValsVecItr++;
1138  int startT = std::get<0>(timeVal);
1139  int maxT = std::get<1>(timeVal);
1140  int endT = std::get<2>(timeVal);
1141  int widT = endT - startT;
1142  int FinalStartT = startT;
1143  int FinalEndT = endT;
1144 
1145  int NFluctuations = EstimateFluctuations(signalVec, startT, maxT, endT);
1146 
1147  peakVals.emplace_back(maxT, widT, startT, endT);
1148 
1149  // See if we want to merge pulses together
1150  // First check if we have more than one pulse on the wire
1151  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1152 
1153  // Loop until no more merged pulses (or candidates in this ROI)
1154  while (checkNextHit) {
1155  // group hits if start time of the next pulse is < end time + fGroupMaxDistance of current pulse and if intermediate signal between two pulses doesn't go below fMinBinToGroup and if the hits are not all above the merge max adc limit
1156  int NextStartT = std::get<0>(*timeValsVecItr);
1157 
1158  double MinADC = signalVec[endT];
1159  for (int i = endT; i <= NextStartT; i++) {
1160  if (signalVec[i] < MinADC) { MinADC = signalVec[i]; }
1161  }
1162 
1163  // Group peaks (grouped peaks are fitted together and can be merged)
1164  if (MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance) {
1165  int CurrentStartT = startT;
1166  int CurrentMaxT = maxT;
1167  int CurrentEndT = endT;
1168 
1169  timeVal = *timeValsVecItr++;
1170  int NextMaxT = std::get<1>(timeVal);
1171  int NextEndT = std::get<2>(timeVal);
1172  int NextWidT = NextEndT - NextStartT;
1173  FinalEndT = NextEndT;
1174 
1175  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1176 
1177  int CurrentSumADC = 0;
1178  for (int i = CurrentStartT; i <= CurrentEndT; i++) {
1179  CurrentSumADC += signalVec[i];
1180  }
1181 
1182  int NextSumADC = 0;
1183  for (int i = NextStartT; i <= NextEndT; i++) {
1184  NextSumADC += signalVec[i];
1185  }
1186 
1187  //Merge peaks within a group
1188  if (fDoMergePeaks) {
1189  if (signalVec[NextMaxT] <= signalVec[CurrentMaxT] &&
1190  ((signalVec[NextMaxT] < fMergeMaxADCThreshold * signalVec[CurrentMaxT] &&
1191  NextSumADC < fMergeADCSumThreshold * CurrentSumADC) ||
1192  (signalVec[NextMaxT] - signalVec[NextStartT]) <
1193  fMinRelativePeakHeightRight * signalVec[CurrentMaxT]) &&
1194  (signalVec[NextMaxT] <= fMergeMaxADCLimit)) {
1195  maxT = CurrentMaxT;
1196  startT = CurrentStartT;
1197  endT = NextEndT;
1198  widT = endT - startT;
1199  peakVals.pop_back();
1200  peakVals.emplace_back(maxT, widT, startT, endT);
1201  }
1202  else if (signalVec[NextMaxT] > signalVec[CurrentMaxT] &&
1203  ((signalVec[CurrentMaxT] < fMergeMaxADCThreshold * signalVec[NextMaxT] &&
1204  CurrentSumADC < fMergeADCSumThreshold * NextSumADC) ||
1205  (signalVec[CurrentMaxT] - signalVec[CurrentEndT]) <
1206  fMinRelativePeakHeightLeft * signalVec[NextMaxT]) &&
1207  (signalVec[CurrentMaxT] <= fMergeMaxADCLimit)) {
1208  maxT = NextMaxT;
1209  startT = CurrentStartT;
1210  endT = NextEndT;
1211  widT = endT - startT;
1212  peakVals.pop_back();
1213  peakVals.emplace_back(maxT, widT, startT, endT);
1214  }
1215  else {
1216  maxT = NextMaxT;
1217  startT = NextStartT;
1218  endT = NextEndT;
1219  widT = NextWidT;
1220  peakVals.emplace_back(maxT, widT, startT, endT);
1221  PeaksInThisMergedPeak++;
1222  }
1223  }
1224  else {
1225  maxT = NextMaxT;
1226  startT = NextStartT;
1227  endT = NextEndT;
1228  widT = NextWidT;
1229  peakVals.emplace_back(maxT, widT, startT, endT);
1230  PeaksInThisMergedPeak++;
1231  }
1232  checkNextHit = timeValsVecItr != timeValsVec.end();
1233  } //<---Checking adjacent pulses
1234  else {
1235  checkNextHit = false;
1236  PeaksInThisMergedPeak = 0;
1237  }
1238 
1239  } //<---End checking if there is more than one pulse on the wire
1240  // Add these to our merged vector
1241  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1242  }
1243  return;
1244  }
1245 
1246  // ----------------------------------------------------------------------------------------------
1247  // Estimate fluctuations for a group of peaks to identify hits from particles in drift direction
1248  // ----------------------------------------------------------------------------------------------
1249  int hit::DPRawHitFinder::EstimateFluctuations(const std::vector<float> fsignalVec,
1250  int peakStart,
1251  int peakMean,
1252  int peakEnd)
1253  {
1254  int NFluctuations = 0;
1255 
1256  for (int j = peakMean - 1; j >= peakStart; j--) {
1257  if (fsignalVec[j] < 5)
1258  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1259 
1260  if (fsignalVec[j] > fsignalVec[j + 1]) { NFluctuations++; }
1261  }
1262 
1263  for (int j = peakMean + 1; j <= peakEnd; j++) {
1264  if (fsignalVec[j] < 5)
1265  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1266 
1267  if (fsignalVec[j] > fsignalVec[j - 1]) { NFluctuations++; }
1268  }
1269 
1270  return NFluctuations;
1271  }
1272 
1273  // --------------------------------------------------------------------------------------------
1274  // Fit Exponentials
1275  // --------------------------------------------------------------------------------------------
1276  void hit::DPRawHitFinder::FitExponentials(const std::vector<float> fSignalVector,
1277  const PeakTimeWidVec fPeakVals,
1278  int fStartTime,
1279  int fEndTime,
1280  ParameterVec& fparamVec,
1281  double& fchi2PerNDF,
1282  int& fNDF,
1283  bool fSameShape)
1284  {
1285  int size = fEndTime - fStartTime + 1;
1286  int NPeaks = fPeakVals.size();
1287 
1288  // #############################################
1289  // ### If size < 0 then set the size to zero ###
1290  // #############################################
1291  if (fEndTime - fStartTime < 0) { size = 0; }
1292 
1293  // --- TH1D HitSignal ---
1294  TH1F hitSignal("hitSignal", "", std::max(size, 1), fStartTime, fEndTime + 1);
1295  hitSignal.Sumw2();
1296 
1297  // #############################
1298  // ### Filling the histogram ###
1299  // #############################
1300  for (int i = fStartTime; i < fEndTime + 1; i++) {
1301  hitSignal.Fill(i, fSignalVector[i]);
1302  hitSignal.SetBinError(i, 0.288675); //1/sqrt(12)
1303  }
1304 
1305  // ################################################
1306  // ### Building TFormula for basic Exponentials ###
1307  // ################################################
1308  std::string eqn = CreateFitFunction(NPeaks, fSameShape);
1309 
1310  // -------------------------------------
1311  // --- TF1 function for Exponentials ---
1312  // -------------------------------------
1313 
1314  TF1 Exponentials("Exponentials", eqn.c_str(), fStartTime, fEndTime + 1);
1315 
1316  if (fLogLevel >= 4) {
1317  std::cout << std::endl;
1318  std::cout << "--- Preparing fit ---" << std::endl;
1319  std::cout << "--- Lower limits, seed, upper limit:" << std::endl;
1320  }
1321 
1322  if (fSameShape) {
1323  Exponentials.SetParameter(0, 0.5);
1324  Exponentials.SetParameter(1, 0.5);
1325  Exponentials.SetParLimits(0, fMinTau, fMaxTau);
1326  Exponentials.SetParLimits(1, fMinTau, fMaxTau);
1327  double amplitude = 0;
1328  double peakMean = 0;
1329 
1330  double peakMeanShift = 2;
1331  double peakMeanSeed = 0;
1332  double peakMeanRangeLow = 0;
1333  double peakMeanRangeHi = 0;
1334  double peakStart = 0;
1335  double peakEnd = 0;
1336 
1337  for (int i = 0; i < NPeaks; i++) {
1338  peakMean = std::get<0>(fPeakVals.at(i));
1339  peakStart = std::get<2>(fPeakVals.at(i));
1340  peakEnd = std::get<3>(fPeakVals.at(i));
1341  peakMeanSeed = peakMean - peakMeanShift;
1342  peakMeanRangeLow = std::max(peakStart - peakMeanShift, peakMeanSeed - fFitPeakMeanRange);
1343  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed + fFitPeakMeanRange);
1344  amplitude = fSignalVector[peakMean];
1345 
1346  Exponentials.SetParameter(2 * (i + 1), 1.65 * amplitude);
1347  Exponentials.SetParLimits(2 * (i + 1), 0.3 * 1.65 * amplitude, 2 * 1.65 * amplitude);
1348  Exponentials.SetParameter(2 * (i + 1) + 1, peakMeanSeed);
1349 
1350  if (NPeaks == 1) {
1351  Exponentials.SetParLimits(2 * (i + 1) + 1, peakMeanRangeLow, peakMeanRangeHi);
1352  }
1353  else if (NPeaks >= 2 && i == 0) {
1354  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1355  Exponentials.SetParLimits(
1356  2 * (i + 1) + 1,
1357  peakMeanRangeLow,
1358  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1359  }
1360  else if (NPeaks >= 2 && i == NPeaks - 1) {
1361  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1362  Exponentials.SetParLimits(
1363  2 * (i + 1) + 1,
1364  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1365  peakMeanRangeHi);
1366  }
1367  else {
1368  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1369  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1370  Exponentials.SetParLimits(
1371  2 * (i + 1) + 1,
1372  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1373  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1374  }
1375 
1376  if (fLogLevel >= 4) {
1377  double t0low, t0high;
1378  Exponentials.GetParLimits(2 * (i + 1) + 1, t0low, t0high);
1379  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3 * 1.65 * amplitude << " , "
1380  << 1.65 * amplitude << " , " << 2 * 1.65 * amplitude << std::endl;
1381  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed
1382  << " , " << t0high << std::endl;
1383  }
1384  }
1385  }
1386  else {
1387  double amplitude = 0;
1388  double peakMean = 0;
1389 
1390  double peakMeanShift = 2;
1391  double peakMeanSeed = 0;
1392  double peakMeanRangeLow = 0;
1393  double peakMeanRangeHi = 0;
1394  double peakStart = 0;
1395  double peakEnd = 0;
1396 
1397  for (int i = 0; i < NPeaks; i++) {
1398  Exponentials.SetParameter(4 * i, 0.5);
1399  Exponentials.SetParameter(4 * i + 1, 0.5);
1400  Exponentials.SetParLimits(4 * i, fMinTau, fMaxTau);
1401  Exponentials.SetParLimits(4 * i + 1, fMinTau, fMaxTau);
1402 
1403  peakMean = std::get<0>(fPeakVals.at(i));
1404  peakStart = std::get<2>(fPeakVals.at(i));
1405  peakEnd = std::get<3>(fPeakVals.at(i));
1406  peakMeanSeed = peakMean - peakMeanShift;
1407  peakMeanRangeLow = std::max(peakStart - peakMeanShift, peakMeanSeed - fFitPeakMeanRange);
1408  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed + fFitPeakMeanRange);
1409  amplitude = fSignalVector[peakMean];
1410 
1411  Exponentials.SetParameter(4 * i + 2, 1.65 * amplitude);
1412  Exponentials.SetParLimits(4 * i + 2, 0.3 * 1.65 * amplitude, 2 * 1.65 * amplitude);
1413  Exponentials.SetParameter(4 * i + 3, peakMeanSeed);
1414 
1415  if (NPeaks == 1) {
1416  Exponentials.SetParLimits(4 * i + 3, peakMeanRangeLow, peakMeanRangeHi);
1417  }
1418  else if (NPeaks >= 2 && i == 0) {
1419  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1420  Exponentials.SetParLimits(
1421  4 * i + 3,
1422  peakMeanRangeLow,
1423  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1424  }
1425  else if (NPeaks >= 2 && i == NPeaks - 1) {
1426  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1427  Exponentials.SetParLimits(
1428  4 * i + 3,
1429  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1430  peakMeanRangeHi);
1431  }
1432  else {
1433  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1434  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1435  Exponentials.SetParLimits(
1436  4 * i + 3,
1437  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1438  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1439  }
1440 
1441  if (fLogLevel >= 4) {
1442  double t0low, t0high;
1443  Exponentials.GetParLimits(4 * i + 3, t0low, t0high);
1444  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3 * 1.65 * amplitude << " , "
1445  << 1.65 * amplitude << " , " << 2 * 1.65 * amplitude << std::endl;
1446  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed
1447  << " , " << t0high << std::endl;
1448  }
1449  }
1450  }
1451 
1452  // ###########################################
1453  // ### PERFORMING THE TOTAL FIT OF THE HIT ###
1454  // ###########################################
1455  try {
1456  hitSignal.Fit(&Exponentials, "QNRWM", "", fStartTime, fEndTime + 1);
1457  }
1458  catch (...) {
1459  mf::LogWarning("DPRawHitFinder") << "Fitter failed finding a hit";
1460  }
1461 
1462  // ##################################################
1463  // ### Getting the fitted parameters from the fit ###
1464  // ##################################################
1465  fchi2PerNDF = (Exponentials.GetChisquare() / Exponentials.GetNDF());
1466  fNDF = Exponentials.GetNDF();
1467 
1468  if (fSameShape) {
1469  fparamVec.emplace_back(Exponentials.GetParameter(0), Exponentials.GetParError(0));
1470  fparamVec.emplace_back(Exponentials.GetParameter(1), Exponentials.GetParError(1));
1471 
1472  for (int i = 0; i < NPeaks; i++) {
1473  fparamVec.emplace_back(Exponentials.GetParameter(2 * (i + 1)),
1474  Exponentials.GetParError(2 * (i + 1)));
1475  fparamVec.emplace_back(Exponentials.GetParameter(2 * (i + 1) + 1),
1476  Exponentials.GetParError(2 * (i + 1) + 1));
1477  }
1478  }
1479  else {
1480  for (int i = 0; i < NPeaks; i++) {
1481  fparamVec.emplace_back(Exponentials.GetParameter(4 * i), Exponentials.GetParError(4 * i));
1482  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 1),
1483  Exponentials.GetParError(4 * i + 1));
1484  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 2),
1485  Exponentials.GetParError(4 * i + 2));
1486  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 3),
1487  Exponentials.GetParError(4 * i + 3));
1488  }
1489  }
1490  Exponentials.Delete();
1491  hitSignal.Delete();
1492  } //<----End FitExponentials
1493 
1494  //---------------------------------------------------------------------------------------------
1495  void hit::DPRawHitFinder::FindPeakWithMaxDeviation(const std::vector<float> fSignalVector,
1496  int fNPeaks,
1497  int fStartTime,
1498  int fEndTime,
1499  bool fSameShape,
1500  ParameterVec fparamVec,
1501  PeakTimeWidVec fpeakVals,
1502  PeakDevVec& fPeakDev)
1503  {
1504  std::string eqn =
1505  CreateFitFunction(fNPeaks, fSameShape); // string for equation of Exponentials fit
1506 
1507  TF1 Exponentials("Exponentials", eqn.c_str(), fStartTime, fEndTime + 1);
1508 
1509  for (size_t i = 0; i < fparamVec.size(); i++) {
1510  Exponentials.SetParameter(i, fparamVec[i].first);
1511  }
1512 
1513  // ##########################################################################
1514  // ### Finding the peak with the max chi2 fit and signal ###
1515  // ##########################################################################
1516  double Chi2PerNDFPeak;
1517  double MaxPosDeviation;
1518  double MaxNegDeviation;
1519  int BinMaxPosDeviation;
1520  int BinMaxNegDeviation;
1521  for (int i = 0; i < fNPeaks; i++) {
1522  Chi2PerNDFPeak = 0.;
1523  MaxPosDeviation = 0.;
1524  MaxNegDeviation = 0.;
1525  BinMaxPosDeviation = 0;
1526  BinMaxNegDeviation = 0;
1527 
1528  for (int j = std::get<2>(fpeakVals.at(i)); j < std::get<3>(fpeakVals.at(i)) + 1; j++) {
1529  if ((Exponentials(j + 0.5) - fSignalVector[j]) > MaxPosDeviation &&
1530  j != std::get<0>(fpeakVals.at(i))) {
1531  MaxPosDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1532  BinMaxPosDeviation = j;
1533  }
1534  if ((Exponentials(j + 0.5) - fSignalVector[j]) < MaxNegDeviation &&
1535  j != std::get<0>(fpeakVals.at(i))) {
1536  MaxNegDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1537  BinMaxNegDeviation = j;
1538  }
1539  Chi2PerNDFPeak +=
1540  pow((Exponentials(j + 0.5) - fSignalVector[j]) / sqrt(fSignalVector[j]), 2);
1541  }
1542 
1543  if (BinMaxNegDeviation != 0) {
1544  Chi2PerNDFPeak /=
1545  static_cast<double>((std::get<3>(fpeakVals.at(i)) - std::get<2>(fpeakVals.at(i))));
1546  fPeakDev.emplace_back(Chi2PerNDFPeak, i, BinMaxNegDeviation, BinMaxPosDeviation);
1547  }
1548  }
1549 
1550  std::sort(
1551  fPeakDev.begin(),
1552  fPeakDev.end(),
1553  [](std::tuple<double, int, int, int> const& t1, std::tuple<double, int, int, int> const& t2) {
1554  return std::get<0>(t1) > std::get<0>(t2);
1555  });
1556  Exponentials.Delete();
1557  }
1558 
1559  //---------------------------------------------------------------------------------------------
1560  std::string hit::DPRawHitFinder::CreateFitFunction(int fNPeaks, bool fSameShape)
1561  {
1562  std::string feqn = ""; // string holding fit formula
1563  std::stringstream numConv;
1564 
1565  if (fSameShape) {
1566  for (int i = 0; i < fNPeaks; i++) {
1567  feqn.append("+( [");
1568  numConv.str("");
1569  numConv << 2 * (i + 1);
1570  feqn.append(numConv.str());
1571  feqn.append("] * exp(0.4*(x-[");
1572  numConv.str("");
1573  numConv << 2 * (i + 1) + 1;
1574  feqn.append(numConv.str());
1575  feqn.append("])/[");
1576  numConv.str("");
1577  numConv << 0;
1578  feqn.append(numConv.str());
1579  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1580  numConv.str("");
1581  numConv << 2 * (i + 1) + 1;
1582  feqn.append(numConv.str());
1583  feqn.append("])/[");
1584  numConv.str("");
1585  numConv << 1;
1586  feqn.append(numConv.str());
1587  feqn.append("]) ) )");
1588  }
1589  }
1590  else {
1591  for (int i = 0; i < fNPeaks; i++) {
1592  feqn.append("+( [");
1593  numConv.str("");
1594  numConv << 4 * i + 2;
1595  feqn.append(numConv.str());
1596  feqn.append("] * exp(0.4*(x-[");
1597  numConv.str("");
1598  numConv << 4 * i + 3;
1599  feqn.append(numConv.str());
1600  feqn.append("])/[");
1601  numConv.str("");
1602  numConv << 4 * i;
1603  feqn.append(numConv.str());
1604  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1605  numConv.str("");
1606  numConv << 2 * (i + 1) + 1;
1607  feqn.append(numConv.str());
1608  feqn.append("])/[");
1609  numConv.str("");
1610  numConv << 4 * i + 1;
1611  feqn.append(numConv.str());
1612  feqn.append("]) ) )");
1613  }
1614  }
1615  return feqn;
1616  }
1617 
1618  //---------------------------------------------------------------------------------------------
1619  void hit::DPRawHitFinder::AddPeak(std::tuple<double, int, int, int> fPeakDevCand,
1620  PeakTimeWidVec& fpeakValsTemp)
1621  {
1622  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1623  int NewPeakMax = std::get<2>(fPeakDevCand);
1624  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1625  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1626  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1627 
1628  int NewPeakStart = 0;
1629  int NewPeakEnd = 0;
1630  int OldPeakNewStart = 0;
1631  int OldPeakNewEnd = 0;
1632  int DistanceBwOldAndNewPeak = 0;
1633 
1634  if (NewPeakMax < OldPeakMax) {
1635  NewPeakStart = OldPeakOldStart;
1636  OldPeakNewEnd = OldPeakOldEnd;
1637  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1638  NewPeakEnd = NewPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1639  if (DistanceBwOldAndNewPeak % 2 == 0) NewPeakEnd -= 1;
1640  OldPeakNewStart = NewPeakEnd + 1;
1641  }
1642  else if (OldPeakMax < NewPeakMax) {
1643  NewPeakEnd = OldPeakOldEnd;
1644  OldPeakNewStart = OldPeakOldStart;
1645  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1646  OldPeakNewEnd = OldPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1647  if (DistanceBwOldAndNewPeak % 2 == 0) OldPeakNewEnd -= 1;
1648  NewPeakStart = OldPeakNewEnd + 1;
1649  }
1650  else if (OldPeakMax == NewPeakMax) {
1651  return;
1652  } //This shouldn't happen
1653 
1654  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1655  std::make_tuple(OldPeakMax, 0, OldPeakNewStart, OldPeakNewEnd);
1656  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1657  std::sort(
1658  fpeakValsTemp.begin(),
1659  fpeakValsTemp.end(),
1660  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1661  return std::get<0>(t1) < std::get<0>(t2);
1662  });
1663 
1664  return;
1665  }
1666 
1667  //---------------------------------------------------------------------------------------------
1668  void hit::DPRawHitFinder::SplitPeak(std::tuple<double, int, int, int> fPeakDevCand,
1669  PeakTimeWidVec& fpeakValsTemp)
1670  {
1671  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1672  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1673  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1674  int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1675 
1676  if (WidthOldPeakOld < 3) { return; } //can't split peaks with widths < 3
1677 
1678  int NewPeakMax = 0;
1679  int NewPeakStart = 0;
1680  int NewPeakEnd = 0;
1681  int OldPeakNewMax = 0;
1682  int OldPeakNewStart = 0;
1683  int OldPeakNewEnd = 0;
1684 
1685  OldPeakNewStart = OldPeakOldStart;
1686  NewPeakEnd = OldPeakOldEnd;
1687 
1688  OldPeakNewEnd = OldPeakNewStart + 0.5 * (WidthOldPeakOld + (WidthOldPeakOld % 2));
1689  NewPeakStart = OldPeakNewEnd + 1;
1690 
1691  int WidthOldPeakNew = OldPeakNewEnd - OldPeakNewStart;
1692  int WidthNewPeak = NewPeakEnd - NewPeakStart;
1693 
1694  OldPeakNewMax = OldPeakNewStart + 0.5 * (WidthOldPeakNew - (WidthOldPeakNew % 2));
1695  NewPeakMax = NewPeakStart + 0.5 * (WidthNewPeak - (WidthNewPeak % 2));
1696 
1697  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1698  std::make_tuple(OldPeakNewMax, 0, OldPeakNewStart, OldPeakNewEnd);
1699  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1700  std::sort(
1701  fpeakValsTemp.begin(),
1702  fpeakValsTemp.end(),
1703  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1704  return std::get<0>(t1) < std::get<0>(t2);
1705  });
1706 
1707  return;
1708  }
1709 
1710  //---------------------------------------------------------------------------------------------
1711  double hit::DPRawHitFinder::WidthFunc(double fPeakMean,
1712  double fPeakAmp,
1713  double fPeakTau1,
1714  double fPeakTau2,
1715  double fStartTime,
1716  double fEndTime,
1717  double fPeakMeanTrue)
1718  {
1719  double MaxValue = (fPeakAmp * exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau1)) /
1720  (1 + exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau2));
1721  double FuncValue = 0.;
1722  double HalfMaxLeftTime = 0.;
1723  double HalfMaxRightTime = 0.;
1724  double ReBin = 10.;
1725 
1726  //First iteration (+- 1 bin)
1727  for (double x = fPeakMeanTrue; x > fStartTime - 1000.; x--) {
1728  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1729  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1730  if (FuncValue < 0.5 * MaxValue) {
1731  HalfMaxLeftTime = x;
1732  break;
1733  }
1734  }
1735 
1736  for (double x = fPeakMeanTrue; x < fEndTime + 1000.; x++) {
1737  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1738  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1739  if (FuncValue < 0.5 * MaxValue) {
1740  HalfMaxRightTime = x;
1741  break;
1742  }
1743  }
1744 
1745  //Second iteration (+- 0.1 bin)
1746  for (double x = HalfMaxLeftTime + 1; x > HalfMaxLeftTime; x -= (1 / ReBin)) {
1747  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1748  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1749  if (FuncValue < 0.5 * MaxValue) {
1750  HalfMaxLeftTime = x;
1751  break;
1752  }
1753  }
1754 
1755  for (double x = HalfMaxRightTime - 1; x < HalfMaxRightTime; x += (1 / ReBin)) {
1756  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1757  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1758  if (FuncValue < 0.5 * MaxValue) {
1759  HalfMaxRightTime = x;
1760  break;
1761  }
1762  }
1763 
1764  //Third iteration (+- 0.01 bin)
1765  for (double x = HalfMaxLeftTime + 1 / ReBin; x > HalfMaxLeftTime; x -= 1 / (ReBin * ReBin)) {
1766  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1767  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1768  if (FuncValue < 0.5 * MaxValue) {
1769  HalfMaxLeftTime = x;
1770  break;
1771  }
1772  }
1773 
1774  for (double x = HalfMaxRightTime - 1 / ReBin; x < HalfMaxRightTime; x += 1 / (ReBin * ReBin)) {
1775  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1776  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1777  if (FuncValue < 0.5 * MaxValue) {
1778  HalfMaxRightTime = x;
1779  break;
1780  }
1781  }
1782 
1783  return HalfMaxRightTime - HalfMaxLeftTime;
1784  }
1785 
1786  //---------------------------------------------------------------------------------------------
1787  double hit::DPRawHitFinder::ChargeFunc(double fPeakMean,
1788  double fPeakAmp,
1789  double fPeakTau1,
1790  double fPeakTau2,
1791  double fChargeNormFactor,
1792  double fPeakMeanTrue)
1793 
1794  {
1795  double ChargeSum = 0.;
1796  double Charge = 0.;
1797  double ReBin = 10.;
1798 
1799  bool ChargeBigEnough = true;
1800  for (double x = fPeakMeanTrue - 1 / ReBin; ChargeBigEnough && x > fPeakMeanTrue - 1000.;
1801  x -= 1.) {
1802  for (double i = 0.; i > -1.; i -= (1 / ReBin)) {
1803  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1804  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1805  ChargeSum += Charge;
1806  }
1807  if (Charge < 0.01) ChargeBigEnough = false;
1808  }
1809 
1810  ChargeBigEnough = true;
1811  for (double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue + 1000.; x += 1.) {
1812  for (double i = 0.; i < 1.; i += (1 / ReBin)) {
1813  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1814  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1815  ChargeSum += Charge;
1816  }
1817  if (Charge < 0.01) ChargeBigEnough = false;
1818  }
1819 
1820  return ChargeSum * fChargeNormFactor / ReBin;
1821  }
1822 
1823  //---------------------------------------------------------------------------------------------
1824  void hit::DPRawHitFinder::doBinAverage(const std::vector<float>& inputVec,
1825  std::vector<float>& outputVec,
1826  size_t binsToAverage) const
1827  {
1828  size_t halfBinsToAverage(binsToAverage / 2);
1829 
1830  float runningSum(0.);
1831 
1832  for (size_t idx = 0; idx < halfBinsToAverage; idx++)
1833  runningSum += inputVec[idx];
1834 
1835  outputVec.resize(inputVec.size());
1836  std::vector<float>::iterator outputVecItr = outputVec.begin();
1837 
1838  // First pass through to build the erosion vector
1839  for (std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end();
1840  inputItr++) {
1841  size_t startOffset = std::distance(inputVec.begin(), inputItr);
1842  size_t stopOffset = std::distance(inputItr, inputVec.end());
1843  size_t count =
1844  std::min(2 * halfBinsToAverage,
1845  std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1846 
1847  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1848  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1849 
1850  *outputVecItr++ = runningSum / float(count);
1851  }
1852 
1853  return;
1854  }
1855 
1856  //---------------------------------------------------------------------------------------------
1857  void hit::DPRawHitFinder::reBin(const std::vector<float>& inputVec,
1858  std::vector<float>& outputVec,
1859  size_t nBinsToCombine) const
1860  {
1861  size_t nNewBins = inputVec.size() / nBinsToCombine;
1862 
1863  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1864 
1865  outputVec.resize(nNewBins, 0.);
1866 
1867  size_t outputBin = 0;
1868 
1869  for (size_t inputIdx = 0; inputIdx < inputVec.size();) {
1870  outputVec[outputBin] += inputVec[inputIdx++];
1871 
1872  if (inputIdx % nBinsToCombine == 0) outputBin++;
1873 
1874  if (outputBin > outputVec.size()) {
1875  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin
1876  << ", inputIdx = " << inputIdx << std::endl;
1877  break;
1878  }
1879  }
1880 
1881  return;
1882  }
1883 
1884  DEFINE_ART_MODULE(DPRawHitFinder)
1885 
1886 } // end of hit namespace
x
Float_t x
Definition: compare.C:6
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
hit::DPRawHitFinder::findCandidatePeaks
void findCandidatePeaks(std::vector< float >::const_iterator startItr, std::vector< float >::const_iterator stopItr, TimeValsVec &timeValsVec, float &PeakMin, int firstTick) const
Definition: DPRawHitFinder_module.cc:1047
t1
TTree * t1
Definition: plottest35.C:26
hit::DPRawHitFinder::mergeCandidatePeaks
void mergeCandidatePeaks(const std::vector< float > signalVec, TimeValsVec, MergedTimeWidVec &)
Definition: DPRawHitFinder_module.cc:1123
hit::DPRawHitFinder::Comp::operator()
bool operator()(std::tuple< int, int, int, int > p, int s) const
Definition: DPRawHitFinder_module.cc:157
Hit.h
Declaration of signal hit object.
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
hit::DPRawHitFinder::produce
void produce(art::Event &evt) override
Definition: DPRawHitFinder_module.cc:292
hit::DPRawHitFinder::FindPeakWithMaxDeviation
void FindPeakWithMaxDeviation(const std::vector< float > fSignalVector, int fNPeaks, int fStartTime, int fEndTime, bool fSameShape, ParameterVec fparamVec, PeakTimeWidVec fpeakVals, PeakDevVec &fPeakDev)
Definition: DPRawHitFinder_module.cc:1495
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
lar::sparse_vector::get_ranges
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
Definition: sparse_vector.h:739
Get
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
hit::DPRawHitFinder::SplitPeak
void SplitPeak(std::tuple< double, int, int, int > fPeakDevCand, PeakTimeWidVec &fpeakValsTemp)
Definition: DPRawHitFinder_module.cc:1668
recob::HitAndAssociationsWriterBase::declare_products
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
hit::DPRawHitFinder::DPRawHitFinder
DPRawHitFinder(fhicl::ParameterSet const &pset)
Definition: DPRawHitFinder_module.cc:207
anab::FVectorWriter::addVector
void addVector(FVector_ID id, std::array< float, N > const &values)
Definition: MVAWriter.h:96
hit::DPRawHitFinder::TimeValsVec
std::vector< std::tuple< int, int, int >> TimeValsVec
Definition: DPRawHitFinder_module.cc:77
raw::TDCtick_t
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
recob::HitCreator
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
hit::DPRawHitFinder::FillOutHitParameterVector
void FillOutHitParameterVector(const std::vector< double > &input, std::vector< double > &output)
Definition: DPRawHitFinder_module.cc:260
HitCreator.h
Helper functions to create a hit.
anab::FVectorWriter::initOutputs
FVector_ID initOutputs(std::string const &dataTag, size_t dataSize, std::vector< std::string > const &names=std::vector< std::string >(N,""))
raw::InvalidChannelID
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:31
hit::DPRawHitFinder::AddPeak
void AddPeak(std::tuple< double, int, int, int > fPeakDevCand, PeakTimeWidVec &fpeakValsTemp)
Definition: DPRawHitFinder_module.cc:1619
recob::HitCollectionCreator
A class handling a collection of hits and its associations.
Definition: HitCreator.h:489
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
recob::Wire::Channel
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:223
hit::DPRawHitFinder::fHitParamWriter
anab::FVectorWriter< 4 > fHitParamWriter
Definition: DPRawHitFinder_module.cc:198
hit::DPRawHitFinder::ParameterVec
std::vector< std::pair< double, double >> ParameterVec
Definition: DPRawHitFinder_module.cc:86
Event.h
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1560
hit::DPRawHitFinder::PeakTimeWidVec
std::vector< std::tuple< int, int, int, int >> PeakTimeWidVec
Definition: DPRawHitFinder_module.cc:79
hit::DPRawHitFinder::PeakDevVec
std::vector< std::tuple< double, int, int, int >> PeakDevVec
Definition: DPRawHitFinder_module.cc:85
util::quantities::tick
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
recob::Wire::SignalROI
const RegionsOfInterest_t & SignalROI() const
Returns the list of regions of interest.
Definition: Wire.h:211
hit::DPRawHitFinder::ChargeFunc
double ChargeFunc(double fPeakMean, double fPeakAmp, double fPeakTau1, double fPeakTau2, double fChargeNormFactor, double fPeakMeanTrue)
Definition: DPRawHitFinder_module.cc:1787
recob::HitCollectionCreator::emplace_back
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
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
recob::HitCollectionCreator::put_into
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:622
anab::FVectorWriter::saveOutputs
void saveOutputs(art::Event &evt)
Check consistency and save all the results in the event.
Definition: MVAWriter.h:404
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
t2
TTree * t2
Definition: plottest35.C:36
art::ProductRegistryHelper::producesCollector
ProducesCollector & producesCollector() noexcept
Definition: ProductRegistryHelper.h:80
hit::DPRawHitFinder::FitExponentials
void FitExponentials(const std::vector< float > fSignalVector, const PeakTimeWidVec fPeakVals, int fStartTime, int fEndTime, ParameterVec &fparamVec, double &fchi2PerNDF, int &fNDF, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1276
lar::sparse_vector< float >
art::ProductRetriever::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: ProductRetriever.h:444
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
Wire.h
Declaration of basic channel signal object.
hit::DPRawHitFinder::EstimateFluctuations
int EstimateFluctuations(const std::vector< float > fsignalVec, int peakStart, int peakMean, int peakEnd)
Definition: DPRawHitFinder_module.cc:1249
hit::DPRawHitFinder::doBinAverage
void doBinAverage(const std::vector< float > &inputVec, std::vector< float > &outputVec, size_t binsToAverage) const
Definition: DPRawHitFinder_module.cc:1824
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
hit::DPRawHitFinder::MergedTimeWidVec
std::vector< std::tuple< int, int, PeakTimeWidVec, int >> MergedTimeWidVec
Definition: DPRawHitFinder_module.cc:84
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
recob::HitCreator::move
recob::Hit && move()
Prepares the constructed hit to be moved away.
Definition: HitCreator.h:338
art::Ptr
Definition: fwd.h:26
hit::DPRawHitFinder::WidthFunc
double WidthFunc(double fPeakMean, double fPeakAmp, double fPeakTau1, double fPeakTau2, double fStartTime, double fEndTime, double fPeakMeanTrue)
Definition: DPRawHitFinder_module.cc:1711
hit::DPRawHitFinder::reBin
void reBin(const std::vector< float > &inputVec, std::vector< float > &outputVec, size_t nBinsToCombine) const
Definition: DPRawHitFinder_module.cc:1857
hit::DPRawHitFinder::Comp::operator()
bool operator()(int s, std::tuple< int, int, int, int > p) const
Definition: DPRawHitFinder_module.cc:158