ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyPolarizedCompton.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LowEnergyPolarizedCompton.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 // GEANT 4 class implementation file
30 // CERN Geneva Switzerland
31 //
32 
33 // --------- G4LowEnergyPolarizedCompton class -----
34 //
35 // by G.Depaola & F.Longo (21 may 2001)
36 //
37 // 21 May 2001 - MGP Modified to inherit from G4VDiscreteProcess
38 // Applies same algorithm as LowEnergyCompton
39 // if the incoming photon is not polarised
40 // Temporary protection to avoid crash in the case
41 // of polarisation || incident photon direction
42 //
43 // 17 October 2001 - F.Longo - Revised according to a design iteration
44 //
45 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola
46 // - better description of parallelism
47 // - system of ref change method improved
48 // 22 January 2003 - V.Ivanchenko Cut per region
49 // 24 April 2003 - V.Ivanchenko Cut per region mfpt
50 //
51 //
52 // ************************************************************
53 //
54 // Corrections by Rui Curado da Silva (2000)
55 // New Implementation by G.Depaola & F.Longo
56 //
57 // - sampling of phi
58 // - polarization of scattered photon
59 //
60 // --------------------------------------------------------------
61 
63 #include "Randomize.hh"
64 #include "G4PhysicalConstants.hh"
65 #include "G4SystemOfUnits.hh"
66 #include "G4ParticleDefinition.hh"
67 #include "G4Track.hh"
68 #include "G4Step.hh"
69 #include "G4ForceCondition.hh"
70 #include "G4Gamma.hh"
71 #include "G4Electron.hh"
72 #include "G4DynamicParticle.hh"
73 #include "G4VParticleChange.hh"
74 #include "G4ThreeVector.hh"
77 #include "G4RDVEMDataSet.hh"
79 #include "G4RDVDataSetAlgorithm.hh"
81 #include "G4RDVRangeTest.hh"
82 #include "G4RDRangeTest.hh"
83 #include "G4MaterialCutsCouple.hh"
84 
85 // constructor
86 
88  : G4VDiscreteProcess(processName),
89  lowEnergyLimit (250*eV), // initialization
90  highEnergyLimit(100*GeV),
91  intrinsicLowEnergyLimit(10*eV),
92  intrinsicHighEnergyLimit(100*GeV)
93 {
96  {
97  G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
98  "OutOfRange", FatalException,
99  "Energy outside intrinsic process validity range!");
100  }
101 
103 
104 
105  G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
106  G4String scatterFile = "comp/ce-sf-";
107  scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.);
108  scatterFunctionData->LoadData(scatterFile);
109 
110  meanFreePathTable = 0;
111 
112  rangeTest = new G4RDRangeTest;
113 
114  // For Doppler broadening
116 
117 
118  if (verboseLevel > 0)
119  {
120  G4cout << GetProcessName() << " is created " << G4endl
121  << "Energy range: "
122  << lowEnergyLimit / keV << " keV - "
123  << highEnergyLimit / GeV << " GeV"
124  << G4endl;
125  }
126 }
127 
128 // destructor
129 
131 {
132  delete meanFreePathTable;
133  delete crossSectionHandler;
134  delete scatterFunctionData;
135  delete rangeTest;
136 }
137 
138 
140 {
141 
143  G4String crossSectionFile = "comp/ce-cs-";
144  crossSectionHandler->LoadData(crossSectionFile);
145  delete meanFreePathTable;
147 
148  // For Doppler broadening
149  G4String file = "/doppler/shell-doppler";
150  shellData.LoadData(file);
151 
152 }
153 
155  const G4Step& aStep)
156 {
157  // The scattered gamma energy is sampled according to Klein - Nishina formula.
158  // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
159  // GEANT4 internal units
160  //
161  // Note : Effects due to binding of atomic electrons are negliged.
162 
163  aParticleChange.Initialize(aTrack);
164 
165  // Dynamic particle quantities
166  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
167  G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
168  G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
169 
170  // gammaPolarization0 = gammaPolarization0.unit(); //
171 
172  // Protection: a polarisation parallel to the
173  // direction causes problems;
174  // in that case find a random polarization
175 
176  G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
177  // ---- MGP ---- Next two lines commented out to remove compilation warnings
178  // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
179  // G4double angle = gammaPolarization0.angle(gammaDirection0);
180 
181  // Make sure that the polarization vector is perpendicular to the
182  // gamma direction. If not
183 
184  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185  { // only for testing now
186  gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187  }
188  else
189  {
190  if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191  {
192  gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193  }
194  }
195 
196  // End of Protection
197 
198  // Within energy limit?
199 
200  if(gammaEnergy0 <= lowEnergyLimit)
201  {
205  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
206  }
207 
208  G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
209 
210  // Select randomly one element in the current material
211 
212  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
213  G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
214 
215  // Sample the energy and the polarization of the scattered photon
216 
217  G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
218 
219  G4double epsilon0 = 1./(1. + 2*E0_m);
220  G4double epsilon0Sq = epsilon0*epsilon0;
221  G4double alpha1 = - std::log(epsilon0);
222  G4double alpha2 = 0.5*(1.- epsilon0Sq);
223 
224  G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225  G4double gammaEnergy1;
226  G4ThreeVector gammaDirection1;
227 
228  do {
229  if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230  {
231  epsilon = std::exp(-alpha1*G4UniformRand());
232  epsilonSq = epsilon*epsilon;
233  }
234  else
235  {
236  epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
237  epsilon = std::sqrt(epsilonSq);
238  }
239 
240  onecost = (1.- epsilon)/(epsilon*E0_m);
241  sinThetaSqr = onecost*(2.-onecost);
242 
243  // Protection
244  if (sinThetaSqr > 1.)
245  {
246  if (verboseLevel>0) G4cout
247  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
248  << "sin(theta)**2 = "
249  << sinThetaSqr
250  << "; set to 1"
251  << G4endl;
252  sinThetaSqr = 1.;
253  }
254  if (sinThetaSqr < 0.)
255  {
256  if (verboseLevel>0) G4cout
257  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
258  << "sin(theta)**2 = "
259  << sinThetaSqr
260  << "; set to 0"
261  << G4endl;
262  sinThetaSqr = 0.;
263  }
264  // End protection
265 
266  G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
267  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268  greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
269 
270  } while(greject < G4UniformRand()*Z);
271 
272 
273  // ****************************************************
274  // Phi determination
275  // ****************************************************
276 
277  G4double phi = SetPhi(epsilon,sinThetaSqr);
278 
279  //
280  // scattered gamma angles. ( Z - axis along the parent gamma)
281  //
282 
283  G4double cosTheta = 1. - onecost;
284 
285  // Protection
286 
287  if (cosTheta > 1.)
288  {
289  if (verboseLevel>0) G4cout
290  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
291  << "cosTheta = "
292  << cosTheta
293  << "; set to 1"
294  << G4endl;
295  cosTheta = 1.;
296  }
297  if (cosTheta < -1.)
298  {
299  if (verboseLevel>0) G4cout
300  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
301  << "cosTheta = "
302  << cosTheta
303  << "; set to -1"
304  << G4endl;
305  cosTheta = -1.;
306  }
307  // End protection
308 
309 
310  G4double sinTheta = std::sqrt (sinThetaSqr);
311 
312  // Protection
313  if (sinTheta > 1.)
314  {
315  if (verboseLevel>0) G4cout
316  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
317  << "sinTheta = "
318  << sinTheta
319  << "; set to 1"
320  << G4endl;
321  sinTheta = 1.;
322  }
323  if (sinTheta < -1.)
324  {
325  if (verboseLevel>0) G4cout
326  << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
327  << "sinTheta = "
328  << sinTheta
329  << "; set to -1"
330  << G4endl;
331  sinTheta = -1.;
332  }
333  // End protection
334 
335 
336  G4double dirx = sinTheta*std::cos(phi);
337  G4double diry = sinTheta*std::sin(phi);
338  G4double dirz = cosTheta ;
339 
340 
341  // oneCosT , eom
342 
343 
344 
345  // Doppler broadening - Method based on:
346  // Y. Namito, S. Ban and H. Hirayama,
347  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
348  // NIM A 349, pp. 489-494, 1994
349 
350  // Maximum number of sampling iterations
351 
352  G4int maxDopplerIterations = 1000;
353  G4double bindingE = 0.;
354  G4double photonEoriginal = epsilon * gammaEnergy0;
355  G4double photonE = -1.;
356  G4int iteration = 0;
357  G4double eMax = gammaEnergy0;
358 
359  do
360  {
361  iteration++;
362  // Select shell based on shell occupancy
363  G4int shell = shellData.SelectRandomShell(Z);
364  bindingE = shellData.BindingEnergy(Z,shell);
365 
366  eMax = gammaEnergy0 - bindingE;
367 
368  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
369  G4double pSample = profileData.RandomSelectMomentum(Z,shell);
370  // Rescale from atomic units
371  G4double pDoppler = pSample * fine_structure_const;
372  G4double pDoppler2 = pDoppler * pDoppler;
373  G4double var2 = 1. + onecost * E0_m;
374  G4double var3 = var2*var2 - pDoppler2;
375  G4double var4 = var2 - pDoppler2 * cosTheta;
376  G4double var = var4*var4 - var3 + pDoppler2 * var3;
377  if (var > 0.)
378  {
379  G4double varSqrt = std::sqrt(var);
380  G4double scale = gammaEnergy0 / var3;
381  // Random select either root
382  if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
383  else photonE = (var4 + varSqrt) * scale;
384  }
385  else
386  {
387  photonE = -1.;
388  }
389  } while ( iteration <= maxDopplerIterations &&
390  (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
391 
392  // End of recalculation of photon energy with Doppler broadening
393  // Revert to original if maximum number of iterations threshold has been reached
394  if (iteration >= maxDopplerIterations)
395  {
396  photonE = photonEoriginal;
397  bindingE = 0.;
398  }
399 
400  gammaEnergy1 = photonE;
401 
402  // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
403 
404 
406 
407 
408 
409 
410  //
411  // update G4VParticleChange for the scattered photon
412  //
413 
414  // gammaEnergy1 = epsilon*gammaEnergy0;
415 
416 
417  // New polarization
418 
419  G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
420  sinThetaSqr,
421  phi,
422  cosTheta);
423 
424  // Set new direction
425  G4ThreeVector tmpDirection1( dirx,diry,dirz );
426  gammaDirection1 = tmpDirection1;
427 
428  // Change reference frame.
429 
430  SystemOfRefChange(gammaDirection0,gammaDirection1,
431  gammaPolarization0,gammaPolarization1);
432 
433  if (gammaEnergy1 > 0.)
434  {
435  aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
436  aParticleChange.ProposeMomentumDirection( gammaDirection1 );
437  aParticleChange.ProposePolarization( gammaPolarization1 );
438  }
439  else
440  {
443  }
444 
445  //
446  // kinematic of the scattered electron
447  //
448 
449  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
450 
451 
452  // Generate the electron only if with large enough range w.r.t. cuts and safety
453 
454  G4double safety = aStep.GetPostStepPoint()->GetSafety();
455 
456 
457  if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
458  {
459  G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
460  G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
461  gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
462  G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
464  aParticleChange.AddSecondary(electron);
465  // aParticleChange.ProposeLocalEnergyDeposit(0.);
467  }
468  else
469  {
471  aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE);
472  }
473 
474  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
475 
476 }
477 
478 
480  G4double sinSqrTh)
481 {
482  G4double rand1;
483  G4double rand2;
484  G4double phiProbability;
485  G4double phi;
486  G4double a, b;
487 
488  do
489  {
490  rand1 = G4UniformRand();
491  rand2 = G4UniformRand();
492  phiProbability=0.;
493  phi = twopi*rand1;
494 
495  a = 2*sinSqrTh;
496  b = energyRate + 1/energyRate;
497 
498  phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
499 
500 
501 
502  }
503  while ( rand2 > phiProbability );
504  return phi;
505 }
506 
507 
509 {
510  G4double dx = a.x();
511  G4double dy = a.y();
512  G4double dz = a.z();
513  G4double x = dx < 0.0 ? -dx : dx;
514  G4double y = dy < 0.0 ? -dy : dy;
515  G4double z = dz < 0.0 ? -dz : dz;
516  if (x < y) {
517  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
518  }else{
519  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
520  }
521 }
522 
524 {
525  G4ThreeVector d0 = direction0.unit();
526  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
527  G4ThreeVector a0 = a1.unit(); // unit vector
528 
530 
531  G4double angle = twopi*rand1; // random polar angle
532  G4ThreeVector b0 = d0.cross(a0); // cross product
533 
535 
536  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
537  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
538  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
539 
540  G4ThreeVector c0 = c.unit();
541 
542  return c0;
543 
544 }
545 
546 
548 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
549 {
550 
551  //
552  // The polarization of a photon is always perpendicular to its momentum direction.
553  // Therefore this function removes those vector component of gammaPolarization, which
554  // points in direction of gammaDirection
555  //
556  // Mathematically we search the projection of the vector a on the plane E, where n is the
557  // plains normal vector.
558  // The basic equation can be found in each geometry book (e.g. Bronstein):
559  // p = a - (a o n)/(n o n)*n
560 
561  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
562 }
563 
564 
566  G4double sinSqrTh,
567  G4double phi,
568  G4double costheta)
569 {
570  G4double rand1;
571  G4double rand2;
572  G4double cosPhi = std::cos(phi);
573  G4double sinPhi = std::sin(phi);
574  G4double sinTheta = std::sqrt(sinSqrTh);
575  G4double cosSqrPhi = cosPhi*cosPhi;
576  // G4double cossqrth = 1.-sinSqrTh;
577  // G4double sinsqrphi = sinPhi*sinPhi;
578  G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
579 
580 
581  // Determination of Theta
582 
583  // ---- MGP ---- Commented out the following 3 lines to avoid compilation
584  // warnings (unused variables)
585  // G4double thetaProbability;
586  G4double theta;
587  // G4double a, b;
588  // G4double cosTheta;
589 
590  /*
591 
592  depaola method
593 
594  do
595  {
596  rand1 = G4UniformRand();
597  rand2 = G4UniformRand();
598  thetaProbability=0.;
599  theta = twopi*rand1;
600  a = 4*normalisation*normalisation;
601  b = (epsilon + 1/epsilon) - 2;
602  thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
603  cosTheta = std::cos(theta);
604  }
605  while ( rand2 > thetaProbability );
606 
607  G4double cosBeta = cosTheta;
608 
609  */
610 
611 
612  // Dan Xu method (IEEE TNS, 52, 1160 (2005))
613 
614  rand1 = G4UniformRand();
615  rand2 = G4UniformRand();
616 
617  if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
618  {
619  if (rand2<0.5)
620  theta = pi/2.0;
621  else
622  theta = 3.0*pi/2.0;
623  }
624  else
625  {
626  if (rand2<0.5)
627  theta = 0;
628  else
629  theta = pi;
630  }
631  G4double cosBeta = std::cos(theta);
632  G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
633 
634  G4ThreeVector gammaPolarization1;
635 
636  G4double xParallel = normalisation*cosBeta;
637  G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
638  G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
639  G4double xPerpendicular = 0.;
640  G4double yPerpendicular = (costheta)*sinBeta/normalisation;
641  G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
642 
643  G4double xTotal = (xParallel + xPerpendicular);
644  G4double yTotal = (yParallel + yPerpendicular);
645  G4double zTotal = (zParallel + zPerpendicular);
646 
647  gammaPolarization1.setX(xTotal);
648  gammaPolarization1.setY(yTotal);
649  gammaPolarization1.setZ(zTotal);
650 
651  return gammaPolarization1;
652 
653 }
654 
655 
657  G4ThreeVector& direction1,
658  G4ThreeVector& polarization0,
659  G4ThreeVector& polarization1)
660 {
661  // direction0 is the original photon direction ---> z
662  // polarization0 is the original photon polarization ---> x
663  // need to specify y axis in the real reference frame ---> y
664  G4ThreeVector Axis_Z0 = direction0.unit();
665  G4ThreeVector Axis_X0 = polarization0.unit();
666  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
667 
668  G4double direction_x = direction1.getX();
669  G4double direction_y = direction1.getY();
670  G4double direction_z = direction1.getZ();
671 
672  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
673  G4double polarization_x = polarization1.getX();
674  G4double polarization_y = polarization1.getY();
675  G4double polarization_z = polarization1.getZ();
676 
677  polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
678 
679 }
680 
681 
683 {
684  return ( &particle == G4Gamma::Gamma() );
685 }
686 
687 
689  G4double,
691 {
693  G4double energy = photon->GetKineticEnergy();
694  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
695  size_t materialIndex = couple->GetIndex();
696  G4double meanFreePath;
697  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
698  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
699  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
700  return meanFreePath;
701 }
702 
703 
704 
705 
706 
707 
708 
709 
710 
711 
712 
713 
714 
715