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 #include "G4LivermoreGammaConversionModelRC.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043
00044
00045
00046 using namespace std;
00047
00048
00049
00050 G4LivermoreGammaConversionModelRC::G4LivermoreGammaConversionModelRC(const G4ParticleDefinition*,
00051 const G4String& nam)
00052 :G4VEmModel(nam),fParticleChange(0),smallEnergy(2.*MeV),isInitialised(false),
00053 crossSectionHandler(0),meanFreePathTable(0)
00054 {
00055 lowEnergyLimit = 2.0*electron_mass_c2;
00056 highEnergyLimit = 100 * GeV;
00057 SetHighEnergyLimit(highEnergyLimit);
00058
00059 verboseLevel= 0;
00060
00061
00062
00063
00064
00065
00066
00067 if(verboseLevel > 0) {
00068 G4cout << "Livermore Gamma conversion is constructed " << G4endl
00069 << "Energy range: "
00070 << lowEnergyLimit / MeV << " MeV - "
00071 << highEnergyLimit / GeV << " GeV"
00072 << G4endl;
00073 }
00074 }
00075
00076
00077
00078 G4LivermoreGammaConversionModelRC::~G4LivermoreGammaConversionModelRC()
00079 {
00080 if (crossSectionHandler) delete crossSectionHandler;
00081 }
00082
00083
00084
00085 void
00086 G4LivermoreGammaConversionModelRC::Initialise(const G4ParticleDefinition*,
00087 const G4DataVector&)
00088 {
00089 if (verboseLevel > 3)
00090 G4cout << "Calling G4LivermoreGammaConversionModelRC::Initialise()" << G4endl;
00091
00092 if (crossSectionHandler)
00093 {
00094 crossSectionHandler->Clear();
00095 delete crossSectionHandler;
00096 }
00097
00098
00099
00100 crossSectionHandler = new G4CrossSectionHandler();
00101 crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
00102 G4String crossSectionFile = "pair/pp-cs-";
00103 crossSectionHandler->LoadData(crossSectionFile);
00104
00105
00106
00107 if (verboseLevel > 2)
00108 G4cout << "Loaded cross section files for Livermore Gamma Conversion model RC" << G4endl;
00109
00110 if (verboseLevel > 0) {
00111 G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
00112 << "Energy range: "
00113 << LowEnergyLimit() / MeV << " MeV - "
00114 << HighEnergyLimit() / GeV << " GeV"
00115 << G4endl;
00116 }
00117
00118 if(isInitialised) return;
00119 fParticleChange = GetParticleChangeForGamma();
00120 isInitialised = true;
00121 }
00122
00123
00124
00125 G4double
00126 G4LivermoreGammaConversionModelRC::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00127 G4double GammaEnergy,
00128 G4double Z, G4double,
00129 G4double, G4double)
00130 {
00131 if (verboseLevel > 3) {
00132 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC"
00133 << G4endl;
00134 }
00135 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
00136
00137 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00138 return cs;
00139 }
00140
00141
00142
00143 void G4LivermoreGammaConversionModelRC::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00144 const G4MaterialCutsCouple* couple,
00145 const G4DynamicParticle* aDynamicGamma,
00146 G4double,
00147 G4double)
00148 {
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 if (verboseLevel > 3)
00161 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC" << G4endl;
00162
00163 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00164 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00165
00166 G4double epsilon ;
00167 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00168 G4double electronTotEnergy;
00169 G4double positronTotEnergy;
00170
00171
00172
00173 if (photonEnergy < smallEnergy )
00174 {
00175 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00176
00177 if (G4int(2*G4UniformRand()))
00178 {
00179 electronTotEnergy = (1. - epsilon) * photonEnergy;
00180 positronTotEnergy = epsilon * photonEnergy;
00181 }
00182 else
00183 {
00184 positronTotEnergy = (1. - epsilon) * photonEnergy;
00185 electronTotEnergy = epsilon * photonEnergy;
00186 }
00187 }
00188 else
00189 {
00190
00191
00192 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00193 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00194 G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries" << G4endl;
00195
00196 if (element == 0)
00197 {
00198 G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
00199 << G4endl;
00200 return;
00201 }
00202 G4IonisParamElm* ionisation = element->GetIonisation();
00203 if (ionisation == 0)
00204 {
00205 G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
00206 << G4endl;
00207 return;
00208 }
00209
00210
00211 G4double fZ = 8. * (ionisation->GetlogZ3());
00212 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00213
00214
00215 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00216 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
00217 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00218
00219
00220 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
00221 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00222 G4double epsilonRange = 0.5 - epsilonMin ;
00223
00224
00225 G4double screen;
00226 G4double gReject ;
00227
00228 G4double f10 = ScreenFunction1(screenMin) - fZ;
00229 G4double f20 = ScreenFunction2(screenMin) - fZ;
00230 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00231 G4double normF2 = std::max(1.5 * f20,0.);
00232 G4double a=393.3750918, b=115.3070201, c=810.6428451, d=19.96497475, e=1016.874592, f=1.936685510,
00233 gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
00234 G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287, ff=-117602.840,
00235 gg=-721455.467, hh=693957.8635, ii=156266.1085, jj=533209.9347;
00236 G4double Rechazo = 0.;
00237 G4double logepsMin = log(epsilonMin);
00238 G4double NormaRC = a + b*logepsMin + c/logepsMin + d*pow(logepsMin,2.) + e/pow(logepsMin,2.) + f*pow(logepsMin,3.) +
00239 gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
00240 k/pow(logepsMin,5.);
00241
00242 do {
00243 do {
00244 if (normF1 / (normF1 + normF2) > G4UniformRand() )
00245 {
00246 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
00247 screen = screenFactor / (epsilon * (1. - epsilon));
00248 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00249 }
00250 else
00251 {
00252 epsilon = epsilonMin + epsilonRange * G4UniformRand();
00253 screen = screenFactor / (epsilon * (1 - epsilon));
00254 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00255 }
00256 } while ( gReject < G4UniformRand() );
00257
00258 if (G4int(2*G4UniformRand())) epsilon = (1. - epsilon);
00259
00260 G4double logepsilon = log(epsilon);
00261 G4double deltaP_R1 = 1. + (a + b*logepsilon + c/logepsilon + d*pow(logepsilon,2.) + e/pow(logepsilon,2.) +
00262 f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
00263 j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
00264 G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
00265 / (1. + bb*logepsilon + dd*pow(logepsilon,2.) + ff*pow(logepsilon,3.) + hh*pow(logepsilon,4.)
00266 + jj*pow(logepsilon,5.) ))/100.;
00267
00268 if (epsilon <= 0.5)
00269 {
00270 Rechazo = deltaP_R1/NormaRC;
00271 }
00272 else
00273 {
00274 Rechazo = deltaP_R2/NormaRC;
00275 }
00276 G4cout << Rechazo << " " << NormaRC << " " << epsilon << G4endl;
00277 } while (Rechazo < G4UniformRand() );
00278
00279 electronTotEnergy = (1. - epsilon) * photonEnergy;
00280 positronTotEnergy = epsilon * photonEnergy;
00281
00282 }
00283
00284
00285
00286
00287
00288
00289
00290 G4double u;
00291 const G4double a1 = 0.625;
00292 G4double a2 = 3. * a1;
00293
00294
00295
00296 if (0.25 > G4UniformRand())
00297 {
00298 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
00299 }
00300 else
00301 {
00302 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
00303 }
00304
00305 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
00306 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
00307 G4double phi = twopi * G4UniformRand();
00308
00309 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00310 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00311
00312
00313
00314
00315
00316
00317 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00318
00319
00320
00321 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00322 electronDirection.rotateUz(photonDirection);
00323
00324 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00325 electronDirection,
00326 electronKineEnergy);
00327
00328
00329 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00330
00331
00332
00333 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00334 positronDirection.rotateUz(photonDirection);
00335
00336
00337 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00338 positronDirection, positronKineEnergy);
00339
00340
00341 fvect->push_back(particle1);
00342 fvect->push_back(particle2);
00343
00344
00345 fParticleChange->SetProposedKineticEnergy(0.);
00346 fParticleChange->ProposeTrackStatus(fStopAndKill);
00347
00348 }
00349
00350
00351
00352 G4double G4LivermoreGammaConversionModelRC::ScreenFunction1(G4double screenVariable)
00353 {
00354
00355
00356 G4double value;
00357
00358 if (screenVariable > 1.)
00359 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00360 else
00361 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00362
00363 return value;
00364 }
00365
00366
00367
00368 G4double G4LivermoreGammaConversionModelRC::ScreenFunction2(G4double screenVariable)
00369 {
00370
00371
00372 G4double value;
00373
00374 if (screenVariable > 1.)
00375 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00376 else
00377 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00378
00379 return value;
00380 }
00381