240 TH1::AddDirectory(kFALSE);
284 std::function<double (double,double,double,double,int,int)> chargeFunc = [](
double peakMean,
double peakAmp,
double peakWidth,
double areaNorm,
int low,
int hi){
return std::sqrt(2*TMath::Pi())*peakAmp*peakWidth/areaNorm;};
290 chargeFunc = [](
double peakMean,
double peakAmp,
double peakWidth,
double areaNorm,
int low,
int hi)
293 for(
int sigPos = low; sigPos < hi; sigPos++)
294 charge += peakAmp * TMath::Gaus(sigPos,peakMean,peakWidth);
301 for(
size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
313 std::vector<geo::WireID> wids = geom->
ChannelToWire(channel);
329 for(
const auto& range : signalROI.
get_ranges())
336 const std::vector<float>& signal = range.data();
354 fHitFinderToolVec.at(plane)->MergeHitCandidates(signal, hitCandidateVec, mergedCandidateHitVec);
360 for(
auto& mergedCands : mergedCandidateHitVec)
362 int startT= mergedCands.front().startTick;
363 int endT = mergedCands.back().stopTick;
368 if (endT - startT < 5)
continue;
375 int nGausForFit = mergedCands.size();
380 double chi2PerNDF(0.);
389 fPeakFitterTool->findPeakParameters(signal, mergedCands, peakParamsVec, chi2PerNDF, NDF);
392 if (!(chi2PerNDF < std::numeric_limits<double>::infinity()))
410 int nHitsThisPulse = (endT - startT) / longPulseWidth;
415 longPulseWidth = (endT - startT) / nHitsThisPulse;
418 if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++;
420 int firstTick = startT;
421 int lastTick = firstTick +
std::min(endT,longPulseWidth);
423 peakParamsVec.clear();
424 nGausForFit = nHitsThisPulse;
426 chi2PerNDF = chi2PerNDF >
fChi2NDF ? chi2PerNDF : -1.;
428 for(
int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++)
431 double sumADC = std::accumulate(signal.begin() + firstTick, signal.begin() + lastTick, 0.);
432 double peakSigma = (lastTick - firstTick) / 3.;
433 double peakAmp = 0.3989 * sumADC / peakSigma;
434 double peakMean = (firstTick + lastTick) / 2.;
439 peakParams.peakCenter = peakMean;
440 peakParams.peakCenterError = 0.1 * peakMean;
441 peakParams.peakSigma = peakSigma;
442 peakParams.peakSigmaError = 0.1 * peakSigma;
443 peakParams.peakAmplitude = peakAmp;
444 peakParams.peakAmplitudeError = 0.1 * peakAmp;
446 peakParamsVec.push_back(peakParams);
449 firstTick = lastTick;
450 lastTick =
std::min(lastTick + longPulseWidth, endT);
461 std::vector<recob::Hit> filteredHitVec;
463 for(
const auto& peakParams : peakParamsVec)
466 float peakAmp = peakParams.peakAmplitude;
467 float peakMean = peakParams.peakCenter;
468 float peakWidth = peakParams.peakSigma;
471 if (std::isnan(peakAmp))
473 std::cout <<
"**** hit peak amplitude is a nan! Channel: " << channel <<
", start tick: " << startT << std::endl;
478 float peakAmpErr = peakParams.peakAmplitudeError;
479 float peakMeanErr = peakParams.peakCenterError;
480 float peakWidthErr = peakParams.peakSigmaError;
483 float charge = chargeFunc(peakMean, peakAmp, peakWidth,
fAreaNormsVec[plane],startT,endT);;
484 float chargeErr = std::sqrt(TMath::Pi()) * (peakAmpErr*peakWidthErr + peakWidthErr*peakAmpErr);
491 double sumADC = std::accumulate(sumStartItr, sumEndItr, 0.);
496 startT+roiFirstBinTick,
497 endT+roiFirstBinTick,
499 peakMean+roiFirstBinTick,
512 filteredHitVec.push_back(hitcreator.copy());
517 allHitCol.emplace_back(std::move(
hit), wire, rawdigits);
522 if (filteredHitCol && !filteredHitVec.empty())
532 std::sort(filteredHitVec.begin(),filteredHitVec.end(),[](
const auto&
left,
const auto&
right){
return left.PeakAmplitude() >
right.PeakAmplitude();});
535 if (filteredHitVec.front().PeakAmplitude() <
fPulseHeightCuts.at(plane) || filteredHitVec.front().RMS() <
fPulseWidthCuts.at(plane)) filteredHitVec.clear();
538 if (filteredHitVec.size() > 1)
541 float largestPH = filteredHitVec.front().PeakAmplitude();
546 std::vector<recob::Hit>::iterator smallHitItr = std::find_if(filteredHitVec.begin(),filteredHitVec.end(),[largestPH,threshold](
const auto&
hit){
return hit.PeakAmplitude() < 8. &&
hit.PeakAmplitude() / largestPH < threshold;});
549 if (smallHitItr != filteredHitVec.end()) filteredHitVec.resize(std::distance(filteredHitVec.begin(),smallHitItr));
552 std::sort(filteredHitVec.begin(),filteredHitVec.end(),[](
const auto&
left,
const auto&
right){
return left.PeakTime() <
right.PeakTime();});
556 for(
const auto& filteredHit : filteredHitVec)
560 fChi2->Fill(chi2PerNDF);
572 if ( filteredHitCol ){
583 allHitCol.put_into(evt);
587 allHitCol.put_into(evt);
std::vector< double > fAreaNormsVec
factors for converting area to same units as peak height
size_t fMaxMultiHit
maximum hits for multi fit
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
std::unique_ptr< HitFilterAlg > fHitFilterAlg
algorithm used to filter out noise hits
ChannelID_t Channel() const
DAQ channel this raw data was read from.
std::vector< float > fPulseRatioCuts
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
std::vector< int > fLongMaxHitsVec
Maximum number hits on a really long pulse train.
std::unique_ptr< reco_tool::IPeakFitter > fPeakFitterTool
Perform fit to candidate peaks.
int TDCtick_t
Type representing a TDC tick.
Class managing the creation of a new recob::Hit object.
unsigned int PlaneID_t
Type for the ID number.
int fAreaMethod
Type of area calculation.
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
A class handling a collection of hits and its associations.
std::vector< float > fPulseWidthCuts
std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
PlaneID_t Plane
Index of the plane within its TPC.
void put_into(art::Event &)
Moves the data into an event.
Detector simulation of raw signals on wires.
std::vector< float > fPulseHeightCuts
maximum Chisquared / NDF allowed for a hit to be saved
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
std::string fCalDataModuleLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
2D representation of charge deposited in the TDC/wire plane
unsigned int ChannelID_t
Type representing the ID of a readout channel.
std::string fAllHitsInstanceName
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
For finding candidate hits.
1.8.11 
