| 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 // 24.11.08 V. Grichine - first implementation 26 // 24.11.08 V. Grichine - first implementation
27 // 27 //
28 // 04.09.18 V. Ivantchenko Major revision of i <<
29 // 27.05.19 V. Ivantchenko Removed obsolete me <<
30 28
31 #include "G4ComponentGGNuclNuclXsc.hh" 29 #include "G4ComponentGGNuclNuclXsc.hh"
32 30
33 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
35 #include "G4NucleiProperties.hh" << 33 #include "G4ParticleTable.hh"
>> 34 #include "G4IonTable.hh"
36 #include "G4ParticleDefinition.hh" 35 #include "G4ParticleDefinition.hh"
>> 36 #include "G4HadTmpUtil.hh"
37 #include "G4HadronNucleonXsc.hh" 37 #include "G4HadronNucleonXsc.hh"
38 #include "G4ComponentGGHadronNucleusXsc.hh" <<
39 #include "G4NuclearRadii.hh" <<
40 #include "G4Pow.hh" <<
41 38
42 namespace <<
43 { <<
44 const G4double inve = 1./CLHEP::eplus; <<
45 } <<
46 39
47 G4ComponentGGNuclNuclXsc::G4ComponentGGNuclNuc 40 G4ComponentGGNuclNuclXsc::G4ComponentGGNuclNuclXsc()
48 : G4VComponentCrossSection("Glauber-Gribov Nu << 41 : G4VComponentCrossSection("Glauber-Gribov nucleus nucleus"),
>> 42 // fUpperLimit(100000*GeV),
>> 43 fLowerLimit(0.1*MeV),
>> 44 fRadiusConst(1.08*fermi), // 1.1, 1.3 ?
>> 45 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
>> 46 fDiffractionXsc(0.0),
>> 47 cacheDP(G4Proton::Proton(),G4ParticleMomentum(1.,0,0),0),
>> 48 dProton(G4Proton::Proton(),G4ParticleMomentum(1.,0,0),0),
>> 49 dNeutron(G4Neutron::Neutron(),G4ParticleMomentum(1.,0,0),0)
>> 50 // , fHadronNucleonXsc(0.0)
49 { 51 {
50 theProton = G4Proton::Proton(); 52 theProton = G4Proton::Proton();
51 theNeutron = G4Neutron::Neutron(); 53 theNeutron = G4Neutron::Neutron();
52 theLambda = G4Lambda::Lambda(); << 54 hnXsc = new G4HadronNucleonXsc();
53 fHNXsc = new G4HadronNucleonXsc(); <<
54 fHadrNucl = new G4ComponentGGHadronNucleusXs <<
55 } 55 }
56 56
>> 57
57 G4ComponentGGNuclNuclXsc::~G4ComponentGGNuclNu 58 G4ComponentGGNuclNuclXsc::~G4ComponentGGNuclNuclXsc()
58 { 59 {
59 delete fHNXsc; << 60 delete hnXsc;
60 } 61 }
61 62
62 ////////////////////////////////////////////// << 63 ////////////////////////////////////////////////////////////////////
63 64
64 G4double G4ComponentGGNuclNuclXsc::GetTotalEle << 65 G4double G4ComponentGGNuclNuclXsc::GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
65 const G4ParticleDefinition* aParticle << 66 G4double kinEnergy,
66 G4int Z, G4double A) << 67 G4int Z, G4int A)
67 { << 68 {
68 ComputeCrossSections(aParticle, kinEnergy, Z << 69 cacheDP.SetDefinition(aParticle);
>> 70 cacheDP.SetKineticEnergy(kinEnergy);
>> 71 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
69 return fTotalXsc; 72 return fTotalXsc;
70 } 73 }
71 74
72 ////////////////////////////////////////////// << 75 //////////////////////////////////////////////////////////////////////
73 76
74 G4double G4ComponentGGNuclNuclXsc::GetTotalIso << 77 G4double G4ComponentGGNuclNuclXsc::GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
75 const G4ParticleDefinition* aParticle << 78 G4double kinEnergy,
76 G4int Z, G4int A) << 79 G4int Z, G4double A)
77 { << 80 {
78 ComputeCrossSections(aParticle, kinEnergy, Z << 81 cacheDP.SetDefinition(aParticle);
>> 82 cacheDP.SetKineticEnergy(kinEnergy);
>> 83 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
79 return fTotalXsc; 84 return fTotalXsc;
80 } 85 }
81 86
82 ////////////////////////////////////////////// << 87 ////////////////////////////////////////////////////////////////////
83 88
84 G4double G4ComponentGGNuclNuclXsc::GetInelasti << 89 G4double G4ComponentGGNuclNuclXsc::GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
85 const G4ParticleDefinition* aParticle << 90 G4double kinEnergy,
86 G4int Z, G4double A) << 91 G4int Z, G4int A)
87 { << 92 {
88 ComputeCrossSections(aParticle, kinEnergy, Z << 93 cacheDP.SetDefinition(aParticle);
>> 94 cacheDP.SetKineticEnergy(kinEnergy);
>> 95 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
89 return fInelasticXsc; 96 return fInelasticXsc;
90 } 97 }
91 98
92 ////////////////////////////////////////////// << 99 /////////////////////////////////////////////////////////////////////
93 100
94 G4double G4ComponentGGNuclNuclXsc::GetInelasti << 101 G4double G4ComponentGGNuclNuclXsc::GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
95 const G4ParticleDefinition* aParticle << 102 G4double kinEnergy,
96 G4int Z, G4int A) << 103 G4int Z, G4double A)
97 { << 104 {
98 ComputeCrossSections(aParticle, kinEnergy, Z << 105 cacheDP.SetDefinition(aParticle);
>> 106 cacheDP.SetKineticEnergy(kinEnergy);
>> 107 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
99 return fInelasticXsc; 108 return fInelasticXsc;
100 } 109 }
101 110
102 ////////////////////////////////////////////// 111 //////////////////////////////////////////////////////////////////
103 112
104 G4double G4ComponentGGNuclNuclXsc::GetElasticE << 113 G4double G4ComponentGGNuclNuclXsc::GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
105 const G4ParticleDefinition* aParticle << 114 G4double kinEnergy,
106 G4int Z, G4double A) << 115 G4int Z, G4double A)
107 { << 116 {
108 ComputeCrossSections(aParticle, kinEnergy, Z << 117 cacheDP.SetDefinition(aParticle);
>> 118 cacheDP.SetKineticEnergy(kinEnergy);
>> 119 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
109 return fElasticXsc; 120 return fElasticXsc;
110 } 121 }
111 122
112 ////////////////////////////////////////////// 123 ///////////////////////////////////////////////////////////////////
113 124
114 G4double G4ComponentGGNuclNuclXsc::GetElasticI << 125 G4double G4ComponentGGNuclNuclXsc::GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
115 const G4ParticleDefinition* aParticle << 126 G4double kinEnergy,
116 G4int Z, G4int A) << 127 G4int Z, G4int A)
117 { << 128 {
118 ComputeCrossSections(aParticle, kinEnergy, Z << 129 cacheDP.SetDefinition(aParticle);
>> 130 cacheDP.SetKineticEnergy(kinEnergy);
>> 131 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
119 return fElasticXsc; 132 return fElasticXsc;
120 } 133 }
121 134
122 ////////////////////////////////////////////// 135 ////////////////////////////////////////////////////////////////
123 136
124 G4double G4ComponentGGNuclNuclXsc::ComputeQuas << 137 G4double G4ComponentGGNuclNuclXsc::ComputeQuasiElasticRatio(const G4ParticleDefinition* aParticle,
125 const G4ParticleDefinition* aParticle << 138 G4double kinEnergy,
126 G4int Z, G4int A) << 139 G4int Z, G4int A)
127 { << 140 {
128 ComputeCrossSections(aParticle, kinEnergy, Z << 141 cacheDP.SetDefinition(aParticle);
129 return (fInelasticXsc > fProductionXsc) << 142 cacheDP.SetKineticEnergy(kinEnergy);
130 ? (fInelasticXsc - fProductionXsc)/fInelas << 143 fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
131 } << 144 G4double ratio = 0.;
132 145
133 ////////////////////////////////////////////// << 146 if(fInelasticXsc > 0.)
>> 147 {
>> 148 ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
>> 149 if(ratio < 0.) ratio = 0.;
>> 150 }
>> 151 return ratio;
>> 152 }
134 153
135 void G4ComponentGGNuclNuclXsc::BuildPhysicsTab <<
136 {} <<
137 <<
138 ////////////////////////////////////////////// 154 //////////////////////////////////////////////////////////////////////
139 155
140 void G4ComponentGGNuclNuclXsc::DumpPhysicsTabl << 156 void
>> 157 G4ComponentGGNuclNuclXsc::CrossSectionDescription(std::ostream& outFile) const
141 { 158 {
142 G4cout << "G4ComponentGGNuclNuclXsc: uses Gl << 159 outFile << "G4ComponentGGNuclNuclXsc calculates total, inelastic and\n"
>> 160 << "elastic cross sections for nucleus-nucleus collisions using\n"
>> 161 << "the Glauber model with Gribov corrections. It is valid for\n"
>> 162 << "all incident energies above 100 keV./n";
143 } 163 }
144 164
145 ////////////////////////////////////////////// << 165 /////////////////////////////////////////////////////////////////////
146 166
147 void G4ComponentGGNuclNuclXsc::Description(std << 167 G4bool
>> 168 G4ComponentGGNuclNuclXsc::IsElementApplicable(const G4DynamicParticle* aDP,
>> 169 G4int Z, const G4Material*)
148 { 170 {
149 outFile << "G4ComponentGGNuclNuclXsc calcula << 171 G4bool applicable = false;
150 << "elastic cross sections for nucle << 172 G4double kineticEnergy = aDP->GetKineticEnergy();
151 << "the Glauber model with Gribov co << 173
152 << "all incident energies above 100 << 174 if (kineticEnergy >= fLowerLimit && Z > 1) applicable = true;
153 << "For the hydrogen target G4HadronNucleo << 175 return applicable;
>> 176 }
>> 177
>> 178 ///////////////////////////////////////////////////////////////////////////////
>> 179 //
>> 180 // Calculates total and inelastic Xsc, derives elastic as total - inelastic
>> 181 // accordong to Glauber model with Gribov correction calculated in the dipole
>> 182 // approximation on light cone. Gaussian density helps to calculate rest
>> 183 // integrals of the model. [1] B.Z. Kopeliovich, nucl-th/0306044
>> 184
>> 185
>> 186 G4double G4ComponentGGNuclNuclXsc::
>> 187 GetElementCrossSection(const G4DynamicParticle* aParticle, G4int Z,
>> 188 const G4Material*)
>> 189 {
>> 190 G4int A = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
>> 191 return GetZandACrossSection(aParticle, Z, A);
154 } 192 }
155 193
156 ////////////////////////////////////////////// 194 ///////////////////////////////////////////////////////////////////////////////
157 // 195 //
158 // Calculates total and inelastic Xsc, derives 196 // Calculates total and inelastic Xsc, derives elastic as total - inelastic
159 // accordong to Glauber model with Gribov corr 197 // accordong to Glauber model with Gribov correction calculated in the dipole
160 // approximation on light cone. Gaussian densi 198 // approximation on light cone. Gaussian density of point-like nucleons helps
161 // to calculate rest integrals of the model. [ 199 // to calculate rest integrals of the model. [1] B.Z. Kopeliovich,
162 // nucl-th/0306044 + simplification above 200 // nucl-th/0306044 + simplification above
163 201
164 void G4ComponentGGNuclNuclXsc::ComputeCrossSec <<
165 const G4ParticleDefinition* aParticle, G4 <<
166 G4int Z, G4int A) <<
167 { <<
168 // check cache <<
169 if(aParticle == fParticle && fZ == Z && fA = <<
170 { return; } <<
171 fParticle = aParticle; <<
172 fZ = Z; <<
173 fA = A; <<
174 fEnergy = kinEnergy; <<
175 G4Pow* pG4Pow = G4Pow::GetInstance(); <<
176 <<
177 G4int pZ = G4lrint(aParticle->GetPDGCharge() <<
178 G4int pA = aParticle->GetBaryonNumber(); <<
179 G4int pL = aParticle->GetNumberOfLambdasInHy <<
180 G4bool pHN = aParticle->IsHypernucleus(); <<
181 G4double cHN(0.88); <<
182 <<
183 // hydrogen <<
184 if(1 == Z && 1 == A) { <<
185 G4double e = kinEnergy*CLHEP::proton_mass_ <<
186 fHadrNucl->ComputeCrossSections( theProton <<
187 fTotalXsc = fHadrNucl->GetTotalGlauberGrib <<
188 fElasticXsc = fHadrNucl->GetElasticGlauber <<
189 fInelasticXsc = fHadrNucl->GetInelasticGla <<
190 fProductionXsc = fHadrNucl->GetProductionG <<
191 fDiffractionXsc = fHadrNucl->GetDiffractio <<
192 return; <<
193 } <<
194 static const G4double cofInelastic = 2.4; <<
195 static const G4double cofTotal = 2.0; <<
196 202
197 G4double pTkin = kinEnergy/(G4double)pA; << 203 G4double G4ComponentGGNuclNuclXsc::
>> 204 GetZandACrossSection(const G4DynamicParticle* aParticle,
>> 205 G4int tZ, G4int tA)
>> 206 {
>> 207 G4double xsection;
>> 208 G4double sigma;
>> 209 G4double cofInelastic = 2.4;
>> 210 G4double cofTotal = 2.0;
>> 211 G4double nucleusSquare;
>> 212 G4double cB;
>> 213 G4double ratio;
>> 214
>> 215 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
>> 216 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
>> 217
>> 218 G4double pTkin = aParticle->GetKineticEnergy();
>> 219 pTkin /= pA;
>> 220
>> 221 G4double pN = pA - pZ;
>> 222 if( pN < 0. ) pN = 0.;
198 223
199 G4int pN = pA - pZ; << 224 G4double tN = tA - tZ;
200 G4int tN = A - Z; << 225 if( tN < 0. ) tN = 0.;
201 226
202 G4double tR = G4NuclearRadii::Radius(Z, A); << 227 G4double tR = GetNucleusRadius( G4double(tZ),G4double(tA) );
203 G4double pR = G4NuclearRadii::Radius(pZ, pA) << 228 G4double pR = GetNucleusRadius(pZ,pA);
204 << 229
205 if(pHN) << 230 cB = GetCoulombBarier(aParticle, G4double(tZ), G4double(tA), pR, tR);
206 pR *= std::sqrt( pG4Pow->Z23( pA - pL ) + <<
207 <<
208 G4double cB = ComputeCoulombBarier(aParticle <<
209 231
210 if ( cB > 0. ) 232 if ( cB > 0. )
211 { 233 {
212 G4double sigma = (pZ*Z+pN*tN)*fHNXsc->Hadr << 234 dProton.SetKineticEnergy(pTkin);
213 if(pHN) sigma += pL*A*fHNXsc->HadronNucleo << 235 dNeutron.SetKineticEnergy(pTkin);
214 G4double ppInXsc = fHNXsc->GetInelasticHad << 236
215 << 237 sigma = (pZ*tZ+pN*tN)*hnXsc->GetHadronNucleonXscNS(&dProton, theProton);
216 sigma += (pZ*tN+pN*Z)*fHNXsc->HadronNucleo << 238
217 G4double npInXsc = fHNXsc->GetInelasticHad << 239 G4double ppInXsc = hnXsc->GetInelasticHadronNucleonXsc();
218 << 240
219 G4double nucleusSquare = cofTotal*CLHEP::p << 241 sigma += (pZ*tN+pN*tZ)*hnXsc->GetHadronNucleonXscNS(&dNeutron, theProton);
220 << 242
221 G4double ratio= sigma/nucleusSquare; << 243 G4double npInXsc = hnXsc->GetInelasticHadronNucleonXsc();
222 fTotalXsc = nucleusSquare*G4Log( 1. + << 244
223 fInelasticXsc = nucleusSquare*G4Log( 1. + << 245 // G4cout<<"ppInXsc = "<<ppInXsc/millibarn<<"; npInXsc = "<<npInXsc/millibarn<<G4endl;
224 fElasticXsc = std::max(fTotalXsc - fInel << 246 // G4cout<<"npTotXsc = "<<hnXsc->GetTotalHadronNucleonXsc()/millibarn<<"; npElXsc = "
>> 247 // <<hnXsc->GetElasticHadronNucleonXsc()/millibarn<<G4endl;
225 248
>> 249 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
>> 250
>> 251 ratio = sigma/nucleusSquare;
>> 252 xsection = nucleusSquare*std::log( 1. + ratio );
>> 253 fTotalXsc = xsection;
>> 254 fTotalXsc *= cB;
>> 255
>> 256 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
>> 257
>> 258 fInelasticXsc *= cB;
>> 259 fElasticXsc = fTotalXsc - fInelasticXsc;
>> 260
>> 261 // if (fElasticXsc < DBL_MIN) fElasticXsc = DBL_MIN;
>> 262 /*
226 G4double difratio = ratio/(1.+ratio); 263 G4double difratio = ratio/(1.+ratio);
227 fDiffractionXsc = 0.5*nucleusSquare*( difr <<
228 264
229 G4double xratio= ((pZ*Z+pN*tN)*ppInXsc + ( << 265 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
230 fProductionXsc = nucleusSquare*G4Log( 1. + << 266 */
231 fProductionXsc = std::min(fProductionXsc, << 267 // production to be checked !!! edit MK xsc
>> 268
>> 269 //sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscMK(theProton, pTkin, theProton) +
>> 270 // (pZ*tN+pN*tZ)*GetHadronNucleonXscMK(theProton, pTkin, theNeutron);
>> 271
>> 272 sigma = (pZ*tZ+pN*tN)*ppInXsc + (pZ*tN+pN*tZ)*npInXsc;
>> 273
>> 274 ratio = sigma/nucleusSquare;
>> 275 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
>> 276
>> 277 if (fElasticXsc < 0.) fElasticXsc = 0.;
232 } 278 }
233 else 279 else
234 { 280 {
235 fInelasticXsc = fTotalXsc = fElasticXsc = << 281 fInelasticXsc = 0.;
>> 282 fTotalXsc = 0.;
>> 283 fElasticXsc = 0.;
>> 284 fProductionXsc = 0.;
236 } 285 }
>> 286 return fInelasticXsc; // xsection;
237 } 287 }
238 288
239 ////////////////////////////////////////////// 289 ///////////////////////////////////////////////////////////////////////////////
>> 290 //
>> 291 //
240 292
241 G4double G4ComponentGGNuclNuclXsc::ComputeCoul << 293 G4double G4ComponentGGNuclNuclXsc::
242 const G4ParticleDefinitio << 294 GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA,
243 G4double pTkin, G4int Z, G4int A, << 295 G4double pR, G4double tR)
244 G4double pR, G4double tR) << 296 {
245 { << 297 G4double ratio;
246 G4int pZ = aParticle->GetPDGCharge()*inve; << 298 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
247 G4double pM = aParticle->GetPDGMass(); << 299
248 G4double tM = G4NucleiProperties::GetNuclear << 300 G4double pTkin = aParticle->GetKineticEnergy();
>> 301 // G4double pPlab = aParticle->GetTotalMomentum();
>> 302 G4double pM = aParticle->GetDefinition()->GetPDGMass();
>> 303 // G4double tM = tZ*proton_mass_c2 + (tA-tZ)*neutron_mass_c2; // ~ 1% accuracy
>> 304 G4double tM = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( G4int(tZ), G4int(tA) );
249 G4double pElab = pTkin + pM; 305 G4double pElab = pTkin + pM;
250 G4double totEcm = std::sqrt(pM*pM + tM*tM + << 306 G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
251 G4double totTcm = totEcm - pM -tM; << 307 // G4double pPcm = pPlab*tM/totEcm;
>> 308 // G4double pTcm = std::sqrt(pM*pM + pPcm*pPcm) - pM;
>> 309 G4double totTcm = totEcm - pM -tM;
252 310
253 // 0.5 defines shape of Cross section correc << 311 G4double bC = fine_structure_const*hbarc*pZ*tZ;
254 // at cB = totTcm it become zero << 312 bC /= pR + tR;
255 static const G4double qfact = 0.5*CLHEP::elm << 313 bC /= 2.; // 4., 2. parametrisation cof ??? vmg
256 G4double bC = qfact*pZ*Z/(pR + tR); << 314
257 << 315 // G4cout<<"pTkin = "<<pTkin/GeV<<"; pPlab = "
258 G4double ratio = (totTcm <= bC) ? 0. : 1. - << 316 // <<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "<<pTcm/GeV<<G4endl;
259 << 317
260 #ifdef G4VERBOSE << 318 if( totTcm <= bC ) ratio = 0.;
261 if (GetVerboseLevel() > 1) { << 319 else ratio = 1. - bC/totTcm;
262 G4cout << "G4ComponentGGNuclNuclXsc::Compu << 320
263 << "; pTkin(GeV)=" << pTkin/CLHEP::MeV << 321 // if(ratio < DBL_MIN) ratio = DBL_MIN;
264 << " totTcm= " << totTcm/CLHEP::MeV<< "; << 322 if( ratio < 0.) ratio = 0.;
265 << G4endl; << 323
266 } << 324 // G4cout <<"ratio = "<<ratio<<G4endl;
267 #endif <<
268 return ratio; 325 return ratio;
269 } 326 }
270 327
>> 328
271 ////////////////////////////////////////////// 329 //////////////////////////////////////////////////////////////////////////
272 // 330 //
273 // Return single-diffraction/inelastic cross-s 331 // Return single-diffraction/inelastic cross-section ratio
274 332
275 G4double G4ComponentGGNuclNuclXsc::GetRatioSD( << 333 G4double G4ComponentGGNuclNuclXsc::
276 const G4DynamicParticle* aParticle, G << 334 GetRatioSD(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
277 { 335 {
278 ComputeCrossSections(aParticle->GetDefinitio << 336 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
279 aParticle->GetKineticEn << 337
280 G4lrint(tZ), G4lrint(tA)); << 338 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
>> 339 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
>> 340
>> 341 G4double pTkin = aParticle->GetKineticEnergy();
>> 342 pTkin /= pA;
>> 343
>> 344 G4double pN = pA - pZ;
>> 345 if( pN < 0. ) pN = 0.;
>> 346
>> 347 G4double tN = tA - tZ;
>> 348 if( tN < 0. ) tN = 0.;
>> 349
>> 350 G4double tR = GetNucleusRadius(tZ,tA);
>> 351 G4double pR = GetNucleusRadius(pZ,pA);
281 352
282 return (fInelasticXsc > 0.0) ? fDiffractionX << 353 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
>> 354 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
>> 355
>> 356 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
>> 357 ratio = sigma/nucleusSquare;
>> 358 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
>> 359 G4double difratio = ratio/(1.+ratio);
>> 360
>> 361 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
>> 362
>> 363 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
>> 364 else ratio = 0.;
>> 365
>> 366 return ratio;
283 } 367 }
284 368
285 ////////////////////////////////////////////// 369 //////////////////////////////////////////////////////////////////////////
286 // 370 //
287 // Return quasi-elastic/inelastic cross-sectio 371 // Return quasi-elastic/inelastic cross-section ratio
288 372
289 G4double G4ComponentGGNuclNuclXsc::GetRatioQE( << 373 G4double G4ComponentGGNuclNuclXsc::
290 const G4DynamicParticle* aParticle, G << 374 GetRatioQE(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
>> 375 {
>> 376 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
>> 377
>> 378 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
>> 379 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
>> 380
>> 381 G4double pTkin = aParticle->GetKineticEnergy();
>> 382 pTkin /= pA;
>> 383
>> 384 G4double pN = pA - pZ;
>> 385 if( pN < 0. ) pN = 0.;
>> 386
>> 387 G4double tN = tA - tZ;
>> 388 if( tN < 0. ) tN = 0.;
>> 389
>> 390 G4double tR = GetNucleusRadius(tZ,tA);
>> 391 G4double pR = GetNucleusRadius(pZ,pA);
>> 392
>> 393 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
>> 394 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
>> 395
>> 396 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
>> 397 ratio = sigma/nucleusSquare;
>> 398 fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
>> 399
>> 400 // sigma = GetHNinelasticXsc(aParticle, tA, tZ);
>> 401 ratio = sigma/nucleusSquare;
>> 402 fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
>> 403
>> 404 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
>> 405 else ratio = 0.;
>> 406 if ( ratio < 0. ) ratio = 0.;
>> 407
>> 408 return ratio;
>> 409 }
>> 410
>> 411 ///////////////////////////////////////////////////////////////////////////////
>> 412 //
>> 413 // Returns hadron-nucleon Xsc according to differnt parametrisations:
>> 414 // [2] E. Levin, hep-ph/9710546
>> 415 // [3] U. Dersch, et al, hep-ex/9910052
>> 416 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
>> 417
>> 418 G4double
>> 419 G4ComponentGGNuclNuclXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
>> 420 const G4Element* anElement)
>> 421 {
>> 422 G4int At = G4lrint(anElement->GetN()); // number of nucleons
>> 423 G4int Zt = G4lrint(anElement->GetZ()); // number of protons
>> 424 return GetHadronNucleonXsc(aParticle, At, Zt);
>> 425 }
>> 426
>> 427
>> 428
>> 429
>> 430 ///////////////////////////////////////////////////////////////////////////////
>> 431 //
>> 432 // Returns hadron-nucleon Xsc according to differnt parametrisations:
>> 433 // [2] E. Levin, hep-ph/9710546
>> 434 // [3] U. Dersch, et al, hep-ex/9910052
>> 435 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
>> 436
>> 437 G4double
>> 438 G4ComponentGGNuclNuclXsc::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
>> 439 G4int At, G4int Zt)
>> 440 {
>> 441 G4double xsection = 0.;
>> 442
>> 443 G4double targ_mass = G4ParticleTable::GetParticleTable()->
>> 444 GetIonTable()->GetIonMass(Zt, At);
>> 445 targ_mass = 0.939*GeV; // ~mean neutron and proton ???
>> 446
>> 447 G4double proj_mass = aParticle->GetMass();
>> 448 G4double proj_momentum = aParticle->GetMomentum().mag();
>> 449 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 450
>> 451 sMand /= GeV*GeV; // in GeV for parametrisation
>> 452 proj_momentum /= GeV;
>> 453 const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
>> 454
>> 455 if(pParticle == theNeutron) // as proton ???
>> 456 {
>> 457 xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
>> 458 }
>> 459 else if(pParticle == theProton)
>> 460 {
>> 461 xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
>> 462 }
>> 463
>> 464 xsection *= millibarn;
>> 465 return xsection;
>> 466 }
>> 467
>> 468
>> 469
>> 470 ///////////////////////////////////////////////////////////////////////////////
>> 471 //
>> 472 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
>> 473 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
>> 474 // At = number of nucleons, Zt = number of protons
>> 475
>> 476 G4double
>> 477 G4ComponentGGNuclNuclXsc::GetHadronNucleonXscPDG(const G4ParticleDefinition* pParticle,
>> 478 G4double sMand,
>> 479 const G4ParticleDefinition* tParticle)
>> 480 {
>> 481 G4double xsection = 0.;
>> 482 // G4bool pORn = (tParticle == theProton || nucleon == theNeutron );
>> 483 G4bool proton = (tParticle == theProton);
>> 484 G4bool neutron = (tParticle == theNeutron);
>> 485
>> 486 // General PDG fit constants
>> 487
>> 488 G4double s0 = 5.38*5.38; // in Gev^2
>> 489 G4double eta1 = 0.458;
>> 490 G4double eta2 = 0.458;
>> 491 G4double B = 0.308;
>> 492
>> 493 // const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
>> 494
>> 495 if(pParticle == theNeutron) // proton-neutron fit
>> 496 {
>> 497 if ( proton )
>> 498 {
>> 499 xsection = ( 35.80 + B*std::pow(std::log(sMand/s0),2.)
>> 500 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
>> 501 }
>> 502 if ( neutron )
>> 503 {
>> 504 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
>> 505 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
>> 506 }
>> 507 }
>> 508 else if(pParticle == theProton)
>> 509 {
>> 510 if ( proton )
>> 511 {
>> 512 xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
>> 513 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
>> 514
>> 515 }
>> 516 if ( neutron )
>> 517 {
>> 518 xsection = (35.80 + B*std::pow(std::log(sMand/s0),2.)
>> 519 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
>> 520 }
>> 521 }
>> 522 xsection *= millibarn; // parametrised in mb
>> 523 return xsection;
>> 524 }
>> 525
>> 526
>> 527 ///////////////////////////////////////////////////////////////////////////////
>> 528 //
>> 529 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
>> 530 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
>> 531 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
>> 532
>> 533 G4double
>> 534 G4ComponentGGNuclNuclXsc::GetHadronNucleonXscNS(const G4ParticleDefinition* pParticle,
>> 535 G4double pTkin,
>> 536 const G4ParticleDefinition* tParticle)
>> 537 {
>> 538 G4double xsection(0);
>> 539 // G4double Delta; DHW 19 May 2011: variable set but not used
>> 540 G4double A0, B0;
>> 541 G4double hpXscv(0);
>> 542 G4double hnXscv(0);
>> 543
>> 544 G4double targ_mass = tParticle->GetPDGMass();
>> 545 G4double proj_mass = pParticle->GetPDGMass();
>> 546
>> 547 G4double proj_energy = proj_mass + pTkin;
>> 548 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
>> 549
>> 550 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
>> 551
>> 552 sMand /= GeV*GeV; // in GeV for parametrisation
>> 553 proj_momentum /= GeV;
>> 554 proj_energy /= GeV;
>> 555 proj_mass /= GeV;
>> 556
>> 557 // General PDG fit constants
>> 558
>> 559 // G4double s0 = 5.38*5.38; // in Gev^2
>> 560 // G4double eta1 = 0.458;
>> 561 // G4double eta2 = 0.458;
>> 562 // G4double B = 0.308;
>> 563
>> 564 if( proj_momentum >= 373.)
>> 565 {
>> 566 return GetHadronNucleonXscPDG(pParticle,sMand,tParticle);
>> 567 }
>> 568 else if( proj_momentum >= 10. ) // high energy: pp = nn = np
>> 569 // if( proj_momentum >= 2.)
>> 570 {
>> 571 // Delta = 1.; // DHW 19 May 2011: variable set but not used
>> 572 // if (proj_energy < 40.) Delta = 0.916+0.0021*proj_energy;
>> 573
>> 574 if (proj_momentum >= 10.) {
>> 575 B0 = 7.5;
>> 576 A0 = 100. - B0*std::log(3.0e7);
>> 577
>> 578 xsection = A0 + B0*std::log(proj_energy) - 11
>> 579 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
>> 580 0.93827*0.93827,-0.165); // mb
>> 581 }
>> 582 }
>> 583 else // low energy pp = nn != np
>> 584 {
>> 585 if(pParticle == tParticle) // pp or nn // nn to be pp
>> 586 {
>> 587 if( proj_momentum < 0.73 )
>> 588 {
>> 589 hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
>> 590 }
>> 591 else if( proj_momentum < 1.05 )
>> 592 {
>> 593 hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
>> 594 (std::log(proj_momentum/0.73));
>> 595 }
>> 596 else // if( proj_momentum < 10. )
>> 597 {
>> 598 hnXscv = 39.0 +
>> 599 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
>> 600 }
>> 601 xsection = hnXscv;
>> 602 }
>> 603 else // pn to be np
>> 604 {
>> 605 if( proj_momentum < 0.8 )
>> 606 {
>> 607 hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
>> 608 }
>> 609 else if( proj_momentum < 1.4 )
>> 610 {
>> 611 hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
>> 612 }
>> 613 else // if( proj_momentum < 10. )
>> 614 {
>> 615 hpXscv = 33.3+
>> 616 20.8*(std::pow(proj_momentum,2.0)-1.35)/
>> 617 (std::pow(proj_momentum,2.50)+0.95);
>> 618 }
>> 619 xsection = hpXscv;
>> 620 }
>> 621 }
>> 622 xsection *= millibarn; // parametrised in mb
>> 623 return xsection;
>> 624 }
>> 625
>> 626 /////////////////////////////////////////////////////////////////////////////////
>> 627 //
>> 628 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
>> 629
>> 630 G4double
>> 631 G4ComponentGGNuclNuclXsc::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
>> 632 G4int At, G4int Zt)
>> 633 {
>> 634 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
>> 635 G4int absPDGcode = std::abs(PDGcode);
>> 636 G4double Elab = aParticle->GetTotalEnergy();
>> 637 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
>> 638 G4double Plab = aParticle->GetMomentum().mag();
>> 639 // std::sqrt(Elab * Elab - 0.88);
>> 640
>> 641 Elab /= GeV;
>> 642 Plab /= GeV;
>> 643
>> 644 G4double LogPlab = std::log( Plab );
>> 645 G4double sqrLogPlab = LogPlab * LogPlab;
>> 646
>> 647 //G4cout<<"Plab = "<<Plab<<G4endl;
>> 648
>> 649 G4double NumberOfTargetProtons = Zt;
>> 650 G4double NumberOfTargetNucleons = At;
>> 651 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
>> 652
>> 653 if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
>> 654
>> 655 G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
>> 656
>> 657 if( absPDGcode > 1000 ) //------Projectile is baryon --------
>> 658 {
>> 659 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
>> 660 0.522*sqrLogPlab - 4.51*LogPlab;
>> 661
>> 662 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
>> 663 0.513*sqrLogPlab - 4.27*LogPlab;
>> 664
>> 665 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
>> 666 0.169*sqrLogPlab - 1.85*LogPlab;
>> 667
>> 668 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
>> 669 0.169*sqrLogPlab - 1.85*LogPlab;
>> 670
>> 671 Xtotal = ( NumberOfTargetProtons * XtotPP +
>> 672 NumberOfTargetNeutrons * XtotPN );
>> 673
>> 674 Xelastic = ( NumberOfTargetProtons * XelPP +
>> 675 NumberOfTargetNeutrons * XelPN );
>> 676 }
>> 677
>> 678 Xinelastic = Xtotal - Xelastic;
>> 679 if(Xinelastic < 0.) Xinelastic = 0.;
>> 680
>> 681 return Xinelastic*= millibarn;
>> 682 }
>> 683
>> 684 ///////////////////////////////////////////////////////////////////////////////
>> 685 //
>> 686 //
>> 687
>> 688 G4double
>> 689 G4ComponentGGNuclNuclXsc::GetNucleusRadius(const G4DynamicParticle* ,
>> 690 const G4Element* anElement)
>> 691 {
>> 692 G4double At = anElement->GetN();
>> 693 G4double oneThird = 1.0/3.0;
>> 694 G4double cubicrAt = std::pow (At, oneThird);
>> 695
>> 696 G4double R; // = fRadiusConst*cubicrAt;
>> 697 R = fRadiusConst*cubicrAt;
>> 698
>> 699 G4double meanA = 21.;
>> 700 G4double tauA1 = 40.;
>> 701 G4double tauA2 = 10.;
>> 702 G4double tauA3 = 5.;
>> 703
>> 704 G4double a1 = 0.85;
>> 705 G4double b1 = 1. - a1;
>> 706
>> 707 G4double b2 = 0.3;
>> 708 G4double b3 = 4.;
>> 709
>> 710 if (At > 20.) // 20.
>> 711 {
>> 712 R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
>> 713 }
>> 714 else if (At > 3.5)
>> 715 {
>> 716 R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
>> 717 }
>> 718 else
>> 719 {
>> 720 R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
>> 721 }
>> 722
>> 723 return R;
>> 724 }
>> 725
>> 726 ///////////////////////////////////////////////////////////////////////////////
>> 727 //
>> 728 //
>> 729
>> 730 G4double
>> 731 G4ComponentGGNuclNuclXsc::GetNucleusRadius(G4double Zt, G4double At)
>> 732 {
>> 733 G4double R;
>> 734 R = GetNucleusRadiusDE(Zt,At);
>> 735 // R = GetNucleusRadiusRMS(Zt,At);
>> 736
>> 737 return R;
>> 738 }
>> 739
>> 740 ///////////////////////////////////////////////////////////////////
>> 741
>> 742 G4double
>> 743 G4ComponentGGNuclNuclXsc::GetNucleusRadiusGG(G4double At)
>> 744 {
>> 745 G4double oneThird = 1.0/3.0;
>> 746 G4double cubicrAt = std::pow (At, oneThird);
>> 747
>> 748 G4double R; // = fRadiusConst*cubicrAt;
>> 749 R = fRadiusConst*cubicrAt;
>> 750
>> 751 G4double meanA = 20.;
>> 752 G4double tauA = 20.;
>> 753
>> 754 if ( At > 20.) // 20.
>> 755 {
>> 756 R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
>> 757 }
>> 758 else
>> 759 {
>> 760 R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
>> 761 }
>> 762
>> 763 return R;
>> 764 }
>> 765
>> 766 /////////////////////////////////////////////////////////////////////////////
>> 767 //
>> 768 //
>> 769
>> 770 G4double
>> 771 G4ComponentGGNuclNuclXsc::GetNucleusRadiusDE(G4double Z, G4double A)
>> 772 {
>> 773 // algorithm from diffuse-elastic
>> 774
>> 775 G4double R, r0, a11, a12, a13, a2, a3;
>> 776
>> 777 a11 = 1.26; // 1.08, 1.16
>> 778 a12 = 1.; // 1.08, 1.16
>> 779 a13 = 1.12; // 1.08, 1.16
>> 780 a2 = 1.1;
>> 781 a3 = 1.;
>> 782
>> 783 // Special rms radii for light nucleii
>> 784
>> 785 if (A < 50.)
>> 786 {
>> 787 if (std::abs(A-1.) < 0.5) return 0.89*fermi; // p
>> 788 else if(std::abs(A-2.) < 0.5) return 2.13*fermi; // d
>> 789 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi; // t
>> 790
>> 791 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi; // He3
>> 792 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi; // He4
>> 793
>> 794 else if(std::abs(Z-3.) < 0.5) return 2.40*fermi; // Li7
>> 795 else if(std::abs(Z-4.) < 0.5) return 2.51*fermi; // Be9
>> 796
>> 797 else if( 10. < A && A <= 16. ) r0 = a11*( 1 - std::pow(A, -2./3.) )*fermi; // 1.08*fermi;
>> 798 else if( 15. < A && A <= 20. ) r0 = a12*( 1 - std::pow(A, -2./3.) )*fermi;
>> 799 else if( 20. < A && A <= 30. ) r0 = a13*( 1 - std::pow(A, -2./3.) )*fermi;
>> 800 else r0 = a2*fermi;
>> 801
>> 802 R = r0*std::pow( A, 1./3. );
>> 803 }
>> 804 else
>> 805 {
>> 806 r0 = a3*fermi;
>> 807
>> 808 R = r0*std::pow(A, 0.27);
>> 809 }
>> 810 return R;
>> 811 }
>> 812
>> 813
>> 814 /////////////////////////////////////////////////////////////////////////////
>> 815 //
>> 816 // RMS radii from e-A scattering data
>> 817
>> 818 G4double
>> 819 G4ComponentGGNuclNuclXsc::GetNucleusRadiusRMS(G4double Z, G4double A)
291 { 820 {
292 ComputeCrossSections(aParticle->GetDefinitio <<
293 aParticle->GetKineticEn <<
294 G4lrint(tZ), G4lrint(tA)); <<
295 821
296 return (fInelasticXsc > 0.0) ? 1.0 - fProduc << 822 if (std::abs(A-1.) < 0.5) return 0.89*fermi; // p
>> 823 else if(std::abs(A-2.) < 0.5) return 2.13*fermi; // d
>> 824 else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi; // t
>> 825
>> 826 else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi; // He3
>> 827 else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi; // He4
>> 828
>> 829 else if(std::abs(Z-3.) < 0.5) return 2.40*fermi; // Li7
>> 830 else if(std::abs(Z-4.) < 0.5) return 2.51*fermi; // Be9
>> 831
>> 832 else return 1.24*std::pow(A, 0.28 )*fermi; // A > 9
297 } 833 }
298 834
>> 835
>> 836 ///////////////////////////////////////////////////////////////////////////////
>> 837 //
>> 838 //
>> 839
>> 840 G4double G4ComponentGGNuclNuclXsc::CalculateEcmValue(const G4double mp,
>> 841 const G4double mt,
>> 842 const G4double Plab)
>> 843 {
>> 844 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
>> 845 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
>> 846 // G4double Pcm = Plab * mt / Ecm;
>> 847 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
>> 848
>> 849 return Ecm ; // KEcm;
>> 850 }
>> 851
>> 852
>> 853 ///////////////////////////////////////////////////////////////////////////////
>> 854 //
>> 855 //
>> 856
>> 857 G4double G4ComponentGGNuclNuclXsc::CalcMandelstamS(const G4double mp,
>> 858 const G4double mt,
>> 859 const G4double Plab)
>> 860 {
>> 861 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
>> 862 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
>> 863
>> 864 return sMand;
>> 865 }
>> 866
>> 867 //
>> 868 //
299 ////////////////////////////////////////////// 869 ///////////////////////////////////////////////////////////////////////////////
300 870