223 TH1::AddDirectory(kFALSE);
253 tbb::concurrent_vector<hitstruct> hitstruct_vec;
254 tbb::concurrent_vector<hitstruct> filthitstruct_vec;
268 std::function<double(double, double, double, double, int, int)> chargeFunc =
269 [](
double peakMean,
double peakAmp,
double peakWidth,
double areaNorm,
int low,
int hi) {
270 return std::sqrt(2 * TMath::Pi()) * peakAmp * peakWidth / areaNorm;
278 [](
double peakMean,
double peakAmp,
double peakWidth,
double areaNorm,
int low,
int hi) {
280 for (
int sigPos = low; sigPos < hi; sigPos++)
281 charge += peakAmp * TMath::Gaus(sigPos, peakMean, peakWidth);
291 static_cast<std::size_t>(0),
292 wireVecHandle->size(),
293 [&](
size_t& wireIter) {
304 std::vector<geo::WireID> wids = geom->
ChannelToWire(channel);
320 for (
const auto& range : signalROI.
get_ranges()) {
331 fHitFinderToolVec.at(plane)->findHitCandidates(range, 0, channel, count, hitCandidateVec);
333 range, hitCandidateVec, mergedCandidateHitVec);
339 for (
auto& mergedCands : mergedCandidateHitVec) {
340 int startT = mergedCands.front().startTick;
341 int endT = mergedCands.back().stopTick;
346 if (endT - startT < 5)
continue;
353 int nGausForFit = mergedCands.size();
358 double chi2PerNDF(0.);
370 range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF);
373 if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) {
389 int nHitsThisPulse = (endT - startT) / longPulseWidth;
393 longPulseWidth = (endT - startT) / nHitsThisPulse;
396 if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++;
398 int firstTick = startT;
399 int lastTick =
std::min(firstTick + longPulseWidth, endT);
401 peakParamsVec.clear();
402 nGausForFit = nHitsThisPulse;
404 chi2PerNDF = chi2PerNDF >
fChi2NDF ? chi2PerNDF : -1.;
406 for (
int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) {
409 std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.);
410 double peakSigma = (lastTick - firstTick) / 3.;
411 double peakAmp = 0.3989 * sumADC / peakSigma;
412 double peakMean = (firstTick + lastTick) / 2.;
424 peakParamsVec.push_back(peakParams);
427 firstTick = lastTick;
428 lastTick =
std::min(lastTick + longPulseWidth, endT);
439 std::vector<recob::Hit> filteredHitVec;
441 for (
const auto& peakParams : peakParamsVec) {
443 float peakAmp = peakParams.peakAmplitude;
444 float peakMean = peakParams.peakCenter;
445 float peakWidth = peakParams.peakSigma;
448 if (std::isnan(peakAmp)) {
449 std::cout <<
"**** hit peak amplitude is a nan! Channel: " << channel
450 <<
", start tick: " << startT <<
std::endl;
455 float peakAmpErr = peakParams.peakAmplitudeError;
456 float peakMeanErr = peakParams.peakCenterError;
457 float peakWidthErr = peakParams.peakSigmaError;
461 chargeFunc(peakMean, peakAmp, peakWidth,
fAreaNormsVec[plane], startT, endT);
464 std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
471 double sumADC = std::accumulate(sumStartItr, sumEndItr, 0.);
477 startT + roiFirstBinTick,
478 endT + roiFirstBinTick,
480 peakMean + roiFirstBinTick,
493 if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy());
505 if (filteredHitCol && !filteredHitVec.empty()) {
514 std::sort(filteredHitVec.begin(),
515 filteredHitVec.end(),
516 [](
const auto&
left,
const auto&
right) {
517 return left.PeakAmplitude() >
right.PeakAmplitude();
523 filteredHitVec.clear();
526 if (filteredHitVec.size() > 1) {
528 float largestPH = filteredHitVec.front().PeakAmplitude();
534 filteredHitVec.begin(),
535 filteredHitVec.end(),
536 [largestPH, threshold](
const auto&
hit) {
537 return hit.PeakAmplitude() < 8. &&
hit.PeakAmplitude() / largestPH < threshold;
541 if (smallHitItr != filteredHitVec.end())
542 filteredHitVec.resize(
std::distance(filteredHitVec.begin(), smallHitItr));
545 std::sort(filteredHitVec.begin(),
546 filteredHitVec.end(),
547 [](
const auto&
left,
const auto&
right) {
548 return left.PeakTime() <
right.PeakTime();
553 for (
const auto& filteredHit : filteredHitVec)
567 for (
size_t i = 0; i < hitstruct_vec.size(); i++) {
568 allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].wire_tbb);
571 for (
size_t j = 0; j < filthitstruct_vec.size(); j++) {
572 filteredHitCol->
emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].wire_tbb);
578 if (filteredHitCol) {
590 allHitCol.put_into(evt);
594 allHitCol.put_into(evt);
std::unique_ptr< HitFilterAlg > fHitFilterAlg
algorithm used to filter out noise hits
unsigned int PlaneID_t
Type for the ID number.
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.
const std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
std::unique_ptr< reco_tool::IPeakFitter > fPeakFitterTool
Perform fit to candidate peaks.
const int fAreaMethod
Type of area calculation.
const std::string fCalDataModuleLabel
const std::vector< float > fPulseWidthCuts
int TDCtick_t
Type representing a TDC tick.
Class managing the creation of a new recob::Hit object.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
A class handling a collection of hits and its associations.
std::atomic< size_t > fEventCount
const size_t fMaxMultiHit
maximum hits for multi fit
const std::vector< int > fLongMaxHitsVec
Maximum number hits on a really long pulse train.
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.
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
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.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const std::vector< float > fPulseRatioCuts
const std::vector< float > fPulseHeightCuts
maximum Chisquared / NDF allowed for a hit to be saved
const std::string fAllHitsInstanceName
2D representation of charge deposited in the TDC/wire plane
const std::vector< double > fAreaNormsVec
factors for converting area to same units as peak height
unsigned int ChannelID_t
Type representing the ID of a readout channel.
QTextStream & endl(QTextStream &s)
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
For finding candidate hits.