505 for (
size_t iq = 0; iq < 2; iq++) {
507 if (wgqvec[igrp][iq].
size() == 0)
continue;
510 wgqvec[igrp][iq].erase(std::remove_if(wgqvec[igrp][iq].
begin(),wgqvec[igrp][iq].
end(),[](
const int& i){
return i<0;}),wgqvec[igrp][iq].
end());
513 size_t nwq = wgqvec[igrp][iq].size();
514 if (nwq <= 2)
continue;
516 std::vector<float> corValVec;
517 corValVec.resize(fftSize, 0.);
522 for(
size_t itck = 0; itck < fftSize; itck++){
523 std::vector<float> adcValuesVec;
525 for (
size_t i = 0; i < nwq; i++) {
527 size_t iwdx = wgqvec[igrp][iq][i];
530 if (itck < wgcvec[igrp][iwdx].tcka || itck >= wgcvec[igrp][iwdx].tckb)
continue;
532 adcValuesVec.push_back(
float(rawadcgvec[igrp][iwdx][itck]) - wgcvec[igrp][iwdx].truncMean);
535 float medval(-10000);
536 if (!adcValuesVec.empty()) {
537 std::sort(adcValuesVec.begin(),adcValuesVec.end());
538 size_t medidx = adcValuesVec.size() / 2;
539 medval = adcValuesVec[medidx];
540 if (adcValuesVec.size() > 1 && medidx % 2) medval = (medval + adcValuesVec[medidx+1]) / 2;
542 corValVec[itck] = std::max(medval,
float(-10000.));
546 size_t iwdx0 = wgqvec[igrp][iq][0];
547 unsigned int plane = wgcvec[igrp][iwdx0].plane;
553 std::vector<float> localCorValVec = corValVec;
554 std::sort(localCorValVec.begin(),localCorValVec.end());
557 float corValSum = std::accumulate(localCorValVec.begin(),localCorValVec.begin() + nTruncVal,0.);
558 float meanCorVal = corValSum / float(nTruncVal);
560 std::vector<float> diffVec(nTruncVal);
561 std::transform(localCorValVec.begin(),localCorValVec.begin() + nTruncVal, diffVec.begin(),
562 std::bind(std::minus<float>(),std::placeholders::_1,meanCorVal));
564 float rmsValSq = std::inner_product(diffVec.begin(),diffVec.end(),diffVec.begin(),0.);
565 float rmsVal = std::sqrt(rmsValSq /
float(nTruncVal));
568 std::vector<float>::iterator lastGoodItr = corValVec.begin();
569 bool wasOutlier(
false);
570 for(std::vector<float>::iterator corValItr = lastGoodItr+1; corValItr != corValVec.end(); corValItr++){
573 float lastVal = *lastGoodItr;
574 float curVal = *corValItr;
576 float slope = (curVal - lastVal) / numTicks;
577 while(lastGoodItr != corValItr){
578 *lastGoodItr++ = (numTicks -
std::distance(lastGoodItr,corValItr)) * slope + lastVal;
582 lastGoodItr = corValItr;
591 std::vector<std::complex<double>> fftOutputVec(halfFFTSize);
593 lfftw.DoFFT(corValVec, fftOutputVec);
595 std::vector<double> powerVec(halfFFTSize);
596 std::transform(fftOutputVec.begin(), fftOutputVec.begin() + halfFFTSize, powerVec.begin(), [](
const auto& val){
return std::abs(val);});
599 std::vector<double> firstDerivVec(powerVec.size(), 0.);
602 for(
size_t idx = 1; idx < firstDerivVec.size() - 1; idx++)
603 firstDerivVec.at(idx) = 0.5 * (powerVec.at(idx + 1) - powerVec.at(idx - 1));
606 std::vector<std::tuple<size_t,size_t,size_t>> peakTupleVec;
610 if (!peakTupleVec.empty())
612 for(
const auto& peakTuple : peakTupleVec)
614 size_t startTick = std::get<0>(peakTuple);
615 size_t stopTick = std::get<2>(peakTuple);
617 if (stopTick > startTick)
619 std::complex<double> slope = (fftOutputVec[stopTick] - fftOutputVec[startTick]) /
double(stopTick - startTick);
623 std::complex<double> interpVal = fftOutputVec[startTick] + double(
tick - startTick) * slope;
625 fftOutputVec[
tick] = interpVal;
631 std::vector<double> tmpVec(corValVec.size());
633 lfftw.DoInvFFT(fftOutputVec, tmpVec);
635 std::transform(corValVec.begin(),corValVec.end(),tmpVec.begin(),corValVec.begin(),std::minus<double>());
642 for(
size_t itck = 0; itck < fftSize; itck++){
643 for (
size_t i = 0; i < nwq; i++) {
645 size_t iwdx = wgqvec[igrp][iq][i];
649 if (itck < wgcvec[igrp][iwdx].tcka || itck >= wgcvec[igrp][iwdx].tckb) {
652 corVal = corValVec[itck];
655 float newAdcValueFloat = float(rawadcgvec[igrp][iwdx][itck]) - corVal - wgcvec[igrp][iwdx].pedCor;
656 rawadcgvec[igrp][iwdx][itck] = std::round(newAdcValueFloat);
664 for (
size_t iwdx = 0; iwdx < wgcvec[igrp].size(); iwdx++) {
666 unsigned int plane = wgcvec[igrp][iwdx].plane;
675 float pedestal = wgcvec[igrp][iwdx].truncMean;
676 float pedCor = wgcvec[igrp][iwdx].pedCor;
677 float deltaPed = pedestal - pedCor;
685 int irdg = wgcvec[igrp][iwdx].irawdig;
688 filteredRawDigit->at(irdg).SetPedestal(pedestal, rmsVal);
690 mf::LogInfo(
"RawDigitFilterICARUS") <<
"--> Rejecting channel for large rms, channel: "
691 << channel <<
", rmsVal: " << rmsVal <<
", truncMean: " << pedestal
692 <<
", pedestal: " << pedCor << std::endl;
bool fApplyFFTCorrection
Use an FFT to get the correlated noise correction.
Collection of charge vs time digitized from a single readout channel.
float fTruncMeanFraction
Fraction for truncated mean.
raw::RawDigit::ADCvector_t RawDigitVector
std::size_t size(FixedBins< T, C > const &) noexcept
std::vector< float > fRmsRejectionCutHi
Maximum rms for input channels, reject if larger.
caldata::RawDigitCharacterizationAlg fCharacterizationAlg
void findPeaks(typename std::vector< T >::iterator startItr, typename std::vector< T >::iterator stopItr, std::vector< std::tuple< size_t, size_t, size_t >> &peakTupleVec, T threshold, size_t firstTick) const
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
bool fApplyCorSmoothing
Attempt to smooth the correlated noise correction?
bool fApplyTopHatFilter
Apply the top hat filter.
auto end(FixedBins< T, C > const &) noexcept
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
std::vector< double > fFFTMinPowerThreshold
Threshold for trimming FFT power spectrum.
void getTruncatedRMS(const RawDigitVector &rawWaveform, float &pedestal, float &truncRms) const
auto begin(FixedBins< T, C > const &) noexcept
std::unique_ptr< caldata::IRawDigitFilter > fRawDigitFilterTool
std::vector< float > fNumRmsToSmoothVec
"sigma" to smooth correlated correction vec
unsigned int ChannelID_t
Type representing the ID of a readout channel.
1.8.5