ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASmoluchowskiDiffusion.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNASmoluchowskiDiffusion.hh
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:
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 
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);
368  }
369 
370  std::vector<double> fInverse;
371  int fNbins;
372  double fEpsilon;
373 };
374 
375 #endif /* SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_MODELS_G4DNASMOLUCHOWSKIDIFFUSION_HH_ */