| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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::G4DNASmoluchowskiD
51 {
52 fNbins = (int) trunc(1/fEpsilon);
53 // std::cout << "fNbins: " << fNbins << std:
54 #ifdef DNADEV
55 assert(fNbins > 0);
56 #endif
57 fInverse.resize(fNbins+2); // trunc sous-est
58
59 // std::cout << "fInverse.capacity(): "<< fI
60 }
61
62 G4DNASmoluchowskiDiffusion::~G4DNASmoluchowski
63 = default;
64 //#endif
65
66 // --> G4DNASmoluchowskiDiffusion -- DEVELOPME
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[
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_
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 *
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("G4DNASmoluchowskiDi
103 double interval = 1e-5;
104 G4DNASmoluchowskiDiffusion* diff = new G4DNA
105 diff->InitialiseInverseProbability();
106
107 // for(size_t i = 0 ; i < diff->fInverse.size
108 // {
109 // std::cout << i*interval << " "<< diff->f
110 // }
111
112 std::cout << diff->fInverse.size() << std::e
113
114 TCanvas* canvas = new TCanvas();
115 //canvas->SetLogx();
116 //canvas->SetLogy();
117 //
118 // TF1 * f = new TF1("f",diff,&G4DNASmoluchow
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-
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_te
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-
147 TH1D* h_irt_distance = new TH1D("h2", "h2",
148
149 for(size_t i = 0 ; i < N ; ++i)
150 {
151 double x = std::sqrt(2*D*time_test)*root_r
152 double y = std::sqrt(2*D*time_test)*root_r
153 double z = std::sqrt(2*D*time_test)*root_r
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
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
173 h_irt_distance->Scale(1./integral_h_irt_dist
174 scalf = h_irt_distance->GetBinWidth ( 1 ) ;
175 h_irt_distance->Scale(1./scalf);
176 h_irt_distance->GetXaxis()->SetTitle("distan
177 }
178
179
180 TF1 * f2 = new TF1("f2",&Plot,0,1e-6,0,"Plot
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 << st
189 std::cout << "integral h1 = " << h1->Integra
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
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]*
206
207 double t = ((test_distance*test_distance)*
208 h3->Fill(t);
209
210 double t_h4 = diff->GetRandomTime(test_di
211 h4->Fill(t_h4);
212 // std::cout << t << " " << t_h4 << std::e
213
214 double proba = root_random.Rndm();
215 double t_irt = 1./(4*D)*std::pow((test_di
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