Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
hit::DPRawHitFinder Class Reference
Inheritance diagram for hit::DPRawHitFinder:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Classes

struct  Comp
 

Public Member Functions

 DPRawHitFinder (fhicl::ParameterSet const &pset)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Producer & operator= (Producer const &)=delete
 
Producer & operator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Types

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

Private Member Functions

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

Private Attributes

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

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 69 of file DPRawHitFinder_module.cc.

Member Typedef Documentation

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

Definition at line 82 of file DPRawHitFinder_module.cc.

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

Definition at line 84 of file DPRawHitFinder_module.cc.

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

Definition at line 83 of file DPRawHitFinder_module.cc.

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

Definition at line 81 of file DPRawHitFinder_module.cc.

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

Definition at line 80 of file DPRawHitFinder_module.cc.

Constructor & Destructor Documentation

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

Definition at line 211 of file DPRawHitFinder_module.cc.

211  :
212  EDProducer{pset},
213  fNewHitsTag(
214  pset.get<std::string>("module_label"), "",
215  art::ServiceHandle<art::TriggerNamesService const>()->getProcessName()),
216  fHitParamWriter(producesCollector())
217 {
218  fLogLevel = pset.get< int >("LogLevel");
219  fCalDataModuleLabel = pset.get< std::string >("CalDataModuleLabel");
220  fMaxMultiHit = pset.get< int >("MaxMultiHit");
221  fMaxGroupLength = pset.get< int >("MaxGroupLength");
222  fTryNplus1Fits = pset.get< bool >("TryNplus1Fits");
223  fChi2NDFRetry = pset.get< double >("Chi2NDFRetry");
224  fChi2NDFRetryFactorMultiHits = pset.get< double >("Chi2NDFRetryFactorMultiHits");
225  fChi2NDFMax = pset.get< double >("Chi2NDFMax");
226  fChi2NDFMaxFactorMultiHits = pset.get< double >("Chi2NDFMaxFactorMultiHits");
227  fNumBinsToAverage = pset.get< size_t >("NumBinsToAverage");
228  fMinSig = pset.get< float >("MinSig");
229  fMinWidth = pset.get< int >("MinWidth");
230  fMinADCSum = pset.get< double >("MinADCSum");
231  fMinADCSumOverWidth = pset.get< double >("MinADCSumOverWidth");
232  fChargeNorm = pset.get< double >("ChargeNorm");
233  fTicksToStopPeakFinder = pset.get< double >("TicksToStopPeakFinder");
234  fMinTau = pset.get< double >("MinTau");
235  fMaxTau = pset.get< double >("MaxTau");
236  fFitPeakMeanRange = pset.get< double >("FitPeakMeanRange");
237  fGroupMaxDistance = pset.get< int >("GroupMaxDistance");
238  fGroupMinADC = pset.get< int >("GroupMinADC");
239  fSameShape = pset.get< bool >("SameShape");
240  fDoMergePeaks = pset.get< bool >("DoMergePeaks");
241  fMergeADCSumThreshold = pset.get< double >("MergeADCSumThreshold");
242  fMergeMaxADCThreshold = pset.get< double >("MergeMaxADCThreshold");
243  fMinRelativePeakHeightLeft = pset.get< double >("MinRelativePeakHeightLeft");
244  fMinRelativePeakHeightRight = pset.get< double >("MinRelativePeakHeightRight");
245  fMergeMaxADCLimit = pset.get< double >("MergeMaxADCLimit");
246  fWidthNormalization = pset.get< double >("WidthNormalization");
247  fLongMaxHits = pset.get< double >("LongMaxHits");
248  fLongPulseWidth = pset.get< double >("LongPulseWidth");
249  fMaxFluctuations = pset.get< double >("MaxFluctuations");
250 
251  // let HitCollectionCreator declare that we are going to produce
252  // hits and associations with wires and raw digits
253  // (with no particular product label)
254  recob::HitCollectionCreator::declare_products(producesCollector());
255 
256  // declare that data products with feature vectors describing
257  // hits is going to be produced
258  fHitParamWriter.produces_using< recob::Hit >();
259 
260 } // DPRawHitFinder::DPRawHitFinder()
string
std::string string
Definition: nybbler.cc:12
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
recob::HitAndAssociationsWriterBase::declare_products
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
hit::DPRawHitFinder::fHitParamWriter
anab::FVectorWriter< 4 > fHitParamWriter
Definition: DPRawHitFinder_module.cc:201
art::ProductRegistryHelper::producesCollector
ProducesCollector & producesCollector() noexcept
Definition: ProductRegistryHelper.h:80
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48

Member Function Documentation

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

Definition at line 1565 of file DPRawHitFinder_module.cc.

1567 {
1568  int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1569  int NewPeakMax = std::get<2>(fPeakDevCand);
1570  int OldPeakMax = std::get<0>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1571  int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1572  int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1573 
1574  int NewPeakStart=0;
1575  int NewPeakEnd=0;
1576  int OldPeakNewStart=0;
1577  int OldPeakNewEnd=0;
1578  int DistanceBwOldAndNewPeak=0;
1579 
1580  if(NewPeakMax<OldPeakMax)
1581  {
1582  NewPeakStart = OldPeakOldStart;
1583  OldPeakNewEnd = OldPeakOldEnd;
1584  DistanceBwOldAndNewPeak = OldPeakMax - NewPeakMax;
1585  NewPeakEnd = NewPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1586  if(DistanceBwOldAndNewPeak%2 == 0) NewPeakEnd -= 1;
1587  OldPeakNewStart = NewPeakEnd+1;
1588  }
1589  else if(OldPeakMax<NewPeakMax)
1590  {
1591  NewPeakEnd = OldPeakOldEnd;
1592  OldPeakNewStart = OldPeakOldStart;
1593  DistanceBwOldAndNewPeak = NewPeakMax - OldPeakMax;
1594  OldPeakNewEnd = OldPeakMax+0.5*(DistanceBwOldAndNewPeak-(DistanceBwOldAndNewPeak%2));
1595  if(DistanceBwOldAndNewPeak%2 == 0) OldPeakNewEnd -= 1;
1596  NewPeakStart = OldPeakNewEnd+1;
1597  }
1598  else if(OldPeakMax==NewPeakMax){return;} //This shouldn't happen
1599 
1600  fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakMax,0,OldPeakNewStart,OldPeakNewEnd);
1601  fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1602  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);} );
1603 
1604  return;
1605 }
test_gpu_visible.t1
t1
Definition: test_gpu_visible.py:16
void hit::DPRawHitFinder::beginJob ( )
overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 286 of file DPRawHitFinder_module.cc.

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

Definition at line 1728 of file DPRawHitFinder_module.cc.

1735 {
1736 double ChargeSum = 0.;
1737 double Charge = 0.;
1738 double ReBin = 10.;
1739 
1740 bool ChargeBigEnough=true;
1741  for(double x = fPeakMeanTrue - 1/ReBin; ChargeBigEnough && x > fPeakMeanTrue-1000.; x-=1.)
1742  {
1743  for(double i=0.; i > -1.; i-=(1/ReBin))
1744  {
1745  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1746  ChargeSum += Charge;
1747  }
1748  if(Charge < 0.01) ChargeBigEnough = false;
1749  }
1750 
1751 ChargeBigEnough=true;
1752  for(double x = fPeakMeanTrue; ChargeBigEnough && x < fPeakMeanTrue+1000.; x+=1.)
1753  {
1754  for(double i=0.; i < 1.; i+=(1/ReBin))
1755  {
1756  Charge = ( fPeakAmp * exp(0.4*(x+i-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x+i-fPeakMean)/fPeakTau2) );
1757  ChargeSum += Charge;
1758  }
1759  if(Charge < 0.01) ChargeBigEnough = false;
1760  }
1761 
1762 
1763 return ChargeSum*fChargeNormFactor/ReBin;
1764 }
train.x
list x
Definition: train.py:276
std::string hit::DPRawHitFinder::CreateFitFunction ( int  fNPeaks,
bool  fSameShape 
)
private

Definition at line 1501 of file DPRawHitFinder_module.cc.

1502 {
1503  std::string feqn = ""; // string holding fit formula
1504  std::stringstream numConv;
1505 
1506  if(fSameShape)
1507  {
1508  for(int i = 0; i < fNPeaks; i++)
1509  {
1510  feqn.append("+( [");
1511  numConv.str("");
1512  numConv << 2*(i+1);
1513  feqn.append(numConv.str());
1514  feqn.append("] * exp(0.4*(x-[");
1515  numConv.str("");
1516  numConv << 2*(i+1)+1;
1517  feqn.append(numConv.str());
1518  feqn.append("])/[");
1519  numConv.str("");
1520  numConv << 0;
1521  feqn.append(numConv.str());
1522  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1523  numConv.str("");
1524  numConv << 2*(i+1)+1;
1525  feqn.append(numConv.str());
1526  feqn.append("])/[");
1527  numConv.str("");
1528  numConv << 1;
1529  feqn.append(numConv.str());
1530  feqn.append("]) ) )");
1531  }
1532  }
1533  else
1534  {
1535  for(int i = 0; i < fNPeaks; i++)
1536  {
1537  feqn.append("+( [");
1538  numConv.str("");
1539  numConv << 4*i+2;
1540  feqn.append(numConv.str());
1541  feqn.append("] * exp(0.4*(x-[");
1542  numConv.str("");
1543  numConv << 4*i+3;
1544  feqn.append(numConv.str());
1545  feqn.append("])/[");
1546  numConv.str("");
1547  numConv << 4*i;
1548  feqn.append(numConv.str());
1549  feqn.append("]) / ( 1 + exp(0.4*(x-[");
1550  numConv.str("");
1551  numConv << 2*(i+1)+1;
1552  feqn.append(numConv.str());
1553  feqn.append("])/[");
1554  numConv.str("");
1555  numConv << 4*i+1;
1556  feqn.append(numConv.str());
1557  feqn.append("]) ) )");
1558  }
1559  }
1560 return feqn;
1561 }
string
std::string string
Definition: nybbler.cc:12
void hit::DPRawHitFinder::doBinAverage ( const std::vector< float > &  inputVec,
std::vector< float > &  outputVec,
size_t  binsToAverage 
) const
private

Definition at line 1767 of file DPRawHitFinder_module.cc.

1770 {
1771  size_t halfBinsToAverage(binsToAverage/2);
1772 
1773  float runningSum(0.);
1774 
1775  for(size_t idx = 0; idx < halfBinsToAverage; idx++) runningSum += inputVec[idx];
1776 
1777  outputVec.resize(inputVec.size());
1778  std::vector<float>::iterator outputVecItr = outputVec.begin();
1779 
1780  // First pass through to build the erosion vector
1781  for(std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end(); inputItr++)
1782  {
1783  size_t startOffset = std::distance(inputVec.begin(),inputItr);
1784  size_t stopOffset = std::distance(inputItr,inputVec.end());
1785  size_t count = std::min(2 * halfBinsToAverage, std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1786 
1787  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1788  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1789 
1790  *outputVecItr++ = runningSum / float(count);
1791  }
1792 
1793  return;
1794 }
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
reco_momentum_tuples.count
int count
Definition: reco_momentum_tuples.py:202
int hit::DPRawHitFinder::EstimateFluctuations ( const std::vector< float >  fsignalVec,
int  peakStart,
int  peakMean,
int  peakEnd 
)
private

Definition at line 1191 of file DPRawHitFinder_module.cc.

1195 {
1196  int NFluctuations=0;
1197 
1198  for(int j = peakMean-1; j >= peakStart; j--)
1199  {
1200  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1201 
1202  if(fsignalVec[j] > fsignalVec[j+1])
1203  {
1204  NFluctuations++;
1205  }
1206  }
1207 
1208  for(int j = peakMean+1; j <= peakEnd; j++)
1209  {
1210  if(fsignalVec[j] < 5) break; //do not look for fluctuations below 5 ADC counts (this should depend on the noise RMS)
1211 
1212  if(fsignalVec[j] > fsignalVec[j-1])
1213  {
1214  NFluctuations++;
1215  }
1216  }
1217 
1218  return NFluctuations;
1219 }
void hit::DPRawHitFinder::FillOutHitParameterVector ( const std::vector< double > &  input,
std::vector< double > &  output 
)
private

Definition at line 265 of file DPRawHitFinder_module.cc.

267 {
268  if(input.size()==0)
269  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
270 
271  art::ServiceHandle<geo::Geometry const> geom;
272  const unsigned int N_PLANES = geom->Nplanes();
273 
274  if(input.size()==1)
275  output.resize(N_PLANES,input[0]);
276  else if(input.size()==N_PLANES)
277  output = input;
278  else
279  throw std::runtime_error("DPRawHitFinder::FillOutHitParameterVector ERROR! Input config vector size !=1 and !=N_PLANES.");
280 
281 }
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:2848
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 981 of file DPRawHitFinder_module.cc.

986 {
987  // Need a minimum number of ticks to do any work here
988  if (std::distance(startItr,stopItr) > 4)
989  {
990  // Find the highest peak in the range given
991  auto maxItr = std::max_element(startItr, stopItr);
992 
993  float maxValue = *maxItr;
994  int maxTime = std::distance(startItr,maxItr);
995 
996  if (maxValue >= PeakMin)
997  {
998  // backwards to find first bin for this candidate hit
999  auto firstItr = std::distance(startItr,maxItr) > 2 ? maxItr - 1 : startItr;
1000  bool PeakStartIsHere = true;
1001 
1002  while(firstItr != startItr)
1003  {
1004  //Check for inflection point & ADC count <= 0
1005  PeakStartIsHere = true;
1006  for(int i=1; i <= fTicksToStopPeakFinder; i++)
1007  {
1008  if( *firstItr >= *(firstItr-i) )
1009  {
1010  PeakStartIsHere = false;
1011  break;
1012  }
1013 
1014  }
1015  if (*firstItr <= 0 || PeakStartIsHere) break;
1016  firstItr--;
1017  }
1018 
1019  int firstTime = std::distance(startItr,firstItr);
1020 
1021  // Recursive call to find all candidate hits earlier than this peak
1022  findCandidatePeaks(startItr, firstItr - 1, timeValsVec, PeakMin, firstTick);
1023 
1024  // forwards to find last bin for this candidate hit
1025  auto lastItr = std::distance(maxItr,stopItr) > 2 ? maxItr + 1 : stopItr - 1;
1026  bool PeakEndIsHere = true;
1027 
1028  while(lastItr != stopItr)
1029  {
1030  //Check for inflection point & ADC count <= 0
1031  PeakEndIsHere = true;
1032  for(int i=1; i <= fTicksToStopPeakFinder; i++)
1033  {
1034  if( *lastItr >= *(lastItr+i) )
1035  {
1036  PeakEndIsHere = false;
1037  break;
1038  }
1039  }
1040  if (*lastItr <= 0 || PeakEndIsHere) break;
1041  lastItr++;
1042  }
1043 
1044  int lastTime = std::distance(startItr,lastItr);
1045 
1046  // Now save this candidate's start and max time info
1047  timeValsVec.push_back(std::make_tuple(firstTick+firstTime,firstTick+maxTime,firstTick+lastTime));
1048 
1049  // Recursive call to find all candidate hits later than this peak
1050  findCandidatePeaks(lastItr + 1, stopItr, timeValsVec, PeakMin, firstTick + std::distance(startItr,lastItr + 1));
1051  }
1052  }
1053 
1054  return;
1055 }
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:981
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
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 1437 of file DPRawHitFinder_module.cc.

1445 {
1446 // int size = fEndTime - fStartTime + 1;
1447 // if(fEndTime - fStartTime < 0){size = 0;}
1448 
1449  std::string eqn = CreateFitFunction(fNPeaks, fSameShape); // string for equation of Exponentials fit
1450 
1451  TF1 Exponentials("Exponentials",eqn.c_str(),fStartTime,fEndTime+1);
1452 
1453  for(size_t i=0; i < fparamVec.size(); i++)
1454  {
1455  Exponentials.SetParameter(i, fparamVec[i].first);
1456  }
1457 
1458  // ##########################################################################
1459  // ### Finding the peak with the max chi2 fit and signal ###
1460  // ##########################################################################
1461  double Chi2PerNDFPeak;
1462  double MaxPosDeviation;
1463  double MaxNegDeviation;
1464  int BinMaxPosDeviation;
1465  int BinMaxNegDeviation;
1466  for(int i = 0; i < fNPeaks; i++)
1467  {
1468  Chi2PerNDFPeak = 0.;
1469  MaxPosDeviation=0.;
1470  MaxNegDeviation=0.;
1471  BinMaxPosDeviation=0;
1472  BinMaxNegDeviation=0;
1473 
1474  for(int j = std::get<2>(fpeakVals.at(i)); j < std::get<3>(fpeakVals.at(i))+1; j++)
1475  {
1476  if( (Exponentials(j+0.5)-fSignalVector[j]) > MaxPosDeviation && j != std::get<0>(fpeakVals.at(i)) )
1477  {
1478  MaxPosDeviation = Exponentials(j+0.5)-fSignalVector[j];
1479  BinMaxPosDeviation = j;
1480  }
1481  if( (Exponentials(j+0.5)-fSignalVector[j]) < MaxNegDeviation && j != std::get<0>(fpeakVals.at(i)) )
1482  {
1483  MaxNegDeviation = Exponentials(j+0.5)-fSignalVector[j];
1484  BinMaxNegDeviation = j;
1485  }
1486  Chi2PerNDFPeak += pow((Exponentials(j+0.5)-fSignalVector[j])/sqrt(fSignalVector[j]),2);
1487  }
1488 
1489  if(BinMaxNegDeviation != 0)
1490  {
1491  Chi2PerNDFPeak /= static_cast<double>((std::get<3>(fpeakVals.at(i))-std::get<2>(fpeakVals.at(i))));
1492  fPeakDev.emplace_back(Chi2PerNDFPeak,i,BinMaxNegDeviation,BinMaxPosDeviation);
1493  }
1494  }
1495 
1496 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);} );
1497 Exponentials.Delete();
1498 }
string
std::string string
Definition: nybbler.cc:12
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
test_gpu_visible.t1
t1
Definition: test_gpu_visible.py:16
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1501
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 1224 of file DPRawHitFinder_module.cc.

1232 {
1233  int size = fEndTime - fStartTime + 1;
1234  int NPeaks = fPeakVals.size();
1235 
1236  // #############################################
1237  // ### If size < 0 then set the size to zero ###
1238  // #############################################
1239  if(fEndTime - fStartTime < 0){size = 0;}
1240 
1241  // --- TH1D HitSignal ---
1242  TH1F hitSignal("hitSignal","",std::max(size,1),fStartTime,fEndTime+1);
1243  hitSignal.Sumw2();
1244 
1245  // #############################
1246  // ### Filling the histogram ###
1247  // #############################
1248  for(int i = fStartTime; i < fEndTime+1; i++)
1249  {
1250  hitSignal.Fill(i,fSignalVector[i]);
1251  hitSignal.SetBinError(i,0.288675); //1/sqrt(12)
1252  }
1253 
1254  // ################################################
1255  // ### Building TFormula for basic Exponentials ###
1256  // ################################################
1257  std::string eqn = CreateFitFunction(NPeaks, fSameShape);
1258 
1259  // -------------------------------------
1260  // --- TF1 function for Exponentials ---
1261  // -------------------------------------
1262 
1263  TF1 Exponentials("Exponentials",eqn.c_str(),fStartTime,fEndTime+1);
1264 
1265  if(fLogLevel >= 4)
1266  {
1267  std::cout << std::endl;
1268  std::cout << "--- Preparing fit ---" << std::endl;
1269  std::cout << "--- Lower limits, seed, upper limit:" << std::endl;
1270  }
1271 
1272  if(fSameShape)
1273  {
1274  Exponentials.SetParameter(0, 0.5);
1275  Exponentials.SetParameter(1, 0.5);
1276  Exponentials.SetParLimits(0, fMinTau, fMaxTau);
1277  Exponentials.SetParLimits(1, fMinTau, fMaxTau);
1278  double amplitude=0;
1279  double peakMean=0;
1280 
1281  double peakMeanShift=2;
1282  double peakMeanSeed=0;
1283  double peakMeanRangeLow=0;
1284  double peakMeanRangeHi=0;
1285  double peakStart=0;
1286  double peakEnd=0;
1287 
1288  for(int i = 0; i < NPeaks; i++)
1289  {
1290  peakMean = std::get<0>(fPeakVals.at(i));
1291  peakStart = std::get<2>(fPeakVals.at(i));
1292  peakEnd = std::get<3>(fPeakVals.at(i));
1293  peakMeanSeed=peakMean-peakMeanShift;
1294  peakMeanRangeLow = std::max(peakStart-peakMeanShift, peakMeanSeed-fFitPeakMeanRange);
1295  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed+fFitPeakMeanRange);
1296  amplitude = fSignalVector[peakMean];
1297 
1298  Exponentials.SetParameter(2*(i+1), 1.65*amplitude);
1299  Exponentials.SetParLimits(2*(i+1), 0.3*1.65*amplitude, 2*1.65*amplitude);
1300  Exponentials.SetParameter(2*(i+1)+1, peakMeanSeed);
1301 
1302  if(NPeaks == 1)
1303  {
1304  Exponentials.SetParLimits(2*(i+1)+1, peakMeanRangeLow, peakMeanRangeHi);
1305  }
1306  else if(NPeaks >= 2 && i == 0)
1307  {
1308  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1309  Exponentials.SetParLimits( 2*(i+1)+1, peakMeanRangeLow, std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1310  }
1311  else if(NPeaks >= 2 && i == NPeaks-1)
1312  {
1313  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1314  Exponentials.SetParLimits(2*(i+1)+1, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), peakMeanRangeHi );
1315  }
1316  else
1317  {
1318  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1319  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1320  Exponentials.SetParLimits(2*(i+1)+1, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1321  }
1322 
1323  if(fLogLevel >= 4)
1324  {
1325  double t0low, t0high;
1326  Exponentials.GetParLimits(2*(i+1)+1, t0low, t0high);
1327  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3*1.65*amplitude << " , " << 1.65*amplitude << " , " << 2*1.65*amplitude << std::endl;
1328  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed << " , " << t0high << std::endl;
1329  }
1330  }
1331  }
1332  else
1333  {
1334  double amplitude=0;
1335  double peakMean=0;
1336 
1337  double peakMeanShift=2;
1338  double peakMeanSeed=0;
1339  double peakMeanRangeLow=0;
1340  double peakMeanRangeHi=0;
1341  double peakStart=0;
1342  double peakEnd=0;
1343 
1344 
1345  for(int i = 0; i < NPeaks; i++)
1346  {
1347  Exponentials.SetParameter(4*i, 0.5);
1348  Exponentials.SetParameter(4*i+1, 0.5);
1349  Exponentials.SetParLimits(4*i, fMinTau, fMaxTau);
1350  Exponentials.SetParLimits(4*i+1, fMinTau, fMaxTau);
1351 
1352  peakMean = std::get<0>(fPeakVals.at(i));
1353  peakStart = std::get<2>(fPeakVals.at(i));
1354  peakEnd = std::get<3>(fPeakVals.at(i));
1355  peakMeanSeed=peakMean-peakMeanShift;
1356  peakMeanRangeLow = std::max(peakStart-peakMeanShift, peakMeanSeed-fFitPeakMeanRange);
1357  peakMeanRangeHi = std::min(peakEnd, peakMeanSeed+fFitPeakMeanRange);
1358  amplitude = fSignalVector[peakMean];
1359 
1360  Exponentials.SetParameter(4*i+2, 1.65*amplitude);
1361  Exponentials.SetParLimits(4*i+2, 0.3*1.65*amplitude, 2*1.65*amplitude);
1362  Exponentials.SetParameter(4*i+3, peakMeanSeed);
1363 
1364  if(NPeaks == 1)
1365  {
1366  Exponentials.SetParLimits(4*i+3, peakMeanRangeLow, peakMeanRangeHi);
1367  }
1368  else if(NPeaks >= 2 && i == 0)
1369  {
1370  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1371  Exponentials.SetParLimits( 4*i+3, peakMeanRangeLow, std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1372  }
1373  else if(NPeaks >= 2 && i == NPeaks-1)
1374  {
1375  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1376  Exponentials.SetParLimits(4*i+3, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), peakMeanRangeHi );
1377  }
1378  else
1379  {
1380  double HalfDistanceToNextMean = 0.5*(std::get<0>(fPeakVals.at(i+1)) - peakMean);
1381  double HalfDistanceToPrevMean = 0.5*(peakMean - std::get<0>(fPeakVals.at(i-1)));
1382  Exponentials.SetParLimits(4*i+3, std::max(peakMeanRangeLow, peakMeanSeed-HalfDistanceToPrevMean), std::min(peakMeanRangeHi, peakMeanSeed+HalfDistanceToNextMean) );
1383  }
1384 
1385  if(fLogLevel >= 4)
1386  {
1387  double t0low, t0high;
1388  Exponentials.GetParLimits(4*i+3, t0low, t0high);
1389  std::cout << "Peak #" << i << ": A [ADC] = " << 0.3*1.65*amplitude << " , " << 1.65*amplitude << " , " << 2*1.65*amplitude << std::endl;
1390  std::cout << "Peak #" << i << ": t0 [ticks] = " << t0low << " , " << peakMeanSeed << " , " << t0high << std::endl;
1391  }
1392  }
1393  }
1394 
1395 
1396  // ###########################################
1397  // ### PERFORMING THE TOTAL FIT OF THE HIT ###
1398  // ###########################################
1399  try
1400  { hitSignal.Fit(&Exponentials,"QNRWM","", fStartTime, fEndTime+1);}
1401  catch(...)
1402  {mf::LogWarning("DPRawHitFinder") << "Fitter failed finding a hit";}
1403 
1404  // ##################################################
1405  // ### Getting the fitted parameters from the fit ###
1406  // ##################################################
1407  fchi2PerNDF = (Exponentials.GetChisquare() / Exponentials.GetNDF());
1408  fNDF = Exponentials.GetNDF();
1409 
1410  if(fSameShape)
1411  {
1412  fparamVec.emplace_back(Exponentials.GetParameter(0),Exponentials.GetParError(0));
1413  fparamVec.emplace_back(Exponentials.GetParameter(1),Exponentials.GetParError(1));
1414 
1415  for(int i = 0; i < NPeaks; i++)
1416  {
1417  fparamVec.emplace_back(Exponentials.GetParameter(2*(i+1)),Exponentials.GetParError(2*(i+1)));
1418  fparamVec.emplace_back(Exponentials.GetParameter(2*(i+1)+1),Exponentials.GetParError(2*(i+1)+1));
1419  }
1420  }
1421  else
1422  {
1423  for(int i = 0; i < NPeaks; i++)
1424  {
1425  fparamVec.emplace_back(Exponentials.GetParameter(4*i),Exponentials.GetParError(4*i));
1426  fparamVec.emplace_back(Exponentials.GetParameter(4*i+1),Exponentials.GetParError(4*i+1));
1427  fparamVec.emplace_back(Exponentials.GetParameter(4*i+2),Exponentials.GetParError(4*i+2));
1428  fparamVec.emplace_back(Exponentials.GetParameter(4*i+3),Exponentials.GetParError(4*i+3));
1429  }
1430  }
1431  Exponentials.Delete();
1432  hitSignal.Delete();
1433 }//<----End FitExponentials
string
std::string string
Definition: nybbler.cc:12
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
hit::DPRawHitFinder::CreateFitFunction
std::string CreateFitFunction(int fNPeaks, bool fSameShape)
Definition: DPRawHitFinder_module.cc:1501
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void hit::DPRawHitFinder::mergeCandidatePeaks ( const std::vector< float >  signalVec,
TimeValsVec  timeValsVec,
MergedTimeWidVec &  mergedVec 
)
private

Definition at line 1061 of file DPRawHitFinder_module.cc.

1062 {
1063  // ################################################################
1064  // ### Lets loop over the candidate pulses we found in this ROI ###
1065  // ################################################################
1066  auto timeValsVecItr = timeValsVec.begin();
1067  unsigned int PeaksInThisMergedPeak = 0;
1068  //int EndTickOfPreviousMergedPeak=0;
1069  while(timeValsVecItr != timeValsVec.end())
1070  {
1071  PeakTimeWidVec peakVals;
1072 
1073  // Setting the start, peak, and end time of the pulse
1074  auto& timeVal = *timeValsVecItr++;
1075  int startT = std::get<0>(timeVal);
1076  int maxT = std::get<1>(timeVal);
1077  int endT = std::get<2>(timeVal);
1078  int widT = endT - startT;
1079  int FinalStartT=startT;
1080  int FinalEndT=endT;
1081 
1082  int NFluctuations=EstimateFluctuations(signalVec, startT, maxT, endT);
1083 
1084  peakVals.emplace_back(maxT,widT,startT,endT);
1085 
1086  // See if we want to merge pulses together
1087  // First check if we have more than one pulse on the wire
1088  bool checkNextHit = timeValsVecItr != timeValsVec.end();
1089 
1090  // Loop until no more merged pulses (or candidates in this ROI)
1091  while(checkNextHit)
1092  {
1093  // 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
1094  int NextStartT = std::get<0>(*timeValsVecItr);
1095 
1096  double MinADC = signalVec[endT];
1097  for(int i = endT; i <= NextStartT; i++)
1098  {
1099  if(signalVec[i]<MinADC)
1100  {
1101  MinADC = signalVec[i];
1102  }
1103  }
1104 
1105  // Group peaks (grouped peaks are fitted together and can be merged)
1106  if( MinADC >= fGroupMinADC && NextStartT - endT <= fGroupMaxDistance )
1107  {
1108  int CurrentStartT=startT;
1109  int CurrentMaxT=maxT;
1110  int CurrentEndT=endT;
1111  //int CurrentWidT=widT;
1112 
1113  timeVal = *timeValsVecItr++;
1114  int NextMaxT = std::get<1>(timeVal);
1115  int NextEndT = std::get<2>(timeVal);
1116  int NextWidT = NextEndT - NextStartT;
1117  FinalEndT=NextEndT;
1118 
1119  NFluctuations += EstimateFluctuations(signalVec, NextStartT, NextMaxT, NextEndT);
1120 
1121  int CurrentSumADC = 0;
1122  for(int i = CurrentStartT; i<= CurrentEndT; i++)
1123  {
1124  CurrentSumADC+=signalVec[i];
1125  }
1126 
1127  int NextSumADC = 0;
1128  for (int i = NextStartT; i<= NextEndT; i++)
1129  {
1130  NextSumADC+=signalVec[i];
1131  }
1132 
1133  //Merge peaks within a group
1134  if(fDoMergePeaks)
1135  {
1136  if( signalVec[NextMaxT] <= signalVec[CurrentMaxT] && ( ( signalVec[NextMaxT] < fMergeMaxADCThreshold*signalVec[CurrentMaxT] && NextSumADC < fMergeADCSumThreshold*CurrentSumADC ) || (signalVec[NextMaxT]-signalVec[NextStartT]) < fMinRelativePeakHeightRight*signalVec[CurrentMaxT] ) && ( signalVec[NextMaxT] <= fMergeMaxADCLimit ) )
1137  {
1138  maxT=CurrentMaxT;
1139  startT=CurrentStartT;
1140  endT=NextEndT;
1141  widT=endT - startT;
1142  peakVals.pop_back();
1143  peakVals.emplace_back(maxT,widT,startT,endT);
1144  }
1145  else if( signalVec[NextMaxT] > signalVec[CurrentMaxT] && ( ( signalVec[CurrentMaxT] < fMergeMaxADCThreshold*signalVec[NextMaxT] && CurrentSumADC < fMergeADCSumThreshold*NextSumADC ) || (signalVec[CurrentMaxT]-signalVec[CurrentEndT]) < fMinRelativePeakHeightLeft*signalVec[NextMaxT] ) && ( signalVec[CurrentMaxT] <= fMergeMaxADCLimit ) )
1146  {
1147  maxT=NextMaxT;
1148  startT=CurrentStartT;
1149  endT=NextEndT;
1150  widT=endT - startT;
1151  peakVals.pop_back();
1152  peakVals.emplace_back(maxT,widT,startT,endT);
1153  }
1154  else
1155  {
1156  maxT=NextMaxT;
1157  startT=NextStartT;
1158  endT=NextEndT;
1159  widT=NextWidT;
1160  peakVals.emplace_back(maxT,widT,startT,endT);
1161  PeaksInThisMergedPeak++;
1162  }
1163  }
1164  else
1165  {
1166  maxT=NextMaxT;
1167  startT=NextStartT;
1168  endT=NextEndT;
1169  widT=NextWidT;
1170  peakVals.emplace_back(maxT,widT,startT,endT);
1171  PeaksInThisMergedPeak++;
1172  }
1173  checkNextHit = timeValsVecItr != timeValsVec.end();
1174  }//<---Checking adjacent pulses
1175  else
1176  {
1177  checkNextHit = false;
1178  PeaksInThisMergedPeak = 0;
1179  }
1180 
1181  }//<---End checking if there is more than one pulse on the wire
1182  // Add these to our merged vector
1183  mergedVec.emplace_back(FinalStartT, FinalEndT, peakVals, NFluctuations);
1184  }
1185  return;
1186 }
hit::DPRawHitFinder::PeakTimeWidVec
std::vector< std::tuple< int, int, int, int >> PeakTimeWidVec
Definition: DPRawHitFinder_module.cc:81
hit::DPRawHitFinder::EstimateFluctuations
int EstimateFluctuations(const std::vector< float > fsignalVec, int peakStart, int peakMean, int peakEnd)
Definition: DPRawHitFinder_module.cc:1191
void hit::DPRawHitFinder::produce ( art::Event &  evt)
overrideprivatevirtual

Implements art::EDProducer.

Definition at line 298 of file DPRawHitFinder_module.cc.

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

Definition at line 1797 of file DPRawHitFinder_module.cc.

1800 {
1801  size_t nNewBins = inputVec.size() / nBinsToCombine;
1802 
1803  if (inputVec.size() % nBinsToCombine > 0) nNewBins++;
1804 
1805  outputVec.resize(nNewBins, 0.);
1806 
1807  size_t outputBin = 0;
1808 
1809  for(size_t inputIdx = 0; inputIdx < inputVec.size();)
1810  {
1811  outputVec[outputBin] += inputVec[inputIdx++];
1812 
1813  if (inputIdx % nBinsToCombine == 0) outputBin++;
1814 
1815  if (outputBin > outputVec.size())
1816  {
1817  std::cout << "***** DISASTER!!! ****** outputBin: " << outputBin << ", inputIdx = " << inputIdx << std::endl;
1818  break;
1819  }
1820  }
1821 
1822  return;
1823 }
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void hit::DPRawHitFinder::SplitPeak ( std::tuple< double, int, int, int >  fPeakDevCand,
PeakTimeWidVec &  fpeakValsTemp 
)
private

Definition at line 1609 of file DPRawHitFinder_module.cc.

1611 {
1612 int PeakNumberWithNewPeak = std::get<1>(fPeakDevCand);
1613 int OldPeakOldStart = std::get<2>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1614 int OldPeakOldEnd = std::get<3>(fpeakValsTemp.at(PeakNumberWithNewPeak));
1615 int WidthOldPeakOld = OldPeakOldEnd - OldPeakOldStart;
1616 
1617 if(WidthOldPeakOld<3) {return;} //can't split peaks with widths < 3
1618 
1619 int NewPeakMax = 0;
1620 int NewPeakStart = 0;
1621 int NewPeakEnd = 0;
1622 int OldPeakNewMax = 0;
1623 int OldPeakNewStart = 0;
1624 int OldPeakNewEnd = 0;
1625 
1626 
1627 OldPeakNewStart = OldPeakOldStart;
1628 NewPeakEnd = OldPeakOldEnd;
1629 
1630 OldPeakNewEnd = OldPeakNewStart + 0.5*(WidthOldPeakOld + (WidthOldPeakOld%2));
1631 NewPeakStart = OldPeakNewEnd+1;
1632 
1633 int WidthOldPeakNew = OldPeakNewEnd-OldPeakNewStart;
1634 int WidthNewPeak = NewPeakEnd-NewPeakStart;
1635 
1636 OldPeakNewMax = OldPeakNewStart + 0.5*(WidthOldPeakNew - (WidthOldPeakNew%2));
1637 NewPeakMax = NewPeakStart + 0.5*(WidthNewPeak - (WidthNewPeak%2));
1638 
1639 fpeakValsTemp.at(PeakNumberWithNewPeak) = std::make_tuple(OldPeakNewMax,0,OldPeakNewStart,OldPeakNewEnd);
1640 fpeakValsTemp.emplace_back(NewPeakMax,0,NewPeakStart,NewPeakEnd);
1641 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);} );
1642 
1643 return;
1644 }
test_gpu_visible.t1
t1
Definition: test_gpu_visible.py:16
double hit::DPRawHitFinder::WidthFunc ( double  fPeakMean,
double  fPeakAmp,
double  fPeakTau1,
double  fPeakTau2,
double  fStartTime,
double  fEndTime,
double  fPeakMeanTrue 
)
private

Definition at line 1647 of file DPRawHitFinder_module.cc.

1654 {
1655 double MaxValue = ( fPeakAmp * exp(0.4*(fPeakMeanTrue-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(fPeakMeanTrue-fPeakMean)/fPeakTau2) );
1656 double FuncValue = 0.;
1657 double HalfMaxLeftTime = 0.;
1658 double HalfMaxRightTime = 0.;
1659 double ReBin=10.;
1660 
1661  //First iteration (+- 1 bin)
1662  for(double x = fPeakMeanTrue; x > fStartTime-1000.; x--)
1663  {
1664  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1665  if( FuncValue < 0.5*MaxValue )
1666  {
1667  HalfMaxLeftTime = x;
1668  break;
1669  }
1670  }
1671 
1672  for(double x = fPeakMeanTrue; x < fEndTime+1000.; x++)
1673  {
1674  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1675  if( FuncValue < 0.5*MaxValue )
1676  {
1677  HalfMaxRightTime = x;
1678  break;
1679  }
1680  }
1681 
1682  //Second iteration (+- 0.1 bin)
1683  for(double x = HalfMaxLeftTime+1; x > HalfMaxLeftTime; x-=(1/ReBin))
1684  {
1685  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1686  if( FuncValue < 0.5*MaxValue )
1687  {
1688  HalfMaxLeftTime = x;
1689  break;
1690  }
1691  }
1692 
1693  for(double x = HalfMaxRightTime-1; x < HalfMaxRightTime; x+=(1/ReBin))
1694  {
1695  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1696  if( FuncValue < 0.5*MaxValue )
1697  {
1698  HalfMaxRightTime = x;
1699  break;
1700  }
1701  }
1702 
1703  //Third iteration (+- 0.01 bin)
1704  for(double x = HalfMaxLeftTime+1/ReBin; x > HalfMaxLeftTime; x-=1/(ReBin*ReBin))
1705  {
1706  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1707  if( FuncValue < 0.5*MaxValue )
1708  {
1709  HalfMaxLeftTime = x;
1710  break;
1711  }
1712  }
1713 
1714  for(double x = HalfMaxRightTime-1/ReBin; x < HalfMaxRightTime; x+=1/(ReBin*ReBin))
1715  {
1716  FuncValue = ( fPeakAmp * exp(0.4*(x-fPeakMean)/fPeakTau1)) / ( 1 + exp(0.4*(x-fPeakMean)/fPeakTau2) );
1717  if( FuncValue < 0.5*MaxValue )
1718  {
1719  HalfMaxRightTime = x;
1720  break;
1721  }
1722  }
1723 
1724 return HalfMaxRightTime-HalfMaxLeftTime;
1725 }
train.x
list x
Definition: train.py:276

Member Data Documentation

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

Definition at line 165 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChargeNorm
private

Definition at line 176 of file DPRawHitFinder_module.cc.

TH1F* hit::DPRawHitFinder::fChi2
private

Definition at line 204 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFMax
private

Definition at line 180 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFMaxFactorMultiHits
private

Definition at line 181 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFRetry
private

Definition at line 178 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFRetryFactorMultiHits
private

Definition at line 179 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fDoMergePeaks
private

Definition at line 189 of file DPRawHitFinder_module.cc.

TH1F* hit::DPRawHitFinder::fFirstChi2
private

Definition at line 203 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fFitPeakMeanRange
private

Definition at line 185 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fGroupMaxDistance
private

Definition at line 186 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fGroupMinADC
private

Definition at line 187 of file DPRawHitFinder_module.cc.

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

Definition at line 201 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLogLevel
private

Definition at line 168 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLongMaxHits
private

Definition at line 196 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLongPulseWidth
private

Definition at line 197 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxFluctuations
private

Definition at line 198 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxGroupLength
private

Definition at line 175 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxMultiHit
private

Definition at line 174 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMaxTau
private

Definition at line 184 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeADCSumThreshold
private

Definition at line 190 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeMaxADCLimit
private

Definition at line 194 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeMaxADCThreshold
private

Definition at line 191 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinADCSum
private

Definition at line 172 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinADCSumOverWidth
private

Definition at line 173 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinRelativePeakHeightLeft
private

Definition at line 192 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinRelativePeakHeightRight
private

Definition at line 193 of file DPRawHitFinder_module.cc.

float hit::DPRawHitFinder::fMinSig
private

Definition at line 169 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinTau
private

Definition at line 183 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMinWidth
private

Definition at line 171 of file DPRawHitFinder_module.cc.

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

Definition at line 200 of file DPRawHitFinder_module.cc.

size_t hit::DPRawHitFinder::fNumBinsToAverage
private

Definition at line 182 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fSameShape
private

Definition at line 188 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fTicksToStopPeakFinder
private

Definition at line 170 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fTryNplus1Fits
private

Definition at line 177 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fWidthNormalization
private

Definition at line 195 of file DPRawHitFinder_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11