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
00040
00041
00042
00043
00044
00045
00046
00047 #include "G4PenelopeBremsstrahlungAngular.hh"
00048
00049 #include "globals.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4PhysicsFreeVector.hh"
00053 #include "G4PhysicsTable.hh"
00054 #include "G4Material.hh"
00055 #include "Randomize.hh"
00056
00057 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular() :
00058 G4VEmAngularDistribution("Penelope"), theEffectiveZSq(0),
00059 theLorentzTables1(0),theLorentzTables2(0)
00060
00061 {
00062 dataRead = false;
00063 verbosityLevel = 0;
00064 }
00065
00066
00067
00068 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular()
00069 {
00070 ClearTables();
00071 }
00072
00073
00074
00075 void G4PenelopeBremsstrahlungAngular::Initialize()
00076 {
00077 ClearTables();
00078 }
00079
00080
00081
00082 void G4PenelopeBremsstrahlungAngular::ClearTables()
00083 {
00084 std::map<G4double,G4PhysicsTable*>::iterator j;
00085
00086 if (theLorentzTables1)
00087 {
00088 for (j=theLorentzTables1->begin(); j != theLorentzTables1->end(); j++)
00089 {
00090 G4PhysicsTable* tab = j->second;
00091 tab->clearAndDestroy();
00092 delete tab;
00093 }
00094 delete theLorentzTables1;
00095 theLorentzTables1 = 0;
00096 }
00097
00098 if (theLorentzTables2)
00099 {
00100 for (j=theLorentzTables2->begin(); j != theLorentzTables2->end(); j++)
00101 {
00102 G4PhysicsTable* tab = j->second;
00103 tab->clearAndDestroy();
00104 delete tab;
00105 }
00106 delete theLorentzTables2;
00107 theLorentzTables2 = 0;
00108 }
00109 if (theEffectiveZSq)
00110 {
00111 delete theEffectiveZSq;
00112 theEffectiveZSq = 0;
00113 }
00114 }
00115
00116
00117
00118 void G4PenelopeBremsstrahlungAngular::ReadDataFile()
00119 {
00120
00121 char* path = getenv("G4LEDATA");
00122 if (!path)
00123 {
00124 G4String excep =
00125 "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
00126 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
00127 "em0006",FatalException,excep);
00128 return;
00129 }
00130 G4String pathString(path);
00131 G4String pathFile = pathString + "/penelope/bremsstrahlung/pdbrang.p08";
00132 std::ifstream file(pathFile);
00133
00134 if (!file.is_open())
00135 {
00136 G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!";
00137 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
00138 "em0003",FatalException,excep);
00139 return;
00140 }
00141 G4int i=0,j=0,k=0;
00142
00143 for (k=0;k<NumberofKPoints;k++)
00144 for (i=0;i<NumberofZPoints;i++)
00145 for (j=0;j<NumberofEPoints;j++)
00146 {
00147 G4double a1,a2;
00148 G4int ik1,iz1,ie1;
00149 G4double zr,er,kr;
00150 file >> iz1 >> ie1 >> ik1 >> zr >> er >> kr >> a1 >> a2;
00151
00152 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 == j))
00153 {
00154 QQ1[i][j][k]=a1;
00155 QQ2[i][j][k]=a2;
00156 }
00157 else
00158 {
00159 G4ExceptionDescription ed;
00160 ed << "Corrupted data file " << pathFile << "?" << G4endl;
00161 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
00162 "em0005",FatalException,ed);
00163 }
00164 }
00165 file.close();
00166 dataRead = true;
00167 }
00168
00169
00170
00171 void G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables(G4double Zmat)
00172 {
00173
00174 if (!dataRead)
00175 {
00176 ReadDataFile();
00177 if (!dataRead)
00178 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
00179 "em2001",FatalException,"Unable to build interpolation table");
00180 }
00181
00182 const G4int reducedEnergyGrid=21;
00183
00184 G4double betas[NumberofEPoints];
00185
00186 G4double Q1[NumberofEPoints][NumberofKPoints];
00187 G4double Q2[NumberofEPoints][NumberofKPoints];
00188
00189 G4double Q1E[NumberofEPoints][reducedEnergyGrid];
00190 G4double Q2E[NumberofEPoints][reducedEnergyGrid];
00191 G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
00192
00193 G4int i=0,j=0,k=0;
00194
00195 for (i=0;i<NumberofEPoints;i++)
00196 {
00197 for (j=0;j<NumberofKPoints;j++)
00198 {
00199 G4PhysicsFreeVector* QQ1vector = new G4PhysicsFreeVector(NumberofZPoints);
00200 QQ1vector->SetSpline(true);
00201 G4PhysicsFreeVector* QQ2vector = new G4PhysicsFreeVector(NumberofZPoints);
00202 QQ2vector->SetSpline(true);
00203
00204
00205 for (k=0;k<NumberofZPoints;k++)
00206 {
00207 QQ1vector->PutValue(k,pZ[k],std::log(QQ1[k][i][j]));
00208 QQ2vector->PutValue(k,pZ[k],QQ2[k][i][j]);
00209 }
00210
00211 Q1[i][j]= std::exp(QQ1vector->Value(Zmat));
00212 Q2[i][j]=QQ2vector->Value(Zmat);
00213 delete QQ1vector;
00214 delete QQ2vector;
00215 }
00216 }
00217 G4double pE[NumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV,
00218 1.0e-01*MeV,5.0e-01*MeV};
00219 G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
00220 G4double ppK[reducedEnergyGrid];
00221
00222 for(i=0;i<reducedEnergyGrid;i++)
00223 ppK[i]=((G4double) i) * 0.05;
00224
00225
00226 for(i=0;i<NumberofEPoints;i++)
00227 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
00228
00229
00230 for (i=0;i<NumberofEPoints;i++)
00231 {
00232 for (j=0;j<NumberofKPoints;j++)
00233 Q1[i][j]=Q1[i][j]/Zmat;
00234 }
00235
00236
00237 for (i=0;i<NumberofEPoints;i++)
00238 {
00239 G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(NumberofKPoints);
00240 G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(NumberofKPoints);
00241
00242 for (j=0;j<NumberofKPoints;j++)
00243 {
00244 Q1vector->PutValue(j,pK[j],std::log(Q1[i][j]));
00245 Q2vector->PutValue(j,pK[j],Q2[i][j]);
00246 }
00247
00248 for (j=0;j<reducedEnergyGrid;j++)
00249 {
00250 Q1E[i][j]=Q1vector->Value(ppK[j]);
00251 Q2E[i][j]=Q2vector->Value(ppK[j]);
00252 }
00253 delete Q1vector;
00254 delete Q2vector;
00255 }
00256
00257
00258
00259 G4PhysicsTable* theTable1 = new G4PhysicsTable();
00260 G4PhysicsTable* theTable2 = new G4PhysicsTable();
00261
00262
00263
00264
00265
00266
00267 for (j=0;j<reducedEnergyGrid;j++)
00268 {
00269 G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(NumberofEPoints);
00270 thevec->SetSpline(true);
00271 theTable1->push_back(thevec);
00272 G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(NumberofEPoints);
00273 thevec2->SetSpline(true);
00274 theTable2->push_back(thevec2);
00275 }
00276
00277 for (j=0;j<reducedEnergyGrid;j++)
00278 {
00279 G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j];
00280 G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j];
00281 for (i=0;i<NumberofEPoints;i++)
00282 {
00283 thevec->PutValue(i,betas[i],Q1E[i][j]);
00284 thevec2->PutValue(i,betas[i],Q2E[i][j]);
00285 }
00286 }
00287
00288 if (theLorentzTables1 && theLorentzTables2)
00289 {
00290 theLorentzTables1->insert(std::make_pair(Zmat,theTable1));
00291 theLorentzTables2->insert(std::make_pair(Zmat,theTable2));
00292 }
00293 else
00294 {
00295 G4ExceptionDescription ed;
00296 ed << "Unable to create tables of Lorentz coefficients for " << G4endl;
00297 ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl;
00298 delete theTable1;
00299 delete theTable2;
00300 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
00301 "em2005",FatalException,ed);
00302 }
00303
00304 return;
00305 }
00306
00307
00308
00309 G4ThreeVector& G4PenelopeBremsstrahlungAngular::SampleDirection(const G4DynamicParticle* dp,
00310 G4double eGamma,
00311 G4int,
00312 const G4Material* material)
00313 {
00314 if (!material)
00315 {
00316 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
00317 "em2040",FatalException,"The pointer to G4Material* is NULL");
00318 return fLocalDirection;
00319 }
00320
00321 G4double Zmat = GetEffectiveZ(material);
00322 if (verbosityLevel > 0)
00323 {
00324 G4cout << "Effective <Z> for material : " << material->GetName() <<
00325 " = " << Zmat << G4endl;
00326 }
00327
00328 G4double ePrimary = dp->GetKineticEnergy();
00329
00330 G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
00331 (ePrimary+electron_mass_c2);
00332 G4double cdt = 0;
00333 G4double sinTheta = 0;
00334 G4double phi = 0;
00335
00336
00337 if (ePrimary > 500*keV)
00338 {
00339 cdt = 2.0*G4UniformRand() - 1.0;
00340 if (G4UniformRand() > 0.75)
00341 {
00342 if (cdt<0)
00343 cdt = -1.0*std::pow(-cdt,1./3.);
00344 else
00345 cdt = std::pow(cdt,1./3.);
00346 }
00347 cdt = (cdt+beta)/(1.0+beta*cdt);
00348
00349 sinTheta = std::sqrt(1. - cdt*cdt);
00350 phi = twopi * G4UniformRand();
00351 fLocalDirection.set(sinTheta* std::cos(phi),
00352 sinTheta* std::sin(phi),
00353 cdt);
00354
00355 fLocalDirection.rotateUz(dp->GetMomentumDirection());
00356
00357 return fLocalDirection;
00358 }
00359
00360
00361 if (!theLorentzTables1)
00362 theLorentzTables1 = new std::map<G4double,G4PhysicsTable*>;
00363 if (!theLorentzTables2)
00364 theLorentzTables2 = new std::map<G4double,G4PhysicsTable*>;
00365
00366
00367 if (!(theLorentzTables1->count(Zmat)))
00368 PrepareInterpolationTables(Zmat);
00369
00370 if (!(theLorentzTables1->count(Zmat)) || !(theLorentzTables2->count(Zmat)))
00371 {
00372 G4ExceptionDescription ed;
00373 ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl;
00374 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
00375 "em2006",FatalException,ed);
00376 }
00377
00378
00379 G4PhysicsTable* theTable1 = theLorentzTables1->find(Zmat)->second;
00380 G4PhysicsTable* theTable2 = theLorentzTables2->find(Zmat)->second;
00381
00382 G4double RK=20.0*eGamma/ePrimary;
00383 G4int ik=std::min((G4int) RK,19);
00384
00385 G4double P10=0,P11=0,P1=0;
00386 G4double P20=0,P21=0,P2=0;
00387
00388
00389 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik];
00390 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1];
00391 P10 = v1->Value(beta);
00392 P11 = v2->Value(beta);
00393 P1=P10+(RK-(G4double) ik)*(P11-P10);
00394
00395
00396 G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik];
00397 G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1];
00398 P20=v3->Value(beta);
00399 P21=v4->Value(beta);
00400 P2=P20+(RK-(G4double) ik)*(P21-P20);
00401
00402
00403 P1=std::min(std::exp(P1)/beta,1.0);
00404 G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
00405
00406 G4double testf=0;
00407
00408 if (G4UniformRand() < P1)
00409 {
00410 do{
00411 cdt = 2.0*G4UniformRand()-1.0;
00412 testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
00413 }while(testf>0);
00414 }
00415 else
00416 {
00417 do{
00418 cdt = 2.0*G4UniformRand()-1.0;
00419 testf=G4UniformRand()-(1.0-cdt*cdt);
00420 }while(testf>0);
00421 }
00422 cdt = (cdt+betap)/(1.0+betap*cdt);
00423
00424
00425 sinTheta = std::sqrt(1. - cdt*cdt);
00426 phi = twopi * G4UniformRand();
00427 fLocalDirection.set(sinTheta* std::cos(phi),
00428 sinTheta* std::sin(phi),
00429 cdt);
00430
00431 fLocalDirection.rotateUz(dp->GetMomentumDirection());
00432
00433 return fLocalDirection;
00434
00435 }
00436
00437
00438
00439 G4double G4PenelopeBremsstrahlungAngular::PolarAngle(const G4double ,
00440 const G4double ,
00441 const G4int )
00442 {
00443 G4cout << "WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" << G4endl;
00444 G4cout << "Please use the alternative interface SampleDirection()" << G4endl;
00445 G4Exception("G4PenelopeBremsstrahlungAngular::PolarAngle()",
00446 "em0005",FatalException,"Unsupported interface");
00447 return 0;
00448 }
00449
00450
00451
00452 G4double G4PenelopeBremsstrahlungAngular::GetEffectiveZ(const G4Material* material)
00453 {
00454 if (!theEffectiveZSq)
00455 theEffectiveZSq = new std::map<const G4Material*,G4double>;
00456
00457
00458 if (theEffectiveZSq->count(material))
00459 return theEffectiveZSq->find(material)->second;
00460
00461
00462
00463 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
00464 G4int nElements = material->GetNumberOfElements();
00465 const G4ElementVector* elementVector = material->GetElementVector();
00466 const G4double* fractionVector = material->GetFractionVector();
00467 for (G4int i=0;i<nElements;i++)
00468 {
00469 G4double fraction = fractionVector[i];
00470 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
00471 StechiometricFactors->push_back(fraction/atomicWeigth);
00472 }
00473
00474 G4double MaxStechiometricFactor = 0.;
00475 for (G4int i=0;i<nElements;i++)
00476 {
00477 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
00478 MaxStechiometricFactor = (*StechiometricFactors)[i];
00479 }
00480
00481 for (G4int i=0;i<nElements;i++)
00482 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
00483
00484 G4double sumz2 = 0;
00485 G4double sums = 0;
00486 for (G4int i=0;i<nElements;i++)
00487 {
00488 G4double Z = (*elementVector)[i]->GetZ();
00489 sumz2 += (*StechiometricFactors)[i]*Z*Z;
00490 sums += (*StechiometricFactors)[i];
00491 }
00492 delete StechiometricFactors;
00493
00494 G4double ZBR = std::sqrt(sumz2/sums);
00495 theEffectiveZSq->insert(std::make_pair(material,ZBR));
00496
00497 return ZBR;
00498 }