110{
111
113 {
114 return;
115 }
116
120
121 if(fVerboseLevel >= 1)
122 {
123 G4cout <<
"G4PolarizedComptonModel::SampleSecondaries in "
125 }
128
129
132
133
136
137
138
139
141
142
143
144 if(targetIsPolarized)
145 {
146 fTargetPolarization.
rotateUz(gamDirection0);
147 }
148
149
150
151
152
153
155 G4double E0_m = gamEnergy0 / electron_mass_c2;
156
157
164
165 G4double eps0 = 1. / (1. + 2. * E0_m);
169
170 G4double polarization = fBeamPolarization.
p3() * fTargetPolarization.
p3();
171
175
177
178 do
179 {
180 do
181 {
182 ++nloop;
183
184 if(nloop > fLoopLim)
185 {
186 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
187 "too many iterations");
188 return;
189 }
190
191
193
194 if(alpha1 >
alpha2 * rndm[0])
195 {
197 }
198 else
199 {
200 epsilon = std::sqrt(epsilon0sq + (1. - epsilon0sq) * rndm[1]);
201 }
202
204 sint2 = onecost * (2. - onecost);
205
209
210 greject = gdist / gdiced;
211
212 if(greject > 1.0)
213 {
214 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
215 "theta majoranta wrong");
216 }
217
218 } while(greject < rndm[2]);
219
220
221 end = true;
222
223
224 cosTeta = 1. - onecost;
225 sinTeta = std::sqrt(sint2);
226 do
227 {
228 ++nloop;
229
230
232
233
234 Phi = twopi * rndm[0];
235 if(nloop > fLoopLim)
236 {
237 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
238 "too many iterations");
239 return;
240 }
241
243 std::abs(fBeamPolarization.
p3()) *
245 fTargetPolarization.
p3()) +
247 (std::sqrt(
sqr(fTargetPolarization.
p1()) +
248 sqr(fTargetPolarization.
p2())))) +
249 sint2 * (std::sqrt(
sqr(fBeamPolarization.
p1()) +
250 sqr(fBeamPolarization.
p2())));
251
254 fBeamPolarization.
p3() *
257 (std::cos(Phi) * fTargetPolarization.
p1() +
258 std::sin(Phi) * fTargetPolarization.
p2())) -
259 sint2 * (std::cos(2. * Phi) * fBeamPolarization.
p1() +
260 std::sin(2. * Phi) * fBeamPolarization.
p2());
261 greject = gdist / gdiced;
262
263 if(greject > 1.0)
264 {
265 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
266 "phi majoranta wrong");
267 }
268
269 if(greject < 1.e-3)
270 {
271 PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
272 "phi loop ineffective");
273
274 end = false;
275 break;
276 }
277
278
279 } while(greject < rndm[1]);
280 } while(!end);
281 G4double dirx = sinTeta * std::cos(Phi);
282 G4double diry = sinTeta * std::sin(Phi);
284
285
287 gamDirection1.rotateUz(gamDirection0);
289
292 {
295 }
296 else
297 {
300 edep = gamEnergy1;
301 }
302
303
306
307
308
309 if(fVerboseLevel >= 1)
310 {
311 G4cout <<
"========================================" <<
G4endl;
312 G4cout <<
" nInteractionFrame = " << nInteractionFrame <<
G4endl;
313 G4cout <<
" GammaDirection0 = " << gamDirection0 <<
G4endl;
314 G4cout <<
" gammaPolarization = " << fBeamPolarization <<
G4endl;
315 G4cout <<
" electronPolarization = " << fTargetPolarization <<
G4endl;
316 }
317
318 fBeamPolarization.
InvRotateAz(nInteractionFrame, gamDirection0);
319 fTargetPolarization.
InvRotateAz(nInteractionFrame, gamDirection0);
320
321 if(fVerboseLevel >= 1)
322 {
323 G4cout <<
"----------------------------------------" <<
G4endl;
324 G4cout <<
" gammaPolarization = " << fBeamPolarization <<
G4endl;
325 G4cout <<
" electronPolarization = " << fTargetPolarization <<
G4endl;
326 G4cout <<
"----------------------------------------" <<
G4endl;
327 }
328
329
331 fTargetPolarization, 2);
332
334 {
335
336
337 fFinalGammaPolarization = fCrossSectionCalculator->
GetPol2();
338 if(fVerboseLevel >= 1)
339 {
340 G4cout <<
" gammaPolarization1 = " << fFinalGammaPolarization <<
G4endl;
341 }
343
344
345 fFinalGammaPolarization.
RotateAz(nInteractionFrame, gamDirection1);
346 if(fFinalGammaPolarization.
mag() > 1. + 1.e-8)
347 {
349 ed << "ERROR in Polarizaed Compton Scattering !\n";
350 ed << "Polarization of final photon more than 100%.\n";
351 ed << fFinalGammaPolarization
352 <<
" mag = " << fFinalGammaPolarization.
mag() <<
"\n";
353 G4Exception(
"G4PolarizedComptonModel::SampleSecondaries",
"pol033",
355 }
356
358 if(fVerboseLevel >= 1)
359 {
360 G4cout <<
" gammaPolarization1 = " << fFinalGammaPolarization <<
G4endl;
361 G4cout <<
" GammaDirection1 = " << gamDirection1 <<
G4endl;
362 }
363 }
364
365
366 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
367
369 {
371 gamEnergy0 * gamDirection0 - gamEnergy1 * gamDirection1;
372 eDirection = eDirection.
unit();
373
374 finalElectronPolarization = fCrossSectionCalculator->
GetPol3();
375 if(fVerboseLevel >= 1)
376 {
377 G4cout <<
" electronPolarization1 = " << finalElectronPolarization
379 }
380
381 finalElectronPolarization.
RotateAz(nInteractionFrame, eDirection);
382 if(fVerboseLevel >= 1)
383 {
384 G4cout <<
" electronPolarization1 = " << finalElectronPolarization
385 <<
G4endl <<
" ElecDirection = " << eDirection <<
G4endl;
386 }
387
388
391
392 if(finalElectronPolarization.
mag() > 1. + 1.e-8)
393 {
395 ed << "ERROR in Polarized Compton Scattering !\n";
396 ed << "Polarization of final electron more than 100%.\n";
397 ed << finalElectronPolarization
398 <<
" mag = " << finalElectronPolarization.
mag() <<
G4endl;
399 G4Exception(
"G4PolarizedComptonModel::SampleSecondaries",
"pol034",
401 }
403 finalElectronPolarization.
p2(),
404 finalElectronPolarization.
p3());
405 fvect->push_back(aElectron);
406 }
407 else
408 {
409 edep += eKinEnergy;
410 }
411
412 if(edep > 0.0)
413 {
415 }
416}
G4double epsilon(G4double density, G4double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
virtual void flatArray(const int size, double *vect)=0
void SetPolarization(const G4ThreeVector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetPolarization() const
G4double lowestSecondaryEnergy
G4ParticleChangeForGamma * fParticleChange
G4ParticleDefinition * theElectron
const G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposePolarization(const G4ThreeVector &dir)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
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 LowEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
const G4Track * GetCurrentTrack() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4LogicalVolume * GetLogicalVolume() const
virtual void Initialize(G4double, G4double, G4double, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)=0
virtual G4StokesVector GetPol2()
virtual G4StokesVector GetPol3()