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