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