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 #include "G4LivermorePolarizedComptonModel.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044
00045
00046
00047 using namespace std;
00048
00049
00050
00051 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
00052 const G4String& nam)
00053 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00054 meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
00055 {
00056 lowEnergyLimit = 250 * eV;
00057 highEnergyLimit = 100 * GeV;
00058
00059 SetHighEnergyLimit(highEnergyLimit);
00060
00061 verboseLevel= 0;
00062
00063
00064
00065
00066
00067
00068
00069 if( verboseLevel>0 ) {
00070 G4cout << "Livermore Polarized Compton is constructed " << G4endl
00071 << "Energy range: "
00072 << lowEnergyLimit / eV << " eV - "
00073 << highEnergyLimit / GeV << " GeV"
00074 << G4endl;
00075 }
00076 }
00077
00078
00079
00080 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
00081 {
00082 if (meanFreePathTable) delete meanFreePathTable;
00083 if (crossSectionHandler) delete crossSectionHandler;
00084 if (scatterFunctionData) delete scatterFunctionData;
00085 }
00086
00087
00088
00089 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
00090 const G4DataVector& cuts)
00091 {
00092 if (verboseLevel > 3)
00093 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
00094
00095 if (crossSectionHandler)
00096 {
00097 crossSectionHandler->Clear();
00098 delete crossSectionHandler;
00099 }
00100
00101
00102
00103 crossSectionHandler = new G4CrossSectionHandler;
00104 crossSectionHandler->Clear();
00105 G4String crossSectionFile = "comp/ce-cs-";
00106 crossSectionHandler->LoadData(crossSectionFile);
00107
00108 meanFreePathTable = 0;
00109 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
00110
00111 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
00112 G4String scatterFile = "comp/ce-sf-";
00113 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
00114 scatterFunctionData->LoadData(scatterFile);
00115
00116
00117 shellData.SetOccupancyData();
00118 G4String file = "/doppler/shell-doppler";
00119 shellData.LoadData(file);
00120
00121 if (verboseLevel > 2)
00122 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
00123
00124 InitialiseElementSelectors(particle,cuts);
00125
00126 if( verboseLevel>0 ) {
00127 G4cout << "Livermore Polarized Compton model is initialized " << G4endl
00128 << "Energy range: "
00129 << LowEnergyLimit() / eV << " eV - "
00130 << HighEnergyLimit() / GeV << " GeV"
00131 << G4endl;
00132 }
00133
00134
00135
00136 if(isInitialised) return;
00137 fParticleChange = GetParticleChangeForGamma();
00138 isInitialised = true;
00139 }
00140
00141
00142
00143 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
00144 const G4ParticleDefinition*,
00145 G4double GammaEnergy,
00146 G4double Z, G4double,
00147 G4double, G4double)
00148 {
00149 if (verboseLevel > 3)
00150 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
00151
00152 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
00153
00154 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00155 return cs;
00156 }
00157
00158
00159
00160 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00161 const G4MaterialCutsCouple* couple,
00162 const G4DynamicParticle* aDynamicGamma,
00163 G4double,
00164 G4double)
00165 {
00166
00167
00168
00169
00170
00171
00172 if (verboseLevel > 3)
00173 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
00174
00175 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
00176 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
00177
00178
00179
00180
00181
00182 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
00183
00184
00185
00186
00187 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
00188 {
00189 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
00190 }
00191 else
00192 {
00193 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
00194 {
00195 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
00196 }
00197 }
00198
00199
00200
00201
00202
00203 if(gammaEnergy0 <= lowEnergyLimit)
00204 {
00205 fParticleChange->ProposeTrackStatus(fStopAndKill);
00206 fParticleChange->SetProposedKineticEnergy(0.);
00207 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
00208 return;
00209 }
00210
00211 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
00212
00213
00214
00215 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00216 const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
00217 G4int Z = (G4int)elm->GetZ();
00218
00219
00220
00221 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
00222
00223 G4double epsilon0Local = 1./(1. + 2*E0_m);
00224 G4double epsilon0Sq = epsilon0Local*epsilon0Local;
00225 G4double alpha1 = - std::log(epsilon0Local);
00226 G4double alpha2 = 0.5*(1.- epsilon0Sq);
00227
00228 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
00229 G4double gammaEnergy1;
00230 G4ThreeVector gammaDirection1;
00231
00232 do {
00233 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
00234 {
00235 epsilon = std::exp(-alpha1*G4UniformRand());
00236 epsilonSq = epsilon*epsilon;
00237 }
00238 else
00239 {
00240 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
00241 epsilon = std::sqrt(epsilonSq);
00242 }
00243
00244 onecost = (1.- epsilon)/(epsilon*E0_m);
00245 sinThetaSqr = onecost*(2.-onecost);
00246
00247
00248 if (sinThetaSqr > 1.)
00249 {
00250 G4cout
00251 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00252 << "sin(theta)**2 = "
00253 << sinThetaSqr
00254 << "; set to 1"
00255 << G4endl;
00256 sinThetaSqr = 1.;
00257 }
00258 if (sinThetaSqr < 0.)
00259 {
00260 G4cout
00261 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00262 << "sin(theta)**2 = "
00263 << sinThetaSqr
00264 << "; set to 0"
00265 << G4endl;
00266 sinThetaSqr = 0.;
00267 }
00268
00269
00270 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
00271 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
00272 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
00273
00274 } while(greject < G4UniformRand()*Z);
00275
00276
00277
00278
00279
00280
00281 G4double phi = SetPhi(epsilon,sinThetaSqr);
00282
00283
00284
00285
00286
00287 G4double cosTheta = 1. - onecost;
00288
00289
00290
00291 if (cosTheta > 1.)
00292 {
00293 G4cout
00294 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00295 << "cosTheta = "
00296 << cosTheta
00297 << "; set to 1"
00298 << G4endl;
00299 cosTheta = 1.;
00300 }
00301 if (cosTheta < -1.)
00302 {
00303 G4cout
00304 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00305 << "cosTheta = "
00306 << cosTheta
00307 << "; set to -1"
00308 << G4endl;
00309 cosTheta = -1.;
00310 }
00311
00312
00313
00314 G4double sinTheta = std::sqrt (sinThetaSqr);
00315
00316
00317 if (sinTheta > 1.)
00318 {
00319 G4cout
00320 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00321 << "sinTheta = "
00322 << sinTheta
00323 << "; set to 1"
00324 << G4endl;
00325 sinTheta = 1.;
00326 }
00327 if (sinTheta < -1.)
00328 {
00329 G4cout
00330 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00331 << "sinTheta = "
00332 << sinTheta
00333 << "; set to -1"
00334 << G4endl;
00335 sinTheta = -1.;
00336 }
00337
00338
00339
00340 G4double dirx = sinTheta*std::cos(phi);
00341 G4double diry = sinTheta*std::sin(phi);
00342 G4double dirz = cosTheta ;
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 G4int maxDopplerIterations = 1000;
00355 G4double bindingE = 0.;
00356 G4double photonEoriginal = epsilon * gammaEnergy0;
00357 G4double photonE = -1.;
00358 G4int iteration = 0;
00359 G4double eMax = gammaEnergy0;
00360
00361 do
00362 {
00363 iteration++;
00364
00365 G4int shell = shellData.SelectRandomShell(Z);
00366 bindingE = shellData.BindingEnergy(Z,shell);
00367
00368 eMax = gammaEnergy0 - bindingE;
00369
00370
00371 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
00372
00373 G4double pDoppler = pSample * fine_structure_const;
00374 G4double pDoppler2 = pDoppler * pDoppler;
00375 G4double var2 = 1. + onecost * E0_m;
00376 G4double var3 = var2*var2 - pDoppler2;
00377 G4double var4 = var2 - pDoppler2 * cosTheta;
00378 G4double var = var4*var4 - var3 + pDoppler2 * var3;
00379 if (var > 0.)
00380 {
00381 G4double varSqrt = std::sqrt(var);
00382 G4double scale = gammaEnergy0 / var3;
00383
00384 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
00385 else photonE = (var4 + varSqrt) * scale;
00386 }
00387 else
00388 {
00389 photonE = -1.;
00390 }
00391 } while ( iteration <= maxDopplerIterations &&
00392 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
00393
00394
00395
00396 if (iteration >= maxDopplerIterations)
00397 {
00398 photonE = photonEoriginal;
00399 bindingE = 0.;
00400 }
00401
00402 gammaEnergy1 = photonE;
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
00414 sinThetaSqr,
00415 phi,
00416 cosTheta);
00417
00418
00419 G4ThreeVector tmpDirection1( dirx,diry,dirz );
00420 gammaDirection1 = tmpDirection1;
00421
00422
00423
00424 SystemOfRefChange(gammaDirection0,gammaDirection1,
00425 gammaPolarization0,gammaPolarization1);
00426
00427 if (gammaEnergy1 > 0.)
00428 {
00429 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
00430 fParticleChange->ProposeMomentumDirection( gammaDirection1 );
00431 fParticleChange->ProposePolarization( gammaPolarization1 );
00432 }
00433 else
00434 {
00435 gammaEnergy1 = 0.;
00436 fParticleChange->SetProposedKineticEnergy(0.) ;
00437 fParticleChange->ProposeTrackStatus(fStopAndKill);
00438 }
00439
00440
00441
00442
00443
00444 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
00445
00446
00447
00448 if(ElecKineEnergy < 0.0) {
00449 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
00450 return;
00451 }
00452
00453
00454
00455 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
00456
00457 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
00458 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
00459
00460 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
00461
00462 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
00463 fvect->push_back(dp);
00464
00465 }
00466
00467
00468
00469 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
00470 G4double sinSqrTh)
00471 {
00472 G4double rand1;
00473 G4double rand2;
00474 G4double phiProbability;
00475 G4double phi;
00476 G4double a, b;
00477
00478 do
00479 {
00480 rand1 = G4UniformRand();
00481 rand2 = G4UniformRand();
00482 phiProbability=0.;
00483 phi = twopi*rand1;
00484
00485 a = 2*sinSqrTh;
00486 b = energyRate + 1/energyRate;
00487
00488 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
00489
00490
00491
00492 }
00493 while ( rand2 > phiProbability );
00494 return phi;
00495 }
00496
00497
00498
00499
00500 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
00501 {
00502 G4double dx = a.x();
00503 G4double dy = a.y();
00504 G4double dz = a.z();
00505 G4double x = dx < 0.0 ? -dx : dx;
00506 G4double y = dy < 0.0 ? -dy : dy;
00507 G4double z = dz < 0.0 ? -dz : dz;
00508 if (x < y) {
00509 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
00510 }else{
00511 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
00512 }
00513 }
00514
00515
00516
00517 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
00518 {
00519 G4ThreeVector d0 = direction0.unit();
00520 G4ThreeVector a1 = SetPerpendicularVector(d0);
00521 G4ThreeVector a0 = a1.unit();
00522
00523 G4double rand1 = G4UniformRand();
00524
00525 G4double angle = twopi*rand1;
00526 G4ThreeVector b0 = d0.cross(a0);
00527
00528 G4ThreeVector c;
00529
00530 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
00531 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
00532 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
00533
00534 G4ThreeVector c0 = c.unit();
00535
00536 return c0;
00537
00538 }
00539
00540
00541
00542 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
00543 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
00544 {
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
00557 }
00558
00559
00560
00561 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
00562 G4double sinSqrTh,
00563 G4double phi,
00564 G4double costheta)
00565 {
00566 G4double rand1;
00567 G4double rand2;
00568 G4double cosPhi = std::cos(phi);
00569 G4double sinPhi = std::sin(phi);
00570 G4double sinTheta = std::sqrt(sinSqrTh);
00571 G4double cosSqrPhi = cosPhi*cosPhi;
00572
00573
00574 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
00575
00576
00577
00578
00579
00580
00581
00582 G4double theta;
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610 rand1 = G4UniformRand();
00611 rand2 = G4UniformRand();
00612
00613 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
00614 {
00615 if (rand2<0.5)
00616 theta = pi/2.0;
00617 else
00618 theta = 3.0*pi/2.0;
00619 }
00620 else
00621 {
00622 if (rand2<0.5)
00623 theta = 0;
00624 else
00625 theta = pi;
00626 }
00627 G4double cosBeta = std::cos(theta);
00628 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
00629
00630 G4ThreeVector gammaPolarization1;
00631
00632 G4double xParallel = normalisation*cosBeta;
00633 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
00634 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
00635 G4double xPerpendicular = 0.;
00636 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
00637 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
00638
00639 G4double xTotal = (xParallel + xPerpendicular);
00640 G4double yTotal = (yParallel + yPerpendicular);
00641 G4double zTotal = (zParallel + zPerpendicular);
00642
00643 gammaPolarization1.setX(xTotal);
00644 gammaPolarization1.setY(yTotal);
00645 gammaPolarization1.setZ(zTotal);
00646
00647 return gammaPolarization1;
00648
00649 }
00650
00651
00652
00653 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
00654 G4ThreeVector& direction1,
00655 G4ThreeVector& polarization0,
00656 G4ThreeVector& polarization1)
00657 {
00658
00659
00660
00661 G4ThreeVector Axis_Z0 = direction0.unit();
00662 G4ThreeVector Axis_X0 = polarization0.unit();
00663 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit();
00664
00665 G4double direction_x = direction1.getX();
00666 G4double direction_y = direction1.getY();
00667 G4double direction_z = direction1.getZ();
00668
00669 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
00670 G4double polarization_x = polarization1.getX();
00671 G4double polarization_y = polarization1.getY();
00672 G4double polarization_z = polarization1.getZ();
00673
00674 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
00675
00676 }
00677
00678