LArSoft  v07_13_02
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 1599 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1601 {
1602  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1603  int NewPeakMax = std::get<2>(fPeakDevCand);
1604  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1605  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1606  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1607 
1608  int NewPeakStart=0;
1609  int NewPeakEnd=0;
1610  int OldPeakNewStart=0;
1611  int OldPeakNewEnd=0;
1612  int DistanceBwOldAndNewPeak=0;
1613 
1614  if(NewPeakMax<OldPeakMax)
1615  {
1616  NewPeakStart = OldPeakOldStart;
1617  OldPeakNewEnd = OldPeakOldEnd;
1618  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1619  NewPeakEnd = NewPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1620  if(DistanceBwOldAndNewPeak%2 == 0) NewPeakEnd -= 1;
1621  OldPeakNewStart = NewPeakEnd+1;
1622  }
1623  else if(OldPeakMax<NewPeakMax)
1624  {
1625  NewPeakEnd = OldPeakOldEnd;
1626  OldPeakNewStart = OldPeakOldStart;
1627  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1628  OldPeakNewEnd = OldPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1629  if(DistanceBwOldAndNewPeak%2 == 0) OldPeakNewEnd -= 1;
1630  NewPeakStart = OldPeakNewEnd+1;
1631  }
1632  else if(OldPeakMax==NewPeakMax){return;} //This shouldn't happen
1633 
1634  fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakMax,0,OldPeakNewStart,OldPeakNewEnd);
1635  fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1636  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);} );
1637 
1638  return;
1639 }
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 1762 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

1769 {
1770 double ChargeSum = 0.;
1771 double Charge = 0.;
1772 double ReBin = 10.;
1773 
1774 bool ChargeBigEnough=true;
1775  for(double x = fPeakMeanTrue - 1/ReBin; ChargeBigEnough && x > fPeakMeanTrue-1000.; x-=1.)
1776  {
1777  for(double i=0.; i > -1.; i-=(1/ReBin))
1778  {
1779  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1780  ChargeSum += Charge;
1781  }
1782  if(Charge < 0.01) ChargeBigEnough = false;
1783  }
1784 
1785 ChargeBigEnough=true;
1786  for(double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue+1000.; x+=1.)
1787  {
1788  for(double i=0.; i < 1.; i+=(1/ReBin))
1789  {
1790  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1791  ChargeSum += Charge;
1792  }
1793  if(Charge < 0.01) ChargeBigEnough = false;
1794  }
1795 
1796 
1797 return ChargeSum*fChargeNormFactor/ReBin;
1798 }
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 1535 of file DPRawHitFinder_module.cc.

Referenced by FindPeakWithMaxDeviation(), and FitExponentials().

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

References min.

Referenced by produce().

1804 {
1805  size_t halfBinsToAverage(binsToAverage/2);
1806 
1807  float runningSum(0.);
1808 
1809  for(size_t idx = 0; idx < halfBinsToAverage; idx++) runningSum += inputVec[idx];
1810 
1811  outputVec.resize(inputVec.size());
1812  std::vector<float>::iterator outputVecItr = outputVec.begin();
1813 
1814  // First pass through to build the erosion vector
1815  for(std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end(); inputItr++)
1816  {
1817  size_t startOffset = std::distance(inputVec.begin(),inputItr);
1818  size_t stopOffset = std::distance(inputItr,inputVec.end());
1819  size_t count = std::min(2 * halfBinsToAverage, std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1820 
1821  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1822  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1823 
1824  *outputVecItr++ = runningSum / float(count);
1825  }
1826 
1827  return;
1828 }
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 1225 of file DPRawHitFinder_module.cc.

Referenced by mergeCandidatePeaks().

1229 {
1230  int NFluctuations=0;
1231 
1232  for(int j = peakMean-1; j >= peakStart; j--)
1233  {
1234  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1235 
1236  if(fsignalVec[j] > fsignalVec[j+1])
1237  {
1238  NFluctuations++;
1239  }
1240  }
1241 
1242  for(int j = peakMean+1; j <= peakEnd; j++)
1243  {
1244  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1245 
1246  if(fsignalVec[j] > fsignalVec[j-1])
1247  {
1248  NFluctuations++;
1249  }
1250  }
1251 
1252  return NFluctuations;
1253 }
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:2799
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 1015 of file DPRawHitFinder_module.cc.

References fTicksToStopPeakFinder.

Referenced by produce().

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

References CreateFitFunction(), t1, and t2.

Referenced by produce().

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

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

Referenced by produce().

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

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

Referenced by produce().

1096 {
1097  // ################################################################
1098  // ### Lets loop over the candidate pulses we found in this ROI ###
1099  // ################################################################
1100  auto timeValsVecItr = timeValsVec.begin();
1101  unsigned int PeaksInThisMergedPeak = 0;
1102  //int EndTickOfPreviousMergedPeak=0;
1103  while(timeValsVecItr != timeValsVec.end())
1104  {
1105  PeakTimeWidVec peakVals;
1106 
1107  // Setting the start, peak, and end time of the pulse
1108  auto& timeVal = *timeValsVecItr++;
1109  int startT = std::get<0>(timeVal);
1110  int maxT = std::get<1>(timeVal);
1111  int endT = std::get<2>(timeVal);
1112  int widT = endT - startT;
1113  int FinalStartT=startT;
1114  int FinalEndT=endT;
1115 
1116  int NFluctuations=EstimateFluctuations(signalVec, startT, maxT, endT);
1117 
1118  peakVals.emplace_back(maxT,widT,startT,endT);
1119 
1120  // See if we want to merge pulses together
1121  // First check if we have more than one pulse on the wire
1122  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1123 
1124  // Loop until no more merged pulses (or candidates in this ROI)
1125  while(checkNextHit)
1126  {
1127  // 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
1128  int NextStartT = std::get<0>(*timeValsVecItr);
1129 
1130  double MinADC = signalVec[endT];
1131  for(int i = endT; i <= NextStartT; i++)
1132  {
1133  if(signalVec[i]<MinADC)
1134  {
1135  MinADC = signalVec[i];
1136  }
1137  }
1138 
1139  // Group peaks (grouped peaks are fitted together and can be merged)
1140  if( MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance )
1141  {
1142  int CurrentStartT=startT;
1143  int CurrentMaxT=maxT;
1144  int CurrentEndT=endT;
1145  //int CurrentWidT=widT;
1146 
1147  timeVal = *timeValsVecItr++;
1148  int NextMaxT = std::get<1>(timeVal);
1149  int NextEndT = std::get<2>(timeVal);
1150  int NextWidT = NextEndT - NextStartT;
1151  FinalEndT=NextEndT;
1152 
1153  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1154 
1155  int CurrentSumADC = 0;
1156  for(int i = CurrentStartT; i<= CurrentEndT; i++)
1157  {
1158  CurrentSumADC+=signalVec[i];
1159  }
1160 
1161  int NextSumADC = 0;
1162  for (int i = NextStartT; i<= NextEndT; i++)
1163  {
1164  NextSumADC+=signalVec[i];
1165  }
1166 
1167  //Merge peaks within a group
1168  if(fDoMergePeaks)
1169  {
1170  if( signalVec[NextMaxT] <= signalVec[CurrentMaxT] && ( ( signalVec[NextMaxT] < fMergeMaxADCThreshold*signalVec[CurrentMaxT] && NextSumADC < fMergeADCSumThreshold*CurrentSumADC ) || (signalVec[NextMaxT]-signalVec[NextStartT]) < fMinRelativePeakHeightRight*signalVec[CurrentMaxT] ) && ( signalVec[NextMaxT] <= fMergeMaxADCLimit ) )
1171  {
1172  maxT=CurrentMaxT;
1173  startT=CurrentStartT;
1174  endT=NextEndT;
1175  widT=endT - startT;
1176  peakVals.pop_back();
1177  peakVals.emplace_back(maxT,widT,startT,endT);
1178  }
1179  else if( signalVec[NextMaxT] > signalVec[CurrentMaxT] && ( ( signalVec[CurrentMaxT] < fMergeMaxADCThreshold*signalVec[NextMaxT] && CurrentSumADC < fMergeADCSumThreshold*NextSumADC ) || (signalVec[CurrentMaxT]-signalVec[CurrentEndT]) < fMinRelativePeakHeightLeft*signalVec[NextMaxT] ) && ( signalVec[CurrentMaxT] <= fMergeMaxADCLimit ) )
1180  {
1181  maxT=NextMaxT;
1182  startT=CurrentStartT;
1183  endT=NextEndT;
1184  widT=endT - startT;
1185  peakVals.pop_back();
1186  peakVals.emplace_back(maxT,widT,startT,endT);
1187  }
1188  else
1189  {
1190  maxT=NextMaxT;
1191  startT=NextStartT;
1192  endT=NextEndT;
1193  widT=NextWidT;
1194  peakVals.emplace_back(maxT,widT,startT,endT);
1195  PeaksInThisMergedPeak++;
1196  }
1197  }
1198  else
1199  {
1200  maxT=NextMaxT;
1201  startT=NextStartT;
1202  endT=NextEndT;
1203  widT=NextWidT;
1204  peakVals.emplace_back(maxT,widT,startT,endT);
1205  PeaksInThisMergedPeak++;
1206  }
1207  checkNextHit = timeValsVecItr != timeValsVec.end();
1208  }//<---Checking adjacent pulses
1209  else
1210  {
1211  checkNextHit = false;
1212  PeaksInThisMergedPeak = 0;
1213  }
1214 
1215  }//<---End checking if there is more than one pulse on the wire
1216  // Add these to our merged vector
1217  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1218  }
1219  return;
1220 }
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:1225
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 = " << chi2PerNDF << std::endl;
558 
559  if(fSameShape)
560  {
561  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
562  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
563 
564  for(unsigned int i = 0; i < nExponentialsForFit; i++)
565  {
566  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2*(i+1)].first << std::endl;
567  std::cout << "Peak #" << i << ": t0 [ticks] = " << 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] = " << paramVec[4*i+2].first << std::endl;
575  std::cout << "Peak #" << i << ": t0 [ticks] = " << range.offset + paramVec[4*i+3].first << std::endl;
576  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4*i].first << std::endl;
577  std::cout << "Peak #" << i << ": tau2 [mus] = " << 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 = " << chi2PerNDF << std::endl;
683 
684  if(fSameShape)
685  {
686  std::cout << "tau1 [mus] = " << paramVec[0].first << std::endl;
687  std::cout << "tau2 [mus] = " << paramVec[1].first << std::endl;
688 
689  for(unsigned int i = 0; i < nExponentialsForFit; i++)
690  {
691  std::cout << "Peak #" << i << ": A [ADC] = " << paramVec[2*(i+1)].first << std::endl;
692  std::cout << "Peak #" << i << ": t0 [ticks] = " << 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] = " << paramVec[4*i+2].first << std::endl;
700  std::cout << "Peak #" << i << ": t0 [ticks] = " << range.offset + paramVec[4*i+3].first << std::endl;
701  std::cout << "Peak #" << i << ": tau1 [mus] = " << paramVec[4*i].first << std::endl;
702  std::cout << "Peak #" << i << ": tau2 [mus] = " << 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 
755  // Extract fit parameter errors
756  double peakMeanErr;
757 
758  if(fSameShape)
759  {
760  peakMeanErr = paramVec[2*(i+1)+1].second;
761  }
762  else
763  {
764  peakMeanErr = paramVec[4*i+3].second;
765  }
766  double peakWidthErr = 0.1*peakWidth;
767 
768  // ### Charge ###
769  double charge = ChargeFunc(peakMean, peakAmp, peakTau1, peakTau2, fChargeNorm, peakMeanTrue);
770  double chargeErr = std::sqrt(TMath::Pi()) * (peakAmpErr*peakWidthErr + peakWidthErr*peakAmpErr);
771 
772  // ### limits for getting sum of ADC counts
773  int startTthisHit = std::get<2>(peakVals.at(i));
774  int endTthisHit = std::get<3>(peakVals.at(i));
775  std::vector<float>::const_iterator sumStartItr = signal.begin() + startTthisHit;
776  std::vector<float>::const_iterator sumEndItr = signal.begin() + endTthisHit;
777 
778  // ### Sum of ADC counts
779  double sumADC = std::accumulate(sumStartItr, sumEndItr + 1, 0.);
780 
781 
782  //Check if fit returns reasonable values and ich chi2 is below threshold
783  if(peakWidth <= 0 || charge <= 0. || charge != charge || ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) || ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) )
784  {
785  if(fLogLevel >= 1)
786  {
787  std::cout << std::endl;
788  std::cout << "WARNING: For peak #" << i << " in this group:" << std::endl;
789  if( peakWidth <= 0 || charge <= 0. || charge != charge) std::cout << "Fit function returned width < 0 or charge < 0 or charge = nan." << std::endl;
790  if( ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) || ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) )
791  {
792  std::cout << std::endl;
793  std::cout << "WARNING: For fit of this group (" << NumberOfPeaksBeforeFit << " peaks before refit, " << nExponentialsForFit << " peaks after refit): " << std::endl;
794  if ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMax << ")." << std::endl;
795  if ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMaxFactorMultiHits*fChi2NDFMax << ")." << std::endl;
796  }
797  std::cout << "---> DO NOT create hit object from fit parameters but use peak values instead." << std::endl;
798  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;
799  }
800  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
801  peakMeanErr=peakWidth/2;
802  charge = sumADC;
803  peakMeanTrue = std::get<0>(peakVals.at(i));
804 
805  //set the fit values to make it visible in the event display that this fit failed
806  peakMean = peakMeanTrue;
807  peakTau1 = 0.008;
808  peakTau2 = 0.0065;
809  peakAmp = 20.;
810  }
811 
812  // Create the hit
813  recob::HitCreator hitcreator(*wire, // wire reference
814  wid, // wire ID
815  startTthisHit+roiFirstBinTick, // start_tick TODO check
816  endTthisHit+roiFirstBinTick, // end_tick TODO check
817  peakWidth, // rms
818  peakMeanTrue+roiFirstBinTick, // peak_time
819  peakMeanErr, // sigma_peak_time
820  peakAmpTrue, // peak_amplitude
821  peakAmpErr, // sigma_peak_amplitude
822  charge, // hit_integral
823  chargeErr, // hit_sigma_integral
824  sumADC, // summedADC FIXME
825  nExponentialsForFit, // multiplicity
826  numHits, // local_index TODO check that the order is correct
827  chi2PerNDF, // goodness_of_fit
828  NDF // dof
829  );
830 
831  if(fLogLevel >=6)
832  {
833  std::cout << std::endl;
834  std::cout << "- Created hit object for peak #" << i << " in this group with the following parameters (obtained from fit):" << std::endl;
835  std::cout << "HitStartTick: " << startTthisHit+roiFirstBinTick << std::endl;
836  std::cout << "HitEndTick: " << endTthisHit+roiFirstBinTick << std::endl;
837  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
838  std::cout << "HitMeanTick: " << peakMeanTrue+roiFirstBinTick << " +- " << peakMeanErr << std::endl;
839  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr << std::endl;
840  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr << std::endl;
841  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
842  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
843  std::cout << "HitIndex in group: " << numHits << std::endl;
844  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
845  std::cout << "HitNDF: " << NDF << std::endl;
846  }
847 
848  const recob::Hit hit(hitcreator.move());
849 
850  hcol.emplace_back(std::move(hit), wire, rawdigits);
851  // add fit parameters associated to the hit just pushed to the collection
852  std::array<float, 4> fitParams;
853  fitParams[0] = peakMean+roiFirstBinTick;
854  fitParams[1] = peakTau1;
855  fitParams[2] = peakTau2;
856  fitParams[3] = peakAmp;
857  fHitParamWriter.addVector(hitID, fitParams);
858  numHits++;
859  } // <---End loop over Exponentials
860 // } // <---End if chi2 <= chi2Max
861  } // <---End if(NumberOfPeaksBeforeFit <= fMaxMultiHit && width <= fMaxGroupLength), then fit
862 
863  // #######################################################
864  // ### If too large then force alternate solution ###
865  // ### - Make n hits from pulse train where n will ###
866  // ### depend on the fhicl parameter fLongPulseWidth ###
867  // ### Also do this if chi^2 is too large ###
868  // #######################################################
869  if( NumberOfPeaksBeforeFit > fMaxMultiHit || (width > fMaxGroupLength) || NFluctuations > fMaxFluctuations)
870  {
871 
872  int nHitsInThisGroup = (endT - startT + 1) / fLongPulseWidth;
873 
874  if (nHitsInThisGroup > fLongMaxHits)
875  {
876  nHitsInThisGroup = fLongMaxHits;
877  fLongPulseWidth = (endT - startT + 1) / nHitsInThisGroup;
878  }
879 
880  if (nHitsInThisGroup * fLongPulseWidth < (endT - startT + 1) ) nHitsInThisGroup++;
881 
882  int firstTick = startT;
883  int lastTick = std::min(endT,firstTick+fLongPulseWidth-1);
884 
885  if(fLogLevel >= 1)
886  {
887  if( NumberOfPeaksBeforeFit > fMaxMultiHit)
888  {
889  std::cout << std::endl;
890  std::cout << "WARNING: Number of peaks in this group (" << NumberOfPeaksBeforeFit << ") is higher than threshold (" << fMaxMultiHit << ")." << std::endl;
891  std::cout << "---> DO NOT fit. Split group of peaks into hits with equal length instead." << std::endl;
892  }
893  if( width > fMaxGroupLength)
894  {
895  std::cout << std::endl;
896  std::cout << "WARNING: group of peak is longer (" << width << " ticks) than threshold (" << fMaxGroupLength << " ticks)." << std::endl;
897  std::cout << "---> DO NOT fit. Split group of peaks into hits with equal length instead." << std::endl;
898  }
899  if( NFluctuations > fMaxFluctuations)
900  {
901  std::cout << std::endl;
902  std::cout << "WARNING: fluctuations (" << NFluctuations << ") higher than threshold (" << fMaxFluctuations << ")." << std::endl;
903  std::cout << "---> DO NOT fit. Split group of peaks into hits with equal length instead." << std::endl;
904  }
905 /*
906  if( ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) || ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) )
907  {
908  std::cout << std::endl;
909  std::cout << "WARNING: For fit of this group (" << NumberOfPeaksBeforeFit << " peaks before refit, " << nExponentialsForFit << " peaks after refit): " << std::endl;
910  if ( nExponentialsForFit == 1 && chi2PerNDF > fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMax << ")." << std::endl;
911  if ( nExponentialsForFit >= 2 && chi2PerNDF > fChi2NDFMaxFactorMultiHits*fChi2NDFMax ) std::cout << "chi2/ndf of this fit (" << chi2PerNDF << ") is higher than threshold (" << fChi2NDFMaxFactorMultiHits*fChi2NDFMax << ")." << std::endl;
912  std::cout << "---> DO NOT create hit object but split group of peaks into hits with equal length instead." << std::endl;
913  }*/
914  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;
915  }
916 
917 
918  for(int hitIdx = 0; hitIdx < nHitsInThisGroup; hitIdx++)
919  {
920  // This hit parameters
921  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
922  double peakMeanTrue = (firstTick + lastTick) / 2.;
923  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
924  double peakMeanErr = (lastTick - firstTick) / 2.;
925  double sumADC = std::accumulate(signal.begin() + firstTick, signal.begin() + lastTick + 1, 0.);
926  double charge = sumADC;
927  double chargeErr = 0.1*sumADC;
928  double peakAmpTrue = 0;
929 
930  for(int tick = firstTick; tick <= lastTick; tick++)
931  {
932  if(signal[tick] > peakAmpTrue) peakAmpTrue = signal[tick];
933  }
934 
935  double peakAmpErr = 1.;
936  nExponentialsForFit = nHitsInThisGroup;
937  NDF = -1;
938  chi2PerNDF = -1.;
939  //set the fit values to make it visible in the event display that this fit failed
940  double peakMean = peakMeanTrue-2;
941  double peakTau1 = 0.008;
942  double peakTau2 = 0.0065;
943  double peakAmp = 20.;
944 
945  recob::HitCreator hitcreator(*wire, // wire reference
946  wid, // wire ID
947  firstTick+roiFirstBinTick, // start_tick TODO check
948  lastTick+roiFirstBinTick, // end_tick TODO check
949  peakWidth, // rms
950  peakMeanTrue+roiFirstBinTick, // peak_time
951  peakMeanErr, // sigma_peak_time
952  peakAmpTrue, // peak_amplitude
953  peakAmpErr, // sigma_peak_amplitude
954  charge, // hit_integral
955  chargeErr, // hit_sigma_integral
956  sumADC, // summedADC FIXME
957  nExponentialsForFit, // multiplicity
958  hitIdx, // local_index TODO check that the order is correct
959  chi2PerNDF, // goodness_of_fit
960  NDF // dof
961  );
962 
963 
964  if(fLogLevel >=6)
965  {
966  std::cout << std::endl;
967  std::cout << "- Created hit object for peak #" << hitIdx << " in this group with the following parameters (obtained from waveform):" << std::endl;
968  std::cout << "HitStartTick: " << firstTick+roiFirstBinTick << std::endl;
969  std::cout << "HitEndTick: " << lastTick+roiFirstBinTick << std::endl;
970  std::cout << "HitWidthTicks: " << peakWidth << std::endl;
971  std::cout << "HitMeanTick: " << peakMeanTrue+roiFirstBinTick << " +- " << peakMeanErr << std::endl;
972  std::cout << "HitAmplitude [ADC]: " << peakAmpTrue << " +- " << peakAmpErr << std::endl;
973  std::cout << "HitIntegral [ADC*ticks]: " << charge << " +- " << chargeErr << std::endl;
974  std::cout << "HitADCSum [ADC*ticks]: " << sumADC << std::endl;
975  std::cout << "HitMultiplicity: " << nExponentialsForFit << std::endl;
976  std::cout << "HitIndex in group: " << hitIdx << std::endl;
977  std::cout << "Hitchi2/ndf: " << chi2PerNDF << std::endl;
978  std::cout << "HitNDF: " << NDF << std::endl;
979  }
980  const recob::Hit hit(hitcreator.move());
981  hcol.emplace_back(std::move(hit), wire, rawdigits);
982 
983  std::array<float, 4> fitParams;
984  fitParams[0] = peakMean+roiFirstBinTick;
985  fitParams[1] = peakTau1;
986  fitParams[2] = peakTau2;
987  fitParams[3] = peakAmp;
988  fHitParamWriter.addVector(hitID, fitParams);
989 
990  // set for next loop
991  firstTick = lastTick+1;
992  lastTick = std::min(firstTick + fLongPulseWidth - 1, endT);
993 
994  }//<---Hits in this group
995  }//<---End if #peaks > MaxMultiHit
996  fChi2->Fill(chi2PerNDF);
997  }//<---End loop over merged candidate hits
998  } //<---End looping over ROI's
999  }//<---End looping over all the wires
1000 
1001  //==================================================================================================
1002  // End of the event
1003 
1004  // move the hit collection and the associations into the event
1005  hcol.put_into(evt);
1006 
1007  // and put hit fit parameters together with metadata into the event
1008  fHitParamWriter.saveOutputs(evt);
1009 
1010 } // 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:1015
hit::DPRawHitFinder::mergeCandidatePeaks
void mergeCandidatePeaks(const std::vector< float > signalVec, TimeValsVec, MergedTimeWidVec &)
Definition: DPRawHitFinder_module.cc:1095
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:1471
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:1168
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:690
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:1643
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:1599
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:1762
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:1258
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:1801
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:1681
void hit::DPRawHitFinder::reBin ( const std::vector< float > &  inputVec,
std::vector< float > &  outputVec,
size_t  nBinsToCombine 
) const
private

Definition at line 1831 of file DPRawHitFinder_module.cc.

References DEFINE_ART_MODULE.

1834 {
1835  size_t nNewBins = inputVec.size() / nBinsToCombine;
1836 
1837  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1838 
1839  outputVec.resize(nNewBins, 0.);
1840 
1841  size_t outputBin = 0;
1842 
1843  for(size_t inputIdx = 0; inputIdx < inputVec.size();)
1844  {
1845  outputVec[outputBin] += inputVec[inputIdx++];
1846 
1847  if (inputIdx % nBinsToCombine == 0) outputBin++;
1848 
1849  if (outputBin > outputVec.size())
1850  {
1851  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin << ", inputIdx = " << inputIdx << std::endl;
1852  break;
1853  }
1854  }
1855 
1856  return;
1857 }
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 1643 of file DPRawHitFinder_module.cc.

References t1, and t2.

Referenced by produce().

1645 {
1646 int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1647 int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1648 int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1649 int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1650 
1651 if(WidthOldPeakOld<3) {return;} //can't split peaks with widths < 3
1652 
1653 int NewPeakMax = 0;
1654 int NewPeakStart = 0;
1655 int NewPeakEnd = 0;
1656 int OldPeakNewMax = 0;
1657 int OldPeakNewStart = 0;
1658 int OldPeakNewEnd = 0;
1659 
1660 
1661 OldPeakNewStart = OldPeakOldStart;
1662 NewPeakEnd = OldPeakOldEnd;
1663 
1664 OldPeakNewEnd = OldPeakNewStart + 0.5*(WidthOldPeakOld + (WidthOldPeakOld%2));
1665 NewPeakStart = OldPeakNewEnd+1;
1666 
1667 int WidthOldPeakNew = OldPeakNewEnd-OldPeakNewStart;
1668 int WidthNewPeak = NewPeakEnd-NewPeakStart;
1669 
1670 OldPeakNewMax = OldPeakNewStart + 0.5*(WidthOldPeakNew - (WidthOldPeakNew%2));
1671 NewPeakMax = NewPeakStart + 0.5*(WidthNewPeak - (WidthNewPeak%2));
1672 
1673 fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakNewMax,0,OldPeakNewStart,OldPeakNewEnd);
1674 fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1675 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);} );
1676 
1677 return;
1678 }
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 1681 of file DPRawHitFinder_module.cc.

References x.

Referenced by produce().

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