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