33 #ifndef G4DNASMOLUCHOWSKIDIFFUSION_HH_
34 #define G4DNASMOLUCHOWSKIDIFFUSION_HH_
61 return root_random.Rndm();
64 #define G4cout std::cout
65 #define G4endl std::endl
74 static double ComputeS(
double r,
double D,
double t)
76 double sTransform = r / (2. * std::sqrt(D * t));
82 return sTransform * 2. * std::sqrt(D * t);
87 return std::pow(r / sTransform, 2.) / (4. *
D);
123 static double constant = -4./std::sqrt(3.141592653589793);
124 return sTransform*sTransform*std::exp(-sTransform*sTransform)*constant;
129 static double my_pi = 3.141592653589793;
130 static double constant = 4.*my_pi/std::pow(4.*my_pi, 1.5);
131 return r*r/std::pow(D * _time,1.5)*std::exp(-r*r/(4. * D * _time))*constant;
158 fXmax(xmax), fXmin(xmin),
159 fToleranceY(toleranceY),
177 G4cout <<
"fXmin: " << fXmin <<
" | fXmax: " << fXmax <<
G4endl;
181 double proposedProba,
183 double& returnedValue)
193 bool returnFlag =
false;
195 if(proposedProba < nextProba-fToleranceY)
199 if(fIncreasingCumulativeFunction > 0)
201 if(proposedXValue > fXmin)
202 fXmin = proposedXValue;
204 else if(fIncreasingCumulativeFunction < 0)
206 if(proposedXValue < fXmax)
207 fXmax = proposedXValue;
210 returnedValue = (fXmax +
fXmin)/2;
214 else if(proposedProba > nextProba+fToleranceY)
218 if(fIncreasingCumulativeFunction>0)
220 if(proposedXValue < fXmax)
221 fXmax = proposedXValue;
223 else if(fIncreasingCumulativeFunction<0)
225 if(proposedXValue > fXmin)
227 fXmin = proposedXValue;
231 returnedValue = (fXmax +
fXmin)/2;
238 fSum = proposedProba;
241 if(fIncreasingCumulativeFunction<0)
244 fXmax = proposedXValue;
246 else if(fIncreasingCumulativeFunction>0)
248 fXmin = proposedXValue;
278 if(boundingBox.
Propose(x, newProba, nextProba, proposedX))
304 static double constant = 2./std::sqrt(3.141592653589793);
305 return erfc(sTransform) + constant*sTransform*std::exp(-sTransform*sTransform);
310 size_t index_low = (size_t) trunc(proba/
fEpsilon);
320 double tangente = (low_y-up_y)/(low_x - up_x);
322 return low_y + tangente*(proba-low_x);
325 size_t index_up = index_low+1;
339 double tangente = (low_y-up_y)/(low_x - up_x);
341 return low_y + tangente*(proba-low_x);
std::vector< double > fInverse
void InitialiseInverseProbability(double xmax=3e28)
static double ComputeDistance(double sTransform, double D, double t)
double EstimateCrossingTime(double proba, double distance, double D)
double PlotInverse(double *x, double *)
static double GetCumulativeProbability(double sTransform)
void PrepareReverseTable(double xmin, double xmax)
G4DNASmoluchowskiDiffusion(double epsilon=1e-5)
double GetInverseProbability(double proba)
double fIncreasingCumulativeFunction
G4GLOB_DLL std::ostream G4cout
BoundingBox(double xmin, double xmax, double toleranceY)
PreviousAction fPreviousAction
double GetRandomTime(double distance, double D)
static double ComputeTime(double sTransform, double D, double r)
static double ComputeS(double r, double D, double t)
const G4double x[NPOINTSGL]
static double GetDifferential(double sTransform)
static double GetDensityProbability(double r, double _time, double D)
bool Propose(double proposedXValue, double proposedProba, double nextProba, double &returnedValue)
double epsilon(double density, double temperature)
double Plot(double *x, double *)
double GetRandomDistance(double _time, double D)
virtual ~G4DNASmoluchowskiDiffusion()