Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class header file 30 // 31 // 32 // File name: G4eDPWAElasticDCS 33 // 34 // Author: Mihaly Novak 35 // 36 // Creation date: 02.07.2020 37 // 38 // Modifications: 39 // 40 // 41 // ------------------------------------------- 42 43 #include "G4eDPWAElasticDCS.hh" 44 #include "G4EmParameters.hh" 45 #include "G4Physics2DVector.hh" 46 47 #include "zlib.h" 48 49 // 50 // Global variables: 51 // 52 G4bool G4eDPWAElasticDCS::gIsGr 53 G4String G4eDPWAElasticDCS::gData 54 // final values of these variables will be set 55 std::size_t G4eDPWAElasticDCS::gNumE 56 std::size_t G4eDPWAElasticDCS::gIndx 57 std::size_t G4eDPWAElasticDCS::gNumT 58 std::size_t G4eDPWAElasticDCS::gNumT 59 G4double G4eDPWAElasticDCS::gLogM 60 G4double G4eDPWAElasticDCS::gInvD 61 // containers for grids: Ekin, mu(t)=0.5[1-cos 62 std::vector<G4double> G4eDPWAElasticDCS::gTheE 63 std::vector<G4double> G4eDPWAElasticDCS::gTheM 64 std::vector<G4double> G4eDPWAElasticDCS::gTheM 65 std::vector<G4double> G4eDPWAElasticDCS::gTheU 66 std::vector<G4double> G4eDPWAElasticDCS::gTheU 67 // abscissas and weights of an 8 point Gauss-L 68 // for numerical integration on [0,1] 69 const G4double G4eDPWAElasticDCS::gXGL[ 70 1.98550718E-02, 1.01666761E-01, 2.37233795E- 71 5.91717321E-01, 7.62766205E-01, 8.98333239E- 72 }; 73 const G4double G4eDPWAElasticDCS::gWGL[ 74 5.06142681E-02, 1.11190517E-01, 1.56853323E- 75 1.81341892E-01, 1.56853323E-01, 1.11190517E- 76 }; 77 78 79 // - iselectron : data for e- (for e+ otherw 80 // - isrestricted : sampling of angular deflec 81 // required (i.e. in case of 82 G4eDPWAElasticDCS::G4eDPWAElasticDCS(G4bool is 83 : fIsRestrictedSamplingRequired(isrestricted), 84 fDCS.resize(gMaxZ+1, nullptr); 85 fDCSLow.resize(gMaxZ+1, nullptr); 86 fSamplingTables.resize(gMaxZ+1, nullptr); 87 } 88 89 90 // DTR 91 G4eDPWAElasticDCS::~G4eDPWAElasticDCS() { 92 for (std::size_t i=0; i<fDCS.size(); ++i) { 93 if (fDCS[i]) delete fDCS[i]; 94 } 95 for (std::size_t i=0; i<fDCSLow.size(); ++i) 96 if (fDCSLow[i]) delete fDCSLow[i]; 97 } 98 for (std::size_t i=0; i<fSamplingTables.size 99 if (fSamplingTables[i]) delete fSamplingTa 100 } 101 // clear scp correction data 102 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 103 if (fSCPCPerMatCuts[imc]) { 104 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 105 delete fSCPCPerMatCuts[imc]; 106 } 107 } 108 fSCPCPerMatCuts.clear(); 109 } 110 111 112 // initialise for a given 'iz' atomic number: 113 // - nothing happens if it has already been i 114 void G4eDPWAElasticDCS::InitialiseForZ(std::si 115 if (!gIsGridLoaded) { 116 LoadGrid(); 117 } 118 LoadDCSForZ((G4int)iz); 119 BuildSmplingTableForZ((G4int)iz); 120 } 121 122 123 // loads the kinetic energy and theta grids fo 124 // should be called only by the master 125 void G4eDPWAElasticDCS::LoadGrid() { 126 G4String fname = FindDirectoryPath() + "grid 127 std::ifstream infile(fname.c_str()); 128 if (!infile.is_open()) { 129 G4String msg = 130 " Problem while trying to read " + 131 " G4LEDATA version should be G4EML 132 G4Exception("G4eDPWAElasticDCS::ReadCompre 133 FatalException,msg.c_str()); 134 return; 135 } 136 // read size 137 infile >> gNumEnergies; 138 infile >> gNumThetas1; 139 infile >> gNumThetas2; 140 // read the grids 141 // - energy in [MeV] 142 G4double dum = 0.0; 143 gTheEnergies.resize(gNumEnergies); 144 for (std::size_t ie=0; ie<gNumEnergies; ++ie 145 infile >> dum; 146 gTheEnergies[ie] = G4Log(dum*CLHEP::MeV); 147 if (gTheEnergies[ie]<G4Log(2.0*CLHEP::keV) 148 } 149 ++gIndxEnergyLim; 150 // store/set usefull logarithms of the kinet 151 gLogMinEkin = gTheEnergies[0]; 152 gInvDelLogEkin = (gNumEnergies-1)/(gTheEnerg 153 // - theta1 in [deg.] (247): we store mu(the 154 gTheMus1.resize(gNumThetas1); 155 gTheU1.resize(gNumThetas1); 156 const double theA = 0.01; 157 for (std::size_t it=0; it<gNumThetas1; ++it) 158 infile >> dum; 159 gTheMus1[it] = 0.5*(1.0-std::cos(dum*CLHEP 160 gTheU1[it] = (theA+1.0)*gTheMus1[it]/(th 161 } 162 // - theta2 in [deg.] (128): we store mu(the 163 gTheMus2.resize(gNumThetas2); 164 gTheU2.resize(gNumThetas2); 165 for (std::size_t it=0; it<gNumThetas2; ++it) 166 infile >> dum; 167 gTheMus2[it] = 0.5*(1.0-std::cos(dum*CLHEP 168 gTheU2[it] = (theA+1.0)*gTheMus2[it]/(th 169 170 } 171 infile.close(); 172 gIsGridLoaded = true; 173 } 174 175 176 // load DCS data for a given Z 177 void G4eDPWAElasticDCS::LoadDCSForZ(G4int iz) 178 // Check if it has already been done: 179 if (fDCS[iz]) return; 180 // Do it otherwise 181 if (fIsElectron) { 182 // e- 183 // load the high energy part firt: 184 // - with gNumThetas2 theta and gNumEnergi 185 const std::size_t hNumEnergries = gNumEner 186 auto v2DHigh = new G4Physics2DVector(gNumT 187 v2DHigh->SetBicubicInterpolation(true); 188 for (std::size_t it=0; it<gNumThetas2; ++i 189 v2DHigh->PutX(it, gTheMus2[it]); 190 } 191 for (std::size_t ie=0; ie<hNumEnergries; + 192 v2DHigh->PutY(ie, gTheEnergies[gIndxEner 193 } 194 std::ostringstream ossh; 195 ossh << FindDirectoryPath() << "dcss/el/dc 196 std::istringstream finh(std::ios::in); 197 ReadCompressedFile(ossh.str(), finh); 198 G4double dum = 0.0; 199 for (std::size_t it=0; it<gNumThetas2; ++i 200 finh >> dum; 201 for (std::size_t ie=0; ie<hNumEnergries; 202 finh >> dum; 203 v2DHigh->PutValue(it, ie, G4Log(dum*CL 204 } 205 } 206 // load the low energy part: 207 // - with gNumThetas1 theta and gIndxEnerg 208 // for including the firts DCS from the 209 // able to perform interpolation between 210 auto v2DLow = new G4Physics2DVector(gNumTh 211 v2DLow->SetBicubicInterpolation(true); 212 for (std::size_t it=0; it<gNumThetas1; ++i 213 v2DLow->PutX(it, gTheMus1[it]); 214 } 215 for (std::size_t ie=0; ie<gIndxEnergyLim+1 216 v2DLow->PutY(ie, gTheEnergies[ie]); 217 } 218 std::ostringstream ossl; 219 ossl << FindDirectoryPath() << "dcss/el/dc 220 std::istringstream finl(std::ios::in); 221 ReadCompressedFile(ossl.str(), finl); 222 for (std::size_t it=0; it<gNumThetas1; ++i 223 finl >> dum; 224 for (std::size_t ie=0; ie<gIndxEnergyLim 225 finl >> dum; 226 v2DLow->PutValue(it, ie, G4Log(dum*CLH 227 } 228 } 229 // add the +1 part: interpolate the firts 230 std::size_t ix = 0; 231 std::size_t iy = 0; 232 for (std::size_t it=0; it<gNumThetas1; ++i 233 const G4double val = v2DHigh->Value(gThe 234 v2DLow->PutValue(it, gIndxEnergyLim, val 235 } 236 // store 237 fDCSLow[iz] = v2DLow; 238 fDCS[iz] = v2DHigh; 239 } else { 240 // e+ 241 auto v2D= new G4Physics2DVector(gNumThetas 242 v2D->SetBicubicInterpolation(true); 243 for (std::size_t it=0; it<gNumThetas2; ++i 244 v2D->PutX(it, gTheMus2[it]); 245 } 246 for (std::size_t ie=0; ie<gNumEnergies; ++ 247 v2D->PutY(ie, gTheEnergies[ie]); 248 } 249 std::ostringstream oss; 250 oss << FindDirectoryPath() << "dcss/pos/dc 251 std::istringstream fin(std::ios::in); 252 ReadCompressedFile(oss.str(), fin); 253 G4double dum = 0.0; 254 for (std::size_t it=0; it<gNumThetas2; ++i 255 fin >> dum; 256 for (std::size_t ie=0; ie<gNumEnergies; 257 fin >> dum; 258 v2D->PutValue(it, ie, G4Log(dum*CLHEP: 259 } 260 } 261 fDCS[iz]= v2D; 262 } 263 } 264 265 266 267 // Computes the elastic, first and second cros 268 // energy and target atom. 269 // Cross sections are zero ff ekin is below/ab 270 void G4eDPWAElasticDCS::ComputeCSPerAtom(G4int 271 G4dou 272 G4dou 273 // init all cross section values to zero; 274 elcs = 0.0; 275 tr1cs = 0.0; 276 tr2cs = 0.0; 277 // make sure that mu(theta) = 0.5[1-cos(thet 278 mumin = std::max(0.0, std::min(1.0, mumin)); 279 mumax = std::max(0.0, std::min(1.0, mumax)); 280 if (mumin>=mumax) return; 281 // make sure that kin. energy is within the 282 const G4double lekin = std::max(gTheEnergies 283 // if the lower, denser in theta, DCS set sh 284 const G4bool isLowerGrid = (fIsElectron && l 285 const std::vector<G4double>& theMuVector = ( 286 const G4Physics2DVector* the2DDCS = ( 287 // find lower/upper mu bin of integration: 288 // 0.0 <= mumin < 1.0 for sure here 289 const std::size_t iMuStart = (mumin == 0.0) 290 // 0.0 < mumax <= 1.0 for sure here 291 const std::size_t iMuEnd = (mumax == 1.0) 292 // perform numerical integration of the DCS 293 // interval (where mu(theta) = 0.5[1-cos(the 294 std::size_t ix = 0; 295 std::size_t iy = 0; 296 for (std::size_t imu=iMuStart; imu<=iMuEnd; 297 G4double elcsPar = 0.0; 298 G4double tr1csPar = 0.0; 299 G4double tr2csPar = 0.0; 300 const G4double low = (imu==iMuStart) ? mum 301 const G4double del = (imu==iMuEnd) ? mum 302 ix = imu; 303 for (std::size_t igl=0; igl<8; ++igl) { 304 const double mu = low + del*gXGL[igl]; 305 const double dcs = G4Exp(the2DDCS->Value 306 elcsPar += gWGL[igl]*dcs; / 307 tr1csPar += gWGL[igl]*dcs*mu; / 308 tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu); / 309 } 310 elcs += del*elcsPar; 311 tr1cs += del*tr1csPar; 312 tr2cs += del*tr2csPar; 313 } 314 elcs *= 2.0*CLHEP::twopi; 315 tr1cs *= 4.0*CLHEP::twopi; 316 tr2cs *= 12.0*CLHEP::twopi; 317 } 318 319 320 // data structure to store one sampling table: 321 // NOTE: when Alias is used, sampling on a res 322 // However, Alias makes possible faster 323 // of single scattering model while it's 324 // when restricted interval sampling is 325 // the fIsRestrictedSamplingRequired fla 326 struct OneSamplingTable { 327 OneSamplingTable () = default; 328 void SetSize(std::size_t nx, G4bool useAlias 329 fN = nx; 330 // Alias 331 if (useAlias) { 332 fW.resize(nx); 333 fI.resize(nx); 334 } 335 // Ratin 336 fCum.resize(nx); 337 fA.resize(nx); 338 fB.resize(nx); 339 } 340 341 // members 342 std::size_t fN; // # da 343 G4double fScreenParA; // the 344 std::vector<G4double> fW; 345 std::vector<G4double> fCum; 346 std::vector<G4double> fA; 347 std::vector<G4double> fB; 348 std::vector<G4int> fI; 349 }; 350 351 352 // loads sampling table for the given Z over t 353 void G4eDPWAElasticDCS::BuildSmplingTableForZ( 354 // Check if it has already been done: 355 if (fSamplingTables[iz]) return; 356 // Do it otherwise: 357 // allocate space 358 auto sTables = new std::vector<OneSamplingTa 359 // read compressed sampling table data 360 std::ostringstream oss; 361 const G4String fname = fIsElectron ? "stable 362 oss << FindDirectoryPath() << fname << "stab 363 std::istringstream fin(std::ios::in); 364 ReadCompressedFile(oss.str(), fin); 365 std::size_t ndata = 0; 366 for (std::size_t ie=0; ie<gNumEnergies; ++ie 367 OneSamplingTable& aTable = (*sTables)[ie]; 368 // #data in this table 369 fin >> ndata; 370 aTable.SetSize(ndata, !fIsRestrictedSampli 371 // the A screening parameter value used fo 372 fin >> aTable.fScreenParA; 373 // load data: Alias(W,I) + RatIn(Cum, A, B 374 if (!fIsRestrictedSamplingRequired) { 375 for (std::size_t id=0; id<ndata; ++id) { 376 fin >> aTable.fW[id]; 377 } 378 for (std::size_t id=0; id<ndata; ++id) { 379 fin >> aTable.fI[id]; 380 } 381 } 382 for (std::size_t id=0; id<ndata; ++id) { 383 fin >> aTable.fCum[id]; 384 } 385 for (std::size_t id=0; id<ndata; ++id) { 386 fin >> aTable.fA[id]; 387 } 388 for (std::size_t id=0; id<ndata; ++id) { 389 fin >> aTable.fB[id]; 390 } 391 } 392 fSamplingTables[iz] = sTables; 393 } 394 395 396 397 398 399 // samples cos(theta) i.e. cosine of the polar 400 // interaction (Coulomb scattering) of the pro 401 // fIsElectron) with kinetic energy of exp('le 402 // muber of 'iz'. See the 'SampleCosineThetaRe 403 // a restricted inteval. 404 G4double 405 G4eDPWAElasticDCS::SampleCosineTheta(std::size 406 G4double 407 lekin = std::max(gTheEnergies[0], std::min(g 408 // determine the discrete ekin sampling tabl 409 // - statistical interpolation (i.e. linear 410 const G4double rem = (lekin-gLogMinEkin 411 const auto k = (std::size_t)rem; 412 const std::size_t iekin = (r1 < rem-k) ? k+1 413 // sample the mu(t)=0.5(1-cos(t)) 414 const double mu = SampleMu(iz, iekin, r2, 415 return std::max(-1.0, std::min(1.0, 1.0-2.0* 416 } 417 418 419 // samples cos(theta) i.e. cosine of the polar 420 // interaction (Coulomb scattering) of the pro 421 // fIsElectron) with kinetic energy of exp('le 422 // muber of 'iz'. 423 // The cosine theta will be in the [costMin, c 424 // corresponds to a maximum allowed polar scat 425 // costMin corresponds to minimum allowed pola 426 // See the 'SampleCosineTheta' for obtain samp 427 G4double 428 G4eDPWAElasticDCS::SampleCosineThetaRestricted 429 430 431 // costMin corresponds to mu-max while costM 432 lekin = std::max(gTheEnergies[0], std::min(g 433 // determine the discrete ekin sampling tabl 434 // - statistical interpolation (i.e. linear 435 const G4double rem = (lekin-gLogMinEkin 436 const auto k = (size_t)rem; 437 const std::size_t iekin = (r1 < rem-k) ? k : 438 // sample the mu(t)=0.5(1-cos(t)) 439 const G4double mu = SampleMu(iz, iekin, r2, 440 return std::max(-1.0, std::min(1.0, 1.0-2.0* 441 } 442 443 444 445 G4double 446 G4eDPWAElasticDCS::SampleMu(std::size_t izet, 447 OneSamplingTable& rtn = (*fSamplingTables[iz 448 // get the lower index of the bin by using t 449 const G4double rest = r1 * (rtn.fN - 1); 450 auto indxl = (std::size_t)rest; 451 const G4double dum0 = rest - indxl; 452 if (rtn.fW[indxl] < dum0) indxl = rtn.fI[ind 453 // sample value within the selected bin by u 454 const G4double delta = rtn.fCum[indxl + 1] - 455 const G4double aval = r2 * delta; 456 457 const G4double dum1 = (1.0 + rtn.fA[indxl] 458 const G4double dum2 = delta * delta + rtn.f 459 const std::vector<G4double>& theUVect = (fIs 460 const G4double u = theUVect[indxl] + dum 461 // transform back u to mu 462 return rtn.fScreenParA*u/(rtn.fScreenParA+1. 463 } 464 465 466 467 G4double 468 G4eDPWAElasticDCS::FindCumValue(G4double u, co 469 const std::vec 470 const std::size_t iLow = std::distance( uvec 471 const G4double tau = (u-uvect[iLow])/(uv 472 const G4double dum0 = (1.0+stable.fA[iLow 473 const G4double dum1 = 2.0*stable.fB[iLow] 474 const G4double dum2 = 1.0 - std::sqrt(std 475 return std::min(stable.fCum[iLow+1], std::ma 476 } 477 478 // muMin and muMax : no checks on these 479 G4double G4eDPWAElasticDCS::SampleMu(std::size 480 G4double 481 const OneSamplingTable& rtn = (*fSamplingTab 482 const G4double theA = rtn.fScreenParA; 483 // 484 const std::vector<G4double>& theUVect = (fIs 485 const G4double xiMin = (muMin > 0.0) ? Fin 486 const G4double xiMax = (muMax < 1.0) ? Fin 487 // 488 const G4double xi = xiMin+r1*(xiMax-xiM 489 const std::size_t iLow = std::distance( rtn. 490 const G4double delta = rtn.fCum[iLow + 1] 491 const G4double aval = xi - rtn.fCum[iLow] 492 493 const G4double dum1 = (1.0 + rtn.fA[iLow] 494 const G4double dum2 = delta * delta + rtn 495 const G4double u = theUVect[iLow] + du 496 return theA*u/(theA+1.0-u); 497 } 498 499 500 501 // set the DCS data directory path 502 const G4String& G4eDPWAElasticDCS::FindDirecto 503 // check environment variable 504 if (gDataDirectory.empty()) { 505 std::ostringstream ost; 506 ost << G4EmParameters::Instance()->GetDirL 507 gDataDirectory = ost.str(); 508 } 509 return gDataDirectory; 510 } 511 512 513 514 // uncompress one data file into the input str 515 void 516 G4eDPWAElasticDCS::ReadCompressedFile(G4String 517 G4String *dataString = nullptr; 518 G4String compfilename(fname+".z"); 519 // create input stream with binary mode oper 520 std::ifstream in(compfilename, std::ios::bin 521 if (in.good()) { 522 // get current position in the stream (wa 523 std::size_t fileSize = in.tellg(); 524 // set current position being the beginni 525 in.seekg(0,std::ios::beg); 526 // create (zlib) byte buffer for the data 527 Bytef *compdata = new Bytef[fileSize]; 528 while(in) { 529 in.read((char*)compdata, fileSize); 530 } 531 // create (zlib) byte buffer for the unco 532 uLongf complen = (uLongf)(fileSize*4); 533 Bytef *uncompdata = new Bytef[complen]; 534 while (Z_OK!=uncompress(uncompdata, &comp 535 // increase uncompressed byte buffer 536 delete[] uncompdata; 537 complen *= 2; 538 uncompdata = new Bytef[complen]; 539 } 540 // delete the compressed data buffer 541 delete [] compdata; 542 // create a string from the uncompressed 543 dataString = new G4String((char*)uncompda 544 // delete the uncompressed data buffer 545 delete [] uncompdata; 546 } else { 547 G4String msg = 548 " Problem while trying to read " + 549 " G4LEDATA version should be G4EML 550 G4Exception("G4eDPWAElasticDCS::ReadCompre 551 FatalException,msg.c_str()); 552 return; 553 } 554 // create the input string stream from the d 555 if (dataString) { 556 iss.str(*dataString); 557 in.close(); 558 delete dataString; 559 } 560 } 561 562 563 564 565 G4double 566 G4eDPWAElasticDCS::ComputeScatteringPowerCorre 567 568 const G4int imc = matcut->GetIndex(); 569 G4double corFactor = 1.0; 570 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin< 571 return corFactor; 572 } 573 // get the scattering power correction facto 574 const G4double lekin = G4Log(ekin); 575 G4double remaining = (lekin-fSCPCPerMatCut 576 std::size_t lindx = (G4int)remaining; 577 remaining -= lindx; 578 std::size_t imax = fSCPCPerMatCuts[imc]- 579 if (lindx>=imax) { 580 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i 581 } else { 582 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l 583 } 584 return corFactor; 585 } 586 587 588 void G4eDPWAElasticDCS::InitSCPCorrection(G4do 589 G4do 590 // get the material-cuts table 591 G4ProductionCutsTable *thePCTable = G4Produc 592 std::size_t numMatCuts = thePCTab 593 // clear container if any 594 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 595 if (fSCPCPerMatCuts[imc]) { 596 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 597 delete fSCPCPerMatCuts[imc]; 598 fSCPCPerMatCuts[imc] = nullptr; 599 } 600 } 601 // 602 // set size of the container and create the 603 fSCPCPerMatCuts.resize(numMatCuts,nullptr); 604 // loop over the material-cuts and create sc 605 for (G4int imc=0; imc<(G4int)numMatCuts; ++i 606 const G4MaterialCutsCouple *matCut = theP 607 const G4Material* mat = matCut->GetMateri 608 // get e- production cut in the current ma 609 const G4double ecut = (*(thePCTable->Ge 610 const G4double limit = fIsElectron ? 2.0 611 const G4double min = std::max(limit,lo 612 const G4double max = highEnergyLimit; 613 if (min>=max) { 614 fSCPCPerMatCuts[imc] = new SCPCorrection 615 fSCPCPerMatCuts[imc]->fIsUse = false; 616 fSCPCPerMatCuts[imc]->fPrCut = min; 617 continue; 618 } 619 G4int numEbins = fNumSPCEbinPerDec 620 numEbins = std::max(numEbins 621 const G4double lmin = G4Log(min); 622 const G4double ldel = G4Log(max/min)/(n 623 fSCPCPerMatCuts[imc] = new SCPCorrection 624 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi 625 fSCPCPerMatCuts[imc]->fIsUse = true; 626 fSCPCPerMatCuts[imc]->fPrCut = min; 627 fSCPCPerMatCuts[imc]->fLEmin = lmin; 628 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; 629 // compute Moliere material dependet param 630 G4double moliereBc = 0.0; 631 G4double moliereXc2 = 0.0; 632 ComputeMParams(mat, moliereBc, moliereXc2) 633 // compute scattering power correction ove 634 for (G4int ie=0; ie<numEbins; ++ie) { 635 const G4double ekin = G4Exp(lmin+ie*ldel 636 G4double scpCorr = 1.0; 637 // compute correction factor: I.Kawrakow 638 if (ie>0) { 639 const G4double tau = ekin/CLHEP:: 640 const G4double tauCut = ecut/CLHEP:: 641 // Moliere's screening parameter 642 const G4double A = moliereXc2/( 643 const G4double gr = (1.+2.*A)*G4 644 const G4double dum0 = (tau+2.)/(ta 645 const G4double dum1 = tau+1.; 646 G4double gm = G4Log(0.5*tau/tauCut) 647 - 0.25*(tau+2.)*( tau+ 648 G4Log((tau+4.)*(tau- 649 + 0.5*(tau-2*tauCut)*( 650 if (gm<gr) { 651 gm = gm/gr; 652 } else { 653 gm = 1.; 654 } 655 const G4double z0 = matCut->GetMateri 656 scpCorr = 1.-gm*z0/(z0*(z0+1.)); 657 } 658 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo 659 } 660 } 661 } 662 663 664 // compute material dependent Moliere MSC para 665 void G4eDPWAElasticDCS::ComputeMParams(const G 666 G4doubl 667 const G4double const1 = 7821.6; // [ 668 const G4double const2 = 0.1569; // [ 669 const G4double finstrc2 = 5.325135453E-5; / 670 // G4double xi = 1.0; 671 const G4ElementVector* theElemVect = ma 672 const std::size_t numelems = ma 673 // 674 const G4double* theNbAtomsPerVolVect = ma 675 G4double theTotNbAtomsPerVol = ma 676 // 677 G4double zs = 0.0; 678 G4double zx = 0.0; 679 G4double ze = 0.0; 680 G4double sa = 0.0; 681 // 682 for(std::size_t ielem = 0; ielem < numelems 683 const G4double zet = (*theElemVect)[iele 684 const G4double iwa = (*theElemVect)[iele 685 const G4double ipz = theNbAtomsPerVolVec 686 const G4double dum = ipz*zet*(zet+1.0); 687 zs += dum; 688 ze += dum*(-2.0/3.0)*G4Log(zet) 689 zx += dum*G4Log(1.0+3.34*finstr 690 sa += ipz*iwa; 691 } 692 const G4double density = mat->GetDensity()* 693 // 694 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa; 695 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx) 696 697 theBc = const1*density*z0*G4Exp(z1); //[1 698 theXc2 = const2*density*z0; // [MeV2/cm] 699 // change to Geant4 internal units of 1/len 700 theBc *= 1.0/CLHEP::cm; 701 theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm; 702 } 703 704