| 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.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 ""
37 #else
38 #include "G4DNASmoluchowskiDiffusion.hh"
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
50 G4DNASmoluchowskiDiffusion::G4DNASmoluchowskiDiffusion(double epsilon) : fEpsilon(epsilon)
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
62 G4DNASmoluchowskiDiffusion::~G4DNASmoluchowskiDiffusion()
63 = default;
64 //#endif
65
66 // --> G4DNASmoluchowskiDiffusion -- DEVELOPMENT TEST
67 #ifdef DNADEV_TEST
68
69 static G4DNASmoluchowskiDiffusion gDiff;
70
71 double time_test = 1e-6 /*s*/;
72 double D = 4.9e-9 /*m2/s*/;
73 double test_distance = 1e-9; // m
74
75 double Plot(double* x, double* )
76 {
77 double diff = gDiff.GetDensityProbability(x[0], time_test, D);
78 return diff;
79 }
80
81 static double InvErfc(double x)
82 {
83 return TMath::ErfcInverse(x);
84 }
85
86 Axis_t* BinLogX(Int_t bins, Axis_t from, Axis_t to) // en puissance de 10
87 {
88 Axis_t width = (to - from) / bins;
89 Axis_t *new_bins = new Axis_t[bins + 1];
90
91 for (int i = 0; i <= bins; i++) {
92 new_bins[i] = TMath::Power(10, from + i * width);
93 // std::cout << new_bins[i] << std::endl;
94 }
95 return new_bins;
96 }
97
98 int main(int argc, char **argv)
99 {
100 gDiff.InitialiseInverseProbability();
101 // srand (time(NULL));
102 TRint* root = new TRint("G4DNASmoluchowskiDiffusion",&argc, argv);
103 double interval = 1e-5;
104 G4DNASmoluchowskiDiffusion* diff = new G4DNASmoluchowskiDiffusion(interval);
105 diff->InitialiseInverseProbability();
106
107 // for(size_t i = 0 ; i < diff->fInverse.size() ; ++i)
108 // {
109 // std::cout << i*interval << " "<< diff->fInverse[i] << std::endl;
110 // }
111
112 std::cout << diff->fInverse.size() << std::endl;
113
114 TCanvas* canvas = new TCanvas();
115 //canvas->SetLogx();
116 //canvas->SetLogy();
117 //
118 // TF1 * f = new TF1("f",diff,&G4DNASmoluchowskiDiffusion::PlotInverse,0,10,0,"G4DNASmoluchowskiDiffusion","Plot"); // create TF1 class.
119 // f->SetNpx(100000);
120 // f->Draw();
121 // canvas->Draw();
122 //
123 // canvas = new TCanvas();
124 TH1D* h1 = new TH1D("h1", "h1", 100, 0., 1e-6);
125 double distance = -1;
126
127 int N = 100000;
128
129 for(size_t i = 0 ; i < N ; ++i)
130 {
131 distance = diff->GetRandomDistance(time_test,D);
132 h1->Fill(distance);
133 //std::cout << distance << std::endl;
134 }
135
136 double scalf;
137
138 {
139 int integral_h1 = h1->Integral();
140 h1->Scale(1./integral_h1);
141 scalf=h1->GetBinWidth ( 1 ) ;
142 h1->Scale(1./scalf);
143 h1->GetXaxis()->SetTitle("distance");
144 }
145
146 TH1D* h2 = new TH1D("h2", "h2", 100, 0., 1e-6);
147 TH1D* h_irt_distance = new TH1D("h2", "h2", 100, 0., 1e-6);
148
149 for(size_t i = 0 ; i < N ; ++i)
150 {
151 double x = std::sqrt(2*D*time_test)*root_random.Gaus();
152 double y = std::sqrt(2*D*time_test)*root_random.Gaus();
153 double z = std::sqrt(2*D*time_test)*root_random.Gaus();
154
155 distance = std::sqrt(x*x+y*y+z*z);
156 h2->Fill(distance);
157 //std::cout << distance << std::endl;
158
159 double proba = root_random.Rndm();
160 double irt_distance = InvErfc(proba)*2*std::sqrt(D*time_test);
161 h_irt_distance->Fill(irt_distance);
162 }
163
164 {
165 int integral_h2 = h2->Integral();
166 h2->Scale(1./integral_h2);
167 scalf=h2->GetBinWidth ( 1 ) ;
168 h2->Scale(1./scalf);
169 }
170
171 {
172 int integral_h_irt_distance = h_irt_distance->Integral();
173 h_irt_distance->Scale(1./integral_h_irt_distance);
174 scalf = h_irt_distance->GetBinWidth ( 1 ) ;
175 h_irt_distance->Scale(1./scalf);
176 h_irt_distance->GetXaxis()->SetTitle("distance");
177 }
178
179
180 TF1 * f2 = new TF1("f2",&Plot,0,1e-6,0,"Plot"); // create TF1 class.
181 //f2->SetNpx(1000);
182 h1->Draw();
183 // h1->DrawNormalized();
184 f2->Draw("SAME");
185 h2->Draw("SAME");
186 h_irt_distance->Draw("SAME");
187 double integral = f2->Integral(0., 1e-6);
188 std::cout << "integral = " << integral << std::endl;
189 std::cout << "integral h1 = " << h1->Integral() << std::endl;
190 canvas->Draw();
191
192 std::vector<double> rdm(3);
193 int nbins = 100;
194 Axis_t* bins = BinLogX(nbins, -12, -1);
195
196 TH1D* h3 = new TH1D("h3", "h3", 100, bins);
197 TH1D* h4 = new TH1D("h4", "h4", 100, bins);
198 TH1D* h_irt = new TH1D("h_irt", "h_irt", 100, bins);
199
200 for(size_t i = 0 ; i < N ; ++i)
201 {
202 for(size_t j = 0 ; j < 3 ; ++j)
203 rdm[j] = root_random.Gaus();
204
205 double denum = 1./(rdm[0]*rdm[0] + rdm[1]*rdm[1] + rdm[2]*rdm[2]);
206
207 double t = ((test_distance*test_distance)*denum)*1./(2*D);
208 h3->Fill(t);
209
210 double t_h4 = diff->GetRandomTime(test_distance,D);
211 h4->Fill(t_h4);
212 // std::cout << t << " " << t_h4 << std::endl;
213
214 double proba = root_random.Rndm();
215 double t_irt = 1./(4*D)*std::pow((test_distance)/InvErfc(proba),2);
216 h_irt ->Fill(t_irt);
217 }
218
219 {
220 TCanvas* c1 = new TCanvas();
221 c1->SetLogx();
222 int integral_h3 = h3->Integral();
223 h3->Scale(1./integral_h3);
224 scalf=h3->GetBinWidth ( 1 ) ;
225 h3->Scale(1./scalf);
226 h3->SetLineColor(1);
227 h3->GetXaxis()->SetTitle("time");;
228 h3->Draw();
229 }
230
231 {
232 // TCanvas* c1 = new TCanvas();
233 // c1->SetLogx();
234 int integral_h4 = h4->Integral();
235 h4->Scale(1./integral_h4);
236 scalf=h4->GetBinWidth ( 1 ) ;
237 h4->Scale(1./scalf);
238 h4->SetLineColor(6);
239 h4->Draw("SAME");
240 // h4->Draw("SAME");
241 }
242
243 {
244 // TCanvas* c1 = new TCanvas();
245 // c1->SetLogx();
246 int integral_h_irt = h_irt->Integral();
247 h_irt->Scale(1./integral_h_irt);
248 scalf=h_irt->GetBinWidth ( 1 ) ;
249 h_irt->Scale(1./scalf);
250 h_irt->SetLineColor(4);
251 h_irt->Draw("SAME");
252 // h4->Draw("SAME");
253 }
254 root->Run();
255 return 0;
256 }
257 #endif
258