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 #include "G4LivermorePolarizedGammaConversionModel.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036
00037
00038
00039 using namespace std;
00040
00041
00042
00043 G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel(
00044 const G4ParticleDefinition*, const G4String& nam)
00045 :G4VEmModel(nam),fParticleChange(0),
00046 isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
00047 {
00048 lowEnergyLimit = 2*electron_mass_c2;
00049 highEnergyLimit = 100 * GeV;
00050 SetLowEnergyLimit(lowEnergyLimit);
00051 SetHighEnergyLimit(highEnergyLimit);
00052 smallEnergy = 2.*MeV;
00053
00054 Phi=0.;
00055 Psi=0.;
00056
00057 verboseLevel= 0;
00058
00059
00060
00061
00062
00063
00064
00065 if(verboseLevel > 0) {
00066 G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl
00067 << "Energy range: "
00068 << lowEnergyLimit / keV << " keV - "
00069 << highEnergyLimit / GeV << " GeV"
00070 << G4endl;
00071 }
00072
00073 crossSectionHandler = new G4CrossSectionHandler();
00074 crossSectionHandler->Initialise(0,lowEnergyLimit,highEnergyLimit,400);
00075 }
00076
00077
00078
00079 G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel()
00080 {
00081 delete crossSectionHandler;
00082 }
00083
00084
00085
00086
00087
00088 void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
00089 const G4DataVector& cuts)
00090 {
00091 if (verboseLevel > 3)
00092 G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl;
00093
00094 if (crossSectionHandler)
00095 {
00096 crossSectionHandler->Clear();
00097 delete crossSectionHandler;
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 if (HighEnergyLimit() > highEnergyLimit)
00111 {
00112 G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " <<
00113 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
00114
00115
00116 }
00117
00118
00119
00120 crossSectionHandler = new G4CrossSectionHandler;
00121 crossSectionHandler->Clear();
00122 G4String crossSectionFile = "pair/pp-cs-";
00123 crossSectionHandler->LoadData(crossSectionFile);
00124
00125
00126 if (verboseLevel > 2) {
00127 G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model"
00128 << G4endl;
00129 }
00130 InitialiseElementSelectors(particle,cuts);
00131
00132 if(verboseLevel > 0) {
00133 G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl
00134 << "Energy range: "
00135 << LowEnergyLimit() / keV << " keV - "
00136 << HighEnergyLimit() / GeV << " GeV"
00137 << G4endl;
00138 }
00139
00140
00141 if(!isInitialised) {
00142 isInitialised = true;
00143 fParticleChange = GetParticleChangeForGamma();
00144 }
00145 }
00146
00147
00148
00149 G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
00150 const G4ParticleDefinition*,
00151 G4double GammaEnergy,
00152 G4double Z, G4double,
00153 G4double, G4double)
00154 {
00155 if (verboseLevel > 3) {
00156 G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()"
00157 << G4endl;
00158 }
00159 if(Z < 0.9 || GammaEnergy <= lowEnergyLimit) { return 0.0; }
00160 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00161 return cs;
00162 }
00163
00164
00165
00166 void
00167 G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00168 const G4MaterialCutsCouple* couple,
00169 const G4DynamicParticle* aDynamicGamma,
00170 G4double,
00171 G4double)
00172 {
00173
00174
00175
00176
00177
00178
00179 if (verboseLevel > 3)
00180 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
00181
00182 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00183
00184
00185 if(photonEnergy <= lowEnergyLimit)
00186 {
00187 fParticleChange->ProposeTrackStatus(fStopAndKill);
00188 fParticleChange->SetProposedKineticEnergy(0.);
00189 return;
00190 }
00191
00192
00193 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
00194 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
00195
00196
00197
00198
00199 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
00200 {
00201 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
00202 }
00203 else
00204 {
00205 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
00206 {
00207 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
00208 }
00209 }
00210
00211
00212
00213
00214 G4double epsilon ;
00215 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00216
00217
00218
00219 if (photonEnergy < smallEnergy )
00220 {
00221 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00222 }
00223 else
00224 {
00225
00226
00227
00228
00229
00230
00231 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00232 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00233
00234
00235
00236
00237
00238
00239
00240
00241 G4IonisParamElm* ionisation = element->GetIonisation();
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 G4double fZ = 8. * (ionisation->GetlogZ3());
00254 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00255
00256
00257 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00258 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
00259 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00260
00261
00262 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
00263 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00264 G4double epsilonRange = 0.5 - epsilonMin ;
00265
00266
00267 G4double screen;
00268 G4double gReject ;
00269
00270 G4double f10 = ScreenFunction1(screenMin) - fZ;
00271 G4double f20 = ScreenFunction2(screenMin) - fZ;
00272 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00273 G4double normF2 = std::max(1.5 * f20,0.);
00274
00275 do {
00276 if (normF1 / (normF1 + normF2) > G4UniformRand() )
00277 {
00278 epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
00279 screen = screenFactor / (epsilon * (1. - epsilon));
00280 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00281 }
00282 else
00283 {
00284 epsilon = epsilonMin + epsilonRange * G4UniformRand();
00285 screen = screenFactor / (epsilon * (1 - epsilon));
00286 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00287
00288
00289 }
00290 } while ( gReject < G4UniformRand() );
00291
00292 }
00293
00294
00295
00296 G4double electronTotEnergy;
00297 G4double positronTotEnergy;
00298
00299
00300 if (G4int(2*G4UniformRand()))
00301 {
00302 electronTotEnergy = (1. - epsilon) * photonEnergy;
00303 positronTotEnergy = epsilon * photonEnergy;
00304 }
00305 else
00306 {
00307 positronTotEnergy = (1. - epsilon) * photonEnergy;
00308 electronTotEnergy = epsilon * photonEnergy;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 G4double Ene = electronTotEnergy/electron_mass_c2;
00331
00332 G4double cosTheta = 0.;
00333 G4double sinTheta = 0.;
00334
00335 SetTheta(&cosTheta,&sinTheta,Ene);
00336
00337
00338
00339
00340 G4double phi,psi=0.;
00341
00342
00343
00344
00345
00346 phi = SetPhi(photonEnergy);
00347 psi = SetPsi(photonEnergy,phi);
00348
00349
00350
00351
00352
00353
00354
00355 Psi = psi;
00356 Phi = phi;
00357
00358
00359
00360 G4double phie = psi;
00361
00362 G4double dirX = sinTheta*cos(phie);
00363 G4double dirY = sinTheta*sin(phie);
00364 G4double dirZ = cosTheta;
00365 G4ThreeVector electronDirection(dirX,dirY,dirZ);
00366
00367
00368
00369
00370
00371
00372 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00373
00374 SystemOfRefChange(gammaDirection0,electronDirection,
00375 gammaPolarization0);
00376
00377 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00378 electronDirection,
00379 electronKineEnergy);
00380
00381
00382
00383 Ene = positronTotEnergy/electron_mass_c2;
00384
00385 cosTheta = 0.;
00386 sinTheta = 0.;
00387
00388 SetTheta(&cosTheta,&sinTheta,Ene);
00389 G4double phip = phie+phi;
00390
00391 dirX = sinTheta*cos(phip);
00392 dirY = sinTheta*sin(phip);
00393 dirZ = cosTheta;
00394 G4ThreeVector positronDirection(dirX,dirY,dirZ);
00395
00396 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00397 SystemOfRefChange(gammaDirection0,positronDirection,
00398 gammaPolarization0);
00399
00400
00401 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00402 positronDirection, positronKineEnergy);
00403
00404
00405 fvect->push_back(particle1);
00406 fvect->push_back(particle2);
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
00420 fParticleChange->SetProposedKineticEnergy(0.);
00421 fParticleChange->ProposeTrackStatus(fStopAndKill);
00422
00423 }
00424
00425
00426
00427 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable)
00428 {
00429
00430
00431 G4double value;
00432
00433 if (screenVariable > 1.)
00434 value = 42.24 - 8.368 * log(screenVariable + 0.952);
00435 else
00436 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00437
00438 return value;
00439 }
00440
00441
00442
00443 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable)
00444 {
00445
00446
00447 G4double value;
00448
00449 if (screenVariable > 1.)
00450 value = 42.24 - 8.368 * log(screenVariable + 0.952);
00451 else
00452 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00453
00454 return value;
00455 }
00456
00457
00458 void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy)
00459 {
00460
00461
00462
00463
00464 G4double Momentum = sqrt(Energy*Energy -1);
00465 G4double Rand = G4UniformRand();
00466
00467 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
00468 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
00469 }
00470
00471
00472
00473 G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy)
00474 {
00475
00476
00477 G4double value = 0.;
00478 G4double Ene = Energy/MeV;
00479
00480 G4double pl[4];
00481
00482
00483 G4double pt[2];
00484 G4double xi = 0;
00485 G4double xe = 0.;
00486 G4double n1=0.;
00487 G4double n2=0.;
00488
00489
00490 if (Ene>=50.)
00491 {
00492 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
00493 const G4double aw = 0.0151, bw = 10.7, cw = -410.;
00494
00495 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
00496
00497 pl[0] = Fln(ay0,by0,Ene);
00498 pl[1] = aa0 + ba0*(Ene);
00499 pl[2] = Poli(aw,bw,cw,Ene);
00500 pl[3] = Poli(axc,bxc,cxc,Ene);
00501
00502 const G4double abf = 3.1216, bbf = 2.68;
00503 pt[0] = -1.4;
00504 pt[1] = abf + bbf/Ene;
00505
00506
00507
00508
00509
00510 xi = 3.0;
00511 xe = Encu(pl,pt,xi);
00512
00513 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
00514 n2 = Finttan(pt,xe) - Finttan(pt,0.);
00515 }
00516 else
00517 {
00518 const G4double ay0=0.144, by0=0.11;
00519 const G4double aa0=2.7, ba0 = 2.74;
00520 const G4double aw = 0.21, bw = 10.8, cw = -58.;
00521 const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
00522
00523 pl[0] = Fln(ay0, by0, Ene);
00524 pl[1] = Fln(aa0, ba0, Ene);
00525 pl[2] = Poli(aw,bw,cw,Ene);
00526 pl[3] = Poli(axc,bxc,cxc,Ene);
00527
00528
00529
00530 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
00531
00532 }
00533
00534
00535 G4double n=0.;
00536 n = n1+n2;
00537
00538 G4double c1 = 0.;
00539 c1 = Glor(pl, xe);
00540
00541
00542
00543
00544
00545
00546 G4double r1,r2,r3;
00547 G4double xco=0.;
00548
00549 if (Ene>=50.)
00550 {
00551 r1= G4UniformRand();
00552 if( r1>=n2/n)
00553 {
00554 do
00555 {
00556 r2 = G4UniformRand();
00557 value = Finvlor(pl,xe,r2);
00558 xco = Glor(pl,value)/c1;
00559 r3 = G4UniformRand();
00560 } while(r3>=xco);
00561 }
00562 else
00563 {
00564 value = Finvtan(pt,n,r1);
00565 }
00566 }
00567 else
00568 {
00569 do
00570 {
00571 r2 = G4UniformRand();
00572 value = Finvlor(pl,xe,r2);
00573 xco = Glor(pl,value)/c1;
00574 r3 = G4UniformRand();
00575 } while(r3>=xco);
00576 }
00577
00578
00579 return value;
00580 }
00581 G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double PhiLocal)
00582 {
00583
00584 G4double value = 0.;
00585 G4double Ene = Energy/MeV;
00586
00587 G4double p0l[4];
00588 G4double ppml[4];
00589 G4double p0t[2];
00590 G4double ppmt[2];
00591
00592 G4double xi = 0.;
00593 G4double xe0 = 0.;
00594 G4double xepm = 0.;
00595
00596 if (Ene>=50.)
00597 {
00598 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
00599 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
00600 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
00601 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
00602 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
00603 const G4double axcp = 2.8E-3,bxcp = -3.133;
00604 const G4double abf0 = 3.1213, bbf0 = 2.61;
00605 const G4double abfpm = 3.1231, bbfpm = 2.84;
00606
00607 p0l[0] = Fln(ay00, by00, Ene);
00608 p0l[1] = Fln(aa00, ba00, Ene);
00609 p0l[2] = Poli(aw0, bw0, cw0, Ene);
00610 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
00611
00612 ppml[0] = Fln(ay0p, by0p, Ene);
00613 ppml[1] = aap + bap*(Ene);
00614 ppml[2] = Poli(awp, bwp, cwp, Ene);
00615 ppml[3] = Fln(axcp,bxcp,Ene);
00616
00617 p0t[0] = -0.81;
00618 p0t[1] = abf0 + bbf0/Ene;
00619 ppmt[0] = -0.6;
00620 ppmt[1] = abfpm + bbfpm/Ene;
00621
00622
00623
00624
00625 xi = 3.0;
00626 xe0 = Encu(p0l, p0t, xi);
00627
00628 xepm = Encu(ppml, ppmt, xi);
00629
00630 }
00631 else
00632 {
00633 const G4double ay00 = 2.82, by00 = 6.35;
00634 const G4double aa00 = -1.75, ba00 = 0.25;
00635
00636 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
00637 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
00638 const G4double ay0p = 1.56, by0p = 3.6;
00639 const G4double aap = 0.86, bap = 8.3E-3;
00640 const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
00641 const G4double xcp = 3.1486;
00642
00643 p0l[0] = Fln(ay00, by00, Ene);
00644 p0l[1] = aa00+pow(Ene, ba00);
00645 p0l[2] = Poli(aw0, bw0, cw0, Ene);
00646 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
00647 ppml[0] = Fln(ay0p, by0p, Ene);
00648 ppml[1] = aap + bap*(Ene);
00649 ppml[2] = Poli(awp, bwp, cwp, Ene);
00650 ppml[3] = xcp;
00651
00652 }
00653
00654 G4double a,b=0.;
00655
00656 if (Ene>=50.)
00657 {
00658 if (PhiLocal>xepm)
00659 {
00660 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
00661 }
00662 else
00663 {
00664 b = Ftan(ppmt,PhiLocal);
00665 }
00666 if (PhiLocal>xe0)
00667 {
00668 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
00669 }
00670 else
00671 {
00672 a = Ftan(p0t,PhiLocal);
00673 }
00674 }
00675 else
00676 {
00677 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
00678 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
00679 }
00680 G4double nr =0.;
00681
00682 if (b>a)
00683 {
00684 nr = 1./b;
00685 }
00686 else
00687 {
00688 nr = 1./a;
00689 }
00690
00691 G4double r1,r2=0.;
00692 G4double r3 =-1.;
00693 do
00694 {
00695 r1 = G4UniformRand();
00696 r2 = G4UniformRand();
00697 value = r2*pi;
00698 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
00699 }while(r1>r3);
00700
00701 return value;
00702 }
00703
00704
00705 G4double G4LivermorePolarizedGammaConversionModel::Poli
00706 (G4double a, G4double b, G4double c, G4double x)
00707 {
00708 G4double value=0.;
00709 if(x>0.)
00710 {
00711 value =(a + b/x + c/(x*x*x));
00712 }
00713 else
00714 {
00715
00716 }
00717 return value;
00718 }
00719 G4double G4LivermorePolarizedGammaConversionModel::Fln
00720 (G4double a, G4double b, G4double x)
00721 {
00722 G4double value=0.;
00723 if(x>0.)
00724 {
00725 value =(a*log(x)-b);
00726 }
00727 else
00728 {
00729
00730 }
00731 return value;
00732 }
00733
00734
00735 G4double G4LivermorePolarizedGammaConversionModel::Encu
00736 (G4double* p_p1, G4double* p_p2, G4double x0)
00737 {
00738 G4int i=0;
00739 G4double fx = 1.;
00740 G4double x = x0;
00741 const G4double xmax = 3.0;
00742
00743 for(i=0; i<100; ++i)
00744 {
00745 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
00746 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
00747 x -= fx;
00748 if(x > xmax) { return xmax; }
00749
00750
00751
00752
00753 if(std::fabs(fx) <= x*1.0e-6) { break; }
00754 }
00755
00756 if(x < 0.0) { x = 0.0; }
00757 return x;
00758 }
00759
00760
00761 G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x)
00762 {
00763 G4double value =0.;
00764
00765
00766 G4double w = p_p1[2];
00767 G4double xc = p_p1[3];
00768
00769 value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
00770 return value;
00771 }
00772
00773
00774 G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x)
00775 {
00776 G4double value =0.;
00777 G4double y0 = p_p1[0];
00778 G4double A = p_p1[1];
00779 G4double w = p_p1[2];
00780 G4double xc = p_p1[3];
00781
00782 value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
00783 return value;
00784 }
00785
00786
00787 G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x)
00788 {
00789 G4double value =0.;
00790
00791 G4double A = p_p1[1];
00792 G4double w = p_p1[2];
00793 G4double xc = p_p1[3];
00794
00795 value = (-16.*A*w*(x-xc))/
00796 (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
00797 return value;
00798 }
00799
00800
00801 G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x)
00802 {
00803 G4double value =0.;
00804 G4double y0 = p_p1[0];
00805 G4double A = p_p1[1];
00806 G4double w = p_p1[2];
00807 G4double xc = p_p1[3];
00808
00809 value = y0*x + A*atan( 2*(x-xc)/w) / pi;
00810 return value;
00811 }
00812
00813
00814 G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r)
00815 {
00816 G4double value = 0.;
00817 G4double nor = 0.;
00818
00819
00820 G4double w = p_p1[2];
00821 G4double xc = p_p1[3];
00822
00823 nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
00824 value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
00825
00826 return value;
00827 }
00828
00829
00830 G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x)
00831 {
00832 G4double value =0.;
00833 G4double a = p_p1[0];
00834 G4double b = p_p1[1];
00835
00836 value = a /(x-b);
00837 return value;
00838 }
00839
00840
00841 G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x)
00842 {
00843 G4double value =0.;
00844 G4double a = p_p1[0];
00845 G4double b = p_p1[1];
00846
00847 value = -1.*a / ((x-b)*(x-b));
00848 return value;
00849 }
00850
00851
00852 G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x)
00853 {
00854 G4double value =0.;
00855 G4double a = p_p1[0];
00856 G4double b = p_p1[1];
00857
00858
00859 value = a*log(b-x);
00860 return value;
00861 }
00862
00863
00864 G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r)
00865 {
00866 G4double value =0.;
00867 G4double a = p_p1[0];
00868 G4double b = p_p1[1];
00869
00870 value = b*(1-exp(r*cnor/a));
00871
00872 return value;
00873 }
00874
00875
00876
00877
00878
00879
00880 G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a)
00881 {
00882 G4double dx = a.x();
00883 G4double dy = a.y();
00884 G4double dz = a.z();
00885 G4double x = dx < 0.0 ? -dx : dx;
00886 G4double y = dy < 0.0 ? -dy : dy;
00887 G4double z = dz < 0.0 ? -dz : dz;
00888 if (x < y) {
00889 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
00890 }else{
00891 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
00892 }
00893 }
00894
00895
00896
00897 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0)
00898 {
00899 G4ThreeVector d0 = direction0.unit();
00900 G4ThreeVector a1 = SetPerpendicularVector(d0);
00901 G4ThreeVector a0 = a1.unit();
00902
00903 G4double rand1 = G4UniformRand();
00904
00905 G4double angle = twopi*rand1;
00906 G4ThreeVector b0 = d0.cross(a0);
00907
00908 G4ThreeVector c;
00909
00910 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
00911 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
00912 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
00913
00914 G4ThreeVector c0 = c.unit();
00915
00916 return c0;
00917
00918 }
00919
00920
00921
00922 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization
00923 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
00924 {
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
00937 }
00938
00939
00940
00941
00942 void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange
00943 (G4ThreeVector& direction0,G4ThreeVector& direction1,
00944 G4ThreeVector& polarization0)
00945 {
00946
00947
00948
00949 G4ThreeVector Axis_Z0 = direction0.unit();
00950 G4ThreeVector Axis_X0 = polarization0.unit();
00951 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit();
00952
00953 G4double direction_x = direction1.getX();
00954 G4double direction_y = direction1.getY();
00955 G4double direction_z = direction1.getZ();
00956
00957 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
00958
00959 }
00960
00961
00962
00963