| 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 // Geant4 Class file
29 //
30 // File name: G4PolarizedComptonModel
31 //
32 // Author: Andreas Schaelicke
33
34 #include "G4PolarizedComptonModel.hh"
35
36 #include "G4Exp.hh"
37 #include "G4Log.hh"
38 #include "G4ParticleChangeForGamma.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4PolarizationManager.hh"
41 #include "G4PolarizationHelper.hh"
42 #include "G4PolarizedComptonXS.hh"
43 #include "G4StokesVector.hh"
44 #include "G4SystemOfUnits.hh"
45
46 //....oooOO0OOooo........oooOO0OOooo........oo
47 G4PolarizedComptonModel::G4PolarizedComptonMod
48
49 : G4KleinNishinaCompton(nullptr, nam)
50 , fVerboseLevel(0)
51 {
52 fCrossSectionCalculator = new G4PolarizedCom
53 fBeamPolarization = G4StokesVector::ZE
54 fTargetPolarization = G4StokesVector::ZE
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58 G4PolarizedComptonModel::~G4PolarizedComptonMo
59 {
60 delete fCrossSectionCalculator;
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64 G4double G4PolarizedComptonModel::ComputeAsymm
65
66 {
67 G4double asymmetry = 0.0;
68
69 G4double k0 = gammaEnergy / electron_mass_c2
70 G4double k1 = 1. + 2. * k0;
71
72 asymmetry = -k0;
73 asymmetry *=
74 (k0 + 1.) * sqr(k1) * G4Log(k1) - 2. * k0
75 asymmetry /= ((k0 - 2.) * k0 - 2.) * sqr(k1)
76 2. * k0 * (k0 * (k0 + 1.) * (k0
77
78 if(asymmetry > 1.)
79 {
80 G4ExceptionDescription ed;
81 ed << "ERROR in G4PolarizedComptonModel::C
82 << " asymmetry = " << asymmetry << "\n"
83 G4Exception("G4PolarizedComptonModel::Comp
84 JustWarning, ed);
85 }
86
87 return asymmetry;
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oo
91 G4double G4PolarizedComptonModel::ComputeCross
92 const G4ParticleDefinition* pd, G4double kin
93 G4double cut, G4double emax)
94 {
95 G4double xs = G4KleinNishinaCompton::Compute
96 pd, kinEnergy, Z, A, cut, emax);
97 G4double polzz = fBeamPolarization.p3() * fT
98 if(polzz > 0.0)
99 {
100 G4double asym = ComputeAsymmetryPerAtom(ki
101 xs *= (1. + polzz * asym);
102 }
103 return xs;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107 void G4PolarizedComptonModel::SampleSecondarie
108 std::vector<G4DynamicParticle*>* fvect, cons
109 const G4DynamicParticle* aDynamicGamma, G4do
110 {
111 // do nothing below the threshold
112 if(aDynamicGamma->GetKineticEnergy() <= LowE
113 {
114 return;
115 }
116
117 const G4Track* aTrack = fParticleChang
118 G4VPhysicalVolume* aPVolume = aTrack->GetVol
119 G4LogicalVolume* aLVolume = aPVolume->GetL
120
121 if(fVerboseLevel >= 1)
122 {
123 G4cout << "G4PolarizedComptonModel::Sample
124 << aLVolume->GetName() << G4endl;
125 }
126 G4PolarizationManager* polarizationManager =
127 G4PolarizationManager::GetInstance();
128
129 // obtain polarization of the beam
130 fBeamPolarization = G4StokesVector(aDynamicG
131 fBeamPolarization.SetPhoton();
132
133 // obtain polarization of the media
134 G4bool targetIsPolarized = polarizationManag
135 fTargetPolarization = polarizationManager->G
136
137 // if beam is linear polarized or target is
138 // determine the angle to x-axis
139 // (assumes same PRF as in the polarization
140 G4ThreeVector gamDirection0 = aDynamicGamma-
141
142 // transfer fTargetPolarization
143 // into the gamma frame (problem electron is
144 if(targetIsPolarized)
145 {
146 fTargetPolarization.rotateUz(gamDirection0
147 }
148 // The scattered gamma energy is sampled acc
149 // Klein - Nishina formula.
150 // The random number techniques of Butcher &
151 // (Nuc Phys 20(1960),15).
152 // Note : Effects due to binding of atomic e
153
154 G4double gamEnergy0 = aDynamicGamma->GetKine
155 G4double E0_m = gamEnergy0 / electron_
156
157 // sample the energy rate of the scattered g
158 G4double epsilon, sint2;
159 G4double onecost = 0.0;
160 G4double Phi = 0.0;
161 G4double greject = 1.0;
162 G4double cosTeta = 1.0;
163 G4double sinTeta = 0.0;
164
165 G4double eps0 = 1. / (1. + 2. * E0_m);
166 G4double epsilon0sq = eps0 * eps0;
167 G4double alpha1 = -G4Log(eps0);
168 G4double alpha2 = alpha1 + 0.5 * (1. - e
169
170 G4double polarization = fBeamPolarization.p3
171
172 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
173 G4int nloop = 0;
174 G4bool end = fals
175
176 G4double rndm[3];
177
178 do
179 {
180 do
181 {
182 ++nloop;
183 // false interaction if too many iterati
184 if(nloop > fLoopLim)
185 {
186 PrintWarning(aDynamicGamma, nloop, gre
187 "too many iterations");
188 return;
189 }
190
191 // 3 random numbers to sample scattering
192 rndmEngineMod->flatArray(3, rndm);
193
194 if(alpha1 > alpha2 * rndm[0])
195 {
196 epsilon = G4Exp(-alpha1 * rndm[1]);
197 }
198 else
199 {
200 epsilon = std::sqrt(epsilon0sq + (1. -
201 }
202
203 onecost = (1. - epsilon) / (epsilon * E0
204 sint2 = onecost * (2. - onecost);
205
206 G4double gdiced = 2. * (1. / epsilon + e
207 G4double gdist = 1. / epsilon + epsilon
208 polarization * (1. / ep
209
210 greject = gdist / gdiced;
211
212 if(greject > 1.0)
213 {
214 PrintWarning(aDynamicGamma, nloop, gre
215 "theta majoranta wrong");
216 }
217 // Loop checking, 03-Aug-2015, Vladimir
218 } while(greject < rndm[2]);
219
220 // assuming phi loop successful
221 end = true;
222
223 // scattered gamma angles. ( Z - axis alon
224 cosTeta = 1. - onecost;
225 sinTeta = std::sqrt(sint2);
226 do
227 {
228 ++nloop;
229
230 // 2 random numbers to sample scattering
231 rndmEngineMod->flatArray(2, rndm);
232
233 // false interaction if too many iterati
234 Phi = twopi * rndm[0];
235 if(nloop > fLoopLim)
236 {
237 PrintWarning(aDynamicGamma, nloop, gre
238 "too many iterations");
239 return;
240 }
241
242 G4double gdiced = 1. / epsilon + epsilon
243 std::abs(fBeamPolariza
244 (std::abs((1. / epsi
245 fTargetPol
246 (1. - epsilon) * si
247 (std::sqrt(sqr(fT
248 sqr(fT
249 sint2 * (std::sqrt(sqr
250 sqr
251
252 G4double gdist =
253 1. / epsilon + epsilon - sint2 +
254 fBeamPolarization.p3() *
255 ((1. / epsilon - epsilon) * cosTeta
256 (1. - epsilon) * sinTeta *
257 (std::cos(Phi) * fTargetPolarizat
258 std::sin(Phi) * fTargetPolarizat
259 sint2 * (std::cos(2. * Phi) * fBeamPol
260 std::sin(2. * Phi) * fBeamPol
261 greject = gdist / gdiced;
262
263 if(greject > 1.0)
264 {
265 PrintWarning(aDynamicGamma, nloop, gre
266 "phi majoranta wrong");
267 }
268
269 if(greject < 1.e-3)
270 {
271 PrintWarning(aDynamicGamma, nloop, gre
272 "phi loop ineffective");
273 // restart theta loop
274 end = false;
275 break;
276 }
277
278 // Loop checking, 03-Aug-2015, Vladimir
279 } while(greject < rndm[1]);
280 } while(!end);
281 G4double dirx = sinTeta * std::cos(Phi);
282 G4double diry = sinTeta * std::sin(Phi);
283 G4double dirz = cosTeta;
284
285 // update G4VParticleChange for the scattere
286 G4ThreeVector gamDirection1(dirx, diry, dirz
287 gamDirection1.rotateUz(gamDirection0);
288 G4double gamEnergy1 = epsilon * gamEnergy0;
289
290 G4double edep = 0.0;
291 if(gamEnergy1 > lowestSecondaryEnergy)
292 {
293 fParticleChange->ProposeMomentumDirection(
294 fParticleChange->SetProposedKineticEnergy(
295 }
296 else
297 {
298 fParticleChange->ProposeTrackStatus(fStopA
299 fParticleChange->SetProposedKineticEnergy(
300 edep = gamEnergy1;
301 }
302
303 // calculate Stokes vector of final state ph
304 G4ThreeVector nInteractionFrame =
305 G4PolarizationHelper::GetFrame(gamDirectio
306
307 // transfer fBeamPolarization and fTargetPol
308 // into the interaction frame (note electron
309 if(fVerboseLevel >= 1)
310 {
311 G4cout << "===============================
312 G4cout << " nInteractionFrame = " << nInte
313 G4cout << " GammaDirection0 = " << gamDire
314 G4cout << " gammaPolarization = " << fBeam
315 G4cout << " electronPolarization = " << fT
316 }
317
318 fBeamPolarization.InvRotateAz(nInteractionFr
319 fTargetPolarization.InvRotateAz(nInteraction
320
321 if(fVerboseLevel >= 1)
322 {
323 G4cout << "-------------------------------
324 G4cout << " gammaPolarization = " << fBeam
325 G4cout << " electronPolarization = " << fT
326 G4cout << "-------------------------------
327 }
328
329 // initialize the polarization transfer matr
330 fCrossSectionCalculator->Initialize(epsilon,
331 fTargetP
332
333 if(gamEnergy1 > lowestSecondaryEnergy)
334 {
335 // in interaction frame
336 // calculate polarization transfer to the
337 fFinalGammaPolarization = fCrossSectionCal
338 if(fVerboseLevel >= 1)
339 {
340 G4cout << " gammaPolarization1 = " << fF
341 }
342 fFinalGammaPolarization.SetPhoton();
343
344 // translate polarization into particle re
345 fFinalGammaPolarization.RotateAz(nInteract
346 if(fFinalGammaPolarization.mag() > 1. + 1.
347 {
348 G4ExceptionDescription ed;
349 ed << "ERROR in Polarizaed Compton Scatt
350 ed << "Polarization of final photon more
351 ed << fFinalGammaPolarization
352 << " mag = " << fFinalGammaPolarizati
353 G4Exception("G4PolarizedComptonModel::Sa
354 FatalException, ed);
355 }
356 // store polarization vector
357 fParticleChange->ProposePolarization(fFina
358 if(fVerboseLevel >= 1)
359 {
360 G4cout << " gammaPolarization1 = " << fF
361 G4cout << " GammaDirection1 = " << gamDi
362 }
363 }
364
365 // kinematic of the scattered electron
366 G4double eKinEnergy = gamEnergy0 - gamEnergy
367
368 if(eKinEnergy > lowestSecondaryEnergy)
369 {
370 G4ThreeVector eDirection =
371 gamEnergy0 * gamDirection0 - gamEnergy1
372 eDirection = eDirection.unit();
373
374 finalElectronPolarization = fCrossSectionC
375 if(fVerboseLevel >= 1)
376 {
377 G4cout << " electronPolarization1 = " <<
378 << G4endl;
379 }
380 // transfer into particle reference frame
381 finalElectronPolarization.RotateAz(nIntera
382 if(fVerboseLevel >= 1)
383 {
384 G4cout << " electronPolarization1 = " <<
385 << G4endl << " ElecDirection = "
386 }
387
388 // create G4DynamicParticle object for the
389 G4DynamicParticle* aElectron =
390 new G4DynamicParticle(theElectron, eDire
391 // store polarization vector
392 if(finalElectronPolarization.mag() > 1. +
393 {
394 G4ExceptionDescription ed;
395 ed << "ERROR in Polarized Compton Scatte
396 ed << "Polarization of final electron mo
397 ed << finalElectronPolarization
398 << " mag = " << finalElectronPolariza
399 G4Exception("G4PolarizedComptonModel::Sa
400 FatalException, ed);
401 }
402 aElectron->SetPolarization(finalElectronPo
403 finalElectronPo
404 finalElectronPo
405 fvect->push_back(aElectron);
406 }
407 else
408 {
409 edep += eKinEnergy;
410 }
411 // energy balance
412 if(edep > 0.0)
413 {
414 fParticleChange->ProposeLocalEnergyDeposit
415 }
416 }
417
418 //....oooOO0OOooo........oooOO0OOooo........oo
419 void G4PolarizedComptonModel::PrintWarning(con
420 G4i
421 G4d
422 con
423 {
424 G4ExceptionDescription ed;
425 ed << "Problem of scattering sampling: " <<
426 << "Niter= " << nloop << " grej= " << gre
427 << " cos(theta)= " << 1.0 - onecos << " p
428 << "Gamma E(MeV)= " << dp->GetKineticEner
429 << " dir= " << dp->GetMomentumDirection()
430 << " pol= " << dp->GetPolarization();
431 G4Exception("G4PolarizedComptonModel::Sample
432 JustWarning, ed, "");
433 }
434