Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASmoluchowskiDiffusion.cc
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.cc
28  *
29  * Created on: 2 févr. 2015
30  * Author: matkara
31  */
32 
33 //#define DNADEV_TEST
34 
35 #ifdef DNADEV_TEST
36 #include "../include/G4DNASmoluchowskiDiffusion.hh"
37 #else
39 #endif
40 
41 //#if __cplusplus >= 201103L
42 #ifdef DNADEV_TEST
43 #include "TRint.h"
44 #include "TCanvas.h"
45 #include "TH1D.h"
46 #include "TRandom.h"
47 #include "TMath.h"
48 #endif
49 
51 {
52  fNbins = (int) trunc(1/fEpsilon);
53  // std::cout << "fNbins: " << fNbins << std::endl;
54 #ifdef DNADEV
55  assert(fNbins > 0);
56 #endif
57  fInverse.resize(fNbins+2); // trunc sous-estime + borne max a rajouter ==> 2
58 
59  // std::cout << "fInverse.capacity(): "<< fInverse.capacity() << std::endl;
60 }
61 
63 {
64 }
65 //#endif
66 
67 // --> G4DNASmoluchowskiDiffusion -- DEVELOPMENT TEST
68 #ifdef DNADEV_TEST
69 
70 static G4DNASmoluchowskiDiffusion gDiff;
71 
72 double time_test = 1e-6 /*s*/;
73 double D = 4.9e-9 /*m2/s*/;
74 double test_distance = 1e-9; // m
75 
76 double Plot(double* x, double* )
77 {
78  double diff = gDiff.GetDensityProbability(x[0], time_test, D);
79  return diff;
80 }
81 
82 static double InvErfc(double x)
83 {
84  return TMath::ErfcInverse(x);
85 }
86 
87 Axis_t* BinLogX(Int_t bins, Axis_t from, Axis_t to) // en puissance de 10
88 {
89  Axis_t width = (to - from) / bins;
90  Axis_t *new_bins = new Axis_t[bins + 1];
91 
92  for (int i = 0; i <= bins; i++) {
93  new_bins[i] = TMath::Power(10, from + i * width);
94 // std::cout << new_bins[i] << std::endl;
95  }
96  return new_bins;
97 }
98 
99 int main(int argc, char **argv)
100 {
102 // srand (time(NULL));
103  TRint* root = new TRint("G4DNASmoluchowskiDiffusion",&argc, argv);
104  double interval = 1e-5;
107 
108 // for(size_t i = 0 ; i < diff->fInverse.size() ; ++i)
109 // {
110 // std::cout << i*interval << " "<< diff->fInverse[i] << std::endl;
111 // }
112 
113  std::cout << diff->fInverse.size() << std::endl;
114 
115  TCanvas* canvas = new TCanvas();
116  //canvas->SetLogx();
117  //canvas->SetLogy();
118 //
119 // TF1 * f = new TF1("f",diff,&G4DNASmoluchowskiDiffusion::PlotInverse,0,10,0,"G4DNASmoluchowskiDiffusion","Plot"); // create TF1 class.
120 // f->SetNpx(100000);
121 // f->Draw();
122 // canvas->Draw();
123 //
124 // canvas = new TCanvas();
125  TH1D* h1 = new TH1D("h1", "h1", 100, 0., 1e-6);
126  double distance = -1;
127 
128  int N = 100000;
129 
130  for(size_t i = 0 ; i < N ; ++i)
131  {
132  distance = diff->GetRandomDistance(time_test,D);
133  h1->Fill(distance);
134  //std::cout << distance << std::endl;
135  }
136 
137  double scalf;
138 
139  {
140  int integral_h1 = h1->Integral();
141  h1->Scale(1./integral_h1);
142  scalf=h1->GetBinWidth ( 1 ) ;
143  h1->Scale(1./scalf);
144  h1->GetXaxis()->SetTitle("distance");
145  }
146 
147  TH1D* h2 = new TH1D("h2", "h2", 100, 0., 1e-6);
148  TH1D* h_irt_distance = new TH1D("h2", "h2", 100, 0., 1e-6);
149 
150  for(size_t i = 0 ; i < N ; ++i)
151  {
152  double x = std::sqrt(2*D*time_test)*root_random.Gaus();
153  double y = std::sqrt(2*D*time_test)*root_random.Gaus();
154  double z = std::sqrt(2*D*time_test)*root_random.Gaus();
155 
156  distance = std::sqrt(x*x+y*y+z*z);
157  h2->Fill(distance);
158  //std::cout << distance << std::endl;
159 
160  double proba = root_random.Rndm();
161  double irt_distance = InvErfc(proba)*2*std::sqrt(D*time_test);
162  h_irt_distance->Fill(irt_distance);
163  }
164 
165  {
166  int integral_h2 = h2->Integral();
167  h2->Scale(1./integral_h2);
168  scalf=h2->GetBinWidth ( 1 ) ;
169  h2->Scale(1./scalf);
170  }
171 
172  {
173  int integral_h_irt_distance = h_irt_distance->Integral();
174  h_irt_distance->Scale(1./integral_h_irt_distance);
175  scalf = h_irt_distance->GetBinWidth ( 1 ) ;
176  h_irt_distance->Scale(1./scalf);
177  h_irt_distance->GetXaxis()->SetTitle("distance");
178  }
179 
180 
181  TF1 * f2 = new TF1("f2",&Plot,0,1e-6,0,"Plot"); // create TF1 class.
182  //f2->SetNpx(1000);
183  h1->Draw();
184  // h1->DrawNormalized();
185  f2->Draw("SAME");
186  h2->Draw("SAME");
187  h_irt_distance->Draw("SAME");
188  double integral = f2->Integral(0., 1e-6);
189  std::cout << "integral = " << integral << std::endl;
190  std::cout << "integral h1 = " << h1->Integral() << std::endl;
191  canvas->Draw();
192 
193  std::vector<double> rdm(3);
194  int nbins = 100;
195  Axis_t* bins = BinLogX(nbins, -12, -1);
196 
197  TH1D* h3 = new TH1D("h3", "h3", 100, bins);
198  TH1D* h4 = new TH1D("h4", "h4", 100, bins);
199  TH1D* h_irt = new TH1D("h_irt", "h_irt", 100, bins);
200 
201  for(size_t i = 0 ; i < N ; ++i)
202  {
203  for(size_t j = 0 ; j < 3 ; ++j)
204  rdm[j] = root_random.Gaus();
205 
206  double denum = 1./(rdm[0]*rdm[0] + rdm[1]*rdm[1] + rdm[2]*rdm[2]);
207 
208  double t = ((test_distance*test_distance)*denum)*1./(2*D);
209  h3->Fill(t);
210 
211  double t_h4 = diff->GetRandomTime(test_distance,D);
212  h4->Fill(t_h4);
213 // std::cout << t << " " << t_h4 << std::endl;
214 
215  double proba = root_random.Rndm();
216  double t_irt = 1./(4*D)*std::pow((test_distance)/InvErfc(proba),2);
217  h_irt ->Fill(t_irt);
218  }
219 
220  {
221  TCanvas* c1 = new TCanvas();
222  c1->SetLogx();
223  int integral_h3 = h3->Integral();
224  h3->Scale(1./integral_h3);
225  scalf=h3->GetBinWidth ( 1 ) ;
226  h3->Scale(1./scalf);
227  h3->SetLineColor(1);
228  h3->GetXaxis()->SetTitle("time");;
229  h3->Draw();
230  }
231 
232  {
233 // TCanvas* c1 = new TCanvas();
234 // c1->SetLogx();
235  int integral_h4 = h4->Integral();
236  h4->Scale(1./integral_h4);
237  scalf=h4->GetBinWidth ( 1 ) ;
238  h4->Scale(1./scalf);
239  h4->SetLineColor(6);
240  h4->Draw("SAME");
241  // h4->Draw("SAME");
242  }
243 
244  {
245 // TCanvas* c1 = new TCanvas();
246 // c1->SetLogx();
247  int integral_h_irt = h_irt->Integral();
248  h_irt->Scale(1./integral_h_irt);
249  scalf=h_irt->GetBinWidth ( 1 ) ;
250  h_irt->Scale(1./scalf);
251  h_irt->SetLineColor(4);
252  h_irt->Draw("SAME");
253  // h4->Draw("SAME");
254  }
255  root->Run();
256  return 0;
257 }
258 #endif
void InitialiseInverseProbability(double xmax=3e28)
#define width
int main(int argc, char **argv)
Definition: genwindef.cpp:30
tuple x
Definition: test.py:50
static double InvErfc(double x)
G4DNASmoluchowskiDiffusion(double epsilon=1e-5)
G4int Int_t
double GetRandomTime(double distance, double D)
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
double D(double temp)
tuple z
Definition: test.py:28
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
static double GetDensityProbability(double r, double _time, double D)
double epsilon(double density, double temperature)
tuple c1
Definition: plottest35.py:14
double GetRandomDistance(double _time, double D)