| 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 // $Id: G4NuclNuclDiffuseElastic.cc 98826 2016-08-12 12:36:07Z gcosmo $
26 // 27 //
27 // 28 //
28 // Physics model class G4NuclNuclDiffuseElasti 29 // Physics model class G4NuclNuclDiffuseElastic
29 // 30 //
30 // 31 //
31 // G4 Model: optical diffuse elastic scatterin 32 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
32 // 33 //
33 // 24-May-07 V. Grichine 34 // 24-May-07 V. Grichine
34 // 35 //
35 36
36 #include "G4NuclNuclDiffuseElastic.hh" 37 #include "G4NuclNuclDiffuseElastic.hh"
37 #include "G4ParticleTable.hh" 38 #include "G4ParticleTable.hh"
38 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh"
39 #include "G4IonTable.hh" 40 #include "G4IonTable.hh"
40 #include "G4NucleiProperties.hh" 41 #include "G4NucleiProperties.hh"
41 42
42 #include "Randomize.hh" 43 #include "Randomize.hh"
43 #include "G4Integrator.hh" 44 #include "G4Integrator.hh"
44 #include "globals.hh" 45 #include "globals.hh"
45 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
47 48
48 #include "G4Proton.hh" 49 #include "G4Proton.hh"
49 #include "G4Neutron.hh" 50 #include "G4Neutron.hh"
50 #include "G4Deuteron.hh" 51 #include "G4Deuteron.hh"
51 #include "G4Alpha.hh" 52 #include "G4Alpha.hh"
52 #include "G4PionPlus.hh" 53 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh" 54 #include "G4PionMinus.hh"
54 55
55 #include "G4Element.hh" 56 #include "G4Element.hh"
56 #include "G4ElementTable.hh" 57 #include "G4ElementTable.hh"
57 #include "G4NistManager.hh" 58 #include "G4NistManager.hh"
58 #include "G4PhysicsTable.hh" 59 #include "G4PhysicsTable.hh"
59 #include "G4PhysicsLogVector.hh" 60 #include "G4PhysicsLogVector.hh"
60 #include "G4PhysicsFreeVector.hh" 61 #include "G4PhysicsFreeVector.hh"
61 #include "G4HadronicParameters.hh" <<
62 62
63 ////////////////////////////////////////////// 63 /////////////////////////////////////////////////////////////////////////
64 // 64 //
65 // Test Constructor. Just to check xsc 65 // Test Constructor. Just to check xsc
66 66
67 67
68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseEla 68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic()
69 : G4HadronElastic("NNDiffuseElastic"), fPart 69 : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
70 { 70 {
71 SetMinEnergy( 50*MeV ); 71 SetMinEnergy( 50*MeV );
72 SetMaxEnergy( G4HadronicParameters::Instance << 72 SetMaxEnergy( 1.*TeV );
73 verboseLevel = 0; 73 verboseLevel = 0;
74 lowEnergyRecoilLimit = 100.*keV; 74 lowEnergyRecoilLimit = 100.*keV;
75 lowEnergyLimitQ = 0.0*GeV; 75 lowEnergyLimitQ = 0.0*GeV;
76 lowEnergyLimitHE = 0.0*GeV; 76 lowEnergyLimitHE = 0.0*GeV;
77 lowestEnergyLimit= 0.0*keV; 77 lowestEnergyLimit= 0.0*keV;
78 plabLowLimit = 20.0*MeV; 78 plabLowLimit = 20.0*MeV;
79 79
80 theProton = G4Proton::Proton(); 80 theProton = G4Proton::Proton();
81 theNeutron = G4Neutron::Neutron(); 81 theNeutron = G4Neutron::Neutron();
82 theDeuteron = G4Deuteron::Deuteron(); 82 theDeuteron = G4Deuteron::Deuteron();
83 theAlpha = G4Alpha::Alpha(); 83 theAlpha = G4Alpha::Alpha();
84 thePionPlus = G4PionPlus::PionPlus(); 84 thePionPlus = G4PionPlus::PionPlus();
85 thePionMinus= G4PionMinus::PionMinus(); 85 thePionMinus= G4PionMinus::PionMinus();
86 86
87 fEnergyBin = 300; // Increased from the ori << 87 fEnergyBin = 200;
88 fAngleBin = 200; 88 fAngleBin = 200;
89 89
90 fEnergyVector = new G4PhysicsLogVector( the 90 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
91 fAngleTable = 0; 91 fAngleTable = 0;
92 92
93 fParticle = 0; 93 fParticle = 0;
94 fWaveVector = 0.; 94 fWaveVector = 0.;
95 fAtomicWeight = 0.; 95 fAtomicWeight = 0.;
96 fAtomicNumber = 0.; 96 fAtomicNumber = 0.;
97 fNuclearRadius = 0.; 97 fNuclearRadius = 0.;
98 fBeta = 0.; 98 fBeta = 0.;
99 fZommerfeld = 0.; 99 fZommerfeld = 0.;
100 fAm = 0.; 100 fAm = 0.;
101 fAddCoulomb = false; 101 fAddCoulomb = false;
102 // Ranges of angle table relative to current 102 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
103 103
104 // Empirical parameters 104 // Empirical parameters
105 105
106 fCofAlphaMax = 1.5; 106 fCofAlphaMax = 1.5;
107 fCofAlphaCoulomb = 0.5; 107 fCofAlphaCoulomb = 0.5;
108 108
109 fProfileDelta = 1.; 109 fProfileDelta = 1.;
110 fProfileAlpha = 0.5; 110 fProfileAlpha = 0.5;
111 111
112 fCofLambda = 1.0; 112 fCofLambda = 1.0;
113 fCofDelta = 0.04; 113 fCofDelta = 0.04;
114 fCofAlpha = 0.095; 114 fCofAlpha = 0.095;
115 115
116 fNuclearRadius1 = fNuclearRadius2 = fNuclear 116 fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare
117 = fRutherfordRatio = fCoulombPhase0 = fHal 117 = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2
118 = fRutherfordTheta = fProfileLambda = fCof << 118 = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar = fCofAlphaMax
119 = fSumSigma = fEtaRatio = fReZ = 0.0; << 119 = fCofAlphaCoulomb = fSumSigma = fEtaRatio = fReZ = 0.0;
120 fMaxL = 0; 120 fMaxL = 0;
121 121
122 fNuclearRadiusCof = 1.0; 122 fNuclearRadiusCof = 1.0;
123 fCoulombMuC = 0.0; <<
124 } 123 }
125 124
126 ////////////////////////////////////////////// 125 //////////////////////////////////////////////////////////////////////////////
127 // 126 //
128 // Destructor 127 // Destructor
129 128
130 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseEl 129 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic()
131 { 130 {
132 if ( fEnergyVector ) { 131 if ( fEnergyVector ) {
133 delete fEnergyVector; 132 delete fEnergyVector;
134 fEnergyVector = 0; 133 fEnergyVector = 0;
135 } 134 }
136 135
137 for ( std::vector<G4PhysicsTable*>::iterator 136 for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
138 it != fAngleBank.end(); ++it ) { 137 it != fAngleBank.end(); ++it ) {
139 if ( (*it) ) (*it)->clearAndDestroy(); 138 if ( (*it) ) (*it)->clearAndDestroy();
140 delete *it; 139 delete *it;
141 *it = 0; 140 *it = 0;
142 } 141 }
143 fAngleTable = 0; 142 fAngleTable = 0;
144 } 143 }
145 144
146 ////////////////////////////////////////////// 145 //////////////////////////////////////////////////////////////////////////////
147 // 146 //
148 // Initialisation for given particle using ele 147 // Initialisation for given particle using element table of application
149 148
150 void G4NuclNuclDiffuseElastic::Initialise() 149 void G4NuclNuclDiffuseElastic::Initialise()
151 { 150 {
152 151
153 // fEnergyVector = new G4PhysicsLogVector( t 152 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
154 153
155 const G4ElementTable* theElementTable = G4El 154 const G4ElementTable* theElementTable = G4Element::GetElementTable();
156 std::size_t jEl, numOfEl = G4Element::GetNum << 155 size_t jEl, numOfEl = G4Element::GetNumberOfElements();
157 156
158 // projectile radius 157 // projectile radius
159 158
160 G4double A1 = G4double( fParticle->GetBaryon 159 G4double A1 = G4double( fParticle->GetBaryonNumber() );
161 G4double R1 = CalculateNuclearRad(A1); 160 G4double R1 = CalculateNuclearRad(A1);
162 161
163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // appli 162 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
164 { 163 {
165 fAtomicNumber = (*theElementTable)[jEl]->G 164 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
166 fAtomicWeight = G4NistManager::Instance()- 165 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
167 166
168 fNuclearRadius = CalculateNuclearRad(fAtom 167 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
169 fNuclearRadius += R1; 168 fNuclearRadius += R1;
170 169
171 if(verboseLevel > 0) 170 if(verboseLevel > 0)
172 { 171 {
173 G4cout<<"G4NuclNuclDiffuseElastic::Initi 172 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
174 <<(*theElementTable)[jEl]->GetName()<<G4 173 <<(*theElementTable)[jEl]->GetName()<<G4endl;
175 } 174 }
176 fElementNumberVector.push_back(fAtomicNumb 175 fElementNumberVector.push_back(fAtomicNumber);
177 fElementNameVector.push_back((*theElementT 176 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
178 177
179 BuildAngleTable(); 178 BuildAngleTable();
180 fAngleBank.push_back(fAngleTable); 179 fAngleBank.push_back(fAngleTable);
181 } 180 }
182 } 181 }
183 182
184 183
185 ////////////////////////////////////////////// 184 ////////////////////////////////////////////////////////////////////////////
186 // 185 //
187 // return differential elastic cross section d 186 // return differential elastic cross section d(sigma)/d(omega)
188 187
189 G4double 188 G4double
190 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc 189 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
191 G4doub 190 G4double theta,
192 G4double momentum, 191 G4double momentum,
193 G4doub 192 G4double A )
194 { 193 {
195 fParticle = particle; 194 fParticle = particle;
196 fWaveVector = momentum/hbarc; 195 fWaveVector = momentum/hbarc;
197 fAtomicWeight = A; 196 fAtomicWeight = A;
198 fAddCoulomb = false; 197 fAddCoulomb = false;
199 fNuclearRadius = CalculateNuclearRad(A); 198 fNuclearRadius = CalculateNuclearRad(A);
200 199
201 G4double sigma = fNuclearRadius*fNuclearRadi 200 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
202 201
203 return sigma; 202 return sigma;
204 } 203 }
205 204
206 ////////////////////////////////////////////// 205 ////////////////////////////////////////////////////////////////////////////
207 // 206 //
208 // return invariant differential elastic cross 207 // return invariant differential elastic cross section d(sigma)/d(tMand)
209 208
210 G4double 209 G4double
211 G4NuclNuclDiffuseElastic::GetInvElasticXsc( co 210 G4NuclNuclDiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle,
212 G4doub 211 G4double tMand,
213 G4double plab, 212 G4double plab,
214 G4doub 213 G4double A, G4double Z )
215 { 214 {
216 G4double m1 = particle->GetPDGMass(); 215 G4double m1 = particle->GetPDGMass();
217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 216 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
218 217
219 G4int iZ = static_cast<G4int>(Z+0.5); 218 G4int iZ = static_cast<G4int>(Z+0.5);
220 G4int iA = static_cast<G4int>(A+0.5); 219 G4int iA = static_cast<G4int>(A+0.5);
221 G4ParticleDefinition * theDef = 0; 220 G4ParticleDefinition * theDef = 0;
222 221
223 if (iZ == 1 && iA == 1) theDef = thePro 222 if (iZ == 1 && iA == 1) theDef = theProton;
224 else if (iZ == 1 && iA == 2) theDef = theDeu 223 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
225 else if (iZ == 1 && iA == 3) theDef = G4Trit 224 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
226 else if (iZ == 2 && iA == 3) theDef = G4He3: 225 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
227 else if (iZ == 2 && iA == 4) theDef = theAlp 226 else if (iZ == 2 && iA == 4) theDef = theAlpha;
228 else theDef = G4ParticleTable::GetParticleTa 227 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
229 228
230 G4double tmass = theDef->GetPDGMass(); 229 G4double tmass = theDef->GetPDGMass();
231 230
232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 231 G4LorentzVector lv(0.0,0.0,0.0,tmass);
233 lv += lv1; 232 lv += lv1;
234 233
235 G4ThreeVector bst = lv.boostVector(); 234 G4ThreeVector bst = lv.boostVector();
236 lv1.boost(-bst); 235 lv1.boost(-bst);
237 236
238 G4ThreeVector p1 = lv1.vect(); 237 G4ThreeVector p1 = lv1.vect();
239 G4double ptot = p1.mag(); 238 G4double ptot = p1.mag();
240 G4double ptot2 = ptot*ptot; 239 G4double ptot2 = ptot*ptot;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 240 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
242 241
243 if( cost >= 1.0 ) cost = 1.0; 242 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0; 243 else if( cost <= -1.0) cost = -1.0;
245 244
246 G4double thetaCMS = std::acos(cost); 245 G4double thetaCMS = std::acos(cost);
247 246
248 G4double sigma = GetDiffuseElasticXsc( parti 247 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
249 248
250 sigma *= pi/ptot2; 249 sigma *= pi/ptot2;
251 250
252 return sigma; 251 return sigma;
253 } 252 }
254 253
255 ////////////////////////////////////////////// 254 ////////////////////////////////////////////////////////////////////////////
256 // 255 //
257 // return differential elastic cross section d 256 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
258 // correction 257 // correction
259 258
260 G4double 259 G4double
261 G4NuclNuclDiffuseElastic::GetDiffuseElasticSum 260 G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
262 G4doub 261 G4double theta,
263 G4double momentum, 262 G4double momentum,
264 G4doub 263 G4double A, G4double Z )
265 { 264 {
266 fParticle = particle; 265 fParticle = particle;
267 fWaveVector = momentum/hbarc; 266 fWaveVector = momentum/hbarc;
268 fAtomicWeight = A; 267 fAtomicWeight = A;
269 fAtomicNumber = Z; 268 fAtomicNumber = Z;
270 fNuclearRadius = CalculateNuclearRad(A); 269 fNuclearRadius = CalculateNuclearRad(A);
271 fAddCoulomb = false; 270 fAddCoulomb = false;
272 271
273 G4double z = particle->GetPDGCharge(); 272 G4double z = particle->GetPDGCharge();
274 273
275 G4double kRt = fWaveVector*fNuclearRadius* 274 G4double kRt = fWaveVector*fNuclearRadius*theta;
276 G4double kRtC = 1.9; 275 G4double kRtC = 1.9;
277 276
278 if( z && (kRt > kRtC) ) 277 if( z && (kRt > kRtC) )
279 { 278 {
280 fAddCoulomb = true; 279 fAddCoulomb = true;
281 fBeta = CalculateParticleBeta( parti 280 fBeta = CalculateParticleBeta( particle, momentum);
282 fZommerfeld = CalculateZommerfeld( fBeta, 281 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
283 fAm = CalculateAm( momentum, fZomm 282 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
284 } 283 }
285 G4double sigma = fNuclearRadius*fNuclearRadi 284 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
286 285
287 return sigma; 286 return sigma;
288 } 287 }
289 288
290 ////////////////////////////////////////////// 289 ////////////////////////////////////////////////////////////////////////////
291 // 290 //
292 // return invariant differential elastic cross 291 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
293 // correction 292 // correction
294 293
295 G4double 294 G4double
296 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( 295 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle,
297 G4doub 296 G4double tMand,
298 G4double plab, 297 G4double plab,
299 G4doub 298 G4double A, G4double Z )
300 { 299 {
301 G4double m1 = particle->GetPDGMass(); 300 G4double m1 = particle->GetPDGMass();
302 301
303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 302 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
304 303
305 G4int iZ = static_cast<G4int>(Z+0.5); 304 G4int iZ = static_cast<G4int>(Z+0.5);
306 G4int iA = static_cast<G4int>(A+0.5); 305 G4int iA = static_cast<G4int>(A+0.5);
307 306
308 G4ParticleDefinition* theDef = 0; 307 G4ParticleDefinition* theDef = 0;
309 308
310 if (iZ == 1 && iA == 1) theDef = thePro 309 if (iZ == 1 && iA == 1) theDef = theProton;
311 else if (iZ == 1 && iA == 2) theDef = theDeu 310 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
312 else if (iZ == 1 && iA == 3) theDef = G4Trit 311 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
313 else if (iZ == 2 && iA == 3) theDef = G4He3: 312 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
314 else if (iZ == 2 && iA == 4) theDef = theAlp 313 else if (iZ == 2 && iA == 4) theDef = theAlpha;
315 else theDef = G4ParticleTable::GetParticleTa 314 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
316 315
317 G4double tmass = theDef->GetPDGMass(); 316 G4double tmass = theDef->GetPDGMass();
318 317
319 G4LorentzVector lv(0.0,0.0,0.0,tmass); 318 G4LorentzVector lv(0.0,0.0,0.0,tmass);
320 lv += lv1; 319 lv += lv1;
321 320
322 G4ThreeVector bst = lv.boostVector(); 321 G4ThreeVector bst = lv.boostVector();
323 lv1.boost(-bst); 322 lv1.boost(-bst);
324 323
325 G4ThreeVector p1 = lv1.vect(); 324 G4ThreeVector p1 = lv1.vect();
326 G4double ptot = p1.mag(); 325 G4double ptot = p1.mag();
327 G4double ptot2 = ptot*ptot; 326 G4double ptot2 = ptot*ptot;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/ 327 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
329 328
330 if( cost >= 1.0 ) cost = 1.0; 329 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0; 330 else if( cost <= -1.0) cost = -1.0;
332 331
333 G4double thetaCMS = std::acos(cost); 332 G4double thetaCMS = std::acos(cost);
334 333
335 G4double sigma = GetDiffuseElasticSumXsc( pa 334 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
336 335
337 sigma *= pi/ptot2; 336 sigma *= pi/ptot2;
338 337
339 return sigma; 338 return sigma;
340 } 339 }
341 340
342 ////////////////////////////////////////////// 341 ////////////////////////////////////////////////////////////////////////////
343 // 342 //
344 // return invariant differential elastic cross 343 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
345 // correction 344 // correction
346 345
347 G4double 346 G4double
348 G4NuclNuclDiffuseElastic::GetInvCoulombElastic 347 G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
349 G4doub 348 G4double tMand,
350 G4double plab, 349 G4double plab,
351 G4doub 350 G4double A, G4double Z )
352 { 351 {
353 G4double m1 = particle->GetPDGMass(); 352 G4double m1 = particle->GetPDGMass();
354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 353 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
355 354
356 G4int iZ = static_cast<G4int>(Z+0.5); 355 G4int iZ = static_cast<G4int>(Z+0.5);
357 G4int iA = static_cast<G4int>(A+0.5); 356 G4int iA = static_cast<G4int>(A+0.5);
358 G4ParticleDefinition * theDef = 0; 357 G4ParticleDefinition * theDef = 0;
359 358
360 if (iZ == 1 && iA == 1) theDef = thePro 359 if (iZ == 1 && iA == 1) theDef = theProton;
361 else if (iZ == 1 && iA == 2) theDef = theDeu 360 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
362 else if (iZ == 1 && iA == 3) theDef = G4Trit 361 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
363 else if (iZ == 2 && iA == 3) theDef = G4He3: 362 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
364 else if (iZ == 2 && iA == 4) theDef = theAlp 363 else if (iZ == 2 && iA == 4) theDef = theAlpha;
365 else theDef = G4ParticleTable::GetParticleTa 364 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
366 365
367 G4double tmass = theDef->GetPDGMass(); 366 G4double tmass = theDef->GetPDGMass();
368 367
369 G4LorentzVector lv(0.0,0.0,0.0,tmass); 368 G4LorentzVector lv(0.0,0.0,0.0,tmass);
370 lv += lv1; 369 lv += lv1;
371 370
372 G4ThreeVector bst = lv.boostVector(); 371 G4ThreeVector bst = lv.boostVector();
373 lv1.boost(-bst); 372 lv1.boost(-bst);
374 373
375 G4ThreeVector p1 = lv1.vect(); 374 G4ThreeVector p1 = lv1.vect();
376 G4double ptot = p1.mag(); 375 G4double ptot = p1.mag();
377 G4double ptot2 = ptot*ptot; 376 G4double ptot2 = ptot*ptot;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 377 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
379 378
380 if( cost >= 1.0 ) cost = 1.0; 379 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0; 380 else if( cost <= -1.0) cost = -1.0;
382 381
383 G4double thetaCMS = std::acos(cost); 382 G4double thetaCMS = std::acos(cost);
384 383
385 G4double sigma = GetCoulombElasticXsc( parti 384 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
386 385
387 sigma *= pi/ptot2; 386 sigma *= pi/ptot2;
388 387
389 return sigma; 388 return sigma;
390 } 389 }
391 390
392 ////////////////////////////////////////////// 391 ////////////////////////////////////////////////////////////////////////////
393 // 392 //
394 // return differential elastic probability d(p 393 // return differential elastic probability d(probability)/d(omega)
395 394
396 G4double 395 G4double
397 G4NuclNuclDiffuseElastic::GetDiffElasticProb( 396 G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
398 G4doub 397 G4double theta
399 // G4double momentum, 398 // G4double momentum,
400 // G4double A 399 // G4double A
401 ) 400 )
402 { 401 {
403 G4double sigma, bzero, bzero2, bonebyarg, bo 402 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
404 G4double delta, diffuse, gamma; 403 G4double delta, diffuse, gamma;
405 G4double e1, e2, bone, bone2; 404 G4double e1, e2, bone, bone2;
406 405
407 // G4double wavek = momentum/hbarc; // wave 406 // G4double wavek = momentum/hbarc; // wave vector
408 // G4double r0 = 1.08*fermi; 407 // G4double r0 = 1.08*fermi;
409 // G4double rad = r0*G4Pow::GetInstance()- 408 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
410 409
411 G4double kr = fWaveVector*fNuclearRadius; 410 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
412 G4double kr2 = kr*kr; 411 G4double kr2 = kr*kr;
413 G4double krt = kr*theta; 412 G4double krt = kr*theta;
414 413
415 bzero = BesselJzero(krt); 414 bzero = BesselJzero(krt);
416 bzero2 = bzero*bzero; 415 bzero2 = bzero*bzero;
417 bone = BesselJone(krt); 416 bone = BesselJone(krt);
418 bone2 = bone*bone; 417 bone2 = bone*bone;
419 bonebyarg = BesselOneByArg(krt); 418 bonebyarg = BesselOneByArg(krt);
420 bonebyarg2 = bonebyarg*bonebyarg; 419 bonebyarg2 = bonebyarg*bonebyarg;
421 420
422 // VI - Coverity complains 421 // VI - Coverity complains
423 /* 422 /*
424 if (fParticle == theProton) 423 if (fParticle == theProton)
425 { 424 {
426 diffuse = 0.63*fermi; 425 diffuse = 0.63*fermi;
427 gamma = 0.3*fermi; 426 gamma = 0.3*fermi;
428 delta = 0.1*fermi*fermi; 427 delta = 0.1*fermi*fermi;
429 e1 = 0.3*fermi; 428 e1 = 0.3*fermi;
430 e2 = 0.35*fermi; 429 e2 = 0.35*fermi;
431 } 430 }
432 else // as proton, if were not defined 431 else // as proton, if were not defined
433 { 432 {
434 */ 433 */
435 diffuse = 0.63*fermi; 434 diffuse = 0.63*fermi;
436 gamma = 0.3*fermi; 435 gamma = 0.3*fermi;
437 delta = 0.1*fermi*fermi; 436 delta = 0.1*fermi*fermi;
438 e1 = 0.3*fermi; 437 e1 = 0.3*fermi;
439 e2 = 0.35*fermi; 438 e2 = 0.35*fermi;
440 //} 439 //}
441 G4double lambda = 15.; // 15 ok 440 G4double lambda = 15.; // 15 ok
442 441
443 // G4double kgamma = fWaveVector*gamma; 442 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
444 443
445 G4double kgamma = lambda*(1.-G4Exp(-fWave 444 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
446 G4double kgamma2 = kgamma*kgamma; 445 G4double kgamma2 = kgamma*kgamma;
447 446
448 // G4double dk2t = delta*fWaveVector*fWaveV 447 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
449 // G4double dk2t2 = dk2t*dk2t; 448 // G4double dk2t2 = dk2t*dk2t;
450 // G4double pikdt = pi*fWaveVector*diffuse*t 449 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
451 450
452 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 451 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
453 452
454 damp = DampFactor(pikdt); 453 damp = DampFactor(pikdt);
455 damp2 = damp*damp; 454 damp2 = damp*damp;
456 455
457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 456 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
458 G4double e2dk3t = -2.*e2*delta*fWaveVector* 457 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
459 458
460 459
461 sigma = kgamma2; 460 sigma = kgamma2;
462 // sigma += dk2t2; 461 // sigma += dk2t2;
463 sigma *= bzero2; 462 sigma *= bzero2;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 463 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2; 464 sigma += kr2*bonebyarg2;
466 sigma *= damp2; // *rad*rad; 465 sigma *= damp2; // *rad*rad;
467 466
468 return sigma; 467 return sigma;
469 } 468 }
470 469
471 ////////////////////////////////////////////// 470 ////////////////////////////////////////////////////////////////////////////
472 // 471 //
473 // return differential elastic probability d(p 472 // return differential elastic probability d(probability)/d(omega) with
474 // Coulomb correction 473 // Coulomb correction
475 474
476 G4double 475 G4double
477 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 476 G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
478 G4doub 477 G4double theta
479 // G4double momentum, 478 // G4double momentum,
480 // G4double A 479 // G4double A
481 ) 480 )
482 { 481 {
483 G4double sigma, bzero, bzero2, bonebyarg, bo 482 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
484 G4double delta, diffuse, gamma; 483 G4double delta, diffuse, gamma;
485 G4double e1, e2, bone, bone2; 484 G4double e1, e2, bone, bone2;
486 485
487 // G4double wavek = momentum/hbarc; // wave 486 // G4double wavek = momentum/hbarc; // wave vector
488 // G4double r0 = 1.08*fermi; 487 // G4double r0 = 1.08*fermi;
489 // G4double rad = r0*G4Pow::GetInstance()- 488 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
490 489
491 G4double kr = fWaveVector*fNuclearRadius; 490 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
492 G4double kr2 = kr*kr; 491 G4double kr2 = kr*kr;
493 G4double krt = kr*theta; 492 G4double krt = kr*theta;
494 493
495 bzero = BesselJzero(krt); 494 bzero = BesselJzero(krt);
496 bzero2 = bzero*bzero; 495 bzero2 = bzero*bzero;
497 bone = BesselJone(krt); 496 bone = BesselJone(krt);
498 bone2 = bone*bone; 497 bone2 = bone*bone;
499 bonebyarg = BesselOneByArg(krt); 498 bonebyarg = BesselOneByArg(krt);
500 bonebyarg2 = bonebyarg*bonebyarg; 499 bonebyarg2 = bonebyarg*bonebyarg;
501 500
502 if (fParticle == theProton) 501 if (fParticle == theProton)
503 { 502 {
504 diffuse = 0.63*fermi; 503 diffuse = 0.63*fermi;
505 // diffuse = 0.6*fermi; 504 // diffuse = 0.6*fermi;
506 gamma = 0.3*fermi; 505 gamma = 0.3*fermi;
507 delta = 0.1*fermi*fermi; 506 delta = 0.1*fermi*fermi;
508 e1 = 0.3*fermi; 507 e1 = 0.3*fermi;
509 e2 = 0.35*fermi; 508 e2 = 0.35*fermi;
510 } 509 }
511 else // as proton, if were not defined 510 else // as proton, if were not defined
512 { 511 {
513 diffuse = 0.63*fermi; 512 diffuse = 0.63*fermi;
514 gamma = 0.3*fermi; 513 gamma = 0.3*fermi;
515 delta = 0.1*fermi*fermi; 514 delta = 0.1*fermi*fermi;
516 e1 = 0.3*fermi; 515 e1 = 0.3*fermi;
517 e2 = 0.35*fermi; 516 e2 = 0.35*fermi;
518 } 517 }
519 G4double lambda = 15.; // 15 ok 518 G4double lambda = 15.; // 15 ok
520 // G4double kgamma = fWaveVector*gamma; 519 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
521 G4double kgamma = lambda*(1.-G4Exp(-fWave 520 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
522 521
523 // G4cout<<"kgamma = "<<kgamma<<G4endl; 522 // G4cout<<"kgamma = "<<kgamma<<G4endl;
524 523
525 if(fAddCoulomb) // add Coulomb correction 524 if(fAddCoulomb) // add Coulomb correction
526 { 525 {
527 G4double sinHalfTheta = std::sin(0.5*thet 526 G4double sinHalfTheta = std::sin(0.5*theta);
528 G4double sinHalfTheta2 = sinHalfTheta*sinH 527 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
529 528
530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 529 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 530 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
532 } 531 }
533 532
534 G4double kgamma2 = kgamma*kgamma; 533 G4double kgamma2 = kgamma*kgamma;
535 534
536 // G4double dk2t = delta*fWaveVector*fWaveV 535 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
537 // G4cout<<"dk2t = "<<dk2t<<G4endl; 536 // G4cout<<"dk2t = "<<dk2t<<G4endl;
538 // G4double dk2t2 = dk2t*dk2t; 537 // G4double dk2t2 = dk2t*dk2t;
539 // G4double pikdt = pi*fWaveVector*diffuse*t 538 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
540 539
541 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 540 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
542 541
543 // G4cout<<"pikdt = "<<pikdt<<G4endl; 542 // G4cout<<"pikdt = "<<pikdt<<G4endl;
544 543
545 damp = DampFactor(pikdt); 544 damp = DampFactor(pikdt);
546 damp2 = damp*damp; 545 damp2 = damp*damp;
547 546
548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 547 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
549 G4double e2dk3t = -2.*e2*delta*fWaveVector* 548 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
550 549
551 sigma = kgamma2; 550 sigma = kgamma2;
552 // sigma += dk2t2; 551 // sigma += dk2t2;
553 sigma *= bzero2; 552 sigma *= bzero2;
554 sigma += mode2k2*bone2; 553 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone; 554 sigma += e2dk3t*bzero*bone;
556 555
557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 556 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
558 sigma += kr2*bonebyarg2; // correction at J 557 sigma += kr2*bonebyarg2; // correction at J1()/()
559 558
560 sigma *= damp2; // *rad*rad; 559 sigma *= damp2; // *rad*rad;
561 560
562 return sigma; 561 return sigma;
563 } 562 }
564 563
565 564
566 ////////////////////////////////////////////// 565 ////////////////////////////////////////////////////////////////////////////
567 // 566 //
568 // return differential elastic probability d(p 567 // return differential elastic probability d(probability)/d(t) with
569 // Coulomb correction 568 // Coulomb correction
570 569
571 G4double 570 G4double
572 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 571 G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
573 { 572 {
574 G4double theta; 573 G4double theta;
575 574
576 theta = std::sqrt(alpha); 575 theta = std::sqrt(alpha);
577 576
578 // theta = std::acos( 1 - alpha/2. ); 577 // theta = std::acos( 1 - alpha/2. );
579 578
580 G4double sigma, bzero, bzero2, bonebyarg, bo 579 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
581 G4double delta, diffuse, gamma; 580 G4double delta, diffuse, gamma;
582 G4double e1, e2, bone, bone2; 581 G4double e1, e2, bone, bone2;
583 582
584 // G4double wavek = momentum/hbarc; // wave 583 // G4double wavek = momentum/hbarc; // wave vector
585 // G4double r0 = 1.08*fermi; 584 // G4double r0 = 1.08*fermi;
586 // G4double rad = r0*G4Pow::GetInstance()- 585 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
587 586
588 G4double kr = fWaveVector*fNuclearRadius; 587 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
589 G4double kr2 = kr*kr; 588 G4double kr2 = kr*kr;
590 G4double krt = kr*theta; 589 G4double krt = kr*theta;
591 590
592 bzero = BesselJzero(krt); 591 bzero = BesselJzero(krt);
593 bzero2 = bzero*bzero; 592 bzero2 = bzero*bzero;
594 bone = BesselJone(krt); 593 bone = BesselJone(krt);
595 bone2 = bone*bone; 594 bone2 = bone*bone;
596 bonebyarg = BesselOneByArg(krt); 595 bonebyarg = BesselOneByArg(krt);
597 bonebyarg2 = bonebyarg*bonebyarg; 596 bonebyarg2 = bonebyarg*bonebyarg;
598 597
599 if (fParticle == theProton) 598 if (fParticle == theProton)
600 { 599 {
601 diffuse = 0.63*fermi; 600 diffuse = 0.63*fermi;
602 // diffuse = 0.6*fermi; 601 // diffuse = 0.6*fermi;
603 gamma = 0.3*fermi; 602 gamma = 0.3*fermi;
604 delta = 0.1*fermi*fermi; 603 delta = 0.1*fermi*fermi;
605 e1 = 0.3*fermi; 604 e1 = 0.3*fermi;
606 e2 = 0.35*fermi; 605 e2 = 0.35*fermi;
607 } 606 }
608 else // as proton, if were not defined 607 else // as proton, if were not defined
609 { 608 {
610 diffuse = 0.63*fermi; 609 diffuse = 0.63*fermi;
611 gamma = 0.3*fermi; 610 gamma = 0.3*fermi;
612 delta = 0.1*fermi*fermi; 611 delta = 0.1*fermi*fermi;
613 e1 = 0.3*fermi; 612 e1 = 0.3*fermi;
614 e2 = 0.35*fermi; 613 e2 = 0.35*fermi;
615 } 614 }
616 G4double lambda = 15.; // 15 ok 615 G4double lambda = 15.; // 15 ok
617 // G4double kgamma = fWaveVector*gamma; 616 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
618 G4double kgamma = lambda*(1.-G4Exp(-fWave 617 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
619 618
620 // G4cout<<"kgamma = "<<kgamma<<G4endl; 619 // G4cout<<"kgamma = "<<kgamma<<G4endl;
621 620
622 if(fAddCoulomb) // add Coulomb correction 621 if(fAddCoulomb) // add Coulomb correction
623 { 622 {
624 G4double sinHalfTheta = theta*0.5; // std 623 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
625 G4double sinHalfTheta2 = sinHalfTheta*sinH 624 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
626 625
627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 626 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 627 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
629 } 628 }
630 629
631 G4double kgamma2 = kgamma*kgamma; 630 G4double kgamma2 = kgamma*kgamma;
632 631
633 // G4double dk2t = delta*fWaveVector*fWaveV 632 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
634 // G4cout<<"dk2t = "<<dk2t<<G4endl; 633 // G4cout<<"dk2t = "<<dk2t<<G4endl;
635 // G4double dk2t2 = dk2t*dk2t; 634 // G4double dk2t2 = dk2t*dk2t;
636 // G4double pikdt = pi*fWaveVector*diffuse*t 635 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
637 636
638 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 637 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
639 638
640 // G4cout<<"pikdt = "<<pikdt<<G4endl; 639 // G4cout<<"pikdt = "<<pikdt<<G4endl;
641 640
642 damp = DampFactor(pikdt); 641 damp = DampFactor(pikdt);
643 damp2 = damp*damp; 642 damp2 = damp*damp;
644 643
645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 644 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
646 G4double e2dk3t = -2.*e2*delta*fWaveVector* 645 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
647 646
648 sigma = kgamma2; 647 sigma = kgamma2;
649 // sigma += dk2t2; 648 // sigma += dk2t2;
650 sigma *= bzero2; 649 sigma *= bzero2;
651 sigma += mode2k2*bone2; 650 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone; 651 sigma += e2dk3t*bzero*bone;
653 652
654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 653 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
655 sigma += kr2*bonebyarg2; // correction at J 654 sigma += kr2*bonebyarg2; // correction at J1()/()
656 655
657 sigma *= damp2; // *rad*rad; 656 sigma *= damp2; // *rad*rad;
658 657
659 return sigma; 658 return sigma;
660 } 659 }
661 660
662 661
663 ////////////////////////////////////////////// 662 ////////////////////////////////////////////////////////////////////////////
664 // 663 //
665 // return differential elastic probability 2*p 664 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
666 665
667 G4double 666 G4double
668 G4NuclNuclDiffuseElastic::GetIntegrandFunction 667 G4NuclNuclDiffuseElastic::GetIntegrandFunction( G4double alpha )
669 { 668 {
670 G4double result; 669 G4double result;
671 670
672 result = GetDiffElasticSumProbA(alpha); 671 result = GetDiffElasticSumProbA(alpha);
673 672
674 // result *= 2*pi*std::sin(theta); 673 // result *= 2*pi*std::sin(theta);
675 674
676 return result; 675 return result;
677 } 676 }
678 677
679 ////////////////////////////////////////////// 678 ////////////////////////////////////////////////////////////////////////////
680 // 679 //
681 // return integral elastic cross section d(sig 680 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
682 681
683 G4double 682 G4double
684 G4NuclNuclDiffuseElastic::IntegralElasticProb( 683 G4NuclNuclDiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle,
685 G4doub 684 G4double theta,
686 G4double momentum, 685 G4double momentum,
687 G4doub 686 G4double A )
688 { 687 {
689 G4double result; 688 G4double result;
690 fParticle = particle; 689 fParticle = particle;
691 fWaveVector = momentum/hbarc; 690 fWaveVector = momentum/hbarc;
692 fAtomicWeight = A; 691 fAtomicWeight = A;
693 692
694 fNuclearRadius = CalculateNuclearRad(A); 693 fNuclearRadius = CalculateNuclearRad(A);
695 694
696 695
697 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 696 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
698 697
699 // result = integral.Legendre10(this,&G4Nucl 698 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
700 result = integral.Legendre96(this,&G4NuclNuc 699 result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
701 700
702 return result; 701 return result;
703 } 702 }
704 703
705 ////////////////////////////////////////////// 704 ////////////////////////////////////////////////////////////////////////////
706 // 705 //
707 // Return inv momentum transfer -t > 0 706 // Return inv momentum transfer -t > 0
708 707
709 G4double G4NuclNuclDiffuseElastic::SampleT( co 708 G4double G4NuclNuclDiffuseElastic::SampleT( const G4ParticleDefinition* aParticle,
710 G4double p, G4double A) 709 G4double p, G4double A)
711 { 710 {
712 G4double theta = SampleThetaCMS( aParticle, 711 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
713 G4double t = 2*p*p*( 1 - std::cos(theta) 712 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
714 return t; 713 return t;
715 } 714 }
716 715
717 ////////////////////////////////////////////// 716 ////////////////////////////////////////////////////////////////////////////
718 // 717 //
719 // Return scattering angle sampled in cms 718 // Return scattering angle sampled in cms
720 719
721 720
722 G4double 721 G4double
723 G4NuclNuclDiffuseElastic::SampleThetaCMS(const 722 G4NuclNuclDiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle,
724 G4doubl 723 G4double momentum, G4double A)
725 { 724 {
726 G4int i, iMax = 100; 725 G4int i, iMax = 100;
727 G4double norm, theta1, theta2, thetaMax; << 726 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
728 G4double result = 0., sum = 0.; <<
729 727
730 fParticle = particle; 728 fParticle = particle;
731 fWaveVector = momentum/hbarc; 729 fWaveVector = momentum/hbarc;
732 fAtomicWeight = A; 730 fAtomicWeight = A;
733 731
734 fNuclearRadius = CalculateNuclearRad(A); 732 fNuclearRadius = CalculateNuclearRad(A);
735 733
736 thetaMax = 10.174/fWaveVector/fNuclearRadius 734 thetaMax = 10.174/fWaveVector/fNuclearRadius;
737 735
738 if (thetaMax > pi) thetaMax = pi; 736 if (thetaMax > pi) thetaMax = pi;
739 737
740 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 738 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
741 739
742 // result = integral.Legendre10(this,&G4Nucl 740 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
743 norm = integral.Legendre96(this,&G4NuclNuclD 741 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
744 742
745 norm *= G4UniformRand(); 743 norm *= G4UniformRand();
746 744
747 for(i = 1; i <= iMax; i++) 745 for(i = 1; i <= iMax; i++)
748 { 746 {
749 theta1 = (i-1)*thetaMax/iMax; 747 theta1 = (i-1)*thetaMax/iMax;
750 theta2 = i*thetaMax/iMax; 748 theta2 = i*thetaMax/iMax;
751 sum += integral.Legendre10(this,&G4NuclN 749 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
752 750
753 if ( sum >= norm ) 751 if ( sum >= norm )
754 { 752 {
755 result = 0.5*(theta1 + theta2); 753 result = 0.5*(theta1 + theta2);
756 break; 754 break;
757 } 755 }
758 } 756 }
759 if (i > iMax ) result = 0.5*(theta1 + theta2 757 if (i > iMax ) result = 0.5*(theta1 + theta2);
760 758
761 G4double sigma = pi*thetaMax/iMax; 759 G4double sigma = pi*thetaMax/iMax;
762 760
763 result += G4RandGauss::shoot(0.,sigma); 761 result += G4RandGauss::shoot(0.,sigma);
764 762
765 if(result < 0.) result = 0.; 763 if(result < 0.) result = 0.;
766 if(result > thetaMax) result = thetaMax; 764 if(result > thetaMax) result = thetaMax;
767 765
768 return result; 766 return result;
769 } 767 }
770 768
771 ////////////////////////////////////////////// 769 /////////////////////////////////////////////////////////////////////////////
772 ///////////////////// Table preparation and r 770 ///////////////////// Table preparation and reading ////////////////////////
773 771
774 ////////////////////////////////////////////// 772 ////////////////////////////////////////////////////////////////////////////
775 // 773 //
776 // Interface function. Return inv momentum tra << 774 // Return inv momentum transfer -t > 0 from initialisation table
777 775
778 G4double G4NuclNuclDiffuseElastic::SampleInvar 776 G4double G4NuclNuclDiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
779 777 G4int Z, G4int A)
780 { 778 {
781 fParticle = aParticle; 779 fParticle = aParticle;
782 fAtomicNumber = G4double(Z); << 780 G4double m1 = fParticle->GetPDGMass();
783 fAtomicWeight = G4double(A); <<
784 <<
785 G4double t(0.), m1 = fParticle->GetPDGMass() <<
786 G4double totElab = std::sqrt(m1*m1+p*p); 781 G4double totElab = std::sqrt(m1*m1+p*p);
787 G4double mass2 = G4NucleiProperties::GetNucl 782 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
788 G4LorentzVector lv1(p,0.0,0.0,totElab); 783 G4LorentzVector lv1(p,0.0,0.0,totElab);
789 G4LorentzVector lv(0.0,0.0,0.0,mass2); 784 G4LorentzVector lv(0.0,0.0,0.0,mass2);
790 lv += lv1; 785 lv += lv1;
791 786
792 G4ThreeVector bst = lv.boostVector(); 787 G4ThreeVector bst = lv.boostVector();
793 lv1.boost(-bst); 788 lv1.boost(-bst);
794 789
795 G4ThreeVector p1 = lv1.vect(); 790 G4ThreeVector p1 = lv1.vect();
796 G4double momentumCMS = p1.mag(); 791 G4double momentumCMS = p1.mag();
797 792
798 // t = SampleTableT( aParticle, momentumCMS << 793 G4double t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
799 <<
800 t = SampleCoulombMuCMS( aParticle, momentum <<
801 <<
802 <<
803 return t; 794 return t;
804 } 795 }
805 796
806 ////////////////////////////////////////////// 797 ////////////////////////////////////////////////////////////////////////////
807 // 798 //
808 // Return inv momentum transfer -t > 0 as Coul <<
809 <<
810 G4double G4NuclNuclDiffuseElastic::SampleCoulo <<
811 { <<
812 G4double t(0.), rand(0.), mu(0.); <<
813 <<
814 G4double A1 = G4double( aParticle->GetBaryon <<
815 G4double R1 = CalculateNuclearRad(A1); <<
816 <<
817 fNuclearRadius = CalculateNuclearRad(fAtomi <<
818 fNuclearRadius += R1; <<
819 <<
820 InitDynParameters(fParticle, p); <<
821 <<
822 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutT <<
823 <<
824 rand = G4UniformRand(); <<
825 <<
826 // sample (1-cosTheta) in 0 < Theta < ThetaC <<
827 <<
828 mu = fCoulombMuC*rand*fAm; <<
829 mu /= fAm + fCoulombMuC*( 1. - rand ); <<
830 <<
831 t = 4.*p*p*mu; <<
832 <<
833 return t; <<
834 } <<
835 <<
836 <<
837 ////////////////////////////////////////////// <<
838 // <<
839 // Return inv momentum transfer -t > 0 from in 799 // Return inv momentum transfer -t > 0 from initialisation table
840 800
841 G4double G4NuclNuclDiffuseElastic::SampleTable 801 G4double G4NuclNuclDiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
842 802 G4double Z, G4double A)
843 { 803 {
844 G4double alpha = SampleTableThetaCMS( aParti 804 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
845 // G4double t = 2*p*p*( 1 - std::cos(std 805 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
846 G4double t = p*p*alpha; // - 806 G4double t = p*p*alpha; // -t !!!
847 return t; 807 return t;
848 } 808 }
849 809
850 ////////////////////////////////////////////// 810 ////////////////////////////////////////////////////////////////////////////
851 // 811 //
852 // Return scattering angle2 sampled in cms acc 812 // Return scattering angle2 sampled in cms according to precalculated table.
853 813
854 814
855 G4double 815 G4double
856 G4NuclNuclDiffuseElastic::SampleTableThetaCMS( 816 G4NuclNuclDiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
857 G4doubl 817 G4double momentum, G4double Z, G4double A)
858 { 818 {
859 std::size_t iElement; << 819 size_t iElement;
860 G4int iMomentum, iAngle; 820 G4int iMomentum, iAngle;
861 G4double randAngle, position, theta1, theta2 821 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
862 G4double m1 = particle->GetPDGMass(); 822 G4double m1 = particle->GetPDGMass();
863 823
864 for(iElement = 0; iElement < fElementNumberV 824 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
865 { 825 {
866 if( std::fabs(Z - fElementNumberVector[iEl 826 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
867 } 827 }
868 if ( iElement == fElementNumberVector.size() 828 if ( iElement == fElementNumberVector.size() )
869 { 829 {
870 InitialiseOnFly(Z,A); // table preparation 830 InitialiseOnFly(Z,A); // table preparation, if needed
871 831
872 // iElement--; 832 // iElement--;
873 833
874 // G4cout << "G4NuclNuclDiffuseElastic: El 834 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
875 // << " is not found, return zero angle" < 835 // << " is not found, return zero angle" << G4endl;
876 // return 0.; // no table for this element 836 // return 0.; // no table for this element
877 } 837 }
878 // G4cout<<"iElement = "<<iElement<<G4endl; 838 // G4cout<<"iElement = "<<iElement<<G4endl;
879 839
880 fAngleTable = fAngleBank[iElement]; 840 fAngleTable = fAngleBank[iElement];
881 841
882 G4double kinE = std::sqrt(momentum*momentum 842 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
883 843
884 for( iMomentum = 0; iMomentum < fEnergyBin; 844 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
885 { 845 {
886 // G4cout<<iMomentum<<", kinE = "<<kinE/Me 846 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
887 847
888 if( kinE < fEnergyVector->GetLowEdgeEnergy 848 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
889 } 849 }
890 // G4cout<<"iMomentum = "<<iMomentum<<", fEn 850 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
891 851
892 852
893 if ( iMomentum >= fEnergyBin ) iMomentum = f 853 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
894 if ( iMomentum < 0 ) iMomentum = 0 854 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
895 855
896 856
897 if (iMomentum == fEnergyBin -1 || iMomentum 857 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
898 { 858 {
899 position = (*(*fAngleTable)(iMomentum))(fA 859 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
900 860
901 // G4cout<<"position = "<<position<<G4endl 861 // G4cout<<"position = "<<position<<G4endl;
902 862
903 for(iAngle = 0; iAngle < fAngleBin-1; iAng 863 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
904 { 864 {
905 if( position < (*(*fAngleTable)(iMomentu 865 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
906 } 866 }
907 if (iAngle >= fAngleBin-1) iAngle = fAngle 867 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
908 868
909 // G4cout<<"iAngle = "<<iAngle<<G4endl; 869 // G4cout<<"iAngle = "<<iAngle<<G4endl;
910 870
911 randAngle = GetScatteringAngle(iMomentum, 871 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
912 872
913 // G4cout<<"randAngle = "<<randAngle<<G4en 873 // G4cout<<"randAngle = "<<randAngle<<G4endl;
914 } 874 }
915 else // kinE inside between energy table ed 875 else // kinE inside between energy table edges
916 { 876 {
917 // position = (*(*fAngleTable)(iMomentum)) 877 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
918 position = (*(*fAngleTable)(iMomentum))(0) 878 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
919 879
920 // G4cout<<"position = "<<position<<G4endl 880 // G4cout<<"position = "<<position<<G4endl;
921 881
922 for(iAngle = 0; iAngle < fAngleBin-1; iAng 882 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
923 { 883 {
924 // if( position < (*(*fAngleTable)(iMome 884 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
925 if( position > (*(*fAngleTable)(iMomentu 885 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
926 } 886 }
927 if (iAngle >= fAngleBin-1) iAngle = fAngle 887 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
928 888
929 // G4cout<<"iAngle = "<<iAngle<<G4endl; 889 // G4cout<<"iAngle = "<<iAngle<<G4endl;
930 890
931 theta2 = GetScatteringAngle(iMomentum, iA 891 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
932 892
933 // G4cout<<"theta2 = "<<theta2<<G4endl; 893 // G4cout<<"theta2 = "<<theta2<<G4endl;
934 894
935 E2 = fEnergyVector->GetLowEdgeEnergy(iMome 895 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
936 896
937 // G4cout<<"E2 = "<<E2<<G4endl; 897 // G4cout<<"E2 = "<<E2<<G4endl;
938 898
939 iMomentum--; 899 iMomentum--;
940 900
941 // position = (*(*fAngleTable)(iMomentum)) 901 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
942 902
943 // G4cout<<"position = "<<position<<G4endl 903 // G4cout<<"position = "<<position<<G4endl;
944 904
945 for(iAngle = 0; iAngle < fAngleBin-1; iAng 905 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
946 { 906 {
947 // if( position < (*(*fAngleTable)(iMome 907 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
948 if( position > (*(*fAngleTable)(iMomentu 908 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
949 } 909 }
950 if (iAngle >= fAngleBin-1) iAngle = fAngle 910 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
951 911
952 theta1 = GetScatteringAngle(iMomentum, iA 912 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
953 913
954 // G4cout<<"theta1 = "<<theta1<<G4endl; 914 // G4cout<<"theta1 = "<<theta1<<G4endl;
955 915
956 E1 = fEnergyVector->GetLowEdgeEnergy(iMome 916 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
957 917
958 // G4cout<<"E1 = "<<E1<<G4endl; 918 // G4cout<<"E1 = "<<E1<<G4endl;
959 919
960 W = 1.0/(E2 - E1); 920 W = 1.0/(E2 - E1);
961 W1 = (E2 - kinE)*W; 921 W1 = (E2 - kinE)*W;
962 W2 = (kinE - E1)*W; 922 W2 = (kinE - E1)*W;
963 923
964 randAngle = W1*theta1 + W2*theta2; 924 randAngle = W1*theta1 + W2*theta2;
965 925
966 // randAngle = theta2; 926 // randAngle = theta2;
967 } 927 }
968 // G4double angle = randAngle; 928 // G4double angle = randAngle;
969 // if (randAngle > 0.) randAngle /= 2*pi*std 929 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
970 // G4cout<<"randAngle = "<<randAngle/degree< 930 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
971 931
972 return randAngle; 932 return randAngle;
973 } 933 }
974 934
975 ////////////////////////////////////////////// 935 //////////////////////////////////////////////////////////////////////////////
976 // 936 //
977 // Initialisation for given particle on fly us 937 // Initialisation for given particle on fly using new element number
978 938
979 void G4NuclNuclDiffuseElastic::InitialiseOnFly 939 void G4NuclNuclDiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
980 { 940 {
981 fAtomicNumber = Z; // atomic number 941 fAtomicNumber = Z; // atomic number
982 fAtomicWeight = G4NistManager::Instance()-> 942 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
983 943
984 G4double A1 = G4double( fParticle->GetBaryon 944 G4double A1 = G4double( fParticle->GetBaryonNumber() );
985 G4double R1 = CalculateNuclearRad(A1); 945 G4double R1 = CalculateNuclearRad(A1);
986 946
987 fNuclearRadius = CalculateNuclearRad(fAtomic 947 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
988 948
989 if( verboseLevel > 0 ) 949 if( verboseLevel > 0 )
990 { 950 {
991 G4cout<<"G4NuclNuclDiffuseElastic::Initial 951 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
992 <<Z<<"; and A = "<<A<<G4endl; 952 <<Z<<"; and A = "<<A<<G4endl;
993 } 953 }
994 fElementNumberVector.push_back(fAtomicNumber 954 fElementNumberVector.push_back(fAtomicNumber);
995 955
996 BuildAngleTable(); 956 BuildAngleTable();
997 957
998 fAngleBank.push_back(fAngleTable); 958 fAngleBank.push_back(fAngleTable);
999 959
1000 return; 960 return;
1001 } 961 }
1002 962
1003 ///////////////////////////////////////////// 963 ///////////////////////////////////////////////////////////////////////////////
1004 // 964 //
1005 // Build for given particle and element table 965 // Build for given particle and element table of momentum, angle probability.
1006 // For the moment in lab system. 966 // For the moment in lab system.
1007 967
1008 void G4NuclNuclDiffuseElastic::BuildAngleTabl 968 void G4NuclNuclDiffuseElastic::BuildAngleTable()
1009 { 969 {
1010 G4int i, j; 970 G4int i, j;
1011 G4double partMom, kinE, m1 = fParticle->Get 971 G4double partMom, kinE, m1 = fParticle->GetPDGMass();
1012 G4double alpha1, alpha2, alphaMax, alphaCou 972 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1013 973
1014 // G4cout<<"particle z = "<<z<<"; particle 974 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
1015 975
1016 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 976 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1017 977
1018 fAngleTable = new G4PhysicsTable(fEnergyBin 978 fAngleTable = new G4PhysicsTable(fEnergyBin);
1019 979
1020 for( i = 0; i < fEnergyBin; i++) 980 for( i = 0; i < fEnergyBin; i++)
1021 { 981 {
1022 kinE = fEnergyVector->GetLowEdgeEn 982 kinE = fEnergyVector->GetLowEdgeEnergy(i);
1023 983
1024 // G4cout<<G4endl; 984 // G4cout<<G4endl;
1025 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G 985 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
1026 986
1027 partMom = std::sqrt( kinE*(kinE + 2*m 987 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1028 988
1029 InitDynParameters(fParticle, partMom); 989 InitDynParameters(fParticle, partMom);
1030 990
1031 alphaMax = fRutherfordTheta*fCofAlphaMax; 991 alphaMax = fRutherfordTheta*fCofAlphaMax;
1032 992
1033 if(alphaMax > pi) alphaMax = pi; 993 if(alphaMax > pi) alphaMax = pi;
1034 994
1035 // VI: Coverity complain 995 // VI: Coverity complain
1036 //alphaMax = pi2; 996 //alphaMax = pi2;
1037 997
1038 alphaCoulomb = fRutherfordTheta*fCofAlpha 998 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
1039 999
1040 // G4cout<<"alphaCoulomb = "<<alphaCoulom 1000 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1041 1001
1042 G4PhysicsFreeVector* angleVector = new G4 1002 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1043 1003
1044 // G4PhysicsLogVector* angleBins = new G 1004 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1045 1005
1046 G4double delth = (alphaMax-alphaCoulomb)/ 1006 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1047 1007
1048 sum = 0.; 1008 sum = 0.;
1049 1009
1050 // fAddCoulomb = false; 1010 // fAddCoulomb = false;
1051 fAddCoulomb = true; 1011 fAddCoulomb = true;
1052 1012
1053 // for(j = 1; j < fAngleBin; j++) 1013 // for(j = 1; j < fAngleBin; j++)
1054 for(j = fAngleBin-1; j >= 1; j--) 1014 for(j = fAngleBin-1; j >= 1; j--)
1055 { 1015 {
1056 // alpha1 = angleBins->GetLowEdgeEnergy 1016 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1057 // alpha2 = angleBins->GetLowEdgeEnergy 1017 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1058 1018
1059 // alpha1 = alphaMax*(j-1)/fAngleBin; 1019 // alpha1 = alphaMax*(j-1)/fAngleBin;
1060 // alpha2 = alphaMax*( j )/fAngleBin; 1020 // alpha2 = alphaMax*( j )/fAngleBin;
1061 1021
1062 alpha1 = alphaCoulomb + delth*(j-1); 1022 alpha1 = alphaCoulomb + delth*(j-1);
1063 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1023 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1064 alpha2 = alpha1 + delth; 1024 alpha2 = alpha1 + delth;
1065 1025
1066 delta = integral.Legendre10(this, &G4Nu 1026 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1067 // delta = integral.Legendre96(this, &G 1027 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1068 1028
1069 sum += delta; 1029 sum += delta;
1070 1030
1071 angleVector->PutValue( j-1 , alpha1, su 1031 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1072 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = " 1032 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1073 } 1033 }
1074 fAngleTable->insertAt(i,angleVector); 1034 fAngleTable->insertAt(i,angleVector);
1075 1035
1076 // delete[] angleVector; 1036 // delete[] angleVector;
1077 // delete[] angleBins; 1037 // delete[] angleBins;
1078 } 1038 }
1079 return; 1039 return;
1080 } 1040 }
1081 1041
1082 ///////////////////////////////////////////// 1042 /////////////////////////////////////////////////////////////////////////////////
1083 // 1043 //
1084 // 1044 //
1085 1045
1086 G4double 1046 G4double
1087 G4NuclNuclDiffuseElastic:: GetScatteringAngle 1047 G4NuclNuclDiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1088 { 1048 {
1089 G4double x1, x2, y1, y2, randAngle; 1049 G4double x1, x2, y1, y2, randAngle;
1090 1050
1091 if( iAngle == 0 ) 1051 if( iAngle == 0 )
1092 { 1052 {
1093 randAngle = (*fAngleTable)(iMomentum)->Ge 1053 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1094 // iAngle++; 1054 // iAngle++;
1095 } 1055 }
1096 else 1056 else
1097 { 1057 {
1098 if ( iAngle >= G4int((*fAngleTable)(iMome 1058 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1099 { 1059 {
1100 iAngle = G4int((*fAngleTable)(iMomentum << 1060 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1101 } 1061 }
1102 y1 = (*(*fAngleTable)(iMomentum))(iAngle- 1062 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1103 y2 = (*(*fAngleTable)(iMomentum))(iAngle) 1063 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1104 1064
1105 x1 = (*fAngleTable)(iMomentum)->GetLowEdg 1065 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1106 x2 = (*fAngleTable)(iMomentum)->GetLowEdg 1066 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1107 1067
1108 if ( x1 == x2 ) randAngle = x2; 1068 if ( x1 == x2 ) randAngle = x2;
1109 else 1069 else
1110 { 1070 {
1111 if ( y1 == y2 ) randAngle = x1 + ( x2 - 1071 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1112 else 1072 else
1113 { 1073 {
1114 randAngle = x1 + ( position - y1 )*( 1074 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1115 } 1075 }
1116 } 1076 }
1117 } 1077 }
1118 return randAngle; 1078 return randAngle;
1119 } 1079 }
1120 1080
1121 1081
1122 1082
1123 ///////////////////////////////////////////// 1083 ////////////////////////////////////////////////////////////////////////////
1124 // 1084 //
1125 // Return scattering angle sampled in lab sys 1085 // Return scattering angle sampled in lab system (target at rest)
1126 1086
1127 1087
1128 1088
1129 G4double 1089 G4double
1130 G4NuclNuclDiffuseElastic::SampleThetaLab( con 1090 G4NuclNuclDiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
1131 G4dou 1091 G4double tmass, G4double A)
1132 { 1092 {
1133 const G4ParticleDefinition* theParticle = a 1093 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1134 G4double m1 = theParticle->GetPDGMass(); 1094 G4double m1 = theParticle->GetPDGMass();
1135 G4double plab = aParticle->GetTotalMomentum 1095 G4double plab = aParticle->GetTotalMomentum();
1136 G4LorentzVector lv1 = aParticle->Get4Moment 1096 G4LorentzVector lv1 = aParticle->Get4Momentum();
1137 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1097 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1138 lv += lv1; 1098 lv += lv1;
1139 1099
1140 G4ThreeVector bst = lv.boostVector(); 1100 G4ThreeVector bst = lv.boostVector();
1141 lv1.boost(-bst); 1101 lv1.boost(-bst);
1142 1102
1143 G4ThreeVector p1 = lv1.vect(); 1103 G4ThreeVector p1 = lv1.vect();
1144 G4double ptot = p1.mag(); 1104 G4double ptot = p1.mag();
1145 G4double tmax = 4.0*ptot*ptot; 1105 G4double tmax = 4.0*ptot*ptot;
1146 G4double t = 0.0; 1106 G4double t = 0.0;
1147 1107
1148 1108
1149 // 1109 //
1150 // Sample t 1110 // Sample t
1151 // 1111 //
1152 1112
1153 t = SampleT( theParticle, ptot, A); 1113 t = SampleT( theParticle, ptot, A);
1154 1114
1155 // NaN finder 1115 // NaN finder
1156 if(!(t < 0.0 || t >= 0.0)) 1116 if(!(t < 0.0 || t >= 0.0))
1157 { 1117 {
1158 if (verboseLevel > 0) 1118 if (verboseLevel > 0)
1159 { 1119 {
1160 G4cout << "G4NuclNuclDiffuseElastic:WAR 1120 G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1161 << " mom(GeV)= " << plab/GeV 1121 << " mom(GeV)= " << plab/GeV
1162 << " S-wave will be sampled" 1122 << " S-wave will be sampled"
1163 << G4endl; 1123 << G4endl;
1164 } 1124 }
1165 t = G4UniformRand()*tmax; 1125 t = G4UniformRand()*tmax;
1166 } 1126 }
1167 if(verboseLevel>1) 1127 if(verboseLevel>1)
1168 { 1128 {
1169 G4cout <<" t= " << t << " tmax= " << tmax 1129 G4cout <<" t= " << t << " tmax= " << tmax
1170 << " ptot= " << ptot << G4endl; 1130 << " ptot= " << ptot << G4endl;
1171 } 1131 }
1172 // Sampling of angles in CM system 1132 // Sampling of angles in CM system
1173 1133
1174 G4double phi = G4UniformRand()*twopi; 1134 G4double phi = G4UniformRand()*twopi;
1175 G4double cost = 1. - 2.0*t/tmax; 1135 G4double cost = 1. - 2.0*t/tmax;
1176 G4double sint; 1136 G4double sint;
1177 1137
1178 if( cost >= 1.0 ) 1138 if( cost >= 1.0 )
1179 { 1139 {
1180 cost = 1.0; 1140 cost = 1.0;
1181 sint = 0.0; 1141 sint = 0.0;
1182 } 1142 }
1183 else if( cost <= -1.0) 1143 else if( cost <= -1.0)
1184 { 1144 {
1185 cost = -1.0; 1145 cost = -1.0;
1186 sint = 0.0; 1146 sint = 0.0;
1187 } 1147 }
1188 else 1148 else
1189 { 1149 {
1190 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1150 sint = std::sqrt((1.0-cost)*(1.0+cost));
1191 } 1151 }
1192 if (verboseLevel>1) 1152 if (verboseLevel>1)
1193 { 1153 {
1194 G4cout << "cos(t)=" << cost << " std::sin 1154 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1195 } 1155 }
1196 G4ThreeVector v1(sint*std::cos(phi),sint*st 1156 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1197 v1 *= ptot; 1157 v1 *= ptot;
1198 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1158 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1199 1159
1200 nlv1.boost(bst); 1160 nlv1.boost(bst);
1201 1161
1202 G4ThreeVector np1 = nlv1.vect(); 1162 G4ThreeVector np1 = nlv1.vect();
1203 1163
1204 // G4double theta = std::acos( np1.z()/np 1164 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1205 1165
1206 G4double theta = np1.theta(); 1166 G4double theta = np1.theta();
1207 1167
1208 return theta; 1168 return theta;
1209 } 1169 }
1210 1170
1211 ///////////////////////////////////////////// 1171 ////////////////////////////////////////////////////////////////////////////
1212 // 1172 //
1213 // Return scattering angle in lab system (tar 1173 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1214 1174
1215 1175
1216 1176
1217 G4double 1177 G4double
1218 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( 1178 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
1219 G4dou 1179 G4double tmass, G4double thetaCMS)
1220 { 1180 {
1221 const G4ParticleDefinition* theParticle = a 1181 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1222 G4double m1 = theParticle->GetPDGMass(); 1182 G4double m1 = theParticle->GetPDGMass();
1223 // G4double plab = aParticle->GetTotalMomen 1183 // G4double plab = aParticle->GetTotalMomentum();
1224 G4LorentzVector lv1 = aParticle->Get4Moment 1184 G4LorentzVector lv1 = aParticle->Get4Momentum();
1225 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1185 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1226 1186
1227 lv += lv1; 1187 lv += lv1;
1228 1188
1229 G4ThreeVector bst = lv.boostVector(); 1189 G4ThreeVector bst = lv.boostVector();
1230 1190
1231 lv1.boost(-bst); 1191 lv1.boost(-bst);
1232 1192
1233 G4ThreeVector p1 = lv1.vect(); 1193 G4ThreeVector p1 = lv1.vect();
1234 G4double ptot = p1.mag(); 1194 G4double ptot = p1.mag();
1235 1195
1236 G4double phi = G4UniformRand()*twopi; 1196 G4double phi = G4UniformRand()*twopi;
1237 G4double cost = std::cos(thetaCMS); 1197 G4double cost = std::cos(thetaCMS);
1238 G4double sint; 1198 G4double sint;
1239 1199
1240 if( cost >= 1.0 ) 1200 if( cost >= 1.0 )
1241 { 1201 {
1242 cost = 1.0; 1202 cost = 1.0;
1243 sint = 0.0; 1203 sint = 0.0;
1244 } 1204 }
1245 else if( cost <= -1.0) 1205 else if( cost <= -1.0)
1246 { 1206 {
1247 cost = -1.0; 1207 cost = -1.0;
1248 sint = 0.0; 1208 sint = 0.0;
1249 } 1209 }
1250 else 1210 else
1251 { 1211 {
1252 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1212 sint = std::sqrt((1.0-cost)*(1.0+cost));
1253 } 1213 }
1254 if (verboseLevel>1) 1214 if (verboseLevel>1)
1255 { 1215 {
1256 G4cout << "cos(tcms)=" << cost << " std:: 1216 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1257 } 1217 }
1258 G4ThreeVector v1(sint*std::cos(phi),sint*st 1218 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1259 v1 *= ptot; 1219 v1 *= ptot;
1260 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1220 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1261 1221
1262 nlv1.boost(bst); 1222 nlv1.boost(bst);
1263 1223
1264 G4ThreeVector np1 = nlv1.vect(); 1224 G4ThreeVector np1 = nlv1.vect();
1265 1225
1266 1226
1267 G4double thetaLab = np1.theta(); 1227 G4double thetaLab = np1.theta();
1268 1228
1269 return thetaLab; 1229 return thetaLab;
1270 } 1230 }
1271 1231
1272 ///////////////////////////////////////////// 1232 ////////////////////////////////////////////////////////////////////////////
1273 // 1233 //
1274 // Return scattering angle in CMS system (tar 1234 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1275 1235
1276 1236
1277 1237
1278 G4double 1238 G4double
1279 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( 1239 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
1280 G4dou 1240 G4double tmass, G4double thetaLab)
1281 { 1241 {
1282 const G4ParticleDefinition* theParticle = a 1242 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1283 G4double m1 = theParticle->GetPDGMass(); 1243 G4double m1 = theParticle->GetPDGMass();
1284 G4double plab = aParticle->GetTotalMomentum 1244 G4double plab = aParticle->GetTotalMomentum();
1285 G4LorentzVector lv1 = aParticle->Get4Moment 1245 G4LorentzVector lv1 = aParticle->Get4Momentum();
1286 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1246 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1287 1247
1288 lv += lv1; 1248 lv += lv1;
1289 1249
1290 G4ThreeVector bst = lv.boostVector(); 1250 G4ThreeVector bst = lv.boostVector();
1291 1251
1292 // lv1.boost(-bst); 1252 // lv1.boost(-bst);
1293 1253
1294 // G4ThreeVector p1 = lv1.vect(); 1254 // G4ThreeVector p1 = lv1.vect();
1295 // G4double ptot = p1.mag(); 1255 // G4double ptot = p1.mag();
1296 1256
1297 G4double phi = G4UniformRand()*twopi; 1257 G4double phi = G4UniformRand()*twopi;
1298 G4double cost = std::cos(thetaLab); 1258 G4double cost = std::cos(thetaLab);
1299 G4double sint; 1259 G4double sint;
1300 1260
1301 if( cost >= 1.0 ) 1261 if( cost >= 1.0 )
1302 { 1262 {
1303 cost = 1.0; 1263 cost = 1.0;
1304 sint = 0.0; 1264 sint = 0.0;
1305 } 1265 }
1306 else if( cost <= -1.0) 1266 else if( cost <= -1.0)
1307 { 1267 {
1308 cost = -1.0; 1268 cost = -1.0;
1309 sint = 0.0; 1269 sint = 0.0;
1310 } 1270 }
1311 else 1271 else
1312 { 1272 {
1313 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1273 sint = std::sqrt((1.0-cost)*(1.0+cost));
1314 } 1274 }
1315 if (verboseLevel>1) 1275 if (verboseLevel>1)
1316 { 1276 {
1317 G4cout << "cos(tlab)=" << cost << " std:: 1277 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1318 } 1278 }
1319 G4ThreeVector v1(sint*std::cos(phi),sint*st 1279 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1320 v1 *= plab; 1280 v1 *= plab;
1321 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1281 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1322 1282
1323 nlv1.boost(-bst); 1283 nlv1.boost(-bst);
1324 1284
1325 G4ThreeVector np1 = nlv1.vect(); 1285 G4ThreeVector np1 = nlv1.vect();
1326 1286
1327 1287
1328 G4double thetaCMS = np1.theta(); 1288 G4double thetaCMS = np1.theta();
1329 1289
1330 return thetaCMS; 1290 return thetaCMS;
1331 } 1291 }
1332 1292
1333 ///////////////////////////////////////////// 1293 ///////////////////////////////////////////////////////////////////////////////
1334 // 1294 //
1335 // Test for given particle and element table 1295 // Test for given particle and element table of momentum, angle probability.
1336 // For the moment in lab system. 1296 // For the moment in lab system.
1337 1297
1338 void G4NuclNuclDiffuseElastic::TestAngleTable 1298 void G4NuclNuclDiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1339 G 1299 G4double Z, G4double A)
1340 { 1300 {
1341 fAtomicNumber = Z; // atomic number 1301 fAtomicNumber = Z; // atomic number
1342 fAtomicWeight = A; // number of nucleo 1302 fAtomicWeight = A; // number of nucleons
1343 fNuclearRadius = CalculateNuclearRad(fAtomi 1303 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1344 1304
1345 1305
1346 1306
1347 G4cout<<"G4NuclNuclDiffuseElastic::TestAngl 1307 G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1348 <<Z<<"; and A = "<<A<<G4endl; 1308 <<Z<<"; and A = "<<A<<G4endl;
1349 1309
1350 fElementNumberVector.push_back(fAtomicNumbe 1310 fElementNumberVector.push_back(fAtomicNumber);
1351 1311
1352 1312
1353 1313
1354 1314
1355 G4int i=0, j; 1315 G4int i=0, j;
1356 G4double a = 0., z = theParticle->GetPDGCha 1316 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1357 G4double alpha1=0., alpha2=0., alphaMax=0., 1317 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1358 G4double deltaL10 = 0., deltaL96 = 0., delt 1318 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1359 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0. 1319 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1360 G4double epsilon = 0.001; 1320 G4double epsilon = 0.001;
1361 1321
1362 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1322 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1363 1323
1364 fAngleTable = new G4PhysicsTable(fEnergyBin 1324 fAngleTable = new G4PhysicsTable(fEnergyBin);
1365 1325
1366 fWaveVector = partMom/hbarc; 1326 fWaveVector = partMom/hbarc;
1367 1327
1368 G4double kR = fWaveVector*fNuclearRadiu 1328 G4double kR = fWaveVector*fNuclearRadius;
1369 G4double kR2 = kR*kR; 1329 G4double kR2 = kR*kR;
1370 G4double kRmax = 10.6; // 10.6, 18, 10.174 1330 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1371 G4double kRcoul = 1.2; // 1.4, 2.5; // on t 1331 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1372 1332
1373 alphaMax = kRmax*kRmax/kR2; 1333 alphaMax = kRmax*kRmax/kR2;
1374 1334
1375 if (alphaMax > 4.) alphaMax = 4.; // vmg05 1335 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1376 1336
1377 alphaCoulomb = kRcoul*kRcoul/kR2; 1337 alphaCoulomb = kRcoul*kRcoul/kR2;
1378 1338
1379 if( z ) 1339 if( z )
1380 { 1340 {
1381 a = partMom/m1; // beta*gamma 1341 a = partMom/m1; // beta*gamma for m1
1382 fBeta = a/std::sqrt(1+a*a); 1342 fBeta = a/std::sqrt(1+a*a);
1383 fZommerfeld = CalculateZommerfeld( fBet 1343 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1384 fAm = CalculateAm( partMom, fZo 1344 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1385 } 1345 }
1386 G4PhysicsFreeVector* angleVector = new G4Ph 1346 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1387 1347
1388 // G4PhysicsLogVector* angleBins = new G4P 1348 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1389 1349
1390 1350
1391 fAddCoulomb = false; 1351 fAddCoulomb = false;
1392 1352
1393 for(j = 1; j < fAngleBin; j++) 1353 for(j = 1; j < fAngleBin; j++)
1394 { 1354 {
1395 // alpha1 = angleBins->GetLowEdgeEnergy 1355 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1396 // alpha2 = angleBins->GetLowEdgeEnergy 1356 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1397 1357
1398 alpha1 = alphaMax*(j-1)/fAngleBin; 1358 alpha1 = alphaMax*(j-1)/fAngleBin;
1399 alpha2 = alphaMax*( j )/fAngleBin; 1359 alpha2 = alphaMax*( j )/fAngleBin;
1400 1360
1401 if( ( alpha2 > alphaCoulomb ) && z ) fAdd 1361 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1402 1362
1403 deltaL10 = integral.Legendre10(this, &G4N 1363 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1404 deltaL96 = integral.Legendre96(this, &G4N 1364 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1405 deltaAG = integral.AdaptiveGauss(this, & 1365 deltaAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1406 alpha1 1366 alpha1, alpha2,epsilon);
1407 1367
1408 // G4cout<<alpha1<<"\t"<<std::sqrt(alph 1368 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1409 // <<deltaL10<<"\t"<<deltaL96<<"\t" 1369 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1410 1370
1411 sumL10 += deltaL10; 1371 sumL10 += deltaL10;
1412 sumL96 += deltaL96; 1372 sumL96 += deltaL96;
1413 sumAG += deltaAG; 1373 sumAG += deltaAG;
1414 1374
1415 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d 1375 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1416 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1376 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1417 1377
1418 angleVector->PutValue( j-1 , alpha1, sumL 1378 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1419 } 1379 }
1420 fAngleTable->insertAt(i,angleVector); 1380 fAngleTable->insertAt(i,angleVector);
1421 fAngleBank.push_back(fAngleTable); 1381 fAngleBank.push_back(fAngleTable);
1422 1382
1423 /* 1383 /*
1424 // Integral over all angle range - Bad accu 1384 // Integral over all angle range - Bad accuracy !!!
1425 1385
1426 sumL10 = integral.Legendre10(this, &G4NuclN 1386 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1427 sumL96 = integral.Legendre96(this, &G4NuclN 1387 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1428 sumAG = integral.AdaptiveGauss(this, &G4Nu 1388 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1429 0., al 1389 0., alpha2,epsilon);
1430 G4cout<<G4endl; 1390 G4cout<<G4endl;
1431 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg 1391 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1432 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1392 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1433 */ 1393 */
1434 return; 1394 return;
1435 } 1395 }
1436 1396
1437 ///////////////////////////////////////////// 1397 /////////////////////////////////////////////////////////////////
1438 // 1398 //
1439 // 1399 //
1440 1400
1441 G4double G4NuclNuclDiffuseElastic::GetLegend 1401 G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
1442 { 1402 {
1443 G4double legPol, epsilon = 1.e-6; 1403 G4double legPol, epsilon = 1.e-6;
1444 G4double x = std::cos(theta); 1404 G4double x = std::cos(theta);
1445 1405
1446 if ( n < 0 ) legPol = 0.; 1406 if ( n < 0 ) legPol = 0.;
1447 else if( n == 0 ) legPol = 1.; 1407 else if( n == 0 ) legPol = 1.;
1448 else if( n == 1 ) legPol = x; 1408 else if( n == 1 ) legPol = x;
1449 else if( n == 2 ) legPol = (3.*x*x-1.)/2.; 1409 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1450 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/ 1410 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
1451 else if( n == 4 ) legPol = (35.*x*x*x*x-30. 1411 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1452 else if( n == 5 ) legPol = (63.*x*x*x*x*x-7 1412 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
1453 else if( n == 6 ) legPol = (231.*x*x*x*x*x* 1413 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1454 else 1414 else
1455 { 1415 {
1456 // legPol = ( (2*n-1)*x*GetLegendrePol(n- 1416 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
1457 1417
1458 legPol = std::sqrt( 2./(n*CLHEP::pi*std:: 1418 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
1459 } 1419 }
1460 return legPol; 1420 return legPol;
1461 } 1421 }
1462 1422
1463 ///////////////////////////////////////////// 1423 /////////////////////////////////////////////////////////////////
1464 // 1424 //
1465 // 1425 //
1466 1426
1467 G4complex G4NuclNuclDiffuseElastic::GetErfCom 1427 G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
1468 { 1428 {
1469 G4int n; 1429 G4int n;
1470 G4double n2, cofn, shny, chny, fn, gn; 1430 G4double n2, cofn, shny, chny, fn, gn;
1471 1431
1472 G4double x = z.real(); 1432 G4double x = z.real();
1473 G4double y = z.imag(); 1433 G4double y = z.imag();
1474 1434
1475 G4double outRe = 0., outIm = 0.; 1435 G4double outRe = 0., outIm = 0.;
1476 1436
1477 G4double twox = 2.*x; 1437 G4double twox = 2.*x;
1478 G4double twoxy = twox*y; 1438 G4double twoxy = twox*y;
1479 G4double twox2 = twox*twox; 1439 G4double twox2 = twox*twox;
1480 1440
1481 G4double cof1 = G4Exp(-x*x)/CLHEP::pi; 1441 G4double cof1 = G4Exp(-x*x)/CLHEP::pi;
1482 1442
1483 G4double cos2xy = std::cos(twoxy); 1443 G4double cos2xy = std::cos(twoxy);
1484 G4double sin2xy = std::sin(twoxy); 1444 G4double sin2xy = std::sin(twoxy);
1485 1445
1486 G4double twoxcos2xy = twox*cos2xy; 1446 G4double twoxcos2xy = twox*cos2xy;
1487 G4double twoxsin2xy = twox*sin2xy; 1447 G4double twoxsin2xy = twox*sin2xy;
1488 1448
1489 for( n = 1; n <= nMax; n++) 1449 for( n = 1; n <= nMax; n++)
1490 { 1450 {
1491 n2 = n*n; 1451 n2 = n*n;
1492 1452
1493 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n 1453 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
1494 1454
1495 chny = std::cosh(n*y); 1455 chny = std::cosh(n*y);
1496 shny = std::sinh(n*y); 1456 shny = std::sinh(n*y);
1497 1457
1498 fn = twox - twoxcos2xy*chny + n*sin2xy*s 1458 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1499 gn = twoxsin2xy*chny + n*cos2xy*s 1459 gn = twoxsin2xy*chny + n*cos2xy*shny;
1500 1460
1501 fn *= cofn; 1461 fn *= cofn;
1502 gn *= cofn; 1462 gn *= cofn;
1503 1463
1504 outRe += fn; 1464 outRe += fn;
1505 outIm += gn; 1465 outIm += gn;
1506 } 1466 }
1507 outRe *= 2*cof1; 1467 outRe *= 2*cof1;
1508 outIm *= 2*cof1; 1468 outIm *= 2*cof1;
1509 1469
1510 if(std::abs(x) < 0.0001) 1470 if(std::abs(x) < 0.0001)
1511 { 1471 {
1512 outRe += GetErf(x); 1472 outRe += GetErf(x);
1513 outIm += cof1*y; 1473 outIm += cof1*y;
1514 } 1474 }
1515 else 1475 else
1516 { 1476 {
1517 outRe += GetErf(x) + cof1*(1-cos2xy)/twox 1477 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
1518 outIm += cof1*sin2xy/twox; 1478 outIm += cof1*sin2xy/twox;
1519 } 1479 }
1520 return G4complex(outRe, outIm); 1480 return G4complex(outRe, outIm);
1521 } 1481 }
1522 1482
1523 1483
1524 ///////////////////////////////////////////// 1484 /////////////////////////////////////////////////////////////////
1525 // 1485 //
1526 // 1486 //
1527 1487
1528 G4complex G4NuclNuclDiffuseElastic::GetErfInt 1488 G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
1529 { 1489 {
1530 G4double outRe, outIm; 1490 G4double outRe, outIm;
1531 1491
1532 G4double x = z.real(); 1492 G4double x = z.real();
1533 G4double y = z.imag(); 1493 G4double y = z.imag();
1534 fReZ = x; 1494 fReZ = x;
1535 1495
1536 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1496 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1537 1497
1538 outRe = integral.Legendre96(this,&G4NuclNuc 1498 outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
1539 outIm = integral.Legendre96(this,&G4NuclNuc 1499 outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
1540 1500
1541 outRe *= 2./std::sqrt(CLHEP::pi); 1501 outRe *= 2./std::sqrt(CLHEP::pi);
1542 outIm *= 2./std::sqrt(CLHEP::pi); 1502 outIm *= 2./std::sqrt(CLHEP::pi);
1543 1503
1544 outRe += GetErf(x); 1504 outRe += GetErf(x);
1545 1505
1546 return G4complex(outRe, outIm); 1506 return G4complex(outRe, outIm);
1547 } 1507 }
1548 1508
1549 ///////////////////////////////////////////// 1509 /////////////////////////////////////////////////////////////////
1550 // 1510 //
1551 // 1511 //
1552 1512
1553 1513
1554 G4complex G4NuclNuclDiffuseElastic::GammaLess 1514 G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
1555 { 1515 {
1556 G4double sinThetaR = 2.*fHalfRutThetaT 1516 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1557 G4double cosHalfThetaR2 = 1./(1. + fHalfRut 1517 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1558 1518
1559 G4double u = std::sqrt(0.5*fPr 1519 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1560 G4double kappa = u/std::sqrt(CLHEP 1520 G4double kappa = u/std::sqrt(CLHEP::pi);
1561 G4double dTheta = theta - fRutherfo 1521 G4double dTheta = theta - fRutherfordTheta;
1562 u *= dTheta; 1522 u *= dTheta;
1563 G4double u2 = u*u; 1523 G4double u2 = u*u;
1564 G4double u2m2p3 = u2*2./3.; 1524 G4double u2m2p3 = u2*2./3.;
1565 1525
1566 G4complex im = G4complex(0.,1.); 1526 G4complex im = G4complex(0.,1.);
1567 G4complex order = G4complex(u,u); 1527 G4complex order = G4complex(u,u);
1568 order /= std::sqrt(2.); 1528 order /= std::sqrt(2.);
1569 1529
1570 G4complex gamma = CLHEP::pi*kappa*G 1530 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1571 G4complex a0 = 0.5*(1. + 4.*(1.+ 1531 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1572 G4complex a1 = 0.5*(1. + 2.*(1.+ 1532 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1573 G4complex out = gamma*(1. - a1*dT 1533 G4complex out = gamma*(1. - a1*dTheta) - a0;
1574 1534
1575 return out; 1535 return out;
1576 } 1536 }
1577 1537
1578 ///////////////////////////////////////////// 1538 /////////////////////////////////////////////////////////////////
1579 // 1539 //
1580 // 1540 //
1581 1541
1582 G4complex G4NuclNuclDiffuseElastic::GammaMore 1542 G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
1583 { 1543 {
1584 G4double sinThetaR = 2.*fHalfRutThetaT 1544 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1585 G4double cosHalfThetaR2 = 1./(1. + fHalfRut 1545 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1586 1546
1587 G4double u = std::sqrt(0.5*fPr 1547 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1588 G4double kappa = u/std::sqrt(CLHEP 1548 G4double kappa = u/std::sqrt(CLHEP::pi);
1589 G4double dTheta = theta - fRutherfo 1549 G4double dTheta = theta - fRutherfordTheta;
1590 u *= dTheta; 1550 u *= dTheta;
1591 G4double u2 = u*u; 1551 G4double u2 = u*u;
1592 G4double u2m2p3 = u2*2./3.; 1552 G4double u2m2p3 = u2*2./3.;
1593 1553
1594 G4complex im = G4complex(0.,1.); 1554 G4complex im = G4complex(0.,1.);
1595 G4complex order = G4complex(u,u); 1555 G4complex order = G4complex(u,u);
1596 order /= std::sqrt(2.); 1556 order /= std::sqrt(2.);
1597 G4complex gamma = CLHEP::pi*kappa*G 1557 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1598 G4complex a0 = 0.5*(1. + 4.*(1.+ 1558 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1599 G4complex a1 = 0.5*(1. + 2.*(1.+ 1559 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1600 G4complex out = -gamma*(1. - a1*d 1560 G4complex out = -gamma*(1. - a1*dTheta) - a0;
1601 1561
1602 return out; 1562 return out;
1603 } 1563 }
1604 1564
1605 ///////////////////////////////////////////// 1565 /////////////////////////////////////////////////////////////////
1606 // 1566 //
1607 // 1567 //
1608 1568
1609 G4complex G4NuclNuclDiffuseElastic::Amplitud 1569 G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
1610 { 1570 {
1611 G4double kappa = std::sqrt(0.5*fProfileLamb 1571 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1612 G4complex out = G4complex(kappa/fWaveVector 1572 G4complex out = G4complex(kappa/fWaveVector,0.);
1613 1573
1614 out *= PhaseNear(theta); 1574 out *= PhaseNear(theta);
1615 1575
1616 if( theta <= fRutherfordTheta ) 1576 if( theta <= fRutherfordTheta )
1617 { 1577 {
1618 out *= GammaLess(theta) + ProfileNear(the 1578 out *= GammaLess(theta) + ProfileNear(theta);
1619 // out *= GammaMore(theta) + ProfileNear( 1579 // out *= GammaMore(theta) + ProfileNear(theta);
1620 out += CoulombAmplitude(theta); 1580 out += CoulombAmplitude(theta);
1621 } 1581 }
1622 else 1582 else
1623 { 1583 {
1624 out *= GammaMore(theta) + ProfileNear(the 1584 out *= GammaMore(theta) + ProfileNear(theta);
1625 // out *= GammaLess(theta) + ProfileNear( 1585 // out *= GammaLess(theta) + ProfileNear(theta);
1626 } 1586 }
1627 return out; 1587 return out;
1628 } 1588 }
1629 1589
1630 ///////////////////////////////////////////// 1590 /////////////////////////////////////////////////////////////////
1631 // 1591 //
1632 // 1592 //
1633 1593
1634 G4complex G4NuclNuclDiffuseElastic::Amplitude 1594 G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
1635 { 1595 {
1636 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 1596 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1637 G4double dTheta = 0.5*(theta - fRutherf 1597 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1638 G4double sindTheta = std::sin(dTheta); 1598 G4double sindTheta = std::sin(dTheta);
1639 G4double persqrt2 = std::sqrt(0.5); 1599 G4double persqrt2 = std::sqrt(0.5);
1640 1600
1641 G4complex order = G4complex(persqrt2,pe 1601 G4complex order = G4complex(persqrt2,persqrt2);
1642 order *= std::sqrt(0.5*fProfil 1602 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1643 // order *= std::sqrt(0.5*fPro 1603 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1644 1604
1645 G4complex out; 1605 G4complex out;
1646 1606
1647 if ( theta <= fRutherfordTheta ) 1607 if ( theta <= fRutherfordTheta )
1648 { 1608 {
1649 out = 1. - 0.5*GetErfcInt(-order)*Profile 1609 out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1650 } 1610 }
1651 else 1611 else
1652 { 1612 {
1653 out = 0.5*GetErfcInt(order)*ProfileNear(t 1613 out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1654 } 1614 }
1655 1615
1656 out *= CoulombAmplitude(theta); 1616 out *= CoulombAmplitude(theta);
1657 1617
1658 return out; 1618 return out;
1659 } 1619 }
1660 1620
1661 ///////////////////////////////////////////// 1621 /////////////////////////////////////////////////////////////////
1662 // 1622 //
1663 // 1623 //
1664 1624
1665 G4complex G4NuclNuclDiffuseElastic::Amplitud 1625 G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
1666 { 1626 {
1667 G4int n; 1627 G4int n;
1668 G4double T12b, b, b2; // cosTheta = std::co 1628 G4double T12b, b, b2; // cosTheta = std::cos(theta);
1669 G4complex out = G4complex(0.,0.), shiftC, s 1629 G4complex out = G4complex(0.,0.), shiftC, shiftN;
1670 G4complex im = G4complex(0.,1.); 1630 G4complex im = G4complex(0.,1.);
1671 1631
1672 for( n = 0; n < fMaxL; n++) 1632 for( n = 0; n < fMaxL; n++)
1673 { 1633 {
1674 shiftC = std::exp( im*2.*CalculateCoulomb 1634 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1675 // b = ( fZommerfeld + std::sqrt( fZommer 1635 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1676 b = ( std::sqrt( G4double(n*(n+1)) ) )/fW 1636 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1677 b2 = b*b; 1637 b2 = b*b;
1678 T12b = fSumSigma*G4Exp(-b2/fNuclearRadius 1638 T12b = fSumSigma*G4Exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1679 shiftN = std::exp( -0.5*(1.-im*fEtaRatio) 1639 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1680 out += (2.*n+1.)*shiftC*shiftN*GetLegend 1640 out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1681 } 1641 }
1682 out /= 2.*im*fWaveVector; 1642 out /= 2.*im*fWaveVector;
1683 out += CoulombAmplitude(theta); 1643 out += CoulombAmplitude(theta);
1684 return out; 1644 return out;
1685 } 1645 }
1686 1646
1687 1647
1688 ///////////////////////////////////////////// 1648 /////////////////////////////////////////////////////////////////
1689 // 1649 //
1690 // 1650 //
1691 1651
1692 G4complex G4NuclNuclDiffuseElastic::Amplitud 1652 G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
1693 { 1653 {
1694 G4int n; 1654 G4int n;
1695 G4double T12b, a, aTemp, b2, sinThetaH = st 1655 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1696 G4double sinThetaH2 = sinThetaH*sinThetaH; 1656 G4double sinThetaH2 = sinThetaH*sinThetaH;
1697 G4complex out = G4complex(0.,0.); 1657 G4complex out = G4complex(0.,0.);
1698 G4complex im = G4complex(0.,1.); 1658 G4complex im = G4complex(0.,1.);
1699 1659
1700 a = -fSumSigma/CLHEP::twopi/fNuclearRadius 1660 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1701 b2 = fWaveVector*fWaveVector*fNuclearRadius 1661 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1702 1662
1703 aTemp = a; 1663 aTemp = a;
1704 1664
1705 for( n = 1; n < fMaxL; n++) 1665 for( n = 1; n < fMaxL; n++)
1706 { 1666 {
1707 T12b = aTemp*G4Exp(-b2/n)/n; 1667 T12b = aTemp*G4Exp(-b2/n)/n;
1708 aTemp *= a; 1668 aTemp *= a;
1709 out += T12b; 1669 out += T12b;
1710 G4cout<<"out = "<<out<<G4endl; 1670 G4cout<<"out = "<<out<<G4endl;
1711 } 1671 }
1712 out *= -4.*im*fWaveVector/CLHEP::pi; 1672 out *= -4.*im*fWaveVector/CLHEP::pi;
1713 out += CoulombAmplitude(theta); 1673 out += CoulombAmplitude(theta);
1714 return out; 1674 return out;
1715 } 1675 }
1716 1676
1717 1677
1718 ///////////////////////////////////////////// 1678 ///////////////////////////////////////////////////////////////////////////////
1719 // 1679 //
1720 // Test for given particle and element table 1680 // Test for given particle and element table of momentum, angle probability.
1721 // For the partMom in CMS. 1681 // For the partMom in CMS.
1722 1682
1723 void G4NuclNuclDiffuseElastic::InitParameters 1683 void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,
1724 G4d 1684 G4double partMom, G4double Z, G4double A)
1725 { 1685 {
1726 fAtomicNumber = Z; // atomic number 1686 fAtomicNumber = Z; // atomic number
1727 fAtomicWeight = A; // number of nucleo 1687 fAtomicWeight = A; // number of nucleons
1728 1688
1729 fNuclearRadius2 = CalculateNuclearRad(fAtom 1689 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1730 G4double A1 = G4double( theParticle->Ge 1690 G4double A1 = G4double( theParticle->GetBaryonNumber() );
1731 fNuclearRadius1 = CalculateNuclearRad(A1); 1691 fNuclearRadius1 = CalculateNuclearRad(A1);
1732 // fNuclearRadius = std::sqrt(fNuclearRadiu 1692 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1733 fNuclearRadius = fNuclearRadius1 + fNuclear 1693 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1734 1694
1735 G4double a = 0.; 1695 G4double a = 0.;
1736 G4double z = theParticle->GetPDGCharge(); 1696 G4double z = theParticle->GetPDGCharge();
1737 G4double m1 = theParticle->GetPDGMass(); 1697 G4double m1 = theParticle->GetPDGMass();
1738 1698
1739 fWaveVector = partMom/CLHEP::hbarc; 1699 fWaveVector = partMom/CLHEP::hbarc;
1740 1700
1741 G4double lambda = fCofLambda*fWaveVector*fN 1701 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1742 G4cout<<"kR = "<<lambda<<G4endl; 1702 G4cout<<"kR = "<<lambda<<G4endl;
1743 1703
1744 if( z ) 1704 if( z )
1745 { 1705 {
1746 a = partMom/m1; // beta*gamma f 1706 a = partMom/m1; // beta*gamma for m1
1747 fBeta = a/std::sqrt(1+a*a); 1707 fBeta = a/std::sqrt(1+a*a);
1748 fZommerfeld = CalculateZommerfeld( fBeta, 1708 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1749 fRutherfordRatio = fZommerfeld/fWaveVecto 1709 fRutherfordRatio = fZommerfeld/fWaveVector;
1750 fAm = CalculateAm( partMom, fZomm 1710 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1751 } 1711 }
1752 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4en 1712 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1753 fProfileLambda = lambda; // *std::sqrt(1.-2 1713 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1754 G4cout<<"fProfileLambda = "<<fProfileLambda 1714 G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1755 fProfileDelta = fCofDelta*fProfileLambda; 1715 fProfileDelta = fCofDelta*fProfileLambda;
1756 fProfileAlpha = fCofAlpha*fProfileLambda; 1716 fProfileAlpha = fCofAlpha*fProfileLambda;
1757 1717
1758 CalculateCoulombPhaseZero(); 1718 CalculateCoulombPhaseZero();
1759 CalculateRutherfordAnglePar(); 1719 CalculateRutherfordAnglePar();
1760 1720
1761 return; 1721 return;
1762 } 1722 }
1763 ///////////////////////////////////////////// 1723 ///////////////////////////////////////////////////////////////////////////////
1764 // 1724 //
1765 // Test for given particle and element table 1725 // Test for given particle and element table of momentum, angle probability.
1766 // For the partMom in CMS. 1726 // For the partMom in CMS.
1767 1727
1768 void G4NuclNuclDiffuseElastic::InitDynParamet 1728 void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,
1769 G4d 1729 G4double partMom)
1770 { 1730 {
1771 G4double a = 0.; 1731 G4double a = 0.;
1772 G4double z = theParticle->GetPDGCharge(); 1732 G4double z = theParticle->GetPDGCharge();
1773 G4double m1 = theParticle->GetPDGMass(); 1733 G4double m1 = theParticle->GetPDGMass();
1774 1734
1775 fWaveVector = partMom/CLHEP::hbarc; 1735 fWaveVector = partMom/CLHEP::hbarc;
1776 1736
1777 G4double lambda = fCofLambda*fWaveVector*fN 1737 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1778 1738
1779 if( z ) 1739 if( z )
1780 { 1740 {
1781 a = partMom/m1; // beta*gamma f 1741 a = partMom/m1; // beta*gamma for m1
1782 fBeta = a/std::sqrt(1+a*a); 1742 fBeta = a/std::sqrt(1+a*a);
1783 fZommerfeld = CalculateZommerfeld( fBeta, 1743 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1784 fRutherfordRatio = fZommerfeld/fWaveVecto 1744 fRutherfordRatio = fZommerfeld/fWaveVector;
1785 fAm = CalculateAm( partMom, fZomm 1745 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1786 } 1746 }
1787 fProfileLambda = lambda; // *std::sqrt(1.-2 1747 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1788 fProfileDelta = fCofDelta*fProfileLambda; 1748 fProfileDelta = fCofDelta*fProfileLambda;
1789 fProfileAlpha = fCofAlpha*fProfileLambda; 1749 fProfileAlpha = fCofAlpha*fProfileLambda;
1790 1750
1791 CalculateCoulombPhaseZero(); 1751 CalculateCoulombPhaseZero();
1792 CalculateRutherfordAnglePar(); 1752 CalculateRutherfordAnglePar();
1793 1753
1794 return; 1754 return;
1795 } 1755 }
1796 1756
1797 1757
1798 ///////////////////////////////////////////// 1758 ///////////////////////////////////////////////////////////////////////////////
1799 // 1759 //
1800 // Test for given particle and element table 1760 // Test for given particle and element table of momentum, angle probability.
1801 // For the partMom in CMS. 1761 // For the partMom in CMS.
1802 1762
1803 void G4NuclNuclDiffuseElastic::InitParameters 1763 void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,
1804 G4d 1764 G4double partMom, G4double Z, G4double A)
1805 { 1765 {
1806 fAtomicNumber = Z; // target atomic nu 1766 fAtomicNumber = Z; // target atomic number
1807 fAtomicWeight = A; // target number of 1767 fAtomicWeight = A; // target number of nucleons
1808 1768
1809 fNuclearRadius2 = CalculateNuclearRad(fAtom 1769 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1810 G4double A1 = G4double( aParticle->GetD 1770 G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1811 fNuclearRadius1 = CalculateNuclearRad(A1); 1771 fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1812 fNuclearRadiusSquare = fNuclearRadius1*fNuc 1772 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1813 1773
1814 1774
1815 G4double a = 0., kR12; 1775 G4double a = 0., kR12;
1816 G4double z = aParticle->GetDefinition()->G 1776 G4double z = aParticle->GetDefinition()->GetPDGCharge();
1817 G4double m1 = aParticle->GetDefinition()->G 1777 G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1818 1778
1819 fWaveVector = partMom/CLHEP::hbarc; 1779 fWaveVector = partMom/CLHEP::hbarc;
1820 1780
1821 G4double pN = A1 - z; 1781 G4double pN = A1 - z;
1822 if( pN < 0. ) pN = 0.; 1782 if( pN < 0. ) pN = 0.;
1823 1783
1824 G4double tN = A - Z; 1784 G4double tN = A - Z;
1825 if( tN < 0. ) tN = 0.; 1785 if( tN < 0. ) tN = 0.;
1826 1786
1827 G4double pTkin = aParticle->GetKineticEnerg 1787 G4double pTkin = aParticle->GetKineticEnergy();
1828 pTkin /= A1; 1788 pTkin /= A1;
1829 1789
1830 1790
1831 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXsc 1791 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1832 (z*tN+pN*Z)*GetHadronNucleonXsc 1792 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1833 1793
1834 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::mi 1794 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1835 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiu 1795 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
1836 kR12 = fWaveVector*std::sqrt(fNuclearRadius 1796 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1837 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl; 1797 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1838 fMaxL = (G4int(kR12)+1)*4; 1798 fMaxL = (G4int(kR12)+1)*4;
1839 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl; 1799 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1840 1800
1841 if( z ) 1801 if( z )
1842 { 1802 {
1843 a = partMom/m1; // beta*gamma 1803 a = partMom/m1; // beta*gamma for m1
1844 fBeta = a/std::sqrt(1+a*a); 1804 fBeta = a/std::sqrt(1+a*a);
1845 fZommerfeld = CalculateZommerfeld( fBet 1805 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1846 fAm = CalculateAm( partMom, fZo 1806 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1847 } 1807 }
1848 1808
1849 CalculateCoulombPhaseZero(); 1809 CalculateCoulombPhaseZero();
1850 1810
1851 1811
1852 return; 1812 return;
1853 } 1813 }
1854 1814
1855 1815
1856 ///////////////////////////////////////////// 1816 /////////////////////////////////////////////////////////////////////////////////////
1857 // 1817 //
1858 // Returns nucleon-nucleon cross-section base 1818 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1859 // data from mainly http://wwwppds.ihep.su:80 1819 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1860 // projectile nucleon is pParticle with pTkin 1820 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1861 1821
1862 G4double 1822 G4double
1863 G4NuclNuclDiffuseElastic::GetHadronNucleonXsc 1823 G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
1864 1824 G4double pTkin,
1865 1825 G4ParticleDefinition* tParticle)
1866 { 1826 {
1867 G4double xsection(0), /*Delta,*/ A0, B0; 1827 G4double xsection(0), /*Delta,*/ A0, B0;
1868 G4double hpXsc(0); 1828 G4double hpXsc(0);
1869 G4double hnXsc(0); 1829 G4double hnXsc(0);
1870 1830
1871 1831
1872 G4double targ_mass = tParticle->GetPDGM 1832 G4double targ_mass = tParticle->GetPDGMass();
1873 G4double proj_mass = pParticle->GetPDGM 1833 G4double proj_mass = pParticle->GetPDGMass();
1874 1834
1875 G4double proj_energy = proj_mass + pTkin; 1835 G4double proj_energy = proj_mass + pTkin;
1876 G4double proj_momentum = std::sqrt(pTkin*(p 1836 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1877 1837
1878 G4double sMand = CalcMandelstamS ( proj_mas 1838 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1879 1839
1880 sMand /= CLHEP::GeV*CLHEP::GeV; // 1840 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1881 proj_momentum /= CLHEP::GeV; 1841 proj_momentum /= CLHEP::GeV;
1882 proj_energy /= CLHEP::GeV; 1842 proj_energy /= CLHEP::GeV;
1883 proj_mass /= CLHEP::GeV; 1843 proj_mass /= CLHEP::GeV;
1884 G4double logS = G4Log(sMand); 1844 G4double logS = G4Log(sMand);
1885 1845
1886 // General PDG fit constants 1846 // General PDG fit constants
1887 1847
1888 1848
1889 // fEtaRatio=Re[f(0)]/Im[f(0)] 1849 // fEtaRatio=Re[f(0)]/Im[f(0)]
1890 1850
1891 if( proj_momentum >= 1.2 ) 1851 if( proj_momentum >= 1.2 )
1892 { 1852 {
1893 fEtaRatio = 0.13*(logS - 5.8579332)*G4Po 1853 fEtaRatio = 0.13*(logS - 5.8579332)*G4Pow::GetInstance()->powA(sMand,-0.18);
1894 } 1854 }
1895 else if( proj_momentum >= 0.6 ) 1855 else if( proj_momentum >= 0.6 )
1896 { 1856 {
1897 fEtaRatio = -75.5*(G4Pow::GetInstance()-> 1857 fEtaRatio = -75.5*(G4Pow::GetInstance()->powA(proj_momentum,0.25)-0.95)/
1898 (G4Pow::GetInstance()->powA(3*proj_moment 1858 (G4Pow::GetInstance()->powA(3*proj_momentum,2.2)+1);
1899 } 1859 }
1900 else 1860 else
1901 { 1861 {
1902 fEtaRatio = 15.5*proj_momentum/(27*proj_m 1862 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1903 } 1863 }
1904 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl; 1864 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1905 1865
1906 // xsc 1866 // xsc
1907 1867
1908 if( proj_momentum >= 10. ) // high energy: 1868 if( proj_momentum >= 10. ) // high energy: pp = nn = np
1909 // if( proj_momentum >= 2.) 1869 // if( proj_momentum >= 2.)
1910 { 1870 {
1911 //Delta = 1.; 1871 //Delta = 1.;
1912 1872
1913 //if( proj_energy < 40. ) Delta = 0.916+0 1873 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1914 1874
1915 //AR-12Aug2016 if( proj_momentum >= 10.) 1875 //AR-12Aug2016 if( proj_momentum >= 10.)
1916 { 1876 {
1917 B0 = 7.5; 1877 B0 = 7.5;
1918 A0 = 100. - B0*G4Log(3.0e7); 1878 A0 = 100. - B0*G4Log(3.0e7);
1919 1879
1920 xsection = A0 + B0*G4Log(proj_energy) 1880 xsection = A0 + B0*G4Log(proj_energy) - 11
1921 + 103*G4Pow::GetInstance()- 1881 + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
1922 0.93827*0.93827,-0.165); 1882 0.93827*0.93827,-0.165); // mb
1923 } 1883 }
1924 } 1884 }
1925 else // low energy pp = nn != np 1885 else // low energy pp = nn != np
1926 { 1886 {
1927 if(pParticle == tParticle) // pp or nn 1887 if(pParticle == tParticle) // pp or nn // nn to be pp
1928 { 1888 {
1929 if( proj_momentum < 0.73 ) 1889 if( proj_momentum < 0.73 )
1930 { 1890 {
1931 hnXsc = 23 + 50*( G4Pow::GetInstanc 1891 hnXsc = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
1932 } 1892 }
1933 else if( proj_momentum < 1.05 ) 1893 else if( proj_momentum < 1.05 )
1934 { 1894 {
1935 hnXsc = 23 + 40*(G4Log(proj_momentu 1895 hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))*
1936 (G4Log(proj_momentum 1896 (G4Log(proj_momentum/0.73));
1937 } 1897 }
1938 else // if( proj_momentum < 10. ) 1898 else // if( proj_momentum < 10. )
1939 { 1899 {
1940 hnXsc = 39.0 + 1900 hnXsc = 39.0 +
1941 75*(proj_momentum - 1.2)/(G4Pow 1901 75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
1942 } 1902 }
1943 xsection = hnXsc; 1903 xsection = hnXsc;
1944 } 1904 }
1945 else // pn to be np 1905 else // pn to be np
1946 { 1906 {
1947 if( proj_momentum < 0.8 ) 1907 if( proj_momentum < 0.8 )
1948 { 1908 {
1949 hpXsc = 33+30*G4Pow::GetInstance()- 1909 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
1950 } 1910 }
1951 else if( proj_momentum < 1.4 ) 1911 else if( proj_momentum < 1.4 )
1952 { 1912 {
1953 hpXsc = 33+30*G4Pow::GetInstance()- 1913 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
1954 } 1914 }
1955 else // if( proj_momentum < 10. ) 1915 else // if( proj_momentum < 10. )
1956 { 1916 {
1957 hpXsc = 33.3+ 1917 hpXsc = 33.3+
1958 20.8*(G4Pow::GetInstance()->pow 1918 20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
1959 (G4Pow::GetInstance()->powA( 1919 (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
1960 } 1920 }
1961 xsection = hpXsc; 1921 xsection = hpXsc;
1962 } 1922 }
1963 } 1923 }
1964 xsection *= CLHEP::millibarn; // parametris 1924 xsection *= CLHEP::millibarn; // parametrised in mb
1965 G4cout<<"xsection = "<<xsection/CLHEP::mill 1925 G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1966 return xsection; 1926 return xsection;
1967 } 1927 }
1968 1928
1969 ///////////////////////////////////////////// 1929 /////////////////////////////////////////////////////////////////
1970 // 1930 //
1971 // The ratio el/ruth for Fresnel smooth nucle 1931 // The ratio el/ruth for Fresnel smooth nucleus profile
1972 1932
1973 G4double G4NuclNuclDiffuseElastic::GetRatioGe 1933 G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
1974 { 1934 {
1975 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 1935 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1976 G4double dTheta = 0.5*(theta - fRutherf 1936 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1977 G4double sindTheta = std::sin(dTheta); 1937 G4double sindTheta = std::sin(dTheta);
1978 1938
1979 G4double prof = Profile(theta); 1939 G4double prof = Profile(theta);
1980 G4double prof2 = prof*prof; 1940 G4double prof2 = prof*prof;
1981 // G4double profmod = std::abs(prof); 1941 // G4double profmod = std::abs(prof);
1982 G4double order = std::sqrt(fProfileLam 1942 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1983 1943
1984 order = std::abs(order); // since sin chang 1944 order = std::abs(order); // since sin changes sign!
1985 // G4cout<<"order = "<<order<<G4endl; 1945 // G4cout<<"order = "<<order<<G4endl;
1986 1946
1987 G4double cosFresnel = GetCint(order); 1947 G4double cosFresnel = GetCint(order);
1988 G4double sinFresnel = GetSint(order); 1948 G4double sinFresnel = GetSint(order);
1989 1949
1990 G4double out; 1950 G4double out;
1991 1951
1992 if ( theta <= fRutherfordTheta ) 1952 if ( theta <= fRutherfordTheta )
1993 { 1953 {
1994 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-c 1954 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1995 out += ( cosFresnel + sinFresnel - 1. )*p 1955 out += ( cosFresnel + sinFresnel - 1. )*prof;
1996 } 1956 }
1997 else 1957 else
1998 { 1958 {
1999 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFres 1959 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
2000 } 1960 }
2001 1961
2002 return out; 1962 return out;
2003 } 1963 }
2004 1964
2005 ///////////////////////////////////////////// 1965 ///////////////////////////////////////////////////////////////////
2006 // 1966 //
2007 // For the calculation of arg Gamma(z) one ne 1967 // For the calculation of arg Gamma(z) one needs complex extension
2008 // of ln(Gamma(z)) 1968 // of ln(Gamma(z))
2009 1969
2010 G4complex G4NuclNuclDiffuseElastic::GammaLoga 1970 G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
2011 { 1971 {
2012 const G4double cof[6] = { 76.18009172947146 1972 const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
2013 24.0140982408309 1973 24.01409824083091, -1.231739572450155,
2014 0.1208650973866 1974 0.1208650973866179e-2, -0.5395239384953e-5 } ;
2015 G4int j; 1975 G4int j;
2016 G4complex z = zz - 1.0; 1976 G4complex z = zz - 1.0;
2017 G4complex tmp = z + 5.5; 1977 G4complex tmp = z + 5.5;
2018 tmp -= (z + 0.5) * std::log(tmp); 1978 tmp -= (z + 0.5) * std::log(tmp);
2019 G4complex ser = G4complex(1.000000000190015 1979 G4complex ser = G4complex(1.000000000190015,0.);
2020 1980
2021 for ( j = 0; j <= 5; j++ ) 1981 for ( j = 0; j <= 5; j++ )
2022 { 1982 {
2023 z += 1.0; 1983 z += 1.0;
2024 ser += cof[j]/z; 1984 ser += cof[j]/z;
2025 } 1985 }
2026 return -tmp + std::log(2.5066282746310005*s 1986 return -tmp + std::log(2.5066282746310005*ser);
2027 } 1987 }
2028 1988
2029 ///////////////////////////////////////////// 1989 /////////////////////////////////////////////////////////////
2030 // 1990 //
2031 // Bessel J0 function based on rational appro 1991 // Bessel J0 function based on rational approximation from
2032 // J.F. Hart, Computer Approximations, New Yo 1992 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2033 1993
2034 G4double G4NuclNuclDiffuseElastic::BesselJzer 1994 G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
2035 { 1995 {
2036 G4double modvalue, value2, fact1, fact2, ar 1996 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2037 1997
2038 modvalue = std::fabs(value); 1998 modvalue = std::fabs(value);
2039 1999
2040 if ( value < 8.0 && value > -8.0 ) 2000 if ( value < 8.0 && value > -8.0 )
2041 { 2001 {
2042 value2 = value*value; 2002 value2 = value*value;
2043 2003
2044 fact1 = 57568490574.0 + value2*(-1336259 2004 fact1 = 57568490574.0 + value2*(-13362590354.0
2045 + value2*( 6516196 2005 + value2*( 651619640.7
2046 + value2*(-1121442 2006 + value2*(-11214424.18
2047 + value2*( 77392.3 2007 + value2*( 77392.33017
2048 + value2*(-184.905 2008 + value2*(-184.9052456 ) ) ) ) );
2049 2009
2050 fact2 = 57568490411.0 + value2*( 1029532 2010 fact2 = 57568490411.0 + value2*( 1029532985.0
2051 + value2*( 9494680 2011 + value2*( 9494680.718
2052 + value2*(59272.64 2012 + value2*(59272.64853
2053 + value2*(267.8532 2013 + value2*(267.8532712
2054 + value2*1.0 2014 + value2*1.0 ) ) ) );
2055 2015
2056 bessel = fact1/fact2; 2016 bessel = fact1/fact2;
2057 } 2017 }
2058 else 2018 else
2059 { 2019 {
2060 arg = 8.0/modvalue; 2020 arg = 8.0/modvalue;
2061 2021
2062 value2 = arg*arg; 2022 value2 = arg*arg;
2063 2023
2064 shift = modvalue-0.785398164; 2024 shift = modvalue-0.785398164;
2065 2025
2066 fact1 = 1.0 + value2*(-0.1098628627e-2 2026 fact1 = 1.0 + value2*(-0.1098628627e-2
2067 + value2*(0.2734510407e-4 2027 + value2*(0.2734510407e-4
2068 + value2*(-0.2073370639e-5 2028 + value2*(-0.2073370639e-5
2069 + value2*0.2093887211e-6 2029 + value2*0.2093887211e-6 ) ) );
2070 2030
2071 fact2 = -0.1562499995e-1 + value2*(0.143 2031 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2072 + value2*(-0.69 2032 + value2*(-0.6911147651e-5
2073 + value2*(0.762 2033 + value2*(0.7621095161e-6
2074 - value2*0.9349 2034 - value2*0.934945152e-7 ) ) );
2075 2035
2076 bessel = std::sqrt(0.636619772/modvalue)* 2036 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2077 } 2037 }
2078 return bessel; 2038 return bessel;
2079 } 2039 }
2080 2040
2081 ///////////////////////////////////////////// 2041 /////////////////////////////////////////////////////////////
2082 // 2042 //
2083 // Bessel J1 function based on rational appro 2043 // Bessel J1 function based on rational approximation from
2084 // J.F. Hart, Computer Approximations, New Yo 2044 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2085 2045
2086 G4double G4NuclNuclDiffuseElastic::BesselJone 2046 G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
2087 { 2047 {
2088 G4double modvalue, value2, fact1, fact2, ar 2048 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2089 2049
2090 modvalue = std::fabs(value); 2050 modvalue = std::fabs(value);
2091 2051
2092 if ( modvalue < 8.0 ) 2052 if ( modvalue < 8.0 )
2093 { 2053 {
2094 value2 = value*value; 2054 value2 = value*value;
2095 2055
2096 fact1 = value*(72362614232.0 + value2*(- 2056 fact1 = value*(72362614232.0 + value2*(-7895059235.0
2097 + value2*( 2057 + value2*( 242396853.1
2098 + value2*(- 2058 + value2*(-2972611.439
2099 + value2*( 2059 + value2*( 15704.48260
2100 + value2*(- 2060 + value2*(-30.16036606 ) ) ) ) ) );
2101 2061
2102 fact2 = 144725228442.0 + value2*(2300535 2062 fact2 = 144725228442.0 + value2*(2300535178.0
2103 + value2*(1858330 2063 + value2*(18583304.74
2104 + value2*(99447.4 2064 + value2*(99447.43394
2105 + value2*(376.999 2065 + value2*(376.9991397
2106 + value2*1.0 2066 + value2*1.0 ) ) ) );
2107 bessel = fact1/fact2; 2067 bessel = fact1/fact2;
2108 } 2068 }
2109 else 2069 else
2110 { 2070 {
2111 arg = 8.0/modvalue; 2071 arg = 8.0/modvalue;
2112 2072
2113 value2 = arg*arg; 2073 value2 = arg*arg;
2114 2074
2115 shift = modvalue - 2.356194491; 2075 shift = modvalue - 2.356194491;
2116 2076
2117 fact1 = 1.0 + value2*( 0.183105e-2 2077 fact1 = 1.0 + value2*( 0.183105e-2
2118 + value2*(-0.3516396496e-4 2078 + value2*(-0.3516396496e-4
2119 + value2*(0.2457520174e-5 2079 + value2*(0.2457520174e-5
2120 + value2*(-0.240337019e-6 2080 + value2*(-0.240337019e-6 ) ) ) );
2121 2081
2122 fact2 = 0.04687499995 + value2*(-0.200269 2082 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2123 + value2*( 0.844919 2083 + value2*( 0.8449199096e-5
2124 + value2*(-0.882289 2084 + value2*(-0.88228987e-6
2125 + value2*0.10578741 2085 + value2*0.105787412e-6 ) ) );
2126 2086
2127 bessel = std::sqrt( 0.636619772/modvalue) 2087 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2128 2088
2129 if (value < 0.0) bessel = -bessel; 2089 if (value < 0.0) bessel = -bessel;
2130 } 2090 }
2131 return bessel; 2091 return bessel;
2132 } 2092 }
2133 2093
2134 // 2094 //
2135 // 2095 //
2136 ///////////////////////////////////////////// 2096 /////////////////////////////////////////////////////////////////////////////////
2137 2097