Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASmoluchowskiDiffusion.hh
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
27  * G4DNASmoluchowskiDiffusion.hh
28  *
29  * Created on: 2 févr. 2015
30  * Author: matkara
31  */
32 
33 #ifndef G4DNASMOLUCHOWSKIDIFFUSION_HH_
34 #define G4DNASMOLUCHOWSKIDIFFUSION_HH_
35 
36 //#if __cplusplus >= 201103L
37 #include <cstdlib>
38 #include <cmath>
39 #include <vector>
40 #include <iostream>
41 #include <algorithm>
42 
43 //#define DNADEV_TEST
44 
45 #ifdef DNADEV_TEST
46 #include <TF1.h>
47 #endif
48 
49 #include <cassert>
50 
51 #ifndef DNADEV_TEST
52 #include "globals.hh"
53 #include "Randomize.hh"
54 #endif
55 
56 #ifdef DNADEV_TEST
57 #include "TRandom.h"
58 TRandom root_random;
59 double G4UniformRand()
60 {
61  return root_random.Rndm();
62 }
63 
64 #define G4cout std::cout
65 #define G4endl std::endl
66 #endif
67 
68 #include "G4Exp.hh"
69 
71 {
72 public:
73  G4DNASmoluchowskiDiffusion(double epsilon = 1e-5);
75 
76  static double ComputeS(double r, double D, double t)
77  {
78  double sTransform = r / (2. * std::sqrt(D * t));
79  return sTransform;
80  }
81 
82  static double ComputeDistance(double sTransform, double D, double t)
83  {
84  return sTransform * 2. * std::sqrt(D * t);
85  }
86 
87  static double ComputeTime(double sTransform, double D, double r)
88  {
89  return std::pow(r / sTransform, 2.) / (4. * D);
90  }
91 
92  //====================================================
93 
94  double GetRandomDistance(double _time, double D)
95  {
96  double proba = G4UniformRand();
97 // G4cout << "proba = " << proba << G4endl;
98  double sTransform = GetInverseProbability(proba);
99 // G4cout << "sTransform = " << sTransform << G4endl;
100  return ComputeDistance(sTransform, D, _time);
101  }
102 
103  double GetRandomTime(double distance, double D)
104  {
105  double proba = G4UniformRand();
106  double sTransform = GetInverseProbability(proba);
107  return ComputeTime(sTransform, D, distance);
108  }
109 
110  double EstimateCrossingTime(double proba,
111  double distance,
112  double D)
113  {
114  double sTransform = GetInverseProbability(proba);
115  return ComputeTime(sTransform, D, distance);
116  }
117 
118  //====================================================
119  // 1-value transformation
120 
121  // WARNING : this is NOT the differential probability
122  // this is the derivative of the function GetCumulativeProbability
123  static double GetDifferential(double sTransform)
124  {
125  static double constant = -4./std::sqrt(3.141592653589793);
126  return sTransform*sTransform*G4Exp(-sTransform*sTransform)*constant; // -4*sTransform*sTransform*exp(-sTransform*sTransform)/sqrt(3.141592653589793)
127  }
128 
129  static double GetDensityProbability(double r, double _time, double D)
130  {
131  static double my_pi = 3.141592653589793;
132  static double constant = 4.*my_pi/std::pow(4.*my_pi, 1.5);
133  return r*r/std::pow(D * _time,1.5)*G4Exp(-r*r/(4. * D * _time))*constant;
134  }
135 
136  //====================================================
137  // BOUNDING BOX
138  struct BoundingBox
139  {
140  double fXmax;
141  double fXmin;
142  double fXmaxDef;
143  double fXminDef;
144  double fToleranceY;
145  double fSum;
147 
149  {
153  };
154 
156 
157  BoundingBox(double xmin,
158  double xmax,
159  double toleranceY) :
160  fXmax(xmax), fXmin(xmin),
161  fToleranceY(toleranceY),
162  fSum(0)
163  {
164  if(fXmax < fXmin)
165  {
166  double tmp = fXmin;
167  fXmin = fXmax;
168  fXmax = tmp;
169  }
170 
171  fXminDef = fXmin;
172  fXmaxDef = fXmax;
175  }
176 
177  void Print()
178  {
179  G4cout << "fXmin: " << fXmin << " | fXmax: " << fXmax << G4endl;
180  }
181 
182  bool Propose(double proposedXValue,
183  double proposedProba,
184  double nextProba,
185  double& returnedValue)
186  {
187 // G4cout << "---------------------------" << G4endl;
188 // G4cout << "Proposed x value: " << proposedXValue
189 // << "| proposedProba: " << proposedProba
190 // << "| nextProba: " << nextProba
191 // << " | fXmin: " << fXmin << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmin) <<")"
192 // << " | fXmax: " << fXmax << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmax) <<")"
193 // << G4endl;
194 
195  bool returnFlag = false;
196 
197  if(proposedProba < nextProba-fToleranceY) // proba trop petite ==> augmente
198  {
199  // G4cout << "proposedProba < nextProba-fToleranceY" << G4endl;
200 
201  if(fIncreasingCumulativeFunction > 0) // croissant
202  {
203  if(proposedXValue > fXmin)
204  fXmin = proposedXValue;
205  }
206  else if(fIncreasingCumulativeFunction < 0) // decroissant
207  {
208  if(proposedXValue < fXmax)
209  fXmax = proposedXValue;
210  }
211 
212  returnedValue = (fXmax + fXmin)/2;
213  returnFlag = false;
215  }
216  else if(proposedProba > nextProba+fToleranceY) // proba trop grande
217  {
218  // G4cout << "proposedProba > nextProba+fToleranceY" << G4endl;
219 
221  {
222  if(proposedXValue < fXmax)
223  fXmax = proposedXValue;
224  }
226  {
227  if(proposedXValue > fXmin)
228  {
229  fXmin = proposedXValue;
230  }
231  }
232 
233  returnedValue = (fXmax + fXmin)/2;
234  returnFlag = false;
236  }
237  else
238  {
239  // G4cout << "IN THE INTERVAL !! : " << nextProba << G4endl;
240  fSum = proposedProba;
241 
242  // Assuming search for next proba is increasing
244  {
245  fXmin = fXminDef;
246  fXmax = proposedXValue;
247  }
249  {
250  fXmin = proposedXValue;
251  fXmax = fXmaxDef;
252  }
253  returnFlag = true;
255  }
256 
257  return returnFlag;
258  }
259  };
260  // END OF BOUNDING BOX
261  //==============================
262 
263  void PrepareReverseTable(double xmin, double xmax)
264  {
265  double x = xmax;
266  int index = 0;
267  double nextProba = fEpsilon;
268  double proposedX;
269 
270  BoundingBox boundingBox(xmin, xmax, fEpsilon*1e-5);
271 
272  while(index <= fNbins)
273  // in case GetCumulativeProbability is exact (digitally speaking), replace with:
274  // while(index <= fNbins+1)
275  {
276  nextProba = fEpsilon*index;
277 
278  double newProba = GetCumulativeProbability(x);
279 
280  if(boundingBox.Propose(x, newProba, nextProba, proposedX))
281  {
282  fInverse[index] = x;
283  index++;
284  }
285  else
286  {
287  if(x == proposedX)
288  {
289  G4cout << "BREAK : x= " << x << G4endl;
290  G4cout << "index= " << index << G4endl;
291  G4cout << "nextProba= " << nextProba << G4endl;
292  G4cout << "newProba= " << newProba << G4endl;
293  abort();
294  }
295  x = proposedX;
296  }
297  }
298 
299  fInverse[fNbins+1] = 0; // P(1) = 0, because we want it exact !
300  // Tips to improve the exactness: get an better value of pi, get better approximation of erf and exp, use long double instead of double
301 // boundingBox.Print();
302  }
303 
304  static double GetCumulativeProbability(double sTransform)
305  {
306  static double constant = 2./std::sqrt(3.141592653589793);
307  return erfc(sTransform) + constant*sTransform*G4Exp(-sTransform*sTransform);
308  }
309 
310  double GetInverseProbability(double proba) // returns sTransform
311  {
312  size_t index_low = (size_t) trunc(proba/fEpsilon);
313 
314  if(index_low == 0) // assymptote en 0
315  {
316  index_low = 1;
317  size_t index_up = 2;
318  double low_y = fInverse[index_low];
319  double up_y = fInverse[index_up];
320  double low_x = index_low*fEpsilon;
321  double up_x = proba+fEpsilon;
322  double tangente = (low_y-up_y)/(low_x - up_x); // ou utiliser GetDifferential(proba) ?
323  // double tangente = GetDifferential(proba);
324  return low_y + tangente*(proba-low_x);
325  }
326 
327  size_t index_up = index_low+1;
328  if(index_low > fInverse.size()) return fInverse.back();
329  double low_y = fInverse[index_low];
330  double up_y = fInverse[index_up];
331 
332  double low_x = index_low*fEpsilon;
333  double up_x = low_x+fEpsilon;
334 
335  if(up_x > 1) // P(1) = 0
336  {
337  up_x = 1;
338  up_y = 0; // more general : fInverse.back()
339  }
340 
341  double tangente = (low_y-up_y)/(low_x - up_x);
342 
343  return low_y + tangente*(proba-low_x);
344  }
345 
346  double PlotInverse(double* x, double* )
347  {
348  return GetInverseProbability(x[0]);
349  }
350 
351  double Plot(double* x, double* )
352  {
353  return GetDifferential(x[0]);
354  }
355 
356 
357  void InitialiseInverseProbability(double xmax = 3e28)
358  {
359  // x > x'
360  // P'(x) = p(x') = lim(x->x') (P(x') - P(x))/(x'-x)
361  // x'-x = (P(x') - P(x))/p(x')
362  // x = x' - (P(x') - P(x))/p(x')
363 
364  // fInverse initialized in the constructor
365 
366  assert(fNbins !=0);
367  PrepareReverseTable(0,xmax);
368  }
369 
370  std::vector<double> fInverse;
371  int fNbins;
372  double fEpsilon;
373 };
374 
375 #endif /* SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_MODELS_G4DNASMOLUCHOWSKIDIFFUSION_HH_ */
void InitialiseInverseProbability(double xmax=3e28)
static double ComputeDistance(double sTransform, double D, double t)
double EstimateCrossingTime(double proba, double distance, double D)
tuple x
Definition: test.py:50
double PlotInverse(double *x, double *)
static double GetCumulativeProbability(double sTransform)
void PrepareReverseTable(double xmin, double xmax)
G4DNASmoluchowskiDiffusion(double epsilon=1e-5)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
BoundingBox(double xmin, double xmax, double toleranceY)
double GetRandomTime(double distance, double D)
static double ComputeTime(double sTransform, double D, double r)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static double ComputeS(double r, double D, double t)
static double GetDifferential(double sTransform)
double D(double temp)
#define G4endl
Definition: G4ios.hh:61
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)