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 // $Id: G4GoudsmitSaundersonTable.cc 107234 2017-11-06 11:53:54Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ----------------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class implementation file 30 // GEANT4 Class implementation file 30 // 31 // 31 // File name: G4GoudsmitSaundersonTable 32 // File name: G4GoudsmitSaundersonTable 32 // 33 // 33 // Author: Mihaly Novak / (Omrane Kadri 34 // Author: Mihaly Novak / (Omrane Kadri) 34 // 35 // 35 // Creation date: 20.02.2009 36 // Creation date: 20.02.2009 36 // 37 // 37 // Class description: 38 // Class description: 38 // Class to handle multiple scattering angul 39 // Class to handle multiple scattering angular distributions precomputed by 39 // using Kawrakow-Bielajew Goudsmit-Saunders 40 // using Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened 40 // Rutherford DCS for elastic scattering of 41 // Rutherford DCS for elastic scattering of electrons/positrons [1,2]. This 41 // class is used by G4GoudsmitSaundersonMscM 42 // class is used by G4GoudsmitSaundersonMscModel to sample the angular 42 // deflection of electrons/positrons after t 43 // deflection of electrons/positrons after travelling a given path. 43 // 44 // 44 // Modifications: 45 // Modifications: 45 // 04.03.2009 V.Ivanchenko cleanup and format 46 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style 46 // 26.08.2009 O.Kadri: avoiding unuseful calcu 47 // 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root 47 // finding parameter error 48 // finding parameter error's within SampleTheta method 48 // 08.02.2010 O.Kadri: reduce delared variable 49 // 08.02.2010 O.Kadri: reduce delared variables; reduce error of finding root 49 // in secant method 50 // in secant method 50 // 26.03.2010 O.Kadri: minimum of used arrays 51 // 26.03.2010 O.Kadri: minimum of used arrays in computation within the dichotomie 51 // finding method the erro 52 // finding method the error was the lowest value of uvalues 52 // 12.05.2010 O.Kadri: changing of sqrt((b-a)* 53 // 12.05.2010 O.Kadri: changing of sqrt((b-a)*(b-a)) with fabs(b-a) 53 // 18.05.2015 M. Novak This class has been com 54 // 18.05.2015 M. Novak This class has been completely replaced (only the original 54 // class name was kept; class descr 55 // class name was kept; class description was also inserted): 55 // A new version of Kawrakow-Bielaj 56 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model 56 // based on the screened Rutherford 57 // based on the screened Rutherford DCS for elastic scattering of 57 // electrons/positrons has been int 58 // electrons/positrons has been introduced[1,2]. The corresponding MSC 58 // angular distributions over a 2D 59 // angular distributions over a 2D parameter grid have been recomputed 59 // and the CDFs are now stored in a 60 // and the CDFs are now stored in a variable transformed (smooth) form 60 // together with the corresponding 61 // together with the corresponding rational interpolation parameters. 61 // The new version is several times 62 // The new version is several times faster, more robust and accurate 62 // compared to the earlier version 63 // compared to the earlier version (G4GoudsmitSaundersonMscModel class 63 // that use these data has been als 64 // that use these data has been also completely replaced) 64 // 28.04.2017 M. Novak: New representation of 65 // 28.04.2017 M. Novak: New representation of the angular distribution data with 65 // significantly reduced data size. 66 // significantly reduced data size. 66 // 23.08.2017 M. Novak: Added funtionality to 67 // 23.08.2017 M. Novak: Added funtionality to handle Mott-correction to the 67 // base GS angular distributions an 68 // base GS angular distributions and some other factors (screening 68 // parameter, first and second mome 69 // parameter, first and second moments) when Mott-correction is 69 // activated in the GS-MSC model. 70 // activated in the GS-MSC model. 70 // 71 // 71 // References: 72 // References: 72 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-20 73 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-208 73 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 74 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336 74 // 75 // 75 // ------------------------------------------- 76 // ----------------------------------------------------------------------------- 76 77 77 #include "G4GoudsmitSaundersonTable.hh" 78 #include "G4GoudsmitSaundersonTable.hh" 78 79 79 80 80 #include "G4PhysicalConstants.hh" 81 #include "G4PhysicalConstants.hh" 81 #include "Randomize.hh" 82 #include "Randomize.hh" 82 #include "G4Log.hh" 83 #include "G4Log.hh" 83 #include "G4Exp.hh" 84 #include "G4Exp.hh" 84 85 85 #include "G4GSMottCorrection.hh" 86 #include "G4GSMottCorrection.hh" 86 #include "G4MaterialTable.hh" 87 #include "G4MaterialTable.hh" 87 #include "G4Material.hh" 88 #include "G4Material.hh" 88 #include "G4MaterialCutsCouple.hh" 89 #include "G4MaterialCutsCouple.hh" 89 #include "G4ProductionCutsTable.hh" 90 #include "G4ProductionCutsTable.hh" 90 #include "G4EmParameters.hh" << 91 << 92 #include "G4String.hh" << 93 91 94 #include <fstream> 92 #include <fstream> 95 #include <cstdlib> 93 #include <cstdlib> 96 #include <cmath> 94 #include <cmath> 97 95 98 #include <iostream> 96 #include <iostream> 99 #include <iomanip> 97 #include <iomanip> 100 98 101 // perecomputed GS angular distributions, base 99 // perecomputed GS angular distributions, based on the Screened-Rutherford DCS 102 // are the same for e- and e+ so make sure we 100 // are the same for e- and e+ so make sure we load them only onece 103 G4bool G4GoudsmitSaundersonTable::gIsInitialis 101 G4bool G4GoudsmitSaundersonTable::gIsInitialised = false; 104 // 102 // 105 std::vector<G4GoudsmitSaundersonTable::GSMSCAn 103 std::vector<G4GoudsmitSaundersonTable::GSMSCAngularDtr*> G4GoudsmitSaundersonTable::gGSMSCAngularDistributions1; 106 std::vector<G4GoudsmitSaundersonTable::GSMSCAn 104 std::vector<G4GoudsmitSaundersonTable::GSMSCAngularDtr*> G4GoudsmitSaundersonTable::gGSMSCAngularDistributions2; 107 // 105 // 108 std::vector<double> G4GoudsmitSaundersonTable: 106 std::vector<double> G4GoudsmitSaundersonTable::gMoliereBc; 109 std::vector<double> G4GoudsmitSaundersonTable: 107 std::vector<double> G4GoudsmitSaundersonTable::gMoliereXc2; 110 108 111 109 112 G4GoudsmitSaundersonTable::G4GoudsmitSaunderso 110 G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable(G4bool iselectron) { 113 fIsElectron = iselectron; 111 fIsElectron = iselectron; 114 // set initial values: final values will be 112 // set initial values: final values will be set in the Initialize method 115 fLogLambda0 = 0.; // wi 113 fLogLambda0 = 0.; // will be set properly at init. 116 fLogDeltaLambda = 0.; // wi 114 fLogDeltaLambda = 0.; // will be set properly at init. 117 fInvLogDeltaLambda = 0.; // wi 115 fInvLogDeltaLambda = 0.; // will be set properly at init. 118 fInvDeltaQ1 = 0.; // wi 116 fInvDeltaQ1 = 0.; // will be set properly at init. 119 fDeltaQ2 = 0.; // wi 117 fDeltaQ2 = 0.; // will be set properly at init. 120 fInvDeltaQ2 = 0.; // wi 118 fInvDeltaQ2 = 0.; // will be set properly at init. 121 // 119 // 122 fLowEnergyLimit = 0.1*CLHEP::keV; // w 120 fLowEnergyLimit = 0.1*CLHEP::keV; // will be set properly at init. 123 fHighEnergyLimit = 100.0*CLHEP::MeV; // w 121 fHighEnergyLimit = 100.0*CLHEP::MeV; // will be set properly at init. 124 // 122 // 125 fIsMottCorrection = false; // w 123 fIsMottCorrection = false; // will be set properly at init. 126 fIsPWACorrection = false; // w 124 fIsPWACorrection = false; // will be set properly at init. 127 fMottCorrection = nullptr; 125 fMottCorrection = nullptr; 128 // 126 // 129 fNumSPCEbinPerDec = 3; 127 fNumSPCEbinPerDec = 3; 130 } 128 } 131 129 132 G4GoudsmitSaundersonTable::~G4GoudsmitSaunders 130 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable() { 133 for (std::size_t i=0; i<gGSMSCAngularDistrib << 131 for (size_t i=0; i<gGSMSCAngularDistributions1.size(); ++i) { 134 if (gGSMSCAngularDistributions1[i]) { 132 if (gGSMSCAngularDistributions1[i]) { 135 delete [] gGSMSCAngularDistributions1[i] 133 delete [] gGSMSCAngularDistributions1[i]->fUValues; 136 delete [] gGSMSCAngularDistributions1[i] 134 delete [] gGSMSCAngularDistributions1[i]->fParamA; 137 delete [] gGSMSCAngularDistributions1[i] 135 delete [] gGSMSCAngularDistributions1[i]->fParamB; 138 delete gGSMSCAngularDistributions1[i]; 136 delete gGSMSCAngularDistributions1[i]; 139 } 137 } 140 } 138 } 141 gGSMSCAngularDistributions1.clear(); 139 gGSMSCAngularDistributions1.clear(); 142 for (std::size_t i=0; i<gGSMSCAngularDistrib << 140 for (size_t i=0; i<gGSMSCAngularDistributions2.size(); ++i) { 143 if (gGSMSCAngularDistributions2[i]) { 141 if (gGSMSCAngularDistributions2[i]) { 144 delete [] gGSMSCAngularDistributions2[i] 142 delete [] gGSMSCAngularDistributions2[i]->fUValues; 145 delete [] gGSMSCAngularDistributions2[i] 143 delete [] gGSMSCAngularDistributions2[i]->fParamA; 146 delete [] gGSMSCAngularDistributions2[i] 144 delete [] gGSMSCAngularDistributions2[i]->fParamB; 147 delete gGSMSCAngularDistributions2[i]; 145 delete gGSMSCAngularDistributions2[i]; 148 } 146 } 149 } 147 } 150 gGSMSCAngularDistributions2.clear(); 148 gGSMSCAngularDistributions2.clear(); 151 if (fMottCorrection) { 149 if (fMottCorrection) { 152 delete fMottCorrection; 150 delete fMottCorrection; 153 fMottCorrection = nullptr; 151 fMottCorrection = nullptr; 154 } 152 } 155 // clear scp correction data 153 // clear scp correction data 156 for (std::size_t imc=0; imc<fSCPCPerMatCuts. << 154 for (size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) { 157 if (fSCPCPerMatCuts[imc]) { 155 if (fSCPCPerMatCuts[imc]) { 158 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 156 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 159 delete fSCPCPerMatCuts[imc]; 157 delete fSCPCPerMatCuts[imc]; 160 } 158 } 161 } 159 } 162 fSCPCPerMatCuts.clear(); 160 fSCPCPerMatCuts.clear(); 163 // 161 // 164 gIsInitialised = false; 162 gIsInitialised = false; 165 } 163 } 166 164 167 void G4GoudsmitSaundersonTable::Initialise(G4d 165 void G4GoudsmitSaundersonTable::Initialise(G4double lownergylimit, G4double highenergylimit) { 168 fLowEnergyLimit = lownergylimit; 166 fLowEnergyLimit = lownergylimit; 169 fHighEnergyLimit = highenergylimit; 167 fHighEnergyLimit = highenergylimit; 170 G4double lLambdaMin = G4Log(gLAMBMIN); 168 G4double lLambdaMin = G4Log(gLAMBMIN); 171 G4double lLambdaMax = G4Log(gLAMBMAX); 169 G4double lLambdaMax = G4Log(gLAMBMAX); 172 fLogLambda0 = lLambdaMin; 170 fLogLambda0 = lLambdaMin; 173 fLogDeltaLambda = (lLambdaMax-lLambdaMin 171 fLogDeltaLambda = (lLambdaMax-lLambdaMin)/(gLAMBNUM-1.); 174 fInvLogDeltaLambda = 1./fLogDeltaLambda; 172 fInvLogDeltaLambda = 1./fLogDeltaLambda; 175 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(g 173 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(gQNUM1-1.)); 176 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM 174 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM2-1.); 177 fInvDeltaQ2 = 1./fDeltaQ2; 175 fInvDeltaQ2 = 1./fDeltaQ2; 178 // load precomputed angular distributions an 176 // load precomputed angular distributions and set up several values used during the sampling 179 // these are particle independet => they go 177 // these are particle independet => they go to static container: load them only onece 180 if (!gIsInitialised) { 178 if (!gIsInitialised) { 181 // load pre-computed GS angular distributi 179 // load pre-computed GS angular distributions (computed based on Screened-Rutherford DCS) 182 LoadMSCData(); 180 LoadMSCData(); 183 gIsInitialised = true; 181 gIsInitialised = true; 184 } 182 } 185 InitMoliereMSCParams(); 183 InitMoliereMSCParams(); 186 // Mott-correction: particle(e- or e+) depen 184 // Mott-correction: particle(e- or e+) dependet so init them 187 if (fIsMottCorrection) { 185 if (fIsMottCorrection) { 188 if (!fMottCorrection) { 186 if (!fMottCorrection) { 189 fMottCorrection = new G4GSMottCorrection 187 fMottCorrection = new G4GSMottCorrection(fIsElectron); 190 } 188 } 191 fMottCorrection->Initialise(); 189 fMottCorrection->Initialise(); 192 } 190 } 193 // init scattering power correction data; us 191 // init scattering power correction data; used only together with Mott-correction 194 // (Moliere's parameters must be initialised 192 // (Moliere's parameters must be initialised before) 195 if (fMottCorrection) { 193 if (fMottCorrection) { 196 InitSCPCorrection(); 194 InitSCPCorrection(); 197 } 195 } 198 } 196 } 199 197 200 198 201 // samplig multiple scattering angles cos(thet 199 // samplig multiple scattering angles cos(theta) and sin(thata) 202 // - including no-scattering, single, "few" s 200 // - including no-scattering, single, "few" scattering cases as well 203 // - Mott-correction will be included if it w 201 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true) 204 // lambdaval : s/lambda_el 202 // lambdaval : s/lambda_el 205 // qval : s/lambda_el G1 203 // qval : s/lambda_el G1 206 // scra : screening parameter 204 // scra : screening parameter 207 // cost : will be the smapled cos(theta) 205 // cost : will be the smapled cos(theta) 208 // sint : will be the smapled sin(theta) 206 // sint : will be the smapled sin(theta) 209 // lekin : logarithm of the current kineti 207 // lekin : logarithm of the current kinetic energy 210 // beta2 : the corresponding beta square 208 // beta2 : the corresponding beta square 211 // matindx : index of the current material 209 // matindx : index of the current material 212 // returns true if it was msc 210 // returns true if it was msc 213 G4bool G4GoudsmitSaundersonTable::Sampling(G4d 211 G4bool G4GoudsmitSaundersonTable::Sampling(G4double lambdaval, G4double qval, G4double scra, G4double &cost, 214 G4d 212 G4double &sint, G4double lekin, G4double beta2, G4int matindx, 215 GSM 213 GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti, 216 G4d 214 G4double &transfPar, G4bool isfirst) { 217 G4double rand0 = G4UniformRand(); 215 G4double rand0 = G4UniformRand(); 218 G4double expn = G4Exp(-lambdaval); 216 G4double expn = G4Exp(-lambdaval); 219 // 217 // 220 // no scattering case 218 // no scattering case 221 if (rand0<expn) { 219 if (rand0<expn) { 222 cost = 1.0; 220 cost = 1.0; 223 sint = 0.0; 221 sint = 0.0; 224 return false; 222 return false; 225 } 223 } 226 // 224 // 227 // single scattering case : sample from the 225 // single scattering case : sample from the single scattering PDF 228 // - Mott-correction will be included if it 226 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true) 229 if (rand0<(1.+lambdaval)*expn) { 227 if (rand0<(1.+lambdaval)*expn) { 230 // cost is sampled in SingleScattering() 228 // cost is sampled in SingleScattering() 231 cost = SingleScattering(lambdaval, scra, l 229 cost = SingleScattering(lambdaval, scra, lekin, beta2, matindx); 232 // add protections 230 // add protections 233 if (cost<-1.0) cost = -1.0; 231 if (cost<-1.0) cost = -1.0; 234 if (cost>1.0) cost = 1.0; 232 if (cost>1.0) cost = 1.0; 235 // compute sin(theta) from the sampled cos 233 // compute sin(theta) from the sampled cos(theta) 236 G4double dum0 = 1.-cost; 234 G4double dum0 = 1.-cost; 237 sint = std::sqrt(dum0*(2.0-dum0)); 235 sint = std::sqrt(dum0*(2.0-dum0)); 238 return false; 236 return false; 239 } 237 } 240 // 238 // 241 // handle this case: 239 // handle this case: 242 // -lambdaval < 1 i.e. mean #elastic ev 240 // -lambdaval < 1 i.e. mean #elastic events along the step is < 1 but 243 // the currently sampled case is not 0 241 // the currently sampled case is not 0 or 1 scattering. [Our minimal 244 // lambdaval (that we have precomputed 242 // lambdaval (that we have precomputed, transformed angular distributions 245 // stored in a form of equally probabe 243 // stored in a form of equally probabe intervalls together with rational 246 // interp. parameters) is 1.] 244 // interp. parameters) is 1.] 247 // -probability of having n elastic eve 245 // -probability of having n elastic events follows Poisson stat. with 248 // lambdaval parameter. 246 // lambdaval parameter. 249 // -the max. probability (when lambdava 247 // -the max. probability (when lambdaval=1) of having more than one 250 // elastic events is 0.2642411 and the 248 // elastic events is 0.2642411 and the prob of having 2,3,..,n elastic 251 // events decays rapidly with n. So se 249 // events decays rapidly with n. So set a max n to 10. 252 // -sampling of this cases is done in a 250 // -sampling of this cases is done in a one-by-one single elastic event way 253 // where the current #elastic event is 251 // where the current #elastic event is sampled from the Poisson distr. 254 if (lambdaval<1.0) { 252 if (lambdaval<1.0) { 255 G4double prob, cumprob; 253 G4double prob, cumprob; 256 prob = cumprob = expn; 254 prob = cumprob = expn; 257 G4double curcost,cursint; 255 G4double curcost,cursint; 258 // init cos(theta) and sin(theta) to the z 256 // init cos(theta) and sin(theta) to the zero scattering values 259 cost = 1.0; 257 cost = 1.0; 260 sint = 0.0; 258 sint = 0.0; 261 for (G4int iel=1; iel<10; ++iel) { 259 for (G4int iel=1; iel<10; ++iel) { 262 // prob of having iel scattering from Po 260 // prob of having iel scattering from Poisson 263 prob *= lambdaval/(G4double)iel; 261 prob *= lambdaval/(G4double)iel; 264 cumprob += prob; 262 cumprob += prob; 265 // 263 // 266 //sample cos(theta) from the singe scatt 264 //sample cos(theta) from the singe scattering pdf: 267 // - Mott-correction will be included if 265 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true) 268 curcost = SingleScattering(lambdav 266 curcost = SingleScattering(lambdaval, scra, lekin, beta2, matindx); 269 G4double dum0 = 1.-curcost; 267 G4double dum0 = 1.-curcost; 270 cursint = dum0*(2.0-dum0); // sin^ 268 cursint = dum0*(2.0-dum0); // sin^2(theta) 271 // 269 // 272 // if we got current deflection that is 270 // if we got current deflection that is not too small 273 // then update cos(theta) sin(theta) 271 // then update cos(theta) sin(theta) 274 if (cursint>1.0e-20) { 272 if (cursint>1.0e-20) { 275 cursint = std::sqrt(cursint); 273 cursint = std::sqrt(cursint); 276 G4double curphi = CLHEP::twopi*G4Unifo 274 G4double curphi = CLHEP::twopi*G4UniformRand(); 277 cost = cost*curcost-sint*cu 275 cost = cost*curcost-sint*cursint*std::cos(curphi); 278 sint = std::sqrt(std::max(0 276 sint = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost))); 279 } 277 } 280 // 278 // 281 // check if we have done enough scatteri 279 // check if we have done enough scattering i.e. sampling from the Poisson 282 if (rand0<cumprob) { 280 if (rand0<cumprob) { 283 return false; 281 return false; 284 } 282 } 285 } 283 } 286 // if reached the max iter i.e. 10 284 // if reached the max iter i.e. 10 287 return false; 285 return false; 288 } 286 } 289 // 287 // 290 // multiple scattering case with lambdavalue 288 // multiple scattering case with lambdavalue >= 1: 291 // - use the precomputed and transformed G 289 // - use the precomputed and transformed Goudsmit-Saunderson angular 292 // distributions to sample cos(theta) 290 // distributions to sample cos(theta) 293 // - Mott-correction will be included if i 291 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true) 294 cost = SampleCosTheta(lambdaval, qval, scra, 292 cost = SampleCosTheta(lambdaval, qval, scra, lekin, beta2, matindx, gsDtr, mcekini, mcdelti, transfPar, isfirst); 295 // add protections 293 // add protections 296 if (cost<-1.0) cost = -1.0; 294 if (cost<-1.0) cost = -1.0; 297 if (cost> 1.0) cost = 1.0; 295 if (cost> 1.0) cost = 1.0; 298 // compute cos(theta) and sin(theta) from th 296 // compute cos(theta) and sin(theta) from the sampled 1-cos(theta) 299 G4double dum0 = 1.0-cost; 297 G4double dum0 = 1.0-cost; 300 sint = std::sqrt(dum0*(2.0-dum0)); 298 sint = std::sqrt(dum0*(2.0-dum0)); 301 // return true if it was msc 299 // return true if it was msc 302 return true; 300 return true; 303 } 301 } 304 302 305 303 306 G4double G4GoudsmitSaundersonTable::SampleCosT 304 G4double G4GoudsmitSaundersonTable::SampleCosTheta(G4double lambdaval, G4double qval, G4double scra, 307 305 G4double lekin, G4double beta2, G4int matindx, 308 306 GSMSCAngularDtr **gsDtr, G4int &mcekini,G4int &mcdelti, 309 307 G4double &transfPar, G4bool isfirst) { 310 G4double cost = 1.; 308 G4double cost = 1.; 311 // determine the base GS angular distributio 309 // determine the base GS angular distribution if it is the first call (when sub-step sampling is used) 312 if (isfirst) { 310 if (isfirst) { 313 *gsDtr = GetGSAngularDtr(scra, lambdaval, 311 *gsDtr = GetGSAngularDtr(scra, lambdaval, qval, transfPar); 314 } 312 } 315 // sample cost from the GS angular distribut 313 // sample cost from the GS angular distribution (computed based on Screened-Rutherford DCS) 316 cost = SampleGSSRCosTheta(*gsDtr, transfPar) 314 cost = SampleGSSRCosTheta(*gsDtr, transfPar); 317 // Mott-correction if it was requested by th 315 // Mott-correction if it was requested by the user 318 if (fIsMottCorrection && *gsDtr) { // no 316 if (fIsMottCorrection && *gsDtr) { // no Mott-correction in case of izotropic theta 319 static const G4int nlooplim = 1000; 317 static const G4int nlooplim = 1000; 320 G4int nloop = 0 ; // rejection loop 318 G4int nloop = 0 ; // rejection loop counter 321 // G4int ekindx = -1; // evaluate only 319 // G4int ekindx = -1; // evaluate only in the first call 322 // G4int deltindx = -1 ; // evaluate onl 320 // G4int deltindx = -1 ; // evaluate only in the first call 323 G4double val = fMottCorrection->GetMo 321 G4double val = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti); 324 while (G4UniformRand()>val && ++nloop<nloo 322 while (G4UniformRand()>val && ++nloop<nlooplim) { 325 // sampling cos(theta) 323 // sampling cos(theta) 326 cost = SampleGSSRCosTheta(*gsDtr, transf 324 cost = SampleGSSRCosTheta(*gsDtr, transfPar); 327 val = fMottCorrection->GetMottRejection 325 val = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti); 328 }; 326 }; 329 } 327 } 330 return cost; 328 return cost; 331 } 329 } 332 330 333 331 334 // returns with cost sampled from the GS angul 332 // returns with cost sampled from the GS angular distribution computed based on Screened-Rutherford DCS 335 G4double G4GoudsmitSaundersonTable::SampleGSSR 333 G4double G4GoudsmitSaundersonTable::SampleGSSRCosTheta(const GSMSCAngularDtr *gsDtr, G4double transfpar) { 336 // check if isotropic theta (i.e. cost is un 334 // check if isotropic theta (i.e. cost is uniform on [-1:1]) 337 if (!gsDtr) { 335 if (!gsDtr) { 338 return 1.-2.0*G4UniformRand(); 336 return 1.-2.0*G4UniformRand(); 339 } 337 } 340 // 338 // 341 // sampling form the selected distribution 339 // sampling form the selected distribution 342 G4double ndatm1 = gsDtr->fNumData-1.; 340 G4double ndatm1 = gsDtr->fNumData-1.; 343 G4double delta = 1.0/ndatm1; 341 G4double delta = 1.0/ndatm1; 344 // determine lower cumulative bin inidex 342 // determine lower cumulative bin inidex 345 G4double rndm = G4UniformRand(); 343 G4double rndm = G4UniformRand(); 346 G4int indxl = rndm*ndatm1; 344 G4int indxl = rndm*ndatm1; 347 G4double aval = rndm-indxl*delta; 345 G4double aval = rndm-indxl*delta; 348 G4double dum0 = delta*aval; 346 G4double dum0 = delta*aval; 349 347 350 G4double dum1 = (1.0+gsDtr->fParamA[indxl] 348 G4double dum1 = (1.0+gsDtr->fParamA[indxl]+gsDtr->fParamB[indxl])*dum0; 351 G4double dum2 = delta*delta + gsDtr->fPara 349 G4double dum2 = delta*delta + gsDtr->fParamA[indxl]*dum0 + gsDtr->fParamB[indxl]*aval*aval; 352 G4double sample = gsDtr->fUValues[indxl] + 350 G4double sample = gsDtr->fUValues[indxl] + dum1/dum2 *(gsDtr->fUValues[indxl+1]-gsDtr->fUValues[indxl]); 353 // transform back u to cos(theta) : 351 // transform back u to cos(theta) : 354 // this is the sampled cos(theta) = (2.0*par 352 // this is the sampled cos(theta) = (2.0*para*sample)/(1.0-sample+para) 355 return 1.-(2.0*transfpar*sample)/(1.0-sample 353 return 1.-(2.0*transfpar*sample)/(1.0-sample+transfpar); 356 } 354 } 357 355 358 356 359 // determine the GS angular distribution we ne 357 // determine the GS angular distribution we need to sample from: will set other things as well ... 360 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4 358 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4GoudsmitSaundersonTable::GetGSAngularDtr(G4double scra, 361 359 G4double &lambdaval, G4double &qval, G4double &transfpar) { 362 GSMSCAngularDtr *dtr = nullptr; 360 GSMSCAngularDtr *dtr = nullptr; 363 G4bool first = false; 361 G4bool first = false; 364 // isotropic cost above gQMAX2 (i.e. dtr sta 362 // isotropic cost above gQMAX2 (i.e. dtr stays nullptr) 365 if (qval<gQMAX2) { 363 if (qval<gQMAX2) { 366 G4int lamIndx = -1; // lambda value in 364 G4int lamIndx = -1; // lambda value index 367 G4int qIndx = -1; // lambda value in 365 G4int qIndx = -1; // lambda value index 368 // init to second grid Q values 366 // init to second grid Q values 369 G4int numQVal = gQNUM2; 367 G4int numQVal = gQNUM2; 370 G4double minQVal = gQMIN2; 368 G4double minQVal = gQMIN2; 371 G4double invDelQ = fInvDeltaQ2; 369 G4double invDelQ = fInvDeltaQ2; 372 G4double pIndxH = 0.; // probability of 370 G4double pIndxH = 0.; // probability of taking higher index 373 // check if first or second grid needs to 371 // check if first or second grid needs to be used 374 if (qval<gQMIN2) { // first grid 372 if (qval<gQMIN2) { // first grid 375 first = true; 373 first = true; 376 // protect against qval<gQMIN1 374 // protect against qval<gQMIN1 377 if (qval<gQMIN1) { 375 if (qval<gQMIN1) { 378 qval = gQMIN1; 376 qval = gQMIN1; 379 qIndx = 0; 377 qIndx = 0; 380 //pIndxH = 0.; 378 //pIndxH = 0.; 381 } 379 } 382 // set to first grid Q values 380 // set to first grid Q values 383 numQVal = gQNUM1; 381 numQVal = gQNUM1; 384 minQVal = gQMIN1; 382 minQVal = gQMIN1; 385 invDelQ = fInvDeltaQ1; 383 invDelQ = fInvDeltaQ1; 386 } 384 } 387 // make sure that lambda = s/lambda_el is 385 // make sure that lambda = s/lambda_el is in [gLAMBMIN,gLAMBMAX) 388 // lambda<gLAMBMIN=1 is already handeled b 386 // lambda<gLAMBMIN=1 is already handeled before so lambda>= gLAMBMIN for sure 389 if (lambdaval>=gLAMBMAX) { 387 if (lambdaval>=gLAMBMAX) { 390 lambdaval = gLAMBMAX-1.e-8; 388 lambdaval = gLAMBMAX-1.e-8; 391 lamIndx = gLAMBNUM-1; 389 lamIndx = gLAMBNUM-1; 392 } 390 } 393 G4double lLambda = G4Log(lambdaval); 391 G4double lLambda = G4Log(lambdaval); 394 // 392 // 395 // determine lower lambda (=s/lambda_el) i 393 // determine lower lambda (=s/lambda_el) index: linear interp. on log(lambda) scale 396 if (lamIndx<0) { 394 if (lamIndx<0) { 397 pIndxH = (lLambda-fLogLambda0)*fInvLogD 395 pIndxH = (lLambda-fLogLambda0)*fInvLogDeltaLambda; 398 lamIndx = (G4int)(pIndxH); // low 396 lamIndx = (G4int)(pIndxH); // lower index of the lambda bin 399 pIndxH = pIndxH-lamIndx; // proba 397 pIndxH = pIndxH-lamIndx; // probability of taking the higher index distribution 400 if (G4UniformRand()<pIndxH) { 398 if (G4UniformRand()<pIndxH) { 401 ++lamIndx; 399 ++lamIndx; 402 } 400 } 403 } 401 } 404 // 402 // 405 // determine lower Q (=s/lambda_el G1) ind 403 // determine lower Q (=s/lambda_el G1) index: linear interp. on Q 406 if (qIndx<0) { 404 if (qIndx<0) { 407 pIndxH = (qval-minQVal)*invDelQ; 405 pIndxH = (qval-minQVal)*invDelQ; 408 qIndx = (G4int)(pIndxH); // low 406 qIndx = (G4int)(pIndxH); // lower index of the Q bin 409 pIndxH = pIndxH-qIndx; 407 pIndxH = pIndxH-qIndx; 410 if (G4UniformRand()<pIndxH) { 408 if (G4UniformRand()<pIndxH) { 411 ++qIndx; 409 ++qIndx; 412 } 410 } 413 } 411 } 414 // set indx 412 // set indx 415 G4int indx = lamIndx*numQVal+qIndx; 413 G4int indx = lamIndx*numQVal+qIndx; 416 if (first) { 414 if (first) { 417 dtr = gGSMSCAngularDistributions1[indx]; 415 dtr = gGSMSCAngularDistributions1[indx]; 418 } else { 416 } else { 419 dtr = gGSMSCAngularDistributions2[indx]; 417 dtr = gGSMSCAngularDistributions2[indx]; 420 } 418 } 421 // dtr might be nullptr that indicates iso 419 // dtr might be nullptr that indicates isotropic cot distribution because: 422 // - if the selected lamIndx, qIndx corres 420 // - if the selected lamIndx, qIndx correspond to L(=s/lambda_el) and Q(=s/lambda_el G1) such that G1(=Q/L) > 1 423 // G1 should always be < 1 and if G1 is 421 // G1 should always be < 1 and if G1 is ~1 -> the dtr is isotropic (this can only happen in case of the 2. grid) 424 // 422 // 425 // compute the transformation parameter 423 // compute the transformation parameter 426 if (lambdaval>10.0) { 424 if (lambdaval>10.0) { 427 transfpar = 0.5*(-2.77164+lLambda*( 2.94 425 transfpar = 0.5*(-2.77164+lLambda*( 2.94874-lLambda*(0.1535754-lLambda*0.00552888) )); 428 } else { 426 } else { 429 transfpar = 0.5*(1.347+lLambda*(0.209364 427 transfpar = 0.5*(1.347+lLambda*(0.209364-lLambda*(0.45525-lLambda*(0.50142-lLambda*0.081234)))); 430 } 428 } 431 transfpar *= (lambdaval+4.0)*scra; 429 transfpar *= (lambdaval+4.0)*scra; 432 } 430 } 433 // return with the selected GS angular distr 431 // return with the selected GS angular distribution that we need to sample cost from (if nullptr => isotropic cost) 434 return dtr; 432 return dtr; 435 } 433 } 436 434 437 435 438 void G4GoudsmitSaundersonTable::LoadMSCData() 436 void G4GoudsmitSaundersonTable::LoadMSCData() { >> 437 char* path = getenv("G4LEDATA"); >> 438 if (!path) { >> 439 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006", >> 440 FatalException, >> 441 "Environment variable G4LEDATA not defined"); >> 442 return; >> 443 } >> 444 // 439 gGSMSCAngularDistributions1.resize(gLAMBNUM* 445 gGSMSCAngularDistributions1.resize(gLAMBNUM*gQNUM1,nullptr); 440 const G4String str1 = G4EmParameters::Instan << 441 for (G4int il=0; il<gLAMBNUM; ++il) { 446 for (G4int il=0; il<gLAMBNUM; ++il) { 442 G4String fname = str1 + std::to_string(il) << 447 char fname[512]; >> 448 sprintf(fname,"%s/msc_GS/GSGrid_1/gsDistr_%d",path,il); 443 std::ifstream infile(fname,std::ios::in); 449 std::ifstream infile(fname,std::ios::in); 444 if (!infile.is_open()) { 450 if (!infile.is_open()) { 445 G4String msgc = "Cannot open file: " + f << 451 char msgc[512]; >> 452 sprintf(msgc,"Cannot open file: %s .",fname); 446 G4Exception("G4GoudsmitSaundersonTable:: 453 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006", 447 FatalException, msgc.c_str()); << 454 FatalException, msgc); 448 return; 455 return; 449 } 456 } 450 for (G4int iq=0; iq<gQNUM1; ++iq) { 457 for (G4int iq=0; iq<gQNUM1; ++iq) { 451 auto gsd = new GSMSCAngularDtr(); << 458 GSMSCAngularDtr *gsd = new GSMSCAngularDtr(); 452 infile >> gsd->fNumData; 459 infile >> gsd->fNumData; 453 gsd->fUValues = new G4double[gsd->fNumDa 460 gsd->fUValues = new G4double[gsd->fNumData](); 454 gsd->fParamA = new G4double[gsd->fNumDa 461 gsd->fParamA = new G4double[gsd->fNumData](); 455 gsd->fParamB = new G4double[gsd->fNumDa 462 gsd->fParamB = new G4double[gsd->fNumData](); 456 G4double ddummy; 463 G4double ddummy; 457 infile >> ddummy; infile >> ddummy; 464 infile >> ddummy; infile >> ddummy; 458 for (G4int i=0; i<gsd->fNumData; ++i) { 465 for (G4int i=0; i<gsd->fNumData; ++i) { 459 infile >> gsd->fUValues[i]; 466 infile >> gsd->fUValues[i]; 460 infile >> gsd->fParamA[i]; 467 infile >> gsd->fParamA[i]; 461 infile >> gsd->fParamB[i]; 468 infile >> gsd->fParamB[i]; 462 } 469 } 463 gGSMSCAngularDistributions1[il*gQNUM1+iq 470 gGSMSCAngularDistributions1[il*gQNUM1+iq] = gsd; 464 } 471 } 465 infile.close(); 472 infile.close(); 466 } 473 } 467 // 474 // 468 // second grid 475 // second grid 469 gGSMSCAngularDistributions2.resize(gLAMBNUM* 476 gGSMSCAngularDistributions2.resize(gLAMBNUM*gQNUM2,nullptr); 470 const G4String str2 = G4EmParameters::Instan << 471 for (G4int il=0; il<gLAMBNUM; ++il) { 477 for (G4int il=0; il<gLAMBNUM; ++il) { 472 G4String fname = str2 + std::to_string(il) << 478 char fname[512]; >> 479 sprintf(fname,"%s/msc_GS/GSGrid_2/gsDistr_%d",path,il); 473 std::ifstream infile(fname,std::ios::in); 480 std::ifstream infile(fname,std::ios::in); 474 if (!infile.is_open()) { 481 if (!infile.is_open()) { 475 G4String msgc = "Cannot open file: " + f << 482 char msgc[512]; >> 483 sprintf(msgc,"Cannot open file: %s .",fname); 476 G4Exception("G4GoudsmitSaundersonTable:: 484 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006", 477 FatalException, msgc.c_str()); << 485 FatalException, msgc); 478 return; 486 return; 479 } 487 } 480 for (G4int iq=0; iq<gQNUM2; ++iq) { 488 for (G4int iq=0; iq<gQNUM2; ++iq) { 481 G4int numData; 489 G4int numData; 482 infile >> numData; 490 infile >> numData; 483 if (numData>1) { 491 if (numData>1) { 484 auto gsd = new GSMSCAngularDtr(); << 492 GSMSCAngularDtr *gsd = new GSMSCAngularDtr(); 485 gsd->fNumData = numData; 493 gsd->fNumData = numData; 486 gsd->fUValues = new G4double[gsd->fNum 494 gsd->fUValues = new G4double[gsd->fNumData](); 487 gsd->fParamA = new G4double[gsd->fNum 495 gsd->fParamA = new G4double[gsd->fNumData](); 488 gsd->fParamB = new G4double[gsd->fNum 496 gsd->fParamB = new G4double[gsd->fNumData](); 489 double ddummy; 497 double ddummy; 490 infile >> ddummy; infile >> ddummy; 498 infile >> ddummy; infile >> ddummy; 491 for (G4int i=0; i<gsd->fNumData; ++i) 499 for (G4int i=0; i<gsd->fNumData; ++i) { 492 infile >> gsd->fUValues[i]; 500 infile >> gsd->fUValues[i]; 493 infile >> gsd->fParamA[i]; 501 infile >> gsd->fParamA[i]; 494 infile >> gsd->fParamB[i]; 502 infile >> gsd->fParamB[i]; 495 } 503 } 496 gGSMSCAngularDistributions2[il*gQNUM2+ 504 gGSMSCAngularDistributions2[il*gQNUM2+iq] = gsd; 497 } else { 505 } else { 498 gGSMSCAngularDistributions2[il*gQNUM2+ 506 gGSMSCAngularDistributions2[il*gQNUM2+iq] = nullptr; 499 } 507 } 500 } 508 } 501 infile.close(); 509 infile.close(); 502 } 510 } 503 } 511 } 504 512 505 // samples cost in single scattering based on 513 // samples cost in single scattering based on Screened-Rutherford DCS 506 // (with Mott-correction if it was requested) 514 // (with Mott-correction if it was requested) 507 G4double G4GoudsmitSaundersonTable::SingleScat 515 G4double G4GoudsmitSaundersonTable::SingleScattering(G4double /*lambdaval*/, G4double scra, 508 516 G4double lekin, G4double beta2, 509 517 G4int matindx) { 510 G4double rand1 = G4UniformRand(); 518 G4double rand1 = G4UniformRand(); 511 // sample cost from the Screened-Rutherford 519 // sample cost from the Screened-Rutherford DCS 512 G4double cost = 1.-2.0*scra*rand1/(1.0-rand 520 G4double cost = 1.-2.0*scra*rand1/(1.0-rand1+scra); 513 // Mott-correction if it was requested by th 521 // Mott-correction if it was requested by the user 514 if (fIsMottCorrection) { 522 if (fIsMottCorrection) { 515 static const G4int nlooplim = 1000; // re 523 static const G4int nlooplim = 1000; // rejection loop limit 516 G4int nloop = 0 ; // loop counter 524 G4int nloop = 0 ; // loop counter 517 G4int ekindx = -1 ; // evaluate only 525 G4int ekindx = -1 ; // evaluate only in the first call 518 G4int deltindx = 0 ; // single scatter 526 G4int deltindx = 0 ; // single scattering case 519 G4double q1 = 0.; // not used when 527 G4double q1 = 0.; // not used when deltindx = 0; 520 // computing Mott rejection function value 528 // computing Mott rejection function value 521 G4double val = fMottCorrection->GetMo 529 G4double val = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost, 522 530 matindx, ekindx, deltindx); 523 while (G4UniformRand()>val && ++nloop<nloo 531 while (G4UniformRand()>val && ++nloop<nlooplim) { 524 // sampling cos(theta) from the Screened 532 // sampling cos(theta) from the Screened-Rutherford DCS 525 rand1 = G4UniformRand(); 533 rand1 = G4UniformRand(); 526 cost = 1.-2.0*scra*rand1/(1.0-rand1+scr 534 cost = 1.-2.0*scra*rand1/(1.0-rand1+scra); 527 // computing Mott rejection function val 535 // computing Mott rejection function value 528 val = fMottCorrection->GetMottRejectio 536 val = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost, matindx, 529 537 ekindx, deltindx); 530 }; 538 }; 531 } 539 } 532 return cost; 540 return cost; 533 } 541 } 534 542 535 543 536 void G4GoudsmitSaundersonTable::GetMottCorrec 544 void G4GoudsmitSaundersonTable::GetMottCorrectionFactors(G4double logekin, G4double beta2, 537 545 G4int matindx, G4double &mcToScr, 538 546 G4double &mcToQ1, G4double &mcToG2PerG1) { 539 if (fIsMottCorrection) { 547 if (fIsMottCorrection) { 540 fMottCorrection->GetMottCorrectionFactors( 548 fMottCorrection->GetMottCorrectionFactors(logekin, beta2, matindx, mcToScr, mcToQ1, mcToG2PerG1); 541 } 549 } 542 } 550 } 543 551 544 552 545 // compute material dependent Moliere MSC para 553 // compute material dependent Moliere MSC parameters at initialisation 546 void G4GoudsmitSaundersonTable::InitMoliereMSC 554 void G4GoudsmitSaundersonTable::InitMoliereMSCParams() { 547 const G4double const1 = 7821.6; // [ 555 const G4double const1 = 7821.6; // [cm2/g] 548 const G4double const2 = 0.1569; // [ 556 const G4double const2 = 0.1569; // [cm2 MeV2 / g] 549 const G4double finstrc2 = 5.325135453E-5; / 557 const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square 550 558 551 G4MaterialTable* theMaterialTable = G4Mater 559 G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 552 // get number of materials in the table 560 // get number of materials in the table 553 std::size_t numMaterials = theMaterialTable << 561 size_t numMaterials = theMaterialTable->size(); 554 // make sure that we have long enough vecto 562 // make sure that we have long enough vectors 555 if(gMoliereBc.size()<numMaterials) { 563 if(gMoliereBc.size()<numMaterials) { 556 gMoliereBc.resize(numMaterials); 564 gMoliereBc.resize(numMaterials); 557 gMoliereXc2.resize(numMaterials); 565 gMoliereXc2.resize(numMaterials); 558 } 566 } 559 G4double xi = 1.0; 567 G4double xi = 1.0; 560 G4int maxZ = 200; 568 G4int maxZ = 200; 561 if (fIsMottCorrection || fIsPWACorrection) 569 if (fIsMottCorrection || fIsPWACorrection) { 562 // xi = 1.0; <= always set to 1 from n 570 // xi = 1.0; <= always set to 1 from now on 563 maxZ = G4GSMottCorrection::GetMaxZet(); 571 maxZ = G4GSMottCorrection::GetMaxZet(); 564 } 572 } 565 // 573 // 566 for (std::size_t imat=0; imat<numMaterials; << 574 for (size_t imat=0; imat<numMaterials; ++imat) { 567 const G4Material* theMaterial = 575 const G4Material* theMaterial = (*theMaterialTable)[imat]; 568 const G4ElementVector* theElemVect = 576 const G4ElementVector* theElemVect = theMaterial->GetElementVector(); 569 const G4int numelems = << 577 const G4int numelems = theMaterial->GetNumberOfElements(); 570 // 578 // 571 const G4double* theNbAtomsPerVolVe 579 const G4double* theNbAtomsPerVolVect = theMaterial->GetVecNbOfAtomsPerVolume(); 572 G4double theTotNbAtomsPerVo 580 G4double theTotNbAtomsPerVol = theMaterial->GetTotNbOfAtomsPerVolume(); 573 // 581 // 574 G4double zs = 0.0; 582 G4double zs = 0.0; 575 G4double zx = 0.0; 583 G4double zx = 0.0; 576 G4double ze = 0.0; 584 G4double ze = 0.0; 577 G4double sa = 0.0; 585 G4double sa = 0.0; 578 // 586 // 579 for(G4int ielem = 0; ielem < numelems; ie 587 for(G4int ielem = 0; ielem < numelems; ielem++) { 580 G4double zet = (*theElemVect)[ielem]->G 588 G4double zet = (*theElemVect)[ielem]->GetZ(); 581 if (zet>maxZ) { 589 if (zet>maxZ) { 582 zet = (G4double)maxZ; 590 zet = (G4double)maxZ; 583 } 591 } 584 G4double iwa = (*theElemVect)[ielem]-> 592 G4double iwa = (*theElemVect)[ielem]->GetN(); 585 G4double ipz = theNbAtomsPerVolVect[ie 593 G4double ipz = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol; 586 G4double dum = ipz*zet*(zet+xi); 594 G4double dum = ipz*zet*(zet+xi); 587 zs += dum; 595 zs += dum; 588 ze += dum*(-2.0/3.0)*G4Log(ze 596 ze += dum*(-2.0/3.0)*G4Log(zet); 589 zx += dum*G4Log(1.0+3.34*fins 597 zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet); 590 sa += ipz*iwa; 598 sa += ipz*iwa; 591 } 599 } 592 G4double density = theMaterial->GetDensit 600 G4double density = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3] 593 // 601 // 594 gMoliereBc[theMaterial->GetIndex()] = co 602 gMoliereBc[theMaterial->GetIndex()] = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm] 595 gMoliereXc2[theMaterial->GetIndex()] = co 603 gMoliereXc2[theMaterial->GetIndex()] = const2*density*zs/sa; // [MeV2/cm] 596 // change to Geant4 internal units of 1/l 604 // change to Geant4 internal units of 1/length and energ2/length 597 gMoliereBc[theMaterial->GetIndex()] *= 1 605 gMoliereBc[theMaterial->GetIndex()] *= 1.0/CLHEP::cm; 598 gMoliereXc2[theMaterial->GetIndex()] *= C 606 gMoliereXc2[theMaterial->GetIndex()] *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm; 599 } 607 } 600 } 608 } 601 609 602 610 603 // this method is temporary, will be removed/r 611 // this method is temporary, will be removed/replaced with a more effictien solution after 10.3.ref09 604 G4double G4GoudsmitSaundersonTable::ComputeSca 612 G4double G4GoudsmitSaundersonTable::ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut, G4double ekin) { 605 G4int imc = matcut->GetIndex(); 613 G4int imc = matcut->GetIndex(); 606 G4double corFactor = 1.0; 614 G4double corFactor = 1.0; 607 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin< 615 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) { 608 return corFactor; 616 return corFactor; 609 } 617 } 610 // get the scattering power correction facto 618 // get the scattering power correction factor 611 G4double lekin = G4Log(ekin); 619 G4double lekin = G4Log(ekin); 612 G4double remaining = (lekin-fSCPCPerMatCuts 620 G4double remaining = (lekin-fSCPCPerMatCuts[imc]->fLEmin)*fSCPCPerMatCuts[imc]->fILDel; 613 std::size_t lindx = (std::size_t)remaining << 621 G4int lindx = (G4int)remaining; 614 remaining -= lindx; 622 remaining -= lindx; 615 std::size_t imax = fSCPCPerMatCuts[imc]-> << 623 G4int imax = fSCPCPerMatCuts[imc]->fVSCPC.size()-1; 616 if (lindx>=imax) { 624 if (lindx>=imax) { 617 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i 625 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax]; 618 } else { 626 } else { 619 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l 627 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]); 620 } 628 } 621 return corFactor; 629 return corFactor; 622 } 630 } 623 631 624 632 625 void G4GoudsmitSaundersonTable::InitSCPCorrect 633 void G4GoudsmitSaundersonTable::InitSCPCorrection() { 626 // get the material-cuts table 634 // get the material-cuts table 627 G4ProductionCutsTable *thePCTable = G4Produc 635 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable(); 628 std::size_t numMatCuts = thePCTab << 636 size_t numMatCuts = thePCTable->GetTableSize(); 629 // clear container if any 637 // clear container if any 630 for (std::size_t imc=0; imc<fSCPCPerMatCuts. << 638 for (size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) { 631 if (fSCPCPerMatCuts[imc]) { 639 if (fSCPCPerMatCuts[imc]) { 632 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 640 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 633 delete fSCPCPerMatCuts[imc]; 641 delete fSCPCPerMatCuts[imc]; 634 fSCPCPerMatCuts[imc] = nullptr; 642 fSCPCPerMatCuts[imc] = nullptr; 635 } 643 } 636 } 644 } 637 // 645 // 638 // set size of the container and create the 646 // set size of the container and create the corresponding data structures 639 fSCPCPerMatCuts.resize(numMatCuts,nullptr); 647 fSCPCPerMatCuts.resize(numMatCuts,nullptr); 640 // loop over the material-cuts and create sc 648 // loop over the material-cuts and create scattering power correction data structure for each 641 for (G4int imc=0; imc<(G4int)numMatCuts; ++i << 649 for (size_t imc=0; imc<numMatCuts; ++imc) { 642 const G4MaterialCutsCouple *matCut = theP 650 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc); 643 // get e- production cut in the current ma 651 // get e- production cut in the current material-cuts in energy 644 G4double limit; 652 G4double limit; 645 G4double ecut; 653 G4double ecut; 646 if (fIsElectron) { 654 if (fIsElectron) { 647 ecut = (*(thePCTable->GetEnergyCutsVect 655 ecut = (*(thePCTable->GetEnergyCutsVector(idxG4ElectronCut)))[matCut->GetIndex()]; 648 limit = 2.*ecut; 656 limit = 2.*ecut; 649 } else { 657 } else { 650 ecut = (*(thePCTable->GetEnergyCutsVect 658 ecut = (*(thePCTable->GetEnergyCutsVector(idxG4PositronCut)))[matCut->GetIndex()]; 651 limit = ecut; 659 limit = ecut; 652 } 660 } 653 G4double min = std::max(limit,fLowEnergyLi 661 G4double min = std::max(limit,fLowEnergyLimit); 654 G4double max = fHighEnergyLimit; 662 G4double max = fHighEnergyLimit; 655 if (min>=max) { 663 if (min>=max) { 656 fSCPCPerMatCuts[imc] = new SCPCorrection 664 fSCPCPerMatCuts[imc] = new SCPCorrection(); 657 fSCPCPerMatCuts[imc]->fIsUse = false; 665 fSCPCPerMatCuts[imc]->fIsUse = false; 658 fSCPCPerMatCuts[imc]->fPrCut = min; 666 fSCPCPerMatCuts[imc]->fPrCut = min; 659 continue; 667 continue; 660 } 668 } 661 G4int numEbins = fNumSPCEbinPerDec*G 669 G4int numEbins = fNumSPCEbinPerDec*G4lrint(std::log10(max/min)); 662 numEbins = std::max(numEbins,3 670 numEbins = std::max(numEbins,3); 663 G4double lmin = G4Log(min); 671 G4double lmin = G4Log(min); 664 G4double ldel = G4Log(max/min)/(num 672 G4double ldel = G4Log(max/min)/(numEbins-1.0); 665 fSCPCPerMatCuts[imc] = new SCPCorrection() 673 fSCPCPerMatCuts[imc] = new SCPCorrection(); 666 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi 674 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbins,1.0); 667 fSCPCPerMatCuts[imc]->fIsUse = true; 675 fSCPCPerMatCuts[imc]->fIsUse = true; 668 fSCPCPerMatCuts[imc]->fPrCut = min; 676 fSCPCPerMatCuts[imc]->fPrCut = min; 669 fSCPCPerMatCuts[imc]->fLEmin = lmin; 677 fSCPCPerMatCuts[imc]->fLEmin = lmin; 670 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; 678 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; 671 for (G4int ie=0; ie<numEbins; ++ie) { 679 for (G4int ie=0; ie<numEbins; ++ie) { 672 G4double ekin = G4Exp(lmin+ie*ldel); 680 G4double ekin = G4Exp(lmin+ie*ldel); 673 G4double scpCorr = 1.0; 681 G4double scpCorr = 1.0; 674 // compute correction factor: I.Kawrakow 682 // compute correction factor: I.Kawrakow NIMB 114(1996)307-326 (Eqs(32-37)) 675 if (ie>0) { 683 if (ie>0) { 676 G4double tau = ekin/CLHEP::electr 684 G4double tau = ekin/CLHEP::electron_mass_c2; 677 G4double tauCut = ecut/CLHEP::electr 685 G4double tauCut = ecut/CLHEP::electron_mass_c2; 678 // Moliere's screening parameter 686 // Moliere's screening parameter 679 G4int matindx = (G4int)matCut->Get << 687 G4int matindx = matCut->GetMaterial()->GetIndex(); 680 G4double A = GetMoliereXc2(mati 688 G4double A = GetMoliereXc2(matindx)/(4.0*tau*(tau+2.)*GetMoliereBc(matindx)); 681 G4double gr = (1.+2.*A)*G4Log(1. 689 G4double gr = (1.+2.*A)*G4Log(1.+1./A)-2.; 682 G4double dum0 = (tau+2.)/(tau+1.); 690 G4double dum0 = (tau+2.)/(tau+1.); 683 G4double dum1 = tau+1.; 691 G4double dum1 = tau+1.; 684 G4double gm = G4Log(0.5*tau/tauC 692 G4double gm = G4Log(0.5*tau/tauCut) + (1.+dum0*dum0)*G4Log(2.*(tau-tauCut+2.)/(tau+4.)) 685 - 0.25*(tau+2.)*( 693 - 0.25*(tau+2.)*( tau+2.+2.*(2.*tau+1.)/(dum1*dum1))* 686 G4Log((tau+4.)*( 694 G4Log((tau+4.)*(tau-tauCut)/tau/(tau-tauCut+2.)) 687 + 0.5*(tau-2*tauCu 695 + 0.5*(tau-2*tauCut)*(tau+2.)*(1./(tau-tauCut)-1./(dum1*dum1)); 688 if (gm<gr) { 696 if (gm<gr) { 689 gm = gm/gr; 697 gm = gm/gr; 690 } else { 698 } else { 691 gm = 1.; 699 gm = 1.; 692 } 700 } 693 G4double z0 = matCut->GetMaterial()-> 701 G4double z0 = matCut->GetMaterial()->GetIonisation()->GetZeffective(); 694 scpCorr = 1.-gm*z0/(z0*(z0+1.)); 702 scpCorr = 1.-gm*z0/(z0*(z0+1.)); 695 } 703 } 696 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo 704 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr; 697 } 705 } 698 } 706 } 699 } 707 } 700 708