144{
145
147
149
150
152
153
158
159
160 if (targetIsPolarized)
162
163
164
165
167 if(tmin >= tmax) return;
168
169
170 G4double polL = theBeamPolarization.
z()*theTargetPolarization.
z();
171 polL=std::fabs(polL);
172 G4double polT = theBeamPolarization.
x()*theTargetPolarization.
x() +
173 theBeamPolarization.
y()*theTargetPolarization.
y();
174 polT=std::fabs(polT);
175
178 G4double totalMomentum = std::sqrt(kineticEnergy*(energy + electron_mass_c2));
193
194
196
198 G4double H = (
sqr(gam - 1.0)/gamma2)*(1. + polT + polL*((gam + 3.)/(gam - 1.)));
199
200 y = 1.0 - xmax;
201 grej = 1.0 - G*xmax + xmax*xmax*(H + (1.0 - G*y)/(y*y));
202 grej2 = 1.0 - G*xmin + xmin*xmin*(H + (1.0 - G*y)/(y*y));
203 if (grej2 > grej) grej = grej2;
204 G4double prefM = gamma2*classic_electr_radius*classic_electr_radius/(gmo2*(
gam + 1.0));
205 grej *= prefM;
206 do {
208 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
209 if (crossSectionCalculator) {
210 crossSectionCalculator->
Initialize(x,gam,phi,theBeamPolarization,
211 theTargetPolarization,1);
215 if (grej < z) {
216 G4cout<<
"WARNING : error in Moller rejection routine! \n"
217 <<" z = "<<z<<" grej="<<grej<<"\n";
218 }
219 } else {
221 }
222
224
225 } else {
226
227 y = xmax*xmax;
228 grej = 0.;
229 grej += y*y*gmo3*(1. + (polL + polT)*(gam + 3.)/gmo);
230 grej += -2.*xmin*xmin*xmin*
gam*gmo2*(1. - (polL + polT)*(gam + 3.)/gmo);
231 grej += y*y*gmo*(3.*gamma2 + 6.*
gam + 4.)*(1. + (polL*(3.*gam + 1.)*(gamma2 +
gam + 1.) + polT*((gam + 2.)*gamma2 + 1.))/(gmo*(3.*
gam*(
gam + 2.) + 4.)));
232 grej /= gpo3;
233 grej += -xmin*(2.*gamma2 + 4.*
gam + 1.)*(1. - gam*(polL*(2.*gam + 1.) + polT)/(2.*
gam*(
gam + 2.) + 1.))/gpo2;
234 grej += gamma2/(gamma2 - 1.);
235 G4double prefB = classic_electr_radius*classic_electr_radius/(
gam - 1.0);
236 grej *= prefB;
237
238 do {
240 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
241 if (crossSectionCalculator) {
242 crossSectionCalculator->
Initialize(x,gam,phi,theBeamPolarization,
243 theTargetPolarization,1);
247 } else {
249 }
250
251 if(z > grej) {
252 G4cout<<
"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<
G4endl;
253 G4cout <<
"G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<
G4endl
254 << "Majorant " << grej << " < "
255 << z <<
" for x= " << x<<
G4endl
256 <<
" e+e- (Bhabha) scattering"<<
" at KinEnergy "<<kineticEnergy<<
G4endl;
257 G4cout<<
"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<
G4endl;
258 }
259
261 }
262
263
264
265
266
267 if (crossSectionCalculator) {
268
269 grej=xs*2.;
270 do {
272 crossSectionCalculator->
Initialize(x,gam,phi,theBeamPolarization,
273 theTargetPolarization,1);
276 if(xs > grej) {
278 G4cout <<
"G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<
G4endl
279 << "Majorant " << grej << " < "
280 << xs <<
" for phi= " << phi<<
G4endl
281 <<
" e-e- (Moller) scattering"<<
G4endl
283 } else {
284 G4cout <<
"G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<
G4endl
285 << "Majorant " << grej << " < "
286 << xs <<
" for phi= " << phi<<
G4endl
287 <<
" e+e- (Bhabha) scattering"<<
G4endl
289 }
290 }
291
293 }
294
295
296 G4double deltaKinEnergy = x * kineticEnergy;
298 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
300 (deltaMomentum * totalMomentum);
302 if(sint > 0.0) sint = std::sqrt(sint);
303
304
305 G4ThreeVector deltaDirection(-sint*std::cos(phi),-sint*std::sin(phi), cost) ;
306 deltaDirection.rotateUz(direction);
307
308
309 kineticEnergy -= deltaKinEnergy;
311
313 G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
314 direction = dir.
unit();
316 }
317
318
320 vdp->push_back(delta);
321
322
325
326 if (crossSectionCalculator) {
327
328
329 theBeamPolarization.
InvRotateAz(nInteractionFrame,direction);
330 theTargetPolarization.
InvRotateAz(nInteractionFrame,direction);
331 crossSectionCalculator->
Initialize(x,gam,phi,theBeamPolarization,
332 theTargetPolarization,2);
333
334
335 fPositronPolarization=crossSectionCalculator->
GetPol2();
336 fPositronPolarization.
RotateAz(nInteractionFrame,direction);
337
339
340
341 fElectronPolarization=crossSectionCalculator->
GetPol3();
342 fElectronPolarization.
RotateAz(nInteractionFrame,deltaDirection);
344 fElectronPolarization.
y(),
345 fElectronPolarization.
z());
346 }
347 else {
350 }
351}
CLHEP::Hep3Vector G4ThreeVector
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)
const G4Track * GetCurrentTrack() const
void ProposePolarization(const G4ThreeVector &dir)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
bool IsPolarized(G4LogicalVolume *lVol) const
static G4PolarizationManager * GetInstance()
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
G4VPhysicalVolume * GetVolume() const
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4LogicalVolume * GetLogicalVolume() const
virtual G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3)=0
virtual G4StokesVector GetPol3()
virtual G4StokesVector GetPol2()
virtual void Initialize(G4double, G4double, G4double, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)
G4double energy(const ThreeVector &p, const G4double m)