LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Classes | Public Types | Public Member Functions | Static 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::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Classes

struct  Comp
 

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 DPRawHitFinder (fhicl::ParameterSet const &pset)
 
void produce (art::Event &evt) override
 
void beginJob () override
 
void endJob () override
 
void reconfigure (fhicl::ParameterSet const &p)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > consumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > consumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > mayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > mayConsumeView (InputTag const &it)
 
base_engine_t & createEngine (seed_t seed)
 
base_engine_t & createEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_t & createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumer > non_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

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 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 79 of file DPRawHitFinder_module.cc.

Member Typedef Documentation

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

Definition at line 94 of file DPRawHitFinder_module.cc.

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

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

Definition at line 96 of file DPRawHitFinder_module.cc.

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

Definition at line 95 of file DPRawHitFinder_module.cc.

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

Definition at line 93 of file DPRawHitFinder_module.cc.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

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

Definition at line 92 of file DPRawHitFinder_module.cc.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

Definition at line 226 of file DPRawHitFinder_module.cc.

References recob::HitAndAssociationsWriterBase::declare_products(), fHitParamWriter, anab::FVectorWriter< N >::produces_using(), and reconfigure().

226  :
227  fNewHitsTag(
228  pset.get<std::string>("module_label"), "",
229  art::ServiceHandle<art::TriggerNamesService>()->getProcessName()),
230  fHitParamWriter(this)
231 {
232  this->reconfigure(pset);
233 
234  // let HitCollectionCreator declare that we are going to produce
235  // hits and associations with wires and raw digits
236  // (with no particular product label)
237  recob::HitCollectionCreator::declare_products(*this);
238 
239  // declare that data products with feature vectors describing
240  // hits is going to be produced
241  fHitParamWriter.produces_using< recob::Hit >();
242 
243 } // DPRawHitFinder::DPRawHitFinder()
recob::HitAndAssociationsWriterBase::declare_products
static void declare_products(ModuleType &producer, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.h:1117
hit::DPRawHitFinder::fHitParamWriter
anab::FVectorWriter< 4 > fHitParamWriter
Definition: DPRawHitFinder_module.cc:213
hit::DPRawHitFinder::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: DPRawHitFinder_module.cc:269
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49

Member Function Documentation

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

Definition at line 1598 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1600 {
1601  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1602  int NewPeakMax = std::get<2>(fPeakDevCand);
1603  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1604  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1605  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1606 
1607  int NewPeakStart=0;
1608  int NewPeakEnd=0;
1609  int OldPeakNewStart=0;
1610  int OldPeakNewEnd=0;
1611  int DistanceBwOldAndNewPeak=0;
1612 
1613  if(NewPeakMax<OldPeakMax)
1614  {
1615  NewPeakStart = OldPeakOldStart;
1616  OldPeakNewEnd = OldPeakOldEnd;
1617  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1618  NewPeakEnd = NewPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1619  if(DistanceBwOldAndNewPeak%2 == 0) NewPeakEnd -= 1;
1620  OldPeakNewStart = NewPeakEnd+1;
1621  }
1622  else if(OldPeakMax<NewPeakMax)
1623  {
1624  NewPeakEnd = OldPeakOldEnd;
1625  OldPeakNewStart = OldPeakOldStart;
1626  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1627  OldPeakNewEnd = OldPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1628  if(DistanceBwOldAndNewPeak%2 == 0) OldPeakNewEnd -= 1;
1629  NewPeakStart = OldPeakNewEnd+1;
1630  }
1631  else if(OldPeakMax==NewPeakMax){return;} //This shouldn't happen
1632 
1633  fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakMax,0,OldPeakNewStart,OldPeakNewEnd);
1634  fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1635  std::sort(fpeakValsTemp.begin(),fpeakValsTemp.end(), [](std::tuple<int,int,int,int> const &t1, std::tuple<int,int,int,int> const &t2) {return std::get<0>(t1) < std::get<0>(t2);} );
1636 
1637  return;
1638 }
t1
TTree * t1
Definition: plottest35.C:26
t2
TTree * t2
Definition: plottest35.C:36
void hit::DPRawHitFinder::beginJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 308 of file DPRawHitFinder_module.cc.

References fChi2, fFirstChi2, and art::TFileDirectory::make().

309 {
310  // get access to the TFile service
311  art::ServiceHandle<art::TFileService> tfs;
312 
313  // ======================================
314  // === Hit Information for Histograms ===
315  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 5000);
316  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 5000);
317 }
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
double hit::DPRawHitFinder::ChargeFunc ( double  fPeakMean,
double  fPeakAmp,
double  fPeakTau1,
double  fPeakTau2,
double  fChargeNormFactor,
double  fPeakMeanTrue 
)
private

Definition at line 1761 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1768 {
1769 double ChargeSum = 0.;
1770 double Charge = 0.;
1771 double ReBin = 10.;
1772 
1773 bool ChargeBigEnough=true;
1774  for(double x = fPeakMeanTrue - 1/ReBin; ChargeBigEnough && x > fPeakMeanTrue-1000.; x-=1.)
1775  {
1776  for(double i=0.; i > -1.; i-=(1/ReBin))
1777  {
1778  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1779  ChargeSum += Charge;
1780  }
1781  if(Charge < 0.01) ChargeBigEnough = false;
1782  }
1783 
1784 ChargeBigEnough=true;
1785  for(double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue+1000.; x+=1.)
1786  {
1787  for(double i=0.; i < 1.; i+=(1/ReBin))
1788  {
1789  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1790  ChargeSum += Charge;
1791  }
1792  if(Charge < 0.01) ChargeBigEnough = false;
1793  }
1794 
1795 
1796 return ChargeSum*fChargeNormFactor/ReBin;
1797 }
x
Float_t x
Definition: compare.C:6
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

148 {
149  if (!moduleContext_)
150  return ProductToken<T>::invalid();
151 
152  consumables_[BT].emplace_back(ConsumableType::Product,
153  TypeID{typeid(T)},
154  it.label(),
155  it.instance(),
156  it.process());
157  return ProductToken<T>{it};
158 }
art::ProductToken::invalid
static ProductToken< T > invalid()
Definition: ProductToken.h:47
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 162 of file Consumer.h.

163 {
164  if (!moduleContext_)
165  return;
166 
167  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
168 }
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

173 {
174  if (!moduleContext_)
175  return ViewToken<T>::invalid();
176 
177  consumables_[BT].emplace_back(ConsumableType::ViewElement,
178  TypeID{typeid(T)},
179  it.label(),
180  it.instance(),
181  it.process());
182  return ViewToken<T>{it};
183 }
art::ViewToken::invalid
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed)
inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

27 {
28  return rng()->createEngine(placeholder_schedule_id(), seed);
29 }
seed
long seed
Definition: chem4.cc:68
art::EngineCreator::rng
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
Definition: EngineCreator.cc:58
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

34 {
35  return rng()->createEngine(
36  placeholder_schedule_id(), seed, kind_of_engine_to_make);
37 }
seed
long seed
Definition: chem4.cc:68
art::EngineCreator::rng
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
Definition: EngineCreator.cc:58
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

43 {
44  return rng()->createEngine(
45  placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
46 }
seed
long seed
Definition: chem4.cc:68
art::EngineCreator::rng
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
Definition: EngineCreator.cc:58
std::string hit::DPRawHitFinder::CreateFitFunction ( int  fNPeaks,
bool  fSameShape 
)
private

Definition at line 1534 of file DPRawHitFinder_module.cc.

Referenced by FindPeakWithMaxDeviation(), and FitExponentials().

1535 {
1536  std::string feqn = ""; // string holding fit formula
1537  std::stringstream numConv;
1538 
1539  if(fSameShape)
1540  {
1541  for(int i = 0; i < fNPeaks; i++)
1542  {
1543  feqn.append("+( [");
1544  numConv.str("");
1545  numConv << 2*(i+1);
1546  feqn.append(numConv.str());
1547  feqn.append("] * exp(0.4*(x-[");
1548  numConv.str("");
1549  numConv << 2*(i+1)+1;
1550  feqn.append(numConv.str());
1551  feqn.append("])/[");
1552  numConv.str("");
1553  numConv << 0;
1554  feqn.append(numConv.str());
1555  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1556  numConv.str("");
1557  numConv << 2*(i+1)+1;
1558  feqn.append(numConv.str());
1559  feqn.append("])/[");
1560  numConv.str("");
1561  numConv << 1;
1562  feqn.append(numConv.str());
1563  feqn.append("]) ) )");
1564  }
1565  }
1566  else
1567  {
1568  for(int i = 0; i < fNPeaks; i++)
1569  {
1570  feqn.append("+( [");
1571  numConv.str("");
1572  numConv << 4*i+2;
1573  feqn.append(numConv.str());
1574  feqn.append("] * exp(0.4*(x-[");
1575  numConv.str("");
1576  numConv << 4*i+3;
1577  feqn.append(numConv.str());
1578  feqn.append("])/[");
1579  numConv.str("");
1580  numConv << 4*i;
1581  feqn.append(numConv.str());
1582  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1583  numConv.str("");
1584  numConv << 2*(i+1)+1;
1585  feqn.append(numConv.str());
1586  feqn.append("])/[");
1587  numConv.str("");
1588  numConv << 4*i+1;
1589  feqn.append(numConv.str());
1590  feqn.append("]) ) )");
1591  }
1592  }
1593 return feqn;
1594 }
CurrentProcessingContext const * art::EDProducer::currentContext ( ) const
protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

121  {
122  return current_context_.get();
123  }
art::EDProducer::current_context_
CPC_exempt_ptr current_context_
Definition: EDProducer.h:116
void hit::DPRawHitFinder::doBinAverage ( const std::vector< float > &  inputVec,
std::vector< float > &  outputVec,
size_t  binsToAverage 
) const
private

Definition at line 1800 of file DPRawHitFinder_module.cc.

References min.

Referenced by produce().

1803 {
1804  size_t halfBinsToAverage(binsToAverage/2);
1805 
1806  float runningSum(0.);
1807 
1808  for(size_t idx = 0; idx < halfBinsToAverage; idx++) runningSum += inputVec[idx];
1809 
1810  outputVec.resize(inputVec.size());
1811  std::vector<float>::iterator outputVecItr = outputVec.begin();
1812 
1813  // First pass through to build the erosion vector
1814  for(std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end(); inputItr++)
1815  {
1816  size_t startOffset = std::distance(inputVec.begin(),inputItr);
1817  size_t stopOffset = std::distance(inputItr,inputVec.end());
1818  size_t count = std::min(2 * halfBinsToAverage, std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1819 
1820  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1821  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1822 
1823  *outputVecItr++ = runningSum / float(count);
1824  }
1825 
1826  return;
1827 }
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
min
Int_t min
Definition: plot.C:26
void hit::DPRawHitFinder::endJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 321 of file DPRawHitFinder_module.cc.

322 {
323 
324 }
int hit::DPRawHitFinder::EstimateFluctuations ( const std::vector< float >  fsignalVec,
int  peakStart,
int  peakMean,
int  peakEnd 
)
private

Definition at line 1224 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks().

1228 {
1229  int NFluctuations=0;
1230 
1231  for(int j = peakMean-1; j >= peakStart; j--)
1232  {
1233  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1234 
1235  if(fsignalVec[j] > fsignalVec[j+1])
1236  {
1237  NFluctuations++;
1238  }
1239  }
1240 
1241  for(int j = peakMean+1; j <= peakEnd; j++)
1242  {
1243  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1244 
1245  if(fsignalVec[j] > fsignalVec[j-1])
1246  {
1247  NFluctuations++;
1248  }
1249  }
1250 
1251  return NFluctuations;
1252 }
void hit::DPRawHitFinder::FillOutHitParameterVector ( const std::vector< double > &  input,
std::vector< double > &  output 
)
private

Definition at line 248 of file DPRawHitFinder_module.cc.

References geo::GeometryCore::Nplanes().

250 {
251  if(input.size()==0)
252  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
253 
254  art::ServiceHandle<geo::Geometry> geom;
255  const unsigned int N_PLANES = geom->Nplanes();
256 
257  if(input.size()==1)
258  output.resize(N_PLANES,input[0]);
259  else if(input.size()==N_PLANES)
260  output = input;
261  else
262  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector size !=1 and !=N_PLANES.");
263 
264 }
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2792
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 1014 of file DPRawHitFinder_module.cc.

References fTicksToStopPeakFinder.

Referenced by produce().

1019 {
1020  // Need a minimum number of ticks to do any work here
1021  if (std::distance(startItr,stopItr) > 4)
1022  {
1023  // Find the highest peak in the range given
1024  auto maxItr = std::max_element(startItr, stopItr);
1025 
1026  float maxValue = *maxItr;
1027  int maxTime = std::distance(startItr,maxItr);
1028 
1029  if (maxValue >= PeakMin)
1030  {
1031  // backwards to find first bin for this candidate hit
1032  auto firstItr = std::distance(startItr,maxItr) > 2 ? maxItr - 1 : startItr;
1033  bool PeakStartIsHere = true;
1034 
1035  while(firstItr != startItr)
1036  {
1037  //Check for inflection point & ADC count <= 0
1038  PeakStartIsHere = true;
1039  for(int i=1; i <= fTicksToStopPeakFinder; i++)
1040  {
1041  if( *firstItr >= *(firstItr-i) )
1042  {
1043  PeakStartIsHere = false;
1044  break;
1045  }
1046 
1047  }
1048  if (*firstItr <= 0 || PeakStartIsHere) break;
1049  firstItr--;
1050  }
1051 
1052  int firstTime = std::distance(startItr,firstItr);
1053 
1054  // Recursive call to find all candidate hits earlier than this peak
1055  findCandidatePeaks(startItr, firstItr - 1, timeValsVec, PeakMin, firstTick);
1056 
1057  // forwards to find last bin for this candidate hit
1058  auto lastItr = std::distance(maxItr,stopItr) > 2 ? maxItr + 1 : stopItr - 1;
1059  bool PeakEndIsHere = true;
1060 
1061  while(lastItr != stopItr)
1062  {
1063  //Check for inflection point & ADC count <= 0
1064  PeakEndIsHere = true;
1065  for(int i=1; i <= fTicksToStopPeakFinder; i++)
1066  {
1067  if( *lastItr >= *(lastItr+i) )
1068  {
1069  PeakEndIsHere = false;
1070  break;
1071  }
1072  }
1073  if (*lastItr <= 0 || PeakEndIsHere) break;
1074  lastItr++;
1075  }
1076 
1077  int lastTime = std::distance(startItr,lastItr);
1078 
1079  // Now save this candidate's start and max time info
1080  timeValsVec.push_back(std::make_tuple(firstTick+firstTime,firstTick+maxTime,firstTick+lastTime));
1081 
1082  // Recursive call to find all candidate hits later than this peak
1083  findCandidatePeaks(lastItr + 1, stopItr, timeValsVec, PeakMin, firstTick + std::distance(startItr,lastItr + 1));
1084  }
1085  }
1086 
1087  return;
1088 }
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:1014
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 1470 of file DPRawHitFinder_module.cc.

References CreateFitFunction(), t1, and t2.

Referenced by produce().

1478 {
1479 // int size = fEndTime - fStartTime + 1;
1480 // if(fEndTime - fStartTime < 0){size = 0;}
1481 
1482  std::string eqn = CreateFitFunction(fNPeaks, fSameShape); // string for equation of Exponentials fit
1483 
1484  TF1 Exponentials("Exponentials",eqn.c_str(),fStartTime,fEndTime+1);
1485 
1486  for(size_t i=0; i < fparamVec.size(); i++)
1487  {
1488  Exponentials.SetParameter(i, fparamVec[i].first);
1489  }
1490 
1491  // ##########################################################################
1492  // ### Finding the peak with the max chi2 fit and signal ###
1493  // ##########################################################################
1494  double Chi2PerNDFPeak;
1495  double MaxPosDeviation;
1496  double MaxNegDeviation;
1497  int BinMaxPosDeviation;
1498  int BinMaxNegDeviation;
1499  for(int i = 0; i < fNPeaks; i++)
1500  {
1501  Chi2PerNDFPeak = 0.;
1502  MaxPosDeviation=0.;
1503  MaxNegDeviation=0.;
1504  BinMaxPosDeviation=0;
1505  BinMaxNegDeviation=0;
1506 
1507  for(int j = std::get<2>(fpeakVals.at(i)); j < std::get<3>(fpeakVals.at(i))+1; j++)
1508  {
1509  if( (Exponentials(j+0.5)-fSignalVector[j]) > MaxPosDeviation && j != std::get<0>(fpeakVals.at(i)) )
1510  {
1511  MaxPosDeviation = Exponentials(j+0.5)-fSignalVector[j];
1512  BinMaxPosDeviation = j;
1513  }
1514  if( (Exponentials(j+0.5)-fSignalVector[j]) < MaxNegDeviation && j != std::get<0>(fpeakVals.at(i)) )
1515  {
1516  MaxNegDeviation = Exponentials(j+0.5)-fSignalVector[j];
1517  BinMaxNegDeviation = j;
1518  }
1519  Chi2PerNDFPeak += pow((Exponentials(j+0.5)-fSignalVector[j])/sqrt(fSignalVector[j]),2);
1520  }
1521 
1522  if(BinMaxNegDeviation != 0)
1523  {
1524  Chi2PerNDFPeak /= static_cast<double>((std::get<3>(fpeakVals.at(i))-std::get<2>(fpeakVals.at(i))));
1525  fPeakDev.emplace_back(Chi2PerNDFPeak,i,BinMaxNegDeviation,BinMaxPosDeviation);
1526  }
1527  }
1528 
1529 std::sort(fPeakDev.begin(),fPeakDev.end(), [](std::tuple<double,int,int,int> const &t1, std::tuple<double,int,int,int> const &t2) {return std::get<0>(t1) > std::get<0>(t2);} );
1530 Exponentials.Delete();
1531 }
t1
TTree * t1
Definition: plottest35.C:26
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1534
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 1257 of file DPRawHitFinder_module.cc.

References CreateFitFunction(), fFitPeakMeanRange, fLogLevel, fMaxTau, fMinTau, max, and min.

Referenced by produce().

1265 {
1266  int size = fEndTime - fStartTime + 1;
1267  int NPeaks = fPeakVals.size();
1268 
1269  // #############################################
1270  // ### If size < 0 then set the size to zero ###
1271  // #############################################
1272  if(fEndTime - fStartTime < 0){size = 0;}
1273 
1274  // --- TH1D HitSignal ---
1275  TH1F hitSignal("hitSignal","",std::max(size,1),fStartTime,fEndTime+1);
1276  hitSignal.Sumw2();
1277 
1278  // #############################
1279  // ### Filling the histogram ###
1280  // #############################
1281  for(int i = fStartTime; i < fEndTime+1; i++)
1282  {
1283  hitSignal.Fill(i,fSignalVector[i]);
1284  hitSignal.SetBinError(i,0.288675); //1/sqrt(12)
1285  }
1286 
1287  // ################################################
1288  // ### Building TFormula for basic Exponentials ###
1289  // ################################################
1290  std::string eqn = CreateFitFunction(NPeaks, fSameShape);
1291 
1292  // -------------------------------------
1293  // --- TF1 function for Exponentials ---
1294  // -------------------------------------
1295 
1296  TF1 Exponentials("Exponentials",eqn.c_str(),fStartTime,fEndTime+1);
1297 
1298  if(fLogLevel >= 4)
1299  {
1300  std::cout << std::endl;
1301  std::cout << "--- Preparing fit ---" << std::endl;
1302  std::cout << "--- Lower limits, seed, upper limit:" << std::endl;
1303  }
1304 
1305  if(fSameShape)
1306  {
1307  Exponentials.SetParameter(0, 0.5);
1308  Exponentials.SetParameter(1, 0.5);
1309  Exponentials.SetParLimits(0, fMinTau, fMaxTau);
1310  Exponentials.SetParLimits(1, fMinTau, fMaxTau);
1311  double amplitude=0;
1312  double peakMean=0;
1313 
1314  double peakMeanShift=2;
1315  double peakMeanSeed=0;
1316  double peakMeanRangeLow=0;
1317  double peakMeanRangeHi=0;
1318  double peakStart=0;
1319  double peakEnd=0;
1320 
1321  for(int i = 0; i < NPeaks; i++)
1322  {
1323  peakMean = std::get<0>(fPeakVals.at(i));
1324  peakStart = std::get<2>(fPeakVals.at(i));
1325  peakEnd = std::get<3>(fPeakVals.at(i));
1326  peakMeanSeed=peakMean-peakMeanShift;
1327  peakMeanRangeLow = std::max(peakStart-peakMeanShift, peakMeanSeed-fFitPeakMeanRange);
1328  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed+fFitPeakMeanRange);
1329  amplitude = fSignalVector[peakMean];
1330 
1331  Exponentials.SetParameter(2*(i+1), 1.65*amplitude);
1332  Exponentials.SetParLimits(2*(i+1), 0.3*1.65*amplitude, 2*1.65*amplitude);
1333  Exponentials.SetParameter(2*(i+1)+1, peakMeanSeed);
1334 
1335  if(NPeaks == 1)
1336  {
1337  Exponentials.SetParLimits(2*(i+1)+1, peakMeanRangeLow, peakMeanRangeHi);
1338  }
1339  else if(NPeaks >= 2 && i == 0)
1340  {
1341  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1342  Exponentials.SetParLimits( 2*(i+1)+1, peakMeanRangeLow, std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1343  }
1344  else if(NPeaks >= 2 && i == NPeaks-1)
1345  {
1346  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1347  Exponentials.SetParLimits(2*(i+1)+1, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), peakMeanRangeHi );
1348  }
1349  else
1350  {
1351  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1352  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1353  Exponentials.SetParLimits(2*(i+1)+1, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1354  }
1355 
1356  if(fLogLevel >= 4)
1357  {
1358  double t0low, t0high;
1359  Exponentials.GetParLimits(2*(i+1)+1, t0low, t0high);
1360  std::cout << "Peak #" << i << ": A [ADC] = " << std::setprecision(1) << std::fixed << 0.3*1.65*amplitude << " , " << 1.65*amplitude << " , " << 2*1.65*amplitude << std::endl;
1361  std::cout << "Peak #" << i << ": t0 [ticks] = " << std::setprecision(1) << std::fixed << t0low << " , " << peakMeanSeed << " , " << t0high << std::endl;
1362  }
1363  }
1364  }
1365  else
1366  {
1367  double amplitude=0;
1368  double peakMean=0;
1369 
1370  double peakMeanShift=2;
1371  double peakMeanSeed=0;
1372  double peakMeanRangeLow=0;
1373  double peakMeanRangeHi=0;
1374  double peakStart=0;
1375  double peakEnd=0;
1376 
1377 
1378  for(int i = 0; i < NPeaks; i++)
1379  {
1380  Exponentials.SetParameter(4*i, 0.5);
1381  Exponentials.SetParameter(4*i+1, 0.5);
1382  Exponentials.SetParLimits(4*i, fMinTau, fMaxTau);
1383  Exponentials.SetParLimits(4*i+1, fMinTau, fMaxTau);
1384 
1385  peakMean = std::get<0>(fPeakVals.at(i));
1386  peakStart = std::get<2>(fPeakVals.at(i));
1387  peakEnd = std::get<3>(fPeakVals.at(i));
1388  peakMeanSeed=peakMean-peakMeanShift;
1389  peakMeanRangeLow = std::max(peakStart-peakMeanShift, peakMeanSeed-fFitPeakMeanRange);
1390  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed+fFitPeakMeanRange);
1391  amplitude = fSignalVector[peakMean];
1392 
1393  Exponentials.SetParameter(4*i+2, 1.65*amplitude);
1394  Exponentials.SetParLimits(4*i+2, 0.3*1.65*amplitude, 2*1.65*amplitude);
1395  Exponentials.SetParameter(4*i+3, peakMeanSeed);
1396 
1397  if(NPeaks == 1)
1398  {
1399  Exponentials.SetParLimits(4*i+3, peakMeanRangeLow, peakMeanRangeHi);
1400  }
1401  else if(NPeaks >= 2 && i == 0)
1402  {
1403  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1404  Exponentials.SetParLimits( 4*i+3, peakMeanRangeLow, std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1405  }
1406  else if(NPeaks >= 2 && i == NPeaks-1)
1407  {
1408  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1409  Exponentials.SetParLimits(4*i+3, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), peakMeanRangeHi );
1410  }
1411  else
1412  {
1413  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1414  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1415  Exponentials.SetParLimits(4*i+3, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1416  }
1417 
1418  if(fLogLevel >= 4)
1419  {
1420  double t0low, t0high;
1421  Exponentials.GetParLimits(4*i+3, t0low, t0high);
1422  std::cout << "Peak #" << i << ": A [ADC] = " << std::setprecision(1) << std::fixed << 0.3*1.65*amplitude << " , " << 1.65*amplitude << " , " << 2*1.65*amplitude << std::endl;
1423  std::cout << "Peak #" << i << ": t0 [ticks] = " << std::setprecision(1) << std::fixed << t0low << " , " << peakMeanSeed << " , " << t0high << std::endl;
1424  }
1425  }
1426  }
1427 
1428 
1429  // ###########################################
1430  // ### PERFORMING THE TOTAL FIT OF THE HIT ###
1431  // ###########################################
1432  try
1433  { hitSignal.Fit(&Exponentials,"QNRWM","", fStartTime, fEndTime+1);}
1434  catch(...)
1435  {mf::LogWarning("DPRawHitFinder") << "Fitter failed finding a hit";}
1436 
1437  // ##################################################
1438  // ### Getting the fitted parameters from the fit ###
1439  // ##################################################
1440  fchi2PerNDF = (Exponentials.GetChisquare() / Exponentials.GetNDF());
1441  fNDF = Exponentials.GetNDF();
1442 
1443  if(fSameShape)
1444  {
1445  fparamVec.emplace_back(Exponentials.GetParameter(0),Exponentials.GetParError(0));
1446  fparamVec.emplace_back(Exponentials.GetParameter(1),Exponentials.GetParError(1));
1447 
1448  for(int i = 0; i < NPeaks; i++)
1449  {
1450  fparamVec.emplace_back(Exponentials.GetParameter(2*(i+1)),Exponentials.GetParError(2*(i+1)));
1451  fparamVec.emplace_back(Exponentials.GetParameter(2*(i+1)+1),Exponentials.GetParError(2*(i+1)+1));
1452  }
1453  }
1454  else
1455  {
1456  for(int i = 0; i < NPeaks; i++)
1457  {
1458  fparamVec.emplace_back(Exponentials.GetParameter(4*i),Exponentials.GetParError(4*i));
1459  fparamVec.emplace_back(Exponentials.GetParameter(4*i+1),Exponentials.GetParError(4*i+1));
1460  fparamVec.emplace_back(Exponentials.GetParameter(4*i+2),Exponentials.GetParError(4*i+2));
1461  fparamVec.emplace_back(Exponentials.GetParameter(4*i+3),Exponentials.GetParError(4*i+3));
1462  }
1463  }
1464  Exponentials.Delete();
1465  hitSignal.Delete();
1466 }//<----End FitExponentials
max
Int_t max
Definition: plot.C:27
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1534
min
Int_t min
Definition: plot.C:26
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
EngineCreator::seed_t EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

52 {
53  auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
54  return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
55 }
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
art::EDProducer::moduleDescription_
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
B
Int_t B
Definition: plot.C:25
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

191 {
192  if (!moduleContext_)
193  return ProductToken<T>::invalid();
194 
195  consumables_[BT].emplace_back(ConsumableType::Product,
196  TypeID{typeid(T)},
197  it.label(),
198  it.instance(),
199  it.process());
200  return ProductToken<T>{it};
201 }
art::ProductToken::invalid
static ProductToken< T > invalid()
Definition: ProductToken.h:47
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 205 of file Consumer.h.

206 {
207  if (!moduleContext_)
208  return;
209 
210  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
211 }
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

216 {
217  if (!moduleContext_)
218  return ViewToken<T>::invalid();
219 
220  consumables_[BT].emplace_back(ConsumableType::ViewElement,
221  TypeID{typeid(T)},
222  it.label(),
223  it.instance(),
224  it.process());
225  return ViewToken<T>{it};
226 }
art::ViewToken::invalid
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
void hit::DPRawHitFinder::mergeCandidatePeaks ( const std::vector< float >  signalVec,
TimeValsVec  timeValsVec,
MergedTimeWidVec &  mergedVec 
)
private

Definition at line 1094 of file DPRawHitFinder_module.cc.

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

Referenced by produce().

1095 {
1096  // ################################################################
1097  // ### Lets loop over the candidate pulses we found in this ROI ###
1098  // ################################################################
1099  auto timeValsVecItr = timeValsVec.begin();
1100  unsigned int PeaksInThisMergedPeak = 0;
1101  //int EndTickOfPreviousMergedPeak=0;
1102  while(timeValsVecItr != timeValsVec.end())
1103  {
1104  PeakTimeWidVec peakVals;
1105 
1106  // Setting the start, peak, and end time of the pulse
1107  auto& timeVal = *timeValsVecItr++;
1108  int startT = std::get<0>(timeVal);
1109  int maxT = std::get<1>(timeVal);
1110  int endT = std::get<2>(timeVal);
1111  int widT = endT - startT;
1112  int FinalStartT=startT;
1113  int FinalEndT=endT;
1114 
1115  int NFluctuations=EstimateFluctuations(signalVec, startT, maxT, endT);
1116 
1117  peakVals.emplace_back(maxT,widT,startT,endT);
1118 
1119  // See if we want to merge pulses together
1120  // First check if we have more than one pulse on the wire
1121  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1122 
1123  // Loop until no more merged pulses (or candidates in this ROI)
1124  while(checkNextHit)
1125  {
1126  // 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
1127  int NextStartT = std::get<0>(*timeValsVecItr);
1128 
1129  double MinADC = signalVec[endT];
1130  for(int i = endT; i <= NextStartT; i++)
1131  {
1132  if(signalVec[i]<MinADC)
1133  {
1134  MinADC = signalVec[i];
1135  }
1136  }
1137 
1138  // Group peaks (grouped peaks are fitted together and can be merged)
1139  if( MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance )
1140  {
1141  int CurrentStartT=startT;
1142  int CurrentMaxT=maxT;
1143  int CurrentEndT=endT;
1144  //int CurrentWidT=widT;
1145 
1146  timeVal = *timeValsVecItr++;
1147  int NextMaxT = std::get<1>(timeVal);
1148  int NextEndT = std::get<2>(timeVal);
1149  int NextWidT = NextEndT - NextStartT;
1150  FinalEndT=NextEndT;
1151 
1152  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1153 
1154  int CurrentSumADC = 0;
1155  for(int i = CurrentStartT; i<= CurrentEndT; i++)
1156  {
1157  CurrentSumADC+=signalVec[i];
1158  }
1159 
1160  int NextSumADC = 0;
1161  for (int i = NextStartT; i<= NextEndT; i++)
1162  {
1163  NextSumADC+=signalVec[i];
1164  }
1165 
1166  //Merge peaks within a group
1167  if(fDoMergePeaks)
1168  {
1169  if( signalVec[NextMaxT] <= signalVec[CurrentMaxT] && ( ( signalVec[NextMaxT] < fMergeMaxADCThreshold*signalVec[CurrentMaxT] && NextSumADC < fMergeADCSumThreshold*CurrentSumADC ) || (signalVec[NextMaxT]-signalVec[NextStartT]) < fMinRelativePeakHeightRight*signalVec[CurrentMaxT] ) && ( signalVec[NextMaxT] <= fMergeMaxADCLimit ) )
1170  {
1171  maxT=CurrentMaxT;
1172  startT=CurrentStartT;
1173  endT=NextEndT;
1174  widT=endT - startT;
1175  peakVals.pop_back();
1176  peakVals.emplace_back(maxT,widT,startT,endT);
1177  }
1178  else if( signalVec[NextMaxT] > signalVec[CurrentMaxT] && ( ( signalVec[CurrentMaxT] < fMergeMaxADCThreshold*signalVec[NextMaxT] && CurrentSumADC < fMergeADCSumThreshold*NextSumADC ) || (signalVec[CurrentMaxT]-signalVec[CurrentEndT]) < fMinRelativePeakHeightLeft*signalVec[NextMaxT] ) && ( signalVec[CurrentMaxT] <= fMergeMaxADCLimit ) )
1179  {
1180  maxT=NextMaxT;
1181  startT=CurrentStartT;
1182  endT=NextEndT;
1183  widT=endT - startT;
1184  peakVals.pop_back();
1185  peakVals.emplace_back(maxT,widT,startT,endT);
1186  }
1187  else
1188  {
1189  maxT=NextMaxT;
1190  startT=NextStartT;
1191  endT=NextEndT;
1192  widT=NextWidT;
1193  peakVals.emplace_back(maxT,widT,startT,endT);
1194  PeaksInThisMergedPeak++;
1195  }
1196  }
1197  else
1198  {
1199  maxT=NextMaxT;
1200  startT=NextStartT;
1201  endT=NextEndT;
1202  widT=NextWidT;
1203  peakVals.emplace_back(maxT,widT,startT,endT);
1204  PeaksInThisMergedPeak++;
1205  }
1206  checkNextHit = timeValsVecItr != timeValsVec.end();
1207  }//<---Checking adjacent pulses
1208  else
1209  {
1210  checkNextHit = false;
1211  PeaksInThisMergedPeak = 0;
1212  }
1213 
1214  }//<---End checking if there is more than one pulse on the wire
1215  // Add these to our merged vector
1216  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1217  }
1218  return;
1219 }
hit::DPRawHitFinder::PeakTimeWidVec
std::vector< std::tuple< int, int, int, int >> PeakTimeWidVec
Definition: DPRawHitFinder_module.cc:93
hit::DPRawHitFinder::EstimateFluctuations
int EstimateFluctuations(const std::vector< float > fsignalVec, int peakStart, int peakMean, int peakEnd)
Definition: DPRawHitFinder_module.cc:1224
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

41  {
42  return true;
43  }
cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )
staticinherited

Definition at line 76 of file Consumer.cc.

Referenced by art::RootOutput::beginSubRun(), art::OutputModule::doBeginRun(), art::OutputModule::doBeginSubRun(), art::OutputModule::doEndRun(), art::OutputModule::doEndSubRun(), art::ProducingService::doPostReadEvent(), art::ProducingService::doPostReadRun(), art::ProducingService::doPostReadSubRun(), art::OutputModule::doWriteEvent(), art::ProcessPackage< L >::postScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::Run >::End::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::End::postScheduleSignal(), art::ProcessPackage< L >::preScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::preScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::preScheduleSignal(), art::EventProcessor::readEvent(), art::EventProcessor::readRun(), art::EmptyEvent::readRun_(), art::EventProcessor::readSubRun(), and art::EmptyEvent::readSubRun_().

77 {
78  static Consumer invalid{InvalidTag{}};
79  return &invalid;
80 }
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

90 {
91  if (!moduleContext_)
92  return;
93 
94  pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
95  for (auto& consumablesPerBranch : consumables_) {
96  cet::sort_all(consumablesPerBranch);
97  }
98 }
art::Consumer::requireConsumes_
bool requireConsumes_
Definition: Consumer.h:137
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
void hit::DPRawHitFinder::produce ( art::Event &  evt)
overridevirtual

Implements art::EDProducer.

Definition at line 327 of file DPRawHitFinder_module.cc.

References AddPeak(), anab::FVectorWriter< N >::addVector(), recob::Wire::Channel(), geo::GeometryCore::ChannelToWire(), ChargeFunc(), geo::CryostatID::Cryostat, doBinAverage(), recob::HitCollectionCreator::emplace_back(), fCalDataModuleLabel, fChargeNorm, fChi2, fChi2NDFMax, fChi2NDFMaxFactorMultiHits, fChi2NDFRetry, fChi2NDFRetryFactorMultiHits, fFirstChi2, fHitParamWriter, findCandidatePeaks(), FindPeakWithMaxDeviation(), FitExponentials(), fLogLevel, fLongMaxHits, fLongPulseWidth, fMaxFluctuations, fMaxGroupLength, fMaxMultiHit, fMinADCSum, fMinADCSumOverWidth, fMinSig, fMinWidth, fNewHitsTag, fNumBinsToAverage, fSameShape, fTryNplus1Fits, fWidthNormalization, lar::sparse_vector< T >::get_ranges(), art::DataViewImpl::getByLabel(), anab::FVectorWriter< N >::initOutputs(), raw::InvalidChannelID, mergeCandidatePeaks(), min, 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.

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

Definition at line 1830 of file DPRawHitFinder_module.cc.

References DEFINE_ART_MODULE.

1833 {
1834  size_t nNewBins = inputVec.size() / nBinsToCombine;
1835 
1836  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1837 
1838  outputVec.resize(nNewBins, 0.);
1839 
1840  size_t outputBin = 0;
1841 
1842  for(size_t inputIdx = 0; inputIdx < inputVec.size();)
1843  {
1844  outputVec[outputBin] += inputVec[inputIdx++];
1845 
1846  if (inputIdx % nBinsToCombine == 0) outputBin++;
1847 
1848  if (outputBin > outputVec.size())
1849  {
1850  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin << ", inputIdx = " << inputIdx << std::endl;
1851  break;
1852  }
1853  }
1854 
1855  return;
1856 }
void hit::DPRawHitFinder::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 269 of file DPRawHitFinder_module.cc.

References fCalDataModuleLabel, fChargeNorm, fChi2NDFMax, fChi2NDFMaxFactorMultiHits, fChi2NDFRetry, fChi2NDFRetryFactorMultiHits, fDoMergePeaks, fFitPeakMeanRange, fGroupMaxDistance, fGroupMinADC, fLogLevel, fLongMaxHits, fLongPulseWidth, fMaxFluctuations, fMaxGroupLength, fMaxMultiHit, fMaxTau, fMergeADCSumThreshold, fMergeMaxADCLimit, fMergeMaxADCThreshold, fMinADCSum, fMinADCSumOverWidth, fMinRelativePeakHeightLeft, fMinRelativePeakHeightRight, fMinSig, fMinTau, fMinWidth, fNumBinsToAverage, fSameShape, fTicksToStopPeakFinder, fTryNplus1Fits, fWidthNormalization, and fhicl::ParameterSet::get().

Referenced by DPRawHitFinder().

270 {
271  // Implementation of optional member function here.
272  fLogLevel = p.get< int >("LogLevel");
273  fCalDataModuleLabel = p.get< std::string >("CalDataModuleLabel");
274  fMaxMultiHit = p.get< int >("MaxMultiHit");
275  fMaxGroupLength = p.get< int >("MaxGroupLength");
276  fTryNplus1Fits = p.get< bool >("TryNplus1Fits");
277  fChi2NDFRetry = p.get< double >("Chi2NDFRetry");
278  fChi2NDFRetryFactorMultiHits = p.get< double >("Chi2NDFRetryFactorMultiHits");
279  fChi2NDFMax = p.get< double >("Chi2NDFMax");
280  fChi2NDFMaxFactorMultiHits = p.get< double >("Chi2NDFMaxFactorMultiHits");
281  fNumBinsToAverage = p.get< size_t >("NumBinsToAverage");
282  fMinSig = p.get< float >("MinSig");
283  fMinWidth = p.get< int >("MinWidth");
284  fMinADCSum = p.get< double >("MinADCSum");
285  fMinADCSumOverWidth = p.get< double >("MinADCSumOverWidth");
286  fChargeNorm = p.get< double >("ChargeNorm");
287  fTicksToStopPeakFinder = p.get< double >("TicksToStopPeakFinder");
288  fMinTau = p.get< double >("MinTau");
289  fMaxTau = p.get< double >("MaxTau");
290  fFitPeakMeanRange = p.get< double >("FitPeakMeanRange");
291  fGroupMaxDistance = p.get< int >("GroupMaxDistance");
292  fGroupMinADC = p.get< int >("GroupMinADC");
293  fSameShape = p.get< bool >("SameShape");
294  fDoMergePeaks = p.get< bool >("DoMergePeaks");
295  fMergeADCSumThreshold = p.get< double >("MergeADCSumThreshold");
296  fMergeMaxADCThreshold = p.get< double >("MergeMaxADCThreshold");
297  fMinRelativePeakHeightLeft = p.get< double >("MinRelativePeakHeightLeft");
298  fMinRelativePeakHeightRight = p.get< double >("MinRelativePeakHeightRight");
299  fMergeMaxADCLimit = p.get< double >("MergeMaxADCLimit");
300  fWidthNormalization = p.get< double >("WidthNormalization");
301  fLongMaxHits = p.get< double >("LongMaxHits");
302  fLongPulseWidth = p.get< double >("LongPulseWidth");
303  fMaxFluctuations = p.get< double >("MaxFluctuations");
304 }
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

126 {
127  if (!moduleContext_)
128  return;
129 
130  // If none of the branches have missing consumes statements, exit early.
131  if (std::all_of(cbegin(missingConsumes_),
132  cend(missingConsumes_),
133  [](auto const& perBranch) { return perBranch.empty(); }))
134  return;
135 
136  constexpr cet::HorizontalRule rule{60};
137  mf::LogPrint log{"MTdiagnostics"};
138  log << '\n'
139  << rule('=') << '\n'
140  << "The following consumes (or mayConsume) statements are missing from\n"
141  << module_context(moduleDescription_) << '\n'
142  << rule('-') << '\n';
143 
144  cet::for_all_with_index(
145  missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
146  for (auto const& pi : perBranch) {
147  log << " "
148  << assemble_consumes_statement(static_cast<BranchType>(i), pi)
149  << '\n';
150  }
151  });
152  log << rule('=');
153 }
mf::MaybeLogger_
Definition: MessageLogger.h:114
art::Consumer::moduleDescription_
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
util::pi
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
Definition: PhysicalConstants.h:80
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
art::Consumer::missingConsumes_
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139
void hit::DPRawHitFinder::SplitPeak ( std::tuple< double, int, int, int >  fPeakDevCand,
PeakTimeWidVec &  fpeakValsTemp 
)
private

Definition at line 1642 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1644 {
1645 int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1646 int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1647 int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1648 int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1649 
1650 if(WidthOldPeakOld<3) {return;} //can't split peaks with widths < 3
1651 
1652 int NewPeakMax = 0;
1653 int NewPeakStart = 0;
1654 int NewPeakEnd = 0;
1655 int OldPeakNewMax = 0;
1656 int OldPeakNewStart = 0;
1657 int OldPeakNewEnd = 0;
1658 
1659 
1660 OldPeakNewStart = OldPeakOldStart;
1661 NewPeakEnd = OldPeakOldEnd;
1662 
1663 OldPeakNewEnd = OldPeakNewStart + 0.5*(WidthOldPeakOld + (WidthOldPeakOld%2));
1664 NewPeakStart = OldPeakNewEnd+1;
1665 
1666 int WidthOldPeakNew = OldPeakNewEnd-OldPeakNewStart;
1667 int WidthNewPeak = NewPeakEnd-NewPeakStart;
1668 
1669 OldPeakNewMax = OldPeakNewStart + 0.5*(WidthOldPeakNew - (WidthOldPeakNew%2));
1670 NewPeakMax = NewPeakStart + 0.5*(WidthNewPeak - (WidthNewPeak%2));
1671 
1672 fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakNewMax,0,OldPeakNewStart,OldPeakNewEnd);
1673 fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1674 std::sort(fpeakValsTemp.begin(),fpeakValsTemp.end(), [](std::tuple<int,int,int,int> const &t1, std::tuple<int,int,int,int> const &t2) {return std::get<0>(t1) < std::get<0>(t2);} );
1675 
1676 return;
1677 }
t1
TTree * t1
Definition: plottest35.C:26
t2
TTree * t2
Definition: plottest35.C:36
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

103 {
104  // Early exits if consumes tracking has been disabled or if the
105  // consumed product is an allowed consumable.
106  if (!moduleContext_)
107  return;
108 
109  if (cet::binary_search_all(consumables_[bt], pi))
110  return;
111 
112  if (requireConsumes_) {
113  throw Exception(errors::ProductRegistrationFailure,
114  "Consumer: an error occurred during validation of a "
115  "retrieved product\n\n")
116  << "The following consumes (or mayConsume) statement is missing from\n"
117  << module_context(moduleDescription_) << ":\n\n"
118  << " " << assemble_consumes_statement(bt, pi) << "\n\n";
119  }
120 
121  missingConsumes_[bt].insert(pi);
122 }
art::Consumer::moduleDescription_
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
art::Consumer::requireConsumes_
bool requireConsumes_
Definition: Consumer.h:137
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
util::pi
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
Definition: PhysicalConstants.h:80
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:138
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:136
art::Consumer::missingConsumes_
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139
double hit::DPRawHitFinder::WidthFunc ( double  fPeakMean,
double  fPeakAmp,
double  fPeakTau1,
double  fPeakTau2,
double  fStartTime,
double  fEndTime,
double  fPeakMeanTrue 
)
private

Definition at line 1680 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1687 {
1688 double MaxValue = ( fPeakAmp * exp(0.4*(fPeakMeanTrue-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(fPeakMeanTrue-fPeakMean)/fPeakTau2) );
1689 double FuncValue = 0.;
1690 double HalfMaxLeftTime = 0.;
1691 double HalfMaxRightTime = 0.;
1692 double ReBin=10.;
1693 
1694  //First iteration (+- 1 bin)
1695  for(double x = fPeakMeanTrue; x > fStartTime-1000.; x--)
1696  {
1697  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1698  if( FuncValue < 0.5*MaxValue )
1699  {
1700  HalfMaxLeftTime = x;
1701  break;
1702  }
1703  }
1704 
1705  for(double x = fPeakMeanTrue; x < fEndTime+1000.; x++)
1706  {
1707  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1708  if( FuncValue < 0.5*MaxValue )
1709  {
1710  HalfMaxRightTime = x;
1711  break;
1712  }
1713  }
1714 
1715  //Second iteration (+- 0.1 bin)
1716  for(double x = HalfMaxLeftTime+1; x > HalfMaxLeftTime; x-=(1/ReBin))
1717  {
1718  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1719  if( FuncValue < 0.5*MaxValue )
1720  {
1721  HalfMaxLeftTime = x;
1722  break;
1723  }
1724  }
1725 
1726  for(double x = HalfMaxRightTime-1; x < HalfMaxRightTime; x+=(1/ReBin))
1727  {
1728  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1729  if( FuncValue < 0.5*MaxValue )
1730  {
1731  HalfMaxRightTime = x;
1732  break;
1733  }
1734  }
1735 
1736  //Third iteration (+- 0.01 bin)
1737  for(double x = HalfMaxLeftTime+1/ReBin; x > HalfMaxLeftTime; x-=1/(ReBin*ReBin))
1738  {
1739  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1740  if( FuncValue < 0.5*MaxValue )
1741  {
1742  HalfMaxLeftTime = x;
1743  break;
1744  }
1745  }
1746 
1747  for(double x = HalfMaxRightTime-1/ReBin; x < HalfMaxRightTime; x+=1/(ReBin*ReBin))
1748  {
1749  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1750  if( FuncValue < 0.5*MaxValue )
1751  {
1752  HalfMaxRightTime = x;
1753  break;
1754  }
1755  }
1756 
1757 return HalfMaxRightTime-HalfMaxLeftTime;
1758 }
x
Float_t x
Definition: compare.C:6

Member Data Documentation

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

Definition at line 177 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fChargeNorm
private

Definition at line 188 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

TH1F* hit::DPRawHitFinder::fChi2
private

Definition at line 216 of file DPRawHitFinder_module.cc.

Referenced by beginJob(), and produce().

double hit::DPRawHitFinder::fChi2NDFMax
private

Definition at line 192 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fChi2NDFMaxFactorMultiHits
private

Definition at line 193 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fChi2NDFRetry
private

Definition at line 190 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fChi2NDFRetryFactorMultiHits
private

Definition at line 191 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

bool hit::DPRawHitFinder::fDoMergePeaks
private

Definition at line 201 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

TH1F* hit::DPRawHitFinder::fFirstChi2
private

Definition at line 215 of file DPRawHitFinder_module.cc.

Referenced by beginJob(), and produce().

double hit::DPRawHitFinder::fFitPeakMeanRange
private

Definition at line 197 of file DPRawHitFinder_module.cc.

Referenced by FitExponentials(), and reconfigure().

int hit::DPRawHitFinder::fGroupMaxDistance
private

Definition at line 198 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

double hit::DPRawHitFinder::fGroupMinADC
private

Definition at line 199 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

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

Definition at line 213 of file DPRawHitFinder_module.cc.

Referenced by DPRawHitFinder(), and produce().

int hit::DPRawHitFinder::fLogLevel
private

Definition at line 180 of file DPRawHitFinder_module.cc.

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

int hit::DPRawHitFinder::fLongMaxHits
private

Definition at line 208 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

int hit::DPRawHitFinder::fLongPulseWidth
private

Definition at line 209 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

int hit::DPRawHitFinder::fMaxFluctuations
private

Definition at line 210 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

int hit::DPRawHitFinder::fMaxGroupLength
private

Definition at line 187 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

int hit::DPRawHitFinder::fMaxMultiHit
private

Definition at line 186 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fMaxTau
private

Definition at line 196 of file DPRawHitFinder_module.cc.

Referenced by FitExponentials(), and reconfigure().

double hit::DPRawHitFinder::fMergeADCSumThreshold
private

Definition at line 202 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

double hit::DPRawHitFinder::fMergeMaxADCLimit
private

Definition at line 206 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

double hit::DPRawHitFinder::fMergeMaxADCThreshold
private

Definition at line 203 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

double hit::DPRawHitFinder::fMinADCSum
private

Definition at line 184 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fMinADCSumOverWidth
private

Definition at line 185 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fMinRelativePeakHeightLeft
private

Definition at line 204 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

double hit::DPRawHitFinder::fMinRelativePeakHeightRight
private

Definition at line 205 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks(), and reconfigure().

float hit::DPRawHitFinder::fMinSig
private

Definition at line 181 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fMinTau
private

Definition at line 195 of file DPRawHitFinder_module.cc.

Referenced by FitExponentials(), and reconfigure().

int hit::DPRawHitFinder::fMinWidth
private

Definition at line 183 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

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

Definition at line 212 of file DPRawHitFinder_module.cc.

Referenced by produce().

size_t hit::DPRawHitFinder::fNumBinsToAverage
private

Definition at line 194 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

bool hit::DPRawHitFinder::fSameShape
private

Definition at line 200 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

int hit::DPRawHitFinder::fTicksToStopPeakFinder
private

Definition at line 182 of file DPRawHitFinder_module.cc.

Referenced by findCandidatePeaks(), and reconfigure().

bool hit::DPRawHitFinder::fTryNplus1Fits
private

Definition at line 189 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

double hit::DPRawHitFinder::fWidthNormalization
private

Definition at line 207 of file DPRawHitFinder_module.cc.

Referenced by produce(), and reconfigure().

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