| 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 // |
28 // | G4LowEPPolarizedComptonModel-- Geant
29 // | polarised low energy Compton scat
30 // | J. M. C. Brown, Monash Universit
31 // |
32 // |
33 // *******************************************
34 // |
35 // | The following is a Geant4 class to simula
36 // | bound electron Compton scattering. Genera
37 // | based on G4LowEnergyCompton.cc and
38 // | G4LivermorePolarizedComptonModel.cc.
39 // | Algorithms for photon energy, and ejected
40 // | direction taken from:
41 // |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill
43 // | "A low energy bound atomic electron Compt
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014
45 // |
46 // | The author acknowledges the work of the G
47 // | in developing the following algorithms th
48 // | or adapeted for the present software:
49 // |
50 // | # sampling of photon scattering angle,
51 // | # target element selection in composite
52 // | # target shell selection in element,
53 // | # and sampling of bound electron momentu
54 // |
55 // *******************************************
56 // |
57 // | History:
58 // | --------
59 // |
60 // | Jan. 2015 JMCB - 1st Version based
61 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi
62 // | Nov. 2016 JMCB - Polarisation track
63 // | of Dr. Merlin Reyn
64 // | University of Gene
65 // |
66 // *******************************************
67
68 #include "G4LowEPPolarizedComptonModel.hh"
69 #include "G4AutoLock.hh"
70 #include "G4PhysicalConstants.hh"
71 #include "G4SystemOfUnits.hh"
72
73 //********************************************
74
75 using namespace std;
76 namespace { G4Mutex LowEPPolarizedComptonModel
77
78
79 G4PhysicsFreeVector* G4LowEPPolarizedComptonMo
80 G4ShellData* G4LowEPPolarizedComptonMode
81 G4DopplerProfile* G4LowEPPolarizedComptonMode
82
83 static const G4double ln10 = G4Log(10.);
84
85 G4LowEPPolarizedComptonModel::G4LowEPPolarized
86
87 : G4VEmModel(nam),isInitialised(false)
88 {
89 verboseLevel=1 ;
90 // Verbosity scale:
91 // 0 = nothing
92 // 1 = warning for energy non-conservation
93 // 2 = details of energy budget
94 // 3 = calculation of cross sections, file o
95 // 4 = entering in methods
96
97 if( verboseLevel>1 ) {
98 G4cout << "Low energy photon Compton model
99 }
100
101 //Mark this model as "applicable" for atomic
102 SetDeexcitationFlag(true);
103
104 fParticleChange = nullptr;
105 fAtomDeexcitation = nullptr;
106 }
107
108 //********************************************
109
110 G4LowEPPolarizedComptonModel::~G4LowEPPolarize
111 {
112 if(IsMaster()) {
113 delete shellData;
114 shellData = nullptr;
115 delete profileData;
116 profileData = nullptr;
117 }
118 }
119
120 //********************************************
121
122 void G4LowEPPolarizedComptonModel::Initialise(
123 const
124 {
125 if (verboseLevel > 1) {
126 G4cout << "Calling G4LowEPPolarizedCompton
127 }
128
129 // Initialise element selector
130 if(IsMaster()) {
131 // Access to elements
132 const char* path = G4FindDataDir("G4LEDATA
133
134 G4ProductionCutsTable* theCoupleTable =
135 G4ProductionCutsTable::GetProductionCuts
136 G4int numOfCouples = (G4int)theCoupleTable
137
138 for(G4int i=0; i<numOfCouples; ++i) {
139 const G4Material* material =
140 theCoupleTable->GetMaterialCutsCouple(
141 const G4ElementVector* theElementVector
142 std::size_t nelm = material->GetNumberOf
143
144 for (std::size_t j=0; j<nelm; ++j) {
145 G4int Z = G4lrint((*theElementVector)[
146 if(Z < 1) { Z = 1; }
147 else if(Z > maxZ){ Z = maxZ; }
148
149 if( (!data[Z]) ) { ReadData(Z, path);
150 }
151 }
152
153 // For Doppler broadening
154 if(!shellData) {
155 shellData = new G4ShellData();
156 shellData->SetOccupancyData();
157 G4String file = "/doppler/shell-doppler"
158 shellData->LoadData(file);
159 }
160 if(!profileData) { profileData = new G4Dop
161
162 InitialiseElementSelectors(particle, cuts)
163 }
164
165 if (verboseLevel > 2) {
166 G4cout << "Loaded cross section files" <<
167 }
168
169 if( verboseLevel>1 ) {
170 G4cout << "G4LowEPPolarizedComptonModel is
171 << "Energy range: "
172 << LowEnergyLimit() / eV << " eV -
173 << HighEnergyLimit() / GeV << " GeV
174 << G4endl;
175 }
176
177 if(isInitialised) { return; }
178
179 fParticleChange = GetParticleChangeForGamma(
180 fAtomDeexcitation = G4LossTableManager::Ins
181 isInitialised = true;
182 }
183
184 //********************************************
185
186 void G4LowEPPolarizedComptonModel::InitialiseL
187
188 {
189 SetElementSelectors(masterModel->GetElementS
190 }
191
192 //********************************************
193
194 void G4LowEPPolarizedComptonModel::ReadData(st
195 {
196 if (verboseLevel > 1)
197 {
198 G4cout << "G4LowEPPolarizedComptonModel::R
199 << G4endl;
200 }
201 if(data[Z]) { return; }
202 const char* datadir = path;
203 if(!datadir)
204 {
205 datadir = G4FindDataDir("G4LEDATA");
206 if(!datadir)
207 {
208 G4Exception("G4LowEPPolarizedComptonMode
209 "em0006",FatalException,
210 "Environment variable G4LEDA
211 return;
212 }
213 }
214
215 data[Z] = new G4PhysicsFreeVector();
216
217 std::ostringstream ost;
218 ost << datadir << "/livermore/comp/ce-cs-" <
219 std::ifstream fin(ost.str().c_str());
220
221 if( !fin.is_open())
222 {
223 G4ExceptionDescription ed;
224 ed << "G4LowEPPolarizedComptonModel data
225 << "> is not opened!" << G4endl;
226 G4Exception("G4LowEPPolarizedComptonModel:
227 "em0003",FatalException,
228 ed,"G4LEDATA version should be
229 return;
230 } else {
231 if(verboseLevel > 3) {
232 G4cout << "File " << ost.str()
233 << " is opened by G4LowEPPolari
234 }
235 data[Z]->Retrieve(fin, true);
236 data[Z]->ScaleVector(MeV, MeV*barn);
237 }
238 fin.close();
239 }
240
241 //********************************************
242
243 G4double
244 G4LowEPPolarizedComptonModel::ComputeCrossSect
245
246
247
248 {
249 if (verboseLevel > 3) {
250 G4cout << "G4LowEPPolarizedComptonModel::C
251 << G4endl;
252 }
253 G4double cs = 0.0;
254
255 if (GammaEnergy < LowEnergyLimit()) { return
256
257 G4int intZ = G4lrint(Z);
258 if(intZ < 1 || intZ > maxZ) { return cs; }
259
260 G4PhysicsFreeVector* pv = data[intZ];
261
262 // if element was not initialised
263 // do initialisation safely for MT mode
264 if(!pv)
265 {
266 InitialiseForElement(0, intZ);
267 pv = data[intZ];
268 if(!pv) { return cs; }
269 }
270
271 G4int n = G4int(pv->GetVectorLength() - 1);
272 G4double e1 = pv->Energy(0);
273 G4double e2 = pv->Energy(n);
274
275 if(GammaEnergy <= e1) { cs = GammaEnerg
276 else if(GammaEnergy <= e2) { cs = pv->Value(
277 else if(GammaEnergy > e2) { cs = pv->Value(
278
279 return cs;
280 }
281
282 //********************************************
283
284 void G4LowEPPolarizedComptonModel::SampleSecon
285 const G4MaterialCutsCouple* c
286 const G4DynamicParticle* aDyn
287 G4double, G4double)
288 {
289
290 //Determine number of digits (in decimal bas
291 G4double g4d_order = G4double(numeric_limits
292 G4double g4d_limit = std::pow(10.,-g4d_order
293
294 // The scattered gamma energy is sampled acc
295 // then accepted or rejected depending on th
296 // by factor from Klein - Nishina formula.
297 // Expression of the angular distribution as
298 // angular and energy distribution and Scatt
299 // D. E. Cullen "A simple model of photon tr
300 // Phys. Res. B 101 (1995). Method of sampli
301 // data are interpolated while in the articl
302 // Reference to the article is from J. Stepa
303 // and Electron Interaction Data for GEANT i
304 // TeV (draft).
305 // The random number techniques of Butcher &
306 // (Nucl Phys 20(1960),15).
307
308 G4double photonEnergy0 = aDynamicGamma->GetK
309
310 if (verboseLevel > 3) {
311 G4cout << "G4LowEPPolarizedComptonModel::S
312 << photonEnergy0/MeV << " in " << c
313 << G4endl;
314 }
315 // do nothing below the threshold
316 // should never get here because the XS is z
317 if (photonEnergy0 < LowEnergyLimit())
318 return ;
319
320 G4double e0m = photonEnergy0 / electron_mass
321 G4ParticleMomentum photonDirection0 = aDynam
322
323 // Polarisation: check orientation of photon
324 // Fix if needed
325
326 G4ThreeVector photonPolarization0 = aDynamic
327 // Check if polarisation vector is perpendic
328
329 if (!(photonPolarization0.isOrthogonal(photo
330 {
331 photonPolarization0 = GetRandomPolarizat
332 }
333 else
334 {
335 if ((photonPolarization0.howOrthogonal(p
336 {
337 photonPolarization0 = GetPerpendicul
338 }
339 }
340
341 // Select randomly one element in the curren
342 const G4ParticleDefinition* particle = aDyn
343 const G4Element* elm = SelectRandomAtom(coup
344 G4int Z = (G4int)elm->GetZ();
345
346 G4double LowEPPCepsilon0 = 1. / (1. + 2. * e
347 G4double LowEPPCepsilon0Sq = LowEPPCepsilon0
348 G4double alpha1 = -std::log(LowEPPCepsilon0)
349 G4double alpha2 = 0.5 * (1. - LowEPPCepsilon
350
351 G4double wlPhoton = h_Planck*c_light/photonE
352
353 // Sample the energy of the scattered photon
354 G4double LowEPPCepsilon;
355 G4double LowEPPCepsilonSq;
356 G4double oneCosT;
357 G4double sinT2;
358 G4double gReject;
359
360 if (verboseLevel > 3) {
361 G4cout << "Started loop to sample gamma en
362 }
363
364 do
365 {
366 if ( alpha1/(alpha1+alpha2) > G4UniformR
367 {
368 LowEPPCepsilon = G4Exp(-alpha1 * G4U
369 LowEPPCepsilonSq = LowEPPCepsilon *
370 }
371 else
372 {
373 LowEPPCepsilonSq = LowEPPCepsilon0Sq
374 LowEPPCepsilon = std::sqrt(LowEPPCep
375 }
376
377 oneCosT = (1. - LowEPPCepsilon) / ( LowE
378 sinT2 = oneCosT * (2. - oneCosT);
379 G4double x = std::sqrt(oneCosT/2.) / (wl
380 G4double scatteringFunction = ComputeSca
381 gReject = (1. - LowEPPCepsilon * sinT2 /
382
383 } while(gReject < G4UniformRand()*Z);
384
385 G4double cosTheta = 1. - oneCosT;
386 G4double sinTheta = std::sqrt(sinT2);
387 G4double phi = SetPhi(LowEPPCepsilon,sinT2);
388 G4double dirx = sinTheta * std::cos(phi);
389 G4double diry = sinTheta * std::sin(phi);
390 G4double dirz = cosTheta ;
391
392 // Set outgoing photon polarization
393
394 G4ThreeVector photonPolarization1 = SetNewPo
395 sinT2,
396 phi,
397 cosTheta);
398
399 // Scatter photon energy and Compton electro
400 // J. M. C. Brown, M. R. Dimmock, J. E. Gill
401 // "A low energy bound atomic electron Compt
402 // NIMB, Vol. 338, 77-88, 2014.
403
404 // Set constants and initialize scattering p
405 const G4double vel_c = c_light / (m/s);
406 const G4double momentum_au_to_nat = halfpi*
407 const G4double e_mass_kg = electron_mass_c2
408
409 const G4int maxDopplerIterations = 1000;
410 G4double bindingE = 0.;
411 G4double pEIncident = photonEnergy0 ;
412 G4double pERecoil = -1.;
413 G4double eERecoil = -1.;
414 G4double e_alpha =0.;
415 G4double e_beta = 0.;
416
417 G4double CE_emission_flag = 0.;
418 G4double ePAU = -1;
419 G4int shellIdx = 0;
420 G4double u_temp = 0;
421 G4double cosPhiE =0;
422 G4double sinThetaE =0;
423 G4double cosThetaE =0;
424 G4int iteration = 0;
425
426 if (verboseLevel > 3) {
427 G4cout << "Started loop to sample photon e
428 }
429
430 do{
431 // ***************************************
432 // | Determine scatter photon energy
433 // ***************************************
434 do
435 {
436 iteration++;
437
438 // *****************************************
439 // | Sample bound electron information
440 // *****************************************
441
442 // Select shell based on shell occupancy
443 shellIdx = shellData->SelectRandomShell(Z);
444 bindingE = shellData->BindingEnergy(Z,shellI
445
446 // Randomly sample bound electron momentum (
447 ePAU = profileData->RandomSelectMomentum(Z,s
448
449 // Convert to SI units
450 G4double ePSI = ePAU * momentum_au_to_nat;
451
452 //Calculate bound electron velocity and norm
453 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) /
454
455 // Sample incident electron direction, amorp
456 e_alpha = pi*G4UniformRand();
457 e_beta = twopi*G4UniformRand();
458
459 // Total energy of system
460 G4double eEIncident = electron_mass_c2 / sqr
461 G4double systemE = eEIncident + pEIncident;
462
463 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem
464 G4double numerator = gamma_temp*electron_mas
465 G4double subdenom1 = u_temp*cosTheta*std::c
466 G4double subdenom2 = u_temp*sinTheta*std::si
467 G4double denominator = (1.0 - cosTheta) + (
468 pERecoil = (numerator/denominator);
469 eERecoil = systemE - pERecoil;
470 CE_emission_flag = pEIncident - pERecoil;
471 } while ( (iteration <= maxDopplerIterat
472
473 // End of recalculation of photon energy w
474 // ***************************************
475 // | Determine ejected Compton electro
476 // ***************************************
477
478 // Calculate velocity of ejected Compton e
479
480 G4double a_temp = eERecoil / electron_mass
481 G4double u_p_temp = sqrt(1 - (1 / (a_temp*
482
483 // Coefficients and terms from simulatenou
484 G4double sinAlpha = std::sin(e_alpha);
485 G4double cosAlpha = std::cos(e_alpha);
486 G4double sinBeta = std::sin(e_beta);
487 G4double cosBeta = std::cos(e_beta);
488
489 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_
490 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem
491
492 G4double var_A = pERecoil*u_p_temp*sinThet
493 G4double var_B = u_p_temp* (pERecoil*cosTh
494 G4double var_C = (pERecoil-pEIncident) - (
495
496 G4double var_D1 = gamma*electron_mass_c2*p
497 G4double var_D2 = (1 - (u_temp*cosTheta*co
498 G4double var_D3 = ((electron_mass_c2*elect
499 G4double var_D = var_D1*var_D2 + var_D3;
500
501 G4double var_E1 = ((gamma*gamma_p)*(electr
502 G4double var_E2 = gamma_p*electron_mass_c2
503 G4double var_E = var_E1 - var_E2;
504
505 G4double var_F1 = ((gamma*gamma_p)*(electr
506 G4double var_F2 = (gamma_p*electron_mass_c
507 G4double var_F = var_F1 - var_F2;
508
509 G4double var_G = (gamma*gamma_p)*(electron
510
511 // Two equations form a quadratic form of
512 // Coefficents and solution to quadratic
513 G4double var_W1 = (var_F*var_B - var_E*var
514 G4double var_W2 = (var_G*var_G)*(var_A*var
515 G4double var_W = var_W1 + var_W2;
516
517 G4double var_Y = 2.0*(((var_A*var_D-var_F*
518
519 G4double var_Z1 = (var_A*var_D - var_F*var
520 G4double var_Z2 = (var_G*var_G)*(var_C*var
521 G4double var_Z = var_Z1 + var_Z2;
522 G4double diff1 = var_Y*var_Y;
523 G4double diff2 = 4*var_W*var_Z;
524 G4double diff = diff1 - diff2;
525
526 // Check if diff is less than zero, if so
527 //Confirm that diff less than zero is due
528 //than 10^(-g4d_order), then set diff to z
529
530 if ((diff < 0.0) && (abs(diff / diff1) < g
531 {
532 diff = 0.0;
533 }
534
535 // Plus and minus of quadratic
536 G4double X_p = (-var_Y + sqrt (diff))/(2*v
537 G4double X_m = (-var_Y - sqrt (diff))/(2*v
538
539 // Floating point precision protection
540 // Check if X_p and X_m are greater than o
541 // Issue due to propagation of FPE and onl
542 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;}
543 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;}
544 // End of FP protection
545
546 G4double ThetaE = 0.;
547
548 // Randomly sample one of the two possible
549 G4double sol_select = G4UniformRand();
550
551 if (sol_select < 0.5)
552 {
553 ThetaE = std::acos(X_p);
554 }
555 if (sol_select > 0.5)
556 {
557 ThetaE = std::acos(X_m);
558 }
559
560 cosThetaE = std::cos(ThetaE);
561 sinThetaE = std::sin(ThetaE);
562 G4double Theta = std::acos(cosTheta);
563
564 //Calculate electron Phi
565 G4double iSinThetaE = std::sqrt(1+std::tan
566 G4double iSinTheta = std::sqrt(1+std::tan(
567 G4double ivar_A = iSinTheta/ (pERecoil*u_p
568 // Trigs
569 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_
570 // End of calculation of ejection Compton
571 //Fix for floating point errors
572 } while ( (iteration <= maxDopplerIterations
573
574 // Revert to original if maximum number of i
575 if (iteration >= maxDopplerIterations)
576 {
577 pERecoil = photonEnergy0 ;
578 bindingE = 0.;
579 dirx=0.0;
580 diry=0.0;
581 dirz=1.0;
582 }
583
584 // Set "scattered" photon direction and ener
585 G4ThreeVector photonDirection1(dirx,diry,dir
586 SystemOfRefChange(photonDirection0,photonDir
587 photonPolarization0,photonPolarization
588
589 if (pERecoil > 0.)
590 {
591 fParticleChange->SetProposedKineticEnerg
592 fParticleChange->ProposeMomentumDirectio
593 fParticleChange->ProposePolarization(pho
594
595 // Set ejected Compton electron directio
596 G4double PhiE = std::acos(cosPhiE);
597 G4double eDirX = sinThetaE * std::cos(ph
598 G4double eDirY = sinThetaE * std::sin(ph
599 G4double eDirZ = cosThetaE;
600
601 G4double eKineticEnergy = pEIncident - p
602
603 G4ThreeVector eDirection(eDirX,eDirY,eDi
604 SystemOfRefChangeElect(photonDirection0,
605 photonPolarization0);
606
607 G4DynamicParticle* dp = new G4DynamicPar
608 eDirection,eKineticEnergy) ;
609 fvect->push_back(dp);
610 }
611 else
612 {
613 fParticleChange->SetProposedKineticEnerg
614 fParticleChange->ProposeTrackStatus(fSto
615 }
616
617 // sample deexcitation
618 //
619 if (verboseLevel > 3) {
620 G4cout << "Started atomic de-excitation "
621 }
622
623 if(fAtomDeexcitation && iteration < maxDoppl
624 G4int index = couple->GetIndex();
625 if(fAtomDeexcitation->CheckDeexcitationAct
626 std::size_t nbefore = fvect->size();
627 G4AtomicShellEnumerator as = G4AtomicShe
628 const G4AtomicShell* shell = fAtomDeexci
629 fAtomDeexcitation->GenerateParticles(fve
630 std::size_t nafter = fvect->size();
631 if(nafter > nbefore) {
632 for (std::size_t i=nbefore; i<nafter;
633 //Check if there is enough residual
634 if (bindingE >= ((*fvect)[i])->GetKi
635 {
636 //Ok, this is a valid secondary: keep
637 bindingE -= ((*fvect)[i])->GetKineticE
638 }
639 else
640 {
641 //Invalid secondary: not enough energy
642 //Keep its energy in the local deposit
643 delete (*fvect)[i];
644 (*fvect)[i]=nullptr;
645 }
646 }
647 }
648 }
649 }
650
651 //This should never happen
652 if(bindingE < 0.0)
653 G4Exception("G4LowEPPolarizedComptonModel:
654 "em2051",FatalException,"Negative local en
655
656 fParticleChange->ProposeLocalEnergyDeposit(b
657 }
658
659 //********************************************
660
661 G4double
662 G4LowEPPolarizedComptonModel::ComputeScatterin
663 {
664 G4double value = Z;
665 if (x <= ScatFuncFitParam[Z][2]) {
666
667 G4double lgq = G4Log(x)/ln10;
668
669 if (lgq < ScatFuncFitParam[Z][1]) {
670 value = ScatFuncFitParam[Z][3] + lgq*Sca
671 } else {
672 value = ScatFuncFitParam[Z][5] + lgq*Sca
673 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l
674 }
675 value = G4Exp(value*ln10);
676 }
677 return value;
678 }
679
680
681 //********************************************
682
683 void
684 G4LowEPPolarizedComptonModel::InitialiseForEle
685 G4int Z)
686 {
687 G4AutoLock l(&LowEPPolarizedComptonModelMute
688 if(!data[Z]) { ReadData(Z); }
689 l.unlock();
690 }
691
692 //********************************************
693
694 //Fitting data to compute scattering function
695 const G4double G4LowEPPolarizedComptonModel::S
696 { 0, 0., 0., 0., 0.,
697 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -1
698 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -5
699 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -6
700 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -1
701 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -1
702 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -1
703 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -1
704 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -1
705 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -1
706 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -1
707 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -4
708 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -5
709 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -6
710 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -7
711 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -8
712 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -9
713 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -9
714 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -1
715 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -5
716 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -5
717 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -5
718 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -5
719 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -5
720 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -6
721 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -5
722 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -6
723 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -6
724 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -6
725 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -6
726 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -6
727 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -6
728 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -6
729 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -6
730 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -6
731 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -7
732 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -7
733 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -4
734 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -4
735 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -4
736 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -4
737 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -5
738 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -5
739 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -5
740 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -5
741 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -5
742 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -7
743 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -6
744 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -5
745 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -5
746 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -5
747 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -6
748 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -6
749 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -6
750 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -6
751 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -4
752 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -3
753 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -3
754 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -4
755 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -4
756 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -3
757 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -4
758 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -4
759 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -4
760 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -4
761 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -4
762 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -4
763 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -4
764 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -4
765 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -4
766 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -4
767 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -4
768 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -4
769 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -4
770 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -4
771 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -4
772 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -4
773 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -4
774 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -5
775 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -5
776 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -4
777 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -4
778 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -5
779 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -5
780 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -5
781 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -5
782 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -5
783 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -3
784 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -3
785 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -3
786 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -3
787 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -3
788 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -3
789 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -3
790 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -3
791 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -3
792 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -3
793 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -4
794 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -3
795 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -3
796 {100, 6.412, 1.00000E+10, -12.900, 1.993, -3
797 };
798
799 //********************************************
800
801 //Supporting functions for photon polarisation
802 G4double G4LowEPPolarizedComptonModel::SetPhi(
803 G4double sinT2)
804 {
805 G4double rand1;
806 G4double rand2;
807 G4double phiProbability;
808 G4double phi;
809 G4double a, b;
810
811 do
812 {
813 rand1 = G4UniformRand();
814 rand2 = G4UniformRand();
815 phiProbability=0.;
816 phi = twopi*rand1;
817
818 a = 2*sinT2;
819 b = energyRate + 1/energyRate;
820
821 phiProbability = 1 - (a/b)*(std::cos(phi
822 }
823 while ( rand2 > phiProbability );
824 return phi;
825 }
826
827 //********************************************
828
829 G4ThreeVector G4LowEPPolarizedComptonModel::Se
830 {
831 G4double dx = a.x();
832 G4double dy = a.y();
833 G4double dz = a.z();
834 G4double x = dx < 0.0 ? -dx : dx;
835 G4double y = dy < 0.0 ? -dy : dy;
836 G4double z = dz < 0.0 ? -dz : dz;
837 if (x < y) {
838 return x < z ? G4ThreeVector(-dy,dx,0) : G
839 }else{
840 return y < z ? G4ThreeVector(dz,0,-dx) : G
841 }
842 }
843
844 //********************************************
845
846 G4ThreeVector G4LowEPPolarizedComptonModel::Ge
847 {
848 G4ThreeVector d0 = direction0.unit();
849 G4ThreeVector a1 = SetPerpendicularVector(d0
850 G4ThreeVector a0 = a1.unit(); // unit vector
851
852 G4double rand1 = G4UniformRand();
853
854 G4double angle = twopi*rand1; // random pola
855 G4ThreeVector b0 = d0.cross(a0); // cross pr
856
857 G4ThreeVector c;
858
859 c.setX(std::cos(angle)*(a0.x())+std::sin(ang
860 c.setY(std::cos(angle)*(a0.y())+std::sin(ang
861 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang
862
863 G4ThreeVector c0 = c.unit();
864
865 return c0;
866 }
867
868 //********************************************
869
870 G4ThreeVector G4LowEPPolarizedComptonModel::Ge
871 (const G4ThreeVector& photonDirection, const G
872 {
873 //
874 // The polarization of a photon is always pe
875 // Therefore this function removes those vec
876 // points in direction of photonDirection
877 //
878 // Mathematically we search the projection o
879 // plains normal vector.
880 // The basic equation can be found in each g
881 // p = a - (a o n)/(n o n)*n
882
883 return photonPolarization - photonPolarizati
884 }
885
886 //********************************************
887
888 G4ThreeVector G4LowEPPolarizedComptonModel::Se
889 G4double sinT2,
890 G4double phi,
891 G4double costheta)
892 {
893 G4double rand1;
894 G4double rand2;
895 G4double cosPhi = std::cos(phi);
896 G4double sinPhi = std::sin(phi);
897 G4double sinTheta = std::sqrt(sinT2);
898 G4double cosP2 = cosPhi*cosPhi;
899 G4double normalisation = std::sqrt(1. - cosP
900
901 // Method based on:
902 // D. Xu, Z. He and F. Zhang
903 // "Detection of Gamma Ray Polarization Usin
904 // IEEE TNS, Vol. 52(4), 1160-1164, 2005.
905
906 // Determination of Theta
907
908 G4double theta;
909 rand1 = G4UniformRand();
910 rand2 = G4UniformRand();
911
912 if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon
913 {
914 if (rand2<0.5)
915 theta = pi/2.0;
916 else
917 theta = 3.0*pi/2.0;
918 }
919 else
920 {
921 if (rand2<0.5)
922 theta = 0;
923 else
924 theta = pi;
925 }
926 G4double cosBeta = std::cos(theta);
927 G4double sinBeta = std::sqrt(1-cosBeta*cosBe
928
929 G4ThreeVector photonPolarization1;
930
931 G4double xParallel = normalisation*cosBeta;
932 G4double yParallel = -(sinT2*cosPhi*sinPhi)*
933 G4double zParallel = -(costheta*sinTheta*cos
934 G4double xPerpendicular = 0.;
935 G4double yPerpendicular = (costheta)*sinBeta
936 G4double zPerpendicular = -(sinTheta*sinPhi)
937
938 G4double xTotal = (xParallel + xPerpendicula
939 G4double yTotal = (yParallel + yPerpendicula
940 G4double zTotal = (zParallel + zPerpendicula
941
942 photonPolarization1.setX(xTotal);
943 photonPolarization1.setY(yTotal);
944 photonPolarization1.setZ(zTotal);
945
946 return photonPolarization1;
947
948 }
949
950 //********************************************
951 void G4LowEPPolarizedComptonModel::SystemOfRef
952 G4ThreeVector& direction1,
953 G4ThreeVector& polarization0,
954 G4ThreeVector& polarization1)
955 {
956 // direction0 is the original photon directi
957 // polarization0 is the original photon pola
958 // need to specify y axis in the real refere
959 G4ThreeVector Axis_Z0 = direction0.unit();
960 G4ThreeVector Axis_X0 = polarization0.unit()
961 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_
962
963 G4double direction_x = direction1.getX();
964 G4double direction_y = direction1.getY();
965 G4double direction_z = direction1.getZ();
966
967 direction1 = (direction_x*Axis_X0 + directio
968 G4double polarization_x = polarization1.getX
969 G4double polarization_y = polarization1.getY
970 G4double polarization_z = polarization1.getZ
971
972 polarization1 = (polarization_x*Axis_X0 + po
973
974 }
975
976 //********************************************
977 void G4LowEPPolarizedComptonModel::SystemOfRef
978 G4ThreeVector& edirection,
979 G4ThreeVector& ppolarization)
980 {
981 // direction0 is the original photon directi
982 // polarization0 is the original photon pola
983 // need to specify y axis in the real refere
984 G4ThreeVector Axis_Z0 = pdirection.unit();
985 G4ThreeVector Axis_X0 = ppolarization.unit()
986 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_
987
988 G4double direction_x = edirection.getX();
989 G4double direction_y = edirection.getY();
990 G4double direction_z = edirection.getZ();
991
992 edirection = (direction_x*Axis_X0 + directio
993 }
994
995
996