00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4PenelopeBremsstrahlungFS.hh"
00042 #include "G4PhysicsFreeVector.hh"
00043 #include "G4PhysicsLogVector.hh"
00044 #include "G4PhysicsTable.hh"
00045 #include "G4Material.hh"
00046 #include "Randomize.hh"
00047
00048
00049
00050 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS() :
00051 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
00052 thePBcut(0)
00053 {
00054 G4double tempvector[nBinsX] =
00055 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
00056 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
00057 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
00058 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
00059
00060 for (size_t ix=0;ix<nBinsX;ix++)
00061 theXGrid[ix] = tempvector[ix];
00062
00063 for (size_t i=0;i<nBinsE;i++)
00064 theEGrid[i] = 0.;
00065
00066 theElementData = new std::map<G4int,G4DataVector*>;
00067 theTempVec = new G4PhysicsFreeVector(nBinsX);
00068 }
00069
00070
00071
00072 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS()
00073 {
00074 ClearTables();
00075
00076 if (theTempVec)
00077 delete theTempVec;
00078
00079
00080 std::map<G4int,G4DataVector*>::iterator i;
00081 if (theElementData)
00082 {
00083 for (i=theElementData->begin(); i != theElementData->end(); i++)
00084 delete i->second;
00085 delete theElementData;
00086 theElementData = 0;
00087 }
00088
00089 }
00090
00091
00092
00093
00094 void G4PenelopeBremsstrahlungFS::ClearTables()
00095 {
00096 std::map< std::pair<const G4Material*,G4double> ,G4PhysicsTable*>::iterator j;
00097
00098 if (theReducedXSTable)
00099 {
00100 for (j=theReducedXSTable->begin(); j != theReducedXSTable->end(); j++)
00101 {
00102 G4PhysicsTable* tab = j->second;
00103 tab->clearAndDestroy();
00104 delete tab;
00105 }
00106 delete theReducedXSTable;
00107 theReducedXSTable = 0;
00108 }
00109
00110 if (theSamplingTable)
00111 {
00112 for (j=theSamplingTable->begin(); j != theSamplingTable->end(); j++)
00113 {
00114 G4PhysicsTable* tab = j->second;
00115 tab->clearAndDestroy();
00116 delete tab;
00117 }
00118 delete theSamplingTable;
00119 theSamplingTable = 0;
00120 }
00121
00122 std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk;
00123 if (thePBcut)
00124 {
00125 for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)
00126 delete kk->second;
00127 delete thePBcut;
00128 thePBcut = 0;
00129 }
00130
00131
00132 if (theEffectiveZSq)
00133 {
00134 delete theEffectiveZSq;
00135 theEffectiveZSq = 0;
00136 }
00137
00138 return;
00139 }
00140
00141
00142
00143 G4double G4PenelopeBremsstrahlungFS::GetEffectiveZSquared(const G4Material* material)
00144 {
00145 if (!theEffectiveZSq)
00146 {
00147 G4ExceptionDescription ed;
00148 ed << "The container for the <Z^2> values is not initialized" << G4endl;
00149 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
00150 "em2007",FatalException,ed);
00151 return 0;
00152 }
00153
00154 if (theEffectiveZSq->count(material))
00155 return theEffectiveZSq->find(material)->second;
00156 else
00157 {
00158 G4ExceptionDescription ed;
00159 ed << "The value of <Z^2> is not properly set for material " <<
00160 material->GetName() << G4endl;
00161
00162 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
00163 "em2008",FatalException,ed);
00164 }
00165 return 0;
00166 }
00167
00168
00169
00170 void G4PenelopeBremsstrahlungFS::BuildScaledXSTable(const G4Material* material,
00171 G4double cut)
00172 {
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
00184 G4int nElements = material->GetNumberOfElements();
00185 const G4ElementVector* elementVector = material->GetElementVector();
00186 const G4double* fractionVector = material->GetFractionVector();
00187 for (G4int i=0;i<nElements;i++)
00188 {
00189 G4double fraction = fractionVector[i];
00190 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
00191 StechiometricFactors->push_back(fraction/atomicWeigth);
00192 }
00193
00194 G4double MaxStechiometricFactor = 0.;
00195 for (G4int i=0;i<nElements;i++)
00196 {
00197 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
00198 MaxStechiometricFactor = (*StechiometricFactors)[i];
00199 }
00200
00201 for (G4int i=0;i<nElements;i++)
00202 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
00203
00204 G4double sumz2 = 0;
00205 G4double sums = 0;
00206 for (G4int i=0;i<nElements;i++)
00207 {
00208 G4double Z = (*elementVector)[i]->GetZ();
00209 sumz2 += (*StechiometricFactors)[i]*Z*Z;
00210 sums += (*StechiometricFactors)[i];
00211 }
00212 G4double ZBR2 = sumz2/sums;
00213
00214 theEffectiveZSq->insert(std::make_pair(material,ZBR2));
00215
00216
00217
00218
00219
00220 G4DataVector* tempData = new G4DataVector(nBinsE);
00221 G4DataVector* tempMatrix = new G4DataVector(nBinsE*nBinsX,0.);
00222
00223 for (G4int iel=0;iel<nElements;iel++)
00224 {
00225 G4double Z = (*elementVector)[iel]->GetZ();
00226 G4int iZ = (G4int) Z;
00227 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
00228
00229
00230
00231 if (!theElementData->count(iZ))
00232 {
00233 ReadDataFile(iZ);
00234 if (!theElementData->count(iZ))
00235 {
00236 G4ExceptionDescription ed;
00237 ed << "Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" << G4endl;
00238 ed << "Unable to retrieve data for element " << iZ << G4endl;
00239 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
00240 "em2009",FatalException,ed);
00241 }
00242 }
00243
00244 G4DataVector* atomData = theElementData->find(iZ)->second;
00245
00246 for (size_t ie=0;ie<nBinsE;ie++)
00247 {
00248 (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX];
00249 for (size_t ix=0;ix<nBinsX;ix++)
00250 (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
00251 }
00252 }
00253
00254
00255
00256
00257
00258 for (size_t ie=0;ie<nBinsE;ie++)
00259 {
00260
00261 G4double* tempData2 = new G4double[nBinsX];
00262 for (size_t ix=0;ix<nBinsX;ix++)
00263 tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
00264 G4double rsum = GetMomentumIntegral(tempData2,1.0,0);
00265 delete[] tempData2;
00266 G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/
00267 (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2));
00268 G4double fnorm = (*tempData)[ie]/(rsum*fact);
00269 G4double TST = 100.*std::fabs(fnorm-1.0);
00270 if (TST > 1.0)
00271 {
00272 G4ExceptionDescription ed;
00273 ed << "G4PenelopeBremsstrahlungFS. Corrupted data files?" << G4endl;
00274 G4cout << "TST= " << TST << "; fnorm = " << fnorm << G4endl;
00275 G4cout << "rsum = " << rsum << G4endl;
00276 G4cout << "fact = " << fact << G4endl;
00277 G4cout << ie << " " << theEGrid[ie]/keV << " " << (*tempData)[ie]/barn << G4endl;
00278 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
00279 "em2010",FatalException,ed);
00280 }
00281 for (size_t ix=0;ix<nBinsX;ix++)
00282 (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
00283 }
00284
00285
00286
00287
00288 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
00289
00290
00291
00292
00293
00294
00295
00296 for (size_t i=0;i<nBinsX;i++)
00297 thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsE+1));
00298
00299 for (size_t ix=0;ix<nBinsX;ix++)
00300 {
00301 G4PhysicsFreeVector* theVec =
00302 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
00303 for (size_t ie=0;ie<nBinsE;ie++)
00304 {
00305 G4double logene = std::log(theEGrid[ie]);
00306 G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
00307 if (aValue < 1e-20*millibarn)
00308 aValue = 1e-20*millibarn;
00309 theVec->PutValue(ie+1,logene,std::log(aValue));
00310 }
00311
00312
00313 G4double derivative = ((*theVec)[2]-(*theVec)[1])/(theVec->Energy(2) - theVec->Energy(1));
00314 G4double log1eV = std::log(1*eV);
00315 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->Energy(1));
00316
00317 theVec->PutValue(0,log1eV,val1eV);
00318 }
00319 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
00320 theReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable));
00321
00322 delete StechiometricFactors;
00323 delete tempData;
00324 delete tempMatrix;
00325
00326 return;
00327 }
00328
00329
00330
00331 void G4PenelopeBremsstrahlungFS::ReadDataFile(G4int Z)
00332 {
00333
00334 char* path = getenv("G4LEDATA");
00335 if (!path)
00336 {
00337 G4String excep = "G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
00338 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
00339 "em0006",FatalException,excep);
00340 return;
00341 }
00342
00343
00344
00345 std::ostringstream ost;
00346 if (Z>9)
00347 ost << path << "/penelope/bremsstrahlung/pdebr" << Z << ".p08";
00348 else
00349 ost << path << "/penelope/bremsstrahlung/pdebr0" << Z << ".p08";
00350 std::ifstream file(ost.str().c_str());
00351 if (!file.is_open())
00352 {
00353 G4String excep = "G4PenelopeBremsstrahlungFS - data file " +
00354 G4String(ost.str()) + " not found!";
00355 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
00356 "em0003",FatalException,excep);
00357 return;
00358 }
00359
00360 G4int readZ =0;
00361 file >> readZ;
00362
00363
00364 if (readZ != Z)
00365 {
00366 G4ExceptionDescription ed;
00367 ed << "Corrupted data file for Z=" << Z << G4endl;
00368 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
00369 "em0005",FatalException,ed);
00370 return;
00371 }
00372
00373 G4DataVector* theMatrix = new G4DataVector(nBinsE*(nBinsX+1),0.);
00374
00375 for (size_t ie=0;ie<nBinsE;ie++)
00376 {
00377 G4double myDouble = 0;
00378 file >> myDouble;
00379 if (!theEGrid[ie])
00380 theEGrid[ie] = myDouble*eV;
00381
00382 for (size_t ix=0;ix<nBinsX;ix++)
00383 {
00384 file >> myDouble;
00385 (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*millibarn;
00386 }
00387 file >> myDouble;
00388 (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*millibarn;
00389 }
00390
00391 if (theElementData)
00392 theElementData->insert(std::make_pair(Z,theMatrix));
00393 else
00394 delete theMatrix;
00395 file.close();
00396 return;
00397 }
00398
00399
00400 G4double G4PenelopeBremsstrahlungFS::GetMomentumIntegral(G4double* y,
00401 G4double xup,G4int momOrder)
00402
00403 {
00404
00405
00406
00407
00408
00409 size_t size = nBinsX;
00410 const G4double eps = 1e-35;
00411
00412
00413 if (momOrder<-1 || size<2 || theXGrid[0]<0)
00414 {
00415 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
00416 "em2011",FatalException,"Invalid call");
00417 }
00418
00419 for (size_t i=1;i<size;i++)
00420 {
00421 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
00422 {
00423 G4ExceptionDescription ed;
00424 ed << "Invalid call for bin " << i << G4endl;
00425 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
00426 "em2012",FatalException,ed);
00427 }
00428 }
00429
00430
00431 G4double result = 0;
00432 if (xup < theXGrid[0])
00433 return result;
00434 bool loopAgain = true;
00435 G4double xt = std::min(xup,theXGrid[size-1]);
00436 G4double xtc = 0;
00437 for (size_t i=0;i<size-1;i++)
00438 {
00439 G4double x1 = std::max(theXGrid[i],eps);
00440 G4double y1 = y[i];
00441 G4double x2 = std::max(theXGrid[i+1],eps);
00442 G4double y2 = y[i+1];
00443 if (xt < x2)
00444 {
00445 xtc = xt;
00446 loopAgain = false;
00447 }
00448 else
00449 xtc = x2;
00450 G4double dx = x2-x1;
00451 G4double dy = y2-y1;
00452 G4double ds = 0;
00453 if (std::fabs(dx)>1e-14*std::fabs(dy))
00454 {
00455 G4double b=dy/dx;
00456 G4double a=y1-b*x1;
00457 if (momOrder == -1)
00458 ds = a*std::log(xtc/x1)+b*(xtc-x1);
00459 else if (momOrder == 0)
00460 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
00461 else
00462 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1))
00463 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2));
00464 }
00465 else
00466 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
00467 result += ds;
00468 if (!loopAgain)
00469 return result;
00470 }
00471 return result;
00472 }
00473
00474
00475
00476 G4PhysicsTable* G4PenelopeBremsstrahlungFS::GetScaledXSTable(const G4Material* mat,
00477 G4double cut)
00478 {
00479
00480 if (!theReducedXSTable)
00481 theReducedXSTable = new std::map< std::pair<const G4Material*,G4double> ,
00482 G4PhysicsTable*>;
00483 if (!theEffectiveZSq)
00484 theEffectiveZSq = new std::map<const G4Material*,G4double>;
00485
00486
00487 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00488 if (!(theReducedXSTable->count(theKey)))
00489 BuildScaledXSTable(mat,cut);
00490
00491 if (!(theReducedXSTable->count(theKey)))
00492 {
00493 G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
00494 "em2013",FatalException,"Unable to retrieve the cross section table");
00495 }
00496
00497 return theReducedXSTable->find(theKey)->second;
00498 }
00499
00500
00501
00502 void G4PenelopeBremsstrahlungFS::InitializeEnergySampling(const G4Material* material,
00503 G4double cut)
00504 {
00505 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
00506
00507 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
00508
00509
00510 G4PhysicsFreeVector* thePBvec = new G4PhysicsFreeVector(nBinsE);
00511
00512
00513 for (size_t i=0;i<nBinsE;i++)
00514 thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsX));
00515
00516
00517
00518
00519 G4PhysicsTable* theTableReduced = GetScaledXSTable(material,cut);
00520
00521 for (size_t ie=0;ie<nBinsE;ie++)
00522 {
00523 G4PhysicsFreeVector* theVec =
00524 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ie]);
00525
00526 G4double value = 0;
00527 theVec->PutValue(0,theXGrid[0],value);
00528 for (size_t ix=1;ix<nBinsX;ix++)
00529 {
00530
00531
00532 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableReduced)[ix-1];
00533 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
00534
00535 G4double x1=std::max(theXGrid[ix-1],1.0e-35);
00536
00537
00538 G4double y1=std::exp((*v1)[ie+1]);
00539 G4double x2=std::max(theXGrid[ix],1.0e-35);
00540 G4double y2=std::exp((*v2)[ie+1]);
00541 G4double B = (y2-y1)/(x2-x1);
00542 G4double A = y1-B*x1;
00543 G4double dS = A*std::log(x2/x1)+B*(x2-x1);
00544 value += dS;
00545 theVec->PutValue(ix,theXGrid[ix],value);
00546 }
00547
00548
00549 G4double xc = cut/theEGrid[ie];
00550
00551 G4double* tempData = new G4double[nBinsX];
00552 for (size_t ix=0;ix<nBinsX;ix++)
00553 {
00554 G4PhysicsFreeVector* vv = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
00555 tempData[ix] = std::exp((*vv)[ie+1]);
00556 }
00557 G4double pbval = (xc<=1) ?
00558 GetMomentumIntegral(tempData,xc,-1) :
00559 GetMomentumIntegral(tempData,1.0,-1);
00560 thePBvec->PutValue(ie,theEGrid[ie],pbval);
00561 delete[] tempData;
00562 }
00563
00564 theSamplingTable->insert(std::make_pair(theKey,thePhysicsTable));
00565 thePBcut->insert(std::make_pair(theKey,thePBvec));
00566 return;
00567 }
00568
00569
00570
00571 G4double G4PenelopeBremsstrahlungFS::SampleGammaEnergy(G4double energy,const G4Material* mat,
00572 G4double cut)
00573 {
00574 if (!theSamplingTable)
00575 theSamplingTable =
00576 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>;
00577 if (!thePBcut)
00578 thePBcut =
00579 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >;
00580
00581 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00582
00583 if (!(theSamplingTable->count(theKey)))
00584 {
00585 InitializeEnergySampling(mat,cut);
00586 if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
00587 {
00588 G4ExceptionDescription ed;
00589 ed << "Unable to create the SamplingTable: " <<
00590 theSamplingTable->count(theKey) << " " <<
00591 thePBcut->count(theKey) << G4endl;
00592 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00593 "em2014",FatalException,ed);
00594 }
00595 }
00596
00597 G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
00598 G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
00599
00600
00601 size_t eBin = 0;
00602 G4bool firstOrLastBin = false;
00603
00604 if (energy < theEGrid[0])
00605 {
00606 eBin = 0;
00607 firstOrLastBin = true;
00608 }
00609 else if (energy > theEGrid[nBinsE-1])
00610 {
00611 eBin = nBinsE-1;
00612 firstOrLastBin = true;
00613 }
00614 else
00615 {
00616 size_t i=0;
00617 size_t j=nBinsE-1;
00618 while ((j-i)>1)
00619 {
00620 size_t k = (i+j)/2;
00621 if (energy > theEGrid[k])
00622 i = k;
00623 else
00624 j = k;
00625 }
00626 eBin = i;
00627 }
00628
00629
00630 G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin];
00631
00632
00633
00634
00635
00636
00637 if (!firstOrLastBin)
00638 {
00639 G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1];
00640 for (size_t iloop=0;iloop<nBinsX;iloop++)
00641 {
00642 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
00643 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
00644 theTempVec->PutValue(iloop,theXGrid[iloop],val);
00645 }
00646 }
00647 else
00648 {
00649 for (size_t iloop=0;iloop<nBinsX;iloop++)
00650 theTempVec->PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);
00651 }
00652
00653
00654 G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
00655
00656 if (!firstOrLastBin)
00657 {
00658 pbcut = (*(thePBcut->find(theKey)->second))[eBin] +
00659 ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
00660 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
00661 }
00662
00663 G4double pCumulative = (*theTempVec)[nBinsX-1];
00664
00665 G4double eGamma = 0;
00666 do
00667 {
00668 G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
00669
00670
00671 size_t ibin = 0;
00672 if (pt < (*theTempVec)[0])
00673 ibin = 0;
00674 else if (pt > (*theTempVec)[nBinsX-1])
00675 {
00676
00677
00678 G4double delta = pt-(*theTempVec)[nBinsX-1];
00679 if (delta < pt*1e-10)
00680 {
00681 ibin = nBinsX-2;
00682 G4ExceptionDescription ed;
00683 ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
00684 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " <<
00685 (pt-(*theTempVec)[nBinsX-1]) << G4endl;
00686 ed << "Possible symptom of problem with numerical precision" << G4endl;
00687 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00688 "em2015",JustWarning,ed);
00689 }
00690 else
00691 {
00692 G4ExceptionDescription ed;
00693 ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
00694 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " <<
00695 nBinsX << G4endl;
00696 ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" <<
00697 G4endl;
00698 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
00699 "em2015",FatalException,ed);
00700 }
00701 }
00702 else
00703 {
00704 size_t i=0;
00705 size_t j=nBinsX-1;
00706 while ((j-i)>1)
00707 {
00708 size_t k = (i+j)/2;
00709 if (pt > (*theTempVec)[k])
00710 i = k;
00711 else
00712 j = k;
00713 }
00714 ibin = i;
00715 }
00716
00717 G4double w1 = theXGrid[ibin];
00718 G4double w2 = theXGrid[ibin+1];
00719
00720 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
00721 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];
00722
00723
00724 G4double pdf1 = std::exp((*v1)[eBin+1]);
00725 G4double pdf2 = std::exp((*v2)[eBin+1]);
00726 G4double deltaW = w2-w1;
00727 G4double dpdfb = pdf2-pdf1;
00728 G4double B = dpdfb/deltaW;
00729 G4double A = pdf1-B*w1;
00730
00731
00732 G4double wbcut = (cut < energy) ? cut/energy : 1.0;
00733 if (firstOrLastBin)
00734 wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
00735
00736 if (w1 < wbcut)
00737 w1 = wbcut;
00738 if (w2 < w1)
00739 {
00740 G4cout << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
00741 G4cout << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
00742 G4cout << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
00743 G4cout << "cut = " << cut/keV << " keV" << G4endl;
00744 return w1*energy;
00745 }
00746
00747 G4double pmax = std::max(A+B*w1,A+B*w2);
00748 G4bool loopAgain = false;
00749 do
00750 {
00751 loopAgain = false;
00752 eGamma = w1* std::pow((w2/w1),G4UniformRand());
00753 if (G4UniformRand()*pmax > (A+B*eGamma))
00754 loopAgain = true;
00755 }while(loopAgain);
00756 eGamma *= energy;
00757 }while(eGamma < cut);
00758
00759 return eGamma;
00760 }
00761
00762
00763
00764