| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Authors: G.Depaola & F.Longo
28 //
29 //
30
31 #include "G4LivermorePolarizedGammaConversionM
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Electron.hh"
35 #include "G4Positron.hh"
36 #include "G4ParticleChangeForGamma.hh"
37 #include "G4Log.hh"
38 #include "G4AutoLock.hh"
39 #include "G4Exp.hh"
40 #include "G4ProductionCutsTable.hh"
41
42 //....oooOO0OOooo........oooOO0OOooo........oo
43
44 using namespace std;
45 namespace { G4Mutex LivermorePolarizedGammaCon
46
47 G4PhysicsFreeVector* G4LivermorePolarizedGamma
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50
51 G4LivermorePolarizedGammaConversionModel::G4Li
52 const G4ParticleDefinition*, const G4String
53 :G4VEmModel(nam), smallEnergy(2.*MeV), isIni
54 {
55 fParticleChange = nullptr;
56 lowEnergyLimit = 2*electron_mass_c2;
57
58 Phi=0.;
59 Psi=0.;
60
61 verboseLevel= 0;
62 // Verbosity scale:
63 // 0 = nothing
64 // 1 = calculation of cross sections, file o
65 // 2 = entering in methods
66
67 if(verboseLevel > 0) {
68 G4cout << "Livermore Polarized GammaConver
69 << G4endl;
70 }
71
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 G4LivermorePolarizedGammaConversionModel::~G4L
77 {
78 if(IsMaster()) {
79 for(G4int i=0; i<maxZ; ++i) {
80 if(data[i]) {
81 delete data[i];
82 data[i] = nullptr;
83 }
84 }
85 }
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
89
90 void G4LivermorePolarizedGammaConversionModel:
91 const G
92 {
93 if (verboseLevel > 1)
94 {
95 G4cout << "Calling1 G4LivermorePolarized
96 << G4endl
97 << "Energy range: "
98 << LowEnergyLimit() / MeV << " MeV - "
99 << HighEnergyLimit() / GeV << " G
100 << G4endl;
101 }
102
103 if(IsMaster())
104 {
105 // Initialise element selector
106 InitialiseElementSelectors(particle, cut
107
108 // Access to elements
109 const char* path = G4FindDataDir("G4LEDA
110
111 G4ProductionCutsTable* theCoupleTable =
112 G4ProductionCutsTable::GetProductionCutsTabl
113
114 G4int numOfCouples = (G4int)theCoupleTab
115
116 for(G4int i=0; i<numOfCouples; ++i)
117 {
118 const G4Material* material =
119 theCoupleTable->GetMaterialCutsCouple(i)
120 const G4ElementVector* theElementVector =
121 std::size_t nelm = material->GetNumberOfEl
122
123 for (std::size_t j=0; j<nelm; ++j)
124 {
125 G4int Z = (G4int)(*theElementVector)[j
126 if(Z < 1) { Z = 1; }
127 else if(Z > maxZ) { Z = maxZ; }
128 if(!data[Z]) { ReadData(Z, path); }
129 }
130 }
131 }
132 if(isInitialised) { return; }
133 fParticleChange = GetParticleChangeForGamma(
134 isInitialised = true;
135 }
136
137 //....oooOO0OOooo........oooOO0OOooo........oo
138
139 void G4LivermorePolarizedGammaConversionModel:
140 const G4ParticleDefinition*, G4VEmModel*
141 {
142 SetElementSelectors(masterModel->GetElementS
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oo
146
147 G4double G4LivermorePolarizedGammaConversionMo
148 const G4ParticleDefinition*, G4doub
149 {
150 return lowEnergyLimit;
151 }
152
153 //....oooOO0OOooo........oooOO0OOooo........oo
154
155 void G4LivermorePolarizedGammaConversionModel:
156 {
157 if (verboseLevel > 1)
158 {
159 G4cout << "Calling ReadData() of G4Liver
160 << G4endl;
161 }
162
163 if(data[Z]) { return; }
164
165 const char* datadir = path;
166
167 if(!datadir)
168 {
169 datadir = G4FindDataDir("G4LEDATA");
170 if(!datadir)
171 {
172 G4Exception("G4LivermorePolarizedGammaConv
173 "em0006",FatalException,
174 "Environment variable G4LEDATA not d
175 return;
176 }
177 }
178 //
179 data[Z] = new G4PhysicsFreeVector(0,/*spline
180 //
181 std::ostringstream ost;
182 ost << datadir << "/livermore/pair/pp-cs-" <
183 std::ifstream fin(ost.str().c_str());
184
185 if( !fin.is_open())
186 {
187 G4ExceptionDescription ed;
188 ed << "G4LivermorePolarizedGammaConversi
189 << "> is not opened!" << G4endl;
190 G4Exception("G4LivermorePolarizedGammaCo
191 "em0003",FatalException,
192 ed,"G4LEDATA version should be G4EMLOW6.
193 return;
194 }
195 else
196 {
197
198 if(verboseLevel > 3) { G4cout << "File "
199 << " is opened by G4LivermorePolar
200
201 data[Z]->Retrieve(fin, true);
202 }
203
204 // Activation of spline interpolation
205 data[Z]->FillSecondDerivatives();
206
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oo
210
211 G4double G4LivermorePolarizedGammaConversionMo
212 const G
213 G4double GammaEnergy,
214 G4double Z, G4double,
215 G4double, G4double)
216 {
217 if (verboseLevel > 1) {
218 G4cout << "G4LivermorePolarizedGammaConver
219 << G4endl;
220 }
221 if (GammaEnergy < lowEnergyLimit) { return 0
222
223 G4double xs = 0.0;
224
225 G4int intZ=G4int(Z);
226
227 if(intZ < 1 || intZ > maxZ) { return xs; }
228
229 G4PhysicsFreeVector* pv = data[intZ];
230
231 // if element was not initialised
232 // do initialisation safely for MT mode
233 if(!pv)
234 {
235 InitialiseForElement(0, intZ);
236 pv = data[intZ];
237 if(!pv) { return xs; }
238 }
239 // x-section is taken from the table
240 xs = pv->Value(GammaEnergy);
241
242 if(verboseLevel > 0)
243 {
244 G4int n = G4int(pv->GetVectorLength() -
245 G4cout << "****** DEBUG: tcs value for
246 << GammaEnergy/MeV << G4endl;
247 G4cout << " cs (Geant4 internal unit)
248 G4cout << " -> first cs value in EA
249 G4cout << " -> last cs value in EA
250 G4cout << "***************************
251 }
252
253 return xs;
254 }
255
256 //....oooOO0OOooo........oooOO0OOooo........oo
257
258 void
259 G4LivermorePolarizedGammaConversionModel::Samp
260 const G4MaterialCutsCouple*
261 const G4DynamicParticle* aDy
262 G4double,
263 G4double)
264 {
265
266 // Fluorescence generated according to:
267 // J. Stepanek ,"A program to determine the
268 // subshell ionisation by a particle or due
269 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
270 if (verboseLevel > 3)
271 G4cout << "Calling SampleSecondaries() of
272
273 G4double photonEnergy = aDynamicGamma->GetKi
274
275 if(photonEnergy <= lowEnergyLimit)
276 {
277 fParticleChange->ProposeTrackStatus(fSto
278 fParticleChange->SetProposedKineticEnerg
279 return;
280 }
281
282 G4ThreeVector gammaPolarization0 = aDynamicG
283 G4ThreeVector gammaDirection0 = aDynamicGamm
284
285 // Make sure that the polarization vector is
286 // gamma direction. If not
287 if(!(gammaPolarization0.isOrthogonal(gammaDi
288 { // only for testing now
289 gammaPolarization0 = GetRandomPolarizati
290 }
291 else
292 {
293 if ( gammaPolarization0.howOrthogonal(ga
294 {
295 gammaPolarization0 = GetPerpendicularPolar
296 }
297 }
298
299 // End of Protection
300
301 G4double epsilon ;
302 G4double epsilon0Local = electron_mass_c2 /
303
304 // Do it fast if photon energy < 2. MeV
305
306 if (photonEnergy < smallEnergy )
307 {
308 epsilon = epsilon0Local + (0.5 - epsilon
309 }
310 else
311 {
312 // Select randomly one element in the cu
313 const G4ParticleDefinition* particle =
314 const G4Element* element = SelectRandomA
315
316 if (element == nullptr)
317 {
318 G4cout << "G4LivermorePolarizedGamma
319 return;
320 }
321
322
323 G4IonisParamElm* ionisation = element->G
324 if (ionisation == nullptr)
325 {
326 G4cout << "G4LivermorePolarizedGamma
327 return;
328 }
329
330 // Extract Coulomb factor for this Eleme
331 G4double fZ = 8. * (ionisation->GetlogZ3
332 if (photonEnergy > 50. * MeV) fZ += 8. *
333
334 // Limits of the screening variable
335 G4double screenFactor = 136. * epsilon0L
336 G4double screenMax = G4Exp ((42.24 - fZ)
337 G4double screenMin = std::min(4.*screenF
338
339 // Limits of the energy sampling
340 G4double epsilon1 = 0.5 - 0.5 * sqrt(1.
341 G4double epsilonMin = std::max(epsilon0L
342 G4double epsilonRange = 0.5 - epsilonMin
343
344 // Sample the energy rate of the created
345 G4double screen;
346 G4double gReject ;
347
348 G4double f10 = ScreenFunction1(screenMin
349 G4double f20 = ScreenFunction2(screenMin
350 G4double normF1 = std::max(f10 * epsilon
351 G4double normF2 = std::max(1.5 * f20,0.)
352
353 do {
354 if (normF1 / (normF1 + normF2) > G4Uni
355 {
356 epsilon = 0.5 - epsilonRange * pow
357 screen = screenFactor / (epsilon *
358 gReject = (ScreenFunction1(screen)
359 }
360 else
361 {
362 epsilon = epsilonMin + epsilonRang
363 screen = screenFactor / (epsilon *
364 gReject = (ScreenFunction2(screen)
365 }
366 } while ( gReject < G4UniformRand() );
367 } // End of epsilon sampling
368
369 // Fix charges randomly
370 G4double electronTotEnergy;
371 G4double positronTotEnergy;
372
373 if (G4UniformRand() > 0.5)
374 {
375 electronTotEnergy = (1. - epsilon) * pho
376 positronTotEnergy = epsilon * photonEner
377 }
378 else
379 {
380 positronTotEnergy = (1. - epsilon) * pho
381 electronTotEnergy = epsilon * photonEner
382 }
383
384 // Scattered electron (positron) angles. ( Z
385 // Universal distribution suggested by L. Ur
386 // derived from Tsai distribution (Rev. Mod.
387 G4double Ene = electronTotEnergy/electron_ma
388
389 G4double cosTheta = 0.;
390 G4double sinTheta = 0.;
391
392 SetTheta(&cosTheta,&sinTheta,Ene);
393 G4double phi,psi=0.;
394
395 //corrected e+ e- angular angular distributi
396 phi = SetPhi(photonEnergy);
397 psi = SetPsi(photonEnergy,phi);
398 Psi = psi;
399 Phi = phi;
400
401 G4double phie, phip;
402 G4double choice, choice2;
403 choice = G4UniformRand();
404 choice2 = G4UniformRand();
405
406 if (choice2 <= 0.5)
407 {
408 // do nothing
409 // phi = phi;
410 }
411 else
412 {
413 phi = -phi;
414 }
415
416 if (choice <= 0.5)
417 {
418 phie = psi; //azimuthal angle for the el
419 phip = phie+phi; //azimuthal angle for t
420 }
421 else
422 {
423 // opzione 1 phie / phip equivalenti
424 phip = psi; //azimuthal angle for the po
425 phie = phip + phi; //azimuthal angle for
426 }
427
428
429 // Electron Kinematics
430 G4double dirX = sinTheta*cos(phie);
431 G4double dirY = sinTheta*sin(phie);
432 G4double dirZ = cosTheta;
433 G4ThreeVector electronDirection(dirX,dirY,di
434
435 // Kinematics of the created pair:
436 // the electron and positron are assumed to
437 // distribution with respect to the Z axis a
438
439 G4double electronKineEnergy = std::max(0.,el
440
441 SystemOfRefChange(gammaDirection0,electronDi
442 gammaPolarization0);
443
444 G4DynamicParticle* particle1 = new G4Dynamic
445 electronDirection,
446 electronKineEnergy);
447
448 // The e+ is always created (even with kinet
449 Ene = positronTotEnergy/electron_mass_c2; //
450
451 cosTheta = 0.;
452 sinTheta = 0.;
453
454 SetTheta(&cosTheta,&sinTheta,Ene);
455
456 // Positron Kinematics
457 dirX = sinTheta*cos(phip);
458 dirY = sinTheta*sin(phip);
459 dirZ = cosTheta;
460 G4ThreeVector positronDirection(dirX,dirY,di
461
462 G4double positronKineEnergy = std::max(0.,po
463 SystemOfRefChange(gammaDirection0,positronDi
464 gammaPolarization0);
465
466 // Create G4DynamicParticle object for the p
467 G4DynamicParticle* particle2 = new G4Dynamic
468
469 fvect->push_back(particle1);
470 fvect->push_back(particle2);
471
472 // Kill the incident photon
473 fParticleChange->SetProposedKineticEnergy(0.
474 fParticleChange->ProposeTrackStatus(fStopAnd
475 }
476
477 //....oooOO0OOooo........oooOO0OOooo........oo
478
479 G4double G4LivermorePolarizedGammaConversionMo
480 {
481 // Compute the value of the screening functi
482 G4double value;
483 if (screenVariable > 1.)
484 value = 42.24 - 8.368 * log(screenVariable
485 else
486 value = 42.392 - screenVariable * (7.796 -
487
488 return value;
489 }
490
491
492
493 G4double G4LivermorePolarizedGammaConversionMo
494 {
495 // Compute the value of the screening functi
496 G4double value;
497
498 if (screenVariable > 1.)
499 value = 42.24 - 8.368 * log(screenVariable
500 else
501 value = 41.405 - screenVariable * (5.828 -
502
503 return value;
504 }
505
506
507 void G4LivermorePolarizedGammaConversionModel:
508 {
509 // to avoid computational errors since Theta
510 // Energy in Normalized Units (!)
511
512 G4double Momentum = sqrt(Energy*Energy -1);
513 G4double Rand = G4UniformRand();
514
515 *p_cosTheta = (Energy*((2*Rand)- 1) + Moment
516 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momen
517 }
518
519
520
521 G4double G4LivermorePolarizedGammaConversionMo
522 {
523 G4double value = 0.;
524 G4double Ene = Energy/MeV;
525
526 G4double pl[4];
527 G4double pt[2];
528 G4double xi = 0;
529 G4double xe = 0.;
530 G4double n1=0.;
531 G4double n2=0.;
532
533 if (Ene>=50.)
534 {
535 const G4double ay0=5.6, by0=18.6, aa0=2.
536 const G4double aw = 0.0151, bw = 10.7, c
537
538 const G4double axc = 3.1455, bxc = -1.11
539
540 pl[0] = Fln(ay0,by0,Ene);
541 pl[1] = aa0 + ba0*(Ene);
542 pl[2] = Poli(aw,bw,cw,Ene);
543 pl[3] = Poli(axc,bxc,cxc,Ene);
544
545 const G4double abf = 3.1216, bbf = 2.68;
546 pt[0] = -1.4;
547 pt[1] = abf + bbf/Ene;
548
549 xi = 3.0;
550 xe = Encu(pl,pt,xi);
551 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
552 n2 = Finttan(pt,xe) - Finttan(pt,0.);
553 }
554 else
555 {
556 const G4double ay0=0.144, by0=0.11;
557 const G4double aa0=2.7, ba0 = 2.74;
558 const G4double aw = 0.21, bw = 10.8, cw
559 const G4double axc = 3.17, bxc = -0.87,
560
561 pl[0] = Fln(ay0, by0, Ene);
562 pl[1] = Fln(aa0, ba0, Ene);
563 pl[2] = Poli(aw,bw,cw,Ene);
564 pl[3] = Poli(axc,bxc,cxc,Ene);
565
566 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
567 }
568
569
570 G4double n=0.;
571 n = n1+n2;
572
573 G4double c1 = 0.;
574 c1 = Glor(pl, xe);
575
576 G4double r1,r2,r3;
577 G4double xco=0.;
578
579 if (Ene>=50.)
580 {
581 r1= G4UniformRand();
582 if( r1>=n2/n)
583 {
584 do
585 {
586 r2 = G4UniformRand();
587 value = Finvlor(pl,xe,r2);
588 xco = Glor(pl,value)/c1;
589 r3 = G4UniformRand();
590 } while(r3>=xco);
591 }
592 else
593 {
594 value = Finvtan(pt,n,r1);
595 }
596 }
597 else
598 {
599 do
600 {
601 r2 = G4UniformRand();
602 value = Finvlor(pl,xe,r2);
603 xco = Glor(pl,value)/c1;
604 r3 = G4UniformRand();
605 } while(r3>=xco);
606 }
607 return value;
608 }
609
610 //....oooOO0OOooo........oooOO0OOooo........oo
611
612 G4double G4LivermorePolarizedGammaConversionMo
613 {
614 G4double value = 0.;
615 G4double Ene = Energy/MeV;
616
617 G4double p0l[4];
618 G4double ppml[4];
619 G4double p0t[2];
620 G4double ppmt[2];
621
622 G4double xi = 0.;
623 G4double xe0 = 0.;
624 G4double xepm = 0.;
625
626 if (Ene>=50.)
627 {
628 const G4double ay00 = 3.4, by00 = 9.8, a
629 const G4double aw0 = 0.014, bw0 = 9.7, c
630 const G4double axc0 = 3.1423, bxc0 = -2.
631 const G4double ay0p = 1.53, by0p = 3.2,
632 const G4double awp = 6.9E-3, bwp = 12.6,
633 const G4double axcp = 2.8E-3,bxcp = -3.1
634 const G4double abf0 = 3.1213, bbf0 = 2.6
635 const G4double abfpm = 3.1231, bbfpm = 2
636
637 p0l[0] = Fln(ay00, by00, Ene);
638 p0l[1] = Fln(aa00, ba00, Ene);
639 p0l[2] = Poli(aw0, bw0, cw0, Ene);
640 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
641
642 ppml[0] = Fln(ay0p, by0p, Ene);
643 ppml[1] = aap + bap*(Ene);
644 ppml[2] = Poli(awp, bwp, cwp, Ene);
645 ppml[3] = Fln(axcp,bxcp,Ene);
646
647 p0t[0] = -0.81;
648 p0t[1] = abf0 + bbf0/Ene;
649 ppmt[0] = -0.6;
650 ppmt[1] = abfpm + bbfpm/Ene;
651
652 xi = 3.0;
653 xe0 = Encu(p0l, p0t, xi);
654 xepm = Encu(ppml, ppmt, xi);
655 }
656 else
657 {
658 const G4double ay00 = 2.82, by00 = 6.35;
659 const G4double aa00 = -1.75, ba00 = 0.25
660
661 const G4double aw0 = 0.028, bw0 = 5., cw
662 const G4double axc0 = 3.14213, bxc0 = -2
663 const G4double ay0p = 1.56, by0p = 3.6;
664 const G4double aap = 0.86, bap = 8.3E-3;
665 const G4double awp = 0.022, bwp = 7.4, c
666 const G4double xcp = 3.1486;
667
668 p0l[0] = Fln(ay00, by00, Ene);
669 p0l[1] = aa00+pow(Ene, ba00);
670 p0l[2] = Poli(aw0, bw0, cw0, Ene);
671 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
672 ppml[0] = Fln(ay0p, by0p, Ene);
673 ppml[1] = aap + bap*(Ene);
674 ppml[2] = Poli(awp, bwp, cwp, Ene);
675 ppml[3] = xcp;
676 }
677
678 G4double a,b=0.;
679
680 if (Ene>=50.)
681 {
682 if (PhiLocal>xepm)
683 {
684 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(
685 }
686 else
687 {
688 b = Ftan(ppmt,PhiLocal);
689 }
690 if (PhiLocal>xe0)
691 {
692 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l
693 }
694 else
695 {
696 a = Ftan(p0t,PhiLocal);
697 }
698 }
699 else
700 {
701 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml
702 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi
703 }
704 G4double nr =0.;
705
706 if (b>a)
707 {
708 nr = 1./b;
709 }
710 else
711 {
712 nr = 1./a;
713 }
714
715 G4double r1,r2=0.;
716 G4double r3 =-1.;
717 do
718 {
719 r1 = G4UniformRand();
720 r2 = G4UniformRand();
721 //value = r2*pi;
722 value = 2.*r2*pi;
723 r3 = nr*(a*cos(value)*cos(value) + b*sin
724 }while(r1>r3);
725
726 return value;
727 }
728
729 //....oooOO0OOooo........oooOO0OOooo........oo
730
731 G4double G4LivermorePolarizedGammaConversionMo
732 (G4double a, G4double b, G4double c, G4double
733 {
734 G4double value=0.;
735 if(x>0.)
736 {
737 value =(a + b/x + c/(x*x*x));
738 }
739 else
740 {
741 //G4cout << "ERROR in Poli! " << G4endl;
742 }
743 return value;
744 }
745
746 //....oooOO0OOooo........oooOO0OOooo........oo
747
748 G4double G4LivermorePolarizedGammaConversionMo
749 (G4double a, G4double b, G4double x)
750 {
751 G4double value=0.;
752 if(x>0.)
753 {
754 value =(a*log(x)-b);
755 }
756 else
757 {
758 //G4cout << "ERROR in Fln! " << G4endl;
759 }
760 return value;
761 }
762
763 //....oooOO0OOooo........oooOO0OOooo........oo
764
765 G4double G4LivermorePolarizedGammaConversionMo
766 (G4double* p_p1, G4double* p_p2, G4double x0)
767 {
768 G4int i=0;
769 G4double fx = 1.;
770 G4double x = x0;
771 const G4double xmax = 3.0;
772
773 for(i=0; i<100; ++i)
774 {
775 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p
776 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
777 x -= fx;
778 if(x > xmax) { return xmax; }
779 if(std::fabs(fx) <= x*1.0e-6) { break; }
780 }
781
782 if(x < 0.0) { x = 0.0; }
783 return x;
784 }
785
786 //....oooOO0OOooo........oooOO0OOooo........oo
787
788 G4double G4LivermorePolarizedGammaConversionMo
789 {
790 G4double value =0.;
791 G4double w = p_p1[2];
792 G4double xc = p_p1[3];
793
794 value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
795 return value;
796 }
797
798 //....oooOO0OOooo........oooOO0OOooo........oo
799
800 G4double G4LivermorePolarizedGammaConversionMo
801 {
802 G4double value =0.;
803 G4double y0 = p_p1[0];
804 G4double A = p_p1[1];
805 G4double w = p_p1[2];
806 G4double xc = p_p1[3];
807
808 value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) +
809 return value;
810 }
811
812 //....oooOO0OOooo........oooOO0OOooo........oo
813
814 G4double G4LivermorePolarizedGammaConversionMo
815 {
816 G4double value =0.;
817 G4double A = p_p1[1];
818 G4double w = p_p1[2];
819 G4double xc = p_p1[3];
820
821 value = (-16.*A*w*(x-xc))/
822 (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*
823 return value;
824 }
825
826 //....oooOO0OOooo........oooOO0OOooo........oo
827
828 G4double G4LivermorePolarizedGammaConversionMo
829 {
830 G4double value =0.;
831 G4double y0 = p_p1[0];
832 G4double A = p_p1[1];
833 G4double w = p_p1[2];
834 G4double xc = p_p1[3];
835
836 value = y0*x + A*atan( 2*(x-xc)/w) / pi;
837 return value;
838 }
839
840
841 G4double G4LivermorePolarizedGammaConversionMo
842 {
843 G4double value = 0.;
844 G4double nor = 0.;
845 G4double w = p_p1[2];
846 G4double xc = p_p1[3];
847
848 nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.
849 value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(
850
851 return value;
852 }
853
854 //....oooOO0OOooo........oooOO0OOooo........oo
855
856 G4double G4LivermorePolarizedGammaConversionMo
857 {
858 G4double value =0.;
859 G4double a = p_p1[0];
860 G4double b = p_p1[1];
861
862 value = a /(x-b);
863 return value;
864 }
865
866 //....oooOO0OOooo........oooOO0OOooo........oo
867
868 G4double G4LivermorePolarizedGammaConversionMo
869 {
870 G4double value =0.;
871 G4double a = p_p1[0];
872 G4double b = p_p1[1];
873
874 value = -1.*a / ((x-b)*(x-b));
875 return value;
876 }
877
878 //....oooOO0OOooo........oooOO0OOooo........oo
879
880 G4double G4LivermorePolarizedGammaConversionMo
881 {
882 G4double value =0.;
883 G4double a = p_p1[0];
884 G4double b = p_p1[1];
885
886 value = a*log(b-x);
887 return value;
888 }
889
890 //....oooOO0OOooo........oooOO0OOooo........oo
891
892 G4double G4LivermorePolarizedGammaConversionMo
893 {
894 G4double value =0.;
895 G4double a = p_p1[0];
896 G4double b = p_p1[1];
897
898 value = b*(1-G4Exp(r*cnor/a));
899
900 return value;
901 }
902
903 //....oooOO0OOooo........oooOO0OOooo........oo
904
905 G4ThreeVector G4LivermorePolarizedGammaConvers
906 {
907 G4double dx = a.x();
908 G4double dy = a.y();
909 G4double dz = a.z();
910 G4double x = dx < 0.0 ? -dx : dx;
911 G4double y = dy < 0.0 ? -dy : dy;
912 G4double z = dz < 0.0 ? -dz : dz;
913 if (x < y) {
914 return x < z ? G4ThreeVector(-dy,dx,0) : G
915 }else{
916 return y < z ? G4ThreeVector(dz,0,-dx) : G
917 }
918 }
919
920 //....oooOO0OOooo........oooOO0OOooo........oo
921
922 G4ThreeVector G4LivermorePolarizedGammaConvers
923 {
924 G4ThreeVector d0 = direction0.unit();
925 G4ThreeVector a1 = SetPerpendicularVector(d0
926 G4ThreeVector a0 = a1.unit(); // unit vector
927
928 G4double rand1 = G4UniformRand();
929
930 G4double angle = twopi*rand1; // random pola
931 G4ThreeVector b0 = d0.cross(a0); // cross pr
932
933 G4ThreeVector c;
934
935 c.setX(std::cos(angle)*(a0.x())+std::sin(ang
936 c.setY(std::cos(angle)*(a0.y())+std::sin(ang
937 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang
938
939 G4ThreeVector c0 = c.unit();
940
941 return c0;
942 }
943
944 //....oooOO0OOooo........oooOO0OOooo........oo
945
946 G4ThreeVector G4LivermorePolarizedGammaConvers
947 (const G4ThreeVector& gammaDirection, const G4
948 {
949 //
950 // The polarization of a photon is always pe
951 // Therefore this function removes those vec
952 // points in direction of gammaDirection
953 //
954 // Mathematically we search the projection o
955 // plains normal vector.
956 // The basic equation can be found in each g
957 // p = a - (a o n)/(n o n)*n
958
959 return gammaPolarization - gammaPolarization
960 }
961
962 //....oooOO0OOooo........oooOO0OOooo........oo
963
964 void G4LivermorePolarizedGammaConversionModel:
965 (G4ThreeVector& direction0,G4ThreeVector&
966 G4ThreeVector& polarization0)
967 {
968 // direction0 is the original photon directi
969 // polarization0 is the original photon pola
970 // need to specify y axis in the real refere
971 G4ThreeVector Axis_Z0 = direction0.unit();
972 G4ThreeVector Axis_X0 = polarization0.unit()
973 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_
974
975 G4double direction_x = direction1.getX();
976 G4double direction_y = direction1.getY();
977 G4double direction_z = direction1.getZ();
978
979 direction1 = (direction_x*Axis_X0 + directio
980 }
981
982 //....oooOO0OOooo........oooOO0OOooo........oo
983
984 void G4LivermorePolarizedGammaConversionModel:
985 const G4ParticleDefiniti
986 G4int Z)
987 {
988 G4AutoLock l(&LivermorePolarizedGammaConvers
989 if(!data[Z]) { ReadData(Z); }
990 l.unlock();
991 }
992