120{
121
124
126
127
129
130
135
136
137 if(targetIsPolarized)
139
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
153 std::sqrt(kineticEnergy * (energy + electron_mass_c2));
154 G4double xmin = tmin / kineticEnergy;
155 G4double xmax = tmax / kineticEnergy;
168
169
171 {
172
173 G4double G = ((2.0 *
gam - 1.0) / gamma2) * (1. - polT - polL * gam);
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 {
188 x = xmin * xmax / (xmin * (1.0 - q) + xmax * q);
189 if(fCrossSectionCalculator)
190 {
191 fCrossSectionCalculator->
Initialize(x, gam, phi, fBeamPolarization,
192 fTargetPolarization, 1);
195 z = xs *
sqr(x) * 4.;
196 if(grej < z)
197 {
199 ed << "WARNING : error in Moller rejection routine! \n"
200 << " z = " << z << " grej=" << grej << "\n";
201 G4Exception(
"G4PolarizedIonisationModel::SampleSecondaries",
"pol019",
203 }
204 }
205 else
206 {
207 G4cout <<
"No calculator in Moller scattering" <<
G4endl;
208 }
209
211
212 }
213 else
214 {
215
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.);
232 classic_electr_radius * classic_electr_radius / (
gam - 1.0);
233 grej *= prefB;
234
235 do
236 {
238 x = xmin * xmax / (xmin * (1.0 - q) + xmax * q);
239 if(fCrossSectionCalculator)
240 {
241 fCrossSectionCalculator->
Initialize(x, gam, phi, fBeamPolarization,
242 fTargetPolarization, 1);
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 {
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",
261 }
262
264 }
265
266
267 if(fCrossSectionCalculator)
268 {
269 grej = xs * 2.;
270 do
271 {
273 fCrossSectionCalculator->
Initialize(x, gam, phi, fBeamPolarization,
274 fTargetPolarization, 1);
277 if(xs > grej)
278 {
280 {
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",
288 }
289 else
290 {
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",
298 }
299 }
300
302 }
303
304
305 G4double deltaKinEnergy = x * kineticEnergy;
307 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0 * electron_mass_c2));
309 (deltaMomentum * totalMomentum);
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
319 kineticEnergy -= deltaKinEnergy;
321
323 {
325 totalMomentum * direction - deltaMomentum * deltaDirection;
326 direction = dir.
unit();
328 }
329
330
333 vdp->push_back(delta);
334
335
338
339 if(fCrossSectionCalculator)
340 {
341
342 fBeamPolarization.
InvRotateAz(nInteractionFrame, direction);
343 fTargetPolarization.
InvRotateAz(nInteractionFrame, direction);
344 fCrossSectionCalculator->
Initialize(x, gam, phi, fBeamPolarization,
345 fTargetPolarization, 2);
346
347
348 fPositronPolarization = fCrossSectionCalculator->
GetPol2();
349 fPositronPolarization.
RotateAz(nInteractionFrame, direction);
350
352
353
354 fElectronPolarization = fCrossSectionCalculator->
GetPol3();
355 fElectronPolarization.
RotateAz(nInteractionFrame, deltaDirection);
357 fElectronPolarization.
z());
358 }
359 else
360 {
363 }
364}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
Hep3Vector & rotateUz(const Hep3Vector &)
void SetPolarization(const G4ThreeVector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetPolarization() const
G4ParticleChangeForLoss * fParticleChange
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposePolarization(const G4ThreeVector &dir)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
bool IsPolarized(G4LogicalVolume *lVol) const
const G4StokesVector GetVolumePolarization(G4LogicalVolume *lVol) const
static G4PolarizationManager * GetInstance()
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
G4VPhysicalVolume * GetVolume() const
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
const G4Track * GetCurrentTrack() const
G4LogicalVolume * GetLogicalVolume() const
virtual void Initialize(G4double, G4double, G4double, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)=0
virtual G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3)=0
virtual G4StokesVector GetPol2()
virtual G4StokesVector GetPol3()
G4double energy(const ThreeVector &p, const G4double m)