| 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 #include "G4Log.hh"
45 #include "G4Lambda.hh" << 40 #include "G4Exp.hh"
46 #include "G4Pow.hh" 41 #include "G4Pow.hh"
47 #include "G4NuclearRadii.hh" <<
48 42
49 ////////////////////////////////////////////// 43 //////////////////////////////////////////////////////////////////////////////
50 // 44 //
51 45
52 G4ComponentGGHadronNucleusXsc::G4ComponentGGHa 46 G4ComponentGGHadronNucleusXsc::G4ComponentGGHadronNucleusXsc()
53 : G4VComponentCrossSection(Default_Name()), 47 : G4VComponentCrossSection(Default_Name()),
54 fTotalXsc(0.0),fElasticXsc(0.0),fInelasticX << 48 // fUpperLimit(100000*GeV),
55 fDiffractionXsc(0.0),fAxsc2piR2(0.0),fModel << 49 fLowerLimit(10.*MeV),// fLowerLimit(3*GeV),
56 fParticle(nullptr),fZ(0),fA(0), fL(0) << 50 fRadiusConst(1.08*fermi), // 1.1, 1.3 ?
>> 51 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
>> 52 fDiffractionXsc(0.0)
>> 53 // , fHadronNucleonXsc(0.0)
57 { 54 {
58 theGamma = G4Gamma::Gamma(); 55 theGamma = G4Gamma::Gamma();
59 theProton = G4Proton::Proton(); 56 theProton = G4Proton::Proton();
60 theNeutron = G4Neutron::Neutron(); 57 theNeutron = G4Neutron::Neutron();
61 theAProton = G4AntiProton::AntiProton(); 58 theAProton = G4AntiProton::AntiProton();
62 theANeutron = G4AntiNeutron::AntiNeutron(); 59 theANeutron = G4AntiNeutron::AntiNeutron();
63 thePiPlus = G4PionPlus::PionPlus(); 60 thePiPlus = G4PionPlus::PionPlus();
64 thePiMinus = G4PionMinus::PionMinus(); 61 thePiMinus = G4PionMinus::PionMinus();
>> 62 thePiZero = G4PionZero::PionZero();
65 theKPlus = G4KaonPlus::KaonPlus(); 63 theKPlus = G4KaonPlus::KaonPlus();
66 theKMinus = G4KaonMinus::KaonMinus(); 64 theKMinus = G4KaonMinus::KaonMinus();
67 theK0S = G4KaonZeroShort::KaonZeroShort 65 theK0S = G4KaonZeroShort::KaonZeroShort();
68 theK0L = G4KaonZeroLong::KaonZeroLong() 66 theK0L = G4KaonZeroLong::KaonZeroLong();
69 theLambda = G4Lambda::Lambda(); << 67 theL = G4Lambda::Lambda();
70 << 68 theAntiL = G4AntiLambda::AntiLambda();
>> 69 theSPlus = G4SigmaPlus::SigmaPlus();
>> 70 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
>> 71 theSMinus = G4SigmaMinus::SigmaMinus();
>> 72 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
>> 73 theS0 = G4SigmaZero::SigmaZero();
>> 74 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
>> 75 theXiMinus = G4XiMinus::XiMinus();
>> 76 theXi0 = G4XiZero::XiZero();
>> 77 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
>> 78 theAXi0 = G4AntiXiZero::AntiXiZero();
>> 79 theOmega = G4OmegaMinus::OmegaMinus();
>> 80 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
>> 81 theD = G4Deuteron::Deuteron();
>> 82 theT = G4Triton::Triton();
>> 83 theA = G4Alpha::Alpha();
>> 84 theHe3 = G4He3::He3();
>> 85
71 hnXsc = new G4HadronNucleonXsc(); 86 hnXsc = new G4HadronNucleonXsc();
72 } 87 }
73 88
74 ////////////////////////////////////////////// << 89 ///////////////////////////////////////////////////////////////////////////////////////
>> 90 //
>> 91 //
75 92
76 G4ComponentGGHadronNucleusXsc::~G4ComponentGGH 93 G4ComponentGGHadronNucleusXsc::~G4ComponentGGHadronNucleusXsc()
77 { 94 {
78 delete hnXsc; << 95 if (hnXsc) delete hnXsc;
79 } 96 }
80 97
81 ////////////////////////////////////////////// << 98 ////////////////////////////////////////////////////////////////////
>> 99
>> 100 G4double G4ComponentGGHadronNucleusXsc::GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
>> 101 G4double kinEnergy,
>> 102 G4int Z, G4int A)
>> 103 {
>> 104 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 105 kinEnergy);
>> 106 fTotalXsc = GetIsoCrossSection(aDP, Z, A);
>> 107 delete aDP;
82 108
83 G4double G4ComponentGGHadronNucleusXsc::GetTot <<
84 const G4ParticleDefinition <<
85 G4double kinEnergy, G4int <<
86 { <<
87 ComputeCrossSections(aParticle, kinEnergy, Z <<
88 return fTotalXsc; 109 return fTotalXsc;
89 } 110 }
90 111
91 ////////////////////////////////////////////// << 112 //////////////////////////////////////////////////////////////////////
>> 113
>> 114 G4double G4ComponentGGHadronNucleusXsc::GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
>> 115 G4double kinEnergy,
>> 116 G4int Z, G4double A)
>> 117 {
>> 118 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 119 kinEnergy);
>> 120 fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
>> 121 delete aDP;
92 122
93 G4double G4ComponentGGHadronNucleusXsc::GetTot <<
94 const G4ParticleDefinition <<
95 G4double kinEnergy, G4int Z, G4int A) <<
96 { <<
97 ComputeCrossSections(aParticle, kinEnergy, Z <<
98 return fTotalXsc; 123 return fTotalXsc;
99 } 124 }
100 125
101 ////////////////////////////////////////////// << 126 ////////////////////////////////////////////////////////////////////
>> 127
>> 128 G4double G4ComponentGGHadronNucleusXsc::GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
>> 129 G4double kinEnergy,
>> 130 G4int Z, G4int A)
>> 131 {
>> 132 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 133 kinEnergy);
>> 134 fTotalXsc = GetIsoCrossSection(aDP, Z, A);
>> 135 delete aDP;
102 136
103 G4double G4ComponentGGHadronNucleusXsc::GetIne <<
104 const G4ParticleDefinition <<
105 G4double kinEnergy, G4int Z, G4double <<
106 { <<
107 ComputeCrossSections(aParticle, kinEnergy, Z <<
108 return fInelasticXsc; 137 return fInelasticXsc;
109 } 138 }
110 139
111 ////////////////////////////////////////////// 140 ////////////////////////////////////////////////////////////////////
112 141
113 G4double G4ComponentGGHadronNucleusXsc::GetIne << 142 G4double G4ComponentGGHadronNucleusXsc::GetProductionIsotopeCrossSection(const G4ParticleDefinition* aParticle,
114 const G4ParticleDefinition << 143 G4double kinEnergy,
115 G4double kinEnergy, G4int Z, G4int A) << 144 G4int Z, G4int A)
116 { << 145 {
117 ComputeCrossSections(aParticle, kinEnergy, Z << 146 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 147 kinEnergy);
>> 148 fTotalXsc = GetIsoCrossSection(aDP, Z, A);
>> 149 delete aDP;
>> 150
>> 151 return fProductionXsc;
>> 152 }
>> 153
>> 154 /////////////////////////////////////////////////////////////////////
>> 155
>> 156 G4double G4ComponentGGHadronNucleusXsc::GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
>> 157 G4double kinEnergy,
>> 158 G4int Z, G4double A)
>> 159 {
>> 160 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 161 kinEnergy);
>> 162 fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
>> 163 delete aDP;
>> 164
118 return fInelasticXsc; 165 return fInelasticXsc;
119 } 166 }
120 167
>> 168 /////////////////////////////////////////////////////////////////////
>> 169
>> 170 G4double G4ComponentGGHadronNucleusXsc::GetProductionElementCrossSection(const G4ParticleDefinition* aParticle,
>> 171 G4double kinEnergy,
>> 172 G4int Z, G4double A)
>> 173 {
>> 174 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 175 kinEnergy);
>> 176 fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
>> 177 delete aDP;
>> 178
>> 179 return fProductionXsc;
>> 180 }
>> 181
121 ////////////////////////////////////////////// 182 //////////////////////////////////////////////////////////////////
122 183
123 G4double G4ComponentGGHadronNucleusXsc::GetEla << 184 G4double G4ComponentGGHadronNucleusXsc::GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
124 const G4ParticleDefinition << 185 G4double kinEnergy,
125 G4double kinEnergy, G4int Z, G4double << 186 G4int Z, G4double A)
126 { << 187 {
127 ComputeCrossSections(aParticle, kinEnergy, Z << 188 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 189 kinEnergy);
>> 190 fTotalXsc = GetIsoCrossSection(aDP, Z, G4int(A));
>> 191 delete aDP;
>> 192
128 return fElasticXsc; 193 return fElasticXsc;
129 } 194 }
130 195
131 ////////////////////////////////////////////// 196 ///////////////////////////////////////////////////////////////////
132 197
133 G4double G4ComponentGGHadronNucleusXsc::GetEla << 198 G4double G4ComponentGGHadronNucleusXsc::GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
134 const G4ParticleDefinition << 199 G4double kinEnergy,
135 G4double kinEnergy, G4int Z, G4int A) << 200 G4int Z, G4int A)
136 { << 201 {
137 ComputeCrossSections(aParticle, kinEnergy, Z << 202 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
>> 203 kinEnergy);
>> 204 fTotalXsc = GetIsoCrossSection(aDP, Z, A);
>> 205 delete aDP;
>> 206
138 return fElasticXsc; 207 return fElasticXsc;
139 } 208 }
140 209
141 ////////////////////////////////////////////// 210 ////////////////////////////////////////////////////////////////
142 211
143 G4double G4ComponentGGHadronNucleusXsc::Comput << 212 G4double G4ComponentGGHadronNucleusXsc::ComputeQuasiElasticRatio(const G4ParticleDefinition* aParticle,
144 const G4ParticleDefinition << 213 G4double kinEnergy,
145 G4double kinEnergy, G4int Z, G4int A) << 214 G4int Z, G4int A)
146 { << 215 {
147 ComputeCrossSections(aParticle, kinEnergy, Z << 216 G4DynamicParticle* aDP = new G4DynamicParticle(aParticle,G4ParticleMomentum(1.,0.,0.),
148 G4double ratio = (fInelasticXsc > 0.) << 217 kinEnergy);
149 ? (fInelasticXsc - fProductionXsc)/fInelas << 218 fTotalXsc = GetIsoCrossSection(aDP, Z, A);
150 ratio = std::max(ratio, 0.); << 219 delete aDP;
>> 220 G4double ratio = 0.;
>> 221
>> 222 if(fInelasticXsc > 0.)
>> 223 {
>> 224 ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
>> 225 if(ratio < 0.) ratio = 0.;
>> 226 }
151 return ratio; 227 return ratio;
152 } 228 }
>> 229
153 230
154 ////////////////////////////////////////////// <<
155 231
156 G4double G4ComponentGGHadronNucleusXsc::GetPro <<
157 const G4ParticleDefinition <<
158 G4double kinEnergy, G4int Z, G4double <<
159 { <<
160 ComputeCrossSections(aParticle, kinEnergy, Z <<
161 return fProductionXsc; <<
162 } <<
163 232
164 ////////////////////////////////////////////// << 233 ////////////////////////////////////////////////////////////////////////////////////////
165 234
166 G4double G4ComponentGGHadronNucleusXsc::GetPro << 235 G4bool
167 const G4ParticleDefinition << 236 G4ComponentGGHadronNucleusXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
168 G4double kinEnergy, G4int Z, G4int A) << 237 G4int Z, G4int /*A*/,
169 { << 238 const G4Element*,
170 ComputeCrossSections(aParticle, kinEnergy, Z << 239 const G4Material*)
171 return fProductionXsc; << 240 {
>> 241 G4bool applicable = false;
>> 242 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
>> 243 G4double kineticEnergy = aDP->GetKineticEnergy();
>> 244
>> 245 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
>> 246
>> 247 if (
>> 248 Z >= 1 // >= H for kaons
>> 249 &&
>> 250 (
>> 251 kineticEnergy >= fLowerLimit
>> 252 &&
>> 253 // Z > 1 && // >= He
>> 254 (
>> 255 theParticle == theAProton ||
>> 256 theParticle == theGamma ||
>> 257 theParticle == theSMinus ||
>> 258 theParticle == theProton ||
>> 259 theParticle == theNeutron ||
>> 260 theParticle == thePiPlus ||
>> 261 theParticle == thePiMinus
>> 262 )
>> 263 )
>> 264 )
>> 265 applicable = true;
>> 266
>> 267 if (
>> 268 Z >= 1 // >= H for kaons
>> 269 &&
>> 270 (
>> 271 kineticEnergy >= 0.01*fLowerLimit
>> 272 &&
>> 273 (
>> 274 theParticle == theKPlus ||
>> 275 theParticle == theKMinus ||
>> 276 theParticle == theK0L ||
>> 277 theParticle == theK0S
>> 278 )
>> 279 )
>> 280 )
>> 281 applicable = true;
>> 282
>> 283 return applicable;
172 } 284 }
173 285
174 ////////////////////////////////////////////// << 286 ////////////////////////////////////////////////////////////////////////////////////////
175 // 287 //
176 // Calculates total and inelastic Xsc, derives << 288 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
177 // inelastic accordong to Glauber model with G << 289 // Glauber model with Gribov correction calculated in the dipole approximation on
178 // in the dipole approximation on light cone. << 290 // 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 291 // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
181 292
182 void G4ComponentGGHadronNucleusXsc::ComputeCro << 293 G4double
183 const G4ParticleDefinition* aP << 294 G4ComponentGGHadronNucleusXsc::GetIsoCrossSection(const G4DynamicParticle* aParticle,
184 G4double kinEnergy, G4int Z, G << 295 G4int Z, G4int A,
185 { << 296 const G4Isotope*,
186 // check cache << 297 const G4Element*,
187 if(aParticle == fParticle && fZ == Z && fA = << 298 const G4Material*)
188 { return; } << 299 {
189 fParticle = aParticle; << 300 G4double xsection, sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
190 fZ = Z; << 301 G4double hpInXsc(0.), hnInXsc(0.);
191 fA = A; << 302 G4double R = GetNucleusRadius(A);
192 fL = nL; <<
193 fEnergy = kinEnergy; <<
194 G4Pow* pG4Pow=G4Pow::GetInstance(); <<
195 // <<
196 G4double cofInelastic = 2.4; <<
197 static const G4double cofTotal = 2.0; <<
198 G4double sigma(0.0), hpInXsc(0.0), hnInXsc(0 <<
199 303
200 G4int N = std::max(A - Z, 0); // number of << 304 G4int N = A - Z; // number of neutrons
>> 305 if (N < 0) N = 0;
201 306
202 if( aParticle == theKPlus || aParticle == th << 307 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
203 aParticle == theK0S || aParticle == th << 308
>> 309 if( theParticle == theProton ||
>> 310 theParticle == theNeutron ||
>> 311 theParticle == thePiPlus ||
>> 312 theParticle == thePiMinus )
204 { 313 {
205 sigma = (1 == Z) << 314 // sigma = GetHadronNucleonXscNS(aParticle, A, Z);
206 ? hnXsc->KaonNucleonXscNS(aParticle, the << 315
207 : Z*hnXsc->KaonNucleonXscGG(aParticle, t << 316 sigma = Z*hnXsc->GetHadronNucleonXscNS(aParticle, theProton);
>> 317
208 hpInXsc = hnXsc->GetInelasticHadronNucleon 318 hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
209 319
210 if(N > 0) { << 320 sigma += N*hnXsc->GetHadronNucleonXscNS(aParticle, theNeutron);
211 sigma += N*hnXsc->KaonNucleonXscGG(aPart << 321
212 hnInXsc = hnXsc->GetInelasticHadronNucle << 322 hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
213 } << 323
214 R = G4NuclearRadii::RadiusKNGG(A); << 324 cofInelastic = 2.4;
215 cofInelastic = 2.2; << 325 cofTotal = 2.0;
216 } 326 }
217 else << 327 else if( theParticle == theKPlus ||
>> 328 theParticle == theKMinus ||
>> 329 theParticle == theK0S ||
>> 330 theParticle == theK0L )
218 { 331 {
219 sigma = Z*hnXsc->HadronNucleonXsc(aParticl << 332 // sigma = GetKaonNucleonXscVector(aParticle, A, Z);
>> 333
>> 334 sigma = Z*hnXsc->GetKaonNucleonXscGG(aParticle, theProton);
>> 335
220 hpInXsc = hnXsc->GetInelasticHadronNucleon 336 hpInXsc = hnXsc->GetInelasticHadronNucleonXsc();
221 337
222 if(N > 0) { << 338 sigma += N*hnXsc->GetKaonNucleonXscGG(aParticle, theNeutron);
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 } <<
236 339
237 G4double nucleusSquare = cofTotal*pi*R*R; << 340 hnInXsc = hnXsc->GetInelasticHadronNucleonXsc();
238 G4double ratio = sigma/nucleusSquare; << 341
239 G4double difratio = ratio/(1.+ratio); << 342 cofInelastic = 2.2;
240 fDiffractionXsc = 0.5*nucleusSquare*( difrat << 343 cofTotal = 2.0;
>> 344 R = 1.3*fermi;
>> 345 R *= G4Pow::GetInstance()->powA(G4double(A), 0.3333);
>> 346 }
>> 347 else
>> 348 {
>> 349 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
>> 350 cofInelastic = 2.2;
>> 351 cofTotal = 2.0;
>> 352 }
>> 353 // cofInelastic = 2.0;
241 354
242 if( A > 1 ) 355 if( A > 1 )
243 { 356 {
244 fTotalXsc = nucleusSquare*G4Log( 1. + rati << 357 nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
245 *GetParticleBarCorTot(aParticle, Z); << 358 ratio = sigma/nucleusSquare;
246 359
247 // inelastic xsc << 360 xsection = nucleusSquare*G4Log( 1. + ratio );
248 fAxsc2piR2 = cofInelastic*ratio; << 361
249 fModelInLog = G4Log( 1. + fAxsc2piR2 ); << 362 xsection *= GetParticleBarCorTot(theParticle, Z);
250 fInelasticXsc = nucleusSquare*fModelInLog/ << 363
251 G4double barcorr = GetParticleBarCorIn(aPa << 364 fTotalXsc = xsection;
252 fInelasticXsc *= barcorr; << 365
253 fElasticXsc = std::max(fTotalXsc - fInelas << 366
254 << 367
255 G4double xratio = (Z*hpInXsc + N*hnInXsc)/ << 368 fInelasticXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )/cofInelastic;
256 fProductionXsc = << 369
257 nucleusSquare*G4Log(1. + cofInelastic*xr << 370 fInelasticXsc *= GetParticleBarCorIn(theParticle, Z);
258 fProductionXsc = std::min(fProductionXsc, << 371
>> 372 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 373
>> 374 if(fElasticXsc < 0.) fElasticXsc = 0.;
>> 375
>> 376 G4double difratio = ratio/(1.+ratio);
>> 377
>> 378 fDiffractionXsc = 0.5*nucleusSquare*( difratio - G4Log( 1. + difratio ) );
>> 379
>> 380
>> 381 // sigma = GetHNinelasticXsc(aParticle, A, Z);
>> 382
>> 383 sigma = Z*hpInXsc + N*hnInXsc;
>> 384
>> 385 ratio = sigma/nucleusSquare;
>> 386
>> 387 fProductionXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )/cofInelastic;
>> 388
>> 389 fProductionXsc *= GetParticleBarCorIn(theParticle, Z);
>> 390
>> 391 if (fElasticXsc < 0.) fElasticXsc = 0.;
259 } 392 }
260 else // H 393 else // H
261 { 394 {
262 fTotalXsc = sigma; 395 fTotalXsc = sigma;
263 fInelasticXsc = hpInXsc; << 396 xsection = sigma;
264 fElasticXsc = std::max(fTotalXsc - fInel << 397
265 fProductionXsc = fInelasticXsc; << 398 fInelasticXsc = hnXsc->GetInelasticHadronNucleonXsc();
266 fDiffractionXsc = 0.2*fInelasticXsc; << 399
267 // G4double xratio = hpInXsc/nucleusSquare << 400 if ( theParticle != theAProton )
268 // fProductionXsc = nucleusSquare*G4Log(1. << 401 {
269 // fProductionXsc = std::min(fProductionXs << 402 fElasticXsc = hnXsc->GetElasticHadronNucleonXsc();
>> 403
>> 404 // sigma = GetHNinelasticXsc(aParticle, A, Z);
>> 405 // fInelasticXsc = sigma;
>> 406 // fElasticXsc = fTotalXsc - fInelasticXsc;
>> 407 }
>> 408 else if( theParticle == theKPlus ||
>> 409 theParticle == theKMinus ||
>> 410 theParticle == theK0S ||
>> 411 theParticle == theK0L )
>> 412 {
>> 413 fInelasticXsc = hpInXsc;
>> 414 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 415 }
>> 416 else
>> 417 {
>> 418 fInelasticXsc = hpInXsc;
>> 419 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 420 }
>> 421 if (fElasticXsc < 0.) fElasticXsc = 0.;
>> 422
270 } 423 }
271 /* << 424 return xsection;
272 G4cout << "GGXsc: Z= " << Z << " A= " << A < <<
273 << " xtot(b)= " << fTotalXsc/barn <<
274 << " xel(b)= " << fElasticXsc/barn << " xi <<
275 << G4endl; <<
276 */ <<
277 } 425 }
278 426
279 ////////////////////////////////////////////// 427 //////////////////////////////////////////////////////////////////////////
280 // 428 //
281 // Return single-diffraction/inelastic cross-s 429 // Return single-diffraction/inelastic cross-section ratio
282 430
283 G4double G4ComponentGGHadronNucleusXsc::GetRat << 431 G4double G4ComponentGGHadronNucleusXsc::
284 const G4DynamicParticle* aParticle, G << 432 GetRatioSD(const G4DynamicParticle* aParticle, G4int A, G4int Z)
285 { 433 {
286 ComputeCrossSections(aParticle->GetDefinitio << 434 G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
287 aParticle->GetKineticEnergy(), Z, A << 435 G4double R = GetNucleusRadius(A);
288 436
289 return (fInelasticXsc > 0.0) ? fDiffractionX << 437 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 438
>> 439 if( theParticle == theProton ||
>> 440 theParticle == theNeutron ||
>> 441 theParticle == thePiPlus ||
>> 442 theParticle == thePiMinus )
>> 443 {
>> 444 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
>> 445 cofInelastic = 2.4;
>> 446 cofTotal = 2.0;
>> 447 }
>> 448 else
>> 449 {
>> 450 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
>> 451 cofInelastic = 2.2;
>> 452 cofTotal = 2.0;
>> 453 }
>> 454 nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
>> 455 ratio = sigma/nucleusSquare;
>> 456
>> 457 fInelasticXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )/cofInelastic;
>> 458
>> 459 G4double difratio = ratio/(1.+ratio);
>> 460
>> 461 fDiffractionXsc = 0.5*nucleusSquare*( difratio - G4Log( 1. + difratio ) );
>> 462
>> 463 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
>> 464 else ratio = 0.;
>> 465
>> 466 return ratio;
290 } 467 }
291 468
292 ////////////////////////////////////////////// 469 //////////////////////////////////////////////////////////////////////////
293 // 470 //
294 // Return quasi-elastic/inelastic cross-sectio << 471 // Return suasi-elastic/inelastic cross-section ratio
295 472
296 G4double G4ComponentGGHadronNucleusXsc:: 473 G4double G4ComponentGGHadronNucleusXsc::
297 GetRatioQE(const G4DynamicParticle* aParticle, 474 GetRatioQE(const G4DynamicParticle* aParticle, G4int A, G4int Z)
298 { 475 {
299 ComputeCrossSections(aParticle->GetDefinitio << 476 G4double sigma, cofInelastic, cofTotal, nucleusSquare, ratio;
300 aParticle->GetKineticEnergy(), Z, A << 477 G4double R = GetNucleusRadius(A);
>> 478
>> 479 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 480
>> 481 if( theParticle == theProton ||
>> 482 theParticle == theNeutron ||
>> 483 theParticle == thePiPlus ||
>> 484 theParticle == thePiMinus )
>> 485 {
>> 486 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
>> 487 cofInelastic = 2.4;
>> 488 cofTotal = 2.0;
>> 489 }
>> 490 else
>> 491 {
>> 492 sigma = GetHadronNucleonXscNS(aParticle, A, Z);
>> 493 cofInelastic = 2.2;
>> 494 cofTotal = 2.0;
>> 495 }
>> 496 nucleusSquare = cofTotal*pi*R*R; // basically 2piRR
>> 497 ratio = sigma/nucleusSquare;
>> 498
>> 499 fInelasticXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )/cofInelastic;
>> 500
>> 501 sigma = GetHNinelasticXsc(aParticle, A, Z);
>> 502 ratio = sigma/nucleusSquare;
>> 503
>> 504 fProductionXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )/cofInelastic;
301 505
302 return (fInelasticXsc > std::max(fProduction << 506 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
303 ? 1.0 - fProductionXsc/fInelasticXsc : 0.0 << 507 else ratio = 0.;
>> 508 if ( ratio < 0. ) ratio = 0.;
>> 509
>> 510 return ratio;
304 } 511 }
305 512
306 ////////////////////////////////////////////// 513 /////////////////////////////////////////////////////////////////////////////////////
307 // 514 //
308 // Returns hadron-nucleon total Xsc according << 515 // Returns hadron-nucleon Xsc according to differnt parametrisations:
>> 516 // [2] E. Levin, hep-ph/9710546
>> 517 // [3] U. Dersch, et al, hep-ex/9910052
>> 518 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
309 519
310 G4double G4ComponentGGHadronNucleusXsc::GetHad << 520 G4double
311 const G4DynamicParticle* aParticle, c << 521 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
>> 522 const G4Element* anElement)
312 { 523 {
313 G4int At = G4lrint(anElement->GetN()); // n 524 G4int At = G4lrint(anElement->GetN()); // number of nucleons
314 G4int Zt = anElement->GetZasInt(); // numbe << 525 G4int Zt = G4lrint(anElement->GetZ()); // number of protons
315 526
316 return GetHadronNucleonXsc(aParticle, At, Zt 527 return GetHadronNucleonXsc(aParticle, At, Zt);
317 } 528 }
318 529
319 ////////////////////////////////////////////// 530 /////////////////////////////////////////////////////////////////////////////////////
320 // 531 //
321 // Returns hadron-nucleon total Xsc according << 532 // Returns hadron-nucleon Xsc according to differnt parametrisations:
>> 533 // [2] E. Levin, hep-ph/9710546
>> 534 // [3] U. Dersch, et al, hep-ex/9910052
>> 535 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
322 536
323 G4double G4ComponentGGHadronNucleusXsc::GetHad << 537 G4double
324 const G4DynamicParticle* aParticle, G << 538 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
>> 539 G4int At, G4int /*Zt*/)
325 { 540 {
326 return hnXsc->HadronNucleonXscEL(aParticle-> << 541 G4double xsection;
327 aParticle->GetKineticEnergy()); << 542
>> 543 //G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
>> 544
>> 545 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 546
>> 547 G4double proj_mass = aParticle->GetMass();
>> 548 G4double proj_momentum = aParticle->GetMomentum().mag();
>> 549 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 550
>> 551 sMand /= GeV*GeV; // in GeV for parametrisation
>> 552 proj_momentum /= GeV;
>> 553
>> 554 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 555
>> 556 G4double aa = At;
>> 557
>> 558 if(theParticle == theGamma)
>> 559 {
>> 560 xsection = aa*(0.0677*G4Pow::GetInstance()->powA(sMand,0.0808) + 0.129*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 561 }
>> 562 else if(theParticle == theNeutron) // as proton ???
>> 563 {
>> 564 xsection = aa*(21.70*G4Pow::GetInstance()->powA(sMand,0.0808) + 56.08*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 565 }
>> 566 else if(theParticle == theProton)
>> 567 {
>> 568 xsection = aa*(21.70*G4Pow::GetInstance()->powA(sMand,0.0808) + 56.08*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 569 // xsection = At*( 49.51*G4Pow::GetInstance()->powA(sMand,-0.097) + 0.314*G4Log(sMand)*G4Log(sMand) );
>> 570 // xsection = At*( 38.4 + 0.85*std::abs(G4Pow::GetInstance()->powA(log(sMand),1.47)) );
>> 571 }
>> 572 else if(theParticle == theAProton)
>> 573 {
>> 574 xsection = aa*( 21.70*G4Pow::GetInstance()->powA(sMand,0.0808) + 98.39*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 575 }
>> 576 else if(theParticle == thePiPlus)
>> 577 {
>> 578 xsection = aa*(13.63*G4Pow::GetInstance()->powA(sMand,0.0808) + 27.56*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 579 }
>> 580 else if(theParticle == thePiMinus)
>> 581 {
>> 582 // xsection = At*( 55.2*G4Pow::GetInstance()->powA(sMand,-0.255) + 0.346*G4Log(sMand)*G4Log(sMand) );
>> 583 xsection = aa*(13.63*G4Pow::GetInstance()->powA(sMand,0.0808) + 36.02*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 584 }
>> 585 else if(theParticle == theKPlus)
>> 586 {
>> 587 xsection = aa*(11.82*G4Pow::GetInstance()->powA(sMand,0.0808) + 8.15*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 588 }
>> 589 else if(theParticle == theKMinus)
>> 590 {
>> 591 xsection = aa*(11.82*G4Pow::GetInstance()->powA(sMand,0.0808) + 26.36*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 592 }
>> 593 else // as proton ???
>> 594 {
>> 595 xsection = aa*(21.70*G4Pow::GetInstance()->powA(sMand,0.0808) + 56.08*G4Pow::GetInstance()->powA(sMand,-0.4525));
>> 596 }
>> 597 xsection *= millibarn;
>> 598 return xsection;
328 } 599 }
329 600
>> 601
330 ////////////////////////////////////////////// 602 /////////////////////////////////////////////////////////////////////////////////////
331 // 603 //
332 // Returns hadron-nucleon total Xsc according << 604 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
333 // << 605 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
334 606
335 G4double G4ComponentGGHadronNucleusXsc::GetHad << 607 G4double
336 const G4DynamicParticle* aParticle, c << 608 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
>> 609 const G4Element* anElement)
337 { 610 {
338 G4int At = G4lrint(anElement->GetN()); // n 611 G4int At = G4lrint(anElement->GetN()); // number of nucleons
339 G4int Zt = anElement->GetZasInt(); // n << 612 G4int Zt = G4lrint(anElement->GetZ()); // number of protons
340 613
341 return GetHadronNucleonXscPDG(aParticle, At, 614 return GetHadronNucleonXscPDG(aParticle, At, Zt);
342 } 615 }
343 616
>> 617
>> 618
>> 619
344 ////////////////////////////////////////////// 620 /////////////////////////////////////////////////////////////////////////////////////
345 // 621 //
346 // Returns hadron-nucleon total Xsc according << 622 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
347 // << 623 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
>> 624 // At = number of nucleons, Zt = number of protons
348 625
349 G4double G4ComponentGGHadronNucleusXsc::GetHad << 626 G4double
350 const G4DynamicParticle* aParticle, G << 627 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
>> 628 G4int At, G4int Zt)
351 { 629 {
352 G4double res = 0.0; << 630 G4double xsection;
353 if(1 == At && 1 == Zt) { << 631
354 res = hnXsc->HadronNucleonXscPDG(aParticle << 632 G4int Nt = At-Zt; // number of neutrons
355 aParticle->GetKineticEnergy()); << 633 if (Nt < 0) Nt = 0;
356 } else if(1 == At && 0 == Zt) { << 634
357 res = hnXsc->HadronNucleonXscPDG(aParticle << 635 G4double zz = Zt;
358 aParticle->GetKineticEnergy()); << 636 G4double aa = At;
359 } else { << 637 G4double nn = Nt;
360 ComputeCrossSections(aParticle->GetDefinit << 638
361 aParticle->GetKineticEnergy(), Zt, At); << 639 G4double targ_mass = G4ParticleTable::GetParticleTable()->
362 res = fTotalXsc; << 640 GetIonTable()->GetIonMass(Zt, At);
>> 641
>> 642 targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 643
>> 644 G4double proj_mass = aParticle->GetMass();
>> 645 G4double proj_momentum = aParticle->GetMomentum().mag();
>> 646
>> 647 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 648
>> 649 sMand /= GeV*GeV; // in GeV for parametrisation
>> 650
>> 651 // General PDG fit constants
>> 652
>> 653 G4double s0 = 5.38*5.38; // in Gev^2
>> 654 G4double eta1 = 0.458;
>> 655 G4double eta2 = 0.458;
>> 656 G4double B = 0.308;
>> 657
>> 658
>> 659 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 660
>> 661
>> 662 if(theParticle == theNeutron) // proton-neutron fit
>> 663 {
>> 664 xsection = zz*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 665 + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) - 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 666 xsection += nn*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 667 + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) - 33.34*G4Pow::GetInstance()->powA(sMand,-eta2)); // pp for nn
>> 668 }
>> 669 else if(theParticle == theProton)
>> 670 {
>> 671
>> 672 xsection = zz*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 673 + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) - 33.34*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 674
>> 675 xsection += nn*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 676 + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) - 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 677 }
>> 678 else if(theParticle == theAProton)
>> 679 {
>> 680 xsection = zz*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 681 + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) + 33.34*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 682
>> 683 xsection += nn*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 684 + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) + 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 685 }
>> 686 else if(theParticle == thePiPlus)
>> 687 {
>> 688 xsection = aa*( 20.86 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 689 + 19.24*G4Pow::GetInstance()->powA(sMand,-eta1) - 6.03*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 690 }
>> 691 else if(theParticle == thePiMinus)
>> 692 {
>> 693 xsection = aa*( 20.86 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 694 + 19.24*G4Pow::GetInstance()->powA(sMand,-eta1) + 6.03*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 695 }
>> 696 else if(theParticle == theKPlus || theParticle == theK0L )
>> 697 {
>> 698 xsection = zz*( 17.91 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 699 + 7.14*G4Pow::GetInstance()->powA(sMand,-eta1) - 13.45*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 700
>> 701 xsection += nn*( 17.87 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 702 + 5.17*G4Pow::GetInstance()->powA(sMand,-eta1) - 7.23*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 703 }
>> 704 else if(theParticle == theKMinus || theParticle == theK0S )
>> 705 {
>> 706 xsection = zz*( 17.91 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 707 + 7.14*G4Pow::GetInstance()->powA(sMand,-eta1) + 13.45*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 708
>> 709 xsection += nn*( 17.87 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 710 + 5.17*G4Pow::GetInstance()->powA(sMand,-eta1) + 7.23*G4Pow::GetInstance()->powA(sMand,-eta2));
363 } 711 }
364 return res; << 712 else if(theParticle == theSMinus)
>> 713 {
>> 714 xsection = aa*( 35.20 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 715 - 199.*G4Pow::GetInstance()->powA(sMand,-eta1) + 264.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 716 }
>> 717 else if(theParticle == theGamma) // modify later on
>> 718 {
>> 719 xsection = aa*( 0.0 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 720 + 0.032*G4Pow::GetInstance()->powA(sMand,-eta1) - 0.0*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 721
>> 722 }
>> 723 else // as proton ???
>> 724 {
>> 725 xsection = zz*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 726 + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) - 33.34*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 727
>> 728 xsection += nn*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 729 + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) - 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 730 }
>> 731 xsection *= millibarn; // parametrised in mb
>> 732 return xsection;
365 } 733 }
366 734
>> 735
367 ////////////////////////////////////////////// 736 /////////////////////////////////////////////////////////////////////////////////////
368 // 737 //
369 // Returns hadron-nucleon total cross-section << 738 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
370 // data from mainly http://wwwppds.ihep.su:800 739 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
371 740
372 G4double G4ComponentGGHadronNucleusXsc::GetHad << 741 G4double
373 const G4DynamicParticle* aParticle, c << 742 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
>> 743 const G4Element* anElement)
374 { 744 {
375 G4int At = G4lrint(anElement->GetN()); // n 745 G4int At = G4lrint(anElement->GetN()); // number of nucleons
376 G4int Zt = anElement->GetZasInt(); // n << 746 G4int Zt = G4lrint(anElement->GetZ()); // number of protons
377 747
378 return GetHadronNucleonXscNS(aParticle, At, 748 return GetHadronNucleonXscNS(aParticle, At, Zt);
379 } 749 }
380 750
381 ////////////////////////////////////////////// << 751
>> 752
>> 753
>> 754 /////////////////////////////////////////////////////////////////////////////////////
382 // 755 //
383 // Returns hadron-nucleon total cross-section << 756 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
384 // data from mainly http://wwwppds.ihep.su:800 757 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
385 758
386 G4double G4ComponentGGHadronNucleusXsc::GetHad << 759 G4double
387 const G4DynamicParticle* aParticle, G << 760 G4ComponentGGHadronNucleusXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
>> 761 G4int At, G4int Zt)
388 { 762 {
389 G4double res = 0.0; << 763 G4double xsection(0);
390 if(1 == At && 1 == Zt) { << 764 // G4double Delta; DHW 19 May 2011: variable set but not used
391 res = hnXsc->HadronNucleonXscNS(aParticle- << 765 G4double A0, B0;
392 aParticle->GetKineticEnergy()); << 766 G4double hpXscv(0);
393 } else if(1 == At && 0 == Zt) { << 767 G4double hnXscv(0);
394 res = hnXsc->HadronNucleonXscNS(aParticle- << 768
395 aParticle->GetKineticEnergy()); << 769 G4int Nt = At-Zt; // number of neutrons
>> 770 if (Nt < 0) Nt = 0;
>> 771
>> 772 G4double aa = At;
>> 773 G4double zz = Zt;
>> 774 G4double nn = Nt;
>> 775
>> 776 G4double targ_mass = G4ParticleTable::GetParticleTable()->
>> 777 GetIonTable()->GetIonMass(Zt, At);
>> 778
>> 779 targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 780
>> 781 G4double proj_mass = aParticle->GetMass();
>> 782 G4double proj_energy = aParticle->GetTotalEnergy();
>> 783 G4double proj_momentum = aParticle->GetMomentum().mag();
>> 784
>> 785 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 786
>> 787 sMand /= GeV*GeV; // in GeV for parametrisation
>> 788 proj_momentum /= GeV;
>> 789 proj_energy /= GeV;
>> 790 proj_mass /= GeV;
>> 791
>> 792 // General PDG fit constants
>> 793
>> 794 G4double s0 = 5.38*5.38; // in Gev^2
>> 795 G4double eta1 = 0.458;
>> 796 G4double eta2 = 0.458;
>> 797 G4double B = 0.308;
>> 798
>> 799
>> 800 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 801
>> 802
>> 803 if(theParticle == theNeutron)
>> 804 {
>> 805 if( proj_momentum >= 373.)
>> 806 {
>> 807 return GetHadronNucleonXscPDG(aParticle,At,Zt);
>> 808 }
>> 809 else if( proj_momentum >= 10.)
>> 810 // if( proj_momentum >= 2.)
>> 811 {
>> 812 // Delta = 1.; // DHW 19 May 2011: variable set but not used
>> 813 // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
>> 814
>> 815 if(proj_momentum >= 10.)
>> 816 {
>> 817 B0 = 7.5;
>> 818 A0 = 100. - B0*G4Log(3.0e7);
>> 819
>> 820 xsection = A0 + B0*G4Log(proj_energy) - 11
>> 821 + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
>> 822 0.93827*0.93827,-0.165); // mb
>> 823 }
>> 824 xsection *= zz + nn;
>> 825 }
>> 826 else
>> 827 {
>> 828 // nn to be pp
>> 829
>> 830 if( proj_momentum < 0.73 )
>> 831 {
>> 832 hnXscv = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
>> 833 }
>> 834 else if( proj_momentum < 1.05 )
>> 835 {
>> 836 hnXscv = 23 + 40*(G4Log(proj_momentum/0.73))*
>> 837 (G4Log(proj_momentum/0.73));
>> 838 }
>> 839 else // if( proj_momentum < 10. )
>> 840 {
>> 841 hnXscv = 39.0+
>> 842 75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
>> 843 }
>> 844 // pn to be np
>> 845
>> 846 if( proj_momentum < 0.8 )
>> 847 {
>> 848 hpXscv = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
>> 849 }
>> 850 else if( proj_momentum < 1.4 )
>> 851 {
>> 852 hpXscv = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
>> 853 }
>> 854 else // if( proj_momentum < 10. )
>> 855 {
>> 856 hpXscv = 33.3+
>> 857 20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
>> 858 (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
>> 859 }
>> 860 xsection = hpXscv*zz + hnXscv*nn;
>> 861 }
>> 862 }
>> 863 else if(theParticle == theProton)
>> 864 {
>> 865 if( proj_momentum >= 373.)
>> 866 {
>> 867 return GetHadronNucleonXscPDG(aParticle,At,Zt);
>> 868 }
>> 869 else if( proj_momentum >= 10.)
>> 870 // if( proj_momentum >= 2.)
>> 871 {
>> 872 // Delta = 1.; DHW 19 May 2011: variable set but not used
>> 873 // if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
>> 874
>> 875 if(proj_momentum >= 10.)
>> 876 {
>> 877 B0 = 7.5;
>> 878 A0 = 100. - B0*G4Log(3.0e7);
>> 879
>> 880 xsection = A0 + B0*G4Log(proj_energy) - 11
>> 881 + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
>> 882 0.93827*0.93827,-0.165); // mb
>> 883 }
>> 884 xsection *= zz + nn;
>> 885 }
>> 886 else
>> 887 {
>> 888 // pp
>> 889
>> 890 if( proj_momentum < 0.73 )
>> 891 {
>> 892 hpXscv = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
>> 893 }
>> 894 else if( proj_momentum < 1.05 )
>> 895 {
>> 896 hpXscv = 23 + 40*(G4Log(proj_momentum/0.73))*
>> 897 (G4Log(proj_momentum/0.73));
>> 898 }
>> 899 else // if( proj_momentum < 10. )
>> 900 {
>> 901 hpXscv = 39.0+
>> 902 75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
>> 903 }
>> 904 // pn to be np
>> 905
>> 906 if( proj_momentum < 0.8 )
>> 907 {
>> 908 hnXscv = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
>> 909 }
>> 910 else if( proj_momentum < 1.4 )
>> 911 {
>> 912 hnXscv = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
>> 913 }
>> 914 else // if( proj_momentum < 10. )
>> 915 {
>> 916 hnXscv = 33.3+
>> 917 20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
>> 918 (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
>> 919 }
>> 920 xsection = hpXscv*zz + hnXscv*nn;
>> 921 // xsection = hpXscv*(Zt + Nt);
>> 922 // xsection = hnXscv*(Zt + Nt);
>> 923 }
>> 924 // xsection *= 0.95;
>> 925 }
>> 926 else if( theParticle == theAProton )
>> 927 {
>> 928 // xsection = Zt*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 929 // + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) + 33.34*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 930
>> 931 // xsection += Nt*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 932 // + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) + 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 933
>> 934 G4double logP = G4Log(proj_momentum);
>> 935
>> 936 if( proj_momentum <= 1.0 )
>> 937 {
>> 938 xsection = zz*(65.55 + 53.84/(proj_momentum+1.e-6) );
>> 939 }
>> 940 else
>> 941 {
>> 942 xsection = zz*( 41.1 + 77.2*G4Pow::GetInstance()->powA( proj_momentum, -0.68)
>> 943 + 0.293*logP*logP - 1.82*logP );
>> 944 }
>> 945 if ( nn > 0.)
>> 946 {
>> 947 xsection += nn*( 41.9 + 96.2*G4Pow::GetInstance()->powA( proj_momentum, -0.99) - 0.154*logP);
>> 948 }
>> 949 else // H
>> 950 {
>> 951 fInelasticXsc = 38.0 + 38.0*G4Pow::GetInstance()->powA( proj_momentum, -0.96)
>> 952 - 0.169*logP*logP;
>> 953 fInelasticXsc *= millibarn;
>> 954 }
>> 955 }
>> 956 else if( theParticle == thePiPlus )
>> 957 {
>> 958 if(proj_momentum < 0.4)
>> 959 {
>> 960 G4double Ex3 = 180*G4Exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
>> 961 hpXscv = Ex3+20.0;
>> 962 }
>> 963 else if( proj_momentum < 1.15 )
>> 964 {
>> 965 G4double Ex4 = 88*(G4Log(proj_momentum/0.75))*(G4Log(proj_momentum/0.75));
>> 966 hpXscv = Ex4+14.0;
>> 967 }
>> 968 else if(proj_momentum < 3.5)
>> 969 {
>> 970 G4double Ex1 = 3.2*G4Exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
>> 971 G4double Ex2 = 12*G4Exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
>> 972 hpXscv = Ex1+Ex2+27.5;
>> 973 }
>> 974 else // if(proj_momentum > 3.5) // mb
>> 975 {
>> 976 hpXscv = 10.6+2.*G4Log(proj_energy)+25*G4Pow::GetInstance()->powA(proj_energy,-0.43);
>> 977 }
>> 978 // pi+n = pi-p??
>> 979
>> 980 if(proj_momentum < 0.37)
>> 981 {
>> 982 hnXscv = 28.0 + 40*G4Exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
>> 983 }
>> 984 else if(proj_momentum<0.65)
>> 985 {
>> 986 hnXscv = 26+110*(G4Log(proj_momentum/0.48))*(G4Log(proj_momentum/0.48));
>> 987 }
>> 988 else if(proj_momentum<1.3)
>> 989 {
>> 990 hnXscv = 36.1+
>> 991 10*G4Exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
>> 992 24*G4Exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
>> 993 }
>> 994 else if(proj_momentum<3.0)
>> 995 {
>> 996 hnXscv = 36.1+0.079-4.313*G4Log(proj_momentum)+
>> 997 3*G4Exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
>> 998 1.5*G4Exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
>> 999 }
>> 1000 else // mb
>> 1001 {
>> 1002 hnXscv = 10.6+2*G4Log(proj_energy)+30*G4Pow::GetInstance()->powA(proj_energy,-0.43);
>> 1003 }
>> 1004 xsection = hpXscv*zz + hnXscv*nn;
>> 1005 }
>> 1006 else if(theParticle == thePiMinus)
>> 1007 {
>> 1008 // pi-n = pi+p??
>> 1009
>> 1010 if(proj_momentum < 0.4)
>> 1011 {
>> 1012 G4double Ex3 = 180*G4Exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
>> 1013 hnXscv = Ex3+20.0;
>> 1014 }
>> 1015 else if(proj_momentum < 1.15)
>> 1016 {
>> 1017 G4double Ex4 = 88*(G4Log(proj_momentum/0.75))*(G4Log(proj_momentum/0.75));
>> 1018 hnXscv = Ex4+14.0;
>> 1019 }
>> 1020 else if(proj_momentum < 3.5)
>> 1021 {
>> 1022 G4double Ex1 = 3.2*G4Exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
>> 1023 G4double Ex2 = 12*G4Exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
>> 1024 hnXscv = Ex1+Ex2+27.5;
>> 1025 }
>> 1026 else // if(proj_momentum > 3.5) // mb
>> 1027 {
>> 1028 hnXscv = 10.6+2.*G4Log(proj_energy)+25*G4Pow::GetInstance()->powA(proj_energy,-0.43);
>> 1029 }
>> 1030 // pi-p
>> 1031
>> 1032 if(proj_momentum < 0.37)
>> 1033 {
>> 1034 hpXscv = 28.0 + 40*G4Exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
>> 1035 }
>> 1036 else if(proj_momentum<0.65)
>> 1037 {
>> 1038 hpXscv = 26+110*(G4Log(proj_momentum/0.48))*(G4Log(proj_momentum/0.48));
>> 1039 }
>> 1040 else if(proj_momentum<1.3)
>> 1041 {
>> 1042 hpXscv = 36.1+
>> 1043 10*G4Exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
>> 1044 24*G4Exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
>> 1045 }
>> 1046 else if(proj_momentum<3.0)
>> 1047 {
>> 1048 hpXscv = 36.1+0.079-4.313*G4Log(proj_momentum)+
>> 1049 3*G4Exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
>> 1050 1.5*G4Exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
>> 1051 }
>> 1052 else // mb
>> 1053 {
>> 1054 hpXscv = 10.6+2*G4Log(proj_energy)+30*G4Pow::GetInstance()->powA(proj_energy,-0.43);
>> 1055 }
>> 1056 xsection = hpXscv*zz + hnXscv*nn;
>> 1057 }
>> 1058 else if(theParticle == theKPlus)
>> 1059 {
>> 1060 xsection = zz*( 17.91 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1061 + 7.14*G4Pow::GetInstance()->powA(sMand,-eta1) - 13.45*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1062
>> 1063 xsection += nn*( 17.87 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1064 + 5.17*G4Pow::GetInstance()->powA(sMand,-eta1) - 7.23*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1065 }
>> 1066 else if(theParticle == theKMinus)
>> 1067 {
>> 1068 xsection = zz*( 17.91 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1069 + 7.14*G4Pow::GetInstance()->powA(sMand,-eta1) + 13.45*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1070
>> 1071 xsection += nn*( 17.87 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1072 + 5.17*G4Pow::GetInstance()->powA(sMand,-eta1) + 7.23*G4Pow::GetInstance()->powA(sMand,-eta2));
396 } 1073 }
397 return res; << 1074 else if(theParticle == theSMinus)
>> 1075 {
>> 1076 xsection = aa*( 35.20 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1077 - 199.*G4Pow::GetInstance()->powA(sMand,-eta1) + 264.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1078 }
>> 1079 else if(theParticle == theGamma) // modify later on
>> 1080 {
>> 1081 xsection = aa*( 0.0 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1082 + 0.032*G4Pow::GetInstance()->powA(sMand,-eta1) - 0.0*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1083
>> 1084 }
>> 1085 else // as proton ???
>> 1086 {
>> 1087 xsection = zz*( 35.45 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1088 + 42.53*G4Pow::GetInstance()->powA(sMand,-eta1) - 33.34*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1089
>> 1090 xsection += nn*( 35.80 + B*G4Pow::GetInstance()->powA(G4Log(sMand/s0),2.)
>> 1091 + 40.15*G4Pow::GetInstance()->powA(sMand,-eta1) - 30.*G4Pow::GetInstance()->powA(sMand,-eta2));
>> 1092 }
>> 1093 xsection *= millibarn; // parametrised in mb
>> 1094 return xsection;
>> 1095 }
>> 1096
>> 1097 /*
>> 1098 G4double
>> 1099 G4ComponentGGHadronNucleusXsc::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 = G4Log(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;
398 } 1131 }
>> 1132 */
399 1133
400 ////////////////////////////////////////////// 1134 /////////////////////////////////////////////////////////////////////////////////////
401 // 1135 //
402 // Returns hadron-nucleon inelastic cross-sect 1136 // Returns hadron-nucleon inelastic cross-section based on proper parametrisation
403 1137
404 G4double 1138 G4double
405 G4ComponentGGHadronNucleusXsc::GetHNinelasticX 1139 G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
406 << 1140 const G4Element* anElement)
407 { 1141 {
408 G4int At = G4lrint(anElement->GetN()); // n 1142 G4int At = G4lrint(anElement->GetN()); // number of nucleons
409 G4int Zt = anElement->GetZasInt(); // n << 1143 G4int Zt = G4lrint(anElement->GetZ()); // number of protons
410 1144
411 return GetHNinelasticXsc(aParticle, At, Zt); 1145 return GetHNinelasticXsc(aParticle, At, Zt);
412 } 1146 }
413 1147
414 ////////////////////////////////////////////// 1148 /////////////////////////////////////////////////////////////////////////////////////
415 // 1149 //
416 // Returns hadron-nucleon inelastic cross-sect << 1150 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
417 1151
418 G4double G4ComponentGGHadronNucleusXsc::GetHNi << 1152 G4double
419 const G4DynamicParticle* aParticle, G << 1153 G4ComponentGGHadronNucleusXsc::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
>> 1154 G4int At, G4int Zt)
420 { 1155 {
421 const G4ParticleDefinition* hadron = aPartic 1156 const G4ParticleDefinition* hadron = aParticle->GetDefinition();
422 G4double e = aParticle->GetKineticEnergy(); << 1157 G4double sumInelastic;
423 G4int Nt = std::max(At - Zt, 0); << 1158 G4int Nt = At - Zt;
424 << 1159 if(Nt < 0) Nt = 0;
425 hnXsc->HadronNucleonXscNS(hadron, theProton, << 1160
426 G4double sumInelastic = Zt*hnXsc->GetInelast << 1161 if( hadron == theKPlus )
427 if(Nt > 0) { << 1162 {
428 hnXsc->HadronNucleonXscNS(hadron, theNeutr << 1163 sumInelastic = GetHNinelasticXscVU(aParticle, At, Zt);
429 sumInelastic += Nt*hnXsc->GetInelasticHadr <<
430 } 1164 }
>> 1165 else
>> 1166 {
>> 1167 //sumInelastic = Zt*GetHadronNucleonXscMK(aParticle, theProton);
>> 1168 // sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);
>> 1169 sumInelastic = G4double(Zt)*GetHadronNucleonXscNS(aParticle, 1, 1);
>> 1170 sumInelastic += G4double(Nt)*GetHadronNucleonXscNS(aParticle, 1, 0);
>> 1171 }
431 return sumInelastic; 1172 return sumInelastic;
432 } 1173 }
433 1174
>> 1175
434 ////////////////////////////////////////////// 1176 /////////////////////////////////////////////////////////////////////////////////////
435 // 1177 //
436 // Returns hadron-nucleon inelastic cross-sect 1178 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
437 1179
438 G4double G4ComponentGGHadronNucleusXsc::GetHNi << 1180 G4double
439 const G4DynamicParticle* aParticle, G << 1181 G4ComponentGGHadronNucleusXsc::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
>> 1182 G4int At, G4int Zt)
440 { 1183 {
441 const G4ParticleDefinition* hadron = aPartic << 1184 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
442 G4double e = aParticle->GetKineticEnergy(); << 1185 G4int absPDGcode = std::abs(PDGcode);
443 G4int Nt = std::max(At - Zt, 0); << 1186
>> 1187 G4double Elab = aParticle->GetTotalEnergy();
>> 1188 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
>> 1189 G4double Plab = aParticle->GetMomentum().mag();
>> 1190 // std::sqrt(Elab * Elab - 0.88);
>> 1191
>> 1192 Elab /= GeV;
>> 1193 Plab /= GeV;
>> 1194
>> 1195 G4double LogPlab = G4Log( Plab );
>> 1196 G4double sqrLogPlab = LogPlab * LogPlab;
444 1197
445 hnXsc->HadronNucleonXscVU(hadron, theProton, << 1198 //G4cout<<"Plab = "<<Plab<<G4endl;
446 G4double sumInelastic = Zt*hnXsc->GetInelast << 1199
447 if(Nt > 0) { << 1200 G4double NumberOfTargetProtons = G4double(Zt);
448 hnXsc->HadronNucleonXscVU(hadron, theNeutr << 1201 G4double NumberOfTargetNucleons = G4double(At);
449 sumInelastic += Nt*hnXsc->GetInelasticHadr << 1202 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
>> 1203
>> 1204 if(NumberOfTargetNeutrons < 0.0) NumberOfTargetNeutrons = 0.0;
>> 1205
>> 1206 G4double Xtotal, Xelastic, Xinelastic;
>> 1207
>> 1208 if( absPDGcode > 1000 ) //------Projectile is baryon --------
>> 1209 {
>> 1210 G4double XtotPP = 48.0 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1211 0.522*sqrLogPlab - 4.51*LogPlab;
>> 1212
>> 1213 G4double XtotPN = 47.3 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1214 0.513*sqrLogPlab - 4.27*LogPlab;
>> 1215
>> 1216 G4double XelPP = 11.9 + 26.9*G4Pow::GetInstance()->powA(Plab,-1.21) +
>> 1217 0.169*sqrLogPlab - 1.85*LogPlab;
>> 1218
>> 1219 G4double XelPN = 11.9 + 26.9*G4Pow::GetInstance()->powA(Plab,-1.21) +
>> 1220 0.169*sqrLogPlab - 1.85*LogPlab;
>> 1221
>> 1222 Xtotal = (NumberOfTargetProtons * XtotPP +
>> 1223 NumberOfTargetNeutrons * XtotPN);
>> 1224
>> 1225 Xelastic = (NumberOfTargetProtons * XelPP +
>> 1226 NumberOfTargetNeutrons * XelPN);
450 } 1227 }
451 return sumInelastic; << 1228 else if( PDGcode == 211 ) //------Projectile is PionPlus -------
>> 1229 {
>> 1230 G4double XtotPiP = 16.4 + 19.3 *G4Pow::GetInstance()->powA(Plab,-0.42) +
>> 1231 0.19 *sqrLogPlab - 0.0 *LogPlab;
>> 1232
>> 1233 G4double XtotPiN = 33.0 + 14.0 *G4Pow::GetInstance()->powA(Plab,-1.36) +
>> 1234 0.456*sqrLogPlab - 4.03*LogPlab;
>> 1235
>> 1236 G4double XelPiP = 0.0 + 11.4*G4Pow::GetInstance()->powA(Plab,-0.40) +
>> 1237 0.079*sqrLogPlab - 0.0 *LogPlab;
>> 1238
>> 1239 G4double XelPiN = 1.76 + 11.2*G4Pow::GetInstance()->powA(Plab,-0.64) +
>> 1240 0.043*sqrLogPlab - 0.0 *LogPlab;
>> 1241
>> 1242 Xtotal = ( NumberOfTargetProtons * XtotPiP +
>> 1243 NumberOfTargetNeutrons * XtotPiN );
>> 1244
>> 1245 Xelastic = ( NumberOfTargetProtons * XelPiP +
>> 1246 NumberOfTargetNeutrons * XelPiN );
>> 1247 }
>> 1248 else if( PDGcode == -211 ) //------Projectile is PionMinus -------
>> 1249 {
>> 1250 G4double XtotPiP = 33.0 + 14.0 *G4Pow::GetInstance()->powA(Plab,-1.36) +
>> 1251 0.456*sqrLogPlab - 4.03*LogPlab;
>> 1252
>> 1253 G4double XtotPiN = 16.4 + 19.3 *G4Pow::GetInstance()->powA(Plab,-0.42) +
>> 1254 0.19 *sqrLogPlab - 0.0 *LogPlab;
>> 1255
>> 1256 G4double XelPiP = 1.76 + 11.2*G4Pow::GetInstance()->powA(Plab,-0.64) +
>> 1257 0.043*sqrLogPlab - 0.0 *LogPlab;
>> 1258
>> 1259 G4double XelPiN = 0.0 + 11.4*G4Pow::GetInstance()->powA(Plab,-0.40) +
>> 1260 0.079*sqrLogPlab - 0.0 *LogPlab;
>> 1261
>> 1262 Xtotal = ( NumberOfTargetProtons * XtotPiP +
>> 1263 NumberOfTargetNeutrons * XtotPiN );
>> 1264
>> 1265 Xelastic = ( NumberOfTargetProtons * XelPiP +
>> 1266 NumberOfTargetNeutrons * XelPiN );
>> 1267 }
>> 1268 else if( PDGcode == 111 ) //------Projectile is PionZero -------
>> 1269 {
>> 1270 G4double XtotPiP =(16.4 + 19.3 *G4Pow::GetInstance()->powA(Plab,-0.42) +
>> 1271 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
>> 1272 33.0 + 14.0 *G4Pow::GetInstance()->powA(Plab,-1.36) +
>> 1273 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
>> 1274
>> 1275 G4double XtotPiN =(33.0 + 14.0 *G4Pow::GetInstance()->powA(Plab,-1.36) +
>> 1276 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
>> 1277 16.4 + 19.3 *G4Pow::GetInstance()->powA(Plab,-0.42) +
>> 1278 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
>> 1279
>> 1280 G4double XelPiP =( 0.0 + 11.4*G4Pow::GetInstance()->powA(Plab,-0.40) +
>> 1281 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
>> 1282 1.76 + 11.2*G4Pow::GetInstance()->powA(Plab,-0.64) +
>> 1283 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
>> 1284
>> 1285 G4double XelPiN =( 1.76 + 11.2*G4Pow::GetInstance()->powA(Plab,-0.64) +
>> 1286 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
>> 1287 0.0 + 11.4*G4Pow::GetInstance()->powA(Plab,-0.40) +
>> 1288 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
>> 1289
>> 1290 Xtotal = ( NumberOfTargetProtons * XtotPiP +
>> 1291 NumberOfTargetNeutrons * XtotPiN );
>> 1292
>> 1293 Xelastic = ( NumberOfTargetProtons * XelPiP +
>> 1294 NumberOfTargetNeutrons * XelPiN );
>> 1295 }
>> 1296 else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
>> 1297 {
>> 1298 G4double XtotKP = 18.1 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1299 0.26 *sqrLogPlab - 1.0 *LogPlab;
>> 1300 G4double XtotKN = 18.7 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1301 0.21 *sqrLogPlab - 0.89*LogPlab;
>> 1302
>> 1303 G4double XelKP = 5.0 + 8.1*G4Pow::GetInstance()->powA(Plab,-1.8 ) +
>> 1304 0.16 *sqrLogPlab - 1.3 *LogPlab;
>> 1305
>> 1306 G4double XelKN = 7.3 + 0. *G4Pow::GetInstance()->powA(Plab,-0. ) +
>> 1307 0.29 *sqrLogPlab - 2.4 *LogPlab;
>> 1308
>> 1309 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 1310 NumberOfTargetNeutrons * XtotKN );
>> 1311
>> 1312 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 1313 NumberOfTargetNeutrons * XelKN );
>> 1314 }
>> 1315 else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
>> 1316 {
>> 1317 G4double XtotKP = 32.1 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1318 0.66 *sqrLogPlab - 5.6 *LogPlab;
>> 1319 G4double XtotKN = 25.2 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1320 0.38 *sqrLogPlab - 2.9 *LogPlab;
>> 1321
>> 1322 G4double XelKP = 7.3 + 0. *G4Pow::GetInstance()->powA(Plab,-0. ) +
>> 1323 0.29 *sqrLogPlab - 2.4 *LogPlab;
>> 1324
>> 1325 G4double XelKN = 5.0 + 8.1*G4Pow::GetInstance()->powA(Plab,-1.8 ) +
>> 1326 0.16 *sqrLogPlab - 1.3 *LogPlab;
>> 1327
>> 1328 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 1329 NumberOfTargetNeutrons * XtotKN );
>> 1330
>> 1331 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 1332 NumberOfTargetNeutrons * XelKN );
>> 1333 }
>> 1334 else if( PDGcode == 311 ) //------Projectile is KaonZero ------
>> 1335 {
>> 1336 G4double XtotKP = ( 18.1 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1337 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
>> 1338 32.1 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1339 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
>> 1340
>> 1341 G4double XtotKN = ( 18.7 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1342 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
>> 1343 25.2 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1344 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
>> 1345
>> 1346 G4double XelKP = ( 5.0 + 8.1*G4Pow::GetInstance()->powA(Plab,-1.8 )
>> 1347 + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
>> 1348 7.3 + 0. *G4Pow::GetInstance()->powA(Plab,-0. ) +
>> 1349 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
>> 1350
>> 1351 G4double XelKN = ( 7.3 + 0. *G4Pow::GetInstance()->powA(Plab,-0. ) +
>> 1352 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
>> 1353 5.0 + 8.1*G4Pow::GetInstance()->powA(Plab,-1.8 ) +
>> 1354 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
>> 1355
>> 1356 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 1357 NumberOfTargetNeutrons * XtotKN );
>> 1358
>> 1359 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 1360 NumberOfTargetNeutrons * XelKN );
>> 1361 }
>> 1362 else //------Projectile is undefined, Nucleon assumed
>> 1363 {
>> 1364 G4double XtotPP = 48.0 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1365 0.522*sqrLogPlab - 4.51*LogPlab;
>> 1366
>> 1367 G4double XtotPN = 47.3 + 0. *G4Pow::GetInstance()->powA(Plab, 0. ) +
>> 1368 0.513*sqrLogPlab - 4.27*LogPlab;
>> 1369
>> 1370 G4double XelPP = 11.9 + 26.9*G4Pow::GetInstance()->powA(Plab,-1.21) +
>> 1371 0.169*sqrLogPlab - 1.85*LogPlab;
>> 1372 G4double XelPN = 11.9 + 26.9*G4Pow::GetInstance()->powA(Plab,-1.21) +
>> 1373 0.169*sqrLogPlab - 1.85*LogPlab;
>> 1374
>> 1375 Xtotal = ( NumberOfTargetProtons * XtotPP +
>> 1376 NumberOfTargetNeutrons * XtotPN );
>> 1377
>> 1378 Xelastic = ( NumberOfTargetProtons * XelPP +
>> 1379 NumberOfTargetNeutrons * XelPN );
>> 1380 }
>> 1381 Xinelastic = Xtotal - Xelastic;
>> 1382
>> 1383 if( Xinelastic < 0.) Xinelastic = 0.;
>> 1384
>> 1385 return Xinelastic*= millibarn;
452 } 1386 }
453 1387
454 ////////////////////////////////////////////// << 1388 ////////////////////////////////////////////////////////////////////////////////////
455 // 1389 //
456 // 1390 //
457 1391
458 void G4ComponentGGHadronNucleusXsc::Descriptio << 1392 G4double
>> 1393 G4ComponentGGHadronNucleusXsc::GetNucleusRadius(const G4DynamicParticle* ,
>> 1394 const G4Element* anElement)
>> 1395 {
>> 1396 G4int At = G4lrint(anElement->GetN());
>> 1397 G4double oneThird = 1.0/3.0;
>> 1398 G4double cubicrAt = G4Pow::GetInstance()->powA(G4double(At), oneThird);
>> 1399
>> 1400 G4double R; // = fRadiusConst*cubicrAt;
>> 1401 /*
>> 1402 G4double tmp = G4Pow::GetInstance()->powA( cubicrAt-1., 3.);
>> 1403 tmp += At;
>> 1404 tmp *= 0.5;
>> 1405
>> 1406 if (At > 20.) // 20.
>> 1407 {
>> 1408 R = fRadiusConst*G4Pow::GetInstance()->powA (tmp, oneThird);
>> 1409 }
>> 1410 else
>> 1411 {
>> 1412 R = fRadiusConst*cubicrAt;
>> 1413 }
>> 1414 */
>> 1415
>> 1416 R = fRadiusConst*cubicrAt;
>> 1417
>> 1418 G4double meanA = 21.;
>> 1419
>> 1420 G4double tauA1 = 40.;
>> 1421 G4double tauA2 = 10.;
>> 1422 G4double tauA3 = 5.;
>> 1423
>> 1424 G4double a1 = 0.85;
>> 1425 G4double b1 = 1. - a1;
>> 1426
>> 1427 G4double b2 = 0.3;
>> 1428 G4double b3 = 4.;
>> 1429
>> 1430 if (At > 20) // 20.
>> 1431 {
>> 1432 R *= ( a1 + b1*G4Exp( -(At - meanA)/tauA1) );
>> 1433 }
>> 1434 else if (At > 3)
>> 1435 {
>> 1436 R *= ( 1.0 + b2*( 1. - G4Exp( (At - meanA)/tauA2) ) );
>> 1437 }
>> 1438 else
>> 1439 {
>> 1440 R *= ( 1.0 + b3*( 1. - G4Exp( (At - meanA)/tauA3) ) );
>> 1441 }
>> 1442 return R;
>> 1443
>> 1444 }
>> 1445 ////////////////////////////////////////////////////////////////////////////////////
>> 1446 //
>> 1447 //
>> 1448
>> 1449 G4double
>> 1450 G4ComponentGGHadronNucleusXsc::GetNucleusRadius(G4int At)
>> 1451 {
>> 1452 G4double oneThird = 1.0/3.0;
>> 1453 G4double cubicrAt = G4Pow::GetInstance()->powA(G4double(At), oneThird);
>> 1454
>> 1455 G4double R; // = fRadiusConst*cubicrAt;
>> 1456
>> 1457 /*
>> 1458 G4double tmp = G4Pow::GetInstance()->powA( cubicrAt-1., 3.);
>> 1459 tmp += At;
>> 1460 tmp *= 0.5;
>> 1461
>> 1462 if (At > 20.)
>> 1463 {
>> 1464 R = fRadiusConst*G4Pow::GetInstance()->powA (tmp, oneThird);
>> 1465 }
>> 1466 else
>> 1467 {
>> 1468 R = fRadiusConst*cubicrAt;
>> 1469 }
>> 1470 */
>> 1471
>> 1472 R = fRadiusConst*cubicrAt;
>> 1473
>> 1474 G4double meanA = 20.;
>> 1475 G4double tauA = 20.;
>> 1476
>> 1477 if (At > 20) // 20.
>> 1478 {
>> 1479 R *= ( 0.8 + 0.2*G4Exp( -(G4double(At) - meanA)/tauA) );
>> 1480 }
>> 1481 else
>> 1482 {
>> 1483 R *= ( 1.0 + 0.1*( 1. - G4Exp( (G4double(At) - meanA)/tauA) ) );
>> 1484 }
>> 1485
>> 1486 return R;
>> 1487 }
>> 1488
>> 1489 ////////////////////////////////////////////////////////////////////////////////////
>> 1490 //
>> 1491 //
>> 1492
>> 1493 G4double G4ComponentGGHadronNucleusXsc::CalculateEcmValue( const G4double mp ,
>> 1494 const G4double mt ,
>> 1495 const G4double Plab )
>> 1496 {
>> 1497 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
>> 1498 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
>> 1499 // G4double Pcm = Plab * mt / Ecm;
>> 1500 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
>> 1501
>> 1502 return Ecm ; // KEcm;
>> 1503 }
>> 1504
>> 1505 ////////////////////////////////////////////////////////////////////////////////////
>> 1506 //
>> 1507 //
>> 1508
>> 1509 G4double G4ComponentGGHadronNucleusXsc::CalcMandelstamS( const G4double mp ,
>> 1510 const G4double mt ,
>> 1511 const G4double Plab )
>> 1512 {
>> 1513 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
>> 1514 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
>> 1515
>> 1516 return sMand;
>> 1517 }
>> 1518
>> 1519 ////////////////////////////////////////////////////////////////////////////////////
>> 1520 //
>> 1521 //
>> 1522
>> 1523 void G4ComponentGGHadronNucleusXsc::CrossSectionDescription(std::ostream& outFile) const
459 { 1524 {
460 outFile << "G4ComponentGGHadronNucleusXsc ca 1525 outFile << "G4ComponentGGHadronNucleusXsc calculates total, inelastic and\n"
461 << "elastic cross sections for hadro 1526 << "elastic cross sections for hadron-nucleus cross sections using\n"
462 << "the Glauber model with Gribov co 1527 << "the Glauber model with Gribov corrections. It is valid for all\n"
463 << "targets except hydrogen, and for 1528 << "targets except hydrogen, and for incident p, pbar, n, sigma-,\n"
464 << "pi+, pi-, K+, K- and gammas with 1529 << "pi+, pi-, K+, K- and gammas with energies above 3 GeV. This is\n"
465 << "a cross section component which 1530 << "a cross section component which is to be used to build a cross\n"
466 << "data set.\n"; 1531 << "data set.\n";
467 } 1532 }
468 1533
469 1534
470 ////////////////////////////////////////////// 1535 ///////////////////////////////////////////////////////////////////////////////
471 // 1536 //
472 // Correction arrays for GG <-> Bar changea at 1537 // Correction arrays for GG <-> Bar changea at ~ 90 GeV
473 1538
474 const G4double G4ComponentGGHadronNucleusXsc:: 1539 const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionTot[93] = {
475 1540
476 1.0, 1.0, 1.42517e+00, // 1.118517e+00, 1541 1.0, 1.0, 1.42517e+00, // 1.118517e+00,
477 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.06 1542 1.082002e+00, 1.116171e+00, 1.078747e+00, 1.061315e+00,
478 1.058205e+00, 1.082663e+00, 1.068500e+00, 1.07 1543 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 1544 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 1545 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 1546 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 1547 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 1548 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 1549 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 1550 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 1551 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 1552 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 1553 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 1554 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 1555 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 1556 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 1557 1.037075e+00, 1.034721e+00
493 1558
494 }; 1559 };
495 1560
496 const G4double G4ComponentGGHadronNucleusXsc:: 1561 const G4double G4ComponentGGHadronNucleusXsc::fNeutronBarCorrectionIn[93] = {
497 1562
498 1.0, 1.0, 1.167421e+00, 1.156250e+00, 1.20 1563 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 1564 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 1565 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 1566 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, 1567 1.131515e+00, 1.144338e+00, // 1.130338e+00,
503 1.134171e+00, 1.139206e+00, 1.148474e+00, // 1 1568 1.134171e+00, 1.139206e+00, 1.148474e+00, // 1.141474e+00,
504 1.142189e+00, 1569 1.142189e+00,
505 1.140725e+00, 1.140100e+00, 1.139848e+00, 1.13 1570 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 1571 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 1572 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 1573 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 1574 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 1575 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 1576 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 1577 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 1578 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 1579 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 1580 1.046435e+00, 1.042614e+00
516 1581
517 }; 1582 };
518 1583
519 const G4double G4ComponentGGHadronNucleusXsc:: 1584 const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionTot[93] = {
520 1585
521 1.0, 1.0, 1586 1.0, 1.0,
522 1.118515e+00, 1.082000e+00, 1.116169e+00, 1.07 1587 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 1588 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 1589 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 1590 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 1591 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 1592 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 1593 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 1594 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 1595 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 1596 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 1597 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 1598 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 1599 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 1600 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 1601 1.031265e+00, 1.032990e+00, 1.033283e+00, 1.035014e+00, 1.033944e+00, 1.037074e+00,
537 1.034720e+00 1602 1.034720e+00
538 1603
539 }; 1604 };
540 1605
541 const G4double G4ComponentGGHadronNucleusXsc:: 1606 const G4double G4ComponentGGHadronNucleusXsc::fProtonBarCorrectionIn[93] = {
542 1607
543 1.0, 1.0, 1608 1.0, 1.0,
544 1.147419e+00, // 1.167419e+00, 1609 1.147419e+00, // 1.167419e+00,
545 1.156248e+00, 1.205362e+00, 1.154224e+00, 1.12 1610 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 1611 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 1612 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 1613 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, 1614 1.140337e+00, // 1.130337e+00,
550 1615
551 1.134170e+00, 1.139205e+00, 1.151472e+00, // 1616 1.134170e+00, 1.139205e+00, 1.151472e+00, // 1.141472e+00,
552 1.142188e+00, 1.140724e+00, 1617 1.142188e+00, 1.140724e+00,
553 1.140099e+00, 1.139847e+00, 1.137672e+00, 1.13 1618 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 1619 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 1620 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 1621 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 1622 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 1623 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 1624 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 1625 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 1626 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 1627 1.051919e+00, 1.051257e+00, 1.049472e+00, 1.048822e+00, 1.045983e+00, 1.046434e+00,
563 1.042613e+00 1628 1.042613e+00
564 1629
565 }; 1630 };
566 1631
567 1632
568 const G4double G4ComponentGGHadronNucleusXsc:: 1633 const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionTot[93] = {
569 1634
570 1.0, 1.0, 1635 1.0, 1.0,
571 1.075927e+00, 1.074407e+00, 1.126098e+00, 1.10 1636 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 1637 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 1638 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 1639 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 1640 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 1641 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 1642 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 1643 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 1644 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 1645 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 1646 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 1647 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 1648 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 1649 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 1650 1.158462e+00, 1.158014e+00, 1.156519e+00, 1.156081e+00, 1.153602e+00, 1.154190e+00,
586 1.152974e+00 1651 1.152974e+00
587 1652
588 }; 1653 };
589 1654
590 const G4double G4ComponentGGHadronNucleusXsc:: 1655 const G4double G4ComponentGGHadronNucleusXsc::fPionPlusBarCorrectionIn[93] = {
591 1656
592 1.0, 1.0, 1657 1.0, 1.0,
593 1.140246e+00, 1.097872e+00, 1.104301e+00, 1.06 1658 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 1659 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 1660 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 1661 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 1662 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 1663 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 1664 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 1665 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 1666 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 1667 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 1668 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 1669 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 1670 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 1671 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 1672 1.062192e+00, 1.061967e+00, 1.061180e+00, 1.060960e+00, 1.059646e+00, 1.059975e+00,
608 1.059658e+00 1673 1.059658e+00
609 1674
610 }; 1675 };
611 1676
612 1677
613 const G4double G4ComponentGGHadronNucleusXsc:: 1678 const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionTot[93] = {
614 1679
615 1.0, 1.0, 1680 1.0, 1.0,
616 1.3956e+00, 1.077959e+00, 1.129145e+00, 1.1020 1681 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 1682 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 1683 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 1684 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 1685 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 1686 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 1687 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 1688 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 1689 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 1690 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 1691 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 1692 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 1693 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 1694 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 1695 1.162774e+00, 1.162217e+00, 1.160740e+00, 1.160196e+00, 1.157857e+00, 1.158220e+00,
631 1.157267e+00 1696 1.157267e+00
632 }; 1697 };
>> 1698
633 1699
634 const G4double G4ComponentGGHadronNucleusXsc:: 1700 const G4double G4ComponentGGHadronNucleusXsc::fPionMinusBarCorrectionIn[93] = {
635 1701
636 1.0, 1.0, 1702 1.0, 1.0,
637 1.463e+00, 1.100898e+00, 1.106773e+00, 1.07 1703 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 1704 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 1705 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 1706 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 1707 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 1708 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 1709 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 1710 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 1711 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 1712 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 1713 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 1714 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 1715 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 1716 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 1717 1.065241e+00, 1.064939e+00, 1.064166e+00, 1.063872e+00, 1.062659e+00, 1.062828e+00,
652 1.062699e+00 1718 1.062699e+00
653 1719
654 }; 1720 };
655 1721
656 1722
657 // 1723 //
658 // 1724 //
659 ////////////////////////////////////////////// 1725 ///////////////////////////////////////////////////////////////////////////////////////
660 1726