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
00048 #include "G4SeltzerBergerModel.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4Electron.hh"
00052 #include "G4Positron.hh"
00053 #include "G4Gamma.hh"
00054 #include "Randomize.hh"
00055 #include "G4Material.hh"
00056 #include "G4Element.hh"
00057 #include "G4ElementVector.hh"
00058 #include "G4ProductionCutsTable.hh"
00059 #include "G4ParticleChangeForLoss.hh"
00060 #include "G4LossTableManager.hh"
00061 #include "G4ModifiedTsai.hh"
00062
00063 #include "G4Physics2DVector.hh"
00064
00065 #include "G4ios.hh"
00066 #include <fstream>
00067 #include <iomanip>
00068
00069
00070
00071
00072 using namespace std;
00073
00074 G4Physics2DVector* G4SeltzerBergerModel::dataSB[101] = {0};
00075 G4double G4SeltzerBergerModel::ylimit[101] = {0.0};
00076 G4double G4SeltzerBergerModel::expnumlim = -12.;
00077
00078 G4SeltzerBergerModel::G4SeltzerBergerModel(const G4ParticleDefinition* p,
00079 const G4String& nam)
00080 : G4eBremsstrahlungRelModel(p,nam),useBicubicInterpolation(false)
00081 {
00082 SetLowEnergyLimit(0.0);
00083 SetLPMFlag(false);
00084 nwarn = 0;
00085 }
00086
00087
00088
00089 G4SeltzerBergerModel::~G4SeltzerBergerModel()
00090 {
00091 for(size_t i=0; i<101; ++i) {
00092 if(dataSB[i]) {
00093 delete dataSB[i];
00094 dataSB[i] = 0;
00095 }
00096 }
00097 }
00098
00099
00100
00101 void G4SeltzerBergerModel::Initialise(const G4ParticleDefinition* p,
00102 const G4DataVector& cuts)
00103 {
00104
00105
00106 char* path = getenv("G4LEDATA");
00107
00108
00109 const G4ElementTable* theElmTable = G4Element::GetElementTable();
00110 size_t numOfElm = G4Element::GetNumberOfElements();
00111 if(numOfElm > 0) {
00112 for(size_t i=0; i<numOfElm; ++i) {
00113 G4int Z = G4int(((*theElmTable)[i])->GetZ());
00114 if(Z < 1) { Z = 1; }
00115 else if(Z > 100) { Z = 100; }
00116
00117
00118 if(!dataSB[Z]) { ReadData(Z, path); }
00119 }
00120 }
00121
00122 G4eBremsstrahlungRelModel::Initialise(p, cuts);
00123 }
00124
00125
00126
00127 void G4SeltzerBergerModel::ReadData(size_t Z, const char* path)
00128 {
00129
00130
00131
00132 if(dataSB[Z]) { return; }
00133 const char* datadir = path;
00134
00135 if(!datadir) {
00136 datadir = getenv("G4LEDATA");
00137 if(!datadir) {
00138 G4Exception("G4SeltzerBergerModel::ReadData()","em0006",FatalException,
00139 "Environment variable G4LEDATA not defined");
00140 return;
00141 }
00142 }
00143 std::ostringstream ost;
00144 ost << datadir << "/brem_SB/br" << Z;
00145 std::ifstream fin(ost.str().c_str());
00146 if( !fin.is_open()) {
00147 G4ExceptionDescription ed;
00148 ed << "Bremsstrahlung data file <" << ost.str().c_str()
00149 << "> is not opened!" << G4endl;
00150 G4Exception("G4SeltzerBergerModel::ReadData()","em0003",FatalException,
00151 ed,"G4LEDATA version should be G4EMLOW6.23 or later.");
00152 return;
00153 }
00154
00155
00156 G4Physics2DVector* v = new G4Physics2DVector();
00157 const G4double emaxlog = 4*log(10.);
00158 if(v->Retrieve(fin)) {
00159 if(useBicubicInterpolation) { v->SetBicubicInterpolation(true); }
00160 dataSB[Z] = v;
00161 ylimit[Z] = v->Value(0.97, emaxlog);
00162 } else {
00163 G4ExceptionDescription ed;
00164 ed << "Bremsstrahlung data file <" << ost.str().c_str()
00165 << "> is not retrieved!" << G4endl;
00166 G4Exception("G4SeltzerBergerModel::ReadData()","em0005",FatalException,
00167 ed,"G4LEDATA version should be G4EMLOW6.23 or later.");
00168 delete v;
00169 }
00170
00171 }
00172
00173
00174
00175 G4double G4SeltzerBergerModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
00176 {
00177
00178 if(gammaEnergy < 0.0 || kinEnergy <= 0.0) { return 0.0; }
00179 G4double x = gammaEnergy/kinEnergy;
00180 G4double y = log(kinEnergy/MeV);
00181 G4int Z = G4int(currentZ);
00182
00183
00184 if(!dataSB[Z]) { ReadData(Z); }
00185 G4double invb2 = totalEnergy*totalEnergy/(kinEnergy*(kinEnergy + 2*particleMass));
00186 G4double cross = dataSB[Z]->Value(x,y)*invb2*millibarn/bremFactor;
00187
00188 if(!isElectron) {
00189 G4double invbeta1 = sqrt(invb2);
00190 G4double e2 = kinEnergy - gammaEnergy;
00191 if(e2 > 0.0) {
00192 G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
00193 G4double xxx = twopi*fine_structure_const*currentZ*(invbeta1 - invbeta2);
00194 if(xxx < expnumlim) { cross = 0.0; }
00195 else { cross *= exp(xxx); }
00196 } else {
00197 cross = 0.0;
00198 }
00199 }
00200
00201 return cross;
00202 }
00203
00204
00205
00206 void
00207 G4SeltzerBergerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
00208 const G4MaterialCutsCouple* couple,
00209 const G4DynamicParticle* dp,
00210 G4double cutEnergy,
00211 G4double maxEnergy)
00212 {
00213 G4double kineticEnergy = dp->GetKineticEnergy();
00214 G4double cut = std::min(cutEnergy, kineticEnergy);
00215 G4double emax = std::min(maxEnergy, kineticEnergy);
00216 if(cut >= emax) { return; }
00217
00218 SetupForMaterial(particle, couple->GetMaterial(), kineticEnergy);
00219
00220 const G4Element* elm =
00221 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
00222 SetCurrentElement(elm->GetZ());
00223 G4int Z = G4int(currentZ);
00224
00225 totalEnergy = kineticEnergy + particleMass;
00226 densityCorr = densityFactor*totalEnergy*totalEnergy;
00227 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
00228
00229
00230
00231
00232
00233
00234 G4double xmin = log(cut*cut + densityCorr);
00235 G4double xmax = log(emax*emax + densityCorr);
00236 G4double y = log(kineticEnergy/MeV);
00237
00238 G4double gammaEnergy, v;
00239
00240
00241 G4double x0 = cut/kineticEnergy;
00242 G4double vmax = dataSB[Z]->Value(x0, y)*1.02;
00243
00244
00245 const G4double epeaklimit= 300*MeV;
00246 const G4double elowlimit = 10*keV;
00247
00248
00249 if(isElectron && x0 < 0.97 &&
00250 ((kineticEnergy > epeaklimit) || (kineticEnergy < elowlimit))) {
00251 G4double ylim = std::min(ylimit[Z],1.1*dataSB[Z]->Value(0.97, y));
00252 if(ylim > vmax) { vmax = ylim; }
00253 }
00254 if(x0 < 0.05) { vmax *= 1.2; }
00255
00256
00257
00258 do {
00259
00260 G4double x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
00261 if(x < 0.0) { x = 0.0; }
00262 gammaEnergy = sqrt(x);
00263 G4double x1 = gammaEnergy/kineticEnergy;
00264 v = dataSB[Z]->Value(x1, y);
00265
00266
00267 if(!isElectron) {
00268 G4double e1 = kineticEnergy - cut;
00269 G4double invbeta1 = (e1 + particleMass)/sqrt(e1*(e1 + 2*particleMass));
00270 G4double e2 = kineticEnergy - gammaEnergy;
00271 G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
00272 G4double xxx = twopi*fine_structure_const*currentZ*(invbeta1 - invbeta2);
00273
00274 if(xxx < expnumlim) { v = 0.0; }
00275 else { v *= exp(xxx); }
00276 }
00277
00278 if (v > 1.05*vmax && nwarn < 20) {
00279 ++nwarn;
00280 G4cout << "### G4SeltzerBergerModel Warning: Majoranta exceeded! "
00281 << v << " > " << vmax << " by " << v/vmax
00282 << " Egamma(MeV)= " << gammaEnergy
00283 << " Ee(MeV)= " << kineticEnergy
00284 << " Z= " << Z << " " << particle->GetParticleName()
00285
00286 << G4endl;
00287
00288 if ( 20 == nwarn ) {
00289 G4cout << "### G4SeltzerBergerModel Warnings will not be printed anymore"
00290 << G4endl;
00291 }
00292 }
00293 } while (v < vmax*G4UniformRand());
00294
00295
00296
00297
00298
00299
00300 G4ThreeVector gammaDirection =
00301 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00302 Z, couple->GetMaterial());
00303
00304
00305 G4DynamicParticle* gamma =
00306 new G4DynamicParticle(theGamma,gammaDirection,gammaEnergy);
00307 vdp->push_back(gamma);
00308
00309 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
00310 - gammaEnergy*gammaDirection).unit();
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 G4double finalE = kineticEnergy - gammaEnergy;
00321
00322
00323 if(gammaEnergy > SecondaryThreshold()) {
00324 fParticleChange->ProposeTrackStatus(fStopAndKill);
00325 fParticleChange->SetProposedKineticEnergy(0.0);
00326 G4DynamicParticle* el =
00327 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
00328 direction, finalE);
00329 vdp->push_back(el);
00330
00331
00332 } else {
00333 fParticleChange->SetProposedMomentumDirection(direction);
00334 fParticleChange->SetProposedKineticEnergy(finalE);
00335 }
00336 }
00337
00338
00339
00340