WireCell::Gen::EmpiricalNoiseModel Class Reference

#include <EmpiricalNoiseModel.h>

Inheritance diagram for WireCell::Gen::EmpiricalNoiseModel:

Classes
struct	NoiseSpectrum

Public Member Functions
	EmpiricalNoiseModel (const std::string &spectra_file="", const int nsamples=10000, const double period=0.5 *units::us, const double wire_length_scale=1.0 *units::cm, const std::string anode_tn="AnodePlane", const std::string chanstat_tn="StaticChannelStatus")
virtual	~EmpiricalNoiseModel ()
virtual const amplitude_t &	operator() (int chid) const
	IChannelSpectrum. More...
virtual const std::vector< float > &	freq () const
virtual const double	gain (int chid) const
virtual const double	shaping_time (int chid) const
virtual void	configure (const WireCell::Configuration &config)
	IConfigurable. More...
virtual WireCell::Configuration	default_configuration () const
	Optional, override to return a hard-coded default configuration. More...
void	resample (NoiseSpectrum &spectrum) const
void	gen_elec_resp_default ()
amplitude_t	interpolate (int plane, double wire_length) const
int	get_nsamples ()
Public Member Functions inherited from WireCell::IChannelSpectrum
virtual	~IChannelSpectrum ()
Public Member Functions inherited from WireCell::IComponent< IChannelSpectrum >
virtual	~IComponent ()
Public Member Functions inherited from WireCell::Interface
virtual	~Interface ()
Public Member Functions inherited from WireCell::IConfigurable
virtual	~IConfigurable ()
Public Member Functions inherited from WireCell::IComponent< IConfigurable >
virtual	~IComponent ()

Private Types
typedef std::unordered_map< int, amplitude_t >	len_amp_cache_t

Private Attributes
IAnodePlane::pointer	m_anode
IChannelStatus::pointer	m_chanstat
std::string	m_spectra_file
int	m_nsamples
int	m_fft_length
double	m_period
double	m_wlres
std::string	m_anode_tn
std::string	m_chanstat_tn
std::map< int, std::vector< NoiseSpectrum * > >	m_spectral_data
std::unordered_map< int, int >	m_chid_to_intlen
std::vector< len_amp_cache_t >	m_amp_cache
Waveform::realseq_t	m_elec_resp_freq
std::unordered_map< int, Waveform::realseq_t >	m_elec_resp_cache
amplitude_t	comb_amp
Log::logptr_t	log

Additional Inherited Members
Public Types inherited from WireCell::IChannelSpectrum
typedef std::vector< float >	amplitude_t
Public Types inherited from WireCell::IComponent< IChannelSpectrum >
typedef std::shared_ptr< IChannelSpectrum >	pointer
	Access subclass facet by pointer. More...
typedef std::vector< pointer >	vector
	Vector of shared pointers. More...
Public Types inherited from WireCell::Interface
typedef std::shared_ptr< Interface >	pointer
Public Types inherited from WireCell::IComponent< IConfigurable >
typedef std::shared_ptr< IConfigurable >	pointer
	Access subclass facet by pointer. More...
typedef std::vector< pointer >	vector
	Vector of shared pointers. More...

Detailed Description

Definition at line 31 of file EmpiricalNoiseModel.h.

Member Typedef Documentation

typedef std::unordered_map<int, amplitude_t> WireCell::Gen::EmpiricalNoiseModel::len_amp_cache_t

private

Definition at line 106 of file EmpiricalNoiseModel.h.

Constructor & Destructor Documentation

Gen::EmpiricalNoiseModel::EmpiricalNoiseModel	(	const std::string &	spectra_file = "",
		const int	nsamples = 10000,
		const double	period = 0.5*units::us,
		const double	wire_length_scale = 1.0*units::cm,
		const std::string	anode_tn = "AnodePlane",
		const std::string	chanstat_tn = "StaticChannelStatus"
	)

Definition at line 21 of file EmpiricalNoiseModel.cxx.

    : m_spectra_file(spectra_file)
    , m_nsamples(nsamples)
    , m_period(period)
    , m_wlres(wire_length_scale)
    // , m_tres(time_scale)
    // , m_gres(gain_scale)
    // , m_fres(freq_scale)
    , m_anode_tn(anode_tn)
    , m_chanstat_tn(chanstat_tn)
    , m_amp_cache(4)
    , log(Log::logger("sim"))
{
  m_fft_length = fft_best_length(m_nsamples);
  //  m_fft_length = m_nsamples;
  gen_elec_resp_default();
}

Gen::EmpiricalNoiseModel::~EmpiricalNoiseModel ( )

virtual

Definition at line 48 of file EmpiricalNoiseModel.cxx.

{
}

Member Function Documentation

void Gen::EmpiricalNoiseModel::configure ( const WireCell::Configuration &  config )

virtual

IConfigurable.

Implements WireCell::IConfigurable.

Definition at line 150 of file EmpiricalNoiseModel.cxx.

{
    m_anode_tn = get(cfg, "anode", m_anode_tn);
    m_anode = Factory::find_tn<IAnodePlane>(m_anode_tn);

    m_chanstat_tn = get(cfg, "chanstat", m_chanstat_tn);
    if (m_chanstat_tn != "") {  // allow for an empty channel status, no deviation from norm
        m_chanstat = Factory::find_tn<IChannelStatus>(m_chanstat_tn);
    }

    log->debug("EmpiricalNoiseModel: using anode {}, chanstat {}", m_anode_tn, m_chanstat_tn);

    m_spectra_file = get(cfg, "spectra_file", m_spectra_file);
        
    m_nsamples = get(cfg, "nsamples", m_nsamples);
    m_fft_length = fft_best_length(m_nsamples);
    //m_fft_length = m_nsamples;
    m_period = get(cfg, "period", m_period);
    m_wlres = get(cfg, "wire_length_scale", m_wlres);
    // m_tres = get(cfg, "time_scale", m_tres);
    // m_gres = get(cfg, "gain_scale", m_gres);
    // m_fres = get(cfg, "freq_scale", m_fres);

    // Load the data file assuming.  The file should be a list of
    // dictionaries matching the
    // wirecell.sigproc.noise.schema.NoiseSpectrum class.  Fixme:
    // should break out this code separate.
    auto jdat = Persist::load(m_spectra_file);
    const int nentries = jdat.size();
    m_spectral_data.clear();

    gen_elec_resp_default();
    
    for (int ientry=0; ientry<nentries; ++ientry) {
        auto jentry = jdat[ientry];
        NoiseSpectrum *nsptr = new NoiseSpectrum();

        nsptr->plane = jentry["plane"].asInt();
        nsptr->nsamples = jentry["nsamples"].asInt();
        nsptr->period = jentry["period"].asFloat();// * m_tres;  
        nsptr->gain = jentry["gain"].asFloat();// *m_gres;    
        nsptr->shaping = jentry["shaping"].asFloat();// * m_tres; 
        nsptr->wirelen = jentry["wirelen"].asFloat();// * m_wlres; 
        nsptr->constant = jentry["const"].asFloat();

        auto jfreqs = jentry["freqs"];
        const int nfreqs = jfreqs.size();
        nsptr->freqs.resize(nfreqs, 0.0);
        for (int ind=0; ind<nfreqs; ++ind) {
          nsptr->freqs[ind] = jfreqs[ind].asFloat();// * m_fres;
        }
        auto jamps = jentry["amps"];
        const int namps = jamps.size();
        nsptr->amps.resize(namps+1, 0.0);
        for (int ind=0; ind<namps; ++ind) {
            nsptr->amps[ind] = jamps[ind].asFloat();
        }
        // put the constant term at the end of the amplitude ... 
        // nsptr->amps[namps] = nsptr->constant;
        
        resample(*nsptr);
        m_spectral_data[nsptr->plane].push_back(nsptr); // assumes ordered by wire length!
    }        
}

WireCell::Configuration Gen::EmpiricalNoiseModel::default_configuration ( ) const

virtual

Optional, override to return a hard-coded default configuration.

The file holding the spectral data

Reimplemented from WireCell::IConfigurable.

Definition at line 80 of file EmpiricalNoiseModel.cxx.

{
    Configuration cfg;

    /// The file holding the spectral data
    cfg["spectra_file"] = m_spectra_file;
    cfg["nsamples"] = m_nsamples;          // number of samples up to Nyquist frequency
    cfg["period"] = m_period; 
    cfg["wire_length_scale"] = m_wlres;    // 
    // cfg["time_scale"] = m_tres;
    // cfg["gain_scale"] = m_gres;
    // cfg["freq_scale"] = m_fres;
    cfg["anode"] = m_anode_tn;            // name of IAnodePlane component

    return cfg;
}

const std::vector< float > & Gen::EmpiricalNoiseModel::freq ( ) const

virtual

Definition at line 266 of file EmpiricalNoiseModel.cxx.

{
  // assume frequency is universal ... 
  const int iplane = 0;
  auto it = m_spectral_data.find(iplane);
  const auto& spectra = it->second;
  return spectra.front()->freqs;
}

const double Gen::EmpiricalNoiseModel::gain ( int  chid ) const

virtual

Definition at line 284 of file EmpiricalNoiseModel.cxx.

{
  auto wpid = m_anode->resolve(chid);
  const int iplane = wpid.index();
  auto it = m_spectral_data.find(iplane);
  const auto& spectra = it->second;
  return spectra.front()->gain;
}

void Gen::EmpiricalNoiseModel::gen_elec_resp_default

(

)

Definition at line 52 of file EmpiricalNoiseModel.cxx.

                                                  {
  // double shaping[5]={1,1.1,2,2.2,3}; // us
  
  // calculate the frequencies ... 
  m_elec_resp_freq.resize(m_fft_length,0);
  for (unsigned int i=0;i!=m_elec_resp_freq.size();i++){
    if (i<=m_elec_resp_freq.size()/2.){
      m_elec_resp_freq.at(i) = i / (m_elec_resp_freq.size() *1.0) * 1./m_period ; // the second half is useless ... 
    }else{
      m_elec_resp_freq.at(i) = (m_elec_resp_freq.size()-i) / (m_elec_resp_freq.size() *1.0) * 1./m_period ; // the second half is useless ... 
    }
      
    // if (m_elec_resp_freq.at(i) > 1./m_period / 2.){
    //   m_elec_resp_freq.resize(i);
    //   break;
    // }
  }

  // for (int i=0;i!=5;i++){
  //   Response::ColdElec elec_resp(1, shaping[i]); // default at 1 mV/fC
  //   auto sig   =   elec_resp.generate(WireCell::Waveform::Domain(0, m_fft_length*m_period), m_fft_length);
  //   auto filt   = Waveform::dft(sig);
  //   int nconfig = shaping[i]/ 0.1;
  //   auto ele_resp_amp = Waveform::magnitude(filt);
  //   m_elec_resp_cache[nconfig] = ele_resp_amp;
  // }
}

int WireCell::Gen::EmpiricalNoiseModel::get_nsamples ( )

inline

Definition at line 87 of file EmpiricalNoiseModel.h.

{return m_nsamples;};

IChannelSpectrum::amplitude_t Gen::EmpiricalNoiseModel::interpolate	(	int	plane,
		double	wire_length
	)		const

Definition at line 216 of file EmpiricalNoiseModel.cxx.

{
    auto it = m_spectral_data.find(iplane);
    if (it == m_spectral_data.end()) {
        return amplitude_t();
    }

    const auto& spectra = it->second;

    

    // Linearly interpolate spectral data to wire length.

    // Any wire lengths
    // which are out of bounds causes the nearest result to be
    // returned (flat extrapolation)
    const NoiseSpectrum* front = spectra.front();
    if (wire_length <= front->wirelen) {
        return front->amps;
    }
    const NoiseSpectrum* back = spectra.back();
    if (wire_length >= back->wirelen) {
        return back->amps;
    }

    const int nspectra = spectra.size();
    for (int ind=1; ind<nspectra; ++ind) {
        const NoiseSpectrum* hi = spectra[ind];
        if (hi->wirelen < wire_length) {
            continue;
        }
        const NoiseSpectrum* lo = spectra[ind-1];

        const double delta = hi->wirelen - lo->wirelen;
        const double dist = wire_length - lo->wirelen;
        const double mu = dist/delta;

        const int nsamp = lo->amps.size();
        amplitude_t amp(nsamp);
        for (int isamp=0; isamp<nsamp; ++isamp) { // regular amplitude ... 
            amp[isamp] = lo->amps[isamp]*(1-mu) + hi->amps[isamp]*mu;
        }
        // amp[nsamp] = lo->constant*(1-mu) + hi->constant*mu; // do the constant term ...
        return amp;
    }
    // should not get here
    return amplitude_t();

}

const IChannelSpectrum::amplitude_t & Gen::EmpiricalNoiseModel::operator() ( int  chid ) const

virtual

IChannelSpectrum.

Implements WireCell::IChannelSpectrum.

Definition at line 295 of file EmpiricalNoiseModel.cxx.

{
    // get truncated wire length for cache
    int ilen=0;
    auto chlen = m_chid_to_intlen.find(chid);
    if (chlen == m_chid_to_intlen.end()) {  // new channel
        auto wires =  m_anode->wires(chid); // sum up wire length
        double len = 0.0;
        for (auto wire : wires) {
            len += ray_length(wire->ray());
        }
        // cache every cm ...
        ilen = int(len/m_wlres);
        // there might be  aproblem with the wire length's unit
        //ilen = int(len);
        m_chid_to_intlen[chid] = ilen;
    }
    else {
        ilen = chlen->second;
    }


    // saved content
    auto wpid = m_anode->resolve(chid);
    const int iplane = wpid.index();

    auto& amp_cache = m_amp_cache.at(iplane);
    auto lenamp = amp_cache.find(ilen);
    if (lenamp == amp_cache.end()) {
      auto amp = interpolate(iplane, ilen*m_wlres); // interpolate ... 
      amp_cache[ilen] = amp;
    }
    lenamp = amp_cache.find(ilen);
    
    // prepare to scale ... 
    double db_gain = gain(chid);
    double db_shaping = shaping_time(chid);

    // get channel gain
    //float ch_gain = 14 * units::mV/units::fC, ch_shaping = 2.0 * units::us; 
    
    double ch_gain = db_gain, ch_shaping = db_shaping;
    if (m_chanstat) {           // allow for deviation from nominal
        ch_gain = m_chanstat->preamp_gain(chid);
        ch_shaping = m_chanstat->preamp_shaping(chid);
    }
    double constant = lenamp->second.back();
    int nbin = lenamp->second.size()-1;

    //    amplitude_t comb_amp;
    comb_amp = lenamp->second;
    comb_amp.pop_back();
    
    if (fabs(ch_gain - db_gain ) > 0.01 * ch_gain){
      // scale the amplitude, not the constant term ...  
      for (int i=0;i!=nbin;i++){
        comb_amp.at(i) *= ch_gain/db_gain;
      }
    }
    
    if (fabs(ch_shaping-db_shaping) > 0.01*ch_shaping){
      // scale the amplitude by different shaping time ...   
      int nconfig = ch_shaping/units::us/0.1;
      auto resp1 = m_elec_resp_cache.find(nconfig);
      if (resp1 == m_elec_resp_cache.end()){

        Response::ColdElec elec_resp(10, ch_shaping); // default at 1 mV/fC
        auto sig   =   elec_resp.generate(WireCell::Waveform::Domain(0, m_fft_length*m_period), m_fft_length);
        auto filt   = Waveform::dft(sig);
        auto ele_resp_amp = Waveform::magnitude(filt);

        ele_resp_amp.resize(m_elec_resp_freq.size());
        m_elec_resp_cache[nconfig] = ele_resp_amp;

      }
      resp1 = m_elec_resp_cache.find(nconfig);
      
      nconfig = db_shaping/units::us/0.1;
      auto resp2 = m_elec_resp_cache.find(nconfig);
      if (resp2 == m_elec_resp_cache.end()){
        Response::ColdElec elec_resp(10, db_shaping); // default at 1 mV/fC
        auto sig   =   elec_resp.generate(WireCell::Waveform::Domain(0, m_fft_length*m_period), m_fft_length);
        auto filt   = Waveform::dft(sig);
        auto ele_resp_amp = Waveform::magnitude(filt);

        ele_resp_amp.resize(m_elec_resp_freq.size());
        m_elec_resp_cache[nconfig] = ele_resp_amp;
      }
      resp2 = m_elec_resp_cache.find(nconfig);
      
      for (int i=0;i!=nbin;i++){
        comb_amp.at(i) *= resp1->second.at(i)/resp2->second.at(i);
      }
      
    }
    
    
    // add the constant terms ... 
    for (int i=0;i!=nbin;i++){
      comb_amp.at(i) = sqrt(pow(comb_amp.at(i),2) + pow(constant,2)); // units still in mV
    }

    return comb_amp;
}

void Gen::EmpiricalNoiseModel::resample

(

NoiseSpectrum &

spectrum

)

const

Definition at line 97 of file EmpiricalNoiseModel.cxx.

{
  
    if (spectrum.nsamples == m_fft_length && spectrum.period == m_period) {
        return;
    }
    
    // scale the amplitude ... 
    double scale = sqrt(m_fft_length/(spectrum.nsamples*1.0) * spectrum.period / (m_period*1.0));
    spectrum.constant *= scale;
    for (unsigned int ind = 0; ind!=spectrum.amps.size(); ind++){
      spectrum.amps[ind] *= scale;
    }
    
    // interpolate ... 
    
    amplitude_t temp_amplitudes(m_fft_length,0);
    int count_low  = 0;
    int count_high = 1;
    double mu=0;
    for (int i=0;i!=m_fft_length;i++){
      double frequency = m_elec_resp_freq.at(i);
      if (frequency <= spectrum.freqs[0]){
        count_low = 0;
        count_high = 1;
        mu = 0;
      }else if (frequency >= spectrum.freqs.back()){
        count_low = spectrum.freqs.size()-2;
        count_high = spectrum.freqs.size()-1;
        mu = 1;
      }else{
        for (unsigned int j=0;j!=spectrum.freqs.size();j++){
          if (frequency>spectrum.freqs.at(j)){
            count_low = j;
            count_high = j+1;
          }else{
            break;
          }
        }
        mu = (frequency - spectrum.freqs.at(count_low)) / (spectrum.freqs.at(count_high)-spectrum.freqs.at(count_low));
      }
      
      
      temp_amplitudes.at(i) = (1-mu) * spectrum.amps[count_low] + mu * spectrum.amps[count_high];
      
    }

    spectrum.amps = temp_amplitudes;
    spectrum.amps.push_back(spectrum.constant);
    
    return;
}

const double Gen::EmpiricalNoiseModel::shaping_time ( int  chid ) const

virtual

Definition at line 275 of file EmpiricalNoiseModel.cxx.

{
  auto wpid = m_anode->resolve(chid);
  const int iplane = wpid.index();
  auto it = m_spectral_data.find(iplane);
  const auto& spectra = it->second;
  return spectra.front()->shaping;
}

Member Data Documentation

amplitude_t WireCell::Gen::EmpiricalNoiseModel::comb_amp

mutableprivate

Definition at line 112 of file EmpiricalNoiseModel.h.

Log::logptr_t WireCell::Gen::EmpiricalNoiseModel::log

private

Definition at line 114 of file EmpiricalNoiseModel.h.

std::vector<len_amp_cache_t> WireCell::Gen::EmpiricalNoiseModel::m_amp_cache

mutableprivate

Definition at line 107 of file EmpiricalNoiseModel.h.

IAnodePlane::pointer WireCell::Gen::EmpiricalNoiseModel::m_anode

private

Definition at line 87 of file EmpiricalNoiseModel.h.

std::string WireCell::Gen::EmpiricalNoiseModel::m_anode_tn

private

Definition at line 98 of file EmpiricalNoiseModel.h.

IChannelStatus::pointer WireCell::Gen::EmpiricalNoiseModel::m_chanstat

private

Definition at line 91 of file EmpiricalNoiseModel.h.

std::string WireCell::Gen::EmpiricalNoiseModel::m_chanstat_tn

private

Definition at line 98 of file EmpiricalNoiseModel.h.

std::unordered_map<int, int> WireCell::Gen::EmpiricalNoiseModel::m_chid_to_intlen

mutableprivate

Definition at line 104 of file EmpiricalNoiseModel.h.

std::unordered_map<int, Waveform::realseq_t> WireCell::Gen::EmpiricalNoiseModel::m_elec_resp_cache

mutableprivate

Definition at line 111 of file EmpiricalNoiseModel.h.

Waveform::realseq_t WireCell::Gen::EmpiricalNoiseModel::m_elec_resp_freq

private

Definition at line 110 of file EmpiricalNoiseModel.h.

int WireCell::Gen::EmpiricalNoiseModel::m_fft_length

private

Definition at line 95 of file EmpiricalNoiseModel.h.

int WireCell::Gen::EmpiricalNoiseModel::m_nsamples

private

Definition at line 94 of file EmpiricalNoiseModel.h.

double WireCell::Gen::EmpiricalNoiseModel::m_period

private

Definition at line 96 of file EmpiricalNoiseModel.h.

std::string WireCell::Gen::EmpiricalNoiseModel::m_spectra_file

private

Definition at line 93 of file EmpiricalNoiseModel.h.

std::map<int, std::vector<NoiseSpectrum*> > WireCell::Gen::EmpiricalNoiseModel::m_spectral_data

private

Definition at line 101 of file EmpiricalNoiseModel.h.

double WireCell::Gen::EmpiricalNoiseModel::m_wlres

private

Definition at line 96 of file EmpiricalNoiseModel.h.

---

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

• wire-cell-build/gen/inc/WireCellGen/EmpiricalNoiseModel.h
• wire-cell-build/gen/src/EmpiricalNoiseModel.cxx