| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // -------------------------------------------------------------------
27 //
28 // Geant4 Class file
29 //
30 // File name: G4PolarizedIonisationModel
31 //
32 // Author: A.Schaelicke on base of Vladimir Ivanchenko code
33 //
34 // Class Description:
35 // Implementation of energy loss and delta-electron production by e+/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........oooOO0OOooo........oooOO0OOooo....
49 G4PolarizedIonisationModel::G4PolarizedIonisationModel(
50 const G4ParticleDefinition* p, const G4String& nam)
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); // necessary due to wrong order in
61 // G4MollerBhabhaModel constructor!
62
63 if(!isElectron)
64 {
65 G4cout << " buildBhabha cross section " << isElectron << G4endl;
66 fCrossSectionCalculator = new G4PolarizedIonisationBhabhaXS();
67 }
68 else
69 {
70 G4cout << " buildMoller cross section " << isElectron << G4endl;
71 fCrossSectionCalculator = new G4PolarizedIonisationMollerXS();
72 }
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 G4PolarizedIonisationModel::~G4PolarizedIonisationModel()
77 {
78 if(fCrossSectionCalculator)
79 {
80 delete fCrossSectionCalculator;
81 }
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 G4double G4PolarizedIonisationModel::ComputeCrossSectionPerElectron(
86 const G4ParticleDefinition* pd, G4double kinEnergy, G4double cut,
87 G4double emax)
88 {
89 G4double xs = G4MollerBhabhaModel::ComputeCrossSectionPerElectron(
90 pd, kinEnergy, cut, emax);
91 G4double factor = 1.;
92 if(xs != 0.)
93 {
94 G4double tmax = MaxSecondaryEnergy(pd, kinEnergy);
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_mass_c2 + 1.0;
105 G4double crossPol = fCrossSectionCalculator->TotalXSection(
106 xmin, xmax, gam, fBeamPolarization, fTargetPolarization);
107 G4double crossUnpol = fCrossSectionCalculator->TotalXSection(
108 xmin, xmax, gam, G4StokesVector::ZERO, G4StokesVector::ZERO);
109 if(crossUnpol > 0.)
110 factor = crossPol / crossUnpol;
111 }
112 }
113 return xs * factor;
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 void G4PolarizedIonisationModel::SampleSecondaries(
118 std::vector<G4DynamicParticle*>* vdp, const G4MaterialCutsCouple*,
119 const G4DynamicParticle* dp, G4double tmin, G4double maxEnergy)
120 {
121 // *** obtain and save target and beam polarization ***
122 G4PolarizationManager* polarizationManager =
123 G4PolarizationManager::GetInstance();
124
125 const G4Track* aTrack = fParticleChange->GetCurrentTrack();
126
127 // obtain polarization of the beam
128 fBeamPolarization = G4StokesVector(dp->GetPolarization());
129
130 // obtain polarization of the media
131 G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
132 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
133 const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
134 fTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
135
136 // transfer target polarization in interaction frame
137 if(targetIsPolarized)
138 fTargetPolarization.rotateUz(dp->GetMomentumDirection());
139
140 G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp));
141 if(tmin >= tmax)
142 return;
143
144 G4double polL = fBeamPolarization.z() * fTargetPolarization.z();
145 polL = std::fabs(polL);
146 G4double polT = fBeamPolarization.x() * fTargetPolarization.x() +
147 fBeamPolarization.y() * fTargetPolarization.y();
148 polT = std::fabs(polT);
149
150 G4double kineticEnergy = dp->GetKineticEnergy();
151 G4double energy = kineticEnergy + electron_mass_c2;
152 G4double totalMomentum =
153 std::sqrt(kineticEnergy * (energy + electron_mass_c2));
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->GetMomentumDirection();
168
169 //(Polarized) Moller (e-e-) scattering
170 if(isElectron)
171 {
172 // *** dice according to polarized cross section
173 G4double G = ((2.0 * gam - 1.0) / gamma2) * (1. - polT - polL * gam);
174 G4double H = (sqr(gam - 1.0) / gamma2) *
175 (1. + polT + polL * ((gam + 3.) / (gam - 1.)));
176
177 y = 1.0 - xmax;
178 grej = 1.0 - G * xmax + xmax * xmax * (H + (1.0 - G * y) / (y * y));
179 grej2 = 1.0 - G * xmin + xmin * xmin * (H + (1.0 - G * y) / (y * y));
180 if(grej2 > grej)
181 grej = grej2;
182 G4double prefM = gamma2 * classic_electr_radius * classic_electr_radius /
183 (gmo2 * (gam + 1.0));
184 grej *= prefM;
185 do
186 {
187 q = G4UniformRand();
188 x = xmin * xmax / (xmin * (1.0 - q) + xmax * q);
189 if(fCrossSectionCalculator)
190 {
191 fCrossSectionCalculator->Initialize(x, gam, phi, fBeamPolarization,
192 fTargetPolarization, 1);
193 xs = fCrossSectionCalculator->XSection(G4StokesVector::ZERO,
194 G4StokesVector::ZERO);
195 z = xs * sqr(x) * 4.;
196 if(grej < z)
197 {
198 G4ExceptionDescription ed;
199 ed << "WARNING : error in Moller rejection routine! \n"
200 << " z = " << z << " grej=" << grej << "\n";
201 G4Exception("G4PolarizedIonisationModel::SampleSecondaries", "pol019",
202 JustWarning, ed);
203 }
204 }
205 else
206 {
207 G4cout << "No calculator in Moller scattering" << G4endl;
208 }
209 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
210 } while(grej * G4UniformRand() > z);
211 // Bhabha (e+e-) scattering
212 }
213 else
214 {
215 // *** dice according to polarized cross section
216 y = xmax * xmax;
217 grej = 0.;
218 grej += y * y * gmo3 * (1. + (polL + polT) * (gam + 3.) / gmo);
219 grej += -2. * xmin * xmin * xmin * gam * gmo2 *
220 (1. - (polL + polT) * (gam + 3.) / gmo);
221 grej += y * y * gmo * (3. * gamma2 + 6. * gam + 4.) *
222 (1. + (polL * (3. * gam + 1.) * (gamma2 + gam + 1.) +
223 polT * ((gam + 2.) * gamma2 + 1.)) /
224 (gmo * (3. * gam * (gam + 2.) + 4.)));
225 grej /= gpo3;
226 grej += -xmin * (2. * gamma2 + 4. * gam + 1.) *
227 (1. - gam * (polL * (2. * gam + 1.) + polT) /
228 (2. * gam * (gam + 2.) + 1.)) /
229 gpo2;
230 grej += gamma2 / (gamma2 - 1.);
231 G4double prefB =
232 classic_electr_radius * classic_electr_radius / (gam - 1.0);
233 grej *= prefB;
234
235 do
236 {
237 q = G4UniformRand();
238 x = xmin * xmax / (xmin * (1.0 - q) + xmax * q);
239 if(fCrossSectionCalculator)
240 {
241 fCrossSectionCalculator->Initialize(x, gam, phi, fBeamPolarization,
242 fTargetPolarization, 1);
243 xs = fCrossSectionCalculator->XSection(G4StokesVector::ZERO,
244 G4StokesVector::ZERO);
245 z = xs * sqr(x) * 4.;
246 }
247 else
248 {
249 G4cout << "No calculator in Bhabha scattering" << G4endl;
250 }
251
252 if(z > grej)
253 {
254 G4ExceptionDescription ed;
255 ed << "G4PolarizedIonisationModel::SampleSecondaries Warning!\n "
256 << "Majorant " << grej << " < " << z << " for x= " << x << G4endl
257 << " e+e- (Bhabha) scattering"
258 << " at KinEnergy " << kineticEnergy << G4endl;
259 G4Exception("G4PolarizedIonisationModel::SampleSecondaries", "pol020",
260 JustWarning, ed);
261 }
262 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
263 } while(grej * G4UniformRand() > z);
264 }
265
266 // polar asymmetries (due to transverse polarizations)
267 if(fCrossSectionCalculator)
268 {
269 grej = xs * 2.;
270 do
271 {
272 phi = twopi * G4UniformRand();
273 fCrossSectionCalculator->Initialize(x, gam, phi, fBeamPolarization,
274 fTargetPolarization, 1);
275 xs = fCrossSectionCalculator->XSection(G4StokesVector::ZERO,
276 G4StokesVector::ZERO);
277 if(xs > grej)
278 {
279 if(isElectron)
280 {
281 G4ExceptionDescription ed;
282 ed << "Majorant " << grej << " < " << xs << " for phi= " << phi
283 << "\n"
284 << " e-e- (Moller) scattering\n"
285 << "PHI DICING\n";
286 G4Exception("G4PolarizedIonisationModel::SampleSecondaries", "pol021",
287 JustWarning, ed);
288 }
289 else
290 {
291 G4ExceptionDescription ed;
292 ed << "Majorant " << grej << " < " << xs << " for phi= " << phi
293 << "\n"
294 << " e+e- (Bhabha) scattering\n"
295 << "PHI DICING\n";
296 G4Exception("G4PolarizedIonisationModel::SampleSecondaries", "pol022",
297 JustWarning, ed);
298 }
299 }
300 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
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 + 2.0 * electron_mass_c2));
308 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
309 (deltaMomentum * totalMomentum);
310 G4double sint = 1.0 - cost * cost;
311 if(sint > 0.0)
312 sint = std::sqrt(sint);
313
314 G4ThreeVector deltaDirection(-sint * std::cos(phi), -sint * std::sin(phi),
315 cost);
316 deltaDirection.rotateUz(direction);
317
318 // primary change
319 kineticEnergy -= deltaKinEnergy;
320 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
321
322 if(kineticEnergy > DBL_MIN)
323 {
324 G4ThreeVector dir =
325 totalMomentum * direction - deltaMomentum * deltaDirection;
326 direction = dir.unit();
327 fParticleChange->SetProposedMomentumDirection(direction);
328 }
329
330 // create G4DynamicParticle object for delta ray
331 G4DynamicParticle* delta =
332 new G4DynamicParticle(theElectron, deltaDirection, deltaKinEnergy);
333 vdp->push_back(delta);
334
335 // get interaction frame
336 G4ThreeVector nInteractionFrame =
337 G4PolarizationHelper::GetFrame(direction, deltaDirection);
338
339 if(fCrossSectionCalculator)
340 {
341 // calculate mean final state polarizations
342 fBeamPolarization.InvRotateAz(nInteractionFrame, direction);
343 fTargetPolarization.InvRotateAz(nInteractionFrame, direction);
344 fCrossSectionCalculator->Initialize(x, gam, phi, fBeamPolarization,
345 fTargetPolarization, 2);
346
347 // electron/positron
348 fPositronPolarization = fCrossSectionCalculator->GetPol2();
349 fPositronPolarization.RotateAz(nInteractionFrame, direction);
350
351 fParticleChange->ProposePolarization(fPositronPolarization);
352
353 // electron
354 fElectronPolarization = fCrossSectionCalculator->GetPol3();
355 fElectronPolarization.RotateAz(nInteractionFrame, deltaDirection);
356 delta->SetPolarization(fElectronPolarization.x(), fElectronPolarization.y(),
357 fElectronPolarization.z());
358 }
359 else
360 {
361 fPositronPolarization = G4StokesVector::ZERO;
362 fElectronPolarization = G4StokesVector::ZERO;
363 }
364 }
365