Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedComptonModel.cc
Go to the documentation of this file.
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 // $Id: G4PolarizedComptonModel.cc 96114 2016-03-16 18:51:33Z gcosmo $
28 //
29 // -------------------------------------------------------------------
30 //
31 // GEANT4 Class file
32 //
33 //
34 // File name: G4PolarizedComptonModel
35 //
36 // Author: Andreas Schaelicke
37 //
38 // Creation date: 01.05.2005
39 //
40 // Modifications:
41 // 18-07-06 use newly calculated cross sections (P. Starovoitov)
42 // 21-08-05 update interface (A. Schaelicke)
43 //
44 // Class Description:
45 //
46 // -------------------------------------------------------------------
47 //
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
52 #include "G4PhysicalConstants.hh"
53 #include "G4Electron.hh"
54 #include "G4Gamma.hh"
55 #include "Randomize.hh"
56 #include "G4DataVector.hh"
58 
59 #include "G4StokesVector.hh"
60 #include "G4PolarizationManager.hh"
61 #include "G4PolarizationHelper.hh"
63 
64 #include "G4SystemOfUnits.hh"
65 #include "G4Log.hh"
66 #include "G4Exp.hh"
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
70 static const G4int nlooplim = 10000;
71 
73  const G4String& nam)
74  : G4KleinNishinaCompton(nullptr,nam),
75  verboseLevel(0)
76 {
77  crossSectionCalculator = new G4PolarizedComptonCrossSection();
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
83 {
84  delete crossSectionCalculator;
85 }
86 
88  (G4double gammaEnergy, G4double /*Z*/)
89 
90 {
91  G4double asymmetry = 0.0 ;
92 
93  G4double k0 = gammaEnergy / electron_mass_c2 ;
94  G4double k1 = 1. + 2.*k0 ;
95 
96  asymmetry = -k0;
97  asymmetry *= (k0 + 1.)*sqr(k1)*G4Log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
98  asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*G4Log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);
99 
100  // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
101  if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
102 
103  return asymmetry;
104 }
105 
106 
108  const G4ParticleDefinition* pd,
109  G4double kinEnergy,
110  G4double Z,
111  G4double A,
112  G4double cut,
113  G4double emax)
114 {
115  double xs =
117  Z,A,cut,emax);
118  G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
119  if (polzz > 0.0) {
120  G4double asym = ComputeAsymmetryPerAtom(kinEnergy, Z);
121  xs *= (1.+polzz*asym);
122  }
123  return xs;
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127 
129  std::vector<G4DynamicParticle*>* fvect,
130  const G4MaterialCutsCouple*,
131  const G4DynamicParticle* aDynamicGamma,
133 {
134  // do nothing below the threshold
135  if(aDynamicGamma->GetKineticEnergy() <= LowEnergyLimit()) { return; }
136 
137  const G4Track * aTrack = fParticleChange->GetCurrentTrack();
138  G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
139  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
140 
141  if (verboseLevel >= 1) {
142  G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
143  << aLVolume->GetName() <<G4endl;
144  }
145  G4PolarizationManager * polarizationManager =
147 
148  // obtain polarization of the beam
149  theBeamPolarization = aDynamicGamma->GetPolarization();
150  theBeamPolarization.SetPhoton();
151 
152  // obtain polarization of the media
153  G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
154  theTargetPolarization =
155  polarizationManager->GetVolumePolarization(aLVolume);
156 
157  // if beam is linear polarized or target is transversely polarized
158  // determine the angle to x-axis
159  // (assumes same PRF as in the polarization definition)
160 
161  G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
162 
163  // transfere theTargetPolarization
164  // into the gamma frame (problem electron is at rest)
165  if (targetIsPolarized) {
166  theTargetPolarization.rotateUz(gamDirection0);
167  }
168  // The scattered gamma energy is sampled according to
169  // Klein - Nishina formula.
170  // The random number techniques of Butcher & Messel are used
171  // (Nuc Phys 20(1960),15).
172  // Note : Effects due to binding of atomic electrons are negliged.
173 
174  G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
175  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
176 
177  //
178  // sample the energy rate of the scattered gamma
179  //
180 
181  G4double epsilon, sint2;
182  G4double onecost = 0.0;
183  G4double Phi = 0.0;
184  G4double greject = 1.0;
185  G4double cosTeta = 1.0;
186  G4double sinTeta = 0.0;
187 
188  G4double eps0 = 1./(1. + 2.*E0_m);
189  G4double epsilon0sq = eps0*eps0;
190  G4double alpha1 = - G4Log(eps0);
191  G4double alpha2 = alpha1 + 0.5*(1.- epsilon0sq);
192 
193  G4double polarization =
194  theBeamPolarization.p3()*theTargetPolarization.p3();
195 
196  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
197  G4int nloop = 0;
198  G4bool end = false;
199 
200  G4double rndm[3];
201 
202  do {
203  do {
204  ++nloop;
205  // false interaction if too many iterations
206  if(nloop > nlooplim) {
207  PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
208  "too many iterations");
209  return;
210  }
211 
212  // 3 random numbers to sample scattering
213  rndmEngineMod->flatArray(3, rndm);
214 
215  if ( alpha1 > alpha2*rndm[0]) {
216  epsilon = G4Exp(-alpha1*rndm[1]); // epsilon0**r
217  } else {
218  epsilon = std::sqrt(epsilon0sq + (1.- epsilon0sq)*rndm[1]);
219  }
220 
221  onecost = (1.- epsilon)/(epsilon*E0_m);
222  sint2 = onecost*(2.-onecost);
223 
224  G4double gdiced = 2.*(1./epsilon+epsilon);
225  G4double gdist = 1./epsilon + epsilon - sint2
226  - polarization*(1./epsilon-epsilon)*(1.-onecost);
227 
228  greject = gdist/gdiced;
229 
230  if (greject > 1.0) {
231  PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
232  "theta majoranta wrong");
233  }
234  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
235  } while (greject < rndm[2]);
236 
237  // assuming phi loop sucessful
238  end = true;
239 
240  //
241  // scattered gamma angles. ( Z - axis along the parent gamma)
242  //
243  cosTeta = 1. - onecost;
244  sinTeta = std::sqrt(sint2);
245  do {
246  ++nloop;
247 
248  // 2 random numbers to sample scattering
249  rndmEngineMod->flatArray(2, rndm);
250 
251  // false interaction if too many iterations
252  Phi = twopi * rndm[0];
253  if(nloop > nlooplim) {
254  PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
255  "too many iterations");
256  return;
257  }
258 
259  G4double gdiced = 1./epsilon + epsilon - sint2
260  + std::abs(theBeamPolarization.p3())*
261  ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
262  +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1())
263  + sqr(theTargetPolarization.p2()))))
264  +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) +
265  sqr(theBeamPolarization.p2())));
266 
267  G4double gdist = 1./epsilon + epsilon - sint2
268  + theBeamPolarization.p3()*
269  ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
270  +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
271  std::sin(Phi)*theTargetPolarization.p2()))
272  -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
273  +std::sin(2.*Phi)*theBeamPolarization.p2());
274  greject = gdist/gdiced;
275 
276  if (greject > 1.0) {
277  PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
278  "phi majoranta wrong");
279  }
280 
281  if(greject < 1.e-3) {
282  PrintWarning(aDynamicGamma, nloop, greject, onecost, Phi,
283  "phi loop ineffective");
284  // restart theta loop
285  end = false;
286  break;
287  }
288 
289  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
290  } while (greject < rndm[1]);
291  } while(!end);
292  G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi),
293  dirz = cosTeta;
294 
295  //
296  // update G4VParticleChange for the scattered gamma
297  //
298 
299  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
300  gamDirection1.rotateUz(gamDirection0);
301  G4double gamEnergy1 = epsilon*gamEnergy0;
302 
303  G4double edep = 0.0;
304  if(gamEnergy1 > lowestSecondaryEnergy) {
307  } else {
310  edep = gamEnergy1;
311  }
312 
313  //
314  // calculate Stokesvector of final state photon and electron
315  //
316  G4ThreeVector nInteractionFrame =
317  G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
318 
319  // transfere theBeamPolarization and theTargetPolarization
320  // into the interaction frame (note electron is in gamma frame)
321  if (verboseLevel>=1) {
322  G4cout << "========================================\n";
323  G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
324  G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
325  G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
326  G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
327  }
328 
329  theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
330  theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
331 
332  if (verboseLevel>=1) {
333  G4cout << "----------------------------------------\n";
334  G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
335  G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
336  G4cout << "----------------------------------------\n";
337  }
338 
339  // initialize the polarization transfer matrix
340  crossSectionCalculator->Initialize(epsilon,E0_m,0.,
341  theBeamPolarization,
342  theTargetPolarization,2);
343 
344  if(gamEnergy1 > lowestSecondaryEnergy) {
345 
346  // in interaction frame
347  // calculate polarization transfer to the photon (in interaction plane)
348  finalGammaPolarization = crossSectionCalculator->GetPol2();
349  if (verboseLevel>=1) {
350  G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
351  }
352  finalGammaPolarization.SetPhoton();
353 
354  // translate polarization into particle reference frame
355  finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
356  if (finalGammaPolarization.mag() > 1.+1.e-8){
357  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
358  G4cout<<"Polarization of final photon more than 100%"<<G4endl;
359  G4cout<<finalGammaPolarization<<" mag = "
360  <<finalGammaPolarization.mag()<<G4endl;
361  }
362  //store polarization vector
363  fParticleChange->ProposePolarization(finalGammaPolarization);
364  if (verboseLevel>=1) {
365  G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
366  G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
367  }
368  }
369 
370  //
371  // kinematic of the scattered electron
372  //
373  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
374 
375  if (eKinEnergy > lowestSecondaryEnergy) {
376 
377  G4ThreeVector eDirection =
378  gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
379  eDirection = eDirection.unit();
380 
381  finalElectronPolarization = crossSectionCalculator->GetPol3();
382  if (verboseLevel>=1) {
383  G4cout << " electronPolarization1 = "
384  <<finalElectronPolarization<<"\n";
385  }
386  // transfer into particle reference frame
387  finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
388  if (verboseLevel>=1) {
389  G4cout << " electronPolarization1 = "
390  <<finalElectronPolarization<<"\n";
391  G4cout << " ElecDirection = " <<eDirection<<"\n";
392  }
393 
394  // create G4DynamicParticle object for the electron.
395  G4DynamicParticle* aElectron =
396  new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
397  //store polarization vector
398  if (finalElectronPolarization.mag() > 1.+1.e-8){
399  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
400  G4cout<<"Polarization of final electron more than 100%"<<G4endl;
401  G4cout<<finalElectronPolarization<<" mag = "
402  <<finalElectronPolarization.mag()<<G4endl;
403  }
404  aElectron->SetPolarization(finalElectronPolarization.p1(),
405  finalElectronPolarization.p2(),
406  finalElectronPolarization.p3());
407  fvect->push_back(aElectron);
408  } else {
409  edep += eKinEnergy;
410  }
411  // energy balance
412  if(edep > 0.0) {
414  }
415 }
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418 
419 void
420 G4PolarizedComptonModel::PrintWarning(const G4DynamicParticle* dp, G4int nloop,
421  G4double grej, G4double onecos,
422  G4double phi, const G4String sss) const
423 {
425  ed << "Problem of scattering sampling: " << sss << "\n"
426  << "Niter= " << nloop << " grej= " << grej << " cos(theta)= "
427  << 1.0-onecos << " phi= " << phi << "\n"
428  << "Gamma E(MeV)= " << dp->GetKineticEnergy()/MeV
429  << " dir= " << dp->GetMomentumDirection()
430  << " pol= " << dp->GetPolarization();
431  G4Exception("G4PolarizedComptonModel::SampleSecondaries","em0044",
432  JustWarning, ed, "");
433 
434 }
435 
436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
437 
438 
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4ParticleChangeForGamma * fParticleChange
G4PolarizedComptonModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="Polarized-Compton")
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
G4double p2() const
static const G4int nlooplim
static G4PolarizationManager * GetInstance()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax) override
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double z() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double p3() const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void ProposePolarization(const G4ThreeVector &dir)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
G4double p1() const
float electron_mass_c2
Definition: hepunit.py:274
void SetPolarization(G4double polX, G4double polY, G4double polZ)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
G4LogicalVolume * GetLogicalVolume() const
static const char sss[MAX_N_PAR+2]
Definition: Evaluator.cc:64
G4double ComputeAsymmetryPerAtom(G4double gammaEnergy, G4double Z)
Hep3Vector unit() const
bool IsPolarized(G4LogicalVolume *lVol) const
const G4ThreeVector & GetPolarization() const
G4VPhysicalVolume * GetVolume() const
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
virtual void Initialize(G4double eps, G4double X, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0) override
G4ParticleDefinition * theElectron
T sqr(const T &x)
Definition: templates.hh:145
const G4String & GetName() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual void flatArray(const int size, double *vect)=0
double mag() const
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
double epsilon(double density, double temperature)
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax) override