| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 * G4DNASmoluchowskiDiffusion.cc 27 * G4DNASmoluchowskiDiffusion.cc
28 * 28 *
29 * Created on: 2 févr. 2015 29 * Created on: 2 févr. 2015
30 * Author: matkara 30 * Author: matkara
31 */ 31 */
32 32
33 //#define DNADEV_TEST 33 //#define DNADEV_TEST
34 34
35 #ifdef DNADEV_TEST 35 #ifdef DNADEV_TEST
36 #include "" 36 #include ""
37 #else 37 #else
38 #include "G4DNASmoluchowskiDiffusion.hh" 38 #include "G4DNASmoluchowskiDiffusion.hh"
39 #endif 39 #endif
40 40
41 //#if __cplusplus >= 201103L 41 //#if __cplusplus >= 201103L
42 #ifdef DNADEV_TEST 42 #ifdef DNADEV_TEST
43 #include "TRint.h" 43 #include "TRint.h"
44 #include "TCanvas.h" 44 #include "TCanvas.h"
45 #include "TH1D.h" 45 #include "TH1D.h"
46 #include "TRandom.h" 46 #include "TRandom.h"
47 #include "TMath.h" 47 #include "TMath.h"
48 #endif 48 #endif
49 49
50 G4DNASmoluchowskiDiffusion::G4DNASmoluchowskiD 50 G4DNASmoluchowskiDiffusion::G4DNASmoluchowskiDiffusion(double epsilon) : fEpsilon(epsilon)
51 { 51 {
52 fNbins = (int) trunc(1/fEpsilon); 52 fNbins = (int) trunc(1/fEpsilon);
53 // std::cout << "fNbins: " << fNbins << std: 53 // std::cout << "fNbins: " << fNbins << std::endl;
54 #ifdef DNADEV 54 #ifdef DNADEV
55 assert(fNbins > 0); 55 assert(fNbins > 0);
56 #endif 56 #endif
57 fInverse.resize(fNbins+2); // trunc sous-est 57 fInverse.resize(fNbins+2); // trunc sous-estime + borne max a rajouter ==> 2
58 58
59 // std::cout << "fInverse.capacity(): "<< fI 59 // std::cout << "fInverse.capacity(): "<< fInverse.capacity() << std::endl;
60 } 60 }
61 61
62 G4DNASmoluchowskiDiffusion::~G4DNASmoluchowski 62 G4DNASmoluchowskiDiffusion::~G4DNASmoluchowskiDiffusion()
63 = default; << 63 {
>> 64 }
64 //#endif 65 //#endif
65 66
66 // --> G4DNASmoluchowskiDiffusion -- DEVELOPME 67 // --> G4DNASmoluchowskiDiffusion -- DEVELOPMENT TEST
67 #ifdef DNADEV_TEST 68 #ifdef DNADEV_TEST
68 69
69 static G4DNASmoluchowskiDiffusion gDiff; 70 static G4DNASmoluchowskiDiffusion gDiff;
70 71
71 double time_test = 1e-6 /*s*/; 72 double time_test = 1e-6 /*s*/;
72 double D = 4.9e-9 /*m2/s*/; 73 double D = 4.9e-9 /*m2/s*/;
73 double test_distance = 1e-9; // m 74 double test_distance = 1e-9; // m
74 75
75 double Plot(double* x, double* ) 76 double Plot(double* x, double* )
76 { 77 {
77 double diff = gDiff.GetDensityProbability(x[ 78 double diff = gDiff.GetDensityProbability(x[0], time_test, D);
78 return diff; 79 return diff;
79 } 80 }
80 81
81 static double InvErfc(double x) 82 static double InvErfc(double x)
82 { 83 {
83 return TMath::ErfcInverse(x); 84 return TMath::ErfcInverse(x);
84 } 85 }
85 86
86 Axis_t* BinLogX(Int_t bins, Axis_t from, Axis_ 87 Axis_t* BinLogX(Int_t bins, Axis_t from, Axis_t to) // en puissance de 10
87 { 88 {
88 Axis_t width = (to - from) / bins; 89 Axis_t width = (to - from) / bins;
89 Axis_t *new_bins = new Axis_t[bins + 1]; 90 Axis_t *new_bins = new Axis_t[bins + 1];
90 91
91 for (int i = 0; i <= bins; i++) { 92 for (int i = 0; i <= bins; i++) {
92 new_bins[i] = TMath::Power(10, from + i * 93 new_bins[i] = TMath::Power(10, from + i * width);
93 // std::cout << new_bins[i] << std::endl; 94 // std::cout << new_bins[i] << std::endl;
94 } 95 }
95 return new_bins; 96 return new_bins;
96 } 97 }
97 98
98 int main(int argc, char **argv) 99 int main(int argc, char **argv)
99 { 100 {
100 gDiff.InitialiseInverseProbability(); 101 gDiff.InitialiseInverseProbability();
101 // srand (time(NULL)); 102 // srand (time(NULL));
102 TRint* root = new TRint("G4DNASmoluchowskiDi 103 TRint* root = new TRint("G4DNASmoluchowskiDiffusion",&argc, argv);
103 double interval = 1e-5; 104 double interval = 1e-5;
104 G4DNASmoluchowskiDiffusion* diff = new G4DNA 105 G4DNASmoluchowskiDiffusion* diff = new G4DNASmoluchowskiDiffusion(interval);
105 diff->InitialiseInverseProbability(); 106 diff->InitialiseInverseProbability();
106 107
107 // for(size_t i = 0 ; i < diff->fInverse.size 108 // for(size_t i = 0 ; i < diff->fInverse.size() ; ++i)
108 // { 109 // {
109 // std::cout << i*interval << " "<< diff->f 110 // std::cout << i*interval << " "<< diff->fInverse[i] << std::endl;
110 // } 111 // }
111 112
112 std::cout << diff->fInverse.size() << std::e 113 std::cout << diff->fInverse.size() << std::endl;
113 114
114 TCanvas* canvas = new TCanvas(); 115 TCanvas* canvas = new TCanvas();
115 //canvas->SetLogx(); 116 //canvas->SetLogx();
116 //canvas->SetLogy(); 117 //canvas->SetLogy();
117 // 118 //
118 // TF1 * f = new TF1("f",diff,&G4DNASmoluchow 119 // TF1 * f = new TF1("f",diff,&G4DNASmoluchowskiDiffusion::PlotInverse,0,10,0,"G4DNASmoluchowskiDiffusion","Plot"); // create TF1 class.
119 // f->SetNpx(100000); 120 // f->SetNpx(100000);
120 // f->Draw(); 121 // f->Draw();
121 // canvas->Draw(); 122 // canvas->Draw();
122 // 123 //
123 // canvas = new TCanvas(); 124 // canvas = new TCanvas();
124 TH1D* h1 = new TH1D("h1", "h1", 100, 0., 1e- 125 TH1D* h1 = new TH1D("h1", "h1", 100, 0., 1e-6);
125 double distance = -1; 126 double distance = -1;
126 127
127 int N = 100000; 128 int N = 100000;
128 129
129 for(size_t i = 0 ; i < N ; ++i) 130 for(size_t i = 0 ; i < N ; ++i)
130 { 131 {
131 distance = diff->GetRandomDistance(time_te 132 distance = diff->GetRandomDistance(time_test,D);
132 h1->Fill(distance); 133 h1->Fill(distance);
133 //std::cout << distance << std::endl; 134 //std::cout << distance << std::endl;
134 } 135 }
135 136
136 double scalf; 137 double scalf;
137 138
138 { 139 {
139 int integral_h1 = h1->Integral(); 140 int integral_h1 = h1->Integral();
140 h1->Scale(1./integral_h1); 141 h1->Scale(1./integral_h1);
141 scalf=h1->GetBinWidth ( 1 ) ; 142 scalf=h1->GetBinWidth ( 1 ) ;
142 h1->Scale(1./scalf); 143 h1->Scale(1./scalf);
143 h1->GetXaxis()->SetTitle("distance"); 144 h1->GetXaxis()->SetTitle("distance");
144 } 145 }
145 146
146 TH1D* h2 = new TH1D("h2", "h2", 100, 0., 1e- 147 TH1D* h2 = new TH1D("h2", "h2", 100, 0., 1e-6);
147 TH1D* h_irt_distance = new TH1D("h2", "h2", 148 TH1D* h_irt_distance = new TH1D("h2", "h2", 100, 0., 1e-6);
148 149
149 for(size_t i = 0 ; i < N ; ++i) 150 for(size_t i = 0 ; i < N ; ++i)
150 { 151 {
151 double x = std::sqrt(2*D*time_test)*root_r 152 double x = std::sqrt(2*D*time_test)*root_random.Gaus();
152 double y = std::sqrt(2*D*time_test)*root_r 153 double y = std::sqrt(2*D*time_test)*root_random.Gaus();
153 double z = std::sqrt(2*D*time_test)*root_r 154 double z = std::sqrt(2*D*time_test)*root_random.Gaus();
154 155
155 distance = std::sqrt(x*x+y*y+z*z); 156 distance = std::sqrt(x*x+y*y+z*z);
156 h2->Fill(distance); 157 h2->Fill(distance);
157 //std::cout << distance << std::endl; 158 //std::cout << distance << std::endl;
158 159
159 double proba = root_random.Rndm(); 160 double proba = root_random.Rndm();
160 double irt_distance = InvErfc(proba)*2*std 161 double irt_distance = InvErfc(proba)*2*std::sqrt(D*time_test);
161 h_irt_distance->Fill(irt_distance); 162 h_irt_distance->Fill(irt_distance);
162 } 163 }
163 164
164 { 165 {
165 int integral_h2 = h2->Integral(); 166 int integral_h2 = h2->Integral();
166 h2->Scale(1./integral_h2); 167 h2->Scale(1./integral_h2);
167 scalf=h2->GetBinWidth ( 1 ) ; 168 scalf=h2->GetBinWidth ( 1 ) ;
168 h2->Scale(1./scalf); 169 h2->Scale(1./scalf);
169 } 170 }
170 171
171 { 172 {
172 int integral_h_irt_distance = h_irt_distance 173 int integral_h_irt_distance = h_irt_distance->Integral();
173 h_irt_distance->Scale(1./integral_h_irt_dist 174 h_irt_distance->Scale(1./integral_h_irt_distance);
174 scalf = h_irt_distance->GetBinWidth ( 1 ) ; 175 scalf = h_irt_distance->GetBinWidth ( 1 ) ;
175 h_irt_distance->Scale(1./scalf); 176 h_irt_distance->Scale(1./scalf);
176 h_irt_distance->GetXaxis()->SetTitle("distan 177 h_irt_distance->GetXaxis()->SetTitle("distance");
177 } 178 }
178 179
179 180
180 TF1 * f2 = new TF1("f2",&Plot,0,1e-6,0,"Plot 181 TF1 * f2 = new TF1("f2",&Plot,0,1e-6,0,"Plot"); // create TF1 class.
181 //f2->SetNpx(1000); 182 //f2->SetNpx(1000);
182 h1->Draw(); 183 h1->Draw();
183 // h1->DrawNormalized(); 184 // h1->DrawNormalized();
184 f2->Draw("SAME"); 185 f2->Draw("SAME");
185 h2->Draw("SAME"); 186 h2->Draw("SAME");
186 h_irt_distance->Draw("SAME"); 187 h_irt_distance->Draw("SAME");
187 double integral = f2->Integral(0., 1e-6); 188 double integral = f2->Integral(0., 1e-6);
188 std::cout << "integral = " << integral << st 189 std::cout << "integral = " << integral << std::endl;
189 std::cout << "integral h1 = " << h1->Integra 190 std::cout << "integral h1 = " << h1->Integral() << std::endl;
190 canvas->Draw(); 191 canvas->Draw();
191 192
192 std::vector<double> rdm(3); 193 std::vector<double> rdm(3);
193 int nbins = 100; 194 int nbins = 100;
194 Axis_t* bins = BinLogX(nbins, -12, -1); 195 Axis_t* bins = BinLogX(nbins, -12, -1);
195 196
196 TH1D* h3 = new TH1D("h3", "h3", 100, bins); 197 TH1D* h3 = new TH1D("h3", "h3", 100, bins);
197 TH1D* h4 = new TH1D("h4", "h4", 100, bins); 198 TH1D* h4 = new TH1D("h4", "h4", 100, bins);
198 TH1D* h_irt = new TH1D("h_irt", "h_irt", 100 199 TH1D* h_irt = new TH1D("h_irt", "h_irt", 100, bins);
199 200
200 for(size_t i = 0 ; i < N ; ++i) 201 for(size_t i = 0 ; i < N ; ++i)
201 { 202 {
202 for(size_t j = 0 ; j < 3 ; ++j) 203 for(size_t j = 0 ; j < 3 ; ++j)
203 rdm[j] = root_random.Gaus(); 204 rdm[j] = root_random.Gaus();
204 205
205 double denum = 1./(rdm[0]*rdm[0] + rdm[1]* 206 double denum = 1./(rdm[0]*rdm[0] + rdm[1]*rdm[1] + rdm[2]*rdm[2]);
206 207
207 double t = ((test_distance*test_distance)* 208 double t = ((test_distance*test_distance)*denum)*1./(2*D);
208 h3->Fill(t); 209 h3->Fill(t);
209 210
210 double t_h4 = diff->GetRandomTime(test_di 211 double t_h4 = diff->GetRandomTime(test_distance,D);
211 h4->Fill(t_h4); 212 h4->Fill(t_h4);
212 // std::cout << t << " " << t_h4 << std::e 213 // std::cout << t << " " << t_h4 << std::endl;
213 214
214 double proba = root_random.Rndm(); 215 double proba = root_random.Rndm();
215 double t_irt = 1./(4*D)*std::pow((test_di 216 double t_irt = 1./(4*D)*std::pow((test_distance)/InvErfc(proba),2);
216 h_irt ->Fill(t_irt); 217 h_irt ->Fill(t_irt);
217 } 218 }
218 219
219 { 220 {
220 TCanvas* c1 = new TCanvas(); 221 TCanvas* c1 = new TCanvas();
221 c1->SetLogx(); 222 c1->SetLogx();
222 int integral_h3 = h3->Integral(); 223 int integral_h3 = h3->Integral();
223 h3->Scale(1./integral_h3); 224 h3->Scale(1./integral_h3);
224 scalf=h3->GetBinWidth ( 1 ) ; 225 scalf=h3->GetBinWidth ( 1 ) ;
225 h3->Scale(1./scalf); 226 h3->Scale(1./scalf);
226 h3->SetLineColor(1); 227 h3->SetLineColor(1);
227 h3->GetXaxis()->SetTitle("time");; 228 h3->GetXaxis()->SetTitle("time");;
228 h3->Draw(); 229 h3->Draw();
229 } 230 }
230 231
231 { 232 {
232 // TCanvas* c1 = new TCanvas(); 233 // TCanvas* c1 = new TCanvas();
233 // c1->SetLogx(); 234 // c1->SetLogx();
234 int integral_h4 = h4->Integral(); 235 int integral_h4 = h4->Integral();
235 h4->Scale(1./integral_h4); 236 h4->Scale(1./integral_h4);
236 scalf=h4->GetBinWidth ( 1 ) ; 237 scalf=h4->GetBinWidth ( 1 ) ;
237 h4->Scale(1./scalf); 238 h4->Scale(1./scalf);
238 h4->SetLineColor(6); 239 h4->SetLineColor(6);
239 h4->Draw("SAME"); 240 h4->Draw("SAME");
240 // h4->Draw("SAME"); 241 // h4->Draw("SAME");
241 } 242 }
242 243
243 { 244 {
244 // TCanvas* c1 = new TCanvas(); 245 // TCanvas* c1 = new TCanvas();
245 // c1->SetLogx(); 246 // c1->SetLogx();
246 int integral_h_irt = h_irt->Integral(); 247 int integral_h_irt = h_irt->Integral();
247 h_irt->Scale(1./integral_h_irt); 248 h_irt->Scale(1./integral_h_irt);
248 scalf=h_irt->GetBinWidth ( 1 ) ; 249 scalf=h_irt->GetBinWidth ( 1 ) ;
249 h_irt->Scale(1./scalf); 250 h_irt->Scale(1./scalf);
250 h_irt->SetLineColor(4); 251 h_irt->SetLineColor(4);
251 h_irt->Draw("SAME"); 252 h_irt->Draw("SAME");
252 // h4->Draw("SAME"); 253 // h4->Draw("SAME");
253 } 254 }
254 root->Run(); 255 root->Run();
255 return 0; 256 return 0;
256 } 257 }
257 #endif 258 #endif
258 259