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