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