icarus_tool::CorrelatedNoise Class Reference

Inheritance diagram for icarus_tool::CorrelatedNoise:

Public Member Functions
	CorrelatedNoise (const fhicl::ParameterSet &pset)
	~CorrelatedNoise ()
void	configure (const fhicl::ParameterSet &pset) override
void	nextEvent () override
void	generateNoise (CLHEP::HepRandomEngine &noise_engine, CLHEP::HepRandomEngine &cornoise_engine, icarusutil::TimeVec &noise, detinfo::DetectorPropertiesData const &, double noise_factor, const geo::PlaneID &, unsigned int) override

Private Types
using	WaveformTools = icarus_signal_processing::WaveformTools< icarusutil::SigProcPrecision >

Private Member Functions
void	GenerateUncorrelatedNoise (CLHEP::HepRandomEngine &, icarusutil::TimeVec &, double)
void	GenerateCorrelatedNoise (CLHEP::HepRandomEngine &, icarusutil::TimeVec &, double, unsigned int)
void	GenNoise (std::function< void(double[])> &, const icarusutil::TimeVec &, icarusutil::TimeVec &, double)
void	ExtractCorrelatedAmplitude (float &, int) const
void	SelectContinuousSpectrum ()
void	FindPeaks ()
void	makeHistograms ()
Private Member Functions inherited from icarus_tool::IGenNoise
virtual	~IGenNoise () noexcept=default

Private Attributes
size_t	fPlane
int	fMedianNumBins
float	fNoiseRand
long	fCorrelatedSeed
long	fUncorrelatedSeed
float	fIncoherentNoiseFrac
bool	fStoreHistograms
std::string	fInputNoiseHistFileName
std::string	fHistogramName
std::string	fCorrAmpHistFileName
std::string	fCorrAmpHistogramName
WaveformTools	fWaveformTool
icarusutil::TimeVec	fNoiseHistVec
icarusutil::TimeVec	fCoherentNoiseVec
icarusutil::TimeVec	fIncoherentNoiseVec
icarusutil::TimeVec	fCorrAmpDistVec
icarusutil::SigProcPrecision	fIncoherentNoiseRMS
icarusutil::SigProcPrecision	fCoherentNoiseRMS
icarusutil::FrequencyVec	fNoiseFrequencyVec
bool	fNeedFirstSeed =true
TProfile *	fInputNoiseHist
TProfile *	fMedianNoiseHist
TProfile *	fPeakNoiseHist
TProfile *	fCorAmpDistHist
Eigen::FFT< double >	fEigenFFT

Detailed Description

Definition at line 40 of file CorrelatedNoise_tool.cc.

Member Typedef Documentation

using icarus_tool::CorrelatedNoise::WaveformTools = icarus_signal_processing::WaveformTools<icarusutil::SigProcPrecision>

private

Definition at line 81 of file CorrelatedNoise_tool.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 113 of file CorrelatedNoise_tool.cc.

{
    // Recover the configuration of the tool from the input fhicl file and set up
    configure(pset);
    
    // Now break out the coherent from the incoherent using the input overall spectrum
    SelectContinuousSpectrum();
    
    // Output some histograms to catalogue what's been done
    makeHistograms();
}

icarus_tool::CorrelatedNoise::~CorrelatedNoise

(

)

Definition at line 125 of file CorrelatedNoise_tool.cc.

{
}

Member Function Documentation

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

overridevirtual

Implements icarus_tool::IGenNoise.

Definition at line 129 of file CorrelatedNoise_tool.cc.

{
    // Recover the histogram used for noise generation
    fPlane                  = pset.get< size_t      >("Plane");
    fMedianNumBins          = pset.get< int         >("MedianNumBins");
    fNoiseRand              = pset.get< float       >("NoiseRand");
    fCorrelatedSeed         = pset.get< long        >("CorrelatedSeed",1000);
    fUncorrelatedSeed       = pset.get< long        >("UncorrelatedSeed",5000);
    fIncoherentNoiseFrac    = pset.get< float       >("IncoherentNoiseFraction",0.5);
    fStoreHistograms        = pset.get< bool        >("StoreHistograms");
    fInputNoiseHistFileName = pset.get< std::string >("NoiseHistFileName");
    fHistogramName          = pset.get< std::string >("HistogramName");
    fCorrAmpHistFileName    = pset.get< std::string >("CorrAmpHistFileName");
    fCorrAmpHistogramName   = pset.get< std::string >("CorrAmpHistogramName");
    
    // Initialize the work vector
    auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
    fNoiseFrequencyVec.resize(detProp.NumberTimeSamples(),std::complex<float>(0.,0.));

    // Set up to input the histogram with the overall noise spectrum
    std::string fullFileName;
    std::string corrAmpFileName;

    cet::search_path searchPath("FW_SEARCH_PATH");
    searchPath.find_file(fInputNoiseHistFileName, fullFileName);
    
    TFile inputFile(fullFileName.c_str(), "READ");
    
    if (!inputFile.IsOpen())
        throw cet::exception("NoiseFromHist::configure") << "Unable to open input file: " << fInputNoiseHistFileName << std::endl;
    
    TH1D* histPtr = (TH1D*)inputFile.Get(fHistogramName.c_str());
    
    if (!histPtr)
        throw cet::exception("NoiseFromHist::configure") << "Unable to recover desired histogram: " << fHistogramName << std::endl;
    
    // Close the input file
    inputFile.Close();
    
    // Now grab the histogram givng the observed "strenght" of the coherent noise
    searchPath.find_file(fCorrAmpHistFileName, corrAmpFileName);
    
    TFile corrAmpInputFile(corrAmpFileName.c_str(), "READ");
    
    if (!corrAmpInputFile.IsOpen())
        throw cet::exception("CorrelatedNoise::configure") << "Unable to open input file: " << fCorrAmpHistFileName << std::endl;
    
    TH1D* corrAmpHistPtr = (TH1D*)corrAmpInputFile.Get(fCorrAmpHistogramName.c_str());
    
    if (!corrAmpHistPtr)
        throw cet::exception("CorrelatedNoise::configure") << "Unable to recover desired histogram: " << fCorrAmpHistogramName << std::endl;
    
    // Close the input file
    corrAmpInputFile.Close();

    fNoiseHistVec.resize(histPtr->GetNbinsX(), 0.);
    fCorrAmpDistVec.resize(corrAmpHistPtr->GetNbinsX(),0.);
    
    // Recover the bin contents into local vectors
    for(size_t histIdx = 0; histIdx < size_t(histPtr->GetNbinsX()); histIdx++)
        fNoiseHistVec[histIdx] = histPtr->GetBinContent(histIdx+1);
    
    for(size_t histIdx = 0; histIdx < size_t(corrAmpHistPtr->GetNbinsX()); histIdx++)  // was 200
        fCorrAmpDistVec[histIdx] = corrAmpHistPtr->GetBinContent(histIdx+1);
    
    // Should we store hists?
    if (fStoreHistograms)
    {
        // Define histograms
        art::ServiceHandle<art::TFileService> tfs;
    
        art::TFileDirectory* histDirectory = tfs.get();
        
        // Make a directory for these histograms
        art::TFileDirectory dir = histDirectory->mkdir(Form("CorNoisePlane%1zu",fPlane));
        
        auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
        float sampleRate  = sampling_rate(clockData);
        float readOutSize = detProp.ReadOutWindowSize();
        float maxFreq     = 1.e6 / (2. * sampleRate);
        float minFreq     = 1.e6 / (2. * sampleRate * readOutSize);
        int   numSamples  = readOutSize / 2;
        
        fInputNoiseHist   = dir.make<TProfile>("InNoise",   ";freq(kHz)", numSamples, minFreq, maxFreq);
        fMedianNoiseHist  = dir.make<TProfile>("MedNoise",  ";freq(kHz)", numSamples, minFreq, maxFreq);;
        fPeakNoiseHist    = dir.make<TProfile>("PeakNoise", ";freq(kHz)", numSamples, minFreq, maxFreq);;
        
        fCorAmpDistHist   = dir.make<TProfile>("CorAmp",    ";Motherboard", fCorrAmpDistVec.size(),0.,fCorrAmpDistVec.size());
    }
    
    return;
}

void icarus_tool::CorrelatedNoise::ExtractCorrelatedAmplitude ( float &  corrFactor, int  board  ) const

private

Definition at line 430 of file CorrelatedNoise_tool.cc.

{
    CLHEP::RandGeneral amp_corr(fCorrAmpDistVec.data(),fCorrAmpDistVec.size(),0);
    amp_corr.setTheSeed(board+1);
    double rnd_corr[1] = {0.};
    
    amp_corr.fireArray(1,rnd_corr);
    
    float cfmedio=0.2287;
    corrFactor=rnd_corr[0]/cfmedio;
}

void icarus_tool::CorrelatedNoise::FindPeaks ( )

private

Definition at line 334 of file CorrelatedNoise_tool.cc.

{
    float dThr=50;
    std::vector<float> dHeight;
    dHeight.resize(fNoiseHistVec.size(), 0.);

    dHeight[0]=0;
    fCoherentNoiseVec.resize(fNoiseHistVec.size(),0);
    
    for(size_t histIdx = 1; histIdx < fNoiseHistVec.size(); histIdx++) {
        dHeight[histIdx]=fNoiseHistVec[histIdx]-fNoiseHistVec[histIdx-1];
    }
    
    for(size_t histIdx = 1; histIdx < fNoiseHistVec.size(); histIdx++)
        if(dHeight[histIdx]>dThr&&dHeight[histIdx+1]<dThr&&histIdx>5) fCoherentNoiseVec[histIdx]=1;
    
    return;
}

void icarus_tool::CorrelatedNoise::GenerateCorrelatedNoise ( CLHEP::HepRandomEngine &  engine, icarusutil::TimeVec &  noise, double  noise_factor, unsigned int  board  )

private

Definition at line 282 of file CorrelatedNoise_tool.cc.

{
    // The goal here is to produce a waveform with a coherent noise component
    // First check the extra scaling for a given motherboard
    float cf=1.;
    
    ExtractCorrelatedAmplitude(cf,board);
    
    // Only proceed if necessary
    if (cf > 0.)
    {
        // Set the engine seed to the board being considered
        engine.setSeed(fCorrelatedSeed+board,0);
        
        CLHEP::RandFlat noiseGen(engine,0,1);
        
        std::function<void (double[])> randGenFunc = [&noiseGen](double randArray[]){noiseGen.fireArray(2,randArray);};
        
        // Make the fraction the value that would happen if the quadrature sum of the two contributions equaled the input value
        float fraction    = std::sqrt(1. - fIncoherentNoiseFrac * fIncoherentNoiseFrac);
        float scaleFactor = fraction * cf * noise_factor / fCoherentNoiseRMS;
        
        GenNoise(randGenFunc, fCoherentNoiseVec, noise, scaleFactor);
    }
    
    return;
}

void icarus_tool::CorrelatedNoise::generateNoise ( CLHEP::HepRandomEngine &  noise_engine, CLHEP::HepRandomEngine &  cornoise_engine, icarusutil::TimeVec &  noise, detinfo::DetectorPropertiesData const &  , double  noise_factor, const geo::PlaneID &  planeID, unsigned int  board  )

overridevirtual

Implements icarus_tool::IGenNoise.

Definition at line 230 of file CorrelatedNoise_tool.cc.

{
    // Define a couple of vectors to hold intermediate work
    icarusutil::TimeVec noise_unc(noise.size(),0.);
    icarusutil::TimeVec noise_corr(noise.size(),0.);
    
    // Make sure the work vector is size right with the output
    if (fNoiseFrequencyVec.size() != noise.size()) fNoiseFrequencyVec.resize(noise.size(),std::complex<float>(0.,0.));
    
    // If applying incoherent noise call the generator
    if (fIncoherentNoiseFrac > 0.) GenerateUncorrelatedNoise(engine_unc,noise_unc,noise_factor);
    
//    int board=channel/32;

    // If applying coherent noise call the generator
    if (fIncoherentNoiseFrac < 1.) GenerateCorrelatedNoise(engine_corr, noise_corr, noise_factor, board);
    
    // Take the noise as the simple sum of the two contributions
    std::transform(noise_unc.begin(),noise_unc.end(),noise_corr.begin(),noise.begin(),std::plus<float>());
    
    return;
}

void icarus_tool::CorrelatedNoise::GenerateUncorrelatedNoise ( CLHEP::HepRandomEngine &  engine, icarusutil::TimeVec &  noise, double  noise_factor  )

private

Definition at line 259 of file CorrelatedNoise_tool.cc.

{
    // Here we aim to produce a waveform consisting of incoherent noise
    // Note that this is expected to be the dominate noise contribution
    // Check for seed initialization
    if (fNeedFirstSeed)
    {
        engine.setSeed(fUncorrelatedSeed,0);
        fNeedFirstSeed = false;
    }
    
    // Get the generator
    CLHEP::RandFlat noiseGen(engine,0,1);
    
    std::function<void (double[])> randGenFunc = [&noiseGen](double randArray[]){noiseGen.fireArray(2,randArray);};

    double scaleFactor = fIncoherentNoiseFrac * noise_factor / fIncoherentNoiseRMS;
    
    GenNoise(randGenFunc, fIncoherentNoiseVec, noise, scaleFactor);

    return;
}

void icarus_tool::CorrelatedNoise::GenNoise ( std::function< void(double[])> &  gen, const icarusutil::TimeVec &  freqDist, icarusutil::TimeVec &  noise, double  scaleFactor  )

private

Definition at line 310 of file CorrelatedNoise_tool.cc.

{
    double rnd_corr[2] = {0.,0.};
    
    // Build out the frequency vector
    for(size_t i=0; i< noise.size()/2; ++i)
    {
        // exponential noise spectrum
        gen(rnd_corr);
        
        double pval  = freqDist[i] * ((1-fNoiseRand) + 2 * fNoiseRand*rnd_corr[0]) * scaleFactor;
        double phase = rnd_corr[1] * 2. * M_PI;
        
        std::complex<double> tc(pval*cos(phase),pval*sin(phase));
        
        fNoiseFrequencyVec[i] = tc;
    }
    
    // inverse FFT MCSignal
    fEigenFFT.inv(noise, fNoiseFrequencyVec);
    
    return;
}

void icarus_tool::CorrelatedNoise::makeHistograms ( )

private

Definition at line 442 of file CorrelatedNoise_tool.cc.

{
    
    return;
}

void icarus_tool::CorrelatedNoise::nextEvent ( )

overridevirtual

Implements icarus_tool::IGenNoise.

Definition at line 222 of file CorrelatedNoise_tool.cc.

{
    // We update the correlated seed because we want to see different noise event-by-event
    fCorrelatedSeed   = (333 * fCorrelatedSeed) % 900000000;
    
    return;
}

void icarus_tool::CorrelatedNoise::SelectContinuousSpectrum ( )

private

Definition at line 353 of file CorrelatedNoise_tool.cc.

{
    // In the below we take the overall noise spectrum and perform a long bin smoothing to get its
    // general shape. We can then subtract this from the input spectrum to return a spectrum containing
    // the "spikes" which are associated to the coherent noise.
    // Generally, the inchorent noise model is then that with the spikes subtracted,
    // the coherent model is derived from the spikes...
    fIncoherentNoiseVec.resize(fNoiseHistVec.size(), 0.);
    
    // This does a median smoothing based on the number of bins we input
    fWaveformTool.medianSmooth(fNoiseHistVec, fIncoherentNoiseVec, fMedianNumBins);
    
    icarusutil::TimeVec peakVec(fNoiseHistVec.size());
    
    // Subtract the smoothed from the input to get the peaks
    std::transform(fNoiseHistVec.begin(),fNoiseHistVec.end(),fIncoherentNoiseVec.begin(),peakVec.begin(),std::minus<float>());
    
    // Try to clean up the "peak" vector a bit
    icarusutil::TimeVec::iterator maxItr = std::max_element(peakVec.begin(),peakVec.end());
    
    float maxPeak   = *maxItr;
    float threshold = std::min(0.1 * maxPeak, 20.);

    fCoherentNoiseVec.clear();
    fCoherentNoiseVec.resize(peakVec.size(),0.);
    
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
    auto const samplingRate = sampling_rate(clockData);
    for(size_t idx = 0; idx < peakVec.size(); idx++)
    {
        if (peakVec[idx] > threshold) fCoherentNoiseVec[idx] = peakVec[idx];
        
        if (fStoreHistograms)
        {
            float freq = 1.e6 * float(idx)/ (2. * samplingRate * fCoherentNoiseVec.size());
            
            fInputNoiseHist->Fill(freq,fNoiseHistVec.at(idx),1.);
            fMedianNoiseHist->Fill(freq,fIncoherentNoiseVec.at(idx),1.);
            fPeakNoiseHist->Fill(freq,peakVec.at(idx),1.);
        }
    }
    
    // Let's get the rms we expect from the incoherent noise contribution to the waveform
    // A couple of ways to do this, let's basically invert the frequency spectrum to
    // produce a waveform and then get the rms from that
    std::function<void (double[])> randGenFunc = [](double randArray[]){randArray[0]=0.5; randArray[1]=0.5;};
    
    icarusutil::TimeVec waveNoise(fIncoherentNoiseVec.size());
    float               scaleFactor = 1.;
    
    GenNoise(randGenFunc, fIncoherentNoiseVec, waveNoise, scaleFactor);
    
    // Now get the details...
    icarusutil::SigProcPrecision nSig(3.);
    icarusutil::SigProcPrecision mean,rmsTrunc;
    int                          nTrunc;
    int                          range;
    
    // Use the waveform tool to recover the full rms
    fWaveformTool.getTruncatedMean(waveNoise, mean, nTrunc, range);
    fWaveformTool.getTruncatedRMS(waveNoise, nSig, fIncoherentNoiseRMS, rmsTrunc, nTrunc);
    
    // Do the same for the coherent term
    GenNoise(randGenFunc, fCoherentNoiseVec, waveNoise, scaleFactor);
    
    fWaveformTool.getTruncatedMean(waveNoise, mean, nTrunc, range);
    fWaveformTool.getTruncatedRMS(waveNoise, nSig, fCoherentNoiseRMS, rmsTrunc, nTrunc);

    if (fStoreHistograms)
    {
        for(size_t idx = 0; idx < fCorrAmpDistVec.size(); idx++)
            fCorAmpDistHist->Fill(idx, fCorrAmpDistVec[idx], 1.);
    }

    return;
}

Member Data Documentation

icarusutil::SigProcPrecision icarus_tool::CorrelatedNoise::fCoherentNoiseRMS

private

Definition at line 93 of file CorrelatedNoise_tool.cc.

icarusutil::TimeVec icarus_tool::CorrelatedNoise::fCoherentNoiseVec

private

Definition at line 88 of file CorrelatedNoise_tool.cc.

TProfile* icarus_tool::CorrelatedNoise::fCorAmpDistHist

private

Definition at line 105 of file CorrelatedNoise_tool.cc.

icarusutil::TimeVec icarus_tool::CorrelatedNoise::fCorrAmpDistVec

private

Definition at line 90 of file CorrelatedNoise_tool.cc.

std::string icarus_tool::CorrelatedNoise::fCorrAmpHistFileName

private

Definition at line 78 of file CorrelatedNoise_tool.cc.

std::string icarus_tool::CorrelatedNoise::fCorrAmpHistogramName

private

Definition at line 79 of file CorrelatedNoise_tool.cc.

long icarus_tool::CorrelatedNoise::fCorrelatedSeed

private

Definition at line 72 of file CorrelatedNoise_tool.cc.

Eigen::FFT<double> icarus_tool::CorrelatedNoise::fEigenFFT

private

Definition at line 108 of file CorrelatedNoise_tool.cc.

std::string icarus_tool::CorrelatedNoise::fHistogramName

private

Definition at line 77 of file CorrelatedNoise_tool.cc.

float icarus_tool::CorrelatedNoise::fIncoherentNoiseFrac

private

Definition at line 74 of file CorrelatedNoise_tool.cc.

icarusutil::SigProcPrecision icarus_tool::CorrelatedNoise::fIncoherentNoiseRMS

private

Definition at line 92 of file CorrelatedNoise_tool.cc.

icarusutil::TimeVec icarus_tool::CorrelatedNoise::fIncoherentNoiseVec

private

Definition at line 89 of file CorrelatedNoise_tool.cc.

TProfile* icarus_tool::CorrelatedNoise::fInputNoiseHist

private

Definition at line 102 of file CorrelatedNoise_tool.cc.

std::string icarus_tool::CorrelatedNoise::fInputNoiseHistFileName

private

Definition at line 76 of file CorrelatedNoise_tool.cc.

TProfile* icarus_tool::CorrelatedNoise::fMedianNoiseHist

private

Definition at line 103 of file CorrelatedNoise_tool.cc.

int icarus_tool::CorrelatedNoise::fMedianNumBins

private

Definition at line 70 of file CorrelatedNoise_tool.cc.

bool icarus_tool::CorrelatedNoise::fNeedFirstSeed =true

private

Definition at line 99 of file CorrelatedNoise_tool.cc.

icarusutil::FrequencyVec icarus_tool::CorrelatedNoise::fNoiseFrequencyVec

private

Definition at line 96 of file CorrelatedNoise_tool.cc.

icarusutil::TimeVec icarus_tool::CorrelatedNoise::fNoiseHistVec

private

Definition at line 87 of file CorrelatedNoise_tool.cc.

float icarus_tool::CorrelatedNoise::fNoiseRand

private

Definition at line 71 of file CorrelatedNoise_tool.cc.

TProfile* icarus_tool::CorrelatedNoise::fPeakNoiseHist

private

Definition at line 104 of file CorrelatedNoise_tool.cc.

size_t icarus_tool::CorrelatedNoise::fPlane

private

Definition at line 69 of file CorrelatedNoise_tool.cc.

bool icarus_tool::CorrelatedNoise::fStoreHistograms

private

Definition at line 75 of file CorrelatedNoise_tool.cc.

long icarus_tool::CorrelatedNoise::fUncorrelatedSeed

private

Definition at line 73 of file CorrelatedNoise_tool.cc.

WaveformTools icarus_tool::CorrelatedNoise::fWaveformTool

private

Definition at line 83 of file CorrelatedNoise_tool.cc.

---

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

• CorrelatedNoise_tool.cc