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