| 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 // History:
30 // -------
31 //
32 // 05 Apr 2021 J Allison added quantum entan
33 // If the photons have been "tagged" as "quant
34 // G4eplusAnnihilation for annihilation into 2
35 // here if - and only if - both photons suffer
36 // predictions from Pryce and Ward, Nature No
37 // Physical Review 73 (1948) p.440. Experiment
38 // entanglement in the MeV regime and its appl
39 // D. Watts, J. Allison et al., Nature Communi
40 // https://doi.org/10.1038/s41467-021-22907-5.
41 //
42 // 02 May 2009 S Incerti as V. Ivanchenko pr
43 //
44 // Cleanup initialisation and generation of se
45 // - apply internal high-ener
46 // - do not apply low-energy
47 // - remove GetMeanFreePath m
48 // - added protection against
49 // - use G4ElementSelector
50
51 #include "G4LivermorePolarizedComptonModel.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4AutoLock.hh"
55 #include "G4Electron.hh"
56 #include "G4ParticleChangeForGamma.hh"
57 #include "G4LossTableManager.hh"
58 #include "G4VAtomDeexcitation.hh"
59 #include "G4AtomicShell.hh"
60 #include "G4Gamma.hh"
61 #include "G4ShellData.hh"
62 #include "G4DopplerProfile.hh"
63 #include "G4Log.hh"
64 #include "G4Exp.hh"
65 #include "G4Pow.hh"
66 #include "G4LogLogInterpolation.hh"
67 #include "G4PhysicsModelCatalog.hh"
68 #include "G4EntanglementAuxInfo.hh"
69 #include "G4eplusAnnihilationEntanglementClipB
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 using namespace std;
74 namespace { G4Mutex LivermorePolarizedComptonM
75
76
77 G4PhysicsFreeVector* G4LivermorePolarizedCompt
78 G4ShellData* G4LivermorePolarizedCompton
79 G4DopplerProfile* G4LivermorePolarizedCompton
80 G4CompositeEMDataSet* G4LivermorePolarizedComp
81
82 //....oooOO0OOooo........oooOO0OOooo........oo
83
84 G4LivermorePolarizedComptonModel::G4LivermoreP
85 :G4VEmModel(nam),isInitialised(false)
86 {
87 verboseLevel= 1;
88 // Verbosity scale:
89 // 0 = nothing
90 // 1 = warning for energy non-conservation
91 // 2 = details of energy budget
92 // 3 = calculation of cross sections, file o
93 // 4 = entering in methods
94
95 if( verboseLevel>1 )
96 G4cout << "Livermore Polarized Compton is
97
98 //Mark this model as "applicable" for atomic
99 SetDeexcitationFlag(true);
100
101 fParticleChange = nullptr;
102 fAtomDeexcitation = nullptr;
103 fEntanglementModelID = G4PhysicsModelCatalog
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 G4LivermorePolarizedComptonModel::~G4Livermore
109 {
110 if(IsMaster()) {
111 delete shellData;
112 shellData = nullptr;
113 delete profileData;
114 profileData = nullptr;
115 delete scatterFunctionData;
116 scatterFunctionData = nullptr;
117 for(G4int i=0; i<maxZ; ++i) {
118 if(data[i]) {
119 delete data[i];
120 data[i] = nullptr;
121 }
122 }
123 }
124 }
125
126 //....oooOO0OOooo........oooOO0OOooo........oo
127
128 void G4LivermorePolarizedComptonModel::Initial
129 const G
130 {
131 if (verboseLevel > 1)
132 G4cout << "Calling G4LivermorePolarizedCom
133
134 // Initialise element selector
135 if(IsMaster()) {
136 // Access to elements
137 const char* path = G4FindDataDir("G4LEDATA
138
139 G4ProductionCutsTable* theCoupleTable =
140 G4ProductionCutsTable::GetProductionCuts
141
142 G4int numOfCouples = (G4int)theCoupleTable
143
144 for(G4int i=0; i<numOfCouples; ++i) {
145 const G4Material* material =
146 theCoupleTable->GetMaterialCutsCouple(
147 const G4ElementVector* theElementVector
148 std::size_t nelm = material->GetNumberOf
149
150 for (std::size_t j=0; j<nelm; ++j) {
151 G4int Z = G4lrint((*theElementVector)[
152 if(Z < 1) { Z = 1; }
153 else if(Z > maxZ){ Z = maxZ; }
154
155 if( (!data[Z]) ) { ReadData(Z, path);
156 }
157 }
158
159 // For Doppler broadening
160 if(!shellData) {
161 shellData = new G4ShellData();
162 shellData->SetOccupancyData();
163 G4String file = "/doppler/shell-doppler"
164 shellData->LoadData(file);
165 }
166 if(!profileData) { profileData = new G4Dop
167
168 // Scattering Function
169 if(!scatterFunctionData)
170 {
171
172 G4VDataSetAlgorithm* scatterInterpolation =
173 G4String scatterFile = "comp/ce-sf-";
174 scatterFunctionData = new G4CompositeEMDataS
175 scatterFunctionData->LoadData(scatterFile);
176 }
177
178 InitialiseElementSelectors(particle, cuts)
179 }
180
181 if (verboseLevel > 2) {
182 G4cout << "Loaded cross section files" <<
183 }
184
185 if( verboseLevel>1 ) {
186 G4cout << "G4LivermoreComptonModel is init
187 << "Energy range: "
188 << LowEnergyLimit() / eV << " eV - "
189 << HighEnergyLimit() / GeV << " GeV"
190 << G4endl;
191 }
192 //
193 if(isInitialised) { return; }
194
195 fParticleChange = GetParticleChangeForGamma(
196 fAtomDeexcitation = G4LossTableManager::Ins
197 isInitialised = true;
198 }
199
200
201 void G4LivermorePolarizedComptonModel::Initial
202 G4VEmModel* masterModel)
203 {
204 SetElementSelectors(masterModel->GetElementS
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oo
208
209 void G4LivermorePolarizedComptonModel::ReadDat
210 {
211 if (verboseLevel > 1)
212 {
213 G4cout << "G4LivermorePolarizedComptonMo
214 << G4endl;
215 }
216 if(data[Z]) { return; }
217 const char* datadir = path;
218 if(!datadir)
219 {
220 datadir = G4FindDataDir("G4LEDATA");
221 if(!datadir)
222 {
223 G4Exception("G4LivermorePolarizedComptonMo
224 "em0006",FatalException,
225 "Environment variable G4LEDATA not d
226 return;
227 }
228 }
229
230 data[Z] = new G4PhysicsFreeVector();
231
232 std::ostringstream ost;
233 ost << datadir << "/livermore/comp/ce-cs-" <
234 std::ifstream fin(ost.str().c_str());
235
236 if( !fin.is_open())
237 {
238 G4ExceptionDescription ed;
239 ed << "G4LivermorePolarizedComptonModel
240 << "> is not opened!" << G4endl;
241 G4Exception("G4LivermoreComptonModel::Re
242 "em0003",FatalException,
243 ed,"G4LEDATA version should be G4EMLOW8.
244 return;
245 } else {
246 if(verboseLevel > 3) {
247 G4cout << "File " << ost.str()
248 << " is opened by G4LivermorePolarizedC
249 }
250 data[Z]->Retrieve(fin, true);
251 data[Z]->ScaleVector(MeV, MeV*barn);
252 }
253 fin.close();
254 }
255
256 //....oooOO0OOooo........oooOO0OOooo........oo
257
258 G4double G4LivermorePolarizedComptonModel::Com
259 const G
260 G
261 G
262 G
263 {
264 if (verboseLevel > 3)
265 G4cout << "Calling ComputeCrossSectionPerA
266
267 G4double cs = 0.0;
268
269 if (GammaEnergy < LowEnergyLimit())
270 return 0.0;
271
272 G4int intZ = G4lrint(Z);
273 if(intZ < 1 || intZ > maxZ) { return cs; }
274
275 G4PhysicsFreeVector* pv = data[intZ];
276
277 // if element was not initialised
278 // do initialisation safely for MT mode
279 if(!pv)
280 {
281 InitialiseForElement(0, intZ);
282 pv = data[intZ];
283 if(!pv) { return cs; }
284 }
285
286 G4int n = G4int(pv->GetVectorLength() - 1);
287 G4double e1 = pv->Energy(0);
288 G4double e2 = pv->Energy(n);
289
290 if(GammaEnergy <= e1) { cs = GammaEnerg
291 else if(GammaEnergy <= e2) { cs = pv->Value(
292 else if(GammaEnergy > e2) { cs = pv->Value(
293
294 return cs;
295
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oo
299
300 void G4LivermorePolarizedComptonModel::SampleS
301 const G4MaterialCutsCouple* co
302 const G4DynamicParticle* aDyna
303 G4double,
304 G4double)
305 {
306 // The scattered gamma energy is sampled acc
307 // The random number techniques of Butcher &
308 // GEANT4 internal units
309 //
310 // Note : Effects due to binding of atomic e
311
312 if (verboseLevel > 3)
313 G4cout << "Calling SampleSecondaries() of
314
315 G4double gammaEnergy0 = aDynamicGamma->GetKi
316
317 // do nothing below the threshold
318 // should never get here because the XS is z
319 if (gammaEnergy0 < LowEnergyLimit())
320 return ;
321
322 G4ThreeVector gammaPolarization0 = aDynamicG
323
324 // Protection: a polarisation parallel to th
325 // direction causes problems;
326 // in that case find a random polarization
327 G4ThreeVector gammaDirection0 = aDynamicGamm
328
329 // Make sure that the polarization vector is
330 // gamma direction. If not
331 if(!(gammaPolarization0.isOrthogonal(gammaDi
332 { // only for testing now
333 gammaPolarization0 = GetRandomPolarizati
334 }
335 else
336 {
337 if ( gammaPolarization0.howOrthogonal(ga
338 {
339 gammaPolarization0 = GetPerpendicularPolar
340 }
341 }
342 // End of Protection
343
344 G4double E0_m = gammaEnergy0 / electron_mass
345
346 // Select randomly one element in the curren
347 //G4int Z = crossSectionHandler->SelectRando
348 const G4ParticleDefinition* particle = aDyn
349 const G4Element* elm = SelectRandomAtom(coup
350 G4int Z = (G4int)elm->GetZ();
351
352 // Sample the energy and the polarization of
353 G4double epsilon, epsilonSq, onecost, sinThe
354
355 G4double epsilon0Local = 1./(1. + 2*E0_m);
356 G4double epsilon0Sq = epsilon0Local*epsilon0
357 G4double alpha1 = - G4Log(epsilon0Local);
358 G4double alpha2 = 0.5*(1.- epsilon0Sq);
359
360 G4double wlGamma = h_Planck*c_light/gammaEne
361 G4double gammaEnergy1;
362 G4ThreeVector gammaDirection1;
363
364 do {
365 if ( alpha1/(alpha1+alpha2) > G4UniformRan
366 {
367 epsilon = G4Exp(-alpha1*G4UniformRand());
368 epsilonSq = epsilon*epsilon;
369 }
370 else
371 {
372 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4
373 epsilon = std::sqrt(epsilonSq);
374 }
375
376 onecost = (1.- epsilon)/(epsilon*E0_m);
377 sinThetaSqr = onecost*(2.-onecost);
378
379 // Protection
380 if (sinThetaSqr > 1.)
381 {
382 G4cout
383 << " -- Warning -- G4LivermorePolarizedCom
384 << "sin(theta)**2 = "
385 << sinThetaSqr
386 << "; set to 1"
387 << G4endl;
388 sinThetaSqr = 1.;
389 }
390 if (sinThetaSqr < 0.)
391 {
392 G4cout
393 << " -- Warning -- G4LivermorePolarizedCom
394 << "sin(theta)**2 = "
395 << sinThetaSqr
396 << "; set to 0"
397 << G4endl;
398 sinThetaSqr = 0.;
399 }
400 // End protection
401
402 G4double x = std::sqrt(onecost/2.) / (wlG
403 G4double scatteringFunction = scatterFunct
404 greject = (1. - epsilon*sinThetaSqr/(1.+ e
405
406 } while(greject < G4UniformRand()*Z);
407
408
409 // *****************************************
410 // Phi determination
411 // *****************************************
412 G4double phi = SetPhi(epsilon,sinThetaSqr);
413
414 //
415 // scattered gamma angles. ( Z - axis along
416 //
417 G4double cosTheta = 1. - onecost;
418
419 // Protection
420 if (cosTheta > 1.)
421 {
422 G4cout
423 << " -- Warning -- G4LivermorePolarizedCompt
424 << "cosTheta = "
425 << cosTheta
426 << "; set to 1"
427 << G4endl;
428 cosTheta = 1.;
429 }
430 if (cosTheta < -1.)
431 {
432 G4cout
433 << " -- Warning -- G4LivermorePolarizedCompt
434 << "cosTheta = "
435 << cosTheta
436 << "; set to -1"
437 << G4endl;
438 cosTheta = -1.;
439 }
440 // End protection
441
442 G4double sinTheta = std::sqrt (sinThetaSqr);
443
444 // Protection
445 if (sinTheta > 1.)
446 {
447 G4cout
448 << " -- Warning -- G4LivermorePolarizedCompt
449 << "sinTheta = "
450 << sinTheta
451 << "; set to 1"
452 << G4endl;
453 sinTheta = 1.;
454 }
455 if (sinTheta < -1.)
456 {
457 G4cout
458 << " -- Warning -- G4LivermorePolarizedCompt
459 << "sinTheta = "
460 << sinTheta
461 << "; set to -1"
462 << G4endl;
463 sinTheta = -1.;
464 }
465 // End protection
466
467 // Check for entanglement and re-sample phi
468
469 const auto* auxInfo
470 = fParticleChange->GetCurrentTrack()->GetAux
471 if (auxInfo) {
472 const auto* entanglementAuxInfo = dynamic_
473 if (entanglementAuxInfo) {
474 auto* clipBoard = dynamic_cast<G4eplusAn
475 (entanglementAuxInfo->GetEntanglementCli
476 if (clipBoard) {
477 // This is an entangled photon from ep
478 // If this is the first scatter of the
479 // phi on the clipboard.
480 // If this is the first scatter of the
481 // phi of the first scatter of the fir
482 // theta of the second photon, to samp
483 if (clipBoard->IsTrack1Measurement())
484 // Check we have the relevant track.
485 // necessary but I want to be sure t
486 // entangled system are properly pai
487 // Note: the tracking manager pops t
488 // will rely on that. (If not, the
489 // more complicated to ensure we mat
490 // So our track 1 is clipboard track
491 if (clipBoard->GetTrackB() == fParti
492 // This is the first scatter of th
493 // // Debug
494 // auto* track1 = fPart
495 // G4cout
496 // << "This is the firs
497 // << "\nTrack: " << tr
498 // << ", Parent: " << t
499 // << ", Name: " << cli
500 // << G4endl;
501 // // End debug
502 clipBoard->ResetTrack1Measurement(
503 // Store cos(theta),phi of first p
504 clipBoard->SetComptonCosTheta1(cos
505 clipBoard->SetComptonPhi1(phi);
506 }
507 } else if (clipBoard->IsTrack2Measurem
508 // Check we have the relevant track.
509 // Remember our track 2 is clipboard
510 if (clipBoard->GetTrackA() == fParti
511 // This is the first scatter of th
512 // // Debug
513 // auto* track2 = fPart
514 // G4cout
515 // << "This is the firs
516 // << "\nTrack: " << tr
517 // << ", Parent: " << t
518 // << ", Name: " << cli
519 // << G4endl;
520 // // End debug
521 clipBoard->ResetTrack2Measurement(
522
523 // Get cos(theta),phi of first pho
524 const G4double& cosTheta1 = clipBo
525 const G4double& phi1 = clipBoard->
526 // For clarity make aliases for th
527 const G4double& cosTheta2 = cosThe
528 G4double& phi2 = phi;
529 // G4cout << "cosTheta1
530 // G4cout << "cosTheta2
531
532 // Re-sample phi
533 // Draw the difference of azimutha
534 // A + B * cos(2*deltaPhi), or rat
535 // C = A / (A + |B|) and D = B / (
536 const G4double sin2Theta1 = 1.-cos
537 const G4double sin2Theta2 = 1.-cos
538
539 // Pryce and Ward, Nature No 4065
540 auto* g4Pow = G4Pow::GetInstance()
541 const G4double A =
542 ((g4Pow->powN(1.-cosTheta1,3))+2.)
543 ((g4Pow->powN(2.-cosTheta1,3)*g4Po
544 const G4double B = -(sin2Theta1*si
545 ((g4Pow->powN(2.-cosTheta1,2)*g4Po
546
547 // // Snyder et al, Physical Re
548 // // (This is an alternative f
549 // const G4double& k0 = gammaEn
550 // const G4double k1 = k0/(2.-c
551 // const G4double k2 = k0/(2.-c
552 // const G4double gamma1 = k1/k
553 // const G4double gamma2 = k2/k
554 // const G4double A1 = gamma1*g
555 // const G4double B1 = 2.*sin2T
556 // // That's A1 + B1*sin2(delta
557 // const G4double A = A1 + 0.5*
558 // const G4double B = -0.5*B1;
559
560 const G4double maxValue = A + std:
561 const G4double C = A / maxValue;
562 const G4double D = B / maxValue;
563 // G4cout << "A,B,C,D: " << A <
564
565 // Sample delta phi
566 G4double deltaPhi;
567 const G4int maxCount = 999999;
568 G4int iCount = 0;
569 for (; iCount < maxCount; ++iCount
570 deltaPhi = twopi * G4UniformRand
571 if (G4UniformRand() < C + D * co
572 }
573 if (iCount >= maxCount ) {
574 G4cout << "G4LivermorePolarizedC
575 << "Re-sampled delta phi not fou
576 << " tries - carrying on anyway.
577 }
578
579 // Thus, the desired second photon
580 phi2 = deltaPhi - phi1 + halfpi;
581 // The minus sign is in above stat
582 // annihilation photons are in opp
583 // are measured in the opposite di
584 // halfpi is added for the followi
585 // In this function phi is relativ
586 // SystemOfRefChange below. We kno
587 // the polarisations of the two an
588 // to each other, i.e., halfpi dif
589 // Furthermore, only sin(phi) and
590 // need to place any range constra
591 // if (phi2 > pi) {
592 // phi2 -= twopi;
593 // }
594 // if (phi2 < -pi) {
595 // phi2 += twopi;
596 // }
597 }
598 }
599 }
600 }
601 }
602
603 // End of entanglement
604
605 G4double dirx = sinTheta*std::cos(phi);
606 G4double diry = sinTheta*std::sin(phi);
607 G4double dirz = cosTheta ;
608
609 // oneCosT , eom
610
611 // Doppler broadening - Method based on:
612 // Y. Namito, S. Ban and H. Hirayama,
613 // "Implementation of the Doppler Broadening
614 // NIM A 349, pp. 489-494, 1994
615
616 // Maximum number of sampling iterations
617 static G4int maxDopplerIterations = 1000;
618 G4double bindingE = 0.;
619 G4double photonEoriginal = epsilon * gammaEn
620 G4double photonE = -1.;
621 G4int iteration = 0;
622 G4double eMax = gammaEnergy0;
623
624 G4int shellIdx = 0;
625
626 if (verboseLevel > 3) {
627 G4cout << "Started loop to sample broading
628 }
629
630 do
631 {
632 iteration++;
633 // Select shell based on shell occupancy
634 shellIdx = shellData->SelectRandomShell(
635 bindingE = shellData->BindingEnergy(Z,sh
636
637 if (verboseLevel > 3) {
638 G4cout << "Shell ID= " << shellIdx
639 << " Ebind(keV)= " << bindingE/keV <<
640 }
641 eMax = gammaEnergy0 - bindingE;
642
643 // Randomly sample bound electron moment
644 G4double pSample = profileData->RandomSe
645
646 if (verboseLevel > 3) {
647 G4cout << "pSample= " << pSample << G4e
648 }
649 // Rescale from atomic units
650 G4double pDoppler = pSample * fine_struc
651 G4double pDoppler2 = pDoppler * pDoppler
652 G4double var2 = 1. + onecost * E0_m;
653 G4double var3 = var2*var2 - pDoppler2;
654 G4double var4 = var2 - pDoppler2 * cosTh
655 G4double var = var4*var4 - var3 + pDoppl
656 if (var > 0.)
657 {
658 G4double varSqrt = std::sqrt(var);
659 G4double scale = gammaEnergy0 / var3;
660 // Random select either root
661 if (G4UniformRand() < 0.5) photonE = (var4
662 else photonE = (var4 + varSqrt) * scale;
663 }
664 else
665 {
666 photonE = -1.;
667 }
668 } while ( iteration <= maxDopplerIterations
669 (photonE < 0. || photonE > eMax || phot
670
671 // End of recalculation of photon energy wit
672 // Revert to original if maximum number of i
673 if (iteration >= maxDopplerIterations)
674 {
675 photonE = photonEoriginal;
676 bindingE = 0.;
677 }
678
679 gammaEnergy1 = photonE;
680
681 //
682 // update G4VParticleChange for the scattere
683 //
684 // New polarization
685 G4ThreeVector gammaPolarization1 = SetNewPol
686 sinThetaSqr,
687 phi,
688 cosTheta);
689
690 // Set new direction
691 G4ThreeVector tmpDirection1( dirx,diry,dirz
692 gammaDirection1 = tmpDirection1;
693
694 // Change reference frame.
695 SystemOfRefChange(gammaDirection0,gammaDirec
696 gammaPolarization0,gammaPolarization1)
697
698 if (gammaEnergy1 > 0.)
699 {
700 fParticleChange->SetProposedKineticEnerg
701 fParticleChange->ProposeMomentumDirectio
702 fParticleChange->ProposePolarization( ga
703 }
704 else
705 {
706 gammaEnergy1 = 0.;
707 fParticleChange->SetProposedKineticEnerg
708 fParticleChange->ProposeTrackStatus(fSto
709 }
710
711 //
712 // kinematic of the scattered electron
713 //
714 G4double ElecKineEnergy = gammaEnergy0 - gam
715
716 // SI -protection against negative final ene
717 // like in G4LivermoreComptonModel.cc
718 if(ElecKineEnergy < 0.0) {
719 fParticleChange->ProposeLocalEnergyDeposit
720 return;
721 }
722
723 G4double ElecMomentum = std::sqrt(ElecKineEn
724
725 G4ThreeVector ElecDirection((gammaEnergy0 *
726 gammaEnergy1 * gammaDirection1) * (
727
728 G4DynamicParticle* dp =
729 new G4DynamicParticle (G4Electron::Electro
730 fvect->push_back(dp);
731
732 // sample deexcitation
733 //
734 if (verboseLevel > 3) {
735 G4cout << "Started atomic de-excitation "
736 }
737
738 if(fAtomDeexcitation && iteration < maxDoppl
739 G4int index = couple->GetIndex();
740 if(fAtomDeexcitation->CheckDeexcitationAct
741 std::size_t nbefore = fvect->size();
742 G4AtomicShellEnumerator as = G4AtomicShe
743 const G4AtomicShell* shell = fAtomDeexci
744 fAtomDeexcitation->GenerateParticles(fve
745 std::size_t nafter = fvect->size();
746 if(nafter > nbefore) {
747 for (std::size_t i=nbefore; i<nafter; ++i) {
748 //Check if there is enough residual energy
749 if (bindingE >= ((*fvect)[i])->GetKineticE
750 {
751 //Ok, this is a valid secondary:
752 bindingE -= ((*fvect)[i])->GetKineticE
753 }
754 else
755 {
756 //Invalid secondary: not enough energy
757 //Keep its energy in the local deposit
758 delete (*fvect)[i];
759 (*fvect)[i]=0;
760 }
761 }
762 }
763 }
764 }
765 //This should never happen
766 if(bindingE < 0.0)
767 G4Exception("G4LivermoreComptonModel::Samp
768 "em2050",FatalException,"Negative local en
769
770 fParticleChange->ProposeLocalEnergyDeposit(b
771
772 }
773
774 //....oooOO0OOooo........oooOO0OOooo........oo
775
776 G4double G4LivermorePolarizedComptonModel::Set
777 G4double sinSqrTh)
778 {
779 G4double rand1;
780 G4double rand2;
781 G4double phiProbability;
782 G4double phi;
783 G4double a, b;
784
785 do
786 {
787 rand1 = G4UniformRand();
788 rand2 = G4UniformRand();
789 phiProbability=0.;
790 phi = twopi*rand1;
791
792 a = 2*sinSqrTh;
793 b = energyRate + 1/energyRate;
794
795 phiProbability = 1 - (a/b)*(std::cos(phi
796 }
797 while ( rand2 > phiProbability );
798 return phi;
799 }
800
801
802 //....oooOO0OOooo........oooOO0OOooo........oo
803
804 G4ThreeVector G4LivermorePolarizedComptonModel
805 {
806 G4double dx = a.x();
807 G4double dy = a.y();
808 G4double dz = a.z();
809 G4double x = dx < 0.0 ? -dx : dx;
810 G4double y = dy < 0.0 ? -dy : dy;
811 G4double z = dz < 0.0 ? -dz : dz;
812 if (x < y) {
813 return x < z ? G4ThreeVector(-dy,dx,0) : G
814 }else{
815 return y < z ? G4ThreeVector(dz,0,-dx) : G
816 }
817 }
818
819 //....oooOO0OOooo........oooOO0OOooo........oo
820
821 G4ThreeVector G4LivermorePolarizedComptonModel
822 {
823 G4ThreeVector d0 = direction0.unit();
824 G4ThreeVector a1 = SetPerpendicularVector(d0
825 G4ThreeVector a0 = a1.unit(); // unit vector
826
827 G4double rand1 = G4UniformRand();
828
829 G4double angle = twopi*rand1; // random pola
830 G4ThreeVector b0 = d0.cross(a0); // cross pr
831
832 G4ThreeVector c;
833
834 c.setX(std::cos(angle)*(a0.x())+std::sin(ang
835 c.setY(std::cos(angle)*(a0.y())+std::sin(ang
836 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang
837
838 G4ThreeVector c0 = c.unit();
839
840 return c0;
841 }
842
843 //....oooOO0OOooo........oooOO0OOooo........oo
844
845 G4ThreeVector G4LivermorePolarizedComptonModel
846 (const G4ThreeVector& gammaDirection, const G4
847 {
848 //
849 // The polarization of a photon is always pe
850 // Therefore this function removes those vec
851 // points in direction of gammaDirection
852 //
853 // Mathematically we search the projection o
854 // plains normal vector.
855 // The basic equation can be found in each g
856 // p = a - (a o n)/(n o n)*n
857
858 return gammaPolarization - gammaPolarization
859 }
860
861 //....oooOO0OOooo........oooOO0OOooo........oo
862
863 G4ThreeVector G4LivermorePolarizedComptonModel
864 G4double sinSqrTh,
865 G4double phi,
866 G4double costheta)
867 {
868 G4double rand1;
869 G4double rand2;
870 G4double cosPhi = std::cos(phi);
871 G4double sinPhi = std::sin(phi);
872 G4double sinTheta = std::sqrt(sinSqrTh);
873 G4double cosSqrPhi = cosPhi*cosPhi;
874 // G4double cossqrth = 1.-sinSqrTh;
875 // G4double sinsqrphi = sinPhi*sinPhi;
876 G4double normalisation = std::sqrt(1. - cosS
877
878 // Determination of Theta
879 G4double theta;
880
881 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
882 rand1 = G4UniformRand();
883 rand2 = G4UniformRand();
884
885 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi
886 {
887 if (rand2<0.5)
888 theta = pi/2.0;
889 else
890 theta = 3.0*pi/2.0;
891 }
892 else
893 {
894 if (rand2<0.5)
895 theta = 0;
896 else
897 theta = pi;
898 }
899 G4double cosBeta = std::cos(theta);
900 G4double sinBeta = std::sqrt(1-cosBeta*cosBe
901
902 G4ThreeVector gammaPolarization1;
903
904 G4double xParallel = normalisation*cosBeta;
905 G4double yParallel = -(sinSqrTh*cosPhi*sinPh
906 G4double zParallel = -(costheta*sinTheta*cos
907 G4double xPerpendicular = 0.;
908 G4double yPerpendicular = (costheta)*sinBeta
909 G4double zPerpendicular = -(sinTheta*sinPhi)
910
911 G4double xTotal = (xParallel + xPerpendicula
912 G4double yTotal = (yParallel + yPerpendicula
913 G4double zTotal = (zParallel + zPerpendicula
914
915 gammaPolarization1.setX(xTotal);
916 gammaPolarization1.setY(yTotal);
917 gammaPolarization1.setZ(zTotal);
918
919 return gammaPolarization1;
920 }
921
922 //....oooOO0OOooo........oooOO0OOooo........oo
923
924 void G4LivermorePolarizedComptonModel::SystemO
925 G4ThreeVector& direction1,
926 G4ThreeVector& polarization0,
927 G4ThreeVector& polarization1)
928 {
929 // direction0 is the original photon directi
930 // polarization0 is the original photon pola
931 // need to specify y axis in the real refere
932 G4ThreeVector Axis_Z0 = direction0.unit();
933 G4ThreeVector Axis_X0 = polarization0.unit()
934 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_
935
936 G4double direction_x = direction1.getX();
937 G4double direction_y = direction1.getY();
938 G4double direction_z = direction1.getZ();
939
940 direction1 = (direction_x*Axis_X0 + directio
941 G4double polarization_x = polarization1.getX
942 G4double polarization_y = polarization1.getY
943 G4double polarization_z = polarization1.getZ
944
945 polarization1 = (polarization_x*Axis_X0 + po
946
947 }
948
949
950 //....oooOO0OOooo........oooOO0OOooo........oo
951 //....oooOO0OOooo........oooOO0OOooo........oo
952
953 void
954 G4LivermorePolarizedComptonModel::InitialiseFo
955 G4int Z)
956 {
957 G4AutoLock l(&LivermorePolarizedComptonModel
958 if(!data[Z]) { ReadData(Z); }
959 l.unlock();
960 }
961