#include <G4DNASmoluchowskiDiffusion.hh>

Classes
struct	BoundingBox

Public Member Functions
	G4DNASmoluchowskiDiffusion (double epsilon=1e-5)

virtual	~G4DNASmoluchowskiDiffusion ()

double	GetRandomDistance (double _time, double D)

double	GetRandomTime (double distance, double D)

double	EstimateCrossingTime (double proba, double distance, double D)

void	PrepareReverseTable (double xmin, double xmax)

double	GetInverseProbability (double proba)

double	PlotInverse (double x, double )

double	Plot (double x, double )

void	InitialiseInverseProbability (double xmax=3e28)

Static Public Member Functions
static double	ComputeS (double r, double D, double t)

static double	ComputeDistance (double sTransform, double D, double t)

static double	ComputeTime (double sTransform, double D, double r)

static double	GetDifferential (double sTransform)

static double	GetDensityProbability (double r, double _time, double D)

static double	GetCumulativeProbability (double sTransform)

Public Attributes
std::vector< double >	fInverse

int	fNbins

double	fEpsilon

Detailed Description

Definition at line 70 of file G4DNASmoluchowskiDiffusion.hh.

Constructor & Destructor Documentation

◆ G4DNASmoluchowskiDiffusion()

G4DNASmoluchowskiDiffusion::G4DNASmoluchowskiDiffusion ( double epsilon = 1e-5 )

Definition at line 50 of file G4DNASmoluchowskiDiffusion.cc.

                                                                     :  fEpsilon(epsilon)
{
  fNbins = (int) trunc(1/fEpsilon);
  // std::cout << "fNbins: " << fNbins << std::endl;
#ifdef DNADEV
  assert(fNbins > 0);
#endif
  fInverse.resize(fNbins+2); // trunc sous-estime + borne max a rajouter ==> 2
 
  // std::cout << "fInverse.capacity(): "<< fInverse.capacity() << std::endl;
}

◆ ~G4DNASmoluchowskiDiffusion()

G4DNASmoluchowskiDiffusion::~G4DNASmoluchowskiDiffusion ( )

virtualdefault

Member Function Documentation

◆ ComputeDistance()

static double G4DNASmoluchowskiDiffusion::ComputeDistance	(	double	sTransform,
		double	D,
		double	t )

inlinestatic

Definition at line 82 of file G4DNASmoluchowskiDiffusion.hh.

  {
    return sTransform * 2. * std::sqrt(D * t);
  }

Referenced by GetRandomDistance().

◆ ComputeS()

static double G4DNASmoluchowskiDiffusion::ComputeS	(	double	r,
		double	D,
		double	t )

inlinestatic

Definition at line 76 of file G4DNASmoluchowskiDiffusion.hh.

  {
    double sTransform = r / (2. * std::sqrt(D * t));
    return sTransform;
  }

◆ ComputeTime()

static double G4DNASmoluchowskiDiffusion::ComputeTime	(	double	sTransform,
		double	D,
		double	r )

inlinestatic

Definition at line 87 of file G4DNASmoluchowskiDiffusion.hh.

  {
    return std::pow(r / sTransform, 2.) / (4. * D);
  }

Referenced by EstimateCrossingTime(), and GetRandomTime().

◆ EstimateCrossingTime()

double G4DNASmoluchowskiDiffusion::EstimateCrossingTime	(	double	proba,
		double	distance,
		double	D )

inline

Definition at line 110 of file G4DNASmoluchowskiDiffusion.hh.

  {
    double sTransform = GetInverseProbability(proba);
    return ComputeTime(sTransform, D, distance);
  }

◆ GetCumulativeProbability()

static double G4DNASmoluchowskiDiffusion::GetCumulativeProbability ( double sTransform )

inlinestatic

Definition at line 303 of file G4DNASmoluchowskiDiffusion.hh.

  {
    static double constant = 2./std::sqrt(3.141592653589793);
    return erfc(sTransform) + constant*sTransform*G4Exp(-sTransform*sTransform);
  }

Referenced by G4DNASmoluchowskiDiffusion::BoundingBox::BoundingBox(), and PrepareReverseTable().

◆ GetDensityProbability()

static double G4DNASmoluchowskiDiffusion::GetDensityProbability	(	double	r,
		double	_time,
		double	D )

inlinestatic

Definition at line 129 of file G4DNASmoluchowskiDiffusion.hh.

  {
    static double my_pi = 3.141592653589793;
    static double constant = 4.*my_pi/std::pow(4.*my_pi, 1.5);
    return r*r/std::pow(D * _time,1.5)*G4Exp(-r*r/(4. * D * _time))*constant;
  }

◆ GetDifferential()

static double G4DNASmoluchowskiDiffusion::GetDifferential ( double sTransform )

inlinestatic

Definition at line 123 of file G4DNASmoluchowskiDiffusion.hh.

  {
    static double constant = -4./std::sqrt(3.141592653589793);
    return sTransform*sTransform*G4Exp(-sTransform*sTransform)*constant; // -4*sTransform*sTransform*exp(-sTransform*sTransform)/sqrt(3.141592653589793)
  }

Referenced by Plot().

◆ GetInverseProbability()

double G4DNASmoluchowskiDiffusion::GetInverseProbability ( double proba )

inline

Definition at line 309 of file G4DNASmoluchowskiDiffusion.hh.

  {
    auto  index_low = (size_t) trunc(proba/fEpsilon);
    
    if(index_low == 0) // assymptote en 0
    {
      index_low = 1;
      size_t index_up = 2;
      double low_y = fInverse[index_low];
      double up_y = fInverse[index_up];
      double low_x = index_low*fEpsilon;
      double up_x = proba+fEpsilon;
      double tangente = (low_y-up_y)/(low_x - up_x); // ou utiliser GetDifferential(proba) ?
      // double tangente = GetDifferential(proba);
      return low_y + tangente*(proba-low_x);
    }
 
    size_t index_up = index_low+1;
    if(index_low > fInverse.size()) return fInverse.back();
    double low_y = fInverse[index_low];
    double up_y = fInverse[index_up];
 
    double low_x = index_low*fEpsilon;
    double up_x = low_x+fEpsilon;
 
    if(up_x > 1) // P(1) = 0
    {
      up_x = 1;
      up_y = 0; // more general : fInverse.back()
    }
 
    double tangente = (low_y-up_y)/(low_x - up_x);
 
    return low_y + tangente*(proba-low_x);
  }

Referenced by EstimateCrossingTime(), GetRandomDistance(), GetRandomTime(), and PlotInverse().

◆ GetRandomDistance()

double G4DNASmoluchowskiDiffusion::GetRandomDistance	(	double	_time,
		double	D )

inline

Definition at line 94 of file G4DNASmoluchowskiDiffusion.hh.

  {
    double proba = G4UniformRand();
//    G4cout << "proba = " << proba << G4endl;
    double sTransform = GetInverseProbability(proba);
//    G4cout << "sTransform = " << sTransform << G4endl;
    return ComputeDistance(sTransform, D, _time);
  }

◆ GetRandomTime()

double G4DNASmoluchowskiDiffusion::GetRandomTime	(	double	distance,
		double	D )

inline

Definition at line 103 of file G4DNASmoluchowskiDiffusion.hh.

  {
    double proba = G4UniformRand();
    double sTransform = GetInverseProbability(proba);
    return ComputeTime(sTransform, D, distance);
  }

◆ InitialiseInverseProbability()

void G4DNASmoluchowskiDiffusion::InitialiseInverseProbability ( double xmax = 3e28 )

inline

Definition at line 356 of file G4DNASmoluchowskiDiffusion.hh.

  {
    // x > x'
    // P'(x) = p(x') = lim(x->x') (P(x') - P(x))/(x'-x)
    // x'-x = (P(x') - P(x))/p(x')
    // x = x' - (P(x') - P(x))/p(x')
 
    // fInverse initialized in the constructor
 
    assert(fNbins !=0);
    PrepareReverseTable(0,xmax);
  }

◆ Plot()

double G4DNASmoluchowskiDiffusion::Plot	(	double *	x,
		double *	)

inline

Definition at line 350 of file G4DNASmoluchowskiDiffusion.hh.

  {
    return GetDifferential(x[0]);
  }

◆ PlotInverse()

double G4DNASmoluchowskiDiffusion::PlotInverse	(	double *	x,
		double *	)

inline

Definition at line 345 of file G4DNASmoluchowskiDiffusion.hh.

  {
    return GetInverseProbability(x[0]);
  }

◆ PrepareReverseTable()

void G4DNASmoluchowskiDiffusion::PrepareReverseTable	(	double	xmin,
		double	xmax )

inline

Definition at line 262 of file G4DNASmoluchowskiDiffusion.hh.

  {
    double x = xmax;
    int index = 0;
    double nextProba = fEpsilon;
    double proposedX;
 
    BoundingBox boundingBox(xmin, xmax, fEpsilon*1e-5);
 
    while(index <= fNbins)
    // in case GetCumulativeProbability is exact (digitally speaking), replace with:
    // while(index <= fNbins+1)
    {
      nextProba = fEpsilon*index;
 
      double newProba = GetCumulativeProbability(x);
 
      if(boundingBox.Propose(x, newProba, nextProba, proposedX))
      {
        fInverse[index] = x;
        index++;
      }
      else
      {
        if(x == proposedX)
        {
          G4cout << "BREAK : x= " << x << G4endl;
          G4cout << "index= " << index << G4endl;
          G4cout << "nextProba= " << nextProba << G4endl;
          G4cout << "newProba= " << newProba << G4endl;
          abort();
        }
        x = proposedX;
      }
    }
    
    fInverse[fNbins+1] = 0; // P(1) = 0, because we want it exact !
    // Tips to improve the exactness: get an better value of pi, get better approximation of erf and exp, use long double instead of double
//    boundingBox.Print();
  }