All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
hit::DPRawHitFinder Class Reference
Inheritance diagram for hit::DPRawHitFinder:

Classes

struct  Comp
 

Public Member Functions

 DPRawHitFinder (fhicl::ParameterSet const &pset)
 

Private Types

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

Private Member Functions

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

Private Attributes

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

Detailed Description

Definition at line 70 of file DPRawHitFinder_module.cc.

Member Typedef Documentation

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

Definition at line 83 of file DPRawHitFinder_module.cc.

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

Definition at line 85 of file DPRawHitFinder_module.cc.

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

Definition at line 84 of file DPRawHitFinder_module.cc.

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

Definition at line 82 of file DPRawHitFinder_module.cc.

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

Definition at line 81 of file DPRawHitFinder_module.cc.

Constructor & Destructor Documentation

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

Definition at line 212 of file DPRawHitFinder_module.cc.

212  :
213  EDProducer{pset},
214  fNewHitsTag(
215  pset.get<std::string>("module_label"), "",
216  art::ServiceHandle<art::TriggerNamesService const>()->getProcessName()),
217  fHitParamWriter(producesCollector())
218 {
219  fLogLevel = pset.get< int >("LogLevel");
220  fCalDataModuleLabel = pset.get< std::string >("CalDataModuleLabel");
221  fMaxMultiHit = pset.get< int >("MaxMultiHit");
222  fMaxGroupLength = pset.get< int >("MaxGroupLength");
223  fTryNplus1Fits = pset.get< bool >("TryNplus1Fits");
224  fChi2NDFRetry = pset.get< double >("Chi2NDFRetry");
225  fChi2NDFRetryFactorMultiHits = pset.get< double >("Chi2NDFRetryFactorMultiHits");
226  fChi2NDFMax = pset.get< double >("Chi2NDFMax");
227  fChi2NDFMaxFactorMultiHits = pset.get< double >("Chi2NDFMaxFactorMultiHits");
228  fNumBinsToAverage = pset.get< size_t >("NumBinsToAverage");
229  fMinSig = pset.get< float >("MinSig");
230  fMinWidth = pset.get< int >("MinWidth");
231  fMinADCSum = pset.get< double >("MinADCSum");
232  fMinADCSumOverWidth = pset.get< double >("MinADCSumOverWidth");
233  fChargeNorm = pset.get< double >("ChargeNorm");
234  fTicksToStopPeakFinder = pset.get< double >("TicksToStopPeakFinder");
235  fMinTau = pset.get< double >("MinTau");
236  fMaxTau = pset.get< double >("MaxTau");
237  fFitPeakMeanRange = pset.get< double >("FitPeakMeanRange");
238  fGroupMaxDistance = pset.get< int >("GroupMaxDistance");
239  fGroupMinADC = pset.get< int >("GroupMinADC");
240  fSameShape = pset.get< bool >("SameShape");
241  fDoMergePeaks = pset.get< bool >("DoMergePeaks");
242  fMergeADCSumThreshold = pset.get< double >("MergeADCSumThreshold");
243  fMergeMaxADCThreshold = pset.get< double >("MergeMaxADCThreshold");
244  fMinRelativePeakHeightLeft = pset.get< double >("MinRelativePeakHeightLeft");
245  fMinRelativePeakHeightRight = pset.get< double >("MinRelativePeakHeightRight");
246  fMergeMaxADCLimit = pset.get< double >("MergeMaxADCLimit");
247  fWidthNormalization = pset.get< double >("WidthNormalization");
248  fLongMaxHits = pset.get< double >("LongMaxHits");
249  fLongPulseWidth = pset.get< double >("LongPulseWidth");
250  fMaxFluctuations = pset.get< double >("MaxFluctuations");
251 
252  // let HitCollectionCreator declare that we are going to produce
253  // hits and associations with wires and raw digits
254  // (with no particular product label)
255  recob::HitCollectionCreator::declare_products(producesCollector());
256 
257  // declare that data products with feature vectors describing
258  // hits is going to be produced
259  fHitParamWriter.produces_using< recob::Hit >();
260 
261 } // DPRawHitFinder::DPRawHitFinder()
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:202
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 1566 of file DPRawHitFinder_module.cc.

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

Definition at line 287 of file DPRawHitFinder_module.cc.

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

Definition at line 1729 of file DPRawHitFinder_module.cc.

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

Definition at line 1502 of file DPRawHitFinder_module.cc.

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

Definition at line 1768 of file DPRawHitFinder_module.cc.

1771 {
1772  size_t halfBinsToAverage(binsToAverage/2);
1773 
1774  float runningSum(0.);
1775 
1776  for(size_t idx = 0; idx < halfBinsToAverage; idx++) runningSum += inputVec[idx];
1777 
1778  outputVec.resize(inputVec.size());
1779  std::vector<float>::iterator outputVecItr = outputVec.begin();
1780 
1781  // First pass through to build the erosion vector
1782  for(std::vector<float>::const_iterator inputItr = inputVec.begin(); inputItr != inputVec.end(); inputItr++)
1783  {
1784  size_t startOffset = std::distance(inputVec.begin(),inputItr);
1785  size_t stopOffset = std::distance(inputItr,inputVec.end());
1786  size_t count = std::min(2 * halfBinsToAverage, std::min(startOffset + halfBinsToAverage + 1, halfBinsToAverage + stopOffset - 1));
1787 
1788  if (startOffset >= halfBinsToAverage) runningSum -= *(inputItr - halfBinsToAverage);
1789  if (stopOffset > halfBinsToAverage) runningSum += *(inputItr + halfBinsToAverage);
1790 
1791  *outputVecItr++ = runningSum / float(count);
1792  }
1793 
1794  return;
1795 }
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
count
std::size_t count(Cont const &cont)
Definition: IntegerRanges_test.cc:29
int hit::DPRawHitFinder::EstimateFluctuations ( const std::vector< float >  fsignalVec,
int  peakStart,
int  peakMean,
int  peakEnd 
)
private

Definition at line 1192 of file DPRawHitFinder_module.cc.

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

Definition at line 266 of file DPRawHitFinder_module.cc.

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

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

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

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

Definition at line 1062 of file DPRawHitFinder_module.cc.

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

Definition at line 299 of file DPRawHitFinder_module.cc.

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

Definition at line 1798 of file DPRawHitFinder_module.cc.

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

Definition at line 1610 of file DPRawHitFinder_module.cc.

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

Definition at line 1648 of file DPRawHitFinder_module.cc.

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

Member Data Documentation

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

Definition at line 166 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChargeNorm
private

Definition at line 177 of file DPRawHitFinder_module.cc.

TH1F* hit::DPRawHitFinder::fChi2
private

Definition at line 205 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFMax
private

Definition at line 181 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFMaxFactorMultiHits
private

Definition at line 182 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFRetry
private

Definition at line 179 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fChi2NDFRetryFactorMultiHits
private

Definition at line 180 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fDoMergePeaks
private

Definition at line 190 of file DPRawHitFinder_module.cc.

TH1F* hit::DPRawHitFinder::fFirstChi2
private

Definition at line 204 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fFitPeakMeanRange
private

Definition at line 186 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fGroupMaxDistance
private

Definition at line 187 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fGroupMinADC
private

Definition at line 188 of file DPRawHitFinder_module.cc.

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

Definition at line 202 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLogLevel
private

Definition at line 169 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLongMaxHits
private

Definition at line 197 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fLongPulseWidth
private

Definition at line 198 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxFluctuations
private

Definition at line 199 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxGroupLength
private

Definition at line 176 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMaxMultiHit
private

Definition at line 175 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMaxTau
private

Definition at line 185 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeADCSumThreshold
private

Definition at line 191 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeMaxADCLimit
private

Definition at line 195 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMergeMaxADCThreshold
private

Definition at line 192 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinADCSum
private

Definition at line 173 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinADCSumOverWidth
private

Definition at line 174 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinRelativePeakHeightLeft
private

Definition at line 193 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinRelativePeakHeightRight
private

Definition at line 194 of file DPRawHitFinder_module.cc.

float hit::DPRawHitFinder::fMinSig
private

Definition at line 170 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fMinTau
private

Definition at line 184 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fMinWidth
private

Definition at line 172 of file DPRawHitFinder_module.cc.

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

Definition at line 201 of file DPRawHitFinder_module.cc.

size_t hit::DPRawHitFinder::fNumBinsToAverage
private

Definition at line 183 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fSameShape
private

Definition at line 189 of file DPRawHitFinder_module.cc.

int hit::DPRawHitFinder::fTicksToStopPeakFinder
private

Definition at line 171 of file DPRawHitFinder_module.cc.

bool hit::DPRawHitFinder::fTryNplus1Fits
private

Definition at line 178 of file DPRawHitFinder_module.cc.

double hit::DPRawHitFinder::fWidthNormalization
private

Definition at line 196 of file DPRawHitFinder_module.cc.

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