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 #include "G4LivermorePolarizedRayleighModel.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046
00047
00048
00049 using namespace std;
00050
00051
00052
00053 G4LivermorePolarizedRayleighModel::G4LivermorePolarizedRayleighModel(const G4ParticleDefinition*,
00054 const G4String& nam)
00055 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00056 crossSectionHandler(0),formFactorData(0)
00057 {
00058 lowEnergyLimit = 250 * eV;
00059 highEnergyLimit = 100 * GeV;
00060
00061
00062 SetHighEnergyLimit(highEnergyLimit);
00063
00064 verboseLevel= 0;
00065
00066
00067
00068
00069
00070
00071
00072 if(verboseLevel > 0) {
00073 G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
00074 << "Energy range: "
00075 << lowEnergyLimit / eV << " eV - "
00076 << highEnergyLimit / GeV << " GeV"
00077 << G4endl;
00078 }
00079 }
00080
00081
00082
00083 G4LivermorePolarizedRayleighModel::~G4LivermorePolarizedRayleighModel()
00084 {
00085 if (crossSectionHandler) delete crossSectionHandler;
00086 if (formFactorData) delete formFactorData;
00087 }
00088
00089
00090
00091 void G4LivermorePolarizedRayleighModel::Initialise(const G4ParticleDefinition* particle,
00092 const G4DataVector& cuts)
00093 {
00094
00095
00096
00097
00098
00099
00100
00101 if (verboseLevel > 3)
00102 G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
00103
00104 if (crossSectionHandler)
00105 {
00106 crossSectionHandler->Clear();
00107 delete crossSectionHandler;
00108 }
00109
00110
00111
00112 crossSectionHandler = new G4CrossSectionHandler;
00113 crossSectionHandler->Clear();
00114 G4String crossSectionFile = "rayl/re-cs-";
00115 crossSectionHandler->LoadData(crossSectionFile);
00116
00117 G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
00118 G4String formFactorFile = "rayl/re-ff-";
00119 formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
00120 formFactorData->LoadData(formFactorFile);
00121
00122 InitialiseElementSelectors(particle,cuts);
00123
00124
00125 if (verboseLevel > 2)
00126 G4cout << "Loaded cross section files for Livermore Polarized Rayleigh model" << G4endl;
00127
00128 InitialiseElementSelectors(particle,cuts);
00129
00130 if (verboseLevel > 0) {
00131 G4cout << "Livermore Polarized Rayleigh model is initialized " << G4endl
00132 << "Energy range: "
00133 << LowEnergyLimit() / eV << " eV - "
00134 << HighEnergyLimit() / GeV << " GeV"
00135 << G4endl;
00136 }
00137
00138 if(isInitialised) return;
00139 fParticleChange = GetParticleChangeForGamma();
00140 isInitialised = true;
00141 }
00142
00143
00144
00145 G4double G4LivermorePolarizedRayleighModel::ComputeCrossSectionPerAtom(
00146 const G4ParticleDefinition*,
00147 G4double GammaEnergy,
00148 G4double Z, G4double,
00149 G4double, G4double)
00150 {
00151 if (verboseLevel > 3)
00152 G4cout << "Calling CrossSectionPerAtom() of G4LivermorePolarizedRayleighModel" << G4endl;
00153
00154 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
00155
00156 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00157 return cs;
00158 }
00159
00160
00161
00162 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00163 const G4MaterialCutsCouple* couple,
00164 const G4DynamicParticle* aDynamicGamma,
00165 G4double,
00166 G4double)
00167 {
00168 if (verboseLevel > 3)
00169 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl;
00170
00171 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
00172
00173 if (photonEnergy0 <= lowEnergyLimit)
00174 {
00175 fParticleChange->ProposeTrackStatus(fStopAndKill);
00176 fParticleChange->SetProposedKineticEnergy(0.);
00177 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
00178 return ;
00179 }
00180
00181 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
00182
00183
00184
00185 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00186 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
00187 G4int Z = (G4int)elm->GetZ();
00188
00189 G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);
00190 G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta);
00191 G4double beta=GeneratePolarizationAngle();
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit());
00204 G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma));
00205 G4ThreeVector y(z.cross(x));
00206
00207
00208 G4double xDir;
00209 G4double yDir;
00210 G4double zDir;
00211 zDir=outcomingPhotonCosTheta;
00212 xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta);
00213 yDir=xDir;
00214 xDir*=std::cos(outcomingPhotonPhi);
00215 yDir*=std::sin(outcomingPhotonPhi);
00216
00217 G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit());
00218 G4ThreeVector xPrime(x.perpPart(zPrime).unit());
00219 G4ThreeVector yPrime(zPrime.cross(xPrime));
00220
00221
00222 G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta));
00223
00224 fParticleChange->ProposeMomentumDirection(zPrime);
00225 fParticleChange->ProposePolarization(outcomingPhotonPolarization);
00226 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
00227
00228 }
00229
00230
00231
00232 G4double G4LivermorePolarizedRayleighModel::GenerateCosTheta(G4double incomingPhotonEnergy, G4int zAtom) const
00233 {
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light);
00253
00254
00255 G4double cosTheta;
00256 G4double fCosTheta;
00257 G4double x;
00258 G4double fValue;
00259
00260 if (incomingPhotonEnergy > 5.*MeV)
00261 {
00262 cosTheta = 1.;
00263 }
00264 else
00265 {
00266 do
00267 {
00268 do
00269 {
00270 cosTheta = 2.*G4UniformRand()-1.;
00271 fCosTheta = (1.+cosTheta*cosTheta)/2.;
00272 }
00273 while (fCosTheta < G4UniformRand());
00274
00275 x = xFactor*std::sqrt((1.-cosTheta)/2.);
00276
00277 if (x > 1.e+005)
00278 fValue = formFactorData->FindValue(x, zAtom-1);
00279 else
00280 fValue = formFactorData->FindValue(0., zAtom-1);
00281
00282 fValue/=zAtom;
00283 fValue*=fValue;
00284 }
00285 while(fValue < G4UniformRand());
00286 }
00287
00288 return cosTheta;
00289 }
00290
00291
00292
00293 G4double G4LivermorePolarizedRayleighModel::GeneratePhi(G4double cosTheta) const
00294 {
00295
00296
00297
00298
00299
00300
00301 G4double phi;
00302 G4double cosPhi;
00303 G4double phiProbability;
00304 G4double sin2Theta;
00305
00306 sin2Theta=1.-cosTheta*cosTheta;
00307
00308 do
00309 {
00310 phi = twopi * G4UniformRand();
00311 cosPhi = std::cos(phi);
00312 phiProbability= 1. - sin2Theta*cosPhi*cosPhi;
00313 }
00314 while (phiProbability < G4UniformRand());
00315
00316 return phi;
00317 }
00318
00319
00320
00321 G4double G4LivermorePolarizedRayleighModel::GeneratePolarizationAngle(void) const
00322 {
00323
00324
00325 return 0;
00326 }
00327
00328
00329
00330 G4ThreeVector G4LivermorePolarizedRayleighModel::GetPhotonPolarization(const G4DynamicParticle& photon)
00331 {
00332
00333
00334
00335 G4ThreeVector photonMomentumDirection;
00336 G4ThreeVector photonPolarization;
00337
00338 photonPolarization = photon.GetPolarization();
00339 photonMomentumDirection = photon.GetMomentumDirection();
00340
00341 if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.)
00342 {
00343
00344
00345
00346 G4ThreeVector e1(photonMomentumDirection.orthogonal().unit());
00347 G4ThreeVector e2(photonMomentumDirection.cross(e1).unit());
00348
00349 G4double angle(G4UniformRand() * twopi);
00350
00351 e1*=std::cos(angle);
00352 e2*=std::sin(angle);
00353
00354 photonPolarization=e1+e2;
00355 }
00356 else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.)
00357 {
00358
00359
00360
00361 photonPolarization=photonPolarization.perpPart(photonMomentumDirection);
00362 }
00363
00364 return photonPolarization.unit();
00365 }
00366