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