Inheritance diagram for icarus_tool::Response:

Public Member Functions
	Response (const fhicl::ParameterSet &pset)

	~Response ()

void	configure (const fhicl::ParameterSet &pset) override
	Setup routines. More...

void	setResponse (double sampling_rate, double weight) override

void	outputHistograms (double sampling_rate, art::TFileDirectory &) const override

size_t	getPlane () const override
	Return the plane these functions represent. More...

const IFieldResponse *	getFieldResponse () const override
	Recover the individual response elements. More...

const IElectronicsResponse *	getElectronicsResponse () const override

const IFilter *	getFilter () const override

size_t	getNumberTimeSamples () const override
	here recover the combined response elements More...

const icarusutil::TimeVec &	getResponse () const override

const icarusutil::FrequencyVec &	getConvKernel () const override

const icarusutil::FrequencyVec &	getDeconvKernel () const override

double	getTOffset () const override

Private Types
using	IFieldResponsePtr = std::unique_ptr< icarus_tool::IFieldResponse >

using	IElectronicsResponsePtr = std::unique_ptr< icarus_tool::IElectronicsResponse >

using	IFilterPtr = std::unique_ptr< icarus_tool::IFilter >

Private Member Functions
void	calculateResponse (double sampling_rate, double weight)

std::string	numberToString (int number)

Private Member Functions inherited from icarus_tool::IResponse
virtual	~IResponse () noexcept=default

Private Attributes
bool	fResponseHasBeenSet

size_t	fThisPlane

double	f3DCorrection

double	fTimeScaleFactor

IFieldResponsePtr	fFieldResponse

IElectronicsResponsePtr	fElectronicsResponse

IFilterPtr	fFilter

size_t	fNumberTimeSamples

icarusutil::TimeVec	fResponse

icarusutil::FrequencyVec	fConvolutionKernel

icarusutil::FrequencyVec	fDeconvolutionKernel

bool	fUseEmpiricalOffsets
	Use emperical offsets divined from data. More...

double	fT0Offset
	The overall T0 offset for the response function. More...

std::unique_ptr < icarus_signal_processing::ICARUSFFT < double > >	fFFT
	Object to handle thread safe FFT. More...

Detailed Description

Definition at line 34 of file Response_tool.cc.

Member Typedef Documentation

using icarus_tool::Response::IElectronicsResponsePtr = std::unique_ptr<icarus_tool::IElectronicsResponse>

private

Definition at line 74 of file Response_tool.cc.

using icarus_tool::Response::IFieldResponsePtr = std::unique_ptr<icarus_tool::IFieldResponse>

private

Definition at line 73 of file Response_tool.cc.

using icarus_tool::Response::IFilterPtr = std::unique_ptr<icarus_tool::IFilter>

private

Definition at line 75 of file Response_tool.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 96 of file Response_tool.cc.

 {
     configure(pset);
 }

icarus_tool::Response::~Response ( )

inline

Definition at line 39 of file Response_tool.cc.

39 {}

Member Function Documentation

void icarus_tool::Response::calculateResponse	(	double	sampling_rate,
		double	weight
	)

private

Definition at line 183 of file Response_tool.cc.

 {
     // First of all set the field response
     fFieldResponse->setResponse(weight, f3DCorrection, fTimeScaleFactor);
 
     // What kind of response? (0= first induction, 1= middle induction, 2= collection)
     size_t responseType = fFieldResponse->getResponseType();
 
     // handle the electronics response for this plane
     fElectronicsResponse->setResponse(fFieldResponse->getResponseVec().size(), fFieldResponse->getBinWidth());
     
     // Next we perform the convolution of the field and electronics responses and then invert to get the
     // interim response function which will be set up in the time binning/range of the field response
     const icarusutil::FrequencyVec& fieldResponseFFTVec       = fFieldResponse->getResponseFFTVec();
     const icarusutil::FrequencyVec& electronicsResponseFFTVec = fElectronicsResponse->getResponseFFTVec();
     
     icarusutil::FrequencyVec curResponseFFTVec(fieldResponseFFTVec.size());
     
     std::transform(fieldResponseFFTVec.begin(), fieldResponseFFTVec.begin() + fieldResponseFFTVec.size()/2, electronicsResponseFFTVec.begin(), curResponseFFTVec.begin(), std::multiplies<icarusutil::ComplexVal>());
     
     // And now we recover the current response vector which is the convolution of the two above
     // (and still in the units of the original field response)
     icarusutil::TimeVec curResponseVec(fieldResponseFFTVec.size());
     
     // Note that we need a local version of the FFT because our time samples currently don't match what we will have
     icarus_signal_processing::ICARUSFFT<double> locFFT(curResponseVec.size());
 
     locFFT.inverseFFT(curResponseFFTVec, curResponseVec);
     
     // Now set to the task of determing the actual sampling response
     // We have to remember that the bin size for determining the field response probably
     // does not match that for the detector readout so we'll need to "convert"
     // from one to the other.
     fResponse.resize(fNumberTimeSamples, 0.);
     
     // Recover the combined response from above
 //    const std::vector<double>& curResponseVec = fSignalShapingICARUS.Response_save();
     
     double respIntegral = std::accumulate(curResponseVec.begin(),curResponseVec.end(),0.);
     
     mf::LogInfo("Response_tool") << "***** Response for plane: " << fThisPlane << " ******" << "\n"
                                  << "      initial response integral: " << respIntegral << std::endl;
     
     // Need two factors: 1) the detector sampling rate and 2) the response sampling rate
     double samplingRate = sampling_rate;                 // This is input in ns/tick
     double responseRate = fFieldResponse->getBinWidth(); // This is returned in ns
     double rateRatio    = samplingRate / responseRate;   // This gives the ratio of time bins for full readout to response bins
 
     // The idea is to step through each bin of the sampling response vector and then to
     // look up the corresponding bins in the current response vector. Since the two sample
     // rates are not the same there will be some "stretching" between the two. In addition,
     // we want to continue to allow for the possibility for further sample stretching
     double binScaleFactor = rateRatio * f3DCorrection * fTimeScaleFactor;
     
     // ok, do the loop
     for(size_t sampleIdx = 0; sampleIdx < fNumberTimeSamples; sampleIdx++)
     {
         // calculate the index for the response
         size_t responseLowIdx = std::floor(sampleIdx * binScaleFactor);
         
         if (responseLowIdx < curResponseVec.size())
         {
             // Calculate the index for the next bin
             size_t responseHiIdx = std::floor((sampleIdx + 1) * binScaleFactor);
             
             // This can't happen? But protect against zero divides...
             if (responseHiIdx == responseLowIdx) responseHiIdx += 1;
             
             if (responseHiIdx >= curResponseVec.size()) break;
             
             icarusutil::TimeVec::const_iterator curResponseItr = curResponseVec.begin();
             
             std::advance(curResponseItr,responseLowIdx);
             
             int nBins = responseHiIdx - responseLowIdx;
             
             // Obtain the average of these bins
             double aveResponse = std::accumulate(curResponseItr,curResponseItr+nBins,0.)/double(nBins);
             
             // Now interpolate between the two bins to get the sampling response for this bin
 //            double responseSlope = (curResponseVec.at(responseHiIdx) - curResponseVec.at(responseLowIdx)) / (responseHiIdx - responseLowIdx);
 //            double response      = curResponseVec.at(responseLowIdx) + 0.5 * responseSlope * (responseHiIdx - responseLowIdx);
 
             fResponse[sampleIdx] = aveResponse;
         }
     }
     
     // We need to scale by the binScaleFactor to preserve normalization
     std::transform(fResponse.begin(),fResponse.end(),fResponse.begin(),std::bind(std::multiplies<double>(),std::placeholders::_1,binScaleFactor));
     
     respIntegral = std::accumulate(fResponse.begin(),fResponse.end(),0.);
 
     // Now compute the T0 offset for the response function
     std::pair<icarusutil::TimeVec::iterator,icarusutil::TimeVec::iterator> minMaxPair = std::minmax_element(fResponse.begin(),fResponse.end());
 
     // Calculation of the T0 offset depends on the signal type
     int timeBin = std::distance(fResponse.begin(),minMaxPair.first);
 
     if (responseType == 2) timeBin = std::distance(fResponse.begin(),minMaxPair.second);
     else
     {
         std::cout << "*** plane: " << fThisPlane << ", first idx/val: " << std::distance(fResponse.begin(),minMaxPair.first) << "/" << *minMaxPair.first << ", second idx/val: " << std::distance(fResponse.begin(),minMaxPair.second) << "/" << *minMaxPair.second << std::endl;
         double weightSum = abs(*minMaxPair.first) + abs(*minMaxPair.second);
         double sum       = std::distance(fResponse.begin(),minMaxPair.first) * abs(*minMaxPair.first) + std::distance(fResponse.begin(),minMaxPair.second) * abs(*minMaxPair.second);
 
         timeBin = std::round(sum/weightSum);
     }
 
     std::cout << "*** Response tool is setting timeBin = " << timeBin << std::endl;
 
     if (fUseEmpiricalOffsets && fThisPlane == 0) timeBin = 83.;
     if (fUseEmpiricalOffsets && fThisPlane == 1) timeBin = 83.;
     
 //    // Do a backwards search to find the first positive bin
 //    while(1)
 //    {
 //        // Did we go too far?
 //        if (timeBin < 0)
 //            throw cet::exception("Response::configure") << "Cannot find zero-point crossover for induction response! ResponseType: " << responseType << ", plane: " << fThisPlane << std::endl;
 //            
 //        double content = fResponse[timeBin]; 
 //        
 //        if (content >= 0.) break;
 //        
 //        timeBin--;
 //    }
 //
 // 
     fT0Offset = -timeBin;     // Note that this value being returned is in tick units now
 
     std::cout << "** Response tool for plane: " << fThisPlane <<  ", T0: " << fT0Offset << std::endl;
     
     mf::LogInfo("Response_tool")  << "      final response integral: " << respIntegral << ", T0Offset: " << fT0Offset << std::endl;
 
     return;
 }

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

overridevirtual

Setup routines.

Parameters

output	the object containting the art output
clusHitPairVector	List of 3D hits to output as "extreme" space points

Implements icarus_tool::IResponse.

Definition at line 101 of file Response_tool.cc.

 {
     // Start by recovering the parameters
     fThisPlane           = pset.get<size_t>("Plane");
     f3DCorrection        = pset.get<size_t>("Correction3D");
     fTimeScaleFactor     = pset.get<size_t>("TimeScaleFactor");
     fUseEmpiricalOffsets = pset.get<bool  >("UseEmpiricalOffsets");
 
     fResponseHasBeenSet  = false;
 
     // Build out the underlying tools we'll be using
     fFieldResponse       = art::make_tool<icarus_tool::IFieldResponse>(pset.get<fhicl::ParameterSet>("FieldResponse"));
     fElectronicsResponse = art::make_tool<icarus_tool::IElectronicsResponse>(pset.get<fhicl::ParameterSet>("ElectronicsResponse"));
     fFilter              = art::make_tool<icarus_tool::IFilter>(pset.get<fhicl::ParameterSet>("Filter"));
 
     auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
     fNumberTimeSamples   = detProp.NumberTimeSamples();
 
     // Now set up our plans for doing the convolution
     fFFT = std::make_unique<icarus_signal_processing::ICARUSFFT<double>>(fNumberTimeSamples);
 
     return;
 }

const icarusutil::FrequencyVec& icarus_tool::Response::getConvKernel ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 53 of file Response_tool.cc.

53 {return fConvolutionKernel;}

icarus_tool::Response::fConvolutionKernel

icarusutil::FrequencyVec fConvolutionKernel

Definition: Response_tool.cc:85

const icarusutil::FrequencyVec& icarus_tool::Response::getDeconvKernel ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 54 of file Response_tool.cc.

54 {return fDeconvolutionKernel;}

icarus_tool::Response::fDeconvolutionKernel

icarusutil::FrequencyVec fDeconvolutionKernel

Definition: Response_tool.cc:86

const IElectronicsResponse* icarus_tool::Response::getElectronicsResponse ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 48 of file Response_tool.cc.

48 {return fElectronicsResponse.get();}

icarus_tool::Response::fElectronicsResponse

IElectronicsResponsePtr fElectronicsResponse

Definition: Response_tool.cc:79

const IFieldResponse* icarus_tool::Response::getFieldResponse ( ) const

inlineoverridevirtual

Recover the individual response elements.

Implements icarus_tool::IResponse.

Definition at line 47 of file Response_tool.cc.

47 {return fFieldResponse.get();}

icarus_tool::Response::fFieldResponse

IFieldResponsePtr fFieldResponse

Definition: Response_tool.cc:78

const IFilter* icarus_tool::Response::getFilter ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 49 of file Response_tool.cc.

49 {return fFilter.get();}

icarus_tool::Response::fFilter

IFilterPtr fFilter

Definition: Response_tool.cc:80

size_t icarus_tool::Response::getNumberTimeSamples ( ) const

inlineoverridevirtual

here recover the combined response elements

Parameters

output	the object containting the art output
clusHitPairVector	List of 3D hits to output as "extreme" space points

Implements icarus_tool::IResponse.

Definition at line 51 of file Response_tool.cc.

51 {return fNumberTimeSamples;}

icarus_tool::Response::fNumberTimeSamples

size_t fNumberTimeSamples

Definition: Response_tool.cc:83

size_t icarus_tool::Response::getPlane ( ) const

inlineoverridevirtual

Return the plane these functions represent.

Implements icarus_tool::IResponse.

Definition at line 45 of file Response_tool.cc.

45 {return fThisPlane;}

icarus_tool::Response::fThisPlane

size_t fThisPlane

Definition: Response_tool.cc:69

const icarusutil::TimeVec& icarus_tool::Response::getResponse ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 52 of file Response_tool.cc.

52 {return fResponse;}

icarus_tool::Response::fResponse

icarusutil::TimeVec fResponse

Definition: Response_tool.cc:84

double icarus_tool::Response::getTOffset ( ) const

inlineoverridevirtual

Implements icarus_tool::IResponse.

Definition at line 55 of file Response_tool.cc.

55 {return fT0Offset;};

icarus_tool::Response::fT0Offset

double fT0Offset

The overall T0 offset for the response function.

Definition: Response_tool.cc:89

std::string icarus_tool::Response::numberToString ( int number )

private

Definition at line 392 of file Response_tool.cc.

 {
     std::ostringstream string;
     
     string << std::setfill('0') << std::setw(2) << number;
     
     return string.str();
 }

void icarus_tool::Response::outputHistograms	(	double	sampling_rate,
		art::TFileDirectory &	histDir
	)		const

overridevirtual

Implements icarus_tool::IResponse.

Definition at line 321 of file Response_tool.cc.

 {
     // Create a subfolder in which to place the "response" histograms
     std::string thisResponse = "ResponsesPlane_" + std::to_string(fThisPlane);
     
     art::TFileDirectory dir = histDir.mkdir(thisResponse.c_str());
     
     // Do the field response histograms
     fFieldResponse->outputHistograms(dir);
     fElectronicsResponse->outputHistograms(dir);
     fFilter->outputHistograms(dir);
     
     // Now make hists for the full response
     std::string dirName = "Response_" + std::to_string(fThisPlane);
     
     art::TFileDirectory        responesDir  = dir.mkdir(dirName.c_str());
     const icarusutil::TimeVec& responseVec  = fResponse;
     double                     numBins      = responseVec.size();
     double                     samplingRate = sampling_rate; // Sampling time in us
     double                     maxFreq      = 1. / samplingRate;   // Max frequency in MHz
     double                     minFreq      = maxFreq / numBins;
     std::string                histName     = "Response_Plane_" + std::to_string(fThisPlane);
     TProfile*                  hist         = dir.make<TProfile>(histName.c_str(), "Response;Time(us)", numBins, 0., numBins * sampling_rate * 1.e-3);
 
     for(int bin = 0; bin < numBins; bin++)
     {
         hist->Fill((double(bin) + 0.5) * sampling_rate * 1.e-3, responseVec.at(bin), 1.);
     }
     
     icarusutil::TimeVec powerVec;
     
     fFFT->getFFTPower(responseVec, powerVec);
     
     std::string freqName  = "Response_FFTPlane_" + std::to_string(fThisPlane);
     TProfile*   freqHist  = dir.make<TProfile>(freqName.c_str(), "Response;Frequency(kHz)", numBins/2, minFreq, 0.5*maxFreq);
     
     for(size_t idx = 0; idx < numBins/2; idx++)
     {
         double freq = minFreq * (idx + 0.5);
         
         freqHist->Fill(freq, powerVec.at(idx), 1.);
     }
     
     const icarusutil::FrequencyVec& convKernel = fConvolutionKernel;
     
     std::string convKernelName   = "ConvKernel_" + std::to_string(fThisPlane);
     TProfile*   fullResponseHist = dir.make<TProfile>(convKernelName.c_str(), "Convolution Kernel;Frequency(kHz)", numBins/2, minFreq, 0.5*maxFreq);
 
     for(size_t idx = 0; idx < numBins/2; idx++)
     {
         double freq = minFreq * (idx + 0.5);
         
         fullResponseHist->Fill(freq, std::abs(convKernel[idx]), 1.);
     }
     
     const icarusutil::FrequencyVec& deconKernel = fDeconvolutionKernel;
     
     std::string deconName = "DeconKernel_" + std::to_string(fThisPlane);
     TProfile*   deconHist = dir.make<TProfile>(deconName.c_str(), "Deconvolution Kernel;Frequency(kHz)", numBins/2, minFreq, 0.5*maxFreq);
     
     for(size_t idx = 0; idx < numBins/2; idx++)
     {
         double freq = minFreq * (idx + 0.5);
         
         deconHist->Fill(freq, std::abs(deconKernel[idx]), 1.);
     }
 
     return;
 }

void icarus_tool::Response::setResponse	(	double	sampling_rate,
		double	weight
	)

overridevirtual

Implements icarus_tool::IResponse.

Definition at line 125 of file Response_tool.cc.

 {
     // If we have already done the setup then can return
     if (fResponseHasBeenSet) return;
 
     // Calculate the combined field and electronics shaping response
     calculateResponse(sampling_rate, weight);
 
     // Now we compute the convolution kernel which is a straigtforward operation
     fFFT->forwardFFT(fResponse, fConvolutionKernel);
 
     // Set up the filter for use in the deconvolution
     fFilter->setResponse(fNumberTimeSamples, f3DCorrection, fTimeScaleFactor);
 
     // Now compute the deconvolution kernel
     fDeconvolutionKernel = fFilter->getResponseVec();
 
     for(size_t idx = 0; idx < fNumberTimeSamples; idx++)
     {
         if (std::abs(fConvolutionKernel[idx]) < 0.0001) fDeconvolutionKernel[idx]  = 0.;
         else                                            fDeconvolutionKernel[idx] /= fConvolutionKernel[idx];
     }
     
     // The following can be uncommented to do some consistency checks if desired
 //    // Do some consistency/cross checks here
 //    // Check area of convolution function
 //    const std::vector<TComplex>& convKernel = fSignalShapingICARUS.ConvKernel();
 //
 //    double normFactor = std::accumulate(convKernel.begin(),convKernel.end(),0.,[](const auto& val, double sum){return sum + std::abs(val);});
 //
 //    mf::LogInfo("Response_tool")  << "Response for plane: " << fThisPlane << ", convKernel integral: " << normFactor << std::endl;
 //
 //    const std::vector<TComplex>& deconvKernel = fSignalShapingICARUS.DeconvKernel();
 //    std::vector<TComplex>  combKernel(deconvKernel.size());
 //
 //    std::transform(convKernel.begin(),convKernel.end(),deconvKernel.begin(),combKernel.begin(),std::multiplies<TComplex>());
 //
 //    const std::vector<TComplex>& filterKernel = fSignalShapingICARUS.Filter();
 //
 //    int    diffCount(0);
 //    double maxRhoDiff(0.);
 //
 //    for(size_t idx = 0; idx < filterKernel.size(); idx++)
 //    {
 //        double rhoDiff = filterKernel[idx].Rho() - combKernel[idx].Rho();
 //
 //        if (std::abs(rhoDiff) > 0.001) diffCount++;
 //
 //        if (std::abs(rhoDiff) > std::abs(maxRhoDiff)) maxRhoDiff = rhoDiff;
 //    }
 //
 //    mf::LogInfo("Response_tool") << "Checking recovery of the filter, # differences: " << diffCount << ", max diff seen: " << maxRhoDiff << std::endl;
 
     fResponseHasBeenSet = true;
 
     return;
 }

Member Data Documentation

double icarus_tool::Response::f3DCorrection

private

Definition at line 70 of file Response_tool.cc.

icarusutil::FrequencyVec icarus_tool::Response::fConvolutionKernel

private

Definition at line 85 of file Response_tool.cc.

icarusutil::FrequencyVec icarus_tool::Response::fDeconvolutionKernel

private

Definition at line 86 of file Response_tool.cc.

IElectronicsResponsePtr icarus_tool::Response::fElectronicsResponse

private

Definition at line 79 of file Response_tool.cc.

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

private

Object to handle thread safe FFT.

Definition at line 91 of file Response_tool.cc.

IFieldResponsePtr icarus_tool::Response::fFieldResponse

private

Definition at line 78 of file Response_tool.cc.

IFilterPtr icarus_tool::Response::fFilter

private

Definition at line 80 of file Response_tool.cc.

size_t icarus_tool::Response::fNumberTimeSamples

private

Definition at line 83 of file Response_tool.cc.

icarusutil::TimeVec icarus_tool::Response::fResponse

private

Definition at line 84 of file Response_tool.cc.

bool icarus_tool::Response::fResponseHasBeenSet

private

Definition at line 66 of file Response_tool.cc.

double icarus_tool::Response::fT0Offset

private

The overall T0 offset for the response function.

Definition at line 89 of file Response_tool.cc.

size_t icarus_tool::Response::fThisPlane

private

Definition at line 69 of file Response_tool.cc.

double icarus_tool::Response::fTimeScaleFactor

private

Definition at line 71 of file Response_tool.cc.

bool icarus_tool::Response::fUseEmpiricalOffsets

private

Use emperical offsets divined from data.

Definition at line 88 of file Response_tool.cc.

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

Response_tool.cc

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation