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