Geant4  10.01.p03
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 82755 2014-07-08 14:07:29Z 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 
60 #include "G4StokesVector.hh"
61 #include "G4PolarizationManager.hh"
62 #include "G4PolarizationHelper.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68  const G4String& nam)
69  : G4KleinNishinaCompton(0,nam),
70  verboseLevel(0)
71 {
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 {
80 }
81 
83  (G4double gammaEnergy, G4double /*Z*/)
84 
85 {
86  G4double asymmetry = 0.0 ;
87 
88  G4double k0 = gammaEnergy / electron_mass_c2 ;
89  G4double k1 = 1 + 2*k0 ;
90 
91  asymmetry = -k0;
92  asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
93  asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);
94 
95  // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
96  if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
97 
98  return asymmetry;
99 }
100 
101 
103  const G4ParticleDefinition* pd,
104  G4double kinEnergy,
105  G4double Z,
106  G4double A,
107  G4double cut,
108  G4double emax)
109 {
110  double xs =
112  Z,A,cut,emax);
114  if (polzz > 0.0) {
115  G4double asym = ComputeAsymmetryPerAtom(kinEnergy, Z);
116  xs*=(1.+polzz*asym);
117  }
118  return xs;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122 
124  std::vector<G4DynamicParticle*>* fvect,
125  const G4MaterialCutsCouple*,
126  const G4DynamicParticle* aDynamicGamma,
128 {
129  // do nothing below the threshold
130  if(aDynamicGamma->GetKineticEnergy() <= LowEnergyLimit()) { return; }
131 
132  const G4Track * aTrack = fParticleChange->GetCurrentTrack();
133  G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
134  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
135 
136  if (verboseLevel >= 1) {
137  G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
138  << aLVolume->GetName() <<G4endl;
139  }
140  G4PolarizationManager * polarizationManager =
142 
143  // obtain polarization of the beam
144  theBeamPolarization = aDynamicGamma->GetPolarization();
146 
147  // obtain polarization of the media
148  G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
150  polarizationManager->GetVolumePolarization(aLVolume);
151 
152  // if beam is linear polarized or target is transversely polarized
153  // determine the angle to x-axis
154  // (assumes same PRF as in the polarization definition)
155 
156  G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
157 
158  // transfere theTargetPolarization
159  // into the gamma frame (problem electron is at rest)
160  if (targetIsPolarized) {
161  theTargetPolarization.rotateUz(gamDirection0);
162  }
163  // The scattered gamma energy is sampled according to
164  // Klein - Nishina formula.
165  // The random number techniques of Butcher & Messel are used
166  // (Nuc Phys 20(1960),15).
167  // Note : Effects due to binding of atomic electrons are negliged.
168 
169  G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
170  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
171 
172  //
173  // sample the energy rate of the scattered gamma
174  //
175 
176  G4double epsilon, epsilonsq, onecost, sint2, greject ;
177 
178  G4double eps0 = 1./(1. + 2.*E0_m);
179  G4double epsilon0sq = eps0*eps0;
180  G4double alpha1 = - std::log(eps0);
181  G4double alpha2 = 0.5*(1.- epsilon0sq);
182 
183  G4double polarization =
185 
186  G4int nloop = 0;
187  do {
188  ++nloop;
189  // false interaction if too many iterations
190  if(nloop > 1000) { return; }
191 
192  if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
193  epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r
194  epsilonsq = epsilon*epsilon;
195 
196  } else {
197  epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
198  epsilon = std::sqrt(epsilonsq);
199  }
200 
201  onecost = (1.- epsilon)/(epsilon*E0_m);
202  sint2 = onecost*(2.-onecost);
203 
204  G4double gdiced = 2.*(1./epsilon+epsilon);
205  G4double gdist = 1./epsilon + epsilon - sint2
206  - polarization*(1./epsilon-epsilon)*(1.-onecost);
207 
208  greject = gdist/gdiced;
209 
210  if (greject>1) {
211  G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
212  <<" costh rejection does not work properly: "<<greject
213  <<G4endl;
214  }
215  } while (greject < G4UniformRand());
216 
217  //
218  // scattered gamma angles. ( Z - axis along the parent gamma)
219  //
220 
221  G4double cosTeta = 1. - onecost;
222  G4double sinTeta = std::sqrt (sint2);
223  G4double Phi;
224  do {
225  ++nloop;
226  // false interaction if too many iterations
227  if(nloop > 1000) { return; }
228 
229  Phi = twopi * G4UniformRand();
230  G4double gdiced = 1./epsilon + epsilon - sint2
231  + std::abs(theBeamPolarization.p3())*
232  ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
233  +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1())
235  +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) +
237 
238  G4double gdist = 1./epsilon + epsilon - sint2
240  ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
241  +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
242  std::sin(Phi)*theTargetPolarization.p2()))
243  -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
244  +std::sin(2.*Phi)*theBeamPolarization.p2());
245  greject = gdist/gdiced;
246 
247  if (greject>1.+1.e-10 || greject<0) {
248  G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
249  <<" phi rejection does not work properly: "<<greject<<G4endl;
250  }
251  if (greject<1.e-3) {
252  G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
253  <<" phi rejection does not work properly: "<<greject<<"\n";
254  G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n";
255  G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n";
256  G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n";
257  }
258 
259  } while (greject < G4UniformRand());
260  G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi),
261  dirz = cosTeta;
262 
263  //
264  // update G4VParticleChange for the scattered gamma
265  //
266 
267  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
268  gamDirection1.rotateUz(gamDirection0);
269  G4double gamEnergy1 = epsilon*gamEnergy0;
270 
271  G4double edep = 0.0;
272  if(gamEnergy1 > lowestSecondaryEnergy) {
275  } else {
278  edep = gamEnergy1;
279  }
280 
281  //
282  // calculate Stokesvector of final state photon and electron
283  //
284  G4ThreeVector nInteractionFrame =
285  G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
286 
287  // transfere theBeamPolarization and theTargetPolarization
288  // into the interaction frame (note electron is in gamma frame)
289  if (verboseLevel>=1) {
290  G4cout << "========================================\n";
291  G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
292  G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
293  G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
294  G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
295  }
296 
297  theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
298  theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
299 
300  if (verboseLevel>=1) {
301  G4cout << "----------------------------------------\n";
302  G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
303  G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
304  G4cout << "----------------------------------------\n";
305  }
306 
307  // initialize the polarization transfer matrix
308  crossSectionCalculator->Initialize(epsilon,E0_m,0.,
311 
312  if(gamEnergy1 > lowestSecondaryEnergy) {
313 
314  // in interaction frame
315  // calculate polarization transfer to the photon (in interaction plane)
317  if (verboseLevel>=1) {
318  G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
319  }
321 
322  // translate polarization into particle reference frame
323  finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
324  //store polarization vector
326  if (finalGammaPolarization.mag() > 1.+1.e-8){
327  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
328  G4cout<<"Polarization of final photon more than 100%"<<G4endl;
329  G4cout<<finalGammaPolarization<<" mag = "
331  }
332  if (verboseLevel>=1) {
333  G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
334  G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
335  }
336  }
337 
338  //
339  // kinematic of the scattered electron
340  //
341  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
342 
343  if (eKinEnergy > lowestSecondaryEnergy) {
344 
345  G4ThreeVector eDirection =
346  gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
347  eDirection = eDirection.unit();
348 
350  if (verboseLevel>=1) {
351  G4cout << " electronPolarization1 = "
353  }
354  // transfer into particle reference frame
355  finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
356  if (verboseLevel>=1) {
357  G4cout << " electronPolarization1 = "
359  G4cout << " ElecDirection = " <<eDirection<<"\n";
360  }
361 
362  // create G4DynamicParticle object for the electron.
363  G4DynamicParticle* aElectron =
364  new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
365  //store polarization vector
366  if (finalElectronPolarization.mag() > 1.+1.e-8){
367  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
368  G4cout<<"Polarization of final electron more than 100%"<<G4endl;
371  }
375  fvect->push_back(aElectron);
376  } else {
377  edep += eKinEnergy;
378  }
379  // energy balance
380  if(edep > 0.0) {
382  }
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386 
387 
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:623
G4ParticleChangeForGamma * fParticleChange
G4String GetName() const
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double p2() const
static G4PolarizationManager * GetInstance()
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double p3() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax)
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4PolarizedComptonCrossSection * crossSectionCalculator
const G4ThreeVector & GetMomentumDirection() const
void ProposePolarization(const G4ThreeVector &dir)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
G4double p1() const
void SetPolarization(G4double polX, G4double polY, G4double polZ)
static const G4double A[nN]
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
virtual void Initialize(G4double eps, G4double X, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)
G4PolarizedComptonModel(const G4ParticleDefinition *p=0, const G4String &nam="Polarized-Compton")
G4LogicalVolume * GetLogicalVolume() const
G4double ComputeAsymmetryPerAtom(G4double gammaEnergy, G4double Z)
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
G4ParticleDefinition * theElectron
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax)
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)