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 | // 17.07.06 V. Grichine - first implementation |
---|
28 | // 22.01.07 V.Ivanchenko - add interface with Z and A |
---|
29 | // 05.03.07 V.Ivanchenko - add IfZAApplicable |
---|
30 | // |
---|
31 | |
---|
32 | #include "G4GlauberGribovCrossSection.hh" |
---|
33 | |
---|
34 | #include "G4ParticleTable.hh" |
---|
35 | #include "G4IonTable.hh" |
---|
36 | #include "G4ParticleDefinition.hh" |
---|
37 | |
---|
38 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
39 | // |
---|
40 | // |
---|
41 | |
---|
42 | const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionTot[93] = { |
---|
43 | |
---|
44 | 1.0, 1.0, 1.118517e+00, 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00, |
---|
45 | 1.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00, |
---|
46 | 1.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00, |
---|
47 | 1.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00, |
---|
48 | 1.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00, |
---|
49 | 1.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00, |
---|
50 | 1.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00, |
---|
51 | 1.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00, |
---|
52 | 1.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00, |
---|
53 | 1.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00, |
---|
54 | 1.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00, |
---|
55 | 1.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00, |
---|
56 | 1.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00, |
---|
57 | 1.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00, |
---|
58 | 1.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00, |
---|
59 | 1.037075e+00, 1.034721e+00 |
---|
60 | |
---|
61 | }; |
---|
62 | |
---|
63 | const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionIn[93] = { |
---|
64 | |
---|
65 | 1.0, 1.0, 1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00, |
---|
66 | 1.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00, |
---|
67 | 1.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00, |
---|
68 | 1.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00, |
---|
69 | 1.131515e+00, 1.130338e+00, 1.134171e+00, 1.139206e+00, 1.141474e+00, 1.142189e+00, |
---|
70 | 1.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00, |
---|
71 | 1.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00, |
---|
72 | 1.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00, |
---|
73 | 1.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00, |
---|
74 | 1.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00, |
---|
75 | 1.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00, |
---|
76 | 1.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00, |
---|
77 | 1.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00, |
---|
78 | 1.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00, |
---|
79 | 1.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00, |
---|
80 | 1.046435e+00, 1.042614e+00 |
---|
81 | |
---|
82 | }; |
---|
83 | |
---|
84 | const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionTot[93] = { |
---|
85 | |
---|
86 | 1.0, 1.0, |
---|
87 | 1.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00, |
---|
88 | 1.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00, |
---|
89 | 1.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00, |
---|
90 | 1.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00, |
---|
91 | 1.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00, |
---|
92 | 1.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00, |
---|
93 | 1.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00, |
---|
94 | 1.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00, |
---|
95 | 1.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00, |
---|
96 | 1.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00, |
---|
97 | 1.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00, |
---|
98 | 1.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00, |
---|
99 | 1.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00, |
---|
100 | 1.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00, |
---|
101 | 1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00, |
---|
102 | 1.034720e+00 |
---|
103 | |
---|
104 | }; |
---|
105 | |
---|
106 | const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionIn[93] = { |
---|
107 | |
---|
108 | 1.0, 1.0, |
---|
109 | 1.167419e+00, 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00, |
---|
110 | 1.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00, |
---|
111 | 1.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00, |
---|
112 | 1.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00, |
---|
113 | 1.130337e+00, 1.134170e+00, 1.139205e+00, 1.141472e+00, 1.142188e+00, 1.140724e+00, |
---|
114 | 1.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00, |
---|
115 | 1.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00, |
---|
116 | 1.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00, |
---|
117 | 1.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00, |
---|
118 | 1.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00, |
---|
119 | 1.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00, |
---|
120 | 1.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00, |
---|
121 | 1.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00, |
---|
122 | 1.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00, |
---|
123 | 1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00, |
---|
124 | 1.042613e+00 |
---|
125 | |
---|
126 | }; |
---|
127 | |
---|
128 | |
---|
129 | const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionTot[93] = { |
---|
130 | |
---|
131 | 1.0, 1.0, |
---|
132 | 1.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00, |
---|
133 | 1.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00, |
---|
134 | 1.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00, |
---|
135 | 1.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00, |
---|
136 | 1.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00, |
---|
137 | 1.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00, |
---|
138 | 1.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00, |
---|
139 | 1.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00, |
---|
140 | 1.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00, |
---|
141 | 1.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00, |
---|
142 | 1.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00, |
---|
143 | 1.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00, |
---|
144 | 1.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00, |
---|
145 | 1.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00, |
---|
146 | 1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00, |
---|
147 | 1.152974e+00 |
---|
148 | |
---|
149 | }; |
---|
150 | |
---|
151 | const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionIn[93] = { |
---|
152 | |
---|
153 | 1.0, 1.0, |
---|
154 | 1.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.044495e+00, 1.062622e+00, |
---|
155 | 1.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.065100e+00, |
---|
156 | 1.070480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00, |
---|
157 | 1.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00, |
---|
158 | 1.137312e+00, 1.126263e+00, 1.126459e+00, 1.115191e+00, 1.116986e+00, 1.117184e+00, |
---|
159 | 1.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00, |
---|
160 | 1.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00, |
---|
161 | 1.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00, |
---|
162 | 1.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00, |
---|
163 | 1.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00, |
---|
164 | 1.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00, |
---|
165 | 1.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00, |
---|
166 | 1.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00, |
---|
167 | 1.071107e+00, 1.069955e+00, 1.064856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00, |
---|
168 | 1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00, |
---|
169 | 1.059658e+00 |
---|
170 | |
---|
171 | }; |
---|
172 | |
---|
173 | |
---|
174 | const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionTot[93] = { |
---|
175 | |
---|
176 | 1.0, 1.0, |
---|
177 | 1.075927e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00, |
---|
178 | 1.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00, |
---|
179 | 1.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00, |
---|
180 | 1.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00, |
---|
181 | 1.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00, |
---|
182 | 1.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00, |
---|
183 | 1.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00, |
---|
184 | 1.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00, |
---|
185 | 1.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00, |
---|
186 | 1.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00, |
---|
187 | 1.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00, |
---|
188 | 1.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00, |
---|
189 | 1.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00, |
---|
190 | 1.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00, |
---|
191 | 1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00, |
---|
192 | 1.157267e+00 |
---|
193 | }; |
---|
194 | |
---|
195 | |
---|
196 | const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionIn[93] = { |
---|
197 | |
---|
198 | 1.0, 1.0, |
---|
199 | 1.140246e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.044514e+00, 1.062628e+00, |
---|
200 | 1.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00, |
---|
201 | 1.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00, |
---|
202 | 1.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00, |
---|
203 | 1.138502e+00, 1.127678e+00, 1.127244e+00, 1.116634e+00, 1.118347e+00, 1.118988e+00, |
---|
204 | 1.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00, |
---|
205 | 1.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00, |
---|
206 | 1.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00, |
---|
207 | 1.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00, |
---|
208 | 1.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00, |
---|
209 | 1.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00, |
---|
210 | 1.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00, |
---|
211 | 1.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00, |
---|
212 | 1.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00, |
---|
213 | 1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00, |
---|
214 | 1.062699e+00 |
---|
215 | |
---|
216 | }; |
---|
217 | |
---|
218 | |
---|
219 | |
---|
220 | |
---|
221 | //////////////////////////////////////////////////////////////////////////////// |
---|
222 | // |
---|
223 | // |
---|
224 | |
---|
225 | G4GlauberGribovCrossSection::G4GlauberGribovCrossSection() |
---|
226 | : fUpperLimit( 100000 * GeV ), |
---|
227 | fLowerLimit( 3 * GeV ), |
---|
228 | fRadiusConst( 1.08*fermi ) // 1.1, 1.3 ? |
---|
229 | { |
---|
230 | theGamma = G4Gamma::Gamma(); |
---|
231 | theProton = G4Proton::Proton(); |
---|
232 | theNeutron = G4Neutron::Neutron(); |
---|
233 | theAProton = G4AntiProton::AntiProton(); |
---|
234 | theANeutron = G4AntiNeutron::AntiNeutron(); |
---|
235 | thePiPlus = G4PionPlus::PionPlus(); |
---|
236 | thePiMinus = G4PionMinus::PionMinus(); |
---|
237 | thePiZero = G4PionZero::PionZero(); |
---|
238 | theKPlus = G4KaonPlus::KaonPlus(); |
---|
239 | theKMinus = G4KaonMinus::KaonMinus(); |
---|
240 | theK0S = G4KaonZeroShort::KaonZeroShort(); |
---|
241 | theK0L = G4KaonZeroLong::KaonZeroLong(); |
---|
242 | theL = G4Lambda::Lambda(); |
---|
243 | theAntiL = G4AntiLambda::AntiLambda(); |
---|
244 | theSPlus = G4SigmaPlus::SigmaPlus(); |
---|
245 | theASPlus = G4AntiSigmaPlus::AntiSigmaPlus(); |
---|
246 | theSMinus = G4SigmaMinus::SigmaMinus(); |
---|
247 | theASMinus = G4AntiSigmaMinus::AntiSigmaMinus(); |
---|
248 | theS0 = G4SigmaZero::SigmaZero(); |
---|
249 | theAS0 = G4AntiSigmaZero::AntiSigmaZero(); |
---|
250 | theXiMinus = G4XiMinus::XiMinus(); |
---|
251 | theXi0 = G4XiZero::XiZero(); |
---|
252 | theAXiMinus = G4AntiXiMinus::AntiXiMinus(); |
---|
253 | theAXi0 = G4AntiXiZero::AntiXiZero(); |
---|
254 | theOmega = G4OmegaMinus::OmegaMinus(); |
---|
255 | theAOmega = G4AntiOmegaMinus::AntiOmegaMinus(); |
---|
256 | theD = G4Deuteron::Deuteron(); |
---|
257 | theT = G4Triton::Triton(); |
---|
258 | theA = G4Alpha::Alpha(); |
---|
259 | theHe3 = G4He3::He3(); |
---|
260 | } |
---|
261 | |
---|
262 | /////////////////////////////////////////////////////////////////////////////////////// |
---|
263 | // |
---|
264 | // |
---|
265 | |
---|
266 | G4GlauberGribovCrossSection::~G4GlauberGribovCrossSection() |
---|
267 | { |
---|
268 | } |
---|
269 | |
---|
270 | |
---|
271 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
272 | // |
---|
273 | // |
---|
274 | |
---|
275 | |
---|
276 | G4bool |
---|
277 | G4GlauberGribovCrossSection::IsApplicable(const G4DynamicParticle* aDP, |
---|
278 | const G4Element* anElement) |
---|
279 | { |
---|
280 | return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN()); |
---|
281 | } |
---|
282 | |
---|
283 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
284 | // |
---|
285 | // |
---|
286 | |
---|
287 | G4bool |
---|
288 | G4GlauberGribovCrossSection::IsZAApplicable(const G4DynamicParticle* aDP, |
---|
289 | G4double Z, G4double) |
---|
290 | { |
---|
291 | G4bool applicable = false; |
---|
292 | // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber(); |
---|
293 | G4double kineticEnergy = aDP->GetKineticEnergy(); |
---|
294 | |
---|
295 | const G4ParticleDefinition* theParticle = aDP->GetDefinition(); |
---|
296 | |
---|
297 | if ( ( kineticEnergy >= fLowerLimit && |
---|
298 | Z > 1.5 && // >= He |
---|
299 | ( theParticle == theAProton || |
---|
300 | theParticle == theGamma || |
---|
301 | theParticle == theKPlus || |
---|
302 | theParticle == theKMinus || |
---|
303 | theParticle == theSMinus) ) || |
---|
304 | |
---|
305 | ( kineticEnergy >= fLowerLimit && |
---|
306 | Z > 1.5 && // >= He |
---|
307 | ( theParticle == theProton || |
---|
308 | theParticle == theNeutron || |
---|
309 | theParticle == thePiPlus || |
---|
310 | theParticle == thePiMinus ) ) ) applicable = true; |
---|
311 | |
---|
312 | return applicable; |
---|
313 | } |
---|
314 | |
---|
315 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
316 | // |
---|
317 | // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to |
---|
318 | // Glauber model with Gribov correction calculated in the dipole approximation on |
---|
319 | // light cone. Gaussian density helps to calculate rest integrals of the model. |
---|
320 | // [1] B.Z. Kopeliovich, nucl-th/0306044 |
---|
321 | |
---|
322 | |
---|
323 | G4double G4GlauberGribovCrossSection:: |
---|
324 | GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T) |
---|
325 | { |
---|
326 | return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T); |
---|
327 | } |
---|
328 | |
---|
329 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
330 | // |
---|
331 | // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to |
---|
332 | // Glauber model with Gribov correction calculated in the dipole approximation on |
---|
333 | // light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model. |
---|
334 | // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above |
---|
335 | |
---|
336 | |
---|
337 | |
---|
338 | G4double G4GlauberGribovCrossSection:: |
---|
339 | GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double Z, G4double A, G4double) |
---|
340 | { |
---|
341 | G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio; |
---|
342 | G4double R = GetNucleusRadius(A); |
---|
343 | |
---|
344 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
345 | |
---|
346 | if( theParticle == theProton || |
---|
347 | theParticle == theNeutron || |
---|
348 | theParticle == thePiPlus || |
---|
349 | theParticle == thePiMinus ) |
---|
350 | { |
---|
351 | sigma = GetHadronNucleonXscNS(aParticle, A, Z); |
---|
352 | cofInelastic = 2.4; |
---|
353 | cofTotal = 2.0; |
---|
354 | } |
---|
355 | else |
---|
356 | { |
---|
357 | sigma = GetHadronNucleonXscPDG(aParticle, A, Z); |
---|
358 | cofInelastic = 2.2; |
---|
359 | cofTotal = 2.0; |
---|
360 | } |
---|
361 | // cofInelastic = 2.0; |
---|
362 | |
---|
363 | |
---|
364 | nucleusSquare = cofTotal*pi*R*R; // basically 2piRR |
---|
365 | ratio = sigma/nucleusSquare; |
---|
366 | |
---|
367 | xsection = nucleusSquare*std::log( 1. + ratio ); |
---|
368 | |
---|
369 | xsection *= GetParticleBarCorTot(theParticle, Z); |
---|
370 | |
---|
371 | fTotalXsc = xsection; |
---|
372 | |
---|
373 | |
---|
374 | |
---|
375 | fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; |
---|
376 | |
---|
377 | fInelasticXsc *= GetParticleBarCorIn(theParticle, Z); |
---|
378 | |
---|
379 | fElasticXsc = fTotalXsc - fInelasticXsc; |
---|
380 | |
---|
381 | |
---|
382 | G4double difratio = ratio/(1.+ratio); |
---|
383 | |
---|
384 | fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) ); |
---|
385 | |
---|
386 | |
---|
387 | sigma = GetHNinelasticXsc(aParticle, A, Z); |
---|
388 | ratio = sigma/nucleusSquare; |
---|
389 | |
---|
390 | fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; |
---|
391 | |
---|
392 | if (fElasticXsc < 0.) fElasticXsc = 0.; |
---|
393 | |
---|
394 | return xsection; |
---|
395 | } |
---|
396 | |
---|
397 | ////////////////////////////////////////////////////////////////////////// |
---|
398 | // |
---|
399 | // Return single-diffraction/inelastic cross-section ratio |
---|
400 | |
---|
401 | G4double G4GlauberGribovCrossSection:: |
---|
402 | GetRatioSD(const G4DynamicParticle* aParticle, G4double A, G4double Z) |
---|
403 | { |
---|
404 | G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio; |
---|
405 | G4double R = GetNucleusRadius(A); |
---|
406 | |
---|
407 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
408 | |
---|
409 | if( theParticle == theProton || |
---|
410 | theParticle == theNeutron || |
---|
411 | theParticle == thePiPlus || |
---|
412 | theParticle == thePiMinus ) |
---|
413 | { |
---|
414 | sigma = GetHadronNucleonXscNS(aParticle, A, Z); |
---|
415 | cofInelastic = 2.4; |
---|
416 | cofTotal = 2.0; |
---|
417 | } |
---|
418 | else |
---|
419 | { |
---|
420 | sigma = GetHadronNucleonXscPDG(aParticle, A, Z); |
---|
421 | cofInelastic = 2.2; |
---|
422 | cofTotal = 2.0; |
---|
423 | } |
---|
424 | nucleusSquare = cofTotal*pi*R*R; // basically 2piRR |
---|
425 | ratio = sigma/nucleusSquare; |
---|
426 | |
---|
427 | fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; |
---|
428 | |
---|
429 | G4double difratio = ratio/(1.+ratio); |
---|
430 | |
---|
431 | fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) ); |
---|
432 | |
---|
433 | if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc; |
---|
434 | else ratio = 0.; |
---|
435 | |
---|
436 | return ratio; |
---|
437 | } |
---|
438 | |
---|
439 | ////////////////////////////////////////////////////////////////////////// |
---|
440 | // |
---|
441 | // Return suasi-elastic/inelastic cross-section ratio |
---|
442 | |
---|
443 | G4double G4GlauberGribovCrossSection:: |
---|
444 | GetRatioQE(const G4DynamicParticle* aParticle, G4double A, G4double Z) |
---|
445 | { |
---|
446 | G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio; |
---|
447 | G4double R = GetNucleusRadius(A); |
---|
448 | |
---|
449 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
450 | |
---|
451 | if( theParticle == theProton || |
---|
452 | theParticle == theNeutron || |
---|
453 | theParticle == thePiPlus || |
---|
454 | theParticle == thePiMinus ) |
---|
455 | { |
---|
456 | sigma = GetHadronNucleonXscNS(aParticle, A, Z); |
---|
457 | cofInelastic = 2.4; |
---|
458 | cofTotal = 2.0; |
---|
459 | } |
---|
460 | else |
---|
461 | { |
---|
462 | sigma = GetHadronNucleonXscPDG(aParticle, A, Z); |
---|
463 | cofInelastic = 2.2; |
---|
464 | cofTotal = 2.0; |
---|
465 | } |
---|
466 | nucleusSquare = cofTotal*pi*R*R; // basically 2piRR |
---|
467 | ratio = sigma/nucleusSquare; |
---|
468 | |
---|
469 | fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; |
---|
470 | |
---|
471 | sigma = GetHNinelasticXsc(aParticle, A, Z); |
---|
472 | ratio = sigma/nucleusSquare; |
---|
473 | |
---|
474 | fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic; |
---|
475 | |
---|
476 | if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc; |
---|
477 | else ratio = 0.; |
---|
478 | if ( ratio < 0. ) ratio = 0.; |
---|
479 | |
---|
480 | return ratio; |
---|
481 | } |
---|
482 | |
---|
483 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
484 | // |
---|
485 | // Returns hadron-nucleon Xsc according to differnt parametrisations: |
---|
486 | // [2] E. Levin, hep-ph/9710546 |
---|
487 | // [3] U. Dersch, et al, hep-ex/9910052 |
---|
488 | // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 |
---|
489 | |
---|
490 | G4double |
---|
491 | G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, |
---|
492 | const G4Element* anElement ) |
---|
493 | { |
---|
494 | G4double At = anElement->GetN(); // number of nucleons |
---|
495 | G4double Zt = anElement->GetZ(); // number of protons |
---|
496 | |
---|
497 | |
---|
498 | return GetHadronNucleonXsc( aParticle, At, Zt ); |
---|
499 | } |
---|
500 | |
---|
501 | |
---|
502 | |
---|
503 | |
---|
504 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
505 | // |
---|
506 | // Returns hadron-nucleon Xsc according to differnt parametrisations: |
---|
507 | // [2] E. Levin, hep-ph/9710546 |
---|
508 | // [3] U. Dersch, et al, hep-ex/9910052 |
---|
509 | // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 |
---|
510 | |
---|
511 | G4double |
---|
512 | G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, |
---|
513 | G4double At, G4double Zt ) |
---|
514 | { |
---|
515 | G4double xsection; |
---|
516 | |
---|
517 | |
---|
518 | G4double targ_mass = G4ParticleTable::GetParticleTable()-> |
---|
519 | GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) ); |
---|
520 | |
---|
521 | targ_mass = 0.939*GeV; // ~mean neutron and proton ??? |
---|
522 | |
---|
523 | G4double proj_mass = aParticle->GetMass(); |
---|
524 | G4double proj_momentum = aParticle->GetMomentum().mag(); |
---|
525 | G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); |
---|
526 | |
---|
527 | sMand /= GeV*GeV; // in GeV for parametrisation |
---|
528 | proj_momentum /= GeV; |
---|
529 | |
---|
530 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
531 | |
---|
532 | |
---|
533 | if(theParticle == theGamma) |
---|
534 | { |
---|
535 | xsection = At*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525)); |
---|
536 | } |
---|
537 | else if(theParticle == theNeutron) // as proton ??? |
---|
538 | { |
---|
539 | xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); |
---|
540 | } |
---|
541 | else if(theParticle == theProton) |
---|
542 | { |
---|
543 | xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); |
---|
544 | // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) ); |
---|
545 | // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) ); |
---|
546 | } |
---|
547 | else if(theParticle == theAProton) |
---|
548 | { |
---|
549 | xsection = At*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525)); |
---|
550 | } |
---|
551 | else if(theParticle == thePiPlus) |
---|
552 | { |
---|
553 | xsection = At*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525)); |
---|
554 | } |
---|
555 | else if(theParticle == thePiMinus) |
---|
556 | { |
---|
557 | // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) ); |
---|
558 | xsection = At*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525)); |
---|
559 | } |
---|
560 | else if(theParticle == theKPlus) |
---|
561 | { |
---|
562 | xsection = At*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525)); |
---|
563 | } |
---|
564 | else if(theParticle == theKMinus) |
---|
565 | { |
---|
566 | xsection = At*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525)); |
---|
567 | } |
---|
568 | else // as proton ??? |
---|
569 | { |
---|
570 | xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525)); |
---|
571 | } |
---|
572 | xsection *= millibarn; |
---|
573 | return xsection; |
---|
574 | } |
---|
575 | |
---|
576 | |
---|
577 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
578 | // |
---|
579 | // Returns hadron-nucleon Xsc according to PDG parametrisation (2005): |
---|
580 | // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf |
---|
581 | |
---|
582 | G4double |
---|
583 | G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, |
---|
584 | const G4Element* anElement ) |
---|
585 | { |
---|
586 | G4double At = anElement->GetN(); // number of nucleons |
---|
587 | G4double Zt = anElement->GetZ(); // number of protons |
---|
588 | |
---|
589 | |
---|
590 | return GetHadronNucleonXscPDG( aParticle, At, Zt ); |
---|
591 | } |
---|
592 | |
---|
593 | |
---|
594 | |
---|
595 | |
---|
596 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
597 | // |
---|
598 | // Returns hadron-nucleon Xsc according to PDG parametrisation (2005): |
---|
599 | // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf |
---|
600 | // At = number of nucleons, Zt = number of protons |
---|
601 | |
---|
602 | G4double |
---|
603 | G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, |
---|
604 | G4double At, G4double Zt ) |
---|
605 | { |
---|
606 | G4double xsection; |
---|
607 | |
---|
608 | G4double Nt = At-Zt; // number of neutrons |
---|
609 | if (Nt < 0.) Nt = 0.; |
---|
610 | |
---|
611 | |
---|
612 | G4double targ_mass = G4ParticleTable::GetParticleTable()-> |
---|
613 | GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) ); |
---|
614 | |
---|
615 | targ_mass = 0.939*GeV; // ~mean neutron and proton ??? |
---|
616 | |
---|
617 | G4double proj_mass = aParticle->GetMass(); |
---|
618 | G4double proj_momentum = aParticle->GetMomentum().mag(); |
---|
619 | |
---|
620 | G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); |
---|
621 | |
---|
622 | sMand /= GeV*GeV; // in GeV for parametrisation |
---|
623 | |
---|
624 | // General PDG fit constants |
---|
625 | |
---|
626 | G4double s0 = 5.38*5.38; // in Gev^2 |
---|
627 | G4double eta1 = 0.458; |
---|
628 | G4double eta2 = 0.458; |
---|
629 | G4double B = 0.308; |
---|
630 | |
---|
631 | |
---|
632 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
633 | |
---|
634 | |
---|
635 | if(theParticle == theNeutron) // proton-neutron fit |
---|
636 | { |
---|
637 | xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
638 | + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); |
---|
639 | xsection += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
640 | + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn |
---|
641 | } |
---|
642 | else if(theParticle == theProton) |
---|
643 | { |
---|
644 | |
---|
645 | xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
646 | + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); |
---|
647 | |
---|
648 | xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
649 | + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); |
---|
650 | } |
---|
651 | else if(theParticle == theAProton) |
---|
652 | { |
---|
653 | xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
654 | + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); |
---|
655 | |
---|
656 | xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
657 | + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); |
---|
658 | } |
---|
659 | else if(theParticle == thePiPlus) |
---|
660 | { |
---|
661 | xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.) |
---|
662 | + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2)); |
---|
663 | } |
---|
664 | else if(theParticle == thePiMinus) |
---|
665 | { |
---|
666 | xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.) |
---|
667 | + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2)); |
---|
668 | } |
---|
669 | else if(theParticle == theKPlus) |
---|
670 | { |
---|
671 | xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) |
---|
672 | + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2)); |
---|
673 | |
---|
674 | xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) |
---|
675 | + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2)); |
---|
676 | } |
---|
677 | else if(theParticle == theKMinus) |
---|
678 | { |
---|
679 | xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) |
---|
680 | + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2)); |
---|
681 | |
---|
682 | xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) |
---|
683 | + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2)); |
---|
684 | } |
---|
685 | else if(theParticle == theSMinus) |
---|
686 | { |
---|
687 | xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.) |
---|
688 | - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2)); |
---|
689 | } |
---|
690 | else if(theParticle == theGamma) // modify later on |
---|
691 | { |
---|
692 | xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.) |
---|
693 | + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2)); |
---|
694 | |
---|
695 | } |
---|
696 | else // as proton ??? |
---|
697 | { |
---|
698 | xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
699 | + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); |
---|
700 | |
---|
701 | xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
702 | + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); |
---|
703 | } |
---|
704 | xsection *= millibarn; // parametrised in mb |
---|
705 | return xsection; |
---|
706 | } |
---|
707 | |
---|
708 | |
---|
709 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
710 | // |
---|
711 | // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of |
---|
712 | // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database |
---|
713 | |
---|
714 | G4double |
---|
715 | G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle, |
---|
716 | const G4Element* anElement ) |
---|
717 | { |
---|
718 | G4double At = anElement->GetN(); // number of nucleons |
---|
719 | G4double Zt = anElement->GetZ(); // number of protons |
---|
720 | |
---|
721 | |
---|
722 | return GetHadronNucleonXscNS( aParticle, At, Zt ); |
---|
723 | } |
---|
724 | |
---|
725 | |
---|
726 | |
---|
727 | |
---|
728 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
729 | // |
---|
730 | // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of |
---|
731 | // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database |
---|
732 | |
---|
733 | G4double |
---|
734 | G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle, |
---|
735 | G4double At, G4double Zt ) |
---|
736 | { |
---|
737 | G4double xsection(0), Delta, A0, B0; |
---|
738 | G4double hpXsc(0); |
---|
739 | G4double hnXsc(0); |
---|
740 | |
---|
741 | G4double Nt = At-Zt; // number of neutrons |
---|
742 | if (Nt < 0.) Nt = 0.; |
---|
743 | |
---|
744 | |
---|
745 | G4double targ_mass = G4ParticleTable::GetParticleTable()-> |
---|
746 | GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) ); |
---|
747 | |
---|
748 | targ_mass = 0.939*GeV; // ~mean neutron and proton ??? |
---|
749 | |
---|
750 | G4double proj_mass = aParticle->GetMass(); |
---|
751 | G4double proj_energy = aParticle->GetTotalEnergy(); |
---|
752 | G4double proj_momentum = aParticle->GetMomentum().mag(); |
---|
753 | |
---|
754 | G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); |
---|
755 | |
---|
756 | sMand /= GeV*GeV; // in GeV for parametrisation |
---|
757 | proj_momentum /= GeV; |
---|
758 | proj_energy /= GeV; |
---|
759 | proj_mass /= GeV; |
---|
760 | |
---|
761 | // General PDG fit constants |
---|
762 | |
---|
763 | G4double s0 = 5.38*5.38; // in Gev^2 |
---|
764 | G4double eta1 = 0.458; |
---|
765 | G4double eta2 = 0.458; |
---|
766 | G4double B = 0.308; |
---|
767 | |
---|
768 | |
---|
769 | const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); |
---|
770 | |
---|
771 | |
---|
772 | if(theParticle == theNeutron) |
---|
773 | { |
---|
774 | if( proj_momentum >= 10.) |
---|
775 | // if( proj_momentum >= 2.) |
---|
776 | { |
---|
777 | Delta = 1.; |
---|
778 | |
---|
779 | if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy; |
---|
780 | |
---|
781 | if(proj_momentum >= 10.) |
---|
782 | { |
---|
783 | B0 = 7.5; |
---|
784 | A0 = 100. - B0*std::log(3.0e7); |
---|
785 | |
---|
786 | xsection = A0 + B0*std::log(proj_energy) - 11 |
---|
787 | + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+ |
---|
788 | 0.93827*0.93827,-0.165); // mb |
---|
789 | } |
---|
790 | xsection *= Zt + Nt; |
---|
791 | } |
---|
792 | else |
---|
793 | { |
---|
794 | // nn to be pp |
---|
795 | |
---|
796 | if( proj_momentum < 0.73 ) |
---|
797 | { |
---|
798 | hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) ); |
---|
799 | } |
---|
800 | else if( proj_momentum < 1.05 ) |
---|
801 | { |
---|
802 | hnXsc = 23 + 40*(std::log(proj_momentum/0.73))* |
---|
803 | (std::log(proj_momentum/0.73)); |
---|
804 | } |
---|
805 | else // if( proj_momentum < 10. ) |
---|
806 | { |
---|
807 | hnXsc = 39.0+ |
---|
808 | 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15); |
---|
809 | } |
---|
810 | // pn to be np |
---|
811 | |
---|
812 | if( proj_momentum < 0.8 ) |
---|
813 | { |
---|
814 | hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0); |
---|
815 | } |
---|
816 | else if( proj_momentum < 1.4 ) |
---|
817 | { |
---|
818 | hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0); |
---|
819 | } |
---|
820 | else // if( proj_momentum < 10. ) |
---|
821 | { |
---|
822 | hpXsc = 33.3+ |
---|
823 | 20.8*(std::pow(proj_momentum,2.0)-1.35)/ |
---|
824 | (std::pow(proj_momentum,2.50)+0.95); |
---|
825 | } |
---|
826 | xsection = hpXsc*Zt + hnXsc*Nt; |
---|
827 | } |
---|
828 | } |
---|
829 | else if(theParticle == theProton) |
---|
830 | { |
---|
831 | if( proj_momentum >= 10.) |
---|
832 | // if( proj_momentum >= 2.) |
---|
833 | { |
---|
834 | Delta = 1.; |
---|
835 | |
---|
836 | if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy; |
---|
837 | |
---|
838 | if(proj_momentum >= 10.) |
---|
839 | { |
---|
840 | B0 = 7.5; |
---|
841 | A0 = 100. - B0*std::log(3.0e7); |
---|
842 | |
---|
843 | xsection = A0 + B0*std::log(proj_energy) - 11 |
---|
844 | + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+ |
---|
845 | 0.93827*0.93827,-0.165); // mb |
---|
846 | } |
---|
847 | xsection *= Zt + Nt; |
---|
848 | } |
---|
849 | else |
---|
850 | { |
---|
851 | // pp |
---|
852 | |
---|
853 | if( proj_momentum < 0.73 ) |
---|
854 | { |
---|
855 | hpXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) ); |
---|
856 | } |
---|
857 | else if( proj_momentum < 1.05 ) |
---|
858 | { |
---|
859 | hpXsc = 23 + 40*(std::log(proj_momentum/0.73))* |
---|
860 | (std::log(proj_momentum/0.73)); |
---|
861 | } |
---|
862 | else // if( proj_momentum < 10. ) |
---|
863 | { |
---|
864 | hpXsc = 39.0+ |
---|
865 | 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15); |
---|
866 | } |
---|
867 | // pn to be np |
---|
868 | |
---|
869 | if( proj_momentum < 0.8 ) |
---|
870 | { |
---|
871 | hnXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0); |
---|
872 | } |
---|
873 | else if( proj_momentum < 1.4 ) |
---|
874 | { |
---|
875 | hnXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0); |
---|
876 | } |
---|
877 | else // if( proj_momentum < 10. ) |
---|
878 | { |
---|
879 | hnXsc = 33.3+ |
---|
880 | 20.8*(std::pow(proj_momentum,2.0)-1.35)/ |
---|
881 | (std::pow(proj_momentum,2.50)+0.95); |
---|
882 | } |
---|
883 | xsection = hpXsc*Zt + hnXsc*Nt; |
---|
884 | // xsection = hpXsc*(Zt + Nt); |
---|
885 | // xsection = hnXsc*(Zt + Nt); |
---|
886 | } |
---|
887 | // xsection *= 0.95; |
---|
888 | } |
---|
889 | else if(theParticle == theAProton) |
---|
890 | { |
---|
891 | xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
892 | + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2)); |
---|
893 | |
---|
894 | xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
895 | + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2)); |
---|
896 | } |
---|
897 | else if(theParticle == thePiPlus) |
---|
898 | { |
---|
899 | if(proj_momentum < 0.4) |
---|
900 | { |
---|
901 | G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085); |
---|
902 | hpXsc = Ex3+20.0; |
---|
903 | } |
---|
904 | else if(proj_momentum < 1.15) |
---|
905 | { |
---|
906 | G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75)); |
---|
907 | hpXsc = Ex4+14.0; |
---|
908 | } |
---|
909 | else if(proj_momentum < 3.5) |
---|
910 | { |
---|
911 | G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55); |
---|
912 | G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225); |
---|
913 | hpXsc = Ex1+Ex2+27.5; |
---|
914 | } |
---|
915 | else // if(proj_momentum > 3.5) // mb |
---|
916 | { |
---|
917 | hpXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43); |
---|
918 | } |
---|
919 | // pi+n = pi-p?? |
---|
920 | |
---|
921 | if(proj_momentum < 0.37) |
---|
922 | { |
---|
923 | hnXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07); |
---|
924 | } |
---|
925 | else if(proj_momentum<0.65) |
---|
926 | { |
---|
927 | hnXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48)); |
---|
928 | } |
---|
929 | else if(proj_momentum<1.3) |
---|
930 | { |
---|
931 | hnXsc = 36.1+ |
---|
932 | 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+ |
---|
933 | 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075); |
---|
934 | } |
---|
935 | else if(proj_momentum<3.0) |
---|
936 | { |
---|
937 | hnXsc = 36.1+0.079-4.313*std::log(proj_momentum)+ |
---|
938 | 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+ |
---|
939 | 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12); |
---|
940 | } |
---|
941 | else // mb |
---|
942 | { |
---|
943 | hnXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); |
---|
944 | } |
---|
945 | xsection = hpXsc*Zt + hnXsc*Nt; |
---|
946 | } |
---|
947 | else if(theParticle == thePiMinus) |
---|
948 | { |
---|
949 | // pi-n = pi+p?? |
---|
950 | |
---|
951 | if(proj_momentum < 0.4) |
---|
952 | { |
---|
953 | G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085); |
---|
954 | hnXsc = Ex3+20.0; |
---|
955 | } |
---|
956 | else if(proj_momentum < 1.15) |
---|
957 | { |
---|
958 | G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75)); |
---|
959 | hnXsc = Ex4+14.0; |
---|
960 | } |
---|
961 | else if(proj_momentum < 3.5) |
---|
962 | { |
---|
963 | G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55); |
---|
964 | G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225); |
---|
965 | hnXsc = Ex1+Ex2+27.5; |
---|
966 | } |
---|
967 | else // if(proj_momentum > 3.5) // mb |
---|
968 | { |
---|
969 | hnXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43); |
---|
970 | } |
---|
971 | // pi-p |
---|
972 | |
---|
973 | if(proj_momentum < 0.37) |
---|
974 | { |
---|
975 | hpXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07); |
---|
976 | } |
---|
977 | else if(proj_momentum<0.65) |
---|
978 | { |
---|
979 | hpXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48)); |
---|
980 | } |
---|
981 | else if(proj_momentum<1.3) |
---|
982 | { |
---|
983 | hpXsc = 36.1+ |
---|
984 | 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+ |
---|
985 | 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075); |
---|
986 | } |
---|
987 | else if(proj_momentum<3.0) |
---|
988 | { |
---|
989 | hpXsc = 36.1+0.079-4.313*std::log(proj_momentum)+ |
---|
990 | 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+ |
---|
991 | 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12); |
---|
992 | } |
---|
993 | else // mb |
---|
994 | { |
---|
995 | hpXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43); |
---|
996 | } |
---|
997 | xsection = hpXsc*Zt + hnXsc*Nt; |
---|
998 | } |
---|
999 | else if(theParticle == theKPlus) |
---|
1000 | { |
---|
1001 | xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) |
---|
1002 | + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2)); |
---|
1003 | |
---|
1004 | xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) |
---|
1005 | + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2)); |
---|
1006 | } |
---|
1007 | else if(theParticle == theKMinus) |
---|
1008 | { |
---|
1009 | xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.) |
---|
1010 | + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2)); |
---|
1011 | |
---|
1012 | xsection += Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.) |
---|
1013 | + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2)); |
---|
1014 | } |
---|
1015 | else if(theParticle == theSMinus) |
---|
1016 | { |
---|
1017 | xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.) |
---|
1018 | - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2)); |
---|
1019 | } |
---|
1020 | else if(theParticle == theGamma) // modify later on |
---|
1021 | { |
---|
1022 | xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.) |
---|
1023 | + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2)); |
---|
1024 | |
---|
1025 | } |
---|
1026 | else // as proton ??? |
---|
1027 | { |
---|
1028 | xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) |
---|
1029 | + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); |
---|
1030 | |
---|
1031 | xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) |
---|
1032 | + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2)); |
---|
1033 | } |
---|
1034 | xsection *= millibarn; // parametrised in mb |
---|
1035 | return xsection; |
---|
1036 | } |
---|
1037 | |
---|
1038 | |
---|
1039 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1040 | // |
---|
1041 | // Returns hadron-nucleon inelastic cross-section based on proper parametrisation |
---|
1042 | |
---|
1043 | G4double |
---|
1044 | G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle, |
---|
1045 | const G4Element* anElement ) |
---|
1046 | { |
---|
1047 | G4double At = anElement->GetN(); // number of nucleons |
---|
1048 | G4double Zt = anElement->GetZ(); // number of protons |
---|
1049 | |
---|
1050 | |
---|
1051 | return GetHNinelasticXsc( aParticle, At, Zt ); |
---|
1052 | } |
---|
1053 | |
---|
1054 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1055 | // |
---|
1056 | // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation |
---|
1057 | |
---|
1058 | G4double |
---|
1059 | G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle, |
---|
1060 | G4double At, G4double Zt ) |
---|
1061 | { |
---|
1062 | G4ParticleDefinition* hadron = aParticle->GetDefinition(); |
---|
1063 | G4double sumInelastic, Nt = At - Zt; |
---|
1064 | if(Nt < 0.) Nt = 0.; |
---|
1065 | |
---|
1066 | if( hadron == theKPlus ) |
---|
1067 | { |
---|
1068 | sumInelastic = GetHNinelasticXscVU(aParticle, At, Zt); |
---|
1069 | } |
---|
1070 | else |
---|
1071 | { |
---|
1072 | sumInelastic = Zt*GetHadronNucleonXscMK(aParticle, theProton); |
---|
1073 | sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron); |
---|
1074 | } |
---|
1075 | return sumInelastic; |
---|
1076 | } |
---|
1077 | |
---|
1078 | |
---|
1079 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1080 | // |
---|
1081 | // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation |
---|
1082 | |
---|
1083 | G4double |
---|
1084 | G4GlauberGribovCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle, |
---|
1085 | G4double At, G4double Zt ) |
---|
1086 | { |
---|
1087 | G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding(); |
---|
1088 | G4int absPDGcode = std::abs(PDGcode); |
---|
1089 | |
---|
1090 | G4double Elab = aParticle->GetTotalEnergy(); |
---|
1091 | // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV; |
---|
1092 | G4double Plab = aParticle->GetMomentum().mag(); |
---|
1093 | // std::sqrt(Elab * Elab - 0.88); |
---|
1094 | |
---|
1095 | Elab /= GeV; |
---|
1096 | Plab /= GeV; |
---|
1097 | |
---|
1098 | G4double LogPlab = std::log( Plab ); |
---|
1099 | G4double sqrLogPlab = LogPlab * LogPlab; |
---|
1100 | |
---|
1101 | //G4cout<<"Plab = "<<Plab<<G4endl; |
---|
1102 | |
---|
1103 | G4double NumberOfTargetProtons = Zt; |
---|
1104 | G4double NumberOfTargetNucleons = At; |
---|
1105 | G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons; |
---|
1106 | |
---|
1107 | if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.; |
---|
1108 | |
---|
1109 | G4double Xtotal, Xelastic, Xinelastic; |
---|
1110 | |
---|
1111 | if( absPDGcode > 1000 ) //------Projectile is baryon -------- |
---|
1112 | { |
---|
1113 | G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + |
---|
1114 | 0.522*sqrLogPlab - 4.51*LogPlab; |
---|
1115 | |
---|
1116 | G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + |
---|
1117 | 0.513*sqrLogPlab - 4.27*LogPlab; |
---|
1118 | |
---|
1119 | G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + |
---|
1120 | 0.169*sqrLogPlab - 1.85*LogPlab; |
---|
1121 | |
---|
1122 | G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + |
---|
1123 | 0.169*sqrLogPlab - 1.85*LogPlab; |
---|
1124 | |
---|
1125 | Xtotal = ( NumberOfTargetProtons * XtotPP + |
---|
1126 | NumberOfTargetNeutrons * XtotPN ); |
---|
1127 | |
---|
1128 | Xelastic = ( NumberOfTargetProtons * XelPP + |
---|
1129 | NumberOfTargetNeutrons * XelPN ); |
---|
1130 | } |
---|
1131 | else if( PDGcode == 211 ) //------Projectile is PionPlus ------- |
---|
1132 | { |
---|
1133 | G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) + |
---|
1134 | 0.19 *sqrLogPlab - 0.0 *LogPlab; |
---|
1135 | |
---|
1136 | G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) + |
---|
1137 | 0.456*sqrLogPlab - 4.03*LogPlab; |
---|
1138 | |
---|
1139 | G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) + |
---|
1140 | 0.079*sqrLogPlab - 0.0 *LogPlab; |
---|
1141 | |
---|
1142 | G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) + |
---|
1143 | 0.043*sqrLogPlab - 0.0 *LogPlab; |
---|
1144 | |
---|
1145 | Xtotal = ( NumberOfTargetProtons * XtotPiP + |
---|
1146 | NumberOfTargetNeutrons * XtotPiN ); |
---|
1147 | |
---|
1148 | Xelastic = ( NumberOfTargetProtons * XelPiP + |
---|
1149 | NumberOfTargetNeutrons * XelPiN ); |
---|
1150 | } |
---|
1151 | else if( PDGcode == -211 ) //------Projectile is PionMinus ------- |
---|
1152 | { |
---|
1153 | G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) + |
---|
1154 | 0.456*sqrLogPlab - 4.03*LogPlab; |
---|
1155 | |
---|
1156 | G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) + |
---|
1157 | 0.19 *sqrLogPlab - 0.0 *LogPlab; |
---|
1158 | |
---|
1159 | G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) + |
---|
1160 | 0.043*sqrLogPlab - 0.0 *LogPlab; |
---|
1161 | |
---|
1162 | G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) + |
---|
1163 | 0.079*sqrLogPlab - 0.0 *LogPlab; |
---|
1164 | |
---|
1165 | Xtotal = ( NumberOfTargetProtons * XtotPiP + |
---|
1166 | NumberOfTargetNeutrons * XtotPiN ); |
---|
1167 | |
---|
1168 | Xelastic = ( NumberOfTargetProtons * XelPiP + |
---|
1169 | NumberOfTargetNeutrons * XelPiN ); |
---|
1170 | } |
---|
1171 | else if( PDGcode == 111 ) //------Projectile is PionZero ------- |
---|
1172 | { |
---|
1173 | G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) + |
---|
1174 | 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+ |
---|
1175 | 33.0 + 14.0 *std::pow(Plab,-1.36) + |
---|
1176 | 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi- |
---|
1177 | |
---|
1178 | G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) + |
---|
1179 | 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+ |
---|
1180 | 16.4 + 19.3 *std::pow(Plab,-0.42) + |
---|
1181 | 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi- |
---|
1182 | |
---|
1183 | G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) + |
---|
1184 | 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+ |
---|
1185 | 1.76 + 11.2*std::pow(Plab,-0.64) + |
---|
1186 | 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi- |
---|
1187 | |
---|
1188 | G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) + |
---|
1189 | 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+ |
---|
1190 | 0.0 + 11.4*std::pow(Plab,-0.40) + |
---|
1191 | 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi- |
---|
1192 | |
---|
1193 | Xtotal = ( NumberOfTargetProtons * XtotPiP + |
---|
1194 | NumberOfTargetNeutrons * XtotPiN ); |
---|
1195 | |
---|
1196 | Xelastic = ( NumberOfTargetProtons * XelPiP + |
---|
1197 | NumberOfTargetNeutrons * XelPiN ); |
---|
1198 | } |
---|
1199 | else if( PDGcode == 321 ) //------Projectile is KaonPlus ------- |
---|
1200 | { |
---|
1201 | G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) + |
---|
1202 | 0.26 *sqrLogPlab - 1.0 *LogPlab; |
---|
1203 | G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) + |
---|
1204 | 0.21 *sqrLogPlab - 0.89*LogPlab; |
---|
1205 | |
---|
1206 | G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) + |
---|
1207 | 0.16 *sqrLogPlab - 1.3 *LogPlab; |
---|
1208 | |
---|
1209 | G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) + |
---|
1210 | 0.29 *sqrLogPlab - 2.4 *LogPlab; |
---|
1211 | |
---|
1212 | Xtotal = ( NumberOfTargetProtons * XtotKP + |
---|
1213 | NumberOfTargetNeutrons * XtotKN ); |
---|
1214 | |
---|
1215 | Xelastic = ( NumberOfTargetProtons * XelKP + |
---|
1216 | NumberOfTargetNeutrons * XelKN ); |
---|
1217 | } |
---|
1218 | else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------ |
---|
1219 | { |
---|
1220 | G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) + |
---|
1221 | 0.66 *sqrLogPlab - 5.6 *LogPlab; |
---|
1222 | G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) + |
---|
1223 | 0.38 *sqrLogPlab - 2.9 *LogPlab; |
---|
1224 | |
---|
1225 | G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) + |
---|
1226 | 0.29 *sqrLogPlab - 2.4 *LogPlab; |
---|
1227 | |
---|
1228 | G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) + |
---|
1229 | 0.16 *sqrLogPlab - 1.3 *LogPlab; |
---|
1230 | |
---|
1231 | Xtotal = ( NumberOfTargetProtons * XtotKP + |
---|
1232 | NumberOfTargetNeutrons * XtotKN ); |
---|
1233 | |
---|
1234 | Xelastic = ( NumberOfTargetProtons * XelKP + |
---|
1235 | NumberOfTargetNeutrons * XelKN ); |
---|
1236 | } |
---|
1237 | else if( PDGcode == 311 ) //------Projectile is KaonZero ------ |
---|
1238 | { |
---|
1239 | G4double XtotKP = ( 18.1 + 0. *std::pow(Plab, 0. ) + |
---|
1240 | 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+ |
---|
1241 | 32.1 + 0. *std::pow(Plab, 0. ) + |
---|
1242 | 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K- |
---|
1243 | |
---|
1244 | G4double XtotKN = ( 18.7 + 0. *std::pow(Plab, 0. ) + |
---|
1245 | 0.21 *sqrLogPlab - 0.89*LogPlab + //K+ |
---|
1246 | 25.2 + 0. *std::pow(Plab, 0. ) + |
---|
1247 | 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K- |
---|
1248 | |
---|
1249 | G4double XelKP = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) |
---|
1250 | + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+ |
---|
1251 | 7.3 + 0. *std::pow(Plab,-0. ) + |
---|
1252 | 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K- |
---|
1253 | |
---|
1254 | G4double XelKN = ( 7.3 + 0. *std::pow(Plab,-0. ) + |
---|
1255 | 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+ |
---|
1256 | 5.0 + 8.1*std::pow(Plab,-1.8 ) + |
---|
1257 | 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K- |
---|
1258 | |
---|
1259 | Xtotal = ( NumberOfTargetProtons * XtotKP + |
---|
1260 | NumberOfTargetNeutrons * XtotKN ); |
---|
1261 | |
---|
1262 | Xelastic = ( NumberOfTargetProtons * XelKP + |
---|
1263 | NumberOfTargetNeutrons * XelKN ); |
---|
1264 | } |
---|
1265 | else //------Projectile is undefined, Nucleon assumed |
---|
1266 | { |
---|
1267 | G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + |
---|
1268 | 0.522*sqrLogPlab - 4.51*LogPlab; |
---|
1269 | |
---|
1270 | G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + |
---|
1271 | 0.513*sqrLogPlab - 4.27*LogPlab; |
---|
1272 | |
---|
1273 | G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + |
---|
1274 | 0.169*sqrLogPlab - 1.85*LogPlab; |
---|
1275 | G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + |
---|
1276 | 0.169*sqrLogPlab - 1.85*LogPlab; |
---|
1277 | |
---|
1278 | Xtotal = ( NumberOfTargetProtons * XtotPP + |
---|
1279 | NumberOfTargetNeutrons * XtotPN ); |
---|
1280 | |
---|
1281 | Xelastic = ( NumberOfTargetProtons * XelPP + |
---|
1282 | NumberOfTargetNeutrons * XelPN ); |
---|
1283 | } |
---|
1284 | Xinelastic = Xtotal - Xelastic; |
---|
1285 | |
---|
1286 | if(Xinelastic < 0.) Xinelastic = 0.; |
---|
1287 | |
---|
1288 | return Xinelastic*= millibarn; |
---|
1289 | } |
---|
1290 | |
---|
1291 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1292 | // |
---|
1293 | // Returns hadron-nucleon cross-section based on Mikhail Kossov CHIPS parametrisation of |
---|
1294 | // data from G4QuasiFreeRatios class |
---|
1295 | |
---|
1296 | G4double |
---|
1297 | G4GlauberGribovCrossSection::GetHadronNucleonXscMK(const G4DynamicParticle* aParticle, |
---|
1298 | const G4ParticleDefinition* nucleon ) |
---|
1299 | { |
---|
1300 | G4int I = -1; |
---|
1301 | G4int PDG = aParticle->GetDefinition()->GetPDGEncoding(); |
---|
1302 | G4double totalXsc = 0; |
---|
1303 | G4double elasticXsc = 0; |
---|
1304 | G4double inelasticXsc; |
---|
1305 | // G4int absPDG = std::abs(PDG); |
---|
1306 | |
---|
1307 | G4double p = aParticle->GetMomentum().mag()/GeV; |
---|
1308 | |
---|
1309 | G4bool F = false; |
---|
1310 | if(nucleon == theProton) F = true; |
---|
1311 | else if(nucleon == theNeutron) F = false; |
---|
1312 | else |
---|
1313 | { |
---|
1314 | G4cout << "nucleon is not proton or neutron, return xsc for proton" << G4endl; |
---|
1315 | F = true; |
---|
1316 | } |
---|
1317 | |
---|
1318 | G4bool kfl = true; // Flag of K0/aK0 oscillation |
---|
1319 | G4bool kf = false; |
---|
1320 | |
---|
1321 | if( PDG == 130 || PDG == 310 ) |
---|
1322 | { |
---|
1323 | kf = true; |
---|
1324 | if( G4UniformRand() > .5 ) kfl = false; |
---|
1325 | } |
---|
1326 | if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) I = 0; // pp/nn |
---|
1327 | else if( (PDG == 2112 && F) || (PDG == 2212 && !F) ) I = 1; // np/pn |
---|
1328 | |
---|
1329 | else if( (PDG == -211 && F) || (PDG == 211 && !F) ) I = 2; // pimp/pipn |
---|
1330 | else if( (PDG == 211 && F) || (PDG ==-211 && !F) ) I = 3; // pipp/pimn |
---|
1331 | |
---|
1332 | else if( PDG == -321 || PDG == -311 || ( kf && !kfl ) ) I = 4; // KmN/K0N |
---|
1333 | else if( PDG == 321 || PDG == 311 || ( kf && kfl ) ) I = 5; // KpN/aK0N |
---|
1334 | |
---|
1335 | else if( PDG > 3000 && PDG < 3335) I = 6; // @@ for all hyperons - take Lambda |
---|
1336 | else if( PDG < -2000 && PDG > -3335) I = 7; // @@ for all anti-baryons - anti-p/anti-n |
---|
1337 | else |
---|
1338 | { |
---|
1339 | G4cout<<"MK PDG = "<<PDG |
---|
1340 | <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; |
---|
1341 | G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash"); |
---|
1342 | } |
---|
1343 | |
---|
1344 | // Each parameter set can have not more than nPoints = 128 parameters |
---|
1345 | |
---|
1346 | static const G4double lmi = 3.5; // min of (lnP-lmi)^2 parabola |
---|
1347 | static const G4double pbe = .0557; // elastic (lnP-lmi)^2 parabola coefficient |
---|
1348 | static const G4double pbt = .3; // total (lnP-lmi)^2 parabola coefficient |
---|
1349 | static const G4double pmi = .1; // Below that fast LE calculation is made |
---|
1350 | static const G4double pma = 1000.; // Above that fast HE calculation is made |
---|
1351 | |
---|
1352 | if( p <= 0.) |
---|
1353 | { |
---|
1354 | G4cout<<" p = "<<p<<" is zero or negative"<<G4endl; |
---|
1355 | |
---|
1356 | elasticXsc = 0.; |
---|
1357 | inelasticXsc = 0.; |
---|
1358 | totalXsc = 0.; |
---|
1359 | |
---|
1360 | return totalXsc; |
---|
1361 | } |
---|
1362 | if (!I) // pp/nn |
---|
1363 | { |
---|
1364 | if( p < pmi ) |
---|
1365 | { |
---|
1366 | G4double p2 = p*p; |
---|
1367 | elasticXsc = 1./(.00012 + p2*.2); |
---|
1368 | totalXsc = elasticXsc; |
---|
1369 | } |
---|
1370 | else if(p>pma) |
---|
1371 | { |
---|
1372 | G4double lp = std::log(p)-lmi; |
---|
1373 | G4double lp2 = lp*lp; |
---|
1374 | elasticXsc = pbe*lp2 + 6.72; |
---|
1375 | totalXsc = pbt*lp2 + 38.2; |
---|
1376 | } |
---|
1377 | else |
---|
1378 | { |
---|
1379 | G4double p2 = p*p; |
---|
1380 | G4double LE = 1./( .00012 + p2*.2); |
---|
1381 | G4double lp = std::log(p) - lmi; |
---|
1382 | G4double lp2 = lp*lp; |
---|
1383 | G4double rp2 = 1./p2; |
---|
1384 | elasticXsc = LE + ( pbe*lp2 + 6.72+32.6/p)/( 1. + rp2/p); |
---|
1385 | totalXsc = LE + ( pbt*lp2 + 38.2+52.7*rp2)/( 1. + 2.72*rp2*rp2); |
---|
1386 | } |
---|
1387 | } |
---|
1388 | else if( I==1 ) // np/pn |
---|
1389 | { |
---|
1390 | if( p < pmi ) |
---|
1391 | { |
---|
1392 | G4double p2 = p*p; |
---|
1393 | elasticXsc = 1./( .00012 + p2*( .051 + .1*p2)); |
---|
1394 | totalXsc = elasticXsc; |
---|
1395 | } |
---|
1396 | else if( p > pma ) |
---|
1397 | { |
---|
1398 | G4double lp = std::log(p) - lmi; |
---|
1399 | G4double lp2 = lp*lp; |
---|
1400 | elasticXsc = pbe*lp2 + 6.72; |
---|
1401 | totalXsc = pbt*lp2 + 38.2; |
---|
1402 | } |
---|
1403 | else |
---|
1404 | { |
---|
1405 | G4double p2 = p*p; |
---|
1406 | G4double LE = 1./( .00012 + p2*( .051 + .1*p2 ) ); |
---|
1407 | G4double lp = std::log(p) - lmi; |
---|
1408 | G4double lp2 = lp*lp; |
---|
1409 | G4double rp2 = 1./p2; |
---|
1410 | elasticXsc = LE + (pbe*lp2 + 6.72 + 30./p)/( 1. + .49*rp2/p); |
---|
1411 | totalXsc = LE + (pbt*lp2 + 38.2)/( 1. + .54*rp2*rp2); |
---|
1412 | } |
---|
1413 | } |
---|
1414 | else if( I == 2 ) // pimp/pipn |
---|
1415 | { |
---|
1416 | G4double lp = std::log(p); |
---|
1417 | |
---|
1418 | if(p<pmi) |
---|
1419 | { |
---|
1420 | G4double lr = lp + 1.27; |
---|
1421 | elasticXsc = 1.53/( lr*lr + .0676); |
---|
1422 | totalXsc = elasticXsc*3; |
---|
1423 | } |
---|
1424 | else if( p > pma ) |
---|
1425 | { |
---|
1426 | G4double ld = lp - lmi; |
---|
1427 | G4double ld2 = ld*ld; |
---|
1428 | G4double sp = std::sqrt(p); |
---|
1429 | elasticXsc = pbe*ld2 + 2.4 + 7./sp; |
---|
1430 | totalXsc = pbt*ld2 + 22.3 + 12./sp; |
---|
1431 | } |
---|
1432 | else |
---|
1433 | { |
---|
1434 | G4double lr = lp + 1.27; |
---|
1435 | G4double LE = 1.53/( lr*lr + .0676); |
---|
1436 | G4double ld = lp - lmi; |
---|
1437 | G4double ld2 = ld*ld; |
---|
1438 | G4double p2 = p*p; |
---|
1439 | G4double p4 = p2*p2; |
---|
1440 | G4double sp = std::sqrt(p); |
---|
1441 | G4double lm = lp + .36; |
---|
1442 | G4double md = lm*lm + .04; |
---|
1443 | G4double lh = lp - .017; |
---|
1444 | G4double hd = lh*lh + .0025; |
---|
1445 | elasticXsc = LE + (pbe*ld2 + 2.4 + 7./sp)/( 1. + .7/p4) + .6/md + .05/hd; |
---|
1446 | totalXsc = LE*3 + (pbt*ld2 + 22.3 + 12./sp)/(1. + .4/p4) + 1./md + .06/hd; |
---|
1447 | } |
---|
1448 | } |
---|
1449 | else if( I == 3 ) // pipp/pimn |
---|
1450 | { |
---|
1451 | G4double lp = std::log(p); |
---|
1452 | |
---|
1453 | if( p < pmi ) |
---|
1454 | { |
---|
1455 | G4double lr = lp + 1.27; |
---|
1456 | G4double lr2 = lr*lr; |
---|
1457 | elasticXsc = 13./( lr2 + lr2*lr2 + .0676); |
---|
1458 | totalXsc = elasticXsc; |
---|
1459 | } |
---|
1460 | else if( p > pma ) |
---|
1461 | { |
---|
1462 | G4double ld = lp - lmi; |
---|
1463 | G4double ld2 = ld*ld; |
---|
1464 | G4double sp = std::sqrt(p); |
---|
1465 | elasticXsc = pbe*ld2 + 2.4 + 6./sp; |
---|
1466 | totalXsc = pbt*ld2 + 22.3 + 5./sp; |
---|
1467 | } |
---|
1468 | else |
---|
1469 | { |
---|
1470 | G4double lr = lp + 1.27; |
---|
1471 | G4double lr2 = lr*lr; |
---|
1472 | G4double LE = 13./(lr2 + lr2*lr2 + .0676); |
---|
1473 | G4double ld = lp - lmi; |
---|
1474 | G4double ld2 = ld*ld; |
---|
1475 | G4double p2 = p*p; |
---|
1476 | G4double p4 = p2*p2; |
---|
1477 | G4double sp = std::sqrt(p); |
---|
1478 | G4double lm = lp - .32; |
---|
1479 | G4double md = lm*lm + .0576; |
---|
1480 | elasticXsc = LE + (pbe*ld2 + 2.4 + 6./sp)/(1. + 3./p4) + .7/md; |
---|
1481 | totalXsc = LE + (pbt*ld2 + 22.3 + 5./sp)/(1. + 1./p4) + .8/md; |
---|
1482 | } |
---|
1483 | } |
---|
1484 | else if( I == 4 ) // Kmp/Kmn/K0p/K0n |
---|
1485 | { |
---|
1486 | if( p < pmi) |
---|
1487 | { |
---|
1488 | G4double psp = p*std::sqrt(p); |
---|
1489 | elasticXsc = 5.2/psp; |
---|
1490 | totalXsc = 14./psp; |
---|
1491 | } |
---|
1492 | else if( p > pma ) |
---|
1493 | { |
---|
1494 | G4double ld = std::log(p) - lmi; |
---|
1495 | G4double ld2 = ld*ld; |
---|
1496 | elasticXsc = pbe*ld2 + 2.23; |
---|
1497 | totalXsc = pbt*ld2 + 19.5; |
---|
1498 | } |
---|
1499 | else |
---|
1500 | { |
---|
1501 | G4double ld = std::log(p) - lmi; |
---|
1502 | G4double ld2 = ld*ld; |
---|
1503 | G4double sp = std::sqrt(p); |
---|
1504 | G4double psp = p*sp; |
---|
1505 | G4double p2 = p*p; |
---|
1506 | G4double p4 = p2*p2; |
---|
1507 | G4double lm = p - .39; |
---|
1508 | G4double md = lm*lm + .000156; |
---|
1509 | G4double lh = p - 1.; |
---|
1510 | G4double hd = lh*lh + .0156; |
---|
1511 | elasticXsc = 5.2/psp + (pbe*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md + .15/hd; |
---|
1512 | totalXsc = 14./psp + (pbt*ld2 + 19.5)/(1. - .21/sp + .52/p4) + .006/md + .30/hd; |
---|
1513 | } |
---|
1514 | } |
---|
1515 | else if( I == 5 ) // Kpp/Kpn/aKp/aKn |
---|
1516 | { |
---|
1517 | if( p < pmi ) |
---|
1518 | { |
---|
1519 | G4double lr = p - .38; |
---|
1520 | G4double lm = p - 1.; |
---|
1521 | G4double md = lm*lm + .372; |
---|
1522 | elasticXsc = .7/(lr*lr + .0676) + 2./md; |
---|
1523 | totalXsc = elasticXsc + .6/md; |
---|
1524 | } |
---|
1525 | else if( p > pma ) |
---|
1526 | { |
---|
1527 | G4double ld = std::log(p) - lmi; |
---|
1528 | G4double ld2 = ld*ld; |
---|
1529 | elasticXsc = pbe*ld2 + 2.23; |
---|
1530 | totalXsc = pbt*ld2 + 19.5; |
---|
1531 | } |
---|
1532 | else |
---|
1533 | { |
---|
1534 | G4double ld = std::log(p) - lmi; |
---|
1535 | G4double ld2 = ld*ld; |
---|
1536 | G4double lr = p - .38; |
---|
1537 | G4double LE = .7/(lr*lr + .0676); |
---|
1538 | G4double sp = std::sqrt(p); |
---|
1539 | G4double p2 = p*p; |
---|
1540 | G4double p4 = p2*p2; |
---|
1541 | G4double lm = p - 1.; |
---|
1542 | G4double md = lm*lm + .372; |
---|
1543 | elasticXsc = LE + (pbe*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md; |
---|
1544 | totalXsc = LE + (pbt*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md; |
---|
1545 | } |
---|
1546 | } |
---|
1547 | else if( I == 6 ) // hyperon-N |
---|
1548 | { |
---|
1549 | if( p < pmi ) |
---|
1550 | { |
---|
1551 | G4double p2 = p*p; |
---|
1552 | elasticXsc = 1./(.002 + p2*(.12 + p2)); |
---|
1553 | totalXsc = elasticXsc; |
---|
1554 | } |
---|
1555 | else if( p > pma ) |
---|
1556 | { |
---|
1557 | G4double lp = std::log(p) - lmi; |
---|
1558 | G4double lp2 = lp*lp; |
---|
1559 | G4double sp = std::sqrt(p); |
---|
1560 | elasticXsc = (pbe*lp2 + 6.72)/(1. + 2./sp); |
---|
1561 | totalXsc = (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp); |
---|
1562 | } |
---|
1563 | else |
---|
1564 | { |
---|
1565 | G4double p2 = p*p; |
---|
1566 | G4double LE = 1./(.002 + p2*(.12 + p2)); |
---|
1567 | G4double lp = std::log(p) - lmi; |
---|
1568 | G4double lp2 = lp*lp; |
---|
1569 | G4double p4 = p2*p2; |
---|
1570 | G4double sp = std::sqrt(p); |
---|
1571 | elasticXsc = LE + (pbe*lp2 + 6.72 + 99./p2)/(1. + 2./sp + 2./p4); |
---|
1572 | totalXsc = LE + (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp + 3./p4); |
---|
1573 | } |
---|
1574 | } |
---|
1575 | else if( I == 7 ) // antibaryon-N |
---|
1576 | { |
---|
1577 | if( p > pma ) |
---|
1578 | { |
---|
1579 | G4double lp = std::log(p) - lmi; |
---|
1580 | G4double lp2 = lp*lp; |
---|
1581 | elasticXsc = pbe*lp2 + 6.72; |
---|
1582 | totalXsc = pbt*lp2 + 38.2; |
---|
1583 | } |
---|
1584 | else |
---|
1585 | { |
---|
1586 | G4double ye = std::pow(p, 1.25); |
---|
1587 | G4double yt = std::pow(p, .35); |
---|
1588 | G4double lp = std::log(p) - lmi; |
---|
1589 | G4double lp2 = lp*lp; |
---|
1590 | elasticXsc = 80./(ye + 1.) + pbe*lp2 + 6.72; |
---|
1591 | totalXsc = (80./yt + .3)/yt +pbt*lp2 + 38.2; |
---|
1592 | } |
---|
1593 | } |
---|
1594 | else |
---|
1595 | { |
---|
1596 | G4cout<<"PDG incoding = "<<I<<" is not defined (0-7)"<<G4endl; |
---|
1597 | |
---|
1598 | } |
---|
1599 | if( elasticXsc > totalXsc ) elasticXsc = totalXsc; |
---|
1600 | |
---|
1601 | totalXsc *= millibarn; |
---|
1602 | elasticXsc *= millibarn; |
---|
1603 | inelasticXsc = totalXsc - elasticXsc; |
---|
1604 | if (inelasticXsc < 0.) inelasticXsc = 0.; |
---|
1605 | |
---|
1606 | return inelasticXsc; |
---|
1607 | } |
---|
1608 | |
---|
1609 | //////////////////////////////////////////////////////////////////////////////////// |
---|
1610 | // |
---|
1611 | // |
---|
1612 | |
---|
1613 | G4double |
---|
1614 | G4GlauberGribovCrossSection::GetNucleusRadius( const G4DynamicParticle* , |
---|
1615 | const G4Element* anElement) |
---|
1616 | { |
---|
1617 | G4double At = anElement->GetN(); |
---|
1618 | G4double oneThird = 1.0/3.0; |
---|
1619 | G4double cubicrAt = std::pow (At, oneThird); |
---|
1620 | |
---|
1621 | |
---|
1622 | G4double R; // = fRadiusConst*cubicrAt; |
---|
1623 | /* |
---|
1624 | G4double tmp = std::pow( cubicrAt-1., 3.); |
---|
1625 | tmp += At; |
---|
1626 | tmp *= 0.5; |
---|
1627 | |
---|
1628 | if (At > 20.) // 20. |
---|
1629 | { |
---|
1630 | R = fRadiusConst*std::pow (tmp, oneThird); |
---|
1631 | } |
---|
1632 | else |
---|
1633 | { |
---|
1634 | R = fRadiusConst*cubicrAt; |
---|
1635 | } |
---|
1636 | */ |
---|
1637 | |
---|
1638 | R = fRadiusConst*cubicrAt; |
---|
1639 | |
---|
1640 | // return R; // !!!! |
---|
1641 | |
---|
1642 | |
---|
1643 | |
---|
1644 | G4double meanA = 21.; |
---|
1645 | |
---|
1646 | G4double tauA1 = 40.; |
---|
1647 | G4double tauA2 = 10.; |
---|
1648 | G4double tauA3 = 5.; |
---|
1649 | |
---|
1650 | G4double a1 = 0.85; |
---|
1651 | G4double b1 = 1. - a1; |
---|
1652 | |
---|
1653 | G4double b2 = 0.3; |
---|
1654 | G4double b3 = 4.; |
---|
1655 | |
---|
1656 | if (At > 20.) // 20. |
---|
1657 | { |
---|
1658 | R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) ); |
---|
1659 | } |
---|
1660 | else if (At > 3.5) |
---|
1661 | { |
---|
1662 | R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) ); |
---|
1663 | } |
---|
1664 | else |
---|
1665 | { |
---|
1666 | R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) ); |
---|
1667 | } |
---|
1668 | return R; |
---|
1669 | |
---|
1670 | } |
---|
1671 | //////////////////////////////////////////////////////////////////////////////////// |
---|
1672 | // |
---|
1673 | // |
---|
1674 | |
---|
1675 | G4double |
---|
1676 | G4GlauberGribovCrossSection::GetNucleusRadius(G4double At) |
---|
1677 | { |
---|
1678 | G4double oneThird = 1.0/3.0; |
---|
1679 | G4double cubicrAt = std::pow (At, oneThird); |
---|
1680 | |
---|
1681 | |
---|
1682 | G4double R; // = fRadiusConst*cubicrAt; |
---|
1683 | |
---|
1684 | /* |
---|
1685 | G4double tmp = std::pow( cubicrAt-1., 3.); |
---|
1686 | tmp += At; |
---|
1687 | tmp *= 0.5; |
---|
1688 | |
---|
1689 | if (At > 20.) |
---|
1690 | { |
---|
1691 | R = fRadiusConst*std::pow (tmp, oneThird); |
---|
1692 | } |
---|
1693 | else |
---|
1694 | { |
---|
1695 | R = fRadiusConst*cubicrAt; |
---|
1696 | } |
---|
1697 | */ |
---|
1698 | |
---|
1699 | R = fRadiusConst*cubicrAt; |
---|
1700 | |
---|
1701 | G4double meanA = 20.; |
---|
1702 | G4double tauA = 20.; |
---|
1703 | |
---|
1704 | if (At > 20.) // 20. |
---|
1705 | { |
---|
1706 | R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) ); |
---|
1707 | } |
---|
1708 | else |
---|
1709 | { |
---|
1710 | R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) ); |
---|
1711 | } |
---|
1712 | |
---|
1713 | return R; |
---|
1714 | } |
---|
1715 | |
---|
1716 | //////////////////////////////////////////////////////////////////////////////////// |
---|
1717 | // |
---|
1718 | // |
---|
1719 | |
---|
1720 | G4double G4GlauberGribovCrossSection::CalculateEcmValue( const G4double mp , |
---|
1721 | const G4double mt , |
---|
1722 | const G4double Plab ) |
---|
1723 | { |
---|
1724 | G4double Elab = std::sqrt ( mp * mp + Plab * Plab ); |
---|
1725 | G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt ); |
---|
1726 | // G4double Pcm = Plab * mt / Ecm; |
---|
1727 | // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp; |
---|
1728 | |
---|
1729 | return Ecm ; // KEcm; |
---|
1730 | } |
---|
1731 | |
---|
1732 | |
---|
1733 | //////////////////////////////////////////////////////////////////////////////////// |
---|
1734 | // |
---|
1735 | // |
---|
1736 | |
---|
1737 | G4double G4GlauberGribovCrossSection::CalcMandelstamS( const G4double mp , |
---|
1738 | const G4double mt , |
---|
1739 | const G4double Plab ) |
---|
1740 | { |
---|
1741 | G4double Elab = std::sqrt ( mp * mp + Plab * Plab ); |
---|
1742 | G4double sMand = mp*mp + mt*mt + 2*Elab*mt ; |
---|
1743 | |
---|
1744 | return sMand; |
---|
1745 | } |
---|
1746 | |
---|
1747 | |
---|
1748 | // |
---|
1749 | // |
---|
1750 | /////////////////////////////////////////////////////////////////////////////////////// |
---|