Geant4  9.6.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 69847 2013-05-16 09:36:18Z 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 
51 #include "G4PolarizedComptonModel.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4Electron.hh"
54 #include "G4Gamma.hh"
55 #include "Randomize.hh"
56 #include "G4DataVector.hh"
57 #include "G4ParticleChangeForGamma.hh"
58 
59 
60 #include "G4StokesVector.hh"
61 #include "G4PolarizationManager.hh"
62 #include "G4PolarizationHelper.hh"
63 #include "G4PolarizedComptonCrossSection.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*,
68  const G4String& nam)
69  : G4KleinNishinaCompton(0,nam),
70  verboseLevel(0)
71 {
72  crossSectionCalculator=new G4PolarizedComptonCrossSection();
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
77 G4PolarizedComptonModel::~G4PolarizedComptonModel()
78 {
79  if (crossSectionCalculator) delete crossSectionCalculator;
80 }
81 
82 
83 
84 G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom
85  (G4double gammaEnergy, G4double /*Z*/)
86 
87 {
88  G4double asymmetry = 0.0 ;
89 
90  G4double k0 = gammaEnergy / electron_mass_c2 ;
91  G4double k1 = 1 + 2*k0 ;
92 
93  asymmetry = -k0;
94  asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
95  asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);
96 
97  // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
98  if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
99 
100  return asymmetry;
101 }
102 
103 
104 G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom(
105  const G4ParticleDefinition* pd,
106  G4double kinEnergy,
107  G4double Z,
108  G4double A,
109  G4double cut,
110  G4double emax)
111 {
112  double xs =
113  G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy,
114  Z,A,cut,emax);
115  G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
116  if (polzz!=0) {
117  G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z);
118  xs*=(1.+polzz*asym);
119  }
120  return xs;
121 }
122 
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 
126 void G4PolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
127  const G4MaterialCutsCouple*,
128  const G4DynamicParticle* aDynamicGamma,
129  G4double,
130  G4double)
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 = G4PolarizationManager::GetInstance();
141 
142  // obtain polarization of the beam
143  theBeamPolarization = aDynamicGamma->GetPolarization();
144  theBeamPolarization.SetPhoton();
145 
146  // obtain polarization of the media
147  const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
148  theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
149 
150  // if beam is linear polarized or target is transversely polarized
151  // determine the angle to x-axis
152  // (assumes same PRF as in the polarization definition)
153 
154  G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
155 
156  // transfere theTargetPolarization
157  // into the gamma frame (problem electron is at rest)
158  if (targetIsPolarized)
159  theTargetPolarization.rotateUz(gamDirection0);
160 
161  // The scattered gamma energy is sampled according to Klein - Nishina formula.
162  // The random number techniques of Butcher & Messel are used
163  // (Nuc Phys 20(1960),15).
164  // Note : Effects due to binding of atomic electrons are negliged.
165 
166  G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
167  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
168 
169 
170  //
171  // sample the energy rate of the scattered gamma
172  //
173 
174  G4double epsilon, epsilonsq, onecost, sint2, greject ;
175 
176  G4double eps0 = 1./(1. + 2.*E0_m);
177  G4double epsilon0sq = eps0*eps0;
178  G4double alpha1 = - std::log(eps0);
179  G4double alpha2 = 0.5*(1.- epsilon0sq);
180 
181  G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
182  do {
183  if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
184  epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r
185  epsilonsq = epsilon*epsilon;
186 
187  } else {
188  epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
189  epsilon = std::sqrt(epsilonsq);
190  };
191 
192  onecost = (1.- epsilon)/(epsilon*E0_m);
193  sint2 = onecost*(2.-onecost);
194 
195 
196  G4double gdiced = 2.*(1./epsilon+epsilon);
197  G4double gdist = 1./epsilon + epsilon - sint2
198  - polarization*(1./epsilon-epsilon)*(1.-onecost);
199 
200  greject = gdist/gdiced;
201 
202  if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
203  <<" costh rejection does not work properly: "<<greject<<G4endl;
204 
205  } while (greject < G4UniformRand());
206 
207  //
208  // scattered gamma angles. ( Z - axis along the parent gamma)
209  //
210 
211  G4double cosTeta = 1. - onecost;
212  G4double sinTeta = std::sqrt (sint2);
213  G4double Phi;
214  do {
215  Phi = twopi * G4UniformRand();
216  G4double gdiced = 1./epsilon + epsilon - sint2
217  + std::abs(theBeamPolarization.p3())*
218  ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
219  +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1())
220  + sqr(theTargetPolarization.p2()))))
221  +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2())));
222 
223  G4double gdist = 1./epsilon + epsilon - sint2
224  + theBeamPolarization.p3()*
225  ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
226  +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
227  std::sin(Phi)*theTargetPolarization.p2()))
228  -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
229  +std::sin(2.*Phi)*theBeamPolarization.p2());
230  greject = gdist/gdiced;
231 
232  if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
233  <<" phi rejection does not work properly: "<<greject<<G4endl;
234 
235  if (greject<1.e-3) {
236  G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
237  <<" phi rejection does not work properly: "<<greject<<"\n";
238  G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n";
239  G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n";
240  G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n";
241  }
242 
243  } while (greject < G4UniformRand());
244  G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta;
245 
246  //
247  // update G4VParticleChange for the scattered gamma
248  //
249 
250  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
251  gamDirection1.rotateUz(gamDirection0);
252  G4double gamEnergy1 = epsilon*gamEnergy0;
253  fParticleChange->SetProposedKineticEnergy(gamEnergy1);
254 
255 
256  if(gamEnergy1 > lowestGammaEnergy) {
257  fParticleChange->ProposeMomentumDirection(gamDirection1);
258  } else {
259  fParticleChange->ProposeTrackStatus(fStopAndKill);
260  gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
261  fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
262  }
263 
264  //
265  // kinematic of the scattered electron
266  //
267 
268  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
269  G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
270  eDirection = eDirection.unit();
271 
272  //
273  // calculate Stokesvector of final state photon and electron
274  //
275  G4ThreeVector nInteractionFrame;
276  if((gamEnergy1 > lowestGammaEnergy) ||
277  (eKinEnergy > DBL_MIN)) {
278 
279  // determine interaction plane
280 // nInteractionFrame =
281 // G4PolarizationHelper::GetFrame(gamDirection1,eDirection);
282  if (gamEnergy1 > lowestGammaEnergy)
283  nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
284  else
285  nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection);
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.,
309  theBeamPolarization,
310  theTargetPolarization,2);
311  }
312 
313  // if(eKinEnergy > DBL_MIN)
314  {
315  // in interaction frame
316  // calculate polarization transfer to the photon (in interaction plane)
317  finalGammaPolarization = crossSectionCalculator->GetPol2();
318  if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
319  finalGammaPolarization.SetPhoton();
320 
321  // translate polarization into particle reference frame
322  finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
323  //store polarization vector
324  fParticleChange->ProposePolarization(finalGammaPolarization);
325  if (finalGammaPolarization.mag() > 1.+1.e-8){
326  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
327  G4cout<<"Polarization of final photon more than 100%"<<G4endl;
328  G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl;
329  }
330  if (verboseLevel>=1) {
331  G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
332  G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
333  }
334  }
335 
336  // if (ElecKineEnergy > fminimalEnergy) {
337  {
338  finalElectronPolarization = crossSectionCalculator->GetPol3();
339  if (verboseLevel>=1)
340  G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
341 
342  // transfer into particle reference frame
343  finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
344  if (verboseLevel>=1) {
345  G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n";
346  G4cout << " ElecDirection = " <<eDirection<<"\n";
347  }
348  }
349  if (verboseLevel>=1)
350  G4cout << "========================================\n";
351 
352 
353  if(eKinEnergy > DBL_MIN) {
354 
355  // create G4DynamicParticle object for the electron.
356  G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
357  //store polarization vector
358  if (finalElectronPolarization.mag() > 1.+1.e-8){
359  G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
360  G4cout<<"Polarization of final electron more than 100%"<<G4endl;
361  G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl;
362  }
363  aElectron->SetPolarization(finalElectronPolarization.p1(),
364  finalElectronPolarization.p2(),
365  finalElectronPolarization.p3());
366  fvect->push_back(aElectron);
367  }
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
371 
372 
Generated on Sat May 25 2013 14:33:36 for Geant4 by   doxygen 1.8.4