LArSoft  v10_04_05
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:261
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 1619 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1621  {
1622  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1623  int NewPeakMax = std::get<2>(fPeakDevCand);
1624  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1625  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1626  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1627 
1628  int NewPeakStart = 0;
1629  int NewPeakEnd = 0;
1630  int OldPeakNewStart = 0;
1631  int OldPeakNewEnd = 0;
1632  int DistanceBwOldAndNewPeak = 0;
1633 
1634  if (NewPeakMax < OldPeakMax) {
1635  NewPeakStart = OldPeakOldStart;
1636  OldPeakNewEnd = OldPeakOldEnd;
1637  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1638  NewPeakEnd = NewPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1639  if (DistanceBwOldAndNewPeak % 2 == 0) NewPeakEnd -= 1;
1640  OldPeakNewStart = NewPeakEnd + 1;
1641  }
1642  else if (OldPeakMax < NewPeakMax) {
1643  NewPeakEnd = OldPeakOldEnd;
1644  OldPeakNewStart = OldPeakOldStart;
1645  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1646  OldPeakNewEnd = OldPeakMax + 0.5 * (DistanceBwOldAndNewPeak - (DistanceBwOldAndNewPeak % 2));
1647  if (DistanceBwOldAndNewPeak % 2 == 0) OldPeakNewEnd -= 1;
1648  NewPeakStart = OldPeakNewEnd + 1;
1649  }
1650  else if (OldPeakMax == NewPeakMax) {
1651  return;
1652  } //This shouldn't happen
1653 
1654  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1655  std::make_tuple(OldPeakMax, 0, OldPeakNewStart, OldPeakNewEnd);
1656  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1657  std::sort(
1658  fpeakValsTemp.begin(),
1659  fpeakValsTemp.end(),
1660  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1661  return std::get<0>(t1) < std::get<0>(t2);
1662  });
1663 
1664  return;
1665  }
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 280 of file DPRawHitFinder_module.cc.

References fChi2, and fFirstChi2.

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

Definition at line 1787 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1794  {
1795  double ChargeSum = 0.;
1796  double Charge = 0.;
1797  double ReBin = 10.;
1798 
1799  bool ChargeBigEnough = true;
1800  for (double x = fPeakMeanTrue - 1 / ReBin; ChargeBigEnough && x > fPeakMeanTrue - 1000.;
1801  x -= 1.) {
1802  for (double i = 0.; i > -1.; i -= (1 / ReBin)) {
1803  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1804  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1805  ChargeSum += Charge;
1806  }
1807  if (Charge < 0.01) ChargeBigEnough = false;
1808  }
1809 
1810  ChargeBigEnough = true;
1811  for (double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue + 1000.; x += 1.) {
1812  for (double i = 0.; i < 1.; i += (1 / ReBin)) {
1813  Charge = (fPeakAmp * exp(0.4 * (x + i - fPeakMean) / fPeakTau1)) /
1814  (1 + exp(0.4 * (x + i - fPeakMean) / fPeakTau2));
1815  ChargeSum += Charge;
1816  }
1817  if (Charge < 0.01) ChargeBigEnough = false;
1818  }
1819 
1820  return ChargeSum * fChargeNormFactor / ReBin;
1821  }
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 1560 of file DPRawHitFinder_module.cc.

Referenced by FindPeakWithMaxDeviation(), and FitExponentials().

1561  {
1562  std::string feqn = ""; // string holding fit formula
1563  std::stringstream numConv;
1564 
1565  if (fSameShape) {
1566  for (int i = 0; i < fNPeaks; i++) {
1567  feqn.append("+( [");
1568  numConv.str("");
1569  numConv << 2 * (i + 1);
1570  feqn.append(numConv.str());
1571  feqn.append("] * exp(0.4*(x-[");
1572  numConv.str("");
1573  numConv << 2 * (i + 1) + 1;
1574  feqn.append(numConv.str());
1575  feqn.append("])/[");
1576  numConv.str("");
1577  numConv << 0;
1578  feqn.append(numConv.str());
1579  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1580  numConv.str("");
1581  numConv << 2 * (i + 1) + 1;
1582  feqn.append(numConv.str());
1583  feqn.append("])/[");
1584  numConv.str("");
1585  numConv << 1;
1586  feqn.append(numConv.str());
1587  feqn.append("]) ) )");
1588  }
1589  }
1590  else {
1591  for (int i = 0; i < fNPeaks; i++) {
1592  feqn.append("+( [");
1593  numConv.str("");
1594  numConv << 4 * i + 2;
1595  feqn.append(numConv.str());
1596  feqn.append("] * exp(0.4*(x-[");
1597  numConv.str("");
1598  numConv << 4 * i + 3;
1599  feqn.append(numConv.str());
1600  feqn.append("])/[");
1601  numConv.str("");
1602  numConv << 4 * i;
1603  feqn.append(numConv.str());
1604  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1605  numConv.str("");
1606  numConv << 2 * (i + 1) + 1;
1607  feqn.append(numConv.str());
1608  feqn.append("])/[");
1609  numConv.str("");
1610  numConv << 4 * i + 1;
1611  feqn.append(numConv.str());
1612  feqn.append("]) ) )");
1613  }
1614  }
1615  return feqn;
1616  }
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 1824 of file DPRawHitFinder_module.cc.

Referenced by produce().

1827  {
1828  size_t halfBinsToAverage(binsToAverage / 2);
1829 
1830  float runningSum(0.);
1831 
1832  for (size_t idx = 0; idx < halfBinsToAverage; idx++)
1833  runningSum += inputVec[idx];
1834 
1835  outputVec.resize(inputVec.size());
1836  std::vector<float>::iterator outputVecItr = outputVec.begin();
1837 
1838  // First pass through to build the erosion vector
1839  for (std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end();
1840  inputItr++) {
1841  size_t startOffset = std::distance(inputVec.begin(), inputItr);
1842  size_t stopOffset = std::distance(inputItr, inputVec.end());
1843  size_t count =
1844  std::min(2 * halfBinsToAverage,
1845  std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1846 
1847  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1848  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1849 
1850  *outputVecItr++ = runningSum / float(count);
1851  }
1852 
1853  return;
1854  }
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 1249 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks().

1253  {
1254  int NFluctuations = 0;
1255 
1256  for (int j = peakMean - 1; j >= peakStart; j--) {
1257  if (fsignalVec[j] < 5)
1258  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1259 
1260  if (fsignalVec[j] > fsignalVec[j + 1]) { NFluctuations++; }
1261  }
1262 
1263  for (int j = peakMean + 1; j <= peakEnd; j++) {
1264  if (fsignalVec[j] < 5)
1265  break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1266 
1267  if (fsignalVec[j] > fsignalVec[j - 1]) { NFluctuations++; }
1268  }
1269 
1270  return NFluctuations;
1271  }
void hit::DPRawHitFinder::FillOutHitParameterVector ( const std::vector< double > &  input,
std::vector< double > &  output 
)
private

Definition at line 260 of file DPRawHitFinder_module.cc.

References Get.

262  {
263  if (input.size() == 0)
264  throw std::runtime_error(
265  "DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
266 
267  const unsigned int N_PLANES = art::ServiceHandle<geo::WireReadout>()->Get().Nplanes();
268 
269  if (input.size() == 1)
270  output.resize(N_PLANES, input[0]);
271  else if (input.size() == N_PLANES)
272  output = input;
273  else
274  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config "
275  "vector size !=1 and !=N_PLANES.");
276  }
Get
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
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 1047 of file DPRawHitFinder_module.cc.

References fTicksToStopPeakFinder.

Referenced by produce().

1052  {
1053  // Need a minimum number of ticks to do any work here
1054  if (std::distance(startItr, stopItr) > 4) {
1055  // Find the highest peak in the range given
1056  auto maxItr = std::max_element(startItr, stopItr);
1057 
1058  float maxValue = *maxItr;
1059  int maxTime = std::distance(startItr, maxItr);
1060 
1061  if (maxValue >= PeakMin) {
1062  // backwards to find first bin for this candidate hit
1063  auto firstItr = std::distance(startItr, maxItr) > 2 ? maxItr - 1 : startItr;
1064  bool PeakStartIsHere = true;
1065 
1066  while (firstItr != startItr) {
1067  //Check for inflection point & ADC count <= 0
1068  PeakStartIsHere = true;
1069  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1070  if (*firstItr >= *(firstItr - i)) {
1071  PeakStartIsHere = false;
1072  break;
1073  }
1074  }
1075  if (*firstItr <= 0 || PeakStartIsHere) break;
1076  firstItr--;
1077  }
1078 
1079  int firstTime = std::distance(startItr, firstItr);
1080 
1081  // Recursive call to find all candidate hits earlier than this peak
1082  findCandidatePeaks(startItr, firstItr - 1, timeValsVec, PeakMin, firstTick);
1083 
1084  // forwards to find last bin for this candidate hit
1085  auto lastItr = std::distance(maxItr, stopItr) > 2 ? maxItr + 1 : stopItr - 1;
1086  bool PeakEndIsHere = true;
1087 
1088  while (lastItr != stopItr) {
1089  //Check for inflection point & ADC count <= 0
1090  PeakEndIsHere = true;
1091  for (int i = 1; i <= fTicksToStopPeakFinder; i++) {
1092  if (*lastItr >= *(lastItr + i)) {
1093  PeakEndIsHere = false;
1094  break;
1095  }
1096  }
1097  if (*lastItr <= 0 || PeakEndIsHere) break;
1098  lastItr++;
1099  }
1100 
1101  int lastTime = std::distance(startItr, lastItr);
1102 
1103  // Now save this candidate's start and max time info
1104  timeValsVec.push_back(
1105  std::make_tuple(firstTick + firstTime, firstTick + maxTime, firstTick + lastTime));
1106 
1107  // Recursive call to find all candidate hits later than this peak
1108  findCandidatePeaks(lastItr + 1,
1109  stopItr,
1110  timeValsVec,
1111  PeakMin,
1112  firstTick + std::distance(startItr, lastItr + 1));
1113  }
1114  }
1115 
1116  return;
1117  }
hit::DPRawHitFinder::findCandidatePeaks
void findCandidatePeaks(std::vector< float >::const_iterator startItr, std::vector< float >::const_iterator stopItr, TimeValsVec &timeValsVec, float &PeakMin, int firstTick) const
Definition: DPRawHitFinder_module.cc:1047
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 1495 of file DPRawHitFinder_module.cc.

References CreateFitFunction(), t1, and t2.

Referenced by produce().

1503  {
1504  std::string eqn =
1505  CreateFitFunction(fNPeaks, fSameShape); // string for equation of Exponentials fit
1506 
1507  TF1 Exponentials("Exponentials", eqn.c_str(), fStartTime, fEndTime + 1);
1508 
1509  for (size_t i = 0; i < fparamVec.size(); i++) {
1510  Exponentials.SetParameter(i, fparamVec[i].first);
1511  }
1512 
1513  // ##########################################################################
1514  // ### Finding the peak with the max chi2 fit and signal ###
1515  // ##########################################################################
1516  double Chi2PerNDFPeak;
1517  double MaxPosDeviation;
1518  double MaxNegDeviation;
1519  int BinMaxPosDeviation;
1520  int BinMaxNegDeviation;
1521  for (int i = 0; i < fNPeaks; i++) {
1522  Chi2PerNDFPeak = 0.;
1523  MaxPosDeviation = 0.;
1524  MaxNegDeviation = 0.;
1525  BinMaxPosDeviation = 0;
1526  BinMaxNegDeviation = 0;
1527 
1528  for (int j = std::get<2>(fpeakVals.at(i)); j < std::get<3>(fpeakVals.at(i)) + 1; j++) {
1529  if ((Exponentials(j + 0.5) - fSignalVector[j]) > MaxPosDeviation &&
1530  j != std::get<0>(fpeakVals.at(i))) {
1531  MaxPosDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1532  BinMaxPosDeviation = j;
1533  }
1534  if ((Exponentials(j + 0.5) - fSignalVector[j]) < MaxNegDeviation &&
1535  j != std::get<0>(fpeakVals.at(i))) {
1536  MaxNegDeviation = Exponentials(j + 0.5) - fSignalVector[j];
1537  BinMaxNegDeviation = j;
1538  }
1539  Chi2PerNDFPeak +=
1540  pow((Exponentials(j + 0.5) - fSignalVector[j]) / sqrt(fSignalVector[j]), 2);
1541  }
1542 
1543  if (BinMaxNegDeviation != 0) {
1544  Chi2PerNDFPeak /=
1545  static_cast<double>((std::get<3>(fpeakVals.at(i)) - std::get<2>(fpeakVals.at(i))));
1546  fPeakDev.emplace_back(Chi2PerNDFPeak, i, BinMaxNegDeviation, BinMaxPosDeviation);
1547  }
1548  }
1549 
1550  std::sort(
1551  fPeakDev.begin(),
1552  fPeakDev.end(),
1553  [](std::tuple<double, int, int, int> const& t1, std::tuple<double, int, int, int> const& t2) {
1554  return std::get<0>(t1) > std::get<0>(t2);
1555  });
1556  Exponentials.Delete();
1557  }
t1
TTree * t1
Definition: plottest35.C:26
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1560
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 1276 of file DPRawHitFinder_module.cc.

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

Referenced by produce().

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

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

Referenced by produce().

1126  {
1127  // ################################################################
1128  // ### Lets loop over the candidate pulses we found in this ROI ###
1129  // ################################################################
1130  auto timeValsVecItr = timeValsVec.begin();
1131  unsigned int PeaksInThisMergedPeak = 0;
1132  //int EndTickOfPreviousMergedPeak=0;
1133  while (timeValsVecItr != timeValsVec.end()) {
1134  PeakTimeWidVec peakVals;
1135 
1136  // Setting the start, peak, and end time of the pulse
1137  auto& timeVal = *timeValsVecItr++;
1138  int startT = std::get<0>(timeVal);
1139  int maxT = std::get<1>(timeVal);
1140  int endT = std::get<2>(timeVal);
1141  int widT = endT - startT;
1142  int FinalStartT = startT;
1143  int FinalEndT = endT;
1144 
1145  int NFluctuations = EstimateFluctuations(signalVec, startT, maxT, endT);
1146 
1147  peakVals.emplace_back(maxT, widT, startT, endT);
1148 
1149  // See if we want to merge pulses together
1150  // First check if we have more than one pulse on the wire
1151  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1152 
1153  // Loop until no more merged pulses (or candidates in this ROI)
1154  while (checkNextHit) {
1155  // group hits if start time of the next pulse is < end time + fGroupMaxDistance of current pulse and if intermediate signal between two pulses doesn't go below fMinBinToGroup and if the hits are not all above the merge max adc limit
1156  int NextStartT = std::get<0>(*timeValsVecItr);
1157 
1158  double MinADC = signalVec[endT];
1159  for (int i = endT; i <= NextStartT; i++) {
1160  if (signalVec[i] < MinADC) { MinADC = signalVec[i]; }
1161  }
1162 
1163  // Group peaks (grouped peaks are fitted together and can be merged)
1164  if (MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance) {
1165  int CurrentStartT = startT;
1166  int CurrentMaxT = maxT;
1167  int CurrentEndT = endT;
1168 
1169  timeVal = *timeValsVecItr++;
1170  int NextMaxT = std::get<1>(timeVal);
1171  int NextEndT = std::get<2>(timeVal);
1172  int NextWidT = NextEndT - NextStartT;
1173  FinalEndT = NextEndT;
1174 
1175  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1176 
1177  int CurrentSumADC = 0;
1178  for (int i = CurrentStartT; i <= CurrentEndT; i++) {
1179  CurrentSumADC += signalVec[i];
1180  }
1181 
1182  int NextSumADC = 0;
1183  for (int i = NextStartT; i <= NextEndT; i++) {
1184  NextSumADC += signalVec[i];
1185  }
1186 
1187  //Merge peaks within a group
1188  if (fDoMergePeaks) {
1189  if (signalVec[NextMaxT] <= signalVec[CurrentMaxT] &&
1190  ((signalVec[NextMaxT] < fMergeMaxADCThreshold * signalVec[CurrentMaxT] &&
1191  NextSumADC < fMergeADCSumThreshold * CurrentSumADC) ||
1192  (signalVec[NextMaxT] - signalVec[NextStartT]) <
1193  fMinRelativePeakHeightRight * signalVec[CurrentMaxT]) &&
1194  (signalVec[NextMaxT] <= fMergeMaxADCLimit)) {
1195  maxT = CurrentMaxT;
1196  startT = CurrentStartT;
1197  endT = NextEndT;
1198  widT = endT - startT;
1199  peakVals.pop_back();
1200  peakVals.emplace_back(maxT, widT, startT, endT);
1201  }
1202  else if (signalVec[NextMaxT] > signalVec[CurrentMaxT] &&
1203  ((signalVec[CurrentMaxT] < fMergeMaxADCThreshold * signalVec[NextMaxT] &&
1204  CurrentSumADC < fMergeADCSumThreshold * NextSumADC) ||
1205  (signalVec[CurrentMaxT] - signalVec[CurrentEndT]) <
1206  fMinRelativePeakHeightLeft * signalVec[NextMaxT]) &&
1207  (signalVec[CurrentMaxT] <= fMergeMaxADCLimit)) {
1208  maxT = NextMaxT;
1209  startT = CurrentStartT;
1210  endT = NextEndT;
1211  widT = endT - startT;
1212  peakVals.pop_back();
1213  peakVals.emplace_back(maxT, widT, startT, endT);
1214  }
1215  else {
1216  maxT = NextMaxT;
1217  startT = NextStartT;
1218  endT = NextEndT;
1219  widT = NextWidT;
1220  peakVals.emplace_back(maxT, widT, startT, endT);
1221  PeaksInThisMergedPeak++;
1222  }
1223  }
1224  else {
1225  maxT = NextMaxT;
1226  startT = NextStartT;
1227  endT = NextEndT;
1228  widT = NextWidT;
1229  peakVals.emplace_back(maxT, widT, startT, endT);
1230  PeaksInThisMergedPeak++;
1231  }
1232  checkNextHit = timeValsVecItr != timeValsVec.end();
1233  } //<---Checking adjacent pulses
1234  else {
1235  checkNextHit = false;
1236  PeaksInThisMergedPeak = 0;
1237  }
1238 
1239  } //<---End checking if there is more than one pulse on the wire
1240  // Add these to our merged vector
1241  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1242  }
1243  return;
1244  }
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:1249
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 292 of file DPRawHitFinder_module.cc.

References AddPeak(), anab::FVectorWriter< N >::addVector(), recob::Wire::Channel(), 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, Get, 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.

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

Definition at line 1857 of file DPRawHitFinder_module.cc.

References DEFINE_ART_MODULE.

1860  {
1861  size_t nNewBins = inputVec.size() / nBinsToCombine;
1862 
1863  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1864 
1865  outputVec.resize(nNewBins, 0.);
1866 
1867  size_t outputBin = 0;
1868 
1869  for (size_t inputIdx = 0; inputIdx < inputVec.size();) {
1870  outputVec[outputBin] += inputVec[inputIdx++];
1871 
1872  if (inputIdx % nBinsToCombine == 0) outputBin++;
1873 
1874  if (outputBin > outputVec.size()) {
1875  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin
1876  << ", inputIdx = " << inputIdx << std::endl;
1877  break;
1878  }
1879  }
1880 
1881  return;
1882  }
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 1668 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1670  {
1671  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1672  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1673  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1674  int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1675 
1676  if (WidthOldPeakOld < 3) { return; } //can't split peaks with widths < 3
1677 
1678  int NewPeakMax = 0;
1679  int NewPeakStart = 0;
1680  int NewPeakEnd = 0;
1681  int OldPeakNewMax = 0;
1682  int OldPeakNewStart = 0;
1683  int OldPeakNewEnd = 0;
1684 
1685  OldPeakNewStart = OldPeakOldStart;
1686  NewPeakEnd = OldPeakOldEnd;
1687 
1688  OldPeakNewEnd = OldPeakNewStart + 0.5 * (WidthOldPeakOld + (WidthOldPeakOld % 2));
1689  NewPeakStart = OldPeakNewEnd + 1;
1690 
1691  int WidthOldPeakNew = OldPeakNewEnd - OldPeakNewStart;
1692  int WidthNewPeak = NewPeakEnd - NewPeakStart;
1693 
1694  OldPeakNewMax = OldPeakNewStart + 0.5 * (WidthOldPeakNew - (WidthOldPeakNew % 2));
1695  NewPeakMax = NewPeakStart + 0.5 * (WidthNewPeak - (WidthNewPeak % 2));
1696 
1697  fpeakValsTemp.at(PeakNumberWithNewPeak) =
1698  std::make_tuple(OldPeakNewMax, 0, OldPeakNewStart, OldPeakNewEnd);
1699  fpeakValsTemp.emplace_back(NewPeakMax, 0, NewPeakStart, NewPeakEnd);
1700  std::sort(
1701  fpeakValsTemp.begin(),
1702  fpeakValsTemp.end(),
1703  [](std::tuple<int, int, int, int> const& t1, std::tuple<int, int, int, int> const& t2) {
1704  return std::get<0>(t1) < std::get<0>(t2);
1705  });
1706 
1707  return;
1708  }
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 1711 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1718  {
1719  double MaxValue = (fPeakAmp * exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau1)) /
1720  (1 + exp(0.4 * (fPeakMeanTrue - fPeakMean) / fPeakTau2));
1721  double FuncValue = 0.;
1722  double HalfMaxLeftTime = 0.;
1723  double HalfMaxRightTime = 0.;
1724  double ReBin = 10.;
1725 
1726  //First iteration (+- 1 bin)
1727  for (double x = fPeakMeanTrue; x > fStartTime - 1000.; x--) {
1728  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1729  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1730  if (FuncValue < 0.5 * MaxValue) {
1731  HalfMaxLeftTime = x;
1732  break;
1733  }
1734  }
1735 
1736  for (double x = fPeakMeanTrue; x < fEndTime + 1000.; x++) {
1737  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1738  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1739  if (FuncValue < 0.5 * MaxValue) {
1740  HalfMaxRightTime = x;
1741  break;
1742  }
1743  }
1744 
1745  //Second iteration (+- 0.1 bin)
1746  for (double x = HalfMaxLeftTime + 1; x > HalfMaxLeftTime; x -= (1 / ReBin)) {
1747  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1748  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1749  if (FuncValue < 0.5 * MaxValue) {
1750  HalfMaxLeftTime = x;
1751  break;
1752  }
1753  }
1754 
1755  for (double x = HalfMaxRightTime - 1; x < HalfMaxRightTime; x += (1 / ReBin)) {
1756  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1757  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1758  if (FuncValue < 0.5 * MaxValue) {
1759  HalfMaxRightTime = x;
1760  break;
1761  }
1762  }
1763 
1764  //Third iteration (+- 0.01 bin)
1765  for (double x = HalfMaxLeftTime + 1 / ReBin; x > HalfMaxLeftTime; x -= 1 / (ReBin * ReBin)) {
1766  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1767  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1768  if (FuncValue < 0.5 * MaxValue) {
1769  HalfMaxLeftTime = x;
1770  break;
1771  }
1772  }
1773 
1774  for (double x = HalfMaxRightTime - 1 / ReBin; x < HalfMaxRightTime; x += 1 / (ReBin * ReBin)) {
1775  FuncValue = (fPeakAmp * exp(0.4 * (x - fPeakMean) / fPeakTau1)) /
1776  (1 + exp(0.4 * (x - fPeakMean) / fPeakTau2));
1777  if (FuncValue < 0.5 * MaxValue) {
1778  HalfMaxRightTime = x;
1779  break;
1780  }
1781  }
1782 
1783  return HalfMaxRightTime - HalfMaxLeftTime;
1784  }
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: