00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "G4GlauberGribovCrossSection.hh"
00035
00036 #include "G4PhysicalConstants.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4ParticleTable.hh"
00039 #include "G4IonTable.hh"
00040 #include "G4ParticleDefinition.hh"
00041 #include "G4HadronNucleonXsc.hh"
00042
00043
00044 #include "G4CrossSectionFactory.hh"
00045
00046 G4_DECLARE_XS_FACTORY(G4GlauberGribovCrossSection);
00047
00049
00050 const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionTot[93] = {
00051
00052 1.0, 1.0, 1.118517e+00, 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00,
00053 1.058205e+00, 1.082663e+00, 1.068500e+00, 1.076912e+00, 1.083475e+00, 1.079117e+00,
00054 1.071856e+00, 1.071990e+00, 1.073774e+00, 1.079356e+00, 1.081314e+00, 1.082056e+00,
00055 1.090772e+00, 1.096776e+00, 1.095828e+00, 1.097678e+00, 1.099157e+00, 1.103677e+00,
00056 1.105132e+00, 1.109806e+00, 1.110816e+00, 1.117378e+00, 1.115165e+00, 1.115710e+00,
00057 1.111855e+00, 1.110482e+00, 1.110112e+00, 1.106676e+00, 1.108706e+00, 1.105549e+00,
00058 1.106318e+00, 1.106242e+00, 1.107672e+00, 1.107342e+00, 1.108119e+00, 1.106655e+00,
00059 1.102588e+00, 1.096657e+00, 1.092920e+00, 1.086629e+00, 1.083592e+00, 1.076030e+00,
00060 1.083777e+00, 1.089460e+00, 1.086545e+00, 1.079924e+00, 1.082218e+00, 1.077798e+00,
00061 1.077062e+00, 1.072825e+00, 1.072241e+00, 1.072104e+00, 1.072490e+00, 1.069829e+00,
00062 1.070398e+00, 1.065458e+00, 1.064968e+00, 1.060524e+00, 1.060048e+00, 1.057620e+00,
00063 1.056428e+00, 1.055366e+00, 1.055017e+00, 1.052304e+00, 1.051767e+00, 1.049728e+00,
00064 1.048745e+00, 1.047399e+00, 1.045876e+00, 1.042972e+00, 1.041824e+00, 1.039993e+00,
00065 1.039021e+00, 1.036627e+00, 1.034176e+00, 1.032526e+00, 1.033633e+00, 1.036107e+00,
00066 1.037803e+00, 1.031266e+00, 1.032991e+00, 1.033284e+00, 1.035015e+00, 1.033945e+00,
00067 1.037075e+00, 1.034721e+00
00068
00069 };
00070
00071 const G4double G4GlauberGribovCrossSection::fNeutronBarCorrectionIn[93] = {
00072
00073 1.0, 1.0, 1.167421e+00, 1.156250e+00, 1.205364e+00, 1.154225e+00, 1.120391e+00,
00074 1.124632e+00, 1.129460e+00, 1.107863e+00, 1.102152e+00, 1.104593e+00, 1.100285e+00,
00075 1.098450e+00, 1.092677e+00, 1.101124e+00, 1.106461e+00, 1.115049e+00, 1.123903e+00,
00076 1.126661e+00, 1.131259e+00, 1.133949e+00, 1.134185e+00, 1.133767e+00, 1.132813e+00,
00077 1.131515e+00, 1.130338e+00, 1.134171e+00, 1.139206e+00, 1.141474e+00, 1.142189e+00,
00078 1.140725e+00, 1.140100e+00, 1.139848e+00, 1.137674e+00, 1.138645e+00, 1.136339e+00,
00079 1.136439e+00, 1.135946e+00, 1.136431e+00, 1.135702e+00, 1.135703e+00, 1.134113e+00,
00080 1.131935e+00, 1.128381e+00, 1.126373e+00, 1.122453e+00, 1.120908e+00, 1.115953e+00,
00081 1.115947e+00, 1.114426e+00, 1.111749e+00, 1.106207e+00, 1.107494e+00, 1.103622e+00,
00082 1.102576e+00, 1.098816e+00, 1.097889e+00, 1.097306e+00, 1.097130e+00, 1.094578e+00,
00083 1.094552e+00, 1.090222e+00, 1.089358e+00, 1.085409e+00, 1.084560e+00, 1.082182e+00,
00084 1.080773e+00, 1.079464e+00, 1.078724e+00, 1.076121e+00, 1.075235e+00, 1.073159e+00,
00085 1.071920e+00, 1.070395e+00, 1.069503e+00, 1.067525e+00, 1.066919e+00, 1.065779e+00,
00086 1.065319e+00, 1.063730e+00, 1.062092e+00, 1.061085e+00, 1.059908e+00, 1.059815e+00,
00087 1.059109e+00, 1.051920e+00, 1.051258e+00, 1.049473e+00, 1.048823e+00, 1.045984e+00,
00088 1.046435e+00, 1.042614e+00
00089
00090 };
00091
00092 const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionTot[93] = {
00093
00094 1.0, 1.0,
00095 1.118515e+00, 1.082000e+00, 1.116169e+00, 1.078745e+00, 1.061313e+00, 1.058203e+00,
00096 1.082661e+00, 1.068498e+00, 1.076910e+00, 1.083474e+00, 1.079115e+00, 1.071854e+00,
00097 1.071988e+00, 1.073772e+00, 1.079355e+00, 1.081312e+00, 1.082054e+00, 1.090770e+00,
00098 1.096774e+00, 1.095827e+00, 1.097677e+00, 1.099156e+00, 1.103676e+00, 1.105130e+00,
00099 1.109805e+00, 1.110814e+00, 1.117377e+00, 1.115163e+00, 1.115708e+00, 1.111853e+00,
00100 1.110480e+00, 1.110111e+00, 1.106674e+00, 1.108705e+00, 1.105548e+00, 1.106317e+00,
00101 1.106241e+00, 1.107671e+00, 1.107341e+00, 1.108118e+00, 1.106654e+00, 1.102586e+00,
00102 1.096655e+00, 1.092918e+00, 1.086628e+00, 1.083590e+00, 1.076028e+00, 1.083776e+00,
00103 1.089458e+00, 1.086543e+00, 1.079923e+00, 1.082216e+00, 1.077797e+00, 1.077061e+00,
00104 1.072824e+00, 1.072239e+00, 1.072103e+00, 1.072488e+00, 1.069828e+00, 1.070396e+00,
00105 1.065456e+00, 1.064966e+00, 1.060523e+00, 1.060047e+00, 1.057618e+00, 1.056427e+00,
00106 1.055365e+00, 1.055016e+00, 1.052303e+00, 1.051766e+00, 1.049727e+00, 1.048743e+00,
00107 1.047397e+00, 1.045875e+00, 1.042971e+00, 1.041823e+00, 1.039992e+00, 1.039019e+00,
00108 1.036626e+00, 1.034175e+00, 1.032525e+00, 1.033632e+00, 1.036106e+00, 1.037802e+00,
00109 1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00,
00110 1.034720e+00
00111
00112 };
00113
00114 const G4double G4GlauberGribovCrossSection::fProtonBarCorrectionIn[93] = {
00115
00116 1.0, 1.0,
00117 1.167419e+00, 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.120390e+00, 1.124630e+00,
00118 1.129459e+00, 1.107861e+00, 1.102151e+00, 1.104591e+00, 1.100284e+00, 1.098449e+00,
00119 1.092675e+00, 1.101122e+00, 1.106460e+00, 1.115048e+00, 1.123902e+00, 1.126659e+00,
00120 1.131258e+00, 1.133948e+00, 1.134183e+00, 1.133766e+00, 1.132812e+00, 1.131514e+00,
00121 1.130337e+00, 1.134170e+00, 1.139205e+00, 1.141472e+00, 1.142188e+00, 1.140724e+00,
00122 1.140099e+00, 1.139847e+00, 1.137672e+00, 1.138644e+00, 1.136338e+00, 1.136438e+00,
00123 1.135945e+00, 1.136429e+00, 1.135701e+00, 1.135702e+00, 1.134112e+00, 1.131934e+00,
00124 1.128380e+00, 1.126371e+00, 1.122452e+00, 1.120907e+00, 1.115952e+00, 1.115946e+00,
00125 1.114425e+00, 1.111748e+00, 1.106205e+00, 1.107493e+00, 1.103621e+00, 1.102575e+00,
00126 1.098815e+00, 1.097888e+00, 1.097305e+00, 1.097129e+00, 1.094577e+00, 1.094551e+00,
00127 1.090221e+00, 1.089357e+00, 1.085408e+00, 1.084559e+00, 1.082181e+00, 1.080772e+00,
00128 1.079463e+00, 1.078723e+00, 1.076120e+00, 1.075234e+00, 1.073158e+00, 1.071919e+00,
00129 1.070394e+00, 1.069502e+00, 1.067524e+00, 1.066918e+00, 1.065778e+00, 1.065318e+00,
00130 1.063729e+00, 1.062091e+00, 1.061084e+00, 1.059907e+00, 1.059814e+00, 1.059108e+00,
00131 1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00,
00132 1.042613e+00
00133
00134 };
00135
00136
00137 const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionTot[93] = {
00138
00139 1.0, 1.0,
00140 1.075927e+00, 1.074407e+00, 1.126098e+00, 1.100127e+00, 1.089742e+00, 1.083536e+00,
00141 1.089988e+00, 1.103566e+00, 1.096922e+00, 1.126573e+00, 1.132734e+00, 1.136512e+00,
00142 1.136629e+00, 1.133086e+00, 1.132428e+00, 1.129299e+00, 1.125622e+00, 1.126992e+00,
00143 1.127840e+00, 1.162670e+00, 1.160392e+00, 1.157864e+00, 1.157227e+00, 1.154627e+00,
00144 1.192555e+00, 1.197243e+00, 1.197911e+00, 1.200326e+00, 1.220053e+00, 1.215019e+00,
00145 1.211703e+00, 1.209080e+00, 1.204248e+00, 1.203328e+00, 1.198671e+00, 1.196840e+00,
00146 1.194392e+00, 1.193037e+00, 1.190408e+00, 1.188583e+00, 1.206127e+00, 1.210028e+00,
00147 1.206434e+00, 1.204456e+00, 1.200547e+00, 1.199058e+00, 1.200174e+00, 1.200276e+00,
00148 1.198912e+00, 1.213048e+00, 1.207160e+00, 1.208020e+00, 1.203814e+00, 1.202380e+00,
00149 1.198306e+00, 1.197002e+00, 1.196027e+00, 1.195449e+00, 1.192563e+00, 1.192135e+00,
00150 1.187556e+00, 1.186308e+00, 1.182124e+00, 1.180900e+00, 1.178224e+00, 1.176471e+00,
00151 1.174811e+00, 1.173702e+00, 1.170827e+00, 1.169581e+00, 1.167205e+00, 1.165626e+00,
00152 1.180244e+00, 1.177626e+00, 1.175121e+00, 1.173903e+00, 1.172192e+00, 1.171128e+00,
00153 1.168997e+00, 1.166826e+00, 1.164130e+00, 1.165412e+00, 1.165504e+00, 1.165020e+00,
00154 1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00,
00155 1.152974e+00
00156
00157 };
00158
00159 const G4double G4GlauberGribovCrossSection::fPionPlusBarCorrectionIn[93] = {
00160
00161 1.0, 1.0,
00162 1.140246e+00, 1.097872e+00, 1.104301e+00, 1.068722e+00, 1.044495e+00, 1.062622e+00,
00163 1.047987e+00, 1.037032e+00, 1.035686e+00, 1.042870e+00, 1.052222e+00, 1.065100e+00,
00164 1.070480e+00, 1.078286e+00, 1.081488e+00, 1.089713e+00, 1.099105e+00, 1.098003e+00,
00165 1.102175e+00, 1.117707e+00, 1.121734e+00, 1.125229e+00, 1.126457e+00, 1.128905e+00,
00166 1.137312e+00, 1.126263e+00, 1.126459e+00, 1.115191e+00, 1.116986e+00, 1.117184e+00,
00167 1.117037e+00, 1.116777e+00, 1.115858e+00, 1.115745e+00, 1.114489e+00, 1.113993e+00,
00168 1.113226e+00, 1.112818e+00, 1.111890e+00, 1.111238e+00, 1.111209e+00, 1.111775e+00,
00169 1.110256e+00, 1.109414e+00, 1.107647e+00, 1.106980e+00, 1.106096e+00, 1.107331e+00,
00170 1.107849e+00, 1.106407e+00, 1.103426e+00, 1.103896e+00, 1.101756e+00, 1.101031e+00,
00171 1.098915e+00, 1.098260e+00, 1.097768e+00, 1.097487e+00, 1.095964e+00, 1.095773e+00,
00172 1.093348e+00, 1.092687e+00, 1.090465e+00, 1.089821e+00, 1.088394e+00, 1.087462e+00,
00173 1.086571e+00, 1.085997e+00, 1.084451e+00, 1.083798e+00, 1.082513e+00, 1.081670e+00,
00174 1.080735e+00, 1.075659e+00, 1.074341e+00, 1.073689e+00, 1.072787e+00, 1.072237e+00,
00175 1.071107e+00, 1.069955e+00, 1.064856e+00, 1.065873e+00, 1.065938e+00, 1.065694e+00,
00176 1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00,
00177 1.059658e+00
00178
00179 };
00180
00181
00182 const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionTot[93] = {
00183
00184 1.0, 1.0,
00185 1.075927e+00, 1.077959e+00, 1.129145e+00, 1.102088e+00, 1.089765e+00, 1.083542e+00,
00186 1.089995e+00, 1.104895e+00, 1.097154e+00, 1.127663e+00, 1.133063e+00, 1.137425e+00,
00187 1.136724e+00, 1.133859e+00, 1.132498e+00, 1.130276e+00, 1.127896e+00, 1.127656e+00,
00188 1.127905e+00, 1.164210e+00, 1.162259e+00, 1.160075e+00, 1.158978e+00, 1.156649e+00,
00189 1.194157e+00, 1.199177e+00, 1.198983e+00, 1.202325e+00, 1.221967e+00, 1.217548e+00,
00190 1.214389e+00, 1.211760e+00, 1.207335e+00, 1.206081e+00, 1.201766e+00, 1.199779e+00,
00191 1.197283e+00, 1.195706e+00, 1.193071e+00, 1.191115e+00, 1.208838e+00, 1.212681e+00,
00192 1.209235e+00, 1.207163e+00, 1.203451e+00, 1.201807e+00, 1.203283e+00, 1.203388e+00,
00193 1.202244e+00, 1.216509e+00, 1.211066e+00, 1.211504e+00, 1.207539e+00, 1.205991e+00,
00194 1.202143e+00, 1.200724e+00, 1.199595e+00, 1.198815e+00, 1.196025e+00, 1.195390e+00,
00195 1.191137e+00, 1.189791e+00, 1.185888e+00, 1.184575e+00, 1.181996e+00, 1.180229e+00,
00196 1.178545e+00, 1.177355e+00, 1.174616e+00, 1.173312e+00, 1.171016e+00, 1.169424e+00,
00197 1.184120e+00, 1.181478e+00, 1.179085e+00, 1.177817e+00, 1.176124e+00, 1.175003e+00,
00198 1.172947e+00, 1.170858e+00, 1.168170e+00, 1.169397e+00, 1.169304e+00, 1.168706e+00,
00199 1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00,
00200 1.157267e+00
00201 };
00202
00203
00204 const G4double G4GlauberGribovCrossSection::fPionMinusBarCorrectionIn[93] = {
00205
00206 1.0, 1.0,
00207 1.140246e+00, 1.100898e+00, 1.106773e+00, 1.070289e+00, 1.044514e+00, 1.062628e+00,
00208 1.047992e+00, 1.038041e+00, 1.035862e+00, 1.043679e+00, 1.052466e+00, 1.065780e+00,
00209 1.070551e+00, 1.078869e+00, 1.081541e+00, 1.090455e+00, 1.100847e+00, 1.098511e+00,
00210 1.102226e+00, 1.118865e+00, 1.123143e+00, 1.126904e+00, 1.127785e+00, 1.130444e+00,
00211 1.138502e+00, 1.127678e+00, 1.127244e+00, 1.116634e+00, 1.118347e+00, 1.118988e+00,
00212 1.118957e+00, 1.118696e+00, 1.118074e+00, 1.117722e+00, 1.116717e+00, 1.116111e+00,
00213 1.115311e+00, 1.114745e+00, 1.113814e+00, 1.113069e+00, 1.113141e+00, 1.113660e+00,
00214 1.112249e+00, 1.111343e+00, 1.109718e+00, 1.108942e+00, 1.108310e+00, 1.109549e+00,
00215 1.110227e+00, 1.108846e+00, 1.106183e+00, 1.106354e+00, 1.104388e+00, 1.103583e+00,
00216 1.101632e+00, 1.100896e+00, 1.100296e+00, 1.099873e+00, 1.098420e+00, 1.098082e+00,
00217 1.095892e+00, 1.095162e+00, 1.093144e+00, 1.092438e+00, 1.091083e+00, 1.090142e+00,
00218 1.089236e+00, 1.088604e+00, 1.087159e+00, 1.086465e+00, 1.085239e+00, 1.084388e+00,
00219 1.083473e+00, 1.078373e+00, 1.077136e+00, 1.076450e+00, 1.075561e+00, 1.074973e+00,
00220 1.073898e+00, 1.072806e+00, 1.067706e+00, 1.068684e+00, 1.068618e+00, 1.068294e+00,
00221 1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00,
00222 1.062699e+00
00223
00224 };
00225
00226
00228
00229
00230 G4GlauberGribovCrossSection::G4GlauberGribovCrossSection()
00231 : G4VCrossSectionDataSet(Default_Name()),
00232 fUpperLimit(100000*GeV), fLowerLimit(10.*MeV),
00233 fRadiusConst(1.08*fermi),
00234 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
00235 fDiffractionXsc(0.0), fHadronNucleonXsc(0.0)
00236 {
00237 theGamma = G4Gamma::Gamma();
00238 theProton = G4Proton::Proton();
00239 theNeutron = G4Neutron::Neutron();
00240 theAProton = G4AntiProton::AntiProton();
00241 theANeutron = G4AntiNeutron::AntiNeutron();
00242 thePiPlus = G4PionPlus::PionPlus();
00243 thePiMinus = G4PionMinus::PionMinus();
00244 thePiZero = G4PionZero::PionZero();
00245 theKPlus = G4KaonPlus::KaonPlus();
00246 theKMinus = G4KaonMinus::KaonMinus();
00247 theK0S = G4KaonZeroShort::KaonZeroShort();
00248 theK0L = G4KaonZeroLong::KaonZeroLong();
00249 theL = G4Lambda::Lambda();
00250 theAntiL = G4AntiLambda::AntiLambda();
00251 theSPlus = G4SigmaPlus::SigmaPlus();
00252 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
00253 theSMinus = G4SigmaMinus::SigmaMinus();
00254 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
00255 theS0 = G4SigmaZero::SigmaZero();
00256 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
00257 theXiMinus = G4XiMinus::XiMinus();
00258 theXi0 = G4XiZero::XiZero();
00259 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00260 theAXi0 = G4AntiXiZero::AntiXiZero();
00261 theOmega = G4OmegaMinus::OmegaMinus();
00262 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
00263 theD = G4Deuteron::Deuteron();
00264 theT = G4Triton::Triton();
00265 theA = G4Alpha::Alpha();
00266 theHe3 = G4He3::He3();
00267
00268 hnXsc = new G4HadronNucleonXsc();
00269 }
00270
00272
00273
00274
00275 G4GlauberGribovCrossSection::~G4GlauberGribovCrossSection()
00276 {
00277 if (hnXsc) delete hnXsc;
00278 }
00279
00281
00282 G4bool
00283 G4GlauberGribovCrossSection::IsIsoApplicable(const G4DynamicParticle* aDP,
00284 G4int Z, G4int ,
00285 const G4Element*,
00286 const G4Material*)
00287 {
00288 G4bool applicable = false;
00289
00290 G4double kineticEnergy = aDP->GetKineticEnergy();
00291
00292 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
00293
00294 if ( ( kineticEnergy >= fLowerLimit &&
00295 Z > 1 &&
00296 ( theParticle == theAProton ||
00297 theParticle == theGamma ||
00298 theParticle == theKPlus ||
00299 theParticle == theKMinus ||
00300 theParticle == theK0L ||
00301 theParticle == theK0S ||
00302 theParticle == theSMinus ||
00303 theParticle == theProton ||
00304 theParticle == theNeutron ||
00305 theParticle == thePiPlus ||
00306 theParticle == thePiMinus ) ) ) applicable = true;
00307
00308 return applicable;
00309 }
00310
00312
00313
00314
00315
00316
00317
00318 G4double
00319 G4GlauberGribovCrossSection::GetIsoCrossSection(const G4DynamicParticle* aParticle,
00320 G4int Z, G4int A,
00321 const G4Isotope*,
00322 const G4Element*,
00323 const G4Material*)
00324 {
00325 G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00326 G4double hpInXsc(0.), hnInXsc(0.);
00327 G4double R = GetNucleusRadius(A);
00328
00329 G4int N = A - Z;
00330 if (N < 0) N = 0;
00331
00332 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00333
00334 if( theParticle == theProton ||
00335 theParticle == theNeutron ||
00336 theParticle == thePiPlus ||
00337 theParticle == thePiMinus )
00338 {
00339
00340
00341 sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
00342
00343 hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00344
00345 sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
00346
00347 hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
00348
00349 cofInelastic = 2.4;
00350 cofTotal = 2.0;
00351 }
00352 else if( theParticle == theKPlus ||
00353 theParticle == theKMinus ||
00354 theParticle == theK0S ||
00355 theParticle == theK0L )
00356 {
00357 sigma = GetKaonNucleonXscVector(aParticle, A, Z);
00358 cofInelastic = 2.2;
00359 cofTotal = 2.0;
00360 R = 1.3*fermi;
00361 R *= std::pow(G4double(A), 0.3333);
00362 }
00363 else
00364 {
00365 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
00366 cofInelastic = 2.2;
00367 cofTotal = 2.0;
00368 }
00369
00370
00371 if( A > 1 )
00372 {
00373 nucleusSquare = cofTotal*pi*R*R;
00374 ratio = sigma/nucleusSquare;
00375
00376 xsection = nucleusSquare*std::log( 1. + ratio );
00377
00378 xsection *= GetParticleBarCorTot(theParticle, Z);
00379
00380 fTotalXsc = xsection;
00381
00382
00383
00384 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00385
00386 fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
00387
00388 fElasticXsc = fTotalXsc - fInelasticXsc;
00389
00390 if(fElasticXsc < 0.) fElasticXsc = 0.;
00391
00392 G4double difratio = ratio/(1.+ratio);
00393
00394 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00395
00396
00397
00398
00399 sigma = Z*hpInXsc + N*hnInXsc;
00400
00401 ratio = sigma/nucleusSquare;
00402
00403 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00404
00405 if (fElasticXsc < 0.) fElasticXsc = 0.;
00406 }
00407 else
00408 {
00409 fTotalXsc = sigma;
00410 xsection = sigma;
00411
00412 if ( theParticle != theAProton )
00413 {
00414 sigma = GetHNinelasticXsc(aParticle, A, Z);
00415 fInelasticXsc = sigma;
00416 fElasticXsc = fTotalXsc - fInelasticXsc;
00417 }
00418 else
00419 {
00420 fElasticXsc = fTotalXsc - fInelasticXsc;
00421 }
00422 if (fElasticXsc < 0.) fElasticXsc = 0.;
00423
00424 }
00425 return xsection;
00426 }
00427
00429
00430
00431
00432 G4double G4GlauberGribovCrossSection::
00433 GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
00434 {
00435 G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00436 G4double R = GetNucleusRadius(A);
00437
00438 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00439
00440 if( theParticle == theProton ||
00441 theParticle == theNeutron ||
00442 theParticle == thePiPlus ||
00443 theParticle == thePiMinus )
00444 {
00445 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
00446 cofInelastic = 2.4;
00447 cofTotal = 2.0;
00448 }
00449 else
00450 {
00451 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
00452 cofInelastic = 2.2;
00453 cofTotal = 2.0;
00454 }
00455 nucleusSquare = cofTotal*pi*R*R;
00456 ratio = sigma/nucleusSquare;
00457
00458 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00459
00460 G4double difratio = ratio/(1.+ratio);
00461
00462 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
00463
00464 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
00465 else ratio = 0.;
00466
00467 return ratio;
00468 }
00469
00471
00472
00473
00474 G4double G4GlauberGribovCrossSection::
00475 GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
00476 {
00477 G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
00478 G4double R = GetNucleusRadius(A);
00479
00480 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00481
00482 if( theParticle == theProton ||
00483 theParticle == theNeutron ||
00484 theParticle == thePiPlus ||
00485 theParticle == thePiMinus )
00486 {
00487 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
00488 cofInelastic = 2.4;
00489 cofTotal = 2.0;
00490 }
00491 else
00492 {
00493 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
00494 cofInelastic = 2.2;
00495 cofTotal = 2.0;
00496 }
00497 nucleusSquare = cofTotal*pi*R*R;
00498 ratio = sigma/nucleusSquare;
00499
00500 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00501
00502 sigma = GetHNinelasticXsc(aParticle, A, Z);
00503 ratio = sigma/nucleusSquare;
00504
00505 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
00506
00507 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
00508 else ratio = 0.;
00509 if ( ratio < 0. ) ratio = 0.;
00510
00511 return ratio;
00512 }
00513
00515
00516
00517
00518
00519
00520
00521 G4double
00522 G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
00523 const G4Element* anElement)
00524 {
00525 G4int At = G4lrint(anElement->GetN());
00526 G4int Zt = G4lrint(anElement->GetZ());
00527
00528 return GetHadronNucleonXsc(aParticle, At, Zt);
00529 }
00530
00532
00533
00534
00535
00536
00537
00538 G4double
00539 G4GlauberGribovCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
00540 G4int At, G4int )
00541 {
00542 G4double xsection;
00543
00544
00545
00546 G4double targ_mass = 0.939*GeV;
00547
00548 G4double proj_mass = aParticle->GetMass();
00549 G4double proj_momentum = aParticle->GetMomentum().mag();
00550 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00551
00552 sMand /= GeV*GeV;
00553 proj_momentum /= GeV;
00554
00555 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00556
00557 G4double aa = At;
00558
00559 if(theParticle == theGamma)
00560 {
00561 xsection = aa*(0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
00562 }
00563 else if(theParticle == theNeutron)
00564 {
00565 xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00566 }
00567 else if(theParticle == theProton)
00568 {
00569 xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00570
00571
00572 }
00573 else if(theParticle == theAProton)
00574 {
00575 xsection = aa*( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
00576 }
00577 else if(theParticle == thePiPlus)
00578 {
00579 xsection = aa*(13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
00580 }
00581 else if(theParticle == thePiMinus)
00582 {
00583
00584 xsection = aa*(13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
00585 }
00586 else if(theParticle == theKPlus)
00587 {
00588 xsection = aa*(11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
00589 }
00590 else if(theParticle == theKMinus)
00591 {
00592 xsection = aa*(11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
00593 }
00594 else
00595 {
00596 xsection = aa*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
00597 }
00598 xsection *= millibarn;
00599 return xsection;
00600 }
00601
00602
00604
00605
00606
00607
00608 G4double
00609 G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
00610 const G4Element* anElement)
00611 {
00612 G4int At = G4lrint(anElement->GetN());
00613 G4int Zt = G4lrint(anElement->GetZ());
00614
00615 return GetHadronNucleonXscPDG(aParticle, At, Zt);
00616 }
00617
00618
00619
00620
00622
00623
00624
00625
00626
00627 G4double
00628 G4GlauberGribovCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
00629 G4int At, G4int Zt)
00630 {
00631 G4double xsection;
00632
00633 G4int Nt = At-Zt;
00634 if (Nt < 0) Nt = 0;
00635
00636 G4double zz = Zt;
00637 G4double aa = At;
00638 G4double nn = Nt;
00639
00640 G4double targ_mass = G4ParticleTable::GetParticleTable()->
00641 GetIonTable()->GetIonMass(Zt, At);
00642
00643 targ_mass = 0.939*GeV;
00644
00645 G4double proj_mass = aParticle->GetMass();
00646 G4double proj_momentum = aParticle->GetMomentum().mag();
00647
00648 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00649
00650 sMand /= GeV*GeV;
00651
00652
00653
00654 G4double s0 = 5.38*5.38;
00655 G4double eta1 = 0.458;
00656 G4double eta2 = 0.458;
00657 G4double B = 0.308;
00658
00659
00660 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00661
00662
00663 if(theParticle == theNeutron)
00664 {
00665 xsection = zz*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00666 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00667 xsection += nn*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00668 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00669 }
00670 else if(theParticle == theProton)
00671 {
00672
00673 xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00674 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00675
00676 xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00677 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00678 }
00679 else if(theParticle == theAProton)
00680 {
00681 xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00682 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
00683
00684 xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00685 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
00686 }
00687 else if(theParticle == thePiPlus)
00688 {
00689 xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
00690 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
00691 }
00692 else if(theParticle == thePiMinus)
00693 {
00694 xsection = aa*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
00695 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
00696 }
00697 else if(theParticle == theKPlus || theParticle == theK0L )
00698 {
00699 xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
00700 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
00701
00702 xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
00703 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
00704 }
00705 else if(theParticle == theKMinus || theParticle == theK0S )
00706 {
00707 xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
00708 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
00709
00710 xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
00711 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
00712 }
00713 else if(theParticle == theSMinus)
00714 {
00715 xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
00716 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
00717 }
00718 else if(theParticle == theGamma)
00719 {
00720 xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
00721 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
00722
00723 }
00724 else
00725 {
00726 xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
00727 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
00728
00729 xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
00730 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
00731 }
00732 xsection *= millibarn;
00733 return xsection;
00734 }
00735
00736
00738
00739
00740
00741
00742 G4double
00743 G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
00744 const G4Element* anElement)
00745 {
00746 G4int At = G4lrint(anElement->GetN());
00747 G4int Zt = G4lrint(anElement->GetZ());
00748
00749 return GetHadronNucleonXscNS(aParticle, At, Zt);
00750 }
00751
00752
00753
00754
00756
00757
00758
00759
00760 G4double
00761 G4GlauberGribovCrossSection::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
00762 G4int At, G4int Zt)
00763 {
00764 G4double xsection(0);
00765
00766 G4double A0, B0;
00767 G4double hpXscv(0);
00768 G4double hnXscv(0);
00769
00770 G4int Nt = At-Zt;
00771 if (Nt < 0) Nt = 0;
00772
00773 G4double aa = At;
00774 G4double zz = Zt;
00775 G4double nn = Nt;
00776
00777 G4double targ_mass = G4ParticleTable::GetParticleTable()->
00778 GetIonTable()->GetIonMass(Zt, At);
00779
00780 targ_mass = 0.939*GeV;
00781
00782 G4double proj_mass = aParticle->GetMass();
00783 G4double proj_energy = aParticle->GetTotalEnergy();
00784 G4double proj_momentum = aParticle->GetMomentum().mag();
00785
00786 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
00787
00788 sMand /= GeV*GeV;
00789 proj_momentum /= GeV;
00790 proj_energy /= GeV;
00791 proj_mass /= GeV;
00792
00793
00794
00795 G4double s0 = 5.38*5.38;
00796 G4double eta1 = 0.458;
00797 G4double eta2 = 0.458;
00798 G4double B = 0.308;
00799
00800
00801 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00802
00803
00804 if(theParticle == theNeutron)
00805 {
00806 if( proj_momentum >= 373.)
00807 {
00808 return GetHadronNucleonXscPDG(aParticle,At,Zt);
00809 }
00810 else if( proj_momentum >= 10.)
00811
00812 {
00813
00814
00815
00816 if(proj_momentum >= 10.)
00817 {
00818 B0 = 7.5;
00819 A0 = 100. - B0*std::log(3.0e7);
00820
00821 xsection = A0 + B0*std::log(proj_energy) - 11
00822 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00823 0.93827*0.93827,-0.165);
00824 }
00825 xsection *= zz + nn;
00826 }
00827 else
00828 {
00829
00830
00831 if( proj_momentum < 0.73 )
00832 {
00833 hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00834 }
00835 else if( proj_momentum < 1.05 )
00836 {
00837 hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00838 (std::log(proj_momentum/0.73));
00839 }
00840 else
00841 {
00842 hnXscv = 39.0+
00843 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00844 }
00845
00846
00847 if( proj_momentum < 0.8 )
00848 {
00849 hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00850 }
00851 else if( proj_momentum < 1.4 )
00852 {
00853 hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00854 }
00855 else
00856 {
00857 hpXscv = 33.3+
00858 20.8*(std::pow(proj_momentum,2.0)-1.35)/
00859 (std::pow(proj_momentum,2.50)+0.95);
00860 }
00861 xsection = hpXscv*zz + hnXscv*nn;
00862 }
00863 }
00864 else if(theParticle == theProton)
00865 {
00866 if( proj_momentum >= 373.)
00867 {
00868 return GetHadronNucleonXscPDG(aParticle,At,Zt);
00869 }
00870 else if( proj_momentum >= 10.)
00871
00872 {
00873
00874
00875
00876 if(proj_momentum >= 10.)
00877 {
00878 B0 = 7.5;
00879 A0 = 100. - B0*std::log(3.0e7);
00880
00881 xsection = A0 + B0*std::log(proj_energy) - 11
00882 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
00883 0.93827*0.93827,-0.165);
00884 }
00885 xsection *= zz + nn;
00886 }
00887 else
00888 {
00889
00890
00891 if( proj_momentum < 0.73 )
00892 {
00893 hpXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
00894 }
00895 else if( proj_momentum < 1.05 )
00896 {
00897 hpXscv = 23 + 40*(std::log(proj_momentum/0.73))*
00898 (std::log(proj_momentum/0.73));
00899 }
00900 else
00901 {
00902 hpXscv = 39.0+
00903 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
00904 }
00905
00906
00907 if( proj_momentum < 0.8 )
00908 {
00909 hnXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
00910 }
00911 else if( proj_momentum < 1.4 )
00912 {
00913 hnXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
00914 }
00915 else
00916 {
00917 hnXscv = 33.3+
00918 20.8*(std::pow(proj_momentum,2.0)-1.35)/
00919 (std::pow(proj_momentum,2.50)+0.95);
00920 }
00921 xsection = hpXscv*zz + hnXscv*nn;
00922
00923
00924 }
00925
00926 }
00927 else if( theParticle == theAProton )
00928 {
00929
00930
00931
00932
00933
00934
00935 G4double logP = std::log(proj_momentum);
00936
00937 if( proj_momentum <= 1.0 )
00938 {
00939 xsection = zz*(65.55 + 53.84/(proj_momentum+1.e-6) );
00940 }
00941 else
00942 {
00943 xsection = zz*( 41.1 + 77.2*std::pow( proj_momentum, -0.68)
00944 + 0.293*logP*logP - 1.82*logP );
00945 }
00946 if ( nn > 0.)
00947 {
00948 xsection += nn*( 41.9 + 96.2*std::pow( proj_momentum, -0.99) - 0.154*logP);
00949 }
00950 else
00951 {
00952 fInelasticXsc = 38.0 + 38.0*std::pow( proj_momentum, -0.96)
00953 - 0.169*logP*logP;
00954 fInelasticXsc *= millibarn;
00955 }
00956 }
00957 else if( theParticle == thePiPlus )
00958 {
00959 if(proj_momentum < 0.4)
00960 {
00961 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
00962 hpXscv = Ex3+20.0;
00963 }
00964 else if( proj_momentum < 1.15 )
00965 {
00966 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
00967 hpXscv = Ex4+14.0;
00968 }
00969 else if(proj_momentum < 3.5)
00970 {
00971 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
00972 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
00973 hpXscv = Ex1+Ex2+27.5;
00974 }
00975 else
00976 {
00977 hpXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
00978 }
00979
00980
00981 if(proj_momentum < 0.37)
00982 {
00983 hnXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
00984 }
00985 else if(proj_momentum<0.65)
00986 {
00987 hnXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
00988 }
00989 else if(proj_momentum<1.3)
00990 {
00991 hnXscv = 36.1+
00992 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
00993 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
00994 }
00995 else if(proj_momentum<3.0)
00996 {
00997 hnXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
00998 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
00999 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01000 }
01001 else
01002 {
01003 hnXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
01004 }
01005 xsection = hpXscv*zz + hnXscv*nn;
01006 }
01007 else if(theParticle == thePiMinus)
01008 {
01009
01010
01011 if(proj_momentum < 0.4)
01012 {
01013 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
01014 hnXscv = Ex3+20.0;
01015 }
01016 else if(proj_momentum < 1.15)
01017 {
01018 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
01019 hnXscv = Ex4+14.0;
01020 }
01021 else if(proj_momentum < 3.5)
01022 {
01023 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
01024 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
01025 hnXscv = Ex1+Ex2+27.5;
01026 }
01027 else
01028 {
01029 hnXscv = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
01030 }
01031
01032
01033 if(proj_momentum < 0.37)
01034 {
01035 hpXscv = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
01036 }
01037 else if(proj_momentum<0.65)
01038 {
01039 hpXscv = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
01040 }
01041 else if(proj_momentum<1.3)
01042 {
01043 hpXscv = 36.1+
01044 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
01045 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
01046 }
01047 else if(proj_momentum<3.0)
01048 {
01049 hpXscv = 36.1+0.079-4.313*std::log(proj_momentum)+
01050 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
01051 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
01052 }
01053 else
01054 {
01055 hpXscv = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
01056 }
01057 xsection = hpXscv*zz + hnXscv*nn;
01058 }
01059 else if(theParticle == theKPlus)
01060 {
01061 xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
01062 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
01063
01064 xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
01065 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
01066 }
01067 else if(theParticle == theKMinus)
01068 {
01069 xsection = zz*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
01070 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
01071
01072 xsection += nn*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
01073 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
01074 }
01075 else if(theParticle == theSMinus)
01076 {
01077 xsection = aa*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
01078 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
01079 }
01080 else if(theParticle == theGamma)
01081 {
01082 xsection = aa*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
01083 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
01084
01085 }
01086 else
01087 {
01088 xsection = zz*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
01089 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
01090
01091 xsection += nn*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
01092 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
01093 }
01094 xsection *= millibarn;
01095 return xsection;
01096 }
01097
01098 G4double
01099 G4GlauberGribovCrossSection::GetKaonNucleonXscVector(const G4DynamicParticle* aParticle,
01100 G4int At, G4int Zt)
01101 {
01102 G4double Tkin, logTkin, xsc, xscP, xscN;
01103 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01104
01105 G4int Nt = At-Zt;
01106 if (Nt < 0) Nt = 0;
01107
01108 Tkin = aParticle->GetKineticEnergy();
01109
01110 if( Tkin > 70*GeV ) return GetHadronNucleonXscPDG(aParticle,At,Zt);
01111
01112 logTkin = std::log(Tkin);
01113
01114 if( theParticle == theKPlus )
01115 {
01116 xscP = hnXsc->GetKpProtonTotXscVector(logTkin);
01117 xscN = hnXsc->GetKpNeutronTotXscVector(logTkin);
01118 }
01119 else if( theParticle == theKMinus )
01120 {
01121 xscP = hnXsc->GetKmProtonTotXscVector(logTkin);
01122 xscN = hnXsc->GetKmNeutronTotXscVector(logTkin);
01123 }
01124 else
01125 {
01126 xscP = (hnXsc->GetKpProtonTotXscVector(logTkin)+hnXsc->GetKmProtonTotXscVector(logTkin))*0.5;
01127 xscN = (hnXsc->GetKpNeutronTotXscVector(logTkin)+hnXsc->GetKmNeutronTotXscVector(logTkin))*0.5;
01128 }
01129 xsc = xscP*Zt + xscN*Nt;
01130 return xsc;
01131 }
01133
01134
01135
01136 G4double
01137 G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
01138 const G4Element* anElement)
01139 {
01140 G4int At = G4lrint(anElement->GetN());
01141 G4int Zt = G4lrint(anElement->GetZ());
01142
01143 return GetHNinelasticXsc(aParticle, At, Zt);
01144 }
01145
01147
01148
01149
01150 G4double
01151 G4GlauberGribovCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
01152 G4int At, G4int Zt)
01153 {
01154 G4ParticleDefinition* hadron = aParticle->GetDefinition();
01155 G4double sumInelastic;
01156 G4int Nt = At - Zt;
01157 if(Nt < 0) Nt = 0;
01158
01159 if( hadron == theKPlus )
01160 {
01161 sumInelastic = GetHNinelasticXscVU(aParticle, At, Zt);
01162 }
01163 else
01164 {
01165
01166
01167 sumInelastic = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
01168 sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);
01169 }
01170 return sumInelastic;
01171 }
01172
01173
01175
01176
01177
01178 G4double
01179 G4GlauberGribovCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
01180 G4int At, G4int Zt)
01181 {
01182 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
01183 G4int absPDGcode = std::abs(PDGcode);
01184
01185 G4double Elab = aParticle->GetTotalEnergy();
01186
01187 G4double Plab = aParticle->GetMomentum().mag();
01188
01189
01190 Elab /= GeV;
01191 Plab /= GeV;
01192
01193 G4double LogPlab = std::log( Plab );
01194 G4double sqrLogPlab = LogPlab * LogPlab;
01195
01196
01197
01198 G4double NumberOfTargetProtons = G4double(Zt);
01199 G4double NumberOfTargetNucleons = G4double(At);
01200 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
01201
01202 if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
01203
01204 G4double Xtotal, Xelastic, Xinelastic;
01205
01206 if( absPDGcode > 1000 )
01207 {
01208 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
01209 0.522*sqrLogPlab - 4.51*LogPlab;
01210
01211 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
01212 0.513*sqrLogPlab - 4.27*LogPlab;
01213
01214 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
01215 0.169*sqrLogPlab - 1.85*LogPlab;
01216
01217 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
01218 0.169*sqrLogPlab - 1.85*LogPlab;
01219
01220 Xtotal = (NumberOfTargetProtons * XtotPP +
01221 NumberOfTargetNeutrons * XtotPN);
01222
01223 Xelastic = (NumberOfTargetProtons * XelPP +
01224 NumberOfTargetNeutrons * XelPN);
01225 }
01226 else if( PDGcode == 211 )
01227 {
01228 G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01229 0.19 *sqrLogPlab - 0.0 *LogPlab;
01230
01231 G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01232 0.456*sqrLogPlab - 4.03*LogPlab;
01233
01234 G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) +
01235 0.079*sqrLogPlab - 0.0 *LogPlab;
01236
01237 G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) +
01238 0.043*sqrLogPlab - 0.0 *LogPlab;
01239
01240 Xtotal = ( NumberOfTargetProtons * XtotPiP +
01241 NumberOfTargetNeutrons * XtotPiN );
01242
01243 Xelastic = ( NumberOfTargetProtons * XelPiP +
01244 NumberOfTargetNeutrons * XelPiN );
01245 }
01246 else if( PDGcode == -211 )
01247 {
01248 G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) +
01249 0.456*sqrLogPlab - 4.03*LogPlab;
01250
01251 G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) +
01252 0.19 *sqrLogPlab - 0.0 *LogPlab;
01253
01254 G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) +
01255 0.043*sqrLogPlab - 0.0 *LogPlab;
01256
01257 G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) +
01258 0.079*sqrLogPlab - 0.0 *LogPlab;
01259
01260 Xtotal = ( NumberOfTargetProtons * XtotPiP +
01261 NumberOfTargetNeutrons * XtotPiN );
01262
01263 Xelastic = ( NumberOfTargetProtons * XelPiP +
01264 NumberOfTargetNeutrons * XelPiN );
01265 }
01266 else if( PDGcode == 111 )
01267 {
01268 G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) +
01269 0.19 *sqrLogPlab - 0.0 *LogPlab +
01270 33.0 + 14.0 *std::pow(Plab,-1.36) +
01271 0.456*sqrLogPlab - 4.03*LogPlab)/2;
01272
01273 G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) +
01274 0.456*sqrLogPlab - 4.03*LogPlab +
01275 16.4 + 19.3 *std::pow(Plab,-0.42) +
01276 0.19 *sqrLogPlab - 0.0 *LogPlab)/2;
01277
01278 G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) +
01279 0.079*sqrLogPlab - 0.0 *LogPlab +
01280 1.76 + 11.2*std::pow(Plab,-0.64) +
01281 0.043*sqrLogPlab - 0.0 *LogPlab)/2;
01282
01283 G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) +
01284 0.043*sqrLogPlab - 0.0 *LogPlab +
01285 0.0 + 11.4*std::pow(Plab,-0.40) +
01286 0.079*sqrLogPlab - 0.0 *LogPlab)/2;
01287
01288 Xtotal = ( NumberOfTargetProtons * XtotPiP +
01289 NumberOfTargetNeutrons * XtotPiN );
01290
01291 Xelastic = ( NumberOfTargetProtons * XelPiP +
01292 NumberOfTargetNeutrons * XelPiN );
01293 }
01294 else if( PDGcode == 321 )
01295 {
01296 G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) +
01297 0.26 *sqrLogPlab - 1.0 *LogPlab;
01298 G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) +
01299 0.21 *sqrLogPlab - 0.89*LogPlab;
01300
01301 G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
01302 0.16 *sqrLogPlab - 1.3 *LogPlab;
01303
01304 G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) +
01305 0.29 *sqrLogPlab - 2.4 *LogPlab;
01306
01307 Xtotal = ( NumberOfTargetProtons * XtotKP +
01308 NumberOfTargetNeutrons * XtotKN );
01309
01310 Xelastic = ( NumberOfTargetProtons * XelKP +
01311 NumberOfTargetNeutrons * XelKN );
01312 }
01313 else if( PDGcode ==-321 )
01314 {
01315 G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) +
01316 0.66 *sqrLogPlab - 5.6 *LogPlab;
01317 G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) +
01318 0.38 *sqrLogPlab - 2.9 *LogPlab;
01319
01320 G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) +
01321 0.29 *sqrLogPlab - 2.4 *LogPlab;
01322
01323 G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) +
01324 0.16 *sqrLogPlab - 1.3 *LogPlab;
01325
01326 Xtotal = ( NumberOfTargetProtons * XtotKP +
01327 NumberOfTargetNeutrons * XtotKN );
01328
01329 Xelastic = ( NumberOfTargetProtons * XelKP +
01330 NumberOfTargetNeutrons * XelKN );
01331 }
01332 else if( PDGcode == 311 )
01333 {
01334 G4double XtotKP = ( 18.1 + 0. *std::pow(Plab, 0. ) +
01335 0.26 *sqrLogPlab - 1.0 *LogPlab +
01336 32.1 + 0. *std::pow(Plab, 0. ) +
01337 0.66 *sqrLogPlab - 5.6 *LogPlab)/2;
01338
01339 G4double XtotKN = ( 18.7 + 0. *std::pow(Plab, 0. ) +
01340 0.21 *sqrLogPlab - 0.89*LogPlab +
01341 25.2 + 0. *std::pow(Plab, 0. ) +
01342 0.38 *sqrLogPlab - 2.9 *LogPlab)/2;
01343
01344 G4double XelKP = ( 5.0 + 8.1*std::pow(Plab,-1.8 )
01345 + 0.16 *sqrLogPlab - 1.3 *LogPlab +
01346 7.3 + 0. *std::pow(Plab,-0. ) +
01347 0.29 *sqrLogPlab - 2.4 *LogPlab)/2;
01348
01349 G4double XelKN = ( 7.3 + 0. *std::pow(Plab,-0. ) +
01350 0.29 *sqrLogPlab - 2.4 *LogPlab +
01351 5.0 + 8.1*std::pow(Plab,-1.8 ) +
01352 0.16 *sqrLogPlab - 1.3 *LogPlab)/2;
01353
01354 Xtotal = ( NumberOfTargetProtons * XtotKP +
01355 NumberOfTargetNeutrons * XtotKN );
01356
01357 Xelastic = ( NumberOfTargetProtons * XelKP +
01358 NumberOfTargetNeutrons * XelKN );
01359 }
01360 else
01361 {
01362 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
01363 0.522*sqrLogPlab - 4.51*LogPlab;
01364
01365 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
01366 0.513*sqrLogPlab - 4.27*LogPlab;
01367
01368 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
01369 0.169*sqrLogPlab - 1.85*LogPlab;
01370 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
01371 0.169*sqrLogPlab - 1.85*LogPlab;
01372
01373 Xtotal = ( NumberOfTargetProtons * XtotPP +
01374 NumberOfTargetNeutrons * XtotPN );
01375
01376 Xelastic = ( NumberOfTargetProtons * XelPP +
01377 NumberOfTargetNeutrons * XelPN );
01378 }
01379 Xinelastic = Xtotal - Xelastic;
01380
01381 if( Xinelastic < 0.) Xinelastic = 0.;
01382
01383 return Xinelastic*= millibarn;
01384 }
01385
01387
01388
01389
01390 G4double
01391 G4GlauberGribovCrossSection::GetNucleusRadius(const G4DynamicParticle* ,
01392 const G4Element* anElement)
01393 {
01394 G4int At = G4lrint(anElement->GetN());
01395 G4double oneThird = 1.0/3.0;
01396 G4double cubicrAt = std::pow(G4double(At), oneThird);
01397
01398 G4double R;
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414 R = fRadiusConst*cubicrAt;
01415
01416 G4double meanA = 21.;
01417
01418 G4double tauA1 = 40.;
01419 G4double tauA2 = 10.;
01420 G4double tauA3 = 5.;
01421
01422 G4double a1 = 0.85;
01423 G4double b1 = 1. - a1;
01424
01425 G4double b2 = 0.3;
01426 G4double b3 = 4.;
01427
01428 if (At > 20)
01429 {
01430 R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
01431 }
01432 else if (At > 3)
01433 {
01434 R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
01435 }
01436 else
01437 {
01438 R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
01439 }
01440 return R;
01441
01442 }
01444
01445
01446
01447 G4double
01448 G4GlauberGribovCrossSection::GetNucleusRadius(G4int At)
01449 {
01450 G4double oneThird = 1.0/3.0;
01451 G4double cubicrAt = std::pow(G4double(At), oneThird);
01452
01453 G4double R;
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470 R = fRadiusConst*cubicrAt;
01471
01472 G4double meanA = 20.;
01473 G4double tauA = 20.;
01474
01475 if (At > 20)
01476 {
01477 R *= ( 0.8 + 0.2*std::exp( -(G4double(At) - meanA)/tauA) );
01478 }
01479 else
01480 {
01481 R *= ( 1.0 + 0.1*( 1. - std::exp( (G4double(At) - meanA)/tauA) ) );
01482 }
01483
01484 return R;
01485 }
01486
01488
01489
01490
01491 G4double G4GlauberGribovCrossSection::CalculateEcmValue( const G4double mp ,
01492 const G4double mt ,
01493 const G4double Plab )
01494 {
01495 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01496 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
01497
01498
01499
01500 return Ecm ;
01501 }
01502
01504
01505
01506
01507 G4double G4GlauberGribovCrossSection::CalcMandelstamS( const G4double mp ,
01508 const G4double mt ,
01509 const G4double Plab )
01510 {
01511 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
01512 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
01513
01514 return sMand;
01515 }
01516
01518
01519
01520
01521 void G4GlauberGribovCrossSection::CrossSectionDescription(std::ostream& outFile) const
01522 {
01523 outFile << "G4GlauberGribovCrossSection calculates total, inelastic and\n"
01524 << "elastic cross sections for hadron-nucleus cross sections using\n"
01525 << "the Glauber model with Gribov corrections. It is valid for all\n"
01526 << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
01527 << "pi+, pi-, K+, K- and gammas with energies above 3 GeV. This is\n"
01528 << "a cross section component which is to be used to build a cross\n"
01529 << "data set.\n";
01530 }
01531
01532
01533