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