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 #include "G4PenelopeRayleighModel.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4PenelopeSamplingData.hh"
00040 #include "G4ParticleDefinition.hh"
00041 #include "G4MaterialCutsCouple.hh"
00042 #include "G4ProductionCutsTable.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4PhysicsTable.hh"
00045 #include "G4ElementTable.hh"
00046 #include "G4Element.hh"
00047 #include "G4PhysicsFreeVector.hh"
00048
00049
00050
00051
00052
00053 G4PenelopeRayleighModel::G4PenelopeRayleighModel(const G4ParticleDefinition*,
00054 const G4String& nam)
00055 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),logAtomicCrossSection(0),
00056 atomicFormFactor(0),logFormFactorTable(0),pMaxTable(0),samplingTable(0)
00057 {
00058 fIntrinsicLowEnergyLimit = 100.0*eV;
00059 fIntrinsicHighEnergyLimit = 100.0*GeV;
00060
00061 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00062
00063 verboseLevel= 0;
00064
00065
00066
00067
00068
00069
00070
00071
00072 G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.);
00073 G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit);
00074
00075 G4double logtransitionenergy = std::log(160*keV);
00076 G4double logfactor1 = std::log(10.)/250.;
00077 G4double logfactor2 = logfactor1*10;
00078 logEnergyGridPMax.push_back(logenergy);
00079 do{
00080 if (logenergy < logtransitionenergy)
00081 logenergy += logfactor1;
00082 else
00083 logenergy += logfactor2;
00084 logEnergyGridPMax.push_back(logenergy);
00085 }while (logenergy < logmaxenergy);
00086 }
00087
00088
00089
00090 G4PenelopeRayleighModel::~G4PenelopeRayleighModel()
00091 {
00092 std::map <G4int,G4PhysicsFreeVector*>::iterator i;
00093 if (logAtomicCrossSection)
00094 {
00095 for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
00096 if (i->second) delete i->second;
00097 delete logAtomicCrossSection;
00098 }
00099
00100 if (atomicFormFactor)
00101 {
00102 for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++)
00103 if (i->second) delete i->second;
00104 delete atomicFormFactor;
00105 }
00106
00107 ClearTables();
00108 }
00109
00110
00111 void G4PenelopeRayleighModel::ClearTables()
00112 {
00113 std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i;
00114
00115 if (logFormFactorTable)
00116 {
00117 for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++)
00118 if (i->second) delete i->second;
00119 delete logFormFactorTable;
00120 logFormFactorTable = 0;
00121 }
00122
00123 if (pMaxTable)
00124 {
00125 for (i=pMaxTable->begin(); i != pMaxTable->end(); i++)
00126 if (i->second) delete i->second;
00127 delete pMaxTable;
00128 pMaxTable = 0;
00129 }
00130
00131 std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii;
00132 if (samplingTable)
00133 {
00134 for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++)
00135 if (ii->second) delete ii->second;
00136 delete samplingTable;
00137 samplingTable = 0;
00138 }
00139
00140 return;
00141 }
00142
00143
00144
00145 void G4PenelopeRayleighModel::Initialise(const G4ParticleDefinition* ,
00146 const G4DataVector& )
00147 {
00148 if (verboseLevel > 3)
00149 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
00150
00151
00152 ClearTables();
00153
00154
00155
00156
00157
00158 if (!logAtomicCrossSection)
00159 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
00160 if (!atomicFormFactor)
00161 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
00162
00163 if (!logFormFactorTable)
00164 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
00165 if (!pMaxTable)
00166 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
00167 if (!samplingTable)
00168 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
00169
00170
00171 if (verboseLevel > 0) {
00172 G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl
00173 << "Energy range: "
00174 << LowEnergyLimit() / keV << " keV - "
00175 << HighEnergyLimit() / GeV << " GeV"
00176 << G4endl;
00177 }
00178
00179 if(isInitialised) return;
00180 fParticleChange = GetParticleChangeForGamma();
00181 isInitialised = true;
00182 }
00183
00184
00185
00186 G4double G4PenelopeRayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00187 G4double energy,
00188 G4double Z,
00189 G4double,
00190 G4double,
00191 G4double)
00192 {
00193
00194
00195
00196
00197 if (verboseLevel > 3)
00198 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
00199
00200 G4int iZ = (G4int) Z;
00201
00202
00203 if (!logAtomicCrossSection->count(iZ))
00204 ReadDataFile(iZ);
00205
00206 if (!logAtomicCrossSection->count(iZ))
00207 {
00208 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
00209 "em2040",FatalException,"Unable to load the cross section table");
00210 }
00211
00212 G4double cross = 0;
00213
00214 G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
00215 if (!atom)
00216 {
00217 G4ExceptionDescription ed;
00218 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
00219 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
00220 "em2041",FatalException,ed);
00221 return 0;
00222 }
00223 G4double logene = std::log(energy);
00224 G4double logXS = atom->Value(logene);
00225 cross = std::exp(logXS);
00226
00227 if (verboseLevel > 2)
00228 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
00229 " = " << cross/barn << " barn" << G4endl;
00230 return cross;
00231 }
00232
00233
00234
00235 void G4PenelopeRayleighModel::BuildFormFactorTable(const G4Material* material)
00236 {
00237
00238
00239
00240
00241
00242 G4int nElements = material->GetNumberOfElements();
00243 const G4ElementVector* elementVector = material->GetElementVector();
00244 const G4double* fractionVector = material->GetFractionVector();
00245
00246 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
00247 for (G4int i=0;i<nElements;i++)
00248 {
00249 G4double fraction = fractionVector[i];
00250 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
00251 StechiometricFactors->push_back(fraction/atomicWeigth);
00252 }
00253
00254 G4double MaxStechiometricFactor = 0.;
00255 for (G4int i=0;i<nElements;i++)
00256 {
00257 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
00258 MaxStechiometricFactor = (*StechiometricFactors)[i];
00259 }
00260 if (MaxStechiometricFactor<1e-16)
00261 {
00262 G4ExceptionDescription ed;
00263 ed << "Inconsistent data of atomic composition for " <<
00264 material->GetName() << G4endl;
00265 G4Exception("G4PenelopeRayleighModel::BuildFormFactorTable()",
00266 "em2042",FatalException,ed);
00267 }
00268
00269 for (G4int i=0;i<nElements;i++)
00270 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
00271
00272
00273 G4double atomsPerMolecule = 0;
00274 for (G4int i=0;i<nElements;i++)
00275 atomsPerMolecule += (*StechiometricFactors)[i];
00276
00277
00278
00279
00280 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size());
00281 theFFVec->SetSpline(true);
00282
00283 for (size_t k=0;k<logQSquareGrid.size();k++)
00284 {
00285 G4double ff2 = 0;
00286 for (G4int i=0;i<nElements;i++)
00287 {
00288 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
00289 G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
00290 G4double f = (*theAtomVec)[k];
00291 ff2 += f*f*(*StechiometricFactors)[i];
00292 }
00293 if (ff2)
00294 theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2));
00295 }
00296 logFormFactorTable->insert(std::make_pair(material,theFFVec));
00297
00298 delete StechiometricFactors;
00299
00300 return;
00301 }
00302
00303
00304
00305
00306 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00307 const G4MaterialCutsCouple* couple,
00308 const G4DynamicParticle* aDynamicGamma,
00309 G4double,
00310 G4double)
00311 {
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 if (verboseLevel > 3)
00326 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
00327
00328 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
00329
00330 if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
00331 {
00332 fParticleChange->ProposeTrackStatus(fStopAndKill);
00333 fParticleChange->SetProposedKineticEnergy(0.);
00334 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
00335 return ;
00336 }
00337
00338 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
00339
00340 const G4Material* theMat = couple->GetMaterial();
00341
00342
00343 if (!pMaxTable || !samplingTable)
00344 {
00345 G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
00346 "em2043",FatalException,"Invalid model initialization");
00347 return;
00348 }
00349
00350
00351 if (!(samplingTable->count(theMat)))
00352 InitializeSamplingAlgorithm(theMat);
00353 G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
00354
00355
00356 if (!pMaxTable->count(theMat))
00357 GetPMaxTable(theMat);
00358 G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
00359
00360 G4double cosTheta = 1.0;
00361
00362
00363 G4double qmax = 2.0*photonEnergy0/electron_mass_c2;
00364
00365 if (qmax < 1e-10)
00366 {
00367 G4bool loopAgain=false;
00368 do
00369 {
00370 loopAgain = false;
00371 cosTheta = 1.0-2.0*G4UniformRand();
00372 G4double G = 0.5*(1+cosTheta*cosTheta);
00373 if (G4UniformRand()>G)
00374 loopAgain = true;
00375 }while(loopAgain);
00376 }
00377 else
00378 {
00379 size_t nData = theDataTable->GetNumberOfStoredPoints();
00380 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
00381 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
00382
00383 G4bool loopAgain = false;
00384 G4double MaxPValue = thePMax->Value(photonEnergy0);
00385 G4double xx=0;
00386
00387
00388
00389
00390 do{
00391 loopAgain = false;
00392 G4double RandomMax = G4UniformRand()*MaxPValue;
00393 xx = theDataTable->SampleValue(RandomMax);
00394
00395
00396 if (xx > q2max)
00397 loopAgain = true;
00398 cosTheta = 1.0-2.0*xx/q2max;
00399 G4double G = 0.5*(1+cosTheta*cosTheta);
00400 if (G4UniformRand()>G)
00401 loopAgain = true;
00402 }while(loopAgain);
00403 }
00404
00405 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
00406
00407
00408 G4double phi = twopi * G4UniformRand() ;
00409 G4double dirX = sinTheta*std::cos(phi);
00410 G4double dirY = sinTheta*std::sin(phi);
00411 G4double dirZ = cosTheta;
00412
00413
00414 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
00415
00416 photonDirection1.rotateUz(photonDirection0);
00417 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
00418 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
00419
00420 return;
00421 }
00422
00423
00424
00425
00426 void G4PenelopeRayleighModel::ReadDataFile(const G4int Z)
00427 {
00428 if (verboseLevel > 2)
00429 {
00430 G4cout << "G4PenelopeRayleighModel::ReadDataFile()" << G4endl;
00431 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
00432 }
00433
00434 char* path = getenv("G4LEDATA");
00435 if (!path)
00436 {
00437 G4String excep = "G4LEDATA environment variable not set!";
00438 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00439 "em0006",FatalException,excep);
00440 return;
00441 }
00442
00443
00444
00445
00446 std::ostringstream ost;
00447 if (Z>9)
00448 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
00449 else
00450 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
00451 std::ifstream file(ost.str().c_str());
00452 if (!file.is_open())
00453 {
00454 G4String excep = "Data file " + G4String(ost.str()) + " not found!";
00455 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00456 "em0003",FatalException,excep);
00457 }
00458 G4int readZ =0;
00459 size_t nPoints= 0;
00460 file >> readZ >> nPoints;
00461
00462 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
00463 {
00464 G4ExceptionDescription ed;
00465 ed << "Corrupted data file for Z=" << Z << G4endl;
00466 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00467 "em0005",FatalException,ed);
00468 return;
00469 }
00470 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
00471 G4double ene=0,f1=0,f2=0,xs=0;
00472 for (size_t i=0;i<nPoints;i++)
00473 {
00474 file >> ene >> f1 >> f2 >> xs;
00475
00476 ene *= eV;
00477 xs *= cm2;
00478 theVec->PutValue(i,std::log(ene),std::log(xs));
00479 if (file.eof() && i != (nPoints-1))
00480 {
00481 G4ExceptionDescription ed ;
00482 ed << "Corrupted data file for Z=" << Z << G4endl;
00483 ed << "Found less than " << nPoints << "entries " <<G4endl;
00484 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00485 "em0005",FatalException,ed);
00486 }
00487 }
00488 if (!logAtomicCrossSection)
00489 {
00490 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00491 "em2044",FatalException,"Unable to allocate the atomic cross section table");
00492 delete theVec;
00493 return;
00494 }
00495 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
00496 file.close();
00497
00498
00499
00500
00501 std::ostringstream ost2;
00502 if (Z>9)
00503 ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
00504 else
00505 ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
00506 file.open(ost2.str().c_str());
00507 if (!file.is_open())
00508 {
00509 G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
00510 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00511 "em0003",FatalException,excep);
00512 }
00513 file >> readZ >> nPoints;
00514
00515 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
00516 {
00517 G4ExceptionDescription ed;
00518 ed << "Corrupted data file for Z=" << Z << G4endl;
00519 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00520 "em0005",FatalException,ed);
00521 return;
00522 }
00523 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
00524 G4double q=0,ff=0,incoh=0;
00525 G4bool fillQGrid = false;
00526
00527 if (!logQSquareGrid.size())
00528 fillQGrid = true;
00529 for (size_t i=0;i<nPoints;i++)
00530 {
00531 file >> q >> ff >> incoh;
00532
00533 theFFVec->PutValue(i,q,ff);
00534 if (fillQGrid)
00535 {
00536 logQSquareGrid.push_back(2.0*std::log(q));
00537 }
00538 if (file.eof() && i != (nPoints-1))
00539 {
00540 G4ExceptionDescription ed;
00541 ed << "Corrupted data file for Z=" << Z << G4endl;
00542 ed << "Found less than " << nPoints << "entries " <<G4endl;
00543 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00544 "em0005",FatalException,ed);
00545 }
00546 }
00547 if (!atomicFormFactor)
00548 {
00549 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
00550 "em2045",FatalException,
00551 "Unable to allocate the atomicFormFactor data table");
00552 delete theFFVec;
00553 return;
00554 }
00555 atomicFormFactor->insert(std::make_pair(Z,theFFVec));
00556 file.close();
00557 return;
00558 }
00559
00560
00561
00562 G4double G4PenelopeRayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared)
00563 {
00564 G4double f2 = 0;
00565
00566
00567
00568 G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.;
00569
00570 G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
00571
00572 if (!logFormFactorTable->count(mat))
00573 BuildFormFactorTable(mat);
00574
00575
00576 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
00577
00578 if (!theVec)
00579 {
00580 G4ExceptionDescription ed;
00581 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
00582 G4Exception("G4PenelopeRayleighModel::GetFSquared()",
00583 "em2046",FatalException,ed);
00584 return 0;
00585 }
00586 if (logQSquared < -20)
00587 {
00588 G4double logf2 = (*theVec)[0];
00589 f2 = std::exp(logf2);
00590 }
00591 else if (logQSquared > maxlogQ2)
00592 f2 =0;
00593 else
00594 {
00595
00596 G4double logf2 = theVec->Value(logQSquared);
00597 f2 = std::exp(logf2);
00598
00599 }
00600 if (verboseLevel > 3)
00601 {
00602 G4cout << "G4PenelopeRayleighModel::GetFSquared() in " << mat->GetName() << G4endl;
00603 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl;
00604 }
00605 return f2;
00606 }
00607
00608
00609
00610 void G4PenelopeRayleighModel::InitializeSamplingAlgorithm(const G4Material* mat)
00611 {
00612 G4double q2min = 0;
00613 G4double q2max = 0;
00614 const size_t np = 150;
00615 for (size_t i=1;i<logQSquareGrid.size();i++)
00616 {
00617 G4double Q2 = std::exp(logQSquareGrid[i]);
00618 if (GetFSquared(mat,Q2) > 1e-35)
00619 {
00620 q2max = std::exp(logQSquareGrid[i-1]);
00621 }
00622 }
00623
00624 size_t nReducedPoints = np/4;
00625
00626
00627 if (np < 16)
00628 {
00629 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
00630 "em2047",FatalException,
00631 "Too few points to initialize the sampling algorithm");
00632 }
00633 if (q2min > (q2max-1e-10))
00634 {
00635 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
00636 "em2048",FatalException,
00637 "Too narrow grid to initialize the sampling algorithm");
00638 }
00639
00640
00641
00642
00643
00644 G4DataVector* x = new G4DataVector();
00645
00646
00647
00648
00649 size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
00650 const G4int nip = 51;
00651
00652 G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
00653 x->push_back(q2min);
00654 for (size_t i=1;i<NUNIF-1;i++)
00655 {
00656 G4double app = q2min + i*dx;
00657 x->push_back(app);
00658 }
00659 x->push_back(q2max);
00660
00661 if (verboseLevel> 3)
00662 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
00663
00664
00665 G4DataVector* area = new G4DataVector();
00666 G4DataVector* a = new G4DataVector();
00667 G4DataVector* b = new G4DataVector();
00668 G4DataVector* c = new G4DataVector();
00669 G4DataVector* err = new G4DataVector();
00670
00671 for (size_t i=0;i<NUNIF-1;i++)
00672 {
00673
00674 G4DataVector* pdfi = new G4DataVector();
00675 G4DataVector* pdfih = new G4DataVector();
00676 G4DataVector* sumi = new G4DataVector();
00677
00678 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
00679 G4double pdfmax = 0;
00680 for (G4int k=0;k<nip;k++)
00681 {
00682 G4double xik = (*x)[i]+k*dxi;
00683 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
00684 pdfi->push_back(pdfk);
00685 pdfmax = std::max(pdfmax,pdfk);
00686 if (k < (nip-1))
00687 {
00688 G4double xih = xik + 0.5*dxi;
00689 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
00690 pdfih->push_back(pdfIK);
00691 pdfmax = std::max(pdfmax,pdfIK);
00692 }
00693 }
00694
00695
00696 G4double cons = dxi*0.5*(1./3.);
00697 sumi->push_back(0.);
00698 for (G4int k=1;k<nip;k++)
00699 {
00700 G4double previous = (*sumi)[k-1];
00701 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
00702 sumi->push_back(next);
00703 }
00704
00705 G4double lastIntegral = (*sumi)[sumi->size()-1];
00706 area->push_back(lastIntegral);
00707
00708 G4double factor = 1.0/lastIntegral;
00709 for (size_t k=0;k<sumi->size();k++)
00710 (*sumi)[k] *= factor;
00711
00712
00713 if ((*pdfi)[0] < 1e-35)
00714 (*pdfi)[0] = 1e-5*pdfmax;
00715 if ((*pdfi)[pdfi->size()-1] < 1e-35)
00716 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
00717
00718 G4double pli = (*pdfi)[0]*factor;
00719 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
00720 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
00721 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
00722 G4double C_temp = 1.0+A_temp+B_temp;
00723 if (C_temp < 1e-35)
00724 {
00725 a->push_back(0.);
00726 b->push_back(0.);
00727 c->push_back(1.);
00728 }
00729 else
00730 {
00731 a->push_back(A_temp);
00732 b->push_back(B_temp);
00733 c->push_back(C_temp);
00734 }
00735
00736
00737
00738 G4int icase = 1;
00739 G4bool reLoop = false;
00740 err->push_back(0.);
00741 do
00742 {
00743 reLoop = false;
00744 (*err)[i] = 0.;
00745 for (G4int k=0;k<nip;k++)
00746 {
00747 G4double rr = (*sumi)[k];
00748 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
00749 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
00750 if (k == 0 || k == nip-1)
00751 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
00752 else
00753 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
00754 }
00755 (*err)[i] *= dxi;
00756
00757
00758 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
00759 {
00760 (*b)[i] = 0;
00761 (*a)[i] = 0;
00762 (*c)[i] = 1.;
00763 icase = 2;
00764 reLoop = true;
00765 }
00766 }while(reLoop);
00767
00768 delete pdfi;
00769 delete pdfih;
00770 delete sumi;
00771 }
00772
00773
00774 (*x)[x->size()-1] = q2max;
00775 a->push_back(0.);
00776 b->push_back(0.);
00777 c->push_back(0.);
00778 err->push_back(0.);
00779 area->push_back(0.);
00780
00781 if (x->size() != NUNIF || a->size() != NUNIF ||
00782 err->size() != NUNIF || area->size() != NUNIF)
00783 {
00784 G4ExceptionDescription ed;
00785 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
00786 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
00787 "em2049",FatalException,ed);
00788 }
00789
00790
00791
00792
00793
00794 do
00795 {
00796 G4double maxError = 0.0;
00797 size_t iErrMax = 0;
00798 for (size_t i=0;i<err->size()-2;i++)
00799 {
00800
00801 if ((*err)[i] > maxError)
00802 {
00803 maxError = (*err)[i];
00804 iErrMax = i;
00805 }
00806 }
00807
00808
00809 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
00810
00811 x->insert(x->begin()+iErrMax+1,newx);
00812
00813 area->insert(area->begin()+iErrMax+1,0.);
00814 a->insert(a->begin()+iErrMax+1,0.);
00815 b->insert(b->begin()+iErrMax+1,0.);
00816 c->insert(c->begin()+iErrMax+1,0.);
00817 err->insert(err->begin()+iErrMax+1,0.);
00818
00819
00820 for (size_t i=iErrMax;i<=iErrMax+1;i++)
00821 {
00822
00823 G4DataVector* pdfi = new G4DataVector();
00824 G4DataVector* pdfih = new G4DataVector();
00825 G4DataVector* sumi = new G4DataVector();
00826
00827 G4double dxLocal = (*x)[i+1]-(*x)[i];
00828 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
00829 G4double pdfmax = 0;
00830 for (G4int k=0;k<nip;k++)
00831 {
00832 G4double xik = (*x)[i]+k*dxi;
00833 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
00834 pdfi->push_back(pdfk);
00835 pdfmax = std::max(pdfmax,pdfk);
00836 if (k < (nip-1))
00837 {
00838 G4double xih = xik + 0.5*dxi;
00839 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
00840 pdfih->push_back(pdfIK);
00841 pdfmax = std::max(pdfmax,pdfIK);
00842 }
00843 }
00844
00845
00846 G4double cons = dxi*0.5*(1./3.);
00847 sumi->push_back(0.);
00848 for (G4int k=1;k<nip;k++)
00849 {
00850 G4double previous = (*sumi)[k-1];
00851 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
00852 sumi->push_back(next);
00853 }
00854 G4double lastIntegral = (*sumi)[sumi->size()-1];
00855 (*area)[i] = lastIntegral;
00856
00857
00858 G4double factor = 1.0/lastIntegral;
00859 for (size_t k=0;k<sumi->size();k++)
00860 (*sumi)[k] *= factor;
00861
00862
00863 if ((*pdfi)[0] < 1e-35)
00864 (*pdfi)[0] = 1e-5*pdfmax;
00865 if ((*pdfi)[pdfi->size()-1] < 1e-35)
00866 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
00867
00868 G4double pli = (*pdfi)[0]*factor;
00869 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
00870 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
00871 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
00872 G4double C_temp = 1.0+A_temp+B_temp;
00873 if (C_temp < 1e-35)
00874 {
00875 (*a)[i]= 0.;
00876 (*b)[i] = 0.;
00877 (*c)[i] = 1;
00878 }
00879 else
00880 {
00881 (*a)[i]= A_temp;
00882 (*b)[i] = B_temp;
00883 (*c)[i] = C_temp;
00884 }
00885
00886
00887 G4int icase = 1;
00888 G4bool reLoop = false;
00889 do
00890 {
00891 reLoop = false;
00892 (*err)[i] = 0.;
00893 for (G4int k=0;k<nip;k++)
00894 {
00895 G4double rr = (*sumi)[k];
00896 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
00897 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
00898 if (k == 0 || k == nip-1)
00899 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
00900 else
00901 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
00902 }
00903 (*err)[i] *= dxi;
00904
00905
00906 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
00907 {
00908 (*b)[i] = 0;
00909 (*a)[i] = 0;
00910 (*c)[i] = 1.;
00911 icase = 2;
00912 reLoop = true;
00913 }
00914 }while(reLoop);
00915 delete pdfi;
00916 delete pdfih;
00917 delete sumi;
00918 }
00919 }while(x->size() < np);
00920
00921 if (x->size() != np || a->size() != np ||
00922 err->size() != np || area->size() != np)
00923 {
00924 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
00925 "em2050",FatalException,
00926 "Problem in building the extended Table for Sampling: array dimensions do not match ");
00927 }
00928
00929
00930
00931
00932 G4double ws = 0;
00933 for (size_t i=0;i<np-1;i++)
00934 ws += (*area)[i];
00935 ws = 1.0/ws;
00936 G4double errMax = 0;
00937 for (size_t i=0;i<np-1;i++)
00938 {
00939 (*area)[i] *= ws;
00940 (*err)[i] *= ws;
00941 errMax = std::max(errMax,(*err)[i]);
00942 }
00943
00944
00945 G4DataVector* PAC = new G4DataVector();
00946 PAC->push_back(0.);
00947 for (size_t i=0;i<np-1;i++)
00948 {
00949 G4double previous = (*PAC)[i];
00950 PAC->push_back(previous+(*area)[i]);
00951 }
00952 (*PAC)[PAC->size()-1] = 1.;
00953
00954
00955
00956
00957
00958
00959 std::vector<size_t> *ITTL = new std::vector<size_t>;
00960 std::vector<size_t> *ITTU = new std::vector<size_t>;
00961
00962
00963 for (size_t i=0;i<np;i++)
00964 {
00965 ITTL->push_back(0);
00966 ITTU->push_back(0);
00967 }
00968
00969 G4double bin = 1.0/(np-1);
00970 (*ITTL)[0]=0;
00971 for (size_t i=1;i<(np-1);i++)
00972 {
00973 G4double ptst = i*bin;
00974 G4bool found = false;
00975 for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
00976 {
00977 if ((*PAC)[j] > ptst)
00978 {
00979 (*ITTL)[i] = j-1;
00980 (*ITTU)[i-1] = j;
00981 found = true;
00982 }
00983 }
00984 }
00985 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
00986 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
00987 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
00988
00989 if (ITTU->size() != np || ITTU->size() != np)
00990 {
00991 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
00992 "em2051",FatalException,
00993 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
00994 }
00995
00996
00997
00998
00999
01000 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
01001 for (size_t i=0;i<np;i++)
01002 {
01003 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
01004 }
01005
01006 if (verboseLevel > 2)
01007 {
01008 G4cout << "*************************************************************************" <<
01009 G4endl;
01010 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
01011 theTable->DumpTable();
01012 }
01013 samplingTable->insert(std::make_pair(mat,theTable));
01014
01015
01016
01017 delete x;
01018 delete a;
01019 delete b;
01020 delete c;
01021 delete err;
01022 delete area;
01023 delete PAC;
01024 delete ITTL;
01025 delete ITTU;
01026
01027
01028 return;
01029
01030 }
01031
01032
01033
01034 void G4PenelopeRayleighModel::GetPMaxTable(const G4Material* mat)
01035 {
01036 if (!pMaxTable)
01037 {
01038 G4cout << "G4PenelopeRayleighModel::BuildPMaxTable" << G4endl;
01039 G4cout << "Going to instanziate the pMaxTable !" << G4endl;
01040 G4cout << "That should _not_ be here! " << G4endl;
01041 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
01042 }
01043
01044 if (pMaxTable->count(mat))
01045 return;
01046
01047
01048 if (!samplingTable)
01049 {
01050 G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
01051 "em2052",FatalException,
01052 "SamplingTable is not properly instantiated");
01053 return;
01054 }
01055
01056 if (!samplingTable->count(mat))
01057 InitializeSamplingAlgorithm(mat);
01058
01059 G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
01060 size_t tablePoints = theTable->GetNumberOfStoredPoints();
01061
01062 size_t nOfEnergyPoints = logEnergyGridPMax.size();
01063 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
01064
01065 const size_t nip = 51;
01066
01067 for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
01068 {
01069 G4double energy = std::exp(logEnergyGridPMax[ie]);
01070 G4double Qm = 2.0*energy/electron_mass_c2;
01071 G4double Qm2 = Qm*Qm;
01072 G4double firstQ2 = theTable->GetX(0);
01073 G4double lastQ2 = theTable->GetX(tablePoints-1);
01074 G4double thePMax = 0;
01075
01076 if (Qm2 > firstQ2)
01077 {
01078 if (Qm2 < lastQ2)
01079 {
01080
01081 size_t lowerBound = 0;
01082 size_t upperBound = tablePoints-1;
01083 while (lowerBound <= upperBound)
01084 {
01085 size_t midBin = (lowerBound + upperBound)/2;
01086 if( Qm2 < theTable->GetX(midBin))
01087 { upperBound = midBin-1; }
01088 else
01089 { lowerBound = midBin+1; }
01090 }
01091
01092 G4double Q1 = theTable->GetX(upperBound);
01093 G4double Q2 = Qm2;
01094 G4double DQ = (Q2-Q1)/((G4double)(nip-1));
01095 G4double theA = theTable->GetA(upperBound);
01096 G4double theB = theTable->GetB(upperBound);
01097 G4double thePAC = theTable->GetPAC(upperBound);
01098 G4DataVector* fun = new G4DataVector();
01099 for (size_t k=0;k<nip;k++)
01100 {
01101 G4double qi = Q1 + k*DQ;
01102 G4double tau = (qi-Q1)/
01103 (theTable->GetX(upperBound+1)-Q1);
01104 G4double con1 = 2.0*theB*tau;
01105 G4double ci = 1.0+theA+theB;
01106 G4double con2 = ci-theA*tau;
01107 G4double etap = 0;
01108 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
01109 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
01110 else
01111 etap = tau/con2;
01112 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
01113 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
01114 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
01115 fun->push_back(theFun);
01116 }
01117
01118 G4DataVector* sum = new G4DataVector;
01119 G4double CONS = DQ*(1./12.);
01120 G4double HCONS = 0.5*CONS;
01121 sum->push_back(0.);
01122 G4double secondPoint = (*sum)[0] +
01123 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
01124 sum->push_back(secondPoint);
01125 for (size_t hh=2;hh<nip-1;hh++)
01126 {
01127 G4double previous = (*sum)[hh-1];
01128 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
01129 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
01130 sum->push_back(next);
01131 }
01132 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
01133 (*fun)[nip-3])*CONS;
01134 sum->push_back(last);
01135 thePMax = thePAC + (*sum)[sum->size()-1];
01136 delete fun;
01137 delete sum;
01138 }
01139 else
01140 {
01141 thePMax = 1.0;
01142 }
01143 }
01144 else
01145 {
01146 thePMax = theTable->GetPAC(0);
01147 }
01148
01149
01150 theVec->PutValue(ie,energy,thePMax);
01151 }
01152
01153 pMaxTable->insert(std::make_pair(mat,theVec));
01154 return;
01155
01156 }
01157
01158
01159
01160 void G4PenelopeRayleighModel::DumpFormFactorTable(const G4Material* mat)
01161 {
01162 G4cout << "*****************************************************************" << G4endl;
01163 G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
01164
01165 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
01166 G4cout << "*****************************************************************" << G4endl;
01167 if (!logFormFactorTable->count(mat))
01168 BuildFormFactorTable(mat);
01169
01170 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
01171 for (size_t i=0;i<theVec->GetVectorLength();i++)
01172 {
01173 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
01174 G4double Q = std::exp(0.5*logQ2);
01175 G4double logF2 = (*theVec)[i];
01176 G4double F = std::exp(0.5*logF2);
01177 G4cout << Q << " " << F << G4endl;
01178 }
01179
01180 return;
01181 }