LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Classes | Public Types | Public Member Functions | Protected Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
hit::DPRawHitFinder Class Reference
Inheritance diagram for hit::DPRawHitFinder:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Classes

struct  Comp
 

Public Types

using ModuleType = EDProducer
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 DPRawHitFinder (fhicl::ParameterSet const &pset)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
void fillProductDescriptions ()
 
void registerProducts (ProductDescriptions &productsToRegister)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
std::unique_ptr< Worker > makeWorker (WorkerParams const &wp)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Types

using TimeValsVec = std::vector< std::tuple< int, int, int >>
 
using PeakTimeWidVec = std::vector< std::tuple< int, int, int, int >>
 
using MergedTimeWidVec = std::vector< std::tuple< int, int, PeakTimeWidVec, int >>
 
using PeakDevVec = std::vector< std::tuple< double, int, int, int >>
 
using ParameterVec = std::vector< std::pair< double, double >>
 

Private Member Functions

void produce (art::Event &evt) override
 
void beginJob () override
 
void findCandidatePeaks (std::vector< float >::const_iterator startItr, std::vector< float >::const_iterator stopItr, TimeValsVec &timeValsVec, float &PeakMin, int firstTick) const
 
int EstimateFluctuations (const std::vector< float > fsignalVec, int peakStart, int peakMean, int peakEnd)
 
void mergeCandidatePeaks (const std::vector< float > signalVec, TimeValsVec, MergedTimeWidVec &)
 
void FitExponentials (const std::vector< float > fSignalVector, const PeakTimeWidVec fPeakVals, int fStartTime, int fEndTime, ParameterVec &fparamVec, double &fchi2PerNDF, int &fNDF, bool fSameShape)
 
void FindPeakWithMaxDeviation (const std::vector< float > fSignalVector, int fNPeaks, int fStartTime, int fEndTime, bool fSameShape, ParameterVec fparamVec, PeakTimeWidVec fpeakVals, PeakDevVec &fPeakDev)
 
std::string CreateFitFunction (int fNPeaks, bool fSameShape)
 
void AddPeak (std::tuple< double, int, int, int > fPeakDevCand, PeakTimeWidVec &fpeakValsTemp)
 
void SplitPeak (std::tuple< double, int, int, int > fPeakDevCand, PeakTimeWidVec &fpeakValsTemp)
 
double WidthFunc (double fPeakMean, double fPeakAmp, double fPeakTau1, double fPeakTau2, double fStartTime, double fEndTime, double fPeakMeanTrue)
 
double ChargeFunc (double fPeakMean, double fPeakAmp, double fPeakTau1, double fPeakTau2, double fChargeNormFactor, double fPeakMeanTrue)
 
void FillOutHitParameterVector (const std::vector< double > &input, std::vector< double > &output)
 
void doBinAverage (const std::vector< float > &inputVec, std::vector< float > &outputVec, size_t binsToAverage) const
 
void reBin (const std::vector< float > &inputVec, std::vector< float > &outputVec, size_t nBinsToCombine) const
 

Private Attributes

std::string fCalDataModuleLabel
 
int fLogLevel
 
float fMinSig
 
int fTicksToStopPeakFinder
 
int fMinWidth
 
double fMinADCSum
 
double fMinADCSumOverWidth
 
int fMaxMultiHit
 
int fMaxGroupLength
 
double fChargeNorm
 
bool fTryNplus1Fits
 
double fChi2NDFRetry
 
double fChi2NDFRetryFactorMultiHits
 
double fChi2NDFMax
 
double fChi2NDFMaxFactorMultiHits
 
size_t fNumBinsToAverage
 
double fMinTau
 
double fMaxTau
 
double fFitPeakMeanRange
 
int fGroupMaxDistance
 
double fGroupMinADC
 
bool fSameShape
 
bool fDoMergePeaks
 
double fMergeADCSumThreshold
 
double fMergeMaxADCThreshold
 
double fMinRelativePeakHeightLeft
 
double fMinRelativePeakHeightRight
 
double fMergeMaxADCLimit
 
double fWidthNormalization
 
int fLongMaxHits
 
int fLongPulseWidth
 
int fMaxFluctuations
 
art::InputTag fNewHitsTag
 
anab::FVectorWriter< 4 > fHitParamWriter
 
TH1F * fFirstChi2
 
TH1F * fChi2
 

Detailed Description

Definition at line 68 of file DPRawHitFinder_module.cc.

Member Typedef Documentation

using hit::DPRawHitFinder::MergedTimeWidVec = std::vector<std::tuple<int, int, PeakTimeWidVec, int>>
private

Definition at line 84 of file DPRawHitFinder_module.cc.

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 17 of file EDProducer.h.

using hit::DPRawHitFinder::ParameterVec = std::vector<std::pair<double, double>>
private

Definition at line 86 of file DPRawHitFinder_module.cc.

using hit::DPRawHitFinder::PeakDevVec = std::vector<std::tuple<double, int, int, int>>
private

Definition at line 85 of file DPRawHitFinder_module.cc.

using hit::DPRawHitFinder::PeakTimeWidVec = std::vector< std::tuple<int, int, int, int>>
private

Definition at line 79 of file DPRawHitFinder_module.cc.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 26 of file Producer.h.

using hit::DPRawHitFinder::TimeValsVec = std::vector<std::tuple<int, int, int>>
private

Definition at line 77 of file DPRawHitFinder_module.cc.

Constructor & Destructor Documentation

hit::DPRawHitFinder::DPRawHitFinder ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 207 of file DPRawHitFinder_module.cc.

References recob::HitAndAssociationsWriterBase::declare_products(), fCalDataModuleLabel, fChargeNorm, fChi2NDFMax, fChi2NDFMaxFactorMultiHits, fChi2NDFRetry, fChi2NDFRetryFactorMultiHits, fDoMergePeaks, fFitPeakMeanRange, fGroupMaxDistance, fGroupMinADC, fHitParamWriter, fLogLevel, fLongMaxHits, fLongPulseWidth, fMaxFluctuations, fMaxGroupLength, fMaxMultiHit, fMaxTau, fMergeADCSumThreshold, fMergeMaxADCLimit, fMergeMaxADCThreshold, fMinADCSum, fMinADCSumOverWidth, fMinRelativePeakHeightLeft, fMinRelativePeakHeightRight, fMinSig, fMinTau, fMinWidth, fNewHitsTag, fNumBinsToAverage, fSameShape, fTicksToStopPeakFinder, fTryNplus1Fits, fWidthNormalization, anab::FVectorWriter< N >::produces_using(), and art::ProductRegistryHelper::producesCollector().

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()
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
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::fHitParamWriter
anab::FVectorWriter< 4 > fHitParamWriter
Definition: DPRawHitFinder_module.cc:198
art::ProductRegistryHelper::producesCollector
ProducesCollector & producesCollector() noexcept
Definition: ProductRegistryHelper.h:80
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46

Member Function Documentation

void hit::DPRawHitFinder::AddPeak ( std::tuple< double, int, int, int >  fPeakDevCand,
PeakTimeWidVec &  fpeakValsTemp 
)
private

Definition at line 1641 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1643  {
1644  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1645  int NewPeakMax = std::get<2>(fPeakDevCand);
1646  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1647  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1648  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1649 
1650  int NewPeakStart = 0;
1651  int NewPeakEnd = 0;
1652  int OldPeakNewStart = 0;
1653  int OldPeakNewEnd = 0;
1654  int DistanceBwOldAndNewPeak = 0;
1655 
1656  if (NewPeakMax < OldPeakMax) {
1657  NewPeakStart = OldPeakOldStart;
1658  OldPeakNewEnd = OldPeakOldEnd;
1659  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1660  NewPeakEnd = NewPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1661  if (DistanceBwOldAndNewPeak % 2 == 0) NewPeakEnd -= 1;
1662  OldPeakNewStart = NewPeakEnd + 1;
1663  }
1664  else if (OldPeakMax < NewPeakMax) {
1665  NewPeakEnd = OldPeakOldEnd;
1666  OldPeakNewStart = OldPeakOldStart;
1667  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1668  OldPeakNewEnd = OldPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1669  if (DistanceBwOldAndNewPeak % 2 == 0) OldPeakNewEnd -= 1;
1670  NewPeakStart = OldPeakNewEnd + 1;
1671  }
1672  else if (OldPeakMax == NewPeakMax) {
1673  return;
1674  } //This shouldn't happen
1675 
1676  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1677  std::make_tuple(OldPeakMax, 0, OldPeakNewStart, OldPeakNewEnd);
1678  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1679  std::sort(
1680  fpeakValsTemp.begin(),
1681  fpeakValsTemp.end(),
1682  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1683  return std::get<0>(t1) < std::get<0>(t2);
1684  });
1685 
1686  return;
1687  }
t1
TTree * t1
Definition: plottest35.C:26
t2
TTree * t2
Definition: plottest35.C:36
void hit::DPRawHitFinder::beginJob ( )
overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 281 of file DPRawHitFinder_module.cc.

References fChi2, and fFirstChi2.

282  {
283  // get access to the TFile service
284  art::ServiceHandle<art::TFileService const> tfs;
285 
286  // ======================================
287  // === Hit Information for Histograms ===
288  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 5000);
289  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 5000);
290  }
double hit::DPRawHitFinder::ChargeFunc ( double  fPeakMean,
double  fPeakAmp,
double  fPeakTau1,
double  fPeakTau2,
double  fChargeNormFactor,
double  fPeakMeanTrue 
)
private

Definition at line 1809 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1816  {
1817  double ChargeSum = 0.;
1818  double Charge = 0.;
1819  double ReBin = 10.;
1820 
1821  bool ChargeBigEnough = true;
1822  for (double x = fPeakMeanTrue - 1 / ReBin; ChargeBigEnough && x > fPeakMeanTrue - 1000.;
1823  x -= 1.) {
1824  for (double i = 0.; i > -1.; i -= (1 / ReBin)) {
1825  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1826  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1827  ChargeSum += Charge;
1828  }
1829  if (Charge < 0.01) ChargeBigEnough = false;
1830  }
1831 
1832  ChargeBigEnough = true;
1833  for (double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue + 1000.; x += 1.) {
1834  for (double i = 0.; i < 1.; i += (1 / ReBin)) {
1835  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1836  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1837  ChargeSum += Charge;
1838  }
1839  if (Charge < 0.01) ChargeBigEnough = false;
1840  }
1841 
1842  return ChargeSum * fChargeNormFactor / ReBin;
1843  }
x
Float_t x
Definition: compare.C:6
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

62  {
63  return collector_.consumes<T, BT>(tag);
64  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::consumes
ProductToken< T > consumes(InputTag const &)
Definition: ConsumesCollector.h:50
ConsumesCollector & art::ModuleBase::consumesCollector ( )
protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

58  {
59  return collector_;
60  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

76  {
77  collector_.consumesMany<T, BT>();
78  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

69  {
70  return collector_.consumesView<T, BT>(tag);
71  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::consumesView
ViewToken< Element > consumesView(InputTag const &)
std::string hit::DPRawHitFinder::CreateFitFunction ( int  fNPeaks,
bool  fSameShape 
)
private

Definition at line 1582 of file DPRawHitFinder_module.cc.

Referenced by FindPeakWithMaxDeviation(), and FitExponentials().

1583  {
1584  std::string feqn = ""; // string holding fit formula
1585  std::stringstream numConv;
1586 
1587  if (fSameShape) {
1588  for (int i = 0; i < fNPeaks; i++) {
1589  feqn.append("+( [");
1590  numConv.str("");
1591  numConv << 2 * (i + 1);
1592  feqn.append(numConv.str());
1593  feqn.append("] * exp(0.4*(x-[");
1594  numConv.str("");
1595  numConv << 2 * (i + 1) + 1;
1596  feqn.append(numConv.str());
1597  feqn.append("])/[");
1598  numConv.str("");
1599  numConv << 0;
1600  feqn.append(numConv.str());
1601  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1602  numConv.str("");
1603  numConv << 2 * (i + 1) + 1;
1604  feqn.append(numConv.str());
1605  feqn.append("])/[");
1606  numConv.str("");
1607  numConv << 1;
1608  feqn.append(numConv.str());
1609  feqn.append("]) ) )");
1610  }
1611  }
1612  else {
1613  for (int i = 0; i < fNPeaks; i++) {
1614  feqn.append("+( [");
1615  numConv.str("");
1616  numConv << 4 * i + 2;
1617  feqn.append(numConv.str());
1618  feqn.append("] * exp(0.4*(x-[");
1619  numConv.str("");
1620  numConv << 4 * i + 3;
1621  feqn.append(numConv.str());
1622  feqn.append("])/[");
1623  numConv.str("");
1624  numConv << 4 * i;
1625  feqn.append(numConv.str());
1626  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1627  numConv.str("");
1628  numConv << 2 * (i + 1) + 1;
1629  feqn.append(numConv.str());
1630  feqn.append("])/[");
1631  numConv.str("");
1632  numConv << 4 * i + 1;
1633  feqn.append(numConv.str());
1634  feqn.append("]) ) )");
1635  }
1636  }
1637  return feqn;
1638  }
void art::detail::Producer::doBeginJob ( SharedResources const &  resources)
inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

23  {
24  setupQueues(resources);
25  ProcessingFrame const frame{ScheduleID{}};
26  beginJobWithFrame(frame);
27  }
art::detail::Producer::setupQueues
virtual void setupQueues(SharedResources const &)=0
art::detail::Producer::beginJobWithFrame
virtual void beginJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doBeginRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

66  {
67  auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
68  ProcessingFrame const frame{mc.scheduleID()};
69  beginRunWithFrame(r, frame);
70  r.commitProducts();
71  return true;
72  }
r
TRandom r
Definition: spectrum.C:23
art::detail::Producer::beginRunWithFrame
virtual void beginRunWithFrame(Run &, ProcessingFrame const &)=0
art::RangeSet::forRun
static RangeSet forRun(RunID)
Definition: RangeSet.cc:51
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

86  {
87  auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
88  ProcessingFrame const frame{mc.scheduleID()};
89  beginSubRunWithFrame(sr, frame);
90  sr.commitProducts();
91  return true;
92  }
art::detail::Producer::beginSubRunWithFrame
virtual void beginSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
art::RangeSet::forSubRun
static RangeSet forSubRun(SubRunID)
Definition: RangeSet.cc:57
void hit::DPRawHitFinder::doBinAverage ( const std::vector< float > &  inputVec,
std::vector< float > &  outputVec,
size_t  binsToAverage 
) const
private

Definition at line 1846 of file DPRawHitFinder_module.cc.

Referenced by produce().

1849  {
1850  size_t halfBinsToAverage(binsToAverage / 2);
1851 
1852  float runningSum(0.);
1853 
1854  for (size_t idx = 0; idx < halfBinsToAverage; idx++)
1855  runningSum += inputVec[idx];
1856 
1857  outputVec.resize(inputVec.size());
1858  std::vector<float>::iterator outputVecItr = outputVec.begin();
1859 
1860  // First pass through to build the erosion vector
1861  for (std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end();
1862  inputItr++) {
1863  size_t startOffset = std::distance(inputVec.begin(), inputItr);
1864  size_t stopOffset = std::distance(inputItr, inputVec.end());
1865  size_t count =
1866  std::min(2 * halfBinsToAverage,
1867  std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1868 
1869  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1870  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1871 
1872  *outputVecItr++ = runningSum / float(count);
1873  }
1874 
1875  return;
1876  }
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
void art::detail::Producer::doEndJob ( )
inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

31  {
32  ProcessingFrame const frame{ScheduleID{}};
33  endJobWithFrame(frame);
34  }
art::detail::Producer::endJobWithFrame
virtual void endJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doEndRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

76  {
77  auto r = rp.makeRun(mc, rp.seenRanges());
78  ProcessingFrame const frame{mc.scheduleID()};
79  endRunWithFrame(r, frame);
80  r.commitProducts();
81  return true;
82  }
r
TRandom r
Definition: spectrum.C:23
art::detail::Producer::endRunWithFrame
virtual void endRunWithFrame(Run &, ProcessingFrame const &)=0
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

96  {
97  auto sr = srp.makeSubRun(mc, srp.seenRanges());
98  ProcessingFrame const frame{mc.scheduleID()};
99  endSubRunWithFrame(sr, frame);
100  sr.commitProducts();
101  return true;
102  }
art::detail::Producer::endSubRunWithFrame
virtual void endSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
bool art::detail::Producer::doEvent ( EventPrincipal &  ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

110  {
111  auto e = ep.makeEvent(mc);
112  ++counts_run;
113  ProcessingFrame const frame{mc.scheduleID()};
114  produceWithFrame(e, frame);
115  e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
116  ++counts_passed;
117  return true;
118  }
art::detail::Producer::checkPutProducts_
bool const checkPutProducts_
Definition: Producer.h:70
e
Float_t e
Definition: plot.C:35
art::detail::Producer::produceWithFrame
virtual void produceWithFrame(Event &, ProcessingFrame const &)=0
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

45  {
46  ProcessingFrame const frame{ScheduleID{}};
47  respondToCloseInputFileWithFrame(fb, frame);
48  }
art::detail::Producer::respondToCloseInputFileWithFrame
virtual void respondToCloseInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

59  {
60  ProcessingFrame const frame{ScheduleID{}};
61  respondToCloseOutputFilesWithFrame(fb, frame);
62  }
art::detail::Producer::respondToCloseOutputFilesWithFrame
virtual void respondToCloseOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

38  {
39  ProcessingFrame const frame{ScheduleID{}};
40  respondToOpenInputFileWithFrame(fb, frame);
41  }
art::detail::Producer::respondToOpenInputFileWithFrame
virtual void respondToOpenInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

52  {
53  ProcessingFrame const frame{ScheduleID{}};
54  respondToOpenOutputFilesWithFrame(fb, frame);
55  }
art::detail::Producer::respondToOpenOutputFilesWithFrame
virtual void respondToOpenOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
int hit::DPRawHitFinder::EstimateFluctuations ( const std::vector< float >  fsignalVec,
int  peakStart,
int  peakMean,
int  peakEnd 
)
private

Definition at line 1268 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks().

1272  {
1273  int NFluctuations = 0;
1274 
1275  for (int j = peakMean - 1; j >= peakStart; j--) {
1276  if (fsignalVec[j] < 5)
1277  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1278 
1279  if (fsignalVec[j] > fsignalVec[j + 1]) { NFluctuations++; }
1280  }
1281 
1282  for (int j = peakMean + 1; j <= peakEnd; j++) {
1283  if (fsignalVec[j] < 5)
1284  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1285 
1286  if (fsignalVec[j] > fsignalVec[j - 1]) { NFluctuations++; }
1287  }
1288 
1289  return NFluctuations;
1290  }
void hit::DPRawHitFinder::FillOutHitParameterVector ( const std::vector< double > &  input,
std::vector< double > &  output 
)
private

Definition at line 260 of file DPRawHitFinder_module.cc.

References geo::GeometryCore::Nplanes().

262  {
263  if (input.size() == 0)
264  throw std::runtime_error(
265  "DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
266 
267  art::ServiceHandle<geo::Geometry const> geom;
268  const unsigned int N_PLANES = geom->Nplanes();
269 
270  if (input.size() == 1)
271  output.resize(N_PLANES, input[0]);
272  else if (input.size() == N_PLANES)
273  output = input;
274  else
275  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config "
276  "vector size !=1 and !=N_PLANES.");
277  }
geo::GeometryCore::Nplanes
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
void art::Modifier::fillProductDescriptions ( )
inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

11  {
12  ProductRegistryHelper::fillDescriptions(moduleDescription());
13  }
art::ProductRegistryHelper::fillDescriptions
void fillDescriptions(ModuleDescription const &md)
Definition: ProductRegistryHelper.cc:48
art::ModuleBase::moduleDescription
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
void hit::DPRawHitFinder::findCandidatePeaks ( std::vector< float >::const_iterator  startItr,
std::vector< float >::const_iterator  stopItr,
TimeValsVec &  timeValsVec,
float &  PeakMin,
int  firstTick 
) const
private

Definition at line 1065 of file DPRawHitFinder_module.cc.

References fTicksToStopPeakFinder.

Referenced by produce().

1070  {
1071  // Need a minimum number of ticks to do any work here
1072  if (std::distance(startItr, stopItr) > 4) {
1073  // Find the highest peak in the range given
1074  auto maxItr = std::max_element(startItr, stopItr);
1075 
1076  float maxValue = *maxItr;
1077  int maxTime = std::distance(startItr, maxItr);
1078 
1079  if (maxValue >= PeakMin) {
1080  // backwards to find first bin for this candidate hit
1081  auto firstItr = std::distance(startItr, maxItr) > 2 ? maxItr - 1 : startItr;
1082  bool PeakStartIsHere = true;
1083 
1084  while (firstItr != startItr) {
1085  //Check for inflection point & ADC count <= 0
1086  PeakStartIsHere = true;
1087  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1088  if (*firstItr >= *(firstItr - i)) {
1089  PeakStartIsHere = false;
1090  break;
1091  }
1092  }
1093  if (*firstItr <= 0 || PeakStartIsHere) break;
1094  firstItr--;
1095  }
1096 
1097  int firstTime = std::distance(startItr, firstItr);
1098 
1099  // Recursive call to find all candidate hits earlier than this peak
1100  findCandidatePeaks(startItr, firstItr - 1, timeValsVec, PeakMin, firstTick);
1101 
1102  // forwards to find last bin for this candidate hit
1103  auto lastItr = std::distance(maxItr, stopItr) > 2 ? maxItr + 1 : stopItr - 1;
1104  bool PeakEndIsHere = true;
1105 
1106  while (lastItr != stopItr) {
1107  //Check for inflection point & ADC count <= 0
1108  PeakEndIsHere = true;
1109  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1110  if (*lastItr >= *(lastItr + i)) {
1111  PeakEndIsHere = false;
1112  break;
1113  }
1114  }
1115  if (*lastItr <= 0 || PeakEndIsHere) break;
1116  lastItr++;
1117  }
1118 
1119  int lastTime = std::distance(startItr, lastItr);
1120 
1121  // Now save this candidate's start and max time info
1122  timeValsVec.push_back(
1123  std::make_tuple(firstTick + firstTime, firstTick + maxTime, firstTick + lastTime));
1124 
1125  // Recursive call to find all candidate hits later than this peak
1126  findCandidatePeaks(lastItr + 1,
1127  stopItr,
1128  timeValsVec,
1129  PeakMin,
1130  firstTick + std::distance(startItr, lastItr + 1));
1131  }
1132  }
1133 
1134  return;
1135  }
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:1065
void hit::DPRawHitFinder::FindPeakWithMaxDeviation ( const std::vector< float >  fSignalVector,
int  fNPeaks,
int  fStartTime,
int  fEndTime,
bool  fSameShape,
ParameterVec  fparamVec,
PeakTimeWidVec  fpeakVals,
PeakDevVec &  fPeakDev 
)
private

Definition at line 1514 of file DPRawHitFinder_module.cc.

References CreateFitFunction(), t1, and t2.

Referenced by produce().

1522  {
1523  // int size = fEndTime - fStartTime + 1;
1524  // if(fEndTime - fStartTime < 0){size = 0;}
1525 
1526  std::string eqn =
1527  CreateFitFunction(fNPeaks, fSameShape); // string for equation of Exponentials fit
1528 
1529  TF1 Exponentials("Exponentials", eqn.c_str(), fStartTime, fEndTime + 1);
1530 
1531  for (size_t i = 0; i < fparamVec.size(); i++) {
1532  Exponentials.SetParameter(i, fparamVec[i].first);
1533  }
1534 
1535  // ##########################################################################
1536  // ### Finding the peak with the max chi2 fit and signal ###
1537  // ##########################################################################
1538  double Chi2PerNDFPeak;
1539  double MaxPosDeviation;
1540  double MaxNegDeviation;
1541  int BinMaxPosDeviation;
1542  int BinMaxNegDeviation;
1543  for (int i = 0; i < fNPeaks; i++) {
1544  Chi2PerNDFPeak = 0.;
1545  MaxPosDeviation = 0.;
1546  MaxNegDeviation = 0.;
1547  BinMaxPosDeviation = 0;
1548  BinMaxNegDeviation = 0;
1549 
1550  for (int j = std::get<2>(fpeakVals.at(i)); j < std::get<3>(fpeakVals.at(i)) + 1; j++) {
1551  if ((Exponentials(j + 0.5) - fSignalVector[j]) > MaxPosDeviation &&
1552  j != std::get<0>(fpeakVals.at(i))) {
1553  MaxPosDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1554  BinMaxPosDeviation = j;
1555  }
1556  if ((Exponentials(j + 0.5) - fSignalVector[j]) < MaxNegDeviation &&
1557  j != std::get<0>(fpeakVals.at(i))) {
1558  MaxNegDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1559  BinMaxNegDeviation = j;
1560  }
1561  Chi2PerNDFPeak +=
1562  pow((Exponentials(j + 0.5) - fSignalVector[j]) / sqrt(fSignalVector[j]), 2);
1563  }
1564 
1565  if (BinMaxNegDeviation != 0) {
1566  Chi2PerNDFPeak /=
1567  static_cast<double>((std::get<3>(fpeakVals.at(i)) - std::get<2>(fpeakVals.at(i))));
1568  fPeakDev.emplace_back(Chi2PerNDFPeak, i, BinMaxNegDeviation, BinMaxPosDeviation);
1569  }
1570  }
1571 
1572  std::sort(
1573  fPeakDev.begin(),
1574  fPeakDev.end(),
1575  [](std::tuple<double, int, int, int> const& t1, std::tuple<double, int, int, int> const& t2) {
1576  return std::get<0>(t1) > std::get<0>(t2);
1577  });
1578  Exponentials.Delete();
1579  }
t1
TTree * t1
Definition: plottest35.C:26
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1582
t2
TTree * t2
Definition: plottest35.C:36
void hit::DPRawHitFinder::FitExponentials ( const std::vector< float >  fSignalVector,
const PeakTimeWidVec  fPeakVals,
int  fStartTime,
int  fEndTime,
ParameterVec &  fparamVec,
double &  fchi2PerNDF,
int &  fNDF,
bool  fSameShape 
)
private

Definition at line 1295 of file DPRawHitFinder_module.cc.

References CreateFitFunction(), fFitPeakMeanRange, fLogLevel, fMaxTau, fMinTau, and util::size().

Referenced by produce().

1303  {
1304  int size = fEndTime - fStartTime + 1;
1305  int NPeaks = fPeakVals.size();
1306 
1307  // #############################################
1308  // ### If size < 0 then set the size to zero ###
1309  // #############################################
1310  if (fEndTime - fStartTime < 0) { size = 0; }
1311 
1312  // --- TH1D HitSignal ---
1313  TH1F hitSignal("hitSignal", "", std::max(size, 1), fStartTime, fEndTime + 1);
1314  hitSignal.Sumw2();
1315 
1316  // #############################
1317  // ### Filling the histogram ###
1318  // #############################
1319  for (int i = fStartTime; i < fEndTime + 1; i++) {
1320  hitSignal.Fill(i, fSignalVector[i]);
1321  hitSignal.SetBinError(i, 0.288675); //1/sqrt(12)
1322  }
1323 
1324  // ################################################
1325  // ### Building TFormula for basic Exponentials ###
1326  // ################################################
1327  std::string eqn = CreateFitFunction(NPeaks, fSameShape);
1328 
1329  // -------------------------------------
1330  // --- TF1 function for Exponentials ---
1331  // -------------------------------------
1332 
1333  TF1 Exponentials("Exponentials", eqn.c_str(), fStartTime, fEndTime + 1);
1334 
1335  if (fLogLevel >= 4) {
1336  std::cout << std::endl;
1337  std::cout << "--- Preparing fit ---" << std::endl;
1338  std::cout << "--- Lower limits, seed, upper limit:" << std::endl;
1339  }
1340 
1341  if (fSameShape) {
1342  Exponentials.SetParameter(0, 0.5);
1343  Exponentials.SetParameter(1, 0.5);
1344  Exponentials.SetParLimits(0, fMinTau, fMaxTau);
1345  Exponentials.SetParLimits(1, fMinTau, fMaxTau);
1346  double amplitude = 0;
1347  double peakMean = 0;
1348 
1349  double peakMeanShift = 2;
1350  double peakMeanSeed = 0;
1351  double peakMeanRangeLow = 0;
1352  double peakMeanRangeHi = 0;
1353  double peakStart = 0;
1354  double peakEnd = 0;
1355 
1356  for (int i = 0; i < NPeaks; i++) {
1357  peakMean = std::get<0>(fPeakVals.at(i));
1358  peakStart = std::get<2>(fPeakVals.at(i));
1359  peakEnd = std::get<3>(fPeakVals.at(i));
1360  peakMeanSeed = peakMean - peakMeanShift;
1361  peakMeanRangeLow = std::max(peakStart - peakMeanShift, peakMeanSeed - fFitPeakMeanRange);
1362  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed + fFitPeakMeanRange);
1363  amplitude = fSignalVector[peakMean];
1364 
1365  Exponentials.SetParameter(2 * (i + 1), 1.65 * amplitude);
1366  Exponentials.SetParLimits(2 * (i + 1), 0.3 * 1.65 * amplitude, 2 * 1.65 * amplitude);
1367  Exponentials.SetParameter(2 * (i + 1) + 1, peakMeanSeed);
1368 
1369  if (NPeaks == 1) {
1370  Exponentials.SetParLimits(2 * (i + 1) + 1, peakMeanRangeLow, peakMeanRangeHi);
1371  }
1372  else if (NPeaks >= 2 && i == 0) {
1373  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1374  Exponentials.SetParLimits(
1375  2 * (i + 1) + 1,
1376  peakMeanRangeLow,
1377  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1378  }
1379  else if (NPeaks >= 2 && i == NPeaks - 1) {
1380  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1381  Exponentials.SetParLimits(
1382  2 * (i + 1) + 1,
1383  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1384  peakMeanRangeHi);
1385  }
1386  else {
1387  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1388  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1389  Exponentials.SetParLimits(
1390  2 * (i + 1) + 1,
1391  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1392  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1393  }
1394 
1395  if (fLogLevel >= 4) {
1396  double t0low, t0high;
1397  Exponentials.GetParLimits(2 * (i + 1) + 1, t0low, t0high);
1398  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3 * 1.65 * amplitude << " , "
1399  << 1.65 * amplitude << " , " << 2 * 1.65 * amplitude << std::endl;
1400  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed
1401  << " , " << t0high << std::endl;
1402  }
1403  }
1404  }
1405  else {
1406  double amplitude = 0;
1407  double peakMean = 0;
1408 
1409  double peakMeanShift = 2;
1410  double peakMeanSeed = 0;
1411  double peakMeanRangeLow = 0;
1412  double peakMeanRangeHi = 0;
1413  double peakStart = 0;
1414  double peakEnd = 0;
1415 
1416  for (int i = 0; i < NPeaks; i++) {
1417  Exponentials.SetParameter(4 * i, 0.5);
1418  Exponentials.SetParameter(4 * i + 1, 0.5);
1419  Exponentials.SetParLimits(4 * i, fMinTau, fMaxTau);
1420  Exponentials.SetParLimits(4 * i + 1, fMinTau, fMaxTau);
1421 
1422  peakMean = std::get<0>(fPeakVals.at(i));
1423  peakStart = std::get<2>(fPeakVals.at(i));
1424  peakEnd = std::get<3>(fPeakVals.at(i));
1425  peakMeanSeed = peakMean - peakMeanShift;
1426  peakMeanRangeLow = std::max(peakStart - peakMeanShift, peakMeanSeed - fFitPeakMeanRange);
1427  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed + fFitPeakMeanRange);
1428  amplitude = fSignalVector[peakMean];
1429 
1430  Exponentials.SetParameter(4 * i + 2, 1.65 * amplitude);
1431  Exponentials.SetParLimits(4 * i + 2, 0.3 * 1.65 * amplitude, 2 * 1.65 * amplitude);
1432  Exponentials.SetParameter(4 * i + 3, peakMeanSeed);
1433 
1434  if (NPeaks == 1) {
1435  Exponentials.SetParLimits(4 * i + 3, peakMeanRangeLow, peakMeanRangeHi);
1436  }
1437  else if (NPeaks >= 2 && i == 0) {
1438  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1439  Exponentials.SetParLimits(
1440  4 * i + 3,
1441  peakMeanRangeLow,
1442  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1443  }
1444  else if (NPeaks >= 2 && i == NPeaks - 1) {
1445  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1446  Exponentials.SetParLimits(
1447  4 * i + 3,
1448  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1449  peakMeanRangeHi);
1450  }
1451  else {
1452  double HalfDistanceToNextMean = 0.5 * (std::get<0>(fPeakVals.at(i + 1)) - peakMean);
1453  double HalfDistanceToPrevMean = 0.5 * (peakMean - std::get<0>(fPeakVals.at(i - 1)));
1454  Exponentials.SetParLimits(
1455  4 * i + 3,
1456  std::max(peakMeanRangeLow, peakMeanSeed - HalfDistanceToPrevMean),
1457  std::min(peakMeanRangeHi, peakMeanSeed + HalfDistanceToNextMean));
1458  }
1459 
1460  if (fLogLevel >= 4) {
1461  double t0low, t0high;
1462  Exponentials.GetParLimits(4 * i + 3, t0low, t0high);
1463  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3 * 1.65 * amplitude << " , "
1464  << 1.65 * amplitude << " , " << 2 * 1.65 * amplitude << std::endl;
1465  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed
1466  << " , " << t0high << std::endl;
1467  }
1468  }
1469  }
1470 
1471  // ###########################################
1472  // ### PERFORMING THE TOTAL FIT OF THE HIT ###
1473  // ###########################################
1474  try {
1475  hitSignal.Fit(&Exponentials, "QNRWM", "", fStartTime, fEndTime + 1);
1476  }
1477  catch (...) {
1478  mf::LogWarning("DPRawHitFinder") << "Fitter failed finding a hit";
1479  }
1480 
1481  // ##################################################
1482  // ### Getting the fitted parameters from the fit ###
1483  // ##################################################
1484  fchi2PerNDF = (Exponentials.GetChisquare() / Exponentials.GetNDF());
1485  fNDF = Exponentials.GetNDF();
1486 
1487  if (fSameShape) {
1488  fparamVec.emplace_back(Exponentials.GetParameter(0), Exponentials.GetParError(0));
1489  fparamVec.emplace_back(Exponentials.GetParameter(1), Exponentials.GetParError(1));
1490 
1491  for (int i = 0; i < NPeaks; i++) {
1492  fparamVec.emplace_back(Exponentials.GetParameter(2 * (i + 1)),
1493  Exponentials.GetParError(2 * (i + 1)));
1494  fparamVec.emplace_back(Exponentials.GetParameter(2 * (i + 1) + 1),
1495  Exponentials.GetParError(2 * (i + 1) + 1));
1496  }
1497  }
1498  else {
1499  for (int i = 0; i < NPeaks; i++) {
1500  fparamVec.emplace_back(Exponentials.GetParameter(4 * i), Exponentials.GetParError(4 * i));
1501  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 1),
1502  Exponentials.GetParError(4 * i + 1));
1503  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 2),
1504  Exponentials.GetParError(4 * i + 2));
1505  fparamVec.emplace_back(Exponentials.GetParameter(4 * i + 3),
1506  Exponentials.GetParError(4 * i + 3));
1507  }
1508  }
1509  Exponentials.Delete();
1510  hitSignal.Delete();
1511  } //<----End FitExponentials
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1582
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const
inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

44  {
45  return collector_.getConsumables();
46  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::getConsumables
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables() const
Definition: ConsumesCollector.cc:13
std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const &  wp)
inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

38  {
39  return doMakeWorker(wp);
40  }
art::ModuleBase::doMakeWorker
virtual std::unique_ptr< Worker > doMakeWorker(WorkerParams const &wp)=0
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

83  {
84  return collector_.mayConsume<T, BT>(tag);
85  }
art::ConsumesCollector::mayConsume
ProductToken< T > mayConsume(InputTag const &)
Definition: ConsumesCollector.h:75
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

97  {
98  collector_.mayConsumeMany<T, BT>();
99  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

90  {
91  return collector_.mayConsumeView<T, BT>(tag);
92  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::mayConsumeView
ViewToken< Element > mayConsumeView(InputTag const &)
void hit::DPRawHitFinder::mergeCandidatePeaks ( const std::vector< float >  signalVec,
TimeValsVec  timeValsVec,
MergedTimeWidVec &  mergedVec 
)
private

Definition at line 1141 of file DPRawHitFinder_module.cc.

References EstimateFluctuations(), fDoMergePeaks, fGroupMaxDistance, fGroupMinADC, fMergeADCSumThreshold, fMergeMaxADCLimit, fMergeMaxADCThreshold, fMinRelativePeakHeightLeft, and fMinRelativePeakHeightRight.

Referenced by produce().

1144  {
1145  // ################################################################
1146  // ### Lets loop over the candidate pulses we found in this ROI ###
1147  // ################################################################
1148  auto timeValsVecItr = timeValsVec.begin();
1149  unsigned int PeaksInThisMergedPeak = 0;
1150  //int EndTickOfPreviousMergedPeak=0;
1151  while (timeValsVecItr != timeValsVec.end()) {
1152  PeakTimeWidVec peakVals;
1153 
1154  // Setting the start, peak, and end time of the pulse
1155  auto& timeVal = *timeValsVecItr++;
1156  int startT = std::get<0>(timeVal);
1157  int maxT = std::get<1>(timeVal);
1158  int endT = std::get<2>(timeVal);
1159  int widT = endT - startT;
1160  int FinalStartT = startT;
1161  int FinalEndT = endT;
1162 
1163  int NFluctuations = EstimateFluctuations(signalVec, startT, maxT, endT);
1164 
1165  peakVals.emplace_back(maxT, widT, startT, endT);
1166 
1167  // See if we want to merge pulses together
1168  // First check if we have more than one pulse on the wire
1169  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1170 
1171  // Loop until no more merged pulses (or candidates in this ROI)
1172  while (checkNextHit) {
1173  // 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
1174  int NextStartT = std::get<0>(*timeValsVecItr);
1175 
1176  double MinADC = signalVec[endT];
1177  for (int i = endT; i <= NextStartT; i++) {
1178  if (signalVec[i] < MinADC) { MinADC = signalVec[i]; }
1179  }
1180 
1181  // Group peaks (grouped peaks are fitted together and can be merged)
1182  if (MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance) {
1183  int CurrentStartT = startT;
1184  int CurrentMaxT = maxT;
1185  int CurrentEndT = endT;
1186  //int CurrentWidT=widT;
1187 
1188  timeVal = *timeValsVecItr++;
1189  int NextMaxT = std::get<1>(timeVal);
1190  int NextEndT = std::get<2>(timeVal);
1191  int NextWidT = NextEndT - NextStartT;
1192  FinalEndT = NextEndT;
1193 
1194  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1195 
1196  int CurrentSumADC = 0;
1197  for (int i = CurrentStartT; i <= CurrentEndT; i++) {
1198  CurrentSumADC += signalVec[i];
1199  }
1200 
1201  int NextSumADC = 0;
1202  for (int i = NextStartT; i <= NextEndT; i++) {
1203  NextSumADC += signalVec[i];
1204  }
1205 
1206  //Merge peaks within a group
1207  if (fDoMergePeaks) {
1208  if (signalVec[NextMaxT] <= signalVec[CurrentMaxT] &&
1209  ((signalVec[NextMaxT] < fMergeMaxADCThreshold * signalVec[CurrentMaxT] &&
1210  NextSumADC < fMergeADCSumThreshold * CurrentSumADC) ||
1211  (signalVec[NextMaxT] - signalVec[NextStartT]) <
1212  fMinRelativePeakHeightRight * signalVec[CurrentMaxT]) &&
1213  (signalVec[NextMaxT] <= fMergeMaxADCLimit)) {
1214  maxT = CurrentMaxT;
1215  startT = CurrentStartT;
1216  endT = NextEndT;
1217  widT = endT - startT;
1218  peakVals.pop_back();
1219  peakVals.emplace_back(maxT, widT, startT, endT);
1220  }
1221  else if (signalVec[NextMaxT] > signalVec[CurrentMaxT] &&
1222  ((signalVec[CurrentMaxT] < fMergeMaxADCThreshold * signalVec[NextMaxT] &&
1223  CurrentSumADC < fMergeADCSumThreshold * NextSumADC) ||
1224  (signalVec[CurrentMaxT] - signalVec[CurrentEndT]) <
1225  fMinRelativePeakHeightLeft * signalVec[NextMaxT]) &&
1226  (signalVec[CurrentMaxT] <= fMergeMaxADCLimit)) {
1227  maxT = NextMaxT;
1228  startT = CurrentStartT;
1229  endT = NextEndT;
1230  widT = endT - startT;
1231  peakVals.pop_back();
1232  peakVals.emplace_back(maxT, widT, startT, endT);
1233  }
1234  else {
1235  maxT = NextMaxT;
1236  startT = NextStartT;
1237  endT = NextEndT;
1238  widT = NextWidT;
1239  peakVals.emplace_back(maxT, widT, startT, endT);
1240  PeaksInThisMergedPeak++;
1241  }
1242  }
1243  else {
1244  maxT = NextMaxT;
1245  startT = NextStartT;
1246  endT = NextEndT;
1247  widT = NextWidT;
1248  peakVals.emplace_back(maxT, widT, startT, endT);
1249  PeaksInThisMergedPeak++;
1250  }
1251  checkNextHit = timeValsVecItr != timeValsVec.end();
1252  } //<---Checking adjacent pulses
1253  else {
1254  checkNextHit = false;
1255  PeaksInThisMergedPeak = 0;
1256  }
1257 
1258  } //<---End checking if there is more than one pulse on the wire
1259  // Add these to our merged vector
1260  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1261  }
1262  return;
1263  }
hit::DPRawHitFinder::PeakTimeWidVec
std::vector< std::tuple< int, int, int, int >> PeakTimeWidVec
Definition: DPRawHitFinder_module.cc:79
hit::DPRawHitFinder::EstimateFluctuations
int EstimateFluctuations(const std::vector< float > fsignalVec, int peakStart, int peakMean, int peakEnd)
Definition: DPRawHitFinder_module.cc:1268
ModuleDescription const & art::ModuleBase::moduleDescription ( ) const
inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

14  {
15  if (md_.has_value()) {
16  return *md_;
17  }
18 
19  throw Exception{errors::LogicError,
20  "There was an error while calling moduleDescription().\n"}
21  << "The moduleDescription() base-class member function cannot be called\n"
22  "during module construction. To determine which module is "
23  "responsible\n"
24  "for calling it, find the '<module type>:<module "
25  "label>@Construction'\n"
26  "tag in the message prefix above. Please contact artists@fnal.gov\n"
27  "for guidance.\n";
28  }
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
art::ModuleBase::md_
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void hit::DPRawHitFinder::produce ( art::Event &  evt)
overrideprivatevirtual

Implements art::EDProducer.

Definition at line 293 of file DPRawHitFinder_module.cc.

References AddPeak(), anab::FVectorWriter< N >::addVector(), recob::Wire::Channel(), geo::GeometryCore::ChannelToWire(), ChargeFunc(), geo::CryostatID::Cryostat, doBinAverage(), recob::HitCollectionCreator::emplace_back(), fCalDataModuleLabel, fChargeNorm, fChi2, fChi2NDFMax, fChi2NDFMaxFactorMultiHits, fChi2NDFRetry, fChi2NDFRetryFactorMultiHits, fFirstChi2, fHitParamWriter, findCandidatePeaks(), FindPeakWithMaxDeviation(), FitExponentials(), fLogLevel, fLongMaxHits, fLongPulseWidth, fMaxFluctuations, fMaxGroupLength, fMaxMultiHit, fMinADCSum, fMinADCSumOverWidth, fMinSig, fMinWidth, fNewHitsTag, fNumBinsToAverage, fSameShape, fTryNplus1Fits, fWidthNormalization, lar::sparse_vector< T >::get_ranges(), art::ProductRetriever::getByLabel(), anab::FVectorWriter< N >::initOutputs(), raw::InvalidChannelID, mergeCandidatePeaks(), recob::HitCreator::move(), geo::PlaneID::Plane, recob::HitCollectionCreator::put_into(), anab::FVectorWriter< N >::saveOutputs(), recob::Wire::SignalROI(), SplitPeak(), geo::TPCID::TPC, WidthFunc(), and geo::WireID::Wire.

294  {
295  //==================================================================================================
296  TH1::AddDirectory(kFALSE);
297 
298  //Instantiate and Reset a stop watch
299  //TStopwatch StopWatch;
300  //StopWatch.Reset();
301 
302  // ################################
303  // ### Calling Geometry service ###
304  // ################################
305  art::ServiceHandle<geo::Geometry const> geom;
306 
307  // ###############################################
308  // ### Making a ptr vector to put on the event ###
309  // ###############################################
310  // this contains the hit collection
311  // and its associations to wires and raw digits
312  recob::HitCollectionCreator hcol(evt);
313 
314  // start collection of fit parameters, initialize metadata describing it
315  auto hitID =
316  fHitParamWriter.initOutputs<recob::Hit>(fNewHitsTag, {"t0", "tau1", "tau2", "ampl"});
317 
318  // ##########################################
319  // ### Reading in the Wire List object(s) ###
320  // ##########################################
321  art::Handle<std::vector<recob::Wire>> wireVecHandle;
322  evt.getByLabel(fCalDataModuleLabel, wireVecHandle);
323 
324  // #################################################################
325  // ### Reading in the RawDigit associated with these wires, too ###
326  // #################################################################
327  art::FindOneP<raw::RawDigit> RawDigits(wireVecHandle, evt, fCalDataModuleLabel);
328  // Channel Number
329  raw::ChannelID_t channel = raw::InvalidChannelID;
330 
331  //##############################
332  //### Looping over the wires ###
333  //##############################
334  for (size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
335  // ####################################
336  // ### Getting this particular wire ###
337  // ####################################
338  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
339  art::Ptr<raw::RawDigit> rawdigits = RawDigits.at(wireIter);
340  // --- Setting Channel Number and Signal type ---
341  channel = wire->Channel();
342  // get the WireID for this hit
343  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
344  // for now, just take the first option returned from ChannelToWire
345  geo::WireID wid = wids[0];
346 
347  if (fLogLevel >= 1) {
348  std::cout << std::endl;
349  std::cout << std::endl;
350  std::cout << std::endl;
351  std::cout << "-----------------------------------------------------------------------------"
352  "------------------------------"
353  << std::endl;
354  std::cout << "Channel: " << channel << std::endl;
355  std::cout << "Cryostat: " << wid.Cryostat << std::endl;
356  std::cout << "TPC: " << wid.TPC << std::endl;
357  std::cout << "Plane: " << wid.Plane << std::endl;
358  std::cout << "Wire: " << wid.Wire << std::endl;
359  }
360 
361  // #################################################
362  // ### Set up to loop over ROI's for this wire ###
363  // #################################################
364  const recob::Wire::RegionsOfInterest_t& signalROI = wire->SignalROI();
365 
366  int CountROI = 0;
367 
368  for (const auto& range : signalROI.get_ranges()) {
369  // #################################################
370  // ### Getting a vector of signals for this wire ###
371  // #################################################
372  const std::vector<float>& signal = range.data();
373 
374  // ROI start time
375  raw::TDCtick_t roiFirstBinTick = range.begin_index();
376  MergedTimeWidVec mergedVec;
377 
378  // ###########################################################
379  // ### If option set do bin averaging before finding peaks ###
380  // ###########################################################
381 
382  if (fNumBinsToAverage > 1) {
383  std::vector<float> timeAve;
384  doBinAverage(signal, timeAve, fNumBinsToAverage);
385 
386  // ###################################################################
387  // ### Search current averaged ROI for candidate peaks and widths ###
388  // ###################################################################
389  TimeValsVec timeValsVec;
390  findCandidatePeaks(timeAve.begin(), timeAve.end(), timeValsVec, fMinSig, 0);
391 
392  // ####################################################
393  // ### If no startTime hit was found skip this wire ###
394  // ####################################################
395  if (timeValsVec.empty()) continue;
396 
397  // #############################################################
398  // ### Merge potentially overlapping peaks and do multi fit ###
399  // #############################################################
400  mergeCandidatePeaks(timeAve, timeValsVec, mergedVec);
401  }
402 
403  // ###########################################################
404  // ### Otherwise, operate directonly on signal vector ###
405  // ###########################################################
406  else {
407  // ##########################################################
408  // ### Search current ROI for candidate peaks and widths ###
409  // ##########################################################
410  TimeValsVec timeValsVec;
411  findCandidatePeaks(signal.begin(), signal.end(), timeValsVec, fMinSig, 0);
412 
413  if (fLogLevel >= 1) {
414  std::cout << std::endl;
415  std::cout << std::endl;
416  std::cout << "-------------------- ROI #" << CountROI << " -------------------- "
417  << std::endl;
418  if (timeValsVec.size() == 1)
419  std::cout << "ROI #" << CountROI << " (" << timeValsVec.size()
420  << " peak): ROIStartTick: " << range.offset
421  << " ROIEndTick:" << range.offset + range.size() << std::endl;
422  else
423  std::cout << "ROI #" << CountROI << " (" << timeValsVec.size()
424  << " peaks): ROIStartTick: " << range.offset
425  << " ROIEndTick:" << range.offset + range.size() << std::endl;
426  CountROI++;
427  }
428 
429  if (fLogLevel >= 2) {
430  int CountPeak = 0;
431  for (auto const& timeValsTmp : timeValsVec) {
432  std::cout << "Peak #" << CountPeak
433  << ": PeakStartTick: " << range.offset + std::get<0>(timeValsTmp)
434  << " PeakMaxTick: " << range.offset + std::get<1>(timeValsTmp)
435  << " PeakEndTick: " << range.offset + std::get<2>(timeValsTmp)
436  << 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  // ### Creating the parameter vector for the new pulse ###
453  // #######################################################
454  ParameterVec paramVec;
455 
456  // === Number of Exponentials to try ===
457  int NumberOfPeaksBeforeFit = 0;
458  unsigned int nExponentialsForFit = 0;
459  double chi2PerNDF = 0.;
460  int NDF = 0;
461 
462  unsigned int NumberOfMergedVecs = mergedVec.size();
463 
464  // ################################################################
465  // ### Lets loop over the groups of peaks we found on this wire ###
466  // ################################################################
467 
468  for (unsigned int j = 0; j < NumberOfMergedVecs; j++) {
469  int startT = std::get<0>(mergedVec.at(j));
470  int endT = std::get<1>(mergedVec.at(j));
471  int width = endT + 1 - startT;
472  PeakTimeWidVec& peakVals = std::get<2>(mergedVec.at(j));
473 
474  int NFluctuations = std::get<3>(mergedVec.at(j));
475 
476  if (fLogLevel >= 3) {
477  std::cout << std::endl;
478  if (peakVals.size() == 1)
479  std::cout << "- Group #" << j << " (" << peakVals.size()
480  << " peak): GroupStartTick: " << range.offset + startT
481  << " GroupEndTick: " << range.offset + endT << std::endl;
482  else
483  std::cout << "- Group #" << j << " (" << peakVals.size()
484  << " peaks): GroupStartTick: " << range.offset + startT
485  << " GroupEndTick: " << range.offset + endT << std::endl;
486  std::cout << "Fluctuations in this group: " << NFluctuations << std::endl;
487  int CountPeakInGroup = 0;
488  for (auto const& peakValsTmp : peakVals) {
489  std::cout << "Peak #" << CountPeakInGroup << " in group #" << j
490  << ": PeakInGroupStartTick: " << range.offset + std::get<2>(peakValsTmp)
491  << " PeakInGroupMaxTick: " << range.offset + std::get<0>(peakValsTmp)
492  << " PeakInGroupEndTick: " << range.offset + std::get<3>(peakValsTmp)
493  << std::endl;
494  CountPeakInGroup++;
495  }
496  }
497 
498  // ### Getting rid of noise hits ###
499  if (width < fMinWidth ||
500  (double)std::accumulate(signal.begin() + startT, signal.begin() + endT + 1, 0) <
501  fMinADCSum ||
502  (double)std::accumulate(signal.begin() + startT, signal.begin() + endT + 1, 0) /
503  width <
504  fMinADCSumOverWidth) {
505  if (fLogLevel >= 3) {
506  std::cout << "Delete this group of peaks because width, integral or width/intergral "
507  "is too small."
508  << std::endl;
509  }
510  continue;
511  }
512 
513  // #####################################################################################################
514  // ### Only attempt to fit if number of peaks <= fMaxMultiHit and if group length <= fMaxGroupLength ###
515  // #####################################################################################################
516  NumberOfPeaksBeforeFit = peakVals.size();
517  nExponentialsForFit = peakVals.size();
518  chi2PerNDF = 0.;
519  NDF = 0;
520  if (NumberOfPeaksBeforeFit <= fMaxMultiHit && width <= fMaxGroupLength &&
521  NFluctuations <= fMaxFluctuations) {
522  // #####################################################
523  // ### Calling the function for fitting Exponentials ###
524  // #####################################################
525  paramVec.clear();
526  FitExponentials(signal, peakVals, startT, endT, paramVec, chi2PerNDF, NDF, fSameShape);
527 
528  if (fLogLevel >= 4) {
529  std::cout << std::endl;
530  std::cout << "--- First fit ---" << std::endl;
531  if (nExponentialsForFit == 1)
532  std::cout << "- Fitted " << nExponentialsForFit << " peak in group #" << j << ":"
533  << std::endl;
534  else
535  std::cout << "- Fitted " << nExponentialsForFit << " peaks in group #" << j << ":"
536  << std::endl;
537  std::cout << "chi2/ndf = " << chi2PerNDF << std::endl;
538 
539  if (fSameShape) {
540  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
541  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
542 
543  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
544  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2 * (i + 1)].first
545  << std::endl;
546  std::cout << "Peak #" << i
547  << ": t0 [ticks] = " << range.offset + paramVec[2 * (i + 1) + 1].first
548  << std::endl;
549  }
550  }
551  else {
552  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
553  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[4 * i + 2].first
554  << std::endl;
555  std::cout << "Peak #" << i
556  << ": t0 [ticks] = " << range.offset + paramVec[4 * i + 3].first
557  << std::endl;
558  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4 * i].first
559  << std::endl;
560  std::cout << "Peak #" << i << ": tau2 [mus] = " << paramVec[4 * i + 1].first
561  << std::endl;
562  }
563  }
564  }
565 
566  // If the chi2 is infinite then there is a real problem so we bail
567  if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) continue;
568 
569  fFirstChi2->Fill(chi2PerNDF);
570 
571  // ########################################################
572  // ### Trying extra Exponentials for an initial bad fit ###
573  // ########################################################
574 
575  if ((fTryNplus1Fits && nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFRetry) ||
576  (fTryNplus1Fits && nExponentialsForFit > 1 &&
577  chi2PerNDF > fChi2NDFRetryFactorMultiHits * fChi2NDFRetry)) {
578  unsigned int nExponentialsBeforeRefit = nExponentialsForFit;
579  unsigned int nExponentialsAfterRefit = nExponentialsForFit;
580  double oldChi2PerNDF = chi2PerNDF;
581  double chi2PerNDF2;
582  int NDF2;
583  bool RefitSuccess;
584  PeakTimeWidVec peakValsTemp;
585  while ((nExponentialsForFit == 1 &&
586  nExponentialsAfterRefit < 3 * nExponentialsBeforeRefit &&
587  chi2PerNDF > fChi2NDFRetry) ||
588  (nExponentialsForFit > 1 &&
589  nExponentialsAfterRefit < 3 * nExponentialsBeforeRefit &&
590  chi2PerNDF > fChi2NDFRetryFactorMultiHits * fChi2NDFRetry)) {
591  RefitSuccess = false;
592  PeakDevVec PeakDev;
593  FindPeakWithMaxDeviation(signal,
594  nExponentialsForFit,
595  startT,
596  endT,
597  fSameShape,
598  paramVec,
599  peakVals,
600  PeakDev);
601 
602  //Add peak and re-fit
603  for (auto& PeakDevCand : PeakDev) {
604  chi2PerNDF2 = 0.;
605  NDF2 = 0.;
606  ParameterVec paramVecRefit;
607  peakValsTemp = peakVals;
608 
609  AddPeak(PeakDevCand, peakValsTemp);
610  FitExponentials(signal,
611  peakValsTemp,
612  startT,
613  endT,
614  paramVecRefit,
615  chi2PerNDF2,
616  NDF2,
617  fSameShape);
618 
619  if (chi2PerNDF2 < chi2PerNDF) {
620  paramVec = paramVecRefit;
621  peakVals = peakValsTemp;
622  nExponentialsForFit = peakVals.size();
623  chi2PerNDF = chi2PerNDF2;
624  NDF = NDF2;
625  nExponentialsAfterRefit++;
626  RefitSuccess = true;
627  break;
628  }
629  }
630 
631  //Split peak and re-fit
632  if (RefitSuccess == false) {
633  for (auto& PeakDevCand : PeakDev) {
634  chi2PerNDF2 = 0.;
635  NDF2 = 0.;
636  ParameterVec paramVecRefit;
637  peakValsTemp = peakVals;
638 
639  SplitPeak(PeakDevCand, peakValsTemp);
640  FitExponentials(signal,
641  peakValsTemp,
642  startT,
643  endT,
644  paramVecRefit,
645  chi2PerNDF2,
646  NDF2,
647  fSameShape);
648 
649  if (chi2PerNDF2 < chi2PerNDF) {
650  paramVec = paramVecRefit;
651  peakVals = peakValsTemp;
652  nExponentialsForFit = peakVals.size();
653  chi2PerNDF = chi2PerNDF2;
654  NDF = NDF2;
655  nExponentialsAfterRefit++;
656  RefitSuccess = true;
657  break;
658  }
659  }
660  }
661 
662  if (RefitSuccess == false) { break; }
663  }
664 
665  if (fLogLevel >= 5) {
666  std::cout << std::endl;
667  std::cout << "--- Refit ---" << std::endl;
668  if (chi2PerNDF == oldChi2PerNDF)
669  std::cout << "chi2/ndf didn't improve. Keep first fit." << std::endl;
670  else {
671  std::cout << "- Added peaks to group #" << j << ". This group now has "
672  << nExponentialsForFit << " peaks:" << std::endl;
673  std::cout << "- Group #" << j << " (" << peakVals.size()
674  << " peaks): GroupStartTick: " << range.offset + startT
675  << " GroupEndTick: " << range.offset + endT << std::endl;
676 
677  int CountPeakInGroup = 0;
678  for (auto const& peakValsTmp : peakVals) {
679  std::cout << "Peak #" << CountPeakInGroup << " in group #" << j
680  << ": PeakInGroupStartTick: "
681  << range.offset + std::get<2>(peakValsTmp)
682  << " PeakInGroupMaxTick: "
683  << range.offset + std::get<0>(peakValsTmp)
684  << " PeakInGroupEndTick: "
685  << range.offset + std::get<3>(peakValsTmp) << std::endl;
686  CountPeakInGroup++;
687  }
688 
689  std::cout << "chi2/ndf = " << chi2PerNDF << std::endl;
690 
691  if (fSameShape) {
692  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
693  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
694 
695  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
696  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2 * (i + 1)].first
697  << std::endl;
698  std::cout << "Peak #" << i << ": t0 [ticks] = "
699  << range.offset + paramVec[2 * (i + 1) + 1].first << std::endl;
700  }
701  }
702  else {
703  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
704  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[4 * i + 2].first
705  << std::endl;
706  std::cout << "Peak #" << i
707  << ": t0 [ticks] = " << range.offset + paramVec[4 * i + 3].first
708  << std::endl;
709  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4 * i].first
710  << std::endl;
711  std::cout << "Peak #" << i << ": tau2 [mus] = " << paramVec[4 * i + 1].first
712  << std::endl;
713  }
714  }
715  }
716  }
717  }
718 
719  // #######################################################
720  // ### Loop through returned peaks and make recob hits ###
721  // #######################################################
722 
723  int numHits(0);
724  for (unsigned int i = 0; i < nExponentialsForFit; i++) {
725  //Extract fit parameters for this hit
726  double peakTau1;
727  double peakTau2;
728  double peakAmp;
729  double peakMean;
730 
731  if (fSameShape) {
732  peakTau1 = paramVec[0].first;
733  peakTau2 = paramVec[1].first;
734  peakAmp = paramVec[2 * (i + 1)].first;
735  peakMean = paramVec[2 * (i + 1) + 1].first;
736  }
737  else {
738  peakTau1 = paramVec[4 * i].first;
739  peakTau2 = paramVec[4 * i + 1].first;
740  peakAmp = paramVec[4 * i + 2].first;
741  peakMean = paramVec[4 * i + 3].first;
742  }
743 
744  //Highest ADC count in peak = peakAmpTrue
745  double peakAmpTrue = signal[std::get<0>(peakVals.at(i))];
746  double peakAmpErr = 1.;
747 
748  //Determine peak position of fitted function (= peakMeanTrue)
749  TF1 Exponentials("Exponentials",
750  "( [0] * exp(0.4*(x-[1])/[2]) / ( 1 + exp(0.4*(x-[1])/[3]) ) )",
751  startT,
752  endT);
753  Exponentials.SetParameter(0, peakAmp);
754  Exponentials.SetParameter(1, peakMean);
755  Exponentials.SetParameter(2, peakTau1);
756  Exponentials.SetParameter(3, peakTau2);
757  double peakMeanTrue = Exponentials.GetMaximumX(startT, endT);
758  Exponentials.Delete();
759 
760  //Calculate width (=FWHM)
761  double peakWidth =
762  WidthFunc(peakMean, peakAmp, peakTau1, peakTau2, startT, endT, peakMeanTrue);
763  peakWidth /=
764  fWidthNormalization; //from FWHM to "standard deviation": standard deviation = FWHM/(2*sqrt(2*ln(2)))
765 
766  // Extract fit parameter errors
767  double peakMeanErr;
768 
769  if (fSameShape) { peakMeanErr = paramVec[2 * (i + 1) + 1].second; }
770  else {
771  peakMeanErr = paramVec[4 * i + 3].second;
772  }
773  double peakWidthErr = 0.1 * peakWidth;
774 
775  // ### Charge ###
776  double charge =
777  ChargeFunc(peakMean, peakAmp, peakTau1, peakTau2, fChargeNorm, peakMeanTrue);
778  double chargeErr =
779  std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
780 
781  // ### limits for getting sum of ADC counts
782  int startTthisHit = std::get<2>(peakVals.at(i));
783  int endTthisHit = std::get<3>(peakVals.at(i));
784  std::vector<float>::const_iterator sumStartItr = signal.begin() + startTthisHit;
785  std::vector<float>::const_iterator sumEndItr = signal.begin() + endTthisHit;
786 
787  // ### Sum of ADC counts
788  double sumADC = std::accumulate(sumStartItr, sumEndItr + 1, 0.);
789 
790  //Check if fit returns reasonable values and ich chi2 is below threshold
791  if (peakWidth <= 0 || charge <= 0. || charge != charge ||
792  (nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax) ||
793  (nExponentialsForFit >= 2 &&
794  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)) {
795  if (fLogLevel >= 1) {
796  std::cout << std::endl;
797  std::cout << "WARNING: For peak #" << i << " in this group:" << std::endl;
798  if (peakWidth <= 0 || charge <= 0. || charge != charge)
799  std::cout << "Fit function returned width < 0 or charge < 0 or charge = nan."
800  << std::endl;
801  if ((nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax) ||
802  (nExponentialsForFit >= 2 &&
803  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)) {
804  std::cout << std::endl;
805  std::cout << "WARNING: For fit of this group (" << NumberOfPeaksBeforeFit
806  << " peaks before refit, " << nExponentialsForFit
807  << " peaks after refit): " << std::endl;
808  if (nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax)
809  std::cout << "chi2/ndf of this fit (" << chi2PerNDF
810  << ") is higher than threshold (" << fChi2NDFMax << ")."
811  << std::endl;
812  if (nExponentialsForFit >= 2 &&
813  chi2PerNDF > fChi2NDFMaxFactorMultiHits * fChi2NDFMax)
814  std::cout << "chi2/ndf of this fit (" << chi2PerNDF
815  << ") is higher than threshold ("
816  << fChi2NDFMaxFactorMultiHits * fChi2NDFMax << ")." << std::endl;
817  }
818  std::cout << "---> DO NOT create hit object from fit parameters but use peak "
819  "values instead."
820  << std::endl;
821  std::cout << "---> Set fit parameter so that a sharp peak with a width of 1 tick "
822  "is shown in the event display. This indicates that the fit failed."
823  << std::endl;
824  }
825  peakWidth =
826  (((double)endTthisHit - (double)startTthisHit) / 4.) /
827  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
828  peakMeanErr = peakWidth / 2;
829  charge = sumADC;
830  peakMeanTrue = std::get<0>(peakVals.at(i));
831 
832  //set the fit values to make it visible in the event display that this fit failed
833  peakMean = peakMeanTrue;
834  peakTau1 = 0.008;
835  peakTau2 = 0.0065;
836  peakAmp = 20.;
837  }
838 
839  // Create the hit
840  recob::HitCreator hitcreator(
841  *wire, // wire reference
842  wid, // wire ID
843  startTthisHit + roiFirstBinTick, // start_tick TODO check
844  endTthisHit + roiFirstBinTick, // end_tick TODO check
845  peakWidth, // rms
846  peakMeanTrue + roiFirstBinTick, // peak_time
847  peakMeanErr, // sigma_peak_time
848  peakAmpTrue, // peak_amplitude
849  peakAmpErr, // sigma_peak_amplitude
850  charge, // hit_integral
851  chargeErr, // hit_sigma_integral
852  sumADC, // summedADC FIXME
853  nExponentialsForFit, // multiplicity
854  numHits, // local_index TODO check that the order is correct
855  chi2PerNDF, // goodness_of_fit
856  NDF // dof
857  );
858 
859  if (fLogLevel >= 6) {
860  std::cout << std::endl;
861  std::cout << "- Created hit object for peak #" << i
862  << " in this group with the following parameters (obtained from fit):"
863  << std::endl;
864  std::cout << "HitStartTick: " << startTthisHit + roiFirstBinTick << std::endl;
865  std::cout << "HitEndTick: " << endTthisHit + roiFirstBinTick << std::endl;
866  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
867  std::cout << "HitMeanTick: " << peakMeanTrue + roiFirstBinTick << " +- "
868  << peakMeanErr << std::endl;
869  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr
870  << std::endl;
871  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr
872  << std::endl;
873  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
874  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
875  std::cout << "HitIndex in group: " << numHits << std::endl;
876  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
877  std::cout << "HitNDF: " << NDF << std::endl;
878  }
879 
880  const recob::Hit hit(hitcreator.move());
881 
882  hcol.emplace_back(std::move(hit), wire, rawdigits);
883  // add fit parameters associated to the hit just pushed to the collection
884  std::array<float, 4> fitParams;
885  fitParams[0] = peakMean + roiFirstBinTick;
886  fitParams[1] = peakTau1;
887  fitParams[2] = peakTau2;
888  fitParams[3] = peakAmp;
889  fHitParamWriter.addVector(hitID, fitParams);
890  numHits++;
891  } // <---End loop over Exponentials
892  // } // <---End if chi2 <= chi2Max
893  } // <---End if(NumberOfPeaksBeforeFit <= fMaxMultiHit && width <= fMaxGroupLength), then fit
894 
895  // #######################################################
896  // ### If too large then force alternate solution ###
897  // ### - Make n hits from pulse train where n will ###
898  // ### depend on the fhicl parameter fLongPulseWidth ###
899  // ### Also do this if chi^2 is too large ###
900  // #######################################################
901  if (NumberOfPeaksBeforeFit > fMaxMultiHit || (width > fMaxGroupLength) ||
902  NFluctuations > fMaxFluctuations) {
903 
904  int nHitsInThisGroup = (endT - startT + 1) / fLongPulseWidth;
905 
906  if (nHitsInThisGroup > fLongMaxHits) {
907  nHitsInThisGroup = fLongMaxHits;
908  fLongPulseWidth = (endT - startT + 1) / nHitsInThisGroup;
909  }
910 
911  if (nHitsInThisGroup * fLongPulseWidth < (endT - startT + 1)) nHitsInThisGroup++;
912 
913  int firstTick = startT;
914  int lastTick = std::min(endT, firstTick + fLongPulseWidth - 1);
915 
916  if (fLogLevel >= 1) {
917  if (NumberOfPeaksBeforeFit > fMaxMultiHit) {
918  std::cout << std::endl;
919  std::cout << "WARNING: Number of peaks in this group (" << NumberOfPeaksBeforeFit
920  << ") is higher than threshold (" << fMaxMultiHit << ")." << 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 (width > fMaxGroupLength) {
926  std::cout << std::endl;
927  std::cout << "WARNING: group of peak is longer (" << width
928  << " ticks) than threshold (" << fMaxGroupLength << " ticks)."
929  << std::endl;
930  std::cout
931  << "---> DO NOT fit. Split group of peaks into hits with equal length instead."
932  << std::endl;
933  }
934  if (NFluctuations > fMaxFluctuations) {
935  std::cout << std::endl;
936  std::cout << "WARNING: fluctuations (" << NFluctuations
937  << ") higher than threshold (" << fMaxFluctuations << ")." << std::endl;
938  std::cout
939  << "---> DO NOT fit. Split group of peaks into hits with equal length instead."
940  << std::endl;
941  }
942  /*
943  if( ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) || ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) )
944  {
945  std::cout << std::endl;
946  std::cout << "WARNING: For fit of this group (" << NumberOfPeaksBeforeFit << " peaks before refit, " << nExponentialsForFit << " peaks after refit): " << std::endl;
947  if ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMax << ")." << std::endl;
948  if ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMaxFactorMultiHits*fChi2NDFMax << ")." << std::endl;
949  std::cout << "---> DO NOT create hit object but split group of peaks into hits with equal length instead." << std::endl;
950  }*/
951  std::cout << "---> Group goes from tick " << roiFirstBinTick + startT << " to "
952  << roiFirstBinTick + endT << ". Split group into ("
953  << roiFirstBinTick + endT << " - " << roiFirstBinTick + startT << ")/"
954  << fLongPulseWidth << " = " << (endT - startT) << "/" << fLongPulseWidth
955  << " = " << nHitsInThisGroup << " peaks (" << fLongPulseWidth
956  << " = LongPulseWidth), or maximum LongMaxHits = " << fLongMaxHits
957  << " peaks." << std::endl;
958  }
959 
960  for (int hitIdx = 0; hitIdx < nHitsInThisGroup; hitIdx++) {
961  // This hit parameters
962  double peakWidth =
963  ((lastTick - firstTick) / 4.) /
964  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
965  double peakMeanTrue = (firstTick + lastTick) / 2.;
966  if (NumberOfPeaksBeforeFit == 1 && nHitsInThisGroup == 1)
967  peakMeanTrue = std::get<0>(peakVals.at(
968  0)); //if only one peak was found, we want the mean of this peak to be the tick with the max. ADC count
969  double peakMeanErr = (lastTick - firstTick) / 2.;
970  double sumADC =
971  std::accumulate(signal.begin() + firstTick, signal.begin() + lastTick + 1, 0.);
972  double charge = sumADC;
973  double chargeErr = 0.1 * sumADC;
974  double peakAmpTrue = 0;
975 
976  for (int tick = firstTick; tick <= lastTick; tick++) {
977  if (signal[tick] > peakAmpTrue) peakAmpTrue = signal[tick];
978  }
979 
980  double peakAmpErr = 1.;
981  nExponentialsForFit = nHitsInThisGroup;
982  NDF = -1;
983  chi2PerNDF = -1.;
984  //set the fit values to make it visible in the event display that this fit failed
985  double peakMean = peakMeanTrue - 2;
986  double peakTau1 = 0.008;
987  double peakTau2 = 0.0065;
988  double peakAmp = 20.;
989 
990  recob::HitCreator hitcreator(
991  *wire, // wire reference
992  wid, // wire ID
993  firstTick + roiFirstBinTick, // start_tick TODO check
994  lastTick + roiFirstBinTick, // end_tick TODO check
995  peakWidth, // rms
996  peakMeanTrue + roiFirstBinTick, // peak_time
997  peakMeanErr, // sigma_peak_time
998  peakAmpTrue, // peak_amplitude
999  peakAmpErr, // sigma_peak_amplitude
1000  charge, // hit_integral
1001  chargeErr, // hit_sigma_integral
1002  sumADC, // summedADC FIXME
1003  nExponentialsForFit, // multiplicity
1004  hitIdx, // local_index TODO check that the order is correct
1005  chi2PerNDF, // goodness_of_fit
1006  NDF // dof
1007  );
1008 
1009  if (fLogLevel >= 6) {
1010  std::cout << std::endl;
1011  std::cout
1012  << "- Created hit object for peak #" << hitIdx
1013  << " in this group with the following parameters (obtained from waveform):"
1014  << std::endl;
1015  std::cout << "HitStartTick: " << firstTick + roiFirstBinTick << std::endl;
1016  std::cout << "HitEndTick: " << lastTick + roiFirstBinTick << std::endl;
1017  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
1018  std::cout << "HitMeanTick: " << peakMeanTrue + roiFirstBinTick << " +- "
1019  << peakMeanErr << std::endl;
1020  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr
1021  << std::endl;
1022  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr
1023  << std::endl;
1024  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
1025  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
1026  std::cout << "HitIndex in group: " << hitIdx << std::endl;
1027  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
1028  std::cout << "HitNDF: " << NDF << std::endl;
1029  }
1030  const recob::Hit hit(hitcreator.move());
1031  hcol.emplace_back(std::move(hit), wire, rawdigits);
1032 
1033  std::array<float, 4> fitParams;
1034  fitParams[0] = peakMean + roiFirstBinTick;
1035  fitParams[1] = peakTau1;
1036  fitParams[2] = peakTau2;
1037  fitParams[3] = peakAmp;
1038  fHitParamWriter.addVector(hitID, fitParams);
1039 
1040  // set for next loop
1041  firstTick = lastTick + 1;
1042  lastTick = std::min(firstTick + fLongPulseWidth - 1, endT);
1043 
1044  } //<---Hits in this group
1045  } //<---End if #peaks > MaxMultiHit
1046  fChi2->Fill(chi2PerNDF);
1047  } //<---End loop over merged candidate hits
1048  } //<---End looping over ROI's
1049  } //<---End looping over all the wires
1050 
1051  //==================================================================================================
1052  // End of the event
1053 
1054  // move the hit collection and the associations into the event
1055  hcol.put_into(evt);
1056 
1057  // and put hit fit parameters together with metadata into the event
1058  fHitParamWriter.saveOutputs(evt);
1059 
1060  } // End of produce()
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:1065
geo::GeometryCore::ChannelToWire
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1081
hit::DPRawHitFinder::mergeCandidatePeaks
void mergeCandidatePeaks(const std::vector< float > signalVec, TimeValsVec, MergedTimeWidVec &)
Definition: DPRawHitFinder_module.cc:1141
raw::RawDigit::Channel
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:213
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:1514
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
std
STL namespace.
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
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
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:1690
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
recob::HitCreator
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
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:1641
recob::HitCollectionCreator
A class handling a collection of hits and its associations.
Definition: HitCreator.h:481
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
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
hit::DPRawHitFinder::ChargeFunc
double ChargeFunc(double fPeakMean, double fPeakAmp, double fPeakTau1, double fPeakTau2, double fChargeNormFactor, double fPeakMeanTrue)
Definition: DPRawHitFinder_module.cc:1809
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
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
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:1295
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
util::begin
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
hit::DPRawHitFinder::doBinAverage
void doBinAverage(const std::vector< float > &inputVec, std::vector< float > &outputVec, size_t binsToAverage) const
Definition: DPRawHitFinder_module.cc:1846
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
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:399
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
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:1733
void hit::DPRawHitFinder::reBin ( const std::vector< float > &  inputVec,
std::vector< float > &  outputVec,
size_t  nBinsToCombine 
) const
private

Definition at line 1879 of file DPRawHitFinder_module.cc.

References DEFINE_ART_MODULE.

1882  {
1883  size_t nNewBins = inputVec.size() / nBinsToCombine;
1884 
1885  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1886 
1887  outputVec.resize(nNewBins, 0.);
1888 
1889  size_t outputBin = 0;
1890 
1891  for (size_t inputIdx = 0; inputIdx < inputVec.size();) {
1892  outputVec[outputBin] += inputVec[inputIdx++];
1893 
1894  if (inputIdx % nBinsToCombine == 0) outputBin++;
1895 
1896  if (outputBin > outputVec.size()) {
1897  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin
1898  << ", inputIdx = " << inputIdx << std::endl;
1899  break;
1900  }
1901  }
1902 
1903  return;
1904  }
void art::Modifier::registerProducts ( ProductDescriptions &  productsToRegister)
inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

17  {
18  ProductRegistryHelper::registerProducts(productsToRegister,
19  moduleDescription());
20  }
art::ProductRegistryHelper::registerProducts
void registerProducts(ProductDescriptions &productsToRegister, ModuleDescription const &md)
Definition: ProductRegistryHelper.cc:24
art::ModuleBase::moduleDescription
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  md)
inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

32  {
33  md_ = md;
34  }
art::ModuleBase::md_
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

50  {
51  // Now that we know we have seen all the consumes declarations,
52  // sort the results for fast lookup later.
53  collector_.sortConsumables(current_process_name);
54  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::sortConsumables
void sortConsumables(std::string const &current_process_name)
Definition: ConsumesCollector.cc:19
void hit::DPRawHitFinder::SplitPeak ( std::tuple< double, int, int, int >  fPeakDevCand,
PeakTimeWidVec &  fpeakValsTemp 
)
private

Definition at line 1690 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1692  {
1693  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1694  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1695  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1696  int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1697 
1698  if (WidthOldPeakOld < 3) { return; } //can't split peaks with widths < 3
1699 
1700  int NewPeakMax = 0;
1701  int NewPeakStart = 0;
1702  int NewPeakEnd = 0;
1703  int OldPeakNewMax = 0;
1704  int OldPeakNewStart = 0;
1705  int OldPeakNewEnd = 0;
1706 
1707  OldPeakNewStart = OldPeakOldStart;
1708  NewPeakEnd = OldPeakOldEnd;
1709 
1710  OldPeakNewEnd = OldPeakNewStart + 0.5 * (WidthOldPeakOld + (WidthOldPeakOld % 2));
1711  NewPeakStart = OldPeakNewEnd + 1;
1712 
1713  int WidthOldPeakNew = OldPeakNewEnd - OldPeakNewStart;
1714  int WidthNewPeak = NewPeakEnd - NewPeakStart;
1715 
1716  OldPeakNewMax = OldPeakNewStart + 0.5 * (WidthOldPeakNew - (WidthOldPeakNew % 2));
1717  NewPeakMax = NewPeakStart + 0.5 * (WidthNewPeak - (WidthNewPeak % 2));
1718 
1719  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1720  std::make_tuple(OldPeakNewMax, 0, OldPeakNewStart, OldPeakNewEnd);
1721  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1722  std::sort(
1723  fpeakValsTemp.begin(),
1724  fpeakValsTemp.end(),
1725  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1726  return std::get<0>(t1) < std::get<0>(t2);
1727  });
1728 
1729  return;
1730  }
t1
TTree * t1
Definition: plottest35.C:26
t2
TTree * t2
Definition: plottest35.C:36
double hit::DPRawHitFinder::WidthFunc ( double  fPeakMean,
double  fPeakAmp,
double  fPeakTau1,
double  fPeakTau2,
double  fStartTime,
double  fEndTime,
double  fPeakMeanTrue 
)
private

Definition at line 1733 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1740  {
1741  double MaxValue = (fPeakAmp * exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau1)) /
1742  (1 + exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau2));
1743  double FuncValue = 0.;
1744  double HalfMaxLeftTime = 0.;
1745  double HalfMaxRightTime = 0.;
1746  double ReBin = 10.;
1747 
1748  //First iteration (+- 1 bin)
1749  for (double x = fPeakMeanTrue; x > fStartTime - 1000.; x--) {
1750  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1751  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1752  if (FuncValue < 0.5 * MaxValue) {
1753  HalfMaxLeftTime = x;
1754  break;
1755  }
1756  }
1757 
1758  for (double x = fPeakMeanTrue; x < fEndTime + 1000.; x++) {
1759  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1760  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1761  if (FuncValue < 0.5 * MaxValue) {
1762  HalfMaxRightTime = x;
1763  break;
1764  }
1765  }
1766 
1767  //Second iteration (+- 0.1 bin)
1768  for (double x = HalfMaxLeftTime + 1; x > HalfMaxLeftTime; x -= (1 / ReBin)) {
1769  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1770  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1771  if (FuncValue < 0.5 * MaxValue) {
1772  HalfMaxLeftTime = x;
1773  break;
1774  }
1775  }
1776 
1777  for (double x = HalfMaxRightTime - 1; x < HalfMaxRightTime; x += (1 / ReBin)) {
1778  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1779  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1780  if (FuncValue < 0.5 * MaxValue) {
1781  HalfMaxRightTime = x;
1782  break;
1783  }
1784  }
1785 
1786  //Third iteration (+- 0.01 bin)
1787  for (double x = HalfMaxLeftTime + 1 / ReBin; x > HalfMaxLeftTime; x -= 1 / (ReBin * ReBin)) {
1788  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1789  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1790  if (FuncValue < 0.5 * MaxValue) {
1791  HalfMaxLeftTime = x;
1792  break;
1793  }
1794  }
1795 
1796  for (double x = HalfMaxRightTime - 1 / ReBin; x < HalfMaxRightTime; x += 1 / (ReBin * ReBin)) {
1797  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1798  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1799  if (FuncValue < 0.5 * MaxValue) {
1800  HalfMaxRightTime = x;
1801  break;
1802  }
1803  }
1804 
1805  return HalfMaxRightTime - HalfMaxLeftTime;
1806  }
x
Float_t x
Definition: compare.C:6

Member Data Documentation

std::string hit::DPRawHitFinder::fCalDataModuleLabel
private

Definition at line 161 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fChargeNorm
private

Definition at line 172 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

TH1F* hit::DPRawHitFinder::fChi2
private

Definition at line 201 of file DPRawHitFinder_module.cc.

Referenced by beginJob(), and produce().

double hit::DPRawHitFinder::fChi2NDFMax
private

Definition at line 176 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fChi2NDFMaxFactorMultiHits
private

Definition at line 177 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fChi2NDFRetry
private

Definition at line 174 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fChi2NDFRetryFactorMultiHits
private

Definition at line 175 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

bool hit::DPRawHitFinder::fDoMergePeaks
private

Definition at line 185 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

TH1F* hit::DPRawHitFinder::fFirstChi2
private

Definition at line 200 of file DPRawHitFinder_module.cc.

Referenced by beginJob(), and produce().

double hit::DPRawHitFinder::fFitPeakMeanRange
private

Definition at line 181 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and FitExponentials().

int hit::DPRawHitFinder::fGroupMaxDistance
private

Definition at line 182 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

double hit::DPRawHitFinder::fGroupMinADC
private

Definition at line 183 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

anab::FVectorWriter<4> hit::DPRawHitFinder::fHitParamWriter
private

Definition at line 198 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fLogLevel
private

Definition at line 164 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), FitExponentials(), and produce().

int hit::DPRawHitFinder::fLongMaxHits
private

Definition at line 192 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fLongPulseWidth
private

Definition at line 193 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fMaxFluctuations
private

Definition at line 194 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fMaxGroupLength
private

Definition at line 171 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fMaxMultiHit
private

Definition at line 170 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fMaxTau
private

Definition at line 180 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and FitExponentials().

double hit::DPRawHitFinder::fMergeADCSumThreshold
private

Definition at line 186 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

double hit::DPRawHitFinder::fMergeMaxADCLimit
private

Definition at line 190 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

double hit::DPRawHitFinder::fMergeMaxADCThreshold
private

Definition at line 187 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

double hit::DPRawHitFinder::fMinADCSum
private

Definition at line 168 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fMinADCSumOverWidth
private

Definition at line 169 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fMinRelativePeakHeightLeft
private

Definition at line 188 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

double hit::DPRawHitFinder::fMinRelativePeakHeightRight
private

Definition at line 189 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and mergeCandidatePeaks().

float hit::DPRawHitFinder::fMinSig
private

Definition at line 165 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fMinTau
private

Definition at line 179 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and FitExponentials().

int hit::DPRawHitFinder::fMinWidth
private

Definition at line 167 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

art::InputTag hit::DPRawHitFinder::fNewHitsTag
private

Definition at line 197 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

size_t hit::DPRawHitFinder::fNumBinsToAverage
private

Definition at line 178 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

bool hit::DPRawHitFinder::fSameShape
private

Definition at line 184 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fTicksToStopPeakFinder
private

Definition at line 166 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and findCandidatePeaks().

bool hit::DPRawHitFinder::fTryNplus1Fits
private

Definition at line 173 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

double hit::DPRawHitFinder::fWidthNormalization
private

Definition at line 191 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

The documentation for this class was generated from the following file: