| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ScreeningMottCrossSection.cc 26 // G4ScreeningMottCrossSection.cc
27 //-------------------------------------------- 27 //-------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class header file 29 // GEANT4 Class header file
30 // 30 //
31 // File name: G4ScreeningMottCrossSection 31 // File name: G4ScreeningMottCrossSection
32 // 32 //
33 // Author: Cristina Consolandi 33 // Author: Cristina Consolandi
34 // 34 //
35 // Creation date: 20.10.2011 35 // Creation date: 20.10.2011
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 27-05-2012 Added Analytic Fitting to the Mo 38 // 27-05-2012 Added Analytic Fitting to the Mott Cross Section by means of G4MottCoefficients class.
39 // 39 //
40 // 40 //
41 // Class Description: 41 // Class Description:
42 // Computation of electron Coulomb Scattering 42 // Computation of electron Coulomb Scattering Cross Section.
43 // Suitable for high energy electrons and lig 43 // Suitable for high energy electrons and light target materials.
44 // 44 //
45 // Reference: 45 // Reference:
46 // M.J. Boschini et al. 46 // M.J. Boschini et al.
47 // "Non Ionizing Energy Loss induced by El 47 // "Non Ionizing Energy Loss induced by Electrons in the Space Environment"
48 // Proc. of the 13th International Confer 48 // Proc. of the 13th International Conference on Particle Physics and Advanced Technology
49 // (13th ICPPAT, Como 3-7/10/2011), World 49 // (13th ICPPAT, Como 3-7/10/2011), World Scientific (Singapore).
50 // Available at: http://arxiv.org/abs/1111.40 50 // Available at: http://arxiv.org/abs/1111.4042v4
51 // 51 //
52 // 1) Mott Differential Cross Section App 52 // 1) Mott Differential Cross Section Approximation:
53 // For Target material up to Z=92 (U): 53 // For Target material up to Z=92 (U):
54 // As described in http://arxiv.org/ab 54 // As described in http://arxiv.org/abs/1111.4042v4
55 // par. 2.1 , eq. (16)-(17) 55 // par. 2.1 , eq. (16)-(17)
56 // Else (Z>92): 56 // Else (Z>92):
57 // W. A. McKinley and H. Fashbach, Phys. R 57 // W. A. McKinley and H. Fashbach, Phys. Rev. 74, (1948) 1759.
58 // 2) Screening coefficient: 58 // 2) Screening coefficient:
59 // vomn G. Moliere, Z. Naturforsh A2 (194 59 // vomn G. Moliere, Z. Naturforsh A2 (1947), 133-145; A3 (1948), 78.
60 // 3) Nuclear Form Factor: 60 // 3) Nuclear Form Factor:
61 // A.V. Butkevich et al. Nucl. Instr. Met 61 // A.V. Butkevich et al. Nucl. Instr. Meth. A488 (2002), 282-294.
62 // 62 //
63 // ------------------------------------------- 63 // ---------------------------------------------------------------------------
64 // 64 //
65 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 66
67 #include "G4ScreeningMottCrossSection.hh" 67 #include "G4ScreeningMottCrossSection.hh"
68 #include "G4PhysicalConstants.hh" 68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh" 69 #include "G4SystemOfUnits.hh"
70 #include "Randomize.hh" 70 #include "Randomize.hh"
71 #include "G4Proton.hh" 71 #include "G4Proton.hh"
72 #include "G4LossTableManager.hh" 72 #include "G4LossTableManager.hh"
73 #include "G4NucleiProperties.hh" 73 #include "G4NucleiProperties.hh"
74 #include "G4Element.hh" 74 #include "G4Element.hh"
75 #include "G4UnitsTable.hh" 75 #include "G4UnitsTable.hh"
76 #include "G4NistManager.hh" 76 #include "G4NistManager.hh"
77 #include "G4ThreeVector.hh" 77 #include "G4ThreeVector.hh"
78 #include "G4Pow.hh" 78 #include "G4Pow.hh"
79 #include "G4MottData.hh" 79 #include "G4MottData.hh"
80 80
81 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 82
83 static G4double invlog10 = 1./std::log(10.); 83 static G4double invlog10 = 1./std::log(10.);
84 84
85 static const G4double angle[DIMMOTT]={ 85 static const G4double angle[DIMMOTT]={
86 1e-07,1.02329e-07,1.04713e-07,1.07152e-07,1.09 86 1e-07,1.02329e-07,1.04713e-07,1.07152e-07,1.09648e-07,1.12202e-07,1.14815e-07,1.1749e-07,1.20226e-07,1.23027e-07,1.25893e-07,1.28825e-07,1.31826e-07,1.34896e-07,1.38038e-07,1.41254e-07,1.44544e-07,1.47911e-07,1.51356e-07,1.54882e-07,1.58489e-07,1.62181e-07,1.65959e-07,1.69824e-07,1.7378e-07,1.77828e-07,1.8197e-07,1.86209e-07,1.90546e-07,1.94984e-07,1.99526e-07,2.04174e-07,2.0893e-07,2.13796e-07,2.18776e-07,2.23872e-07,2.29087e-07,2.34423e-07,2.39883e-07,2.45471e-07,2.51189e-07,2.5704e-07,2.63027e-07,2.69153e-07,2.75423e-07,2.81838e-07,2.88403e-07,2.95121e-07,3.01995e-07,3.0903e-07,3.16228e-07,3.23594e-07,3.31131e-07,3.38844e-07,3.46737e-07,3.54813e-07,3.63078e-07,3.71535e-07,3.80189e-07,3.89045e-07,3.98107e-07,4.0738e-07,4.16869e-07,4.2658e-07,4.36516e-07,4.46684e-07,4.57088e-07,4.67735e-07,4.7863e-07,4.89779e-07,5.01187e-07,5.12861e-07,5.24807e-07,5.37032e-07,5.49541e-07,5.62341e-07,5.7544e-07,5.88844e-07,6.0256e-07,6.16595e-07,6.30957e-07,6.45654e-07,6.60693e-07,6.76083e-07,6.91831e-07,7.07946e-07,7.24436e-07,7.4131e-07,7.58578e-07,7.76247e-07,7.94328e-07,8.12831e-07,8.31764e-07,8.51138e-07,8.70964e-07,8.91251e-07,9.12011e-07,9.33254e-07,9.54993e-07,9.77237e-07,1e-06,1.02329e-06,1.04713e-06,1.07152e-06,1.09648e-06,1.12202e-06,1.14815e-06,1.1749e-06,1.20226e-06,1.23027e-06,1.25893e-06,1.28825e-06,1.31826e-06,1.34896e-06,1.38038e-06,1.41254e-06,1.44544e-06,1.47911e-06,1.51356e-06,1.54882e-06,1.58489e-06,1.62181e-06,1.65959e-06,1.69824e-06,1.7378e-06,1.77828e-06,1.8197e-06,1.86209e-06,1.90546e-06,1.94984e-06,1.99526e-06,2.04174e-06,2.0893e-06,2.13796e-06,2.18776e-06,2.23872e-06,2.29087e-06,2.34423e-06,2.39883e-06,2.45471e-06,2.51189e-06,2.5704e-06,2.63027e-06,2.69153e-06,2.75423e-06,2.81838e-06,2.88403e-06,2.95121e-06,3.01995e-06,3.0903e-06,3.16228e-06,3.23594e-06,3.31131e-06,3.38844e-06,3.46737e-06,3.54813e-06,3.63078e-06,3.71535e-06,3.80189e-06,3.89045e-06,3.98107e-06,4.0738e-06,4.16869e-06,4.2658e-06,4.36516e-06,4.46684e-06,4.57088e-06,4.67735e-06,4.7863e-06,4.89779e-06,5.01187e-06,5.12861e-06,5.24807e-06,5.37032e-06,5.49541e-06,5.62341e-06,5.7544e-06,5.88844e-06,6.0256e-06,6.16595e-06,6.30957e-06,6.45654e-06,6.60693e-06,6.76083e-06,6.91831e-06,7.07946e-06,7.24436e-06,7.4131e-06,7.58578e-06,7.76247e-06,7.94328e-06,8.12831e-06,8.31764e-06,8.51138e-06,8.70964e-06,8.91251e-06,9.12011e-06,9.33254e-06,9.54993e-06,9.77237e-06,1e-05,1.02329e-05,1.04713e-05,1.07152e-05,1.09648e-05,1.12202e-05,1.14815e-05,1.1749e-05,1.20226e-05,1.23027e-05,1.25893e-05,1.28825e-05,1.31826e-05,1.34896e-05,1.38038e-05,1.41254e-05,1.44544e-05,1.47911e-05,1.51356e-05,1.54882e-05,1.58489e-05,1.62181e-05,1.65959e-05,1.69824e-05,1.7378e-05,1.77828e-05,1.8197e-05,1.86209e-05,1.90546e-05,1.94984e-05,1.99526e-05,2.04174e-05,2.0893e-05,2.13796e-05,2.18776e-05,2.23872e-05,2.29087e-05,2.34423e-05,2.39883e-05,2.45471e-05,2.51189e-05,2.5704e-05,2.63027e-05,2.69153e-05,2.75423e-05,2.81838e-05,2.88403e-05,2.95121e-05,3.01995e-05,3.0903e-05,3.16228e-05,3.23594e-05,3.31131e-05,3.38844e-05,3.46737e-05,3.54813e-05,3.63078e-05,3.71535e-05,3.80189e-05,3.89045e-05,3.98107e-05,4.0738e-05,4.16869e-05,4.2658e-05,4.36516e-05,4.46684e-05,4.57088e-05,4.67735e-05,4.7863e-05,4.89779e-05,5.01187e-05,5.12861e-05,5.24807e-05,5.37032e-05,5.49541e-05,5.62341e-05,5.7544e-05,5.88844e-05,6.0256e-05,6.16595e-05,6.30957e-05,6.45654e-05,6.60693e-05,6.76083e-05,6.91831e-05,7.07946e-05,7.24436e-05,7.4131e-05,7.58578e-05,7.76247e-05,7.94328e-05,8.12831e-05,8.31764e-05,8.51138e-05,8.70964e-05,8.91251e-05,9.12011e-05,9.33254e-05,9.54993e-05,9.77237e-05,0.0001,0.000102329,0.000104713,0.000107152,0.000109648,0.000112202,0.000114815,0.00011749,0.000120226,0.000123027,0.000125893,0.000128825,0.000131826,0.000134896,0.000138038,0.000141254,0.000144544,0.000147911,0.000151356,0.000154882,0.000158489,0.000162181,0.000165959,0.000169824,0.00017378,0.000177828,0.00018197,0.000186209,0.000190546,0.000194984,0.000199526,0.000204174,0.00020893,0.000213796,0.000218776,0.000223872,0.000229087,0.000234423,0.000239883,0.000245471,0.000251189,0.00025704,0.000263027,0.000269153,0.000275423,0.000281838,0.000288403,0.000295121,0.000301995,0.00030903,0.000316228,0.000323594,0.000331131,0.000338844,0.000346737,0.000354813,0.000363078,0.000371535,0.000380189,0.000389045,0.000398107,0.00040738,0.000416869,0.00042658,0.000436516,0.000446684,0.000457088,0.000467735,0.00047863,0.000489779,0.000501187,0.000512861,0.000524807,0.000537032,0.000549541,0.000562341,0.00057544,0.000588844,0.00060256,0.000616595,0.000630957,0.000645654,0.000660693,0.000676083,0.000691831,0.000707946,0.000724436,0.00074131,0.000758578,0.000776247,0.000794328,0.000812831,0.000831764,0.000851138,0.000870964,0.000891251,0.000912011,0.000933254,0.000954993,0.000977237,0.001,0.00102329,0.00104713,0.00107152,0.00109648,0.00112202,0.00114815,0.0011749,0.00120226,0.00123027,0.00125893,0.00128825,0.00131826,0.00134896,0.00138038,0.00141254,0.00144544,0.00147911,0.00151356,0.00154882,0.00158489,0.00162181,0.00165959,0.00169824,0.0017378,0.00177828,0.0018197,0.00186209,0.00190546,0.00194984,0.00199526,0.00204174,0.0020893,0.00213796,0.00218776,0.00223872,0.00229087,0.00234423,0.00239883,0.00245471,0.00251189,0.0025704,0.00263027,0.00269153,0.00275423,0.00281838,0.00288403,0.00295121,0.00301995,0.0030903,0.00316228,0.00323594,0.00331131,0.00338844,0.00346737,0.00354813,0.00363078,0.00371535,0.00380189,0.00389045,0.00398107,0.0040738,0.00416869,0.0042658,0.00436516,0.00446684,0.00457088,0.00467735,0.0047863,0.00489779,0.00501187,0.00512861,0.00524807,0.00537032,0.00549541,0.00562341,0.0057544,0.00588844,0.0060256,0.00616595,0.00630957,0.00645654,0.00660693,0.00676083,0.00691831,0.00707946,0.00724436,0.0074131,0.00758578,0.00776247,0.00794328,0.00812831,0.00831764,0.00851138,0.00870964,0.00891251,0.00912011,0.00933254,0.00954993,0.00977237,0.01,0.0102329,0.0104713,0.0107152,0.0109648,0.0112202,0.0114815,0.011749,0.0120226,0.0123027,0.0125893,0.0128825,0.0131826,0.0134896,0.0138038,0.0141254,0.0144544,0.0147911,0.0151356,0.0154882,0.0158489,0.0162181,0.0165959,0.0169824,0.017378,0.0177828,0.018197,0.0186209,0.0190546,0.0194984,0.0199526,0.0204174,0.020893,0.0213796,0.0218776,0.0223872,0.0229087,0.0234423,0.0239883,0.0245471,0.0251189,0.025704,0.0263027,0.0269153,0.0275423,0.0281838,0.0288403,0.0295121,0.0301995,0.030903,0.0316228,0.0323594,0.0331131,0.0338844,0.0346737,0.0354813,0.0363078,0.0371535,0.0380189,0.0389045,0.0398107,0.040738,0.0416869,0.042658,0.0436516,0.0446684,0.0457088,0.0467735,0.047863,0.0489779,0.0501187,0.0512861,0.0524807,0.0537032,0.0549541,0.0562341,0.057544,0.0588844,0.060256,0.0616595,0.0630957,0.0645654,0.0660693,0.0676083,0.0691831,0.0707946,0.0724436,0.074131,0.0758578,0.0776247,0.0794328,0.0812831,0.0831764,0.0851138,0.0870964,0.0891251,0.0912011,0.0933254,0.0954993,0.0977237,0.1,0.102329,0.104713,0.107152,0.109648,0.112202,0.114815,0.11749,0.120226,0.123027,0.125893,0.128825,0.131826,0.134896,0.138038,0.141254,0.144544,0.147911,0.151356,0.154882,0.158489,0.162181,0.165959,0.169824,0.17378,0.177828,0.18197,0.186209,0.190546,0.194984,0.199526,0.204174,0.20893,0.213796,0.218776,0.223872,0.229087,0.234423,0.239883,0.245471,0.251189,0.25704,0.263027,0.269153,0.275423,0.281838,0.288403,0.295121,0.301995,0.30903,0.316228,0.323594,0.331131,0.338844,0.346737,0.354813,0.363078,0.371535,0.380189,0.389045,0.398107,0.40738,0.416869,0.42658,0.436516,0.446684,0.457088,0.467735,0.47863,0.489779,0.501187,0.512861,0.524807,0.537032,0.549541,0.562341,0.57544,0.588844,0.60256,0.616595,0.630957,0.645654,0.660693,0.676083,0.691831,0.707946,0.724436,0.74131,0.758578,0.776247,0.794328,0.812831,0.831764,0.851138,0.870964,0.891251,0.912011,0.933254,0.954993,0.977237,1,1.02329,1.04713,1.07152,1.09648,1.12202,1.14815,1.1749,1.20226,1.23027,1.25893,1.28825,1.31826,1.34896,1.38038,1.41254,1.44544,1.47911,1.51356,1.54882,1.58489,1.62181,1.65959,1.69824,1.7378,1.77828,1.8197,1.86209,1.90546,1.94984,1.99526,2.04174,2.0893,2.13796,2.18776,2.23872,2.29087,2.34423,2.39883,2.45471,2.51189,2.5704,2.63027,2.69153,2.75423,2.81838,2.88403,2.95121,3.01995,3.0903};
87 87
88 //using namespace std; 88 //using namespace std;
89 89
90 G4ScreeningMottCrossSection::G4ScreeningMottCr 90 G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
91 cosThetaMin(1.0), 91 cosThetaMin(1.0),
92 cosThetaMax(-1.0), 92 cosThetaMax(-1.0),
93 alpha(fine_structure_const), 93 alpha(fine_structure_const),
94 htc2(hbarc_squared), 94 htc2(hbarc_squared),
95 e2(electron_mass_c2*classic_electr_radius) 95 e2(electron_mass_c2*classic_electr_radius)
96 { 96 {
97 fTotalCross=0; 97 fTotalCross=0;
98 98
99 fNistManager = G4NistManager::Instance(); 99 fNistManager = G4NistManager::Instance();
100 fG4pow = G4Pow::GetInstance(); 100 fG4pow = G4Pow::GetInstance();
101 particle=nullptr; 101 particle=nullptr;
102 102
103 spin = mass = mu_rel = 0.0; 103 spin = mass = mu_rel = 0.0;
104 tkinLab = momLab2 = invbetaLab2 = 0.0; 104 tkinLab = momLab2 = invbetaLab2 = 0.0;
105 tkin = mom2 = invbeta2 = beta = gamma = 0.0; 105 tkin = mom2 = invbeta2 = beta = gamma = 0.0;
106 106
107 targetMass = As = etag = ecut = 0.0; 107 targetMass = As = etag = ecut = 0.0;
108 108
109 targetZ = targetA = 0; 109 targetZ = targetA = 0;
110 110
111 cosTetMinNuc = cosTetMaxNuc = 0.0; 111 cosTetMinNuc = cosTetMaxNuc = 0.0;
112 } 112 }
113 113
114 //....Ooooo0ooooo........oooOO0OOooo........oo 114 //....Ooooo0ooooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 115
116 G4ScreeningMottCrossSection::~G4ScreeningMottC << 116 G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
>> 117 {}
117 118
118 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 120
120 void G4ScreeningMottCrossSection::Initialise(c 121 void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
121 G 122 G4double cosThetaLim)
122 { 123 {
123 SetupParticle(p); 124 SetupParticle(p);
124 tkin = mom2 = 0.0; 125 tkin = mom2 = 0.0;
125 ecut = etag = DBL_MAX; 126 ecut = etag = DBL_MAX;
126 particle = p; 127 particle = p;
127 cosThetaMin = cosThetaLim; 128 cosThetaMin = cosThetaLim;
128 } 129 }
129 130
130 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131 132
132 void G4ScreeningMottCrossSection::SetupKinemat 133 void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4int Z)
133 { 134 {
134 //...Target 135 //...Target
135 const G4int iz = std::min(92, Z); 136 const G4int iz = std::min(92, Z);
136 const G4int ia = G4lrint(fNistManager->GetAt 137 const G4int ia = G4lrint(fNistManager->GetAtomicMassAmu(iz));
137 138
138 targetZ = iz; 139 targetZ = iz;
139 targetA = ia; 140 targetA = ia;
140 targetMass = G4NucleiProperties::GetNuclearM 141 targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
141 142
142 //G4cout<<"......... targetA "<< targetA <<G 143 //G4cout<<"......... targetA "<< targetA <<G4endl;
143 //G4cout<<"......... targetMass "<< targetMa 144 //G4cout<<"......... targetMass "<< targetMass/MeV <<G4endl;
144 145
145 // incident particle lab 146 // incident particle lab
146 tkinLab = ekin; 147 tkinLab = ekin;
147 momLab2 = tkinLab*(tkinLab + 2.0*mass); 148 momLab2 = tkinLab*(tkinLab + 2.0*mass);
148 invbetaLab2 = 1.0 + mass*mass/momLab2; 149 invbetaLab2 = 1.0 + mass*mass/momLab2;
149 150
150 const G4double etot = tkinLab + mass; 151 const G4double etot = tkinLab + mass;
151 const G4double ptot = std::sqrt(momLab2); 152 const G4double ptot = std::sqrt(momLab2);
152 const G4double m12 = mass*mass; 153 const G4double m12 = mass*mass;
153 154
154 // relativistic reduced mass from publucatio 155 // relativistic reduced mass from publucation
155 // A.P. Martynenko, R.N. Faustov, Teoret. ma 156 // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
156 157
157 //incident particle & target nucleus 158 //incident particle & target nucleus
158 const G4double Ecm = std::sqrt(m12 + targetM 159 const G4double Ecm = std::sqrt(m12 + targetMass*targetMass + 2.0*etot*targetMass);
159 mu_rel = mass*targetMass/Ecm; 160 mu_rel = mass*targetMass/Ecm;
160 const G4double momCM= ptot*targetMass/Ecm; 161 const G4double momCM= ptot*targetMass/Ecm;
161 // relative system 162 // relative system
162 mom2 = momCM*momCM; 163 mom2 = momCM*momCM;
163 const G4double x = mu_rel*mu_rel/mom2; 164 const G4double x = mu_rel*mu_rel/mom2;
164 invbeta2 = 1.0 + x; 165 invbeta2 = 1.0 + x;
165 tkin = momCM*std::sqrt(invbeta2) - mu_rel;// 166 tkin = momCM*std::sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
166 const G4double beta2 = 1./invbeta2; 167 const G4double beta2 = 1./invbeta2;
167 beta = std::sqrt(beta2) ; 168 beta = std::sqrt(beta2) ;
168 const G4double gamma2= invbeta2/x; 169 const G4double gamma2= invbeta2/x;
169 gamma = std::sqrt(gamma2); 170 gamma = std::sqrt(gamma2);
170 171
171 //Thomas-Fermi screening length 172 //Thomas-Fermi screening length
172 const G4double alpha2 = alpha*alpha; 173 const G4double alpha2 = alpha*alpha;
173 const G4double aU = 0.88534*CLHEP::Bohr_radi 174 const G4double aU = 0.88534*CLHEP::Bohr_radius/fG4pow->Z13(targetZ);
174 const G4double twoR2 = aU*aU; 175 const G4double twoR2 = aU*aU;
175 As = 0.25*htc2*(1.13 + 3.76*targetZ*targetZ* 176 As = 0.25*htc2*(1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2)/(twoR2*mom2);
176 177
177 //Integration Angles definition 178 //Integration Angles definition
178 cosTetMinNuc = cosThetaMin; 179 cosTetMinNuc = cosThetaMin;
179 cosTetMaxNuc = cosThetaMax; 180 cosTetMaxNuc = cosThetaMax;
180 } 181 }
181 182
182 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 184
184 G4double G4ScreeningMottCrossSection::FormFact 185 G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double sin2t2)
185 { 186 {
186 G4double M=targetMass; 187 G4double M=targetMass;
187 G4double E=tkinLab; 188 G4double E=tkinLab;
188 G4double Etot=E+mass; 189 G4double Etot=E+mass;
189 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ 190 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
190 G4double T=Tmax*sin2t2; 191 G4double T=Tmax*sin2t2;
191 G4double q2=T*(T+2.*M); 192 G4double q2=T*(T+2.*M);
192 q2/=htc2;//1/cm2 193 q2/=htc2;//1/cm2
193 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ 194 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targetA)*0.27)*CLHEP::cm;
194 G4double xN= (RN*RN*q2); 195 G4double xN= (RN*RN*q2);
195 G4double den=(1.+xN/12.); 196 G4double den=(1.+xN/12.);
196 G4double FN=1./(den*den); 197 G4double FN=1./(den*den);
197 G4double form2=(FN*FN); 198 G4double form2=(FN*FN);
198 199
199 return form2; 200 return form2;
200 } 201 }
201 202
202 //....oooOO0OOooo........oooOO0OOooo........oo 203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 204
204 G4double G4ScreeningMottCrossSection::FormFact 205 G4double G4ScreeningMottCrossSection::FormFactor2Gauss(G4double sin2t2)
205 { 206 {
206 G4double M=targetMass; 207 G4double M=targetMass;
207 G4double E=tkinLab; 208 G4double E=tkinLab;
208 G4double Etot=E+mass; 209 G4double Etot=E+mass;
209 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ 210 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
210 G4double T=Tmax*sin2t2; 211 G4double T=Tmax*sin2t2;
211 G4double q2=T*(T+2.*M); 212 G4double q2=T*(T+2.*M);
212 q2/=htc2;//1/cm2 213 q2/=htc2;//1/cm2
213 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targ 214 G4double RN=1.27e-13*G4Exp(fG4pow->logZ(targetA)*0.27)*CLHEP::cm;
214 G4double xN= (RN*RN*q2); 215 G4double xN= (RN*RN*q2);
215 216
216 G4double expo=(-xN/6.); 217 G4double expo=(-xN/6.);
217 G4double FN=G4Exp(expo); 218 G4double FN=G4Exp(expo);
218 219
219 G4double form2=(FN*FN); 220 G4double form2=(FN*FN);
220 221
221 return form2; 222 return form2;
222 } 223 }
223 224
224 //....oooOO0OOooo........oooOO0OOooo........oo 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 226
226 G4double G4ScreeningMottCrossSection::FormFact 227 G4double G4ScreeningMottCrossSection::FormFactor2UniformHelm(G4double sin2t2)
227 { 228 {
228 G4double M=targetMass; 229 G4double M=targetMass;
229 G4double E=tkinLab; 230 G4double E=tkinLab;
230 G4double Etot=E+mass; 231 G4double Etot=E+mass;
231 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+ 232 G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
232 G4double T=Tmax*sin2t2; 233 G4double T=Tmax*sin2t2;
233 G4double q2=T*(T+2.*M); 234 G4double q2=T*(T+2.*M);
234 q2=q2/(htc2*0.01);//1/cm2 235 q2=q2/(htc2*0.01);//1/cm2
235 236
236 G4double q=std::sqrt(q2); 237 G4double q=std::sqrt(q2);
237 G4double R0=1.2E-13*fG4pow->Z13(targetA); 238 G4double R0=1.2E-13*fG4pow->Z13(targetA);
238 G4double R1=2.0E-13; 239 G4double R1=2.0E-13;
239 240
240 G4double x0=q*R0; 241 G4double x0=q*R0;
241 G4double F0=(3./fG4pow->powN(x0,3))*(std::si 242 G4double F0=(3./fG4pow->powN(x0,3))*(std::sin(x0)-x0*std::cos(x0));
242 243
243 G4double x1=q*R1; 244 G4double x1=q*R1;
244 G4double F1=(3./fG4pow->powN(x1,3))*(std::si 245 G4double F1=(3./fG4pow->powN(x1,3))*(std::sin(x1)-x1*std::cos(x1));
245 246
246 G4double F=F0*F1; 247 G4double F=F0*F1;
247 248
248 G4double form2=(F*F); 249 G4double form2=(F*F);
249 250
250 return form2; 251 return form2;
251 } 252 }
252 253
253 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 255
255 G4double G4ScreeningMottCrossSection::McFcorre 256 G4double G4ScreeningMottCrossSection::McFcorrection(G4double sin2t2)
256 { 257 {
257 const G4double sint = std::sqrt(sin2t2); 258 const G4double sint = std::sqrt(sin2t2);
258 return 1.-beta*beta*sin2t2 + targetZ*alpha*b 259 return 1.-beta*beta*sin2t2 + targetZ*alpha*beta*pi*sint*(1.-sint);
259 } 260 }
260 261
261 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 263
263 G4double G4ScreeningMottCrossSection::RatioMot 264 G4double G4ScreeningMottCrossSection::RatioMottRutherford(G4double angles)
264 { 265 {
265 return RatioMottRutherfordCosT(std::sqrt(1. 266 return RatioMottRutherfordCosT(std::sqrt(1. - std::cos(angles)));
266 } 267 }
267 268
268 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
269 270
270 G4double G4ScreeningMottCrossSection::RatioMot 271 G4double G4ScreeningMottCrossSection::RatioMottRutherfordCosT(G4double fcost)
271 { 272 {
272 G4double R(0.); 273 G4double R(0.);
273 const G4double shift = 0.7181228; 274 const G4double shift = 0.7181228;
274 G4double beta0 = beta - shift; 275 G4double beta0 = beta - shift;
275 276
276 G4double b0 = 1.0; 277 G4double b0 = 1.0;
277 G4double b[6]; 278 G4double b[6];
278 for(G4int i=0; i<6; ++i) { 279 for(G4int i=0; i<6; ++i) {
279 b[i] = b0; 280 b[i] = b0;
280 b0 *= beta0; 281 b0 *= beta0;
281 } 282 }
282 283
283 b0 = 1.0; 284 b0 = 1.0;
284 for(G4int j=0; j<=4; ++j) { 285 for(G4int j=0; j<=4; ++j) {
285 G4double a = 0.0; 286 G4double a = 0.0;
286 for(G4int k=0; k<=5; ++k){ 287 for(G4int k=0; k<=5; ++k){
287 a += fMottCoef[targetZ][j][k]*b[k]; 288 a += fMottCoef[targetZ][j][k]*b[k];
288 } 289 }
289 R += a*b0; 290 R += a*b0;
290 b0 *= fcost; 291 b0 *= fcost;
291 } 292 }
292 return R; 293 return R;
293 } 294 }
294 295
295 //....oooOO0OOooo........oooOO0OOooo........oo 296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 297
297 G4double G4ScreeningMottCrossSection::GetTrans 298 G4double G4ScreeningMottCrossSection::GetTransitionRandom()
298 { 299 {
299 G4double x = G4Log(tkinLab)*invlog10; 300 G4double x = G4Log(tkinLab)*invlog10;
300 G4double res(0.), x0(1.0); 301 G4double res(0.), x0(1.0);
301 for(G4int i=0; i<11; ++i) { 302 for(G4int i=0; i<11; ++i) {
302 res += fPRM[targetZ][i]*x0; 303 res += fPRM[targetZ][i]*x0;
303 x0 *= x; 304 x0 *= x;
304 } 305 }
305 //G4cout << "Z= " << targetZ << " x0= " << x 306 //G4cout << "Z= " << targetZ << " x0= " << x0 << " res= " << res << G4endl;
306 return res; 307 return res;
307 } 308 }
308 309
309 //....oooOO0OOooo........oooOO0OOooo........oo 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310 311
311 G4double 312 G4double
312 G4ScreeningMottCrossSection::DifferentialXSect 313 G4ScreeningMottCrossSection::DifferentialXSection(G4int idx, G4int form)
313 { 314 {
314 const G4double ang = angle[idx]; 315 const G4double ang = angle[idx];
315 const G4double z1 = 1.0 - std::cos(ang); 316 const G4double z1 = 1.0 - std::cos(ang);
316 G4double step; 317 G4double step;
317 if(0 == idx) { 318 if(0 == idx) {
318 step = 0.5*(angle[1] + angle[0]); 319 step = 0.5*(angle[1] + angle[0]);
319 } else if(DIMMOTT == idx + 1) { 320 } else if(DIMMOTT == idx + 1) {
320 step = CLHEP::pi - 0.5*(angle[DIMMOTT-2] + 321 step = CLHEP::pi - 0.5*(angle[DIMMOTT-2] + angle[DIMMOTT-1]);
321 } else { 322 } else {
322 step = 0.5*(angle[idx+1] - angle[idx-1]); 323 step = 0.5*(angle[idx+1] - angle[idx-1]);
323 } 324 }
324 325
325 G4double F2 = 1.; 326 G4double F2 = 1.;
326 switch (form) { 327 switch (form) {
327 case 1: 328 case 1:
328 F2 = FormFactor2ExpHof(z1*0.5); 329 F2 = FormFactor2ExpHof(z1*0.5);
329 break; 330 break;
330 case 2: 331 case 2:
331 F2 = FormFactor2Gauss(z1*0.5); 332 F2 = FormFactor2Gauss(z1*0.5);
332 break; 333 break;
333 case 3: 334 case 3:
334 F2 = FormFactor2UniformHelm(z1*0.5); 335 F2 = FormFactor2UniformHelm(z1*0.5);
335 break; 336 break;
336 } 337 }
337 338
338 const G4double R = RatioMottRutherfordCosT(s 339 const G4double R = RatioMottRutherfordCosT(std::sqrt(z1));
339 340
340 G4double den = 2.*As + z1; 341 G4double den = 2.*As + z1;
341 G4double func = 1./(den*den); 342 G4double func = 1./(den*den);
342 343
343 G4double fatt = (G4double)targetZ/(mu_rel*ga 344 G4double fatt = (G4double)targetZ/(mu_rel*gamma*beta*beta);
344 G4double sigma= e2*e2*fatt*fatt*func; 345 G4double sigma= e2*e2*fatt*fatt*func;
345 G4double pi2sintet = CLHEP::twopi*std::sqrt( 346 G4double pi2sintet = CLHEP::twopi*std::sqrt(z1*(2. - z1));
346 347
347 G4double Xsec = std::max(pi2sintet*F2*R*sigm 348 G4double Xsec = std::max(pi2sintet*F2*R*sigma*step, 0.);
348 return Xsec; 349 return Xsec;
349 } 350 }
350 351
351 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 353
353 G4double 354 G4double
354 G4ScreeningMottCrossSection::NuclearCrossSecti 355 G4ScreeningMottCrossSection::NuclearCrossSection(G4int form, G4int fast)
355 { 356 {
356 fTotalCross=0.; 357 fTotalCross=0.;
357 if(cosTetMaxNuc >= cosTetMinNuc) return fTot 358 if(cosTetMaxNuc >= cosTetMinNuc) return fTotalCross;
358 if(0 == cross.size()) { cross.resize(DIMMOTT 359 if(0 == cross.size()) { cross.resize(DIMMOTT, 0.0); }
359 360
360 //G4cout<<"MODEL: "<<fast<<G4endl; 361 //G4cout<<"MODEL: "<<fast<<G4endl;
361 //************ PRECISE COMPUTATION 362 //************ PRECISE COMPUTATION
362 if(fast == 0){ 363 if(fast == 0){
363 364
364 for(G4int i=0; i<DIMMOTT; ++i){ 365 for(G4int i=0; i<DIMMOTT; ++i){
365 G4double y = DifferentialXSection(i, for 366 G4double y = DifferentialXSection(i, form);
366 fTotalCross += y; 367 fTotalCross += y;
367 cross[i] = fTotalCross; 368 cross[i] = fTotalCross;
368 if(fTotalCross*1.e-9 > y) { 369 if(fTotalCross*1.e-9 > y) {
369 for(G4int j=i+1; j<DIMMOTT; ++j) { 370 for(G4int j=i+1; j<DIMMOTT; ++j) {
370 cross[j] = fTotalCross; 371 cross[j] = fTotalCross;
371 } 372 }
372 break; 373 break;
373 } 374 }
374 } 375 }
375 //**************** FAST COMPUTATION 376 //**************** FAST COMPUTATION
376 } else if(fast == 1) { 377 } else if(fast == 1) {
377 378
378 G4double p0 = electron_mass_c2*classic_ele 379 G4double p0 = electron_mass_c2*classic_electr_radius;
379 G4double coeff = twopi*p0*p0; 380 G4double coeff = twopi*p0*p0;
380 381
381 G4double fac = coeff*targetZ*targetZ*invbe 382 G4double fac = coeff*targetZ*targetZ*invbeta2/mom2;
382 383
383 G4double x = 1.0 - cosTetMinNuc; 384 G4double x = 1.0 - cosTetMinNuc;
384 G4double x1 = x + 2*As; 385 G4double x1 = x + 2*As;
385 386
386 // scattering with nucleus 387 // scattering with nucleus
387 fTotalCross = fac*(cosTetMinNuc - cosTetMa 388 fTotalCross = fac*(cosTetMinNuc - cosTetMaxNuc)/
388 (x1*(1.0 - cosTetMaxNuc + 2*As)); 389 (x1*(1.0 - cosTetMaxNuc + 2*As));
389 } 390 }
390 /* 391 /*
391 G4cout << "Energy(MeV): " << tkinLab/CLHEP:: 392 G4cout << "Energy(MeV): " << tkinLab/CLHEP::MeV
392 << " Total Cross(mb): " << fTotalCros 393 << " Total Cross(mb): " << fTotalCross
393 << " form= " << form << " fast= " << 394 << " form= " << form << " fast= " << fast << G4endl;
394 */ 395 */
395 return fTotalCross; 396 return fTotalCross;
396 } 397 }
397 398
398 //....oooOO0OOooo........oooOO0OOooo........oo 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399 400
400 G4double 401 G4double
401 G4ScreeningMottCrossSection::GetScatteringAngl 402 G4ScreeningMottCrossSection::GetScatteringAngle(G4int form, G4int fast)
402 { 403 {
403 G4double scattangle = 0.0; 404 G4double scattangle = 0.0;
404 G4double r = G4UniformRand(); 405 G4double r = G4UniformRand();
405 //************ PRECISE COMPUTATION 406 //************ PRECISE COMPUTATION
406 if(fast == 0){ 407 if(fast == 0){
407 //G4cout << "r= " << r << " tot= " << fTot 408 //G4cout << "r= " << r << " tot= " << fTotalCross << G4endl;
408 r *= fTotalCross; 409 r *= fTotalCross;
409 for(G4int i=0; i<DIMMOTT; ++i){ 410 for(G4int i=0; i<DIMMOTT; ++i){
410 //G4cout << i << ". r= " << r << " xs= " 411 //G4cout << i << ". r= " << r << " xs= " << cross[i] << G4endl;
411 if(r <= cross[i]) { 412 if(r <= cross[i]) {
412 scattangle = ComputeAngle(i, r); 413 scattangle = ComputeAngle(i, r);
413 break; 414 break;
414 } 415 }
415 } 416 }
416 417
417 //**************** FAST COMPUTATION 418 //**************** FAST COMPUTATION
418 } else if(fast == 1) { 419 } else if(fast == 1) {
419 420
420 G4double limit = GetTransitionRandom(); 421 G4double limit = GetTransitionRandom();
421 if(limit > 0.0) { 422 if(limit > 0.0) {
422 G4double Sz = 2*As; 423 G4double Sz = 2*As;
423 G4double x = Sz-((Sz*(2+Sz))/(Sz+2*limit 424 G4double x = Sz-((Sz*(2+Sz))/(Sz+2*limit))+1.0;
424 //G4cout << "limit= " << limit << " Sz= 425 //G4cout << "limit= " << limit << " Sz= " << Sz << " x= " << x << G4endl;
425 G4double angle_limit = (std::abs(x) < 1. 426 G4double angle_limit = (std::abs(x) < 1.0) ? std::acos(x) : 0.0;
426 //G4cout<<"ANGLE LIMIT: "<<angle_limit<< 427 //G4cout<<"ANGLE LIMIT: "<<angle_limit<<" x= " << x << G4endl;
427 428
428 if(r > limit && angle_limit != 0.0) { 429 if(r > limit && angle_limit != 0.0) {
429 x = Sz-((Sz*(2+Sz))/(Sz+2*r))+1.0; 430 x = Sz-((Sz*(2+Sz))/(Sz+2*r))+1.0;
430 scattangle = (x >= 1.0) ? 0.0 : ((x <= -1.0) 431 scattangle = (x >= 1.0) ? 0.0 : ((x <= -1.0) ? pi : std::acos(x));
431 //G4cout<<"FAST: scattangle= "<< scattangle 432 //G4cout<<"FAST: scattangle= "<< scattangle <<G4endl;
432 } 433 }
433 } else { 434 } else {
434 //G4cout<<"SLOW: total= "<<fTotalCross<< 435 //G4cout<<"SLOW: total= "<<fTotalCross<< G4endl;
435 r *= fTotalCross; 436 r *= fTotalCross;
436 G4double tot = 0.0; 437 G4double tot = 0.0;
437 for(G4int i=0; i<DIMMOTT; ++i) { 438 for(G4int i=0; i<DIMMOTT; ++i) {
438 G4double xs = DifferentialXSection(i, form); 439 G4double xs = DifferentialXSection(i, form);
439 tot += xs; 440 tot += xs;
440 cross[i] = tot; 441 cross[i] = tot;
441 if(r <= tot) { 442 if(r <= tot) {
442 scattangle = ComputeAngle(i, r); 443 scattangle = ComputeAngle(i, r);
443 break; 444 break;
444 } 445 }
445 } 446 }
446 } 447 }
447 } 448 }
448 //****************************************** 449 //************************************************
449 //G4cout<<"Energy(MeV): "<<tkinLab/MeV<<" SC 450 //G4cout<<"Energy(MeV): "<<tkinLab/MeV<<" SCATTANGLE: "<<scattangle<<G4endl;
450 return scattangle; 451 return scattangle;
451 } 452 }
452 453
453 G4double G4ScreeningMottCrossSection::ComputeA 454 G4double G4ScreeningMottCrossSection::ComputeAngle(G4int i, G4double& r)
454 { 455 {
455 G4double x1, x2, y; 456 G4double x1, x2, y;
456 if(0 == i) { 457 if(0 == i) {
457 x1 = 0.0; 458 x1 = 0.0;
458 x2 = 0.5*(angle[0] + angle[1]); 459 x2 = 0.5*(angle[0] + angle[1]);
459 y = cross[0]; 460 y = cross[0];
460 } else if(i == DIMMOTT-1) { 461 } else if(i == DIMMOTT-1) {
461 x1 = 0.5*(angle[i] + angle[i-1]); 462 x1 = 0.5*(angle[i] + angle[i-1]);
462 x2 = CLHEP::pi; 463 x2 = CLHEP::pi;
463 y = cross[i] - cross[i-1]; 464 y = cross[i] - cross[i-1];
464 r -= cross[i-1]; 465 r -= cross[i-1];
465 } else { 466 } else {
466 x1 = 0.5*(angle[i] + angle[i-1]); 467 x1 = 0.5*(angle[i] + angle[i-1]);
467 x2 = 0.5*(angle[i] + angle[i+1]); 468 x2 = 0.5*(angle[i] + angle[i+1]);
468 y = cross[i] - cross[i-1]; 469 y = cross[i] - cross[i-1];
469 r -= cross[i-1]; 470 r -= cross[i-1];
470 } 471 }
471 //G4cout << i << ". r= " << r << " y= " << y 472 //G4cout << i << ". r= " << r << " y= " << y
472 // << " x1= " << " x2= " << x2 << G4endl; 473 // << " x1= " << " x2= " << x2 << G4endl;
473 return x1 + (x2 - x1)*r/y; 474 return x1 + (x2 - x1)*r/y;
474 } 475 }
475 476