| 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: G4PolarizedIonisationModel
31 //
32 // Author: A.Schaelicke on base of Vlad
33 //
34 // Class Description:
35 // Implementation of energy loss and delta-e
36 // (including polarization effects)
37
38 #include "G4PolarizedIonisationModel.hh"
39
40 #include "G4ParticleChangeForLoss.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4PolarizationHelper.hh"
43 #include "G4PolarizationManager.hh"
44 #include "G4PolarizedIonisationBhabhaXS.hh"
45 #include "G4PolarizedIonisationMollerXS.hh"
46 #include "Randomize.hh"
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49 G4PolarizedIonisationModel::G4PolarizedIonisat
50 const G4ParticleDefinition* p, const G4Strin
51 : G4MollerBhabhaModel(p, nam)
52 , fCrossSectionCalculator(nullptr)
53 {
54 fBeamPolarization = G4StokesVector::ZERO;
55 fTargetPolarization = G4StokesVector::ZERO;
56
57 fPositronPolarization = G4StokesVector::ZERO
58 fElectronPolarization = G4StokesVector::ZERO
59
60 isElectron = (p == theElectron); // necessa
61 // G4Molle
62
63 if(!isElectron)
64 {
65 G4cout << " buildBhabha cross section " <<
66 fCrossSectionCalculator = new G4PolarizedI
67 }
68 else
69 {
70 G4cout << " buildMoller cross section " <<
71 fCrossSectionCalculator = new G4PolarizedI
72 }
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76 G4PolarizedIonisationModel::~G4PolarizedIonisa
77 {
78 if(fCrossSectionCalculator)
79 {
80 delete fCrossSectionCalculator;
81 }
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oo
85 G4double G4PolarizedIonisationModel::ComputeCr
86 const G4ParticleDefinition* pd, G4double kin
87 G4double emax)
88 {
89 G4double xs = G4MollerBhabhaModel::ComputeCr
90 pd, kinEnergy, cut, emax);
91 G4double factor = 1.;
92 if(xs != 0.)
93 {
94 G4double tmax = MaxSecondaryEnergy(pd, kin
95 tmax = std::min(emax, tmax);
96
97 if(std::fabs(cut / emax - 1.) < 1.e-10)
98 return xs;
99
100 if(cut < tmax)
101 {
102 G4double xmin = cut / kinEnergy;
103 G4double xmax = tmax / kinEnergy;
104 G4double gam = kinEnergy / electron
105 G4double crossPol = fCrossSectionCalcula
106 xmin, xmax, gam, fBeamPolarization, fT
107 G4double crossUnpol = fCrossSectionCalcu
108 xmin, xmax, gam, G4StokesVector::ZERO,
109 if(crossUnpol > 0.)
110 factor = crossPol / crossUnpol;
111 }
112 }
113 return xs * factor;
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117 void G4PolarizedIonisationModel::SampleSeconda
118 std::vector<G4DynamicParticle*>* vdp, const
119 const G4DynamicParticle* dp, G4double tmin,
120 {
121 // *** obtain and save target and beam polar
122 G4PolarizationManager* polarizationManager =
123 G4PolarizationManager::GetInstance();
124
125 const G4Track* aTrack = fParticleChange->Get
126
127 // obtain polarization of the beam
128 fBeamPolarization = G4StokesVector(dp->GetPo
129
130 // obtain polarization of the media
131 G4VPhysicalVolume* aPVolume = aTrack->Get
132 G4LogicalVolume* aLVolume = aPVolume->G
133 const G4bool targetIsPolarized = polarizatio
134 fTargetPolarization = polarizationManager->G
135
136 // transfer target polarization in interacti
137 if(targetIsPolarized)
138 fTargetPolarization.rotateUz(dp->GetMoment
139
140 G4double tmax = std::min(maxEnergy, MaxSecon
141 if(tmin >= tmax)
142 return;
143
144 G4double polL = fBeamPolarization.z() * fTar
145 polL = std::fabs(polL);
146 G4double polT = fBeamPolarization.x() * fTar
147 fBeamPolarization.y() * fTar
148 polT = std::fabs(polT);
149
150 G4double kineticEnergy = dp->GetKineticEnerg
151 G4double energy = kineticEnergy + ele
152 G4double totalMomentum =
153 std::sqrt(kineticEnergy * (energy + electr
154 G4double xmin = tmin / kineticEnergy;
155 G4double xmax = tmax / kineticEnergy;
156 G4double gam = energy / electron_mass_c2;
157 G4double gamma2 = gam * gam;
158 G4double gmo = gam - 1.;
159 G4double gmo2 = gmo * gmo;
160 G4double gmo3 = gmo2 * gmo;
161 G4double gpo = gam + 1.;
162 G4double gpo2 = gpo * gpo;
163 G4double gpo3 = gpo2 * gpo;
164 G4double x, y, q, grej, grej2;
165 G4double z = 0.;
166 G4double xs = 0., phi = 0.;
167 G4ThreeVector direction = dp->GetMomentumDir
168
169 //(Polarized) Moller (e-e-) scattering
170 if(isElectron)
171 {
172 // *** dice according to polarized cross s
173 G4double G = ((2.0 * gam - 1.0) / gamma2)
174 G4double H = (sqr(gam - 1.0) / gamma2) *
175 (1. + polT + polL * ((gam + 3
176
177 y = 1.0 - xmax;
178 grej = 1.0 - G * xmax + xmax * xmax * (H
179 grej2 = 1.0 - G * xmin + xmin * xmin * (H
180 if(grej2 > grej)
181 grej = grej2;
182 G4double prefM = gamma2 * classic_electr_r
183 (gmo2 * (gam + 1.0));
184 grej *= prefM;
185 do
186 {
187 q = G4UniformRand();
188 x = xmin * xmax / (xmin * (1.0 - q) + xm
189 if(fCrossSectionCalculator)
190 {
191 fCrossSectionCalculator->Initialize(x,
192 fT
193 xs = fCrossSectionCalculator->XSection
194
195 z = xs * sqr(x) * 4.;
196 if(grej < z)
197 {
198 G4ExceptionDescription ed;
199 ed << "WARNING : error in Moller rej
200 << " z = " << z << " grej=" << gr
201 G4Exception("G4PolarizedIonisationMo
202 JustWarning, ed);
203 }
204 }
205 else
206 {
207 G4cout << "No calculator in Moller sca
208 }
209 // Loop checking, 03-Aug-2015, Vladimir
210 } while(grej * G4UniformRand() > z);
211 // Bhabha (e+e-) scattering
212 }
213 else
214 {
215 // *** dice according to polarized cross s
216 y = xmax * xmax;
217 grej = 0.;
218 grej += y * y * gmo3 * (1. + (polL + polT)
219 grej += -2. * xmin * xmin * xmin * gam * g
220 (1. - (polL + polT) * (gam + 3.) /
221 grej += y * y * gmo * (3. * gamma2 + 6. *
222 (1. + (polL * (3. * gam + 1.) * (g
223 polT * ((gam + 2.) * gamma2
224 (gmo * (3. * gam * (gam +
225 grej /= gpo3;
226 grej += -xmin * (2. * gamma2 + 4. * gam +
227 (1. - gam * (polL * (2. * gam + 1.
228 (2. * gam * (gam + 2.) + 1
229 gpo2;
230 grej += gamma2 / (gamma2 - 1.);
231 G4double prefB =
232 classic_electr_radius * classic_electr_r
233 grej *= prefB;
234
235 do
236 {
237 q = G4UniformRand();
238 x = xmin * xmax / (xmin * (1.0 - q) + xm
239 if(fCrossSectionCalculator)
240 {
241 fCrossSectionCalculator->Initialize(x,
242 fT
243 xs = fCrossSectionCalculator->XSection
244
245 z = xs * sqr(x) * 4.;
246 }
247 else
248 {
249 G4cout << "No calculator in Bhabha sca
250 }
251
252 if(z > grej)
253 {
254 G4ExceptionDescription ed;
255 ed << "G4PolarizedIonisationModel::Sam
256 << "Majorant " << grej << " < " <<
257 << " e+e- (Bhabha) scattering"
258 << " at KinEnergy " << kineticEnerg
259 G4Exception("G4PolarizedIonisationMode
260 JustWarning, ed);
261 }
262 // Loop checking, 03-Aug-2015, Vladimir
263 } while(grej * G4UniformRand() > z);
264 }
265
266 // polar asymmetries (due to transverse pola
267 if(fCrossSectionCalculator)
268 {
269 grej = xs * 2.;
270 do
271 {
272 phi = twopi * G4UniformRand();
273 fCrossSectionCalculator->Initialize(x, g
274 fTar
275 xs = fCrossSectionCalculator->XSection(G
276 G
277 if(xs > grej)
278 {
279 if(isElectron)
280 {
281 G4ExceptionDescription ed;
282 ed << "Majorant " << grej << " < " <
283 << "\n"
284 << " e-e- (Moller) scattering\n"
285 << "PHI DICING\n";
286 G4Exception("G4PolarizedIonisationMo
287 JustWarning, ed);
288 }
289 else
290 {
291 G4ExceptionDescription ed;
292 ed << "Majorant " << grej << " < " <
293 << "\n"
294 << " e+e- (Bhabha) scattering\n"
295 << "PHI DICING\n";
296 G4Exception("G4PolarizedIonisationMo
297 JustWarning, ed);
298 }
299 }
300 // Loop checking, 03-Aug-2015, Vladimir
301 } while(grej * G4UniformRand() > xs);
302 }
303
304 // fix kinematics of delta electron
305 G4double deltaKinEnergy = x * kineticEnergy;
306 G4double deltaMomentum =
307 std::sqrt(deltaKinEnergy * (deltaKinEnergy
308 G4double cost = deltaKinEnergy * (energy + e
309 (deltaMomentum * totalMoment
310 G4double sint = 1.0 - cost * cost;
311 if(sint > 0.0)
312 sint = std::sqrt(sint);
313
314 G4ThreeVector deltaDirection(-sint * std::co
315 cost);
316 deltaDirection.rotateUz(direction);
317
318 // primary change
319 kineticEnergy -= deltaKinEnergy;
320 fParticleChange->SetProposedKineticEnergy(ki
321
322 if(kineticEnergy > DBL_MIN)
323 {
324 G4ThreeVector dir =
325 totalMomentum * direction - deltaMomentu
326 direction = dir.unit();
327 fParticleChange->SetProposedMomentumDirect
328 }
329
330 // create G4DynamicParticle object for delta
331 G4DynamicParticle* delta =
332 new G4DynamicParticle(theElectron, deltaDi
333 vdp->push_back(delta);
334
335 // get interaction frame
336 G4ThreeVector nInteractionFrame =
337 G4PolarizationHelper::GetFrame(direction,
338
339 if(fCrossSectionCalculator)
340 {
341 // calculate mean final state polarization
342 fBeamPolarization.InvRotateAz(nInteraction
343 fTargetPolarization.InvRotateAz(nInteracti
344 fCrossSectionCalculator->Initialize(x, gam
345 fTarge
346
347 // electron/positron
348 fPositronPolarization = fCrossSectionCalcu
349 fPositronPolarization.RotateAz(nInteractio
350
351 fParticleChange->ProposePolarization(fPosi
352
353 // electron
354 fElectronPolarization = fCrossSectionCalcu
355 fElectronPolarization.RotateAz(nInteractio
356 delta->SetPolarization(fElectronPolarizati
357 fElectronPolarizati
358 }
359 else
360 {
361 fPositronPolarization = G4StokesVector::ZE
362 fElectronPolarization = G4StokesVector::ZE
363 }
364 }
365