00001
00002
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
00044
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00069
00070 #include "G4ios.hh"
00071 #include "G4OpProcessSubType.hh"
00072
00073 #include "DsG4OpRayleigh.h"
00074
00076
00078
00080
00082
00083
00084
00085
00086
00088
00090
00091 DsG4OpRayleigh::DsG4OpRayleigh(const G4String& processName, G4ProcessType type)
00092 : G4VDiscreteProcess(processName, type)
00093 {
00094 SetProcessSubType(fOpRayleigh);
00095
00096 thePhysicsTable = 0;
00097
00098 DefaultWater = false;
00099
00100 if (verboseLevel>0) {
00101 G4cout << GetProcessName() << " is created " << G4endl;
00102 }
00103
00104 BuildThePhysicsTable();
00105 }
00106
00107
00108
00109
00110
00112
00114
00115 DsG4OpRayleigh::~DsG4OpRayleigh()
00116 {
00117 if (thePhysicsTable!= 0) {
00118 thePhysicsTable->clearAndDestroy();
00119 delete thePhysicsTable;
00120 }
00121 }
00122
00124
00126
00127
00128
00129
00130 G4VParticleChange*
00131 DsG4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00132 {
00133 aParticleChange.Initialize(aTrack);
00134
00135 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00136
00137 if (verboseLevel>0) {
00138 G4cout << "Scattering Photon!" << G4endl;
00139 G4cout << "Old Momentum Direction: "
00140 << aParticle->GetMomentumDirection() << G4endl;
00141 G4cout << "Old Polarization: "
00142 << aParticle->GetPolarization() << G4endl;
00143 }
00144
00145
00146
00147 G4double rand = G4UniformRand();
00148
00149 G4double CosTheta = std::pow(rand, 1./3.);
00150 G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
00151
00152 if(G4UniformRand() < 0.5)CosTheta = -CosTheta;
00153
00154
00155
00156 rand = G4UniformRand();
00157
00158 G4double Phi = twopi*rand;
00159 G4double SinPhi = std::sin(Phi);
00160 G4double CosPhi = std::cos(Phi);
00161
00162 G4double unit_x = SinTheta * CosPhi;
00163 G4double unit_y = SinTheta * SinPhi;
00164 G4double unit_z = CosTheta;
00165
00166 G4ThreeVector NewPolarization (unit_x,unit_y,unit_z);
00167
00168
00169
00170 G4ThreeVector OldPolarization = aParticle->GetPolarization();
00171 OldPolarization = OldPolarization.unit();
00172
00173 NewPolarization.rotateUz(OldPolarization);
00174 NewPolarization = NewPolarization.unit();
00175
00176
00177
00178
00179
00180 G4ThreeVector NewMomentumDirection =
00181 OldPolarization - NewPolarization * CosTheta;
00182
00183 if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection;
00184 NewMomentumDirection = NewMomentumDirection.unit();
00185
00186 aParticleChange.ProposePolarization(NewPolarization);
00187
00188 aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
00189
00190 if (verboseLevel>0) {
00191 G4cout << "New Polarization: "
00192 << NewPolarization << G4endl;
00193 G4cout << "Polarization Change: "
00194 << *(aParticleChange.GetPolarization()) << G4endl;
00195 G4cout << "New Momentum Direction: "
00196 << NewMomentumDirection << G4endl;
00197 G4cout << "Momentum Change: "
00198 << *(aParticleChange.GetMomentumDirection()) << G4endl;
00199 }
00200
00201 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00202 }
00203
00204
00205
00206
00207 void DsG4OpRayleigh::BuildThePhysicsTable()
00208 {
00209
00210
00211 if (thePhysicsTable) return;
00212
00213 const G4MaterialTable* theMaterialTable=
00214 G4Material::GetMaterialTable();
00215 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00216
00217
00218
00219 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
00220
00221
00222
00223 for (G4int i=0 ; i < numOfMaterials; i++)
00224 {
00225 G4PhysicsOrderedFreeVector* ScatteringLengths =
00226 new G4PhysicsOrderedFreeVector();
00227
00228 G4MaterialPropertiesTable *aMaterialPropertiesTable =
00229 (*theMaterialTable)[i]->GetMaterialPropertiesTable();
00230
00231 if(aMaterialPropertiesTable){
00232
00233 G4MaterialPropertyVector* AttenuationLengthVector =
00234 aMaterialPropertiesTable->GetProperty("RAYLEIGH");
00235
00236 if(!AttenuationLengthVector){
00237
00238 if ((*theMaterialTable)[i]->GetName() == "/dd/Materials/Water" ||
00239 (*theMaterialTable)[i]->GetName() == "/dd/Materials/OwsWater" ||
00240 (*theMaterialTable)[i]->GetName() == "/dd/Materials/IwsWater"
00241 )
00242 {
00243
00244
00245
00246 DefaultWater = true;
00247
00248 ScatteringLengths =
00249 RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
00250 }
00251 }
00252 }
00253
00254 thePhysicsTable->insertAt(i,ScatteringLengths);
00255 }
00256 }
00257
00258
00259
00260
00261 G4double DsG4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
00262 G4double ,
00263 G4ForceCondition* )
00264 {
00265 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00266 const G4Material* aMaterial = aTrack.GetMaterial();
00267
00268 G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00269
00270 G4double AttenuationLength = DBL_MAX;
00271
00272 if ((strcmp(aMaterial->GetName(), "/dd/Materials/Water") == 0 ||
00273 strcmp(aMaterial->GetName(), "/dd/Materials/OwsWater") == 0 ||
00274 strcmp(aMaterial->GetName(), "/dd/Materials/IwsWater") == 0 )
00275 && DefaultWater){
00276
00277 G4bool isOutRange;
00278
00279 AttenuationLength =
00280 (*thePhysicsTable)(aMaterial->GetIndex())->
00281 GetValue(thePhotonEnergy, isOutRange);
00282 }
00283 else {
00284
00285 G4MaterialPropertiesTable* aMaterialPropertyTable =
00286 aMaterial->GetMaterialPropertiesTable();
00287
00288 if(aMaterialPropertyTable){
00289 G4MaterialPropertyVector* AttenuationLengthVector =
00290 aMaterialPropertyTable->GetProperty("RAYLEIGH");
00291 if(AttenuationLengthVector){
00292 AttenuationLength = AttenuationLengthVector ->
00293 GetProperty(thePhotonEnergy);
00294 }
00295 else{
00296
00297 }
00298 }
00299 else{
00300
00301 }
00302 }
00303
00304 return AttenuationLength;
00305 }
00306
00307
00308
00309
00310
00311 G4PhysicsOrderedFreeVector*
00312 DsG4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT)
00313 {
00314
00315
00316
00317 G4double betat = 7.658e-23*m3/MeV;
00318
00319
00320 G4double kboltz = 8.61739e-11*MeV/kelvin;
00321
00322
00323
00324
00325 G4double temp = 283.15*kelvin;
00326
00327
00328
00329
00330 G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
00331
00332 G4double refsq;
00333 G4double e;
00334 G4double xlambda;
00335 G4double c1, c2, c3, c4;
00336 G4double Dist;
00337 G4double refraction_index;
00338
00339 G4PhysicsOrderedFreeVector *RayleighScatteringLengths =
00340 new G4PhysicsOrderedFreeVector();
00341
00342 if (Rindex ) {
00343
00344 Rindex->ResetIterator();
00345
00346 while (++(*Rindex)) {
00347
00348 e = (Rindex->GetPhotonEnergy());
00349
00350 refraction_index = Rindex->GetProperty();
00351 refsq = refraction_index*refraction_index;
00352 xlambda = h_Planck*c_light/e;
00353
00354 if (verboseLevel>0) {
00355 G4cout << Rindex->GetPhotonEnergy() << " MeV\t";
00356 G4cout << xlambda << " mm\t";
00357 }
00358
00359 c1 = 1 / (6.0 * pi);
00360 c2 = std::pow((2.0 * pi / xlambda), 4);
00361 c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
00362 c4 = betat * temp * kboltz;
00363
00364 Dist = 1.0 / (c1*c2*c3*c4);
00365
00366 if (verboseLevel>0) {
00367 G4cout << Dist << " mm" << G4endl;
00368 }
00369 RayleighScatteringLengths->
00370 InsertValues(Rindex->GetPhotonEnergy(), Dist);
00371 }
00372
00373 }
00374
00375 return RayleighScatteringLengths;
00376 }