| Geant4 Cross Reference |
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
70 class G4DNASmoluchowskiDiffusion
71 {
72 public:
73 G4DNASmoluchowskiDiffusion(double epsilon = 1e-5);
74 virtual ~G4DNASmoluchowskiDiffusion();
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{0};
146 double fIncreasingCumulativeFunction;
147
148 enum PreviousAction
149 {
150 IncreaseProba,
151 DecreaseProba,
152 Undefined
153 };
154
155 PreviousAction fPreviousAction;
156
157 BoundingBox(double xmin,
158 double xmax,
159 double toleranceY) :
160 fXmax(xmax), fXmin(xmin),
161 fToleranceY(toleranceY)
162 {
163 if(fXmax < fXmin)
164 {
165 double tmp = fXmin;
166 fXmin = fXmax;
167 fXmax = tmp;
168 }
169
170 fXminDef = fXmin;
171 fXmaxDef = fXmax;
172 fPreviousAction = BoundingBox::Undefined;
173 fIncreasingCumulativeFunction = (GetCumulativeProbability(fXmax) - GetCumulativeProbability(fXmin))/(fXmax-fXmin);
174 }
175
176 void Print()
177 {
178 G4cout << "fXmin: " << fXmin << " | fXmax: " << fXmax << G4endl;
179 }
180
181 bool Propose(double proposedXValue,
182 double proposedProba,
183 double nextProba,
184 double& returnedValue)
185 {
186 // G4cout << "---------------------------" << G4endl;
187 // G4cout << "Proposed x value: " << proposedXValue
188 // << "| proposedProba: " << proposedProba
189 // << "| nextProba: " << nextProba
190 // << " | fXmin: " << fXmin << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmin) <<")"
191 // << " | fXmax: " << fXmax << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmax) <<")"
192 // << G4endl;
193
194 bool returnFlag = false;
195
196 if(proposedProba < nextProba-fToleranceY) // proba trop petite ==> augmente
197 {
198 // G4cout << "proposedProba < nextProba-fToleranceY" << G4endl;
199
200 if(fIncreasingCumulativeFunction > 0) // croissant
201 {
202 if(proposedXValue > fXmin)
203 fXmin = proposedXValue;
204 }
205 else if(fIncreasingCumulativeFunction < 0) // decroissant
206 {
207 if(proposedXValue < fXmax)
208 fXmax = proposedXValue;
209 }
210
211 returnedValue = (fXmax + fXmin)/2;
212 returnFlag = false;
213 fPreviousAction = BoundingBox::IncreaseProba;
214 }
215 else if(proposedProba > nextProba+fToleranceY) // proba trop grande
216 {
217 // G4cout << "proposedProba > nextProba+fToleranceY" << G4endl;
218
219 if(fIncreasingCumulativeFunction>0)
220 {
221 if(proposedXValue < fXmax)
222 fXmax = proposedXValue;
223 }
224 else if(fIncreasingCumulativeFunction<0)
225 {
226 if(proposedXValue > fXmin)
227 {
228 fXmin = proposedXValue;
229 }
230 }
231
232 returnedValue = (fXmax + fXmin)/2;
233 returnFlag = false;
234 fPreviousAction = BoundingBox::DecreaseProba;
235 }
236 else
237 {
238 // G4cout << "IN THE INTERVAL !! : " << nextProba << G4endl;
239 fSum = proposedProba;
240
241 // Assuming search for next proba is increasing
242 if(fIncreasingCumulativeFunction<0)
243 {
244 fXmin = fXminDef;
245 fXmax = proposedXValue;
246 }
247 else if(fIncreasingCumulativeFunction>0)
248 {
249 fXmin = proposedXValue;
250 fXmax = fXmaxDef;
251 }
252 returnFlag = true;
253 fPreviousAction = BoundingBox::Undefined;
254 }
255
256 return returnFlag;
257 }
258 };
259 // END OF BOUNDING BOX
260 //==============================
261
262 void PrepareReverseTable(double xmin, double xmax)
263 {
264 double x = xmax;
265 int index = 0;
266 double nextProba = fEpsilon;
267 double proposedX;
268
269 BoundingBox boundingBox(xmin, xmax, fEpsilon*1e-5);
270
271 while(index <= fNbins)
272 // in case GetCumulativeProbability is exact (digitally speaking), replace with:
273 // while(index <= fNbins+1)
274 {
275 nextProba = fEpsilon*index;
276
277 double newProba = GetCumulativeProbability(x);
278
279 if(boundingBox.Propose(x, newProba, nextProba, proposedX))
280 {
281 fInverse[index] = x;
282 index++;
283 }
284 else
285 {
286 if(x == proposedX)
287 {
288 G4cout << "BREAK : x= " << x << G4endl;
289 G4cout << "index= " << index << G4endl;
290 G4cout << "nextProba= " << nextProba << G4endl;
291 G4cout << "newProba= " << newProba << G4endl;
292 abort();
293 }
294 x = proposedX;
295 }
296 }
297
298 fInverse[fNbins+1] = 0; // P(1) = 0, because we want it exact !
299 // Tips to improve the exactness: get an better value of pi, get better approximation of erf and exp, use long double instead of double
300 // boundingBox.Print();
301 }
302
303 static double GetCumulativeProbability(double sTransform)
304 {
305 static double constant = 2./std::sqrt(3.141592653589793);
306 return erfc(sTransform) + constant*sTransform*G4Exp(-sTransform*sTransform);
307 }
308
309 double GetInverseProbability(double proba) // returns sTransform
310 {
311 auto index_low = (size_t) trunc(proba/fEpsilon);
312
313 if(index_low == 0) // assymptote en 0
314 {
315 index_low = 1;
316 size_t index_up = 2;
317 double low_y = fInverse[index_low];
318 double up_y = fInverse[index_up];
319 double low_x = index_low*fEpsilon;
320 double up_x = proba+fEpsilon;
321 double tangente = (low_y-up_y)/(low_x - up_x); // ou utiliser GetDifferential(proba) ?
322 // double tangente = GetDifferential(proba);
323 return low_y + tangente*(proba-low_x);
324 }
325
326 size_t index_up = index_low+1;
327 if(index_low > fInverse.size()) return fInverse.back();
328 double low_y = fInverse[index_low];
329 double up_y = fInverse[index_up];
330
331 double low_x = index_low*fEpsilon;
332 double up_x = low_x+fEpsilon;
333
334 if(up_x > 1) // P(1) = 0
335 {
336 up_x = 1;
337 up_y = 0; // more general : fInverse.back()
338 }
339
340 double tangente = (low_y-up_y)/(low_x - up_x);
341
342 return low_y + tangente*(proba-low_x);
343 }
344
345 double PlotInverse(double* x, double* )
346 {
347 return GetInverseProbability(x[0]);
348 }
349
350 double Plot(double* x, double* )
351 {
352 return GetDifferential(x[0]);
353 }
354
355
356 void InitialiseInverseProbability(double xmax = 3e28)
357 {
358 // x > x'
359 // P'(x) = p(x') = lim(x->x') (P(x') - P(x))/(x'-x)
360 // x'-x = (P(x') - P(x))/p(x')
361 // x = x' - (P(x') - P(x))/p(x')
362
363 // fInverse initialized in the constructor
364
365 assert(fNbins !=0);
366 PrepareReverseTable(0,xmax);
367 }
368
369 std::vector<double> fInverse;
370 int fNbins;
371 double fEpsilon;
372 };
373
374 #endif /* SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_MODELS_G4DNASMOLUCHOWSKIDIFFUSION_HH_ */
375