icarus_tool::ROIDeconvolution Class Reference

Inheritance diagram for icarus_tool::ROIDeconvolution:

Public Member Functions
	ROIDeconvolution (const fhicl::ParameterSet &pset)
	~ROIDeconvolution ()
void	configure (const fhicl::ParameterSet &pset) override
void	initializeHistograms (art::TFileDirectory &) const override
void	Deconvolve (const IROIFinder::Waveform &, double samplingRate, raw::ChannelID_t, IROIFinder::CandidateROIVec const &, recob::Wire::RegionsOfInterest_t &) const override
Public Member Functions inherited from icarus_tool::IDeconvolution
virtual	~IDeconvolution () noexcept=default

Private Attributes
size_t	fFFTSize
	FFT size for ROI deconvolution. More...
bool	fDodQdxCalib
	Do we apply wire-by-wire calibration? More...
std::string	fdQdxCalibFileName
	Text file for constants to do wire-by-wire calibration. More...
std::map< unsigned int, float >	fdQdxCalib
	number, content is correction factor More...
std::unique_ptr < icarus_tool::IBaseline >	fBaseline
const geo::GeometryCore *	fGeometry = lar::providerFrom<geo::Geometry>()
art::ServiceHandle < icarusutil::SignalShapingICARUSService >	fSignalShaping
std::unique_ptr < icarus_signal_processing::ICARUSFFT < double > >	fFFT
	Object to handle thread safe FFT. More...

Detailed Description

Definition at line 28 of file ROIDeconvolution_tool.cc.

Constructor & Destructor Documentation

icarus_tool::ROIDeconvolution::ROIDeconvolution ( const fhicl::ParameterSet &  pset )

explicit

Definition at line 61 of file ROIDeconvolution_tool.cc.

{
    configure(pset);
}

icarus_tool::ROIDeconvolution::~ROIDeconvolution

(

)

Definition at line 66 of file ROIDeconvolution_tool.cc.

{
}

Member Function Documentation

void icarus_tool::ROIDeconvolution::configure ( const fhicl::ParameterSet &  pset )

overridevirtual

Implements icarus_tool::IDeconvolution.

Definition at line 70 of file ROIDeconvolution_tool.cc.

{
    // Start by recovering the parameters
    fFFTSize    = pset.get< size_t >("FFTSize"                );
    
    //wire-by-wire calibration
    fDodQdxCalib = pset.get< bool >("DodQdxCalib", false);
    
    if (fDodQdxCalib)
    {
        fdQdxCalibFileName = pset.get< std::string >("dQdxCalibFileName");
        std::string fullname;
        cet::search_path sp("FW_SEARCH_PATH");
        sp.find_file(fdQdxCalibFileName, fullname);
        
        if (fullname.empty())
        {
            std::cout << "Input file " << fdQdxCalibFileName << " not found" << std::endl;
            throw cet::exception("File not found");
        }
        else
            std::cout << "Applying wire-by-wire calibration using file " << fdQdxCalibFileName << std::endl;
        
        std::ifstream inFile(fullname, std::ios::in);
        std::string line;
        
        while (std::getline(inFile,line))
        {
            unsigned int channel;
            float        constant;
            std::stringstream linestream(line);
            linestream >> channel >> constant;
            fdQdxCalib[channel] = constant;
            if (channel%1000==0) std::cout<<"Channel "<<channel<<" correction factor "<<fdQdxCalib[channel]<<std::endl;
        }
    }
    
    // Recover the baseline tool
    fBaseline  = art::make_tool<icarus_tool::IBaseline> (pset.get<fhicl::ParameterSet>("Baseline"));
    
    // Get signal shaping service.
    fSignalShaping = art::ServiceHandle<icarusutil::SignalShapingICARUSService>();

    // Now set up our plans for doing the convolution
    auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
    fFFT = std::make_unique<icarus_signal_processing::ICARUSFFT<double>>(detProp.NumberTimeSamples());
    
    return;
}

void icarus_tool::ROIDeconvolution::Deconvolve ( const IROIFinder::Waveform &  waveform, double  samplingRate, raw::ChannelID_t  channel, IROIFinder::CandidateROIVec const &  roiVec, recob::Wire::RegionsOfInterest_t &  ROIVec  ) const

overridevirtual

Implements icarus_tool::IDeconvolution.

Definition at line 119 of file ROIDeconvolution_tool.cc.

{
    double deconNorm = fSignalShaping->GetDeconNorm();

    // And now process them
    for(auto const& roi : roiVec)
    {
        // First up: copy out the relevent ADC bins into the ROI holder
        size_t roiLen = roi.second - roi.first;
        
        // We want the deconvolution buffer size to be a power of 2 in length
        // to facilitate the FFT
        size_t deconSize = fFFTSize;
        
        while(1)
        {
            if (roiLen > deconSize  ) deconSize *= 2;
            else break;
        }
        
        // In theory, most ROI's are around the same size so this should mostly be a noop
        fSignalShaping->SetDecon(samplingRate, deconSize, channel);
        
        deconSize = fFFTSize;

        icarusutil::TimeVec deconVec(deconSize);
        
        // Pad with zeroes if the deconvolution buffer is larger than the input waveform
        if (deconSize > waveform.size()) deconVec.resize(deconSize, 0.);
        
        // Watch for the case where the input ROI is long enough to want an deconvolution buffer that is
        // larger than the input waveform.
        size_t maxActualSize = std::min(deconSize, waveform.size());
        
        // Extend the ROI to accommodate the extra bins for the FFT
        // The idea is to try to center the desired ROI in the buffer used by deconvolution
        size_t halfLeftOver = (maxActualSize - roiLen) / 2;           // Number bins either side of ROI
        int    roiStartInt  = halfLeftOver;                           // Start in the buffer of the ROI
        int    roiStopInt   = halfLeftOver + roiLen;                  // Stop in the buffer of the ROI
        int    firstOffset  = roi.first - halfLeftOver;               // Offset into the ADC vector of buffer start
        int    secondOffset = roi.second + halfLeftOver + roiLen % 2; // Offset into the ADC vector of buffer end
        
        // Check for the two edge conditions - starting before the ADC vector or running off the end
        // In either case we shift the actual roi within the FFT buffer
        // First is the case where we would be starting before the ADC vector
        if (firstOffset < 0)
        {
            roiStartInt  += firstOffset;  // remember that firstOffset is negative
            roiStopInt   += firstOffset;
            secondOffset -= firstOffset;
            firstOffset   = 0;
        }
        // Second is the case where we would overshoot the end
        else if (size_t(secondOffset) > waveform.size())
        {
            size_t overshoot = secondOffset - waveform.size();
            
            roiStartInt  += overshoot;
            roiStopInt   += overshoot;
            firstOffset  -= overshoot;
            secondOffset  = waveform.size();
        }
        
        size_t roiStart(roiStartInt);
        size_t roiStop(roiStopInt);
        size_t holderOffset = 0; //deconSize > waveform.size() ? (deconSize - waveform.size()) / 2 : 0;
        
        // Fill the buffer and do the deconvolution
        std::copy(waveform.begin()+firstOffset, waveform.begin()+secondOffset, deconVec.begin() + holderOffset);
        
        // Deconvolute the raw signal using the channel's nominal response
        fFFT->deconvolute(deconVec, fSignalShaping->GetResponse(channel).getDeconvKernel(), fSignalShaping->ResponseTOffset(channel));

        std::vector<float>  holder(deconVec.size());
        
        // Get rid of the leading and trailing "extra" bins needed to keep the FFT happy
        if (roiStart > 0 || holderOffset > 0) std::copy(deconVec.begin() + holderOffset + roiStart, deconVec.begin() + holderOffset + roiStop, holder.begin());
        
        // Resize the holder to the ROI length
        holder.resize(roiLen);
       
        // "normalize" the vector
        std::transform(holder.begin(),holder.end(),holder.begin(),[deconNorm](auto& deconVal){return deconVal/deconNorm;});
        
        // Now we do the baseline determination and correct the ROI
        //float base = fBaseline->GetBaseline(holder, channel, roiStart, roiLen);
        float base = fBaseline->GetBaseline(holder, channel, 0, roiLen);
        
        std::transform(holder.begin(),holder.end(),holder.begin(),[base](const auto& adcVal){return adcVal - base;});
        
        // apply wire-by-wire calibration
        if (fDodQdxCalib)
        {
            if(fdQdxCalib.find(channel) != fdQdxCalib.end())
            {
                float constant = fdQdxCalib.at(channel);
                
                for (size_t iholder = 0; iholder < holder.size(); ++iholder) holder[iholder] *= constant;
            }
        }

        // add the range into ROIVec
        ROIVec.add_range(roi.first, std::move(holder));
    } // loop over candidate roi's
    
    return;
}

void icarus_tool::ROIDeconvolution::initializeHistograms ( art::TFileDirectory &  histDir ) const

overridevirtual

Implements icarus_tool::IDeconvolution.

Definition at line 233 of file ROIDeconvolution_tool.cc.

{
    // It is assumed that the input TFileDirectory has been set up to group histograms into a common
    // folder at the calling routine's level. Here we create one more level of indirection to keep
    // histograms made by this tool separate.
/*
    std::string dirName = "ROIDeconvolutionPlane_" + std::to_string(fPlane);
 
    art::TFileDirectory dir = histDir.mkdir(dirName.c_str());
 
    auto const* detprop      = lar::providerFrom<detinfo::DetectorPropertiesService>();
    double      samplingRate = detprop->SamplingRate();
    double      numBins      = fROIDeconvolutionVec.size();
    double      maxFreq      = 500. / samplingRate;
    std::string histName     = "ROIDeconvolutionPlane_" + std::to_string(fPlane);
    
    TH1D*       hist         = dir.make<TH1D>(histName.c_str(), "ROIDeconvolution;Frequency(MHz)", numBins, 0., maxFreq);
    
    for(int bin = 0; bin < numBins; bin++)
    {
        double freq = maxFreq * double(bin + 0.5) / double(numBins);
        
        hist->Fill(freq, fROIDeconvolutionVec.at(bin).Re());
    }
*/
    
    return;
}

Member Data Documentation

std::unique_ptr<icarus_tool::IBaseline> icarus_tool::ROIDeconvolution::fBaseline

private

Definition at line 52 of file ROIDeconvolution_tool.cc.

bool icarus_tool::ROIDeconvolution::fDodQdxCalib

private

Do we apply wire-by-wire calibration?

Definition at line 47 of file ROIDeconvolution_tool.cc.

std::map<unsigned int, float> icarus_tool::ROIDeconvolution::fdQdxCalib

private

number, content is correction factor

Map to do wire-by-wire calibration, key is channel

Definition at line 49 of file ROIDeconvolution_tool.cc.

std::string icarus_tool::ROIDeconvolution::fdQdxCalibFileName

private

Text file for constants to do wire-by-wire calibration.

Definition at line 48 of file ROIDeconvolution_tool.cc.

std::unique_ptr<icarus_signal_processing::ICARUSFFT<double> > icarus_tool::ROIDeconvolution::fFFT

private

Object to handle thread safe FFT.

Definition at line 56 of file ROIDeconvolution_tool.cc.

size_t icarus_tool::ROIDeconvolution::fFFTSize

private

FFT size for ROI deconvolution.

Definition at line 46 of file ROIDeconvolution_tool.cc.

const geo::GeometryCore* icarus_tool::ROIDeconvolution::fGeometry = lar::providerFrom<geo::Geometry>()

private

Definition at line 54 of file ROIDeconvolution_tool.cc.

art::ServiceHandle<icarusutil::SignalShapingICARUSService> icarus_tool::ROIDeconvolution::fSignalShaping

private

Definition at line 55 of file ROIDeconvolution_tool.cc.

---

The documentation for this class was generated from the following file:

• ROIDeconvolution_tool.cc