WireCell::Gen::GaussianDiffusion Class Reference

#include <GaussianDiffusion.h>

Public Types
typedef std::shared_ptr< GaussianDiffusion >	pointer
typedef Array::array_xxf	patch_t

Public Member Functions
	GaussianDiffusion (const IDepo::pointer &depo, const GausDesc &time_desc, const GausDesc &pitch_desc)
	GaussianDiffusion. More...
void	set_sampling (const Binning &tbin, const Binning &pbin, double nsigma=3.0, IRandom::pointer fluctuate=nullptr, unsigned int weightstrat=1)
void	clear_sampling ()
const patch_t &	patch () const
const std::vector< double >	weights () const
int	toffset_bin () const
	Return the absolute time bin in the binning corresponding to column 0 of the patch. More...
int	poffset_bin () const
	Return the absolute impact bin in the binning corresponding to column 0 of the patch. More...
IDepo::pointer	depo () const
	Access deposition. More...
double	depo_time () const
double	depo_x () const
const GausDesc	pitch_desc ()
const GausDesc	time_desc ()

Private Attributes
IDepo::pointer	m_deposition
GausDesc	m_time_desc
GausDesc	m_pitch_desc
patch_t	m_patch
std::vector< double >	m_qweights
int	m_toffset_bin
int	m_poffset_bin

Detailed Description

Definition at line 68 of file GaussianDiffusion.h.

Member Typedef Documentation

typedef Array::array_xxf WireCell::Gen::GaussianDiffusion::patch_t

A patch is a 2D array of diffuse charge in (nimpacts X nticks) bins. patch[0] is drifted/diffused charge waveform at impact position 0 (relative to min pitch for the patch). patch[0][0] is charge at this impact position at time = 0 (relative to min time for the patch). See bin().

Definition at line 80 of file GaussianDiffusion.h.

typedef std::shared_ptr<GaussianDiffusion> WireCell::Gen::GaussianDiffusion::pointer

Definition at line 71 of file GaussianDiffusion.h.

Constructor & Destructor Documentation

Gen::GaussianDiffusion::GaussianDiffusion	(	const IDepo::pointer &	depo,
		const GausDesc &	time_desc,
		const GausDesc &	pitch_desc
	)

GaussianDiffusion.

Create a diffused deposition.

Definition at line 113 of file GaussianDiffusion.cxx.

    : m_deposition(depo)
    , m_time_desc(time_desc)
    , m_pitch_desc(pitch_desc)
    , m_toffset_bin(-1)
    , m_poffset_bin(-1)
{
}

Member Function Documentation

void Gen::GaussianDiffusion::clear_sampling

(

)

Definition at line 238 of file GaussianDiffusion.cxx.

                                         {
  m_patch.resize(0,0); 
  m_qweights.clear();
  m_qweights.shrink_to_fit();
}

IDepo::pointer WireCell::Gen::GaussianDiffusion::depo ( ) const

inline

Access deposition.

Definition at line 119 of file GaussianDiffusion.h.

{ return m_deposition; }

double WireCell::Gen::GaussianDiffusion::depo_time ( ) const

inline

Definition at line 121 of file GaussianDiffusion.h.

{ return m_deposition->time();}

double WireCell::Gen::GaussianDiffusion::depo_x ( ) const

inline

Definition at line 122 of file GaussianDiffusion.h.

{ return m_deposition->pos().x();}

const Gen::GaussianDiffusion::patch_t & Gen::GaussianDiffusion::patch

(

)

const

Get the diffusion patch as an array of N_pitch rows X N_time columns. Index as patch(i_pitch, i_time). Call set_sampling() first.

Definition at line 246 of file GaussianDiffusion.cxx.

{
    return m_patch;
}

const GausDesc WireCell::Gen::GaussianDiffusion::pitch_desc ( )

inline

Definition at line 124 of file GaussianDiffusion.h.

{ return m_pitch_desc; }

int WireCell::Gen::GaussianDiffusion::poffset_bin ( ) const

inline

Return the absolute impact bin in the binning corresponding to column 0 of the patch.

Definition at line 116 of file GaussianDiffusion.h.

{ return m_poffset_bin; }

void Gen::GaussianDiffusion::set_sampling	(	const Binning &	tbin,
		const Binning &	pbin,
		double	nsigma = 3.0,
		IRandom::pointer	fluctuate = nullptr,
		unsigned int	weightstrat = 1
	)

This fills the patch once matching the given time and pitch binning. The patch is limited to the 2D sample points that cover the subdomain determined by the number of sigma. If there should be Poisson fluctuations applied pass in an IRandom. Total charge is charge is preserved. Each cell of the patch represents the 2D bin-centered sampling of the Gaussian.

Sample time dimension

Sample pitch dimension.

fixme: for hanyu.

Definition at line 124 of file GaussianDiffusion.cxx.

{
    if (m_patch.size() > 0) {
        return;
    }

    /// Sample time dimension
    auto tval_range = m_time_desc.sigma_range(nsigma);
    auto tbin_range = tbin.sample_bin_range(tval_range.first, tval_range.second);
    const size_t ntss = tbin_range.second - tbin_range.first;
    m_toffset_bin = tbin_range.first;
    //auto tvec =  m_time_desc.sample(tbin.center(m_toffset_bin), tbin.binsize(), ntss);
    auto tvec =  m_time_desc.binint(tbin.edge(m_toffset_bin), tbin.binsize(), ntss);

    if (!ntss) {
        cerr << "Gen::GaussianDiffusion: no time bins for [" << tval_range.first/units::us << "," << tval_range.second/units::us << "] us\n";
        return;
    }

    /// Sample pitch dimension.
    auto pval_range = m_pitch_desc.sigma_range(nsigma);
    auto pbin_range = pbin.sample_bin_range(pval_range.first, pval_range.second);
    const size_t npss = pbin_range.second - pbin_range.first;
    m_poffset_bin = pbin_range.first;
    //auto pvec = m_pitch_desc.sample(pbin.center(m_poffset_bin), pbin.binsize(), npss);
    auto pvec = m_pitch_desc.binint(pbin.edge(m_poffset_bin), pbin.binsize(), npss);
    

    if (!npss) {
        cerr << "No impact bins [" << pval_range.first/units::mm << "," << pval_range.second/units::mm << "] mm\n";
        return;
    }

    // weightstrat = 1;

    // make charge weights for later interpolation.
    /// fixme: for hanyu.
    if(weightstrat == 2){
        auto wvec = m_pitch_desc.weight(pbin.edge(m_poffset_bin), pbin.binsize(), npss, pvec);
        m_qweights = wvec;
    }
    if(weightstrat == 1){
        m_qweights.resize(npss, 0.5);
    }

    // start making the time vs impact patch of charge.
    patch_t ret = patch_t::Zero(npss, ntss);
    double raw_sum=0.0;

    // Convolve the two independent Gaussians
    for (size_t ip = 0; ip < npss; ++ip) {
        for (size_t it = 0; it < ntss; ++it) {
            const double val = pvec[ip]*tvec[it];
            raw_sum += val;
            ret(ip,it) = (float)val;
        }
    }

    // normalize to total charge
    ret *= m_deposition->charge() / raw_sum;

    const double charge_sign = m_deposition->charge() < 0 ? -1 : 1;

    double fluc_sum = 0;
    if (fluctuate) {
        double unfluc_sum = 0;

        for (size_t ip = 0; ip < npss; ++ip) {
            for (size_t it = 0; it < ntss; ++it) {
                const float oldval = ret(ip,it);
                unfluc_sum += oldval;
                // should be a multinomial distribution, n_i follows binomial distribution
                // but n_i, n_j has covariance -n_tot * p_i * p_j
                // normalize later to approximate this multinomial distribution (how precise?)
                // how precise? better than poisson and 10000 total charge corresponds to a <1% level error.
                float number = fluctuate->binomial((int)(std::abs(m_deposition->charge())), oldval/m_deposition->charge());
                //the charge should be netagive -- ionization electrons
                fluc_sum += charge_sign*number;
                ret(ip,it) = charge_sign*number;
            }
        }
        if (fluc_sum == 0) {
           // cerr << "No net charge after fluctuation. Total unfluctuated = "
           //      << unfluc_sum
           //      << " Qdepo = " << m_deposition->charge()
           //      << endl;
           return;
        }
        else {
            ret *= m_deposition->charge() / fluc_sum;
        }
    }


    {                           // debugging
        //double retsum=0.0;
        //for (size_t ip = 0; ip < npss; ++ip) {
        //    for (size_t it = 0; it < ntss; ++it) {
        //        retsum += ret(ip,it);
        //    }
        //}
        // cerr << "GaussianDiffusion: Q in electrons: depo=" << m_deposition->charge()/units::eplus
        //      << " rawsum=" << raw_sum/units::eplus << " flucsum=" << fluc_sum/units::eplus
        //      << " returned=" << retsum/units::eplus << endl;
    }


    m_patch = ret;
}

const GausDesc WireCell::Gen::GaussianDiffusion::time_desc ( )

inline

Definition at line 125 of file GaussianDiffusion.h.

{ return m_time_desc; }

int WireCell::Gen::GaussianDiffusion::toffset_bin ( ) const

inline

Return the absolute time bin in the binning corresponding to column 0 of the patch.

Definition at line 113 of file GaussianDiffusion.h.

{ return m_toffset_bin; }

const std::vector< double > Gen::GaussianDiffusion::weights

(

)

const

Definition at line 251 of file GaussianDiffusion.cxx.

{
    return m_qweights;
}

Member Data Documentation

IDepo::pointer WireCell::Gen::GaussianDiffusion::m_deposition

private

Definition at line 129 of file GaussianDiffusion.h.

patch_t WireCell::Gen::GaussianDiffusion::m_patch

private

Definition at line 133 of file GaussianDiffusion.h.

GausDesc WireCell::Gen::GaussianDiffusion::m_pitch_desc

private

Definition at line 131 of file GaussianDiffusion.h.

int WireCell::Gen::GaussianDiffusion::m_poffset_bin

private

Definition at line 137 of file GaussianDiffusion.h.

std::vector<double> WireCell::Gen::GaussianDiffusion::m_qweights

private

Definition at line 134 of file GaussianDiffusion.h.

GausDesc WireCell::Gen::GaussianDiffusion::m_time_desc

private

Definition at line 131 of file GaussianDiffusion.h.

int WireCell::Gen::GaussianDiffusion::m_toffset_bin

private

Definition at line 136 of file GaussianDiffusion.h.

---

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

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