EIC code displayed by LXR master/​include/​Geant4/​G4DNASmoluchowskiDiffusion.hh	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 09:58:09

 //
 // ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
 // * conditions of the Geant4 Software License,  included in the file *
 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
 // * include a list of copyright holders.                             *
 // *                                                                  *
 // * Neither the authors of this software system, nor their employing *
 // * institutes,nor the agencies providing financial support for this *
 // * work  make  any representation or  warranty, express or implied, *
 // * regarding  this  software system or assume any liability for its *
 // * use.  Please see the license in the file  LICENSE  and URL above *
 // * for the full disclaimer and the limitation of liability.         *
 // *                                                                  *
 // * This  code  implementation is the result of  the  scientific and *
 // * technical work of the GEANT4 collaboration.                      *
 // * By using,  copying,  modifying or  distributing the software (or *
 // * any work based  on the software)  you  agree  to acknowledge its *
 // * use  in  resulting  scientific  publications,  and indicate your *
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
 //
 /*
  * G4DNASmoluchowskiDiffusion.hh
  *
  *  Created on: 2 févr. 2015
  *      Author: matkara
  */
 
 #ifndef G4DNASMOLUCHOWSKIDIFFUSION_HH_
 #define G4DNASMOLUCHOWSKIDIFFUSION_HH_
 
 //#if __cplusplus >= 201103L
 #include <cstdlib>
 #include <cmath>
 #include <vector>
 #include <iostream>
 #include <algorithm>
 
 //#define DNADEV_TEST
 
 #ifdef DNADEV_TEST
 #include <TF1.h>
 #endif
 
 #include <cassert>
 
 #ifndef DNADEV_TEST
 #include "globals.hh"
 #include "Randomize.hh"
 #endif
 
 #ifdef DNADEV_TEST
 #include "TRandom.h"
 TRandom root_random;
 double G4UniformRand()
 {
   return root_random.Rndm();
 }
 
 #define G4cout std::cout
 #define G4endl std::endl
 #endif
 
 #include "G4Exp.hh"
 
 class G4DNASmoluchowskiDiffusion
 {
 public:
   G4DNASmoluchowskiDiffusion(double epsilon = 1e-5);
   virtual ~G4DNASmoluchowskiDiffusion();
 
   static double ComputeS(double r, double D, double t)
   {
     double sTransform = r / (2. * std::sqrt(D * t));
     return sTransform;
   }
 
   static double ComputeDistance(double sTransform, double D, double t)
   {
     return sTransform * 2. * std::sqrt(D * t);
   }
 
   static double ComputeTime(double sTransform, double D, double r)
   {
     return std::pow(r / sTransform, 2.) / (4. * D);
   }
 
   //====================================================
 
   double GetRandomDistance(double _time, double D)
   {
     double proba = G4UniformRand();
 //    G4cout << "proba = " << proba << G4endl;
     double sTransform = GetInverseProbability(proba);
 //    G4cout << "sTransform = " << sTransform << G4endl;
     return ComputeDistance(sTransform, D, _time);
   }
 
   double GetRandomTime(double distance, double D)
   {
     double proba = G4UniformRand();
     double sTransform = GetInverseProbability(proba);
     return ComputeTime(sTransform, D, distance);
   }
 
   double EstimateCrossingTime(double proba,
                               double distance,
                               double D)
   {
     double sTransform = GetInverseProbability(proba);
     return ComputeTime(sTransform, D, distance);
   }
 
   //====================================================
   // 1-value transformation
 
   // WARNING : this is NOT the differential probability
   // this is the derivative of the function GetCumulativeProbability
   static double GetDifferential(double sTransform)
   {
     static double constant = -4./std::sqrt(3.141592653589793);
     return sTransform*sTransform*G4Exp(-sTransform*sTransform)*constant; // -4*sTransform*sTransform*exp(-sTransform*sTransform)/sqrt(3.141592653589793)
   }
 
   static double GetDensityProbability(double r, double _time, double D)
   {
     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;
   }
 
   //====================================================
   // BOUNDING BOX
   struct BoundingBox
   {
     double fXmax;
     double fXmin;
     double fXmaxDef;
     double fXminDef;
     double fToleranceY;
     double fSum{0};
     double    fIncreasingCumulativeFunction;
 
     enum PreviousAction
     {
       IncreaseProba,
       DecreaseProba,
       Undefined
     };
 
     PreviousAction fPreviousAction;
 
     BoundingBox(double xmin,
                 double xmax,
                 double toleranceY) :
      fXmax(xmax), fXmin(xmin),
      fToleranceY(toleranceY) 
     {
       if(fXmax < fXmin)
       {
         double tmp = fXmin;
         fXmin = fXmax;
         fXmax = tmp;
       }
       
       fXminDef = fXmin;
       fXmaxDef = fXmax;
       fPreviousAction = BoundingBox::Undefined;
       fIncreasingCumulativeFunction = (GetCumulativeProbability(fXmax) - GetCumulativeProbability(fXmin))/(fXmax-fXmin);
     }
     
     void Print()
     {
       G4cout << "fXmin: " << fXmin << " | fXmax: " << fXmax << G4endl;
     }
 
     bool Propose(double proposedXValue,
                  double proposedProba,
                  double nextProba,
                  double& returnedValue)
     {
 //      G4cout << "---------------------------" << G4endl;
 //      G4cout << "Proposed x value: " << proposedXValue
 //          << "| proposedProba: " << proposedProba
 //          << "| nextProba: " << nextProba
 //          << " | fXmin: " << fXmin << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmin) <<")"
 //          << " | fXmax: " << fXmax << " (" << G4DNASmoluchowskiDiffusion::GetCumulativeProbability(fXmax) <<")"
 //          << G4endl;
 
       bool returnFlag = false;
       
       if(proposedProba < nextProba-fToleranceY) // proba trop petite ==> augmente
       {
         // G4cout << "proposedProba < nextProba-fToleranceY" << G4endl;
 
         if(fIncreasingCumulativeFunction > 0) // croissant
         {
           if(proposedXValue > fXmin)
             fXmin = proposedXValue;
         }
         else if(fIncreasingCumulativeFunction < 0) // decroissant
         {
           if(proposedXValue < fXmax)
             fXmax = proposedXValue;
         }
         
         returnedValue = (fXmax + fXmin)/2;
         returnFlag = false;
         fPreviousAction = BoundingBox::IncreaseProba;
       }
       else if(proposedProba > nextProba+fToleranceY) // proba trop grande
       {
         // G4cout << "proposedProba > nextProba+fToleranceY" << G4endl;
 
         if(fIncreasingCumulativeFunction>0)
         {
           if(proposedXValue < fXmax)
             fXmax = proposedXValue;
         }
         else if(fIncreasingCumulativeFunction<0)
         {
           if(proposedXValue > fXmin)
           {
             fXmin = proposedXValue;
           }
         }
         
         returnedValue = (fXmax + fXmin)/2;
         returnFlag = false;
         fPreviousAction = BoundingBox::DecreaseProba;
       }
       else
       {
         // G4cout << "IN THE INTERVAL !! : " << nextProba << G4endl;
         fSum = proposedProba;
         
         // Assuming search for next proba is increasing
         if(fIncreasingCumulativeFunction<0)
         {
          fXmin = fXminDef;
          fXmax = proposedXValue;
         }
         else if(fIncreasingCumulativeFunction>0)
         {
           fXmin = proposedXValue;
           fXmax = fXmaxDef;
         }
         returnFlag = true;
         fPreviousAction = BoundingBox::Undefined;
       }
       
       return returnFlag;
     }
   };
   // END OF BOUNDING BOX
   //==============================
   
   void PrepareReverseTable(double xmin, double xmax)
   {
     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();
   }
 
   static double GetCumulativeProbability(double sTransform)
   {
     static double constant = 2./std::sqrt(3.141592653589793);
     return erfc(sTransform) + constant*sTransform*G4Exp(-sTransform*sTransform);
   }
 
   double GetInverseProbability(double proba) // returns sTransform
   {
     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);
   }
 
   double PlotInverse(double* x, double* )
   {
     return GetInverseProbability(x[0]);
   }
 
   double Plot(double* x, double* )
   {
     return GetDifferential(x[0]);
   }
 
 
   void InitialiseInverseProbability(double xmax = 3e28)
   {
     // 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);
   }
 
   std::vector<double> fInverse;
   int fNbins;
   double fEpsilon;
 };
 
 #endif /* SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_MODELS_G4DNASMOLUCHOWSKIDIFFUSION_HH_ */

This page was automatically generated by the 2.3.7 LXR engine. The LXR team