| 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 // 26 //
27 // 27 //
28 // Physics model class G4DiffuseElastic 28 // Physics model class G4DiffuseElastic
29 // 29 //
30 // 30 //
31 // G4 Model: optical diffuse elastic scatterin 31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
32 // 32 //
33 // 24-May-07 V. Grichine 33 // 24-May-07 V. Grichine
34 // 34 //
35 // 21.10.15 V. Grichine 35 // 21.10.15 V. Grichine
36 // Bug fixed in BuildAngleTable, i 36 // Bug fixed in BuildAngleTable, improving accuracy for
37 // angle bins at high energies > 5 37 // angle bins at high energies > 50 GeV for pions.
38 // 38 //
39 39
40 #include "G4DiffuseElastic.hh" 40 #include "G4DiffuseElastic.hh"
41 #include "G4ParticleTable.hh" 41 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh" 42 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh" 43 #include "G4IonTable.hh"
44 #include "G4NucleiProperties.hh" 44 #include "G4NucleiProperties.hh"
45 45
46 #include "Randomize.hh" 46 #include "Randomize.hh"
47 #include "G4Integrator.hh" 47 #include "G4Integrator.hh"
48 #include "globals.hh" 48 #include "globals.hh"
49 #include "G4PhysicalConstants.hh" 49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
51 51
52 #include "G4Proton.hh" 52 #include "G4Proton.hh"
53 #include "G4Neutron.hh" 53 #include "G4Neutron.hh"
54 #include "G4Deuteron.hh" 54 #include "G4Deuteron.hh"
55 #include "G4Alpha.hh" 55 #include "G4Alpha.hh"
56 #include "G4PionPlus.hh" 56 #include "G4PionPlus.hh"
57 #include "G4PionMinus.hh" 57 #include "G4PionMinus.hh"
58 58
59 #include "G4Element.hh" 59 #include "G4Element.hh"
60 #include "G4ElementTable.hh" 60 #include "G4ElementTable.hh"
61 #include "G4NistManager.hh" 61 #include "G4NistManager.hh"
62 #include "G4PhysicsTable.hh" 62 #include "G4PhysicsTable.hh"
63 #include "G4PhysicsLogVector.hh" 63 #include "G4PhysicsLogVector.hh"
64 #include "G4PhysicsFreeVector.hh" 64 #include "G4PhysicsFreeVector.hh"
65 65
66 #include "G4Exp.hh" 66 #include "G4Exp.hh"
67 67
68 #include "G4HadronicParameters.hh" 68 #include "G4HadronicParameters.hh"
69 69
70 ////////////////////////////////////////////// 70 /////////////////////////////////////////////////////////////////////////
71 // 71 //
72 // Test Constructor. Just to check xsc 72 // Test Constructor. Just to check xsc
73 73
74 74
75 G4DiffuseElastic::G4DiffuseElastic() 75 G4DiffuseElastic::G4DiffuseElastic()
76 : G4HadronElastic("DiffuseElastic"), fPartic 76 : G4HadronElastic("DiffuseElastic"), fParticle(0)
77 { 77 {
78 SetMinEnergy( 0.01*MeV ); 78 SetMinEnergy( 0.01*MeV );
79 SetMaxEnergy( G4HadronicParameters::Instance 79 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
80 80
81 verboseLevel = 0; 81 verboseLevel = 0;
82 lowEnergyRecoilLimit = 100.*keV; 82 lowEnergyRecoilLimit = 100.*keV;
83 lowEnergyLimitQ = 0.0*GeV; 83 lowEnergyLimitQ = 0.0*GeV;
84 lowEnergyLimitHE = 0.0*GeV; 84 lowEnergyLimitHE = 0.0*GeV;
85 lowestEnergyLimit = 0.0*keV; 85 lowestEnergyLimit = 0.0*keV;
86 plabLowLimit = 20.0*MeV; 86 plabLowLimit = 20.0*MeV;
87 87
88 theProton = G4Proton::Proton(); 88 theProton = G4Proton::Proton();
89 theNeutron = G4Neutron::Neutron(); 89 theNeutron = G4Neutron::Neutron();
90 theDeuteron = G4Deuteron::Deuteron(); 90 theDeuteron = G4Deuteron::Deuteron();
91 theAlpha = G4Alpha::Alpha(); 91 theAlpha = G4Alpha::Alpha();
92 thePionPlus = G4PionPlus::PionPlus(); 92 thePionPlus = G4PionPlus::PionPlus();
93 thePionMinus = G4PionMinus::PionMinus(); 93 thePionMinus = G4PionMinus::PionMinus();
94 94
95 fEnergyBin = 300; // Increased from the ori 95 fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV
96 fAngleBin = 200; 96 fAngleBin = 200;
97 97
98 fEnergyVector = new G4PhysicsLogVector( the 98 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
99 99
100 fAngleTable = 0; 100 fAngleTable = 0;
101 101
102 fParticle = 0; 102 fParticle = 0;
103 fWaveVector = 0.; 103 fWaveVector = 0.;
104 fAtomicWeight = 0.; 104 fAtomicWeight = 0.;
105 fAtomicNumber = 0.; 105 fAtomicNumber = 0.;
106 fNuclearRadius = 0.; 106 fNuclearRadius = 0.;
107 fBeta = 0.; 107 fBeta = 0.;
108 fZommerfeld = 0.; 108 fZommerfeld = 0.;
109 fAm = 0.; 109 fAm = 0.;
110 fAddCoulomb = false; 110 fAddCoulomb = false;
111 } 111 }
112 112
113 ////////////////////////////////////////////// 113 //////////////////////////////////////////////////////////////////////////////
114 // 114 //
115 // Destructor 115 // Destructor
116 116
117 G4DiffuseElastic::~G4DiffuseElastic() 117 G4DiffuseElastic::~G4DiffuseElastic()
118 { 118 {
119 if ( fEnergyVector ) 119 if ( fEnergyVector )
120 { 120 {
121 delete fEnergyVector; 121 delete fEnergyVector;
122 fEnergyVector = 0; 122 fEnergyVector = 0;
123 } 123 }
124 for ( std::vector<G4PhysicsTable*>::iterator 124 for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
125 it != fAngleBank.end(); ++it ) 125 it != fAngleBank.end(); ++it )
126 { 126 {
127 if ( (*it) ) (*it)->clearAndDestroy(); 127 if ( (*it) ) (*it)->clearAndDestroy();
128 128
129 delete *it; 129 delete *it;
130 *it = 0; 130 *it = 0;
131 } 131 }
132 fAngleTable = 0; 132 fAngleTable = 0;
133 } 133 }
134 134
135 ////////////////////////////////////////////// 135 //////////////////////////////////////////////////////////////////////////////
136 // 136 //
137 // Initialisation for given particle using ele 137 // Initialisation for given particle using element table of application
138 138
139 void G4DiffuseElastic::Initialise() 139 void G4DiffuseElastic::Initialise()
140 { 140 {
141 141
142 // fEnergyVector = new G4PhysicsLogVector( t 142 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
143 143
144 const G4ElementTable* theElementTable = G4El 144 const G4ElementTable* theElementTable = G4Element::GetElementTable();
145 145
146 std::size_t jEl, numOfEl = G4Element::GetNum << 146 size_t jEl, numOfEl = G4Element::GetNumberOfElements();
147 147
148 for( jEl = 0; jEl < numOfEl; ++jEl) // appli 148 for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop
149 { 149 {
150 fAtomicNumber = (*theElementTable)[jEl]->G 150 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
151 fAtomicWeight = G4NistManager::Instance()- 151 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
152 fNuclearRadius = CalculateNuclearRad(fAtom 152 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
153 153
154 if( verboseLevel > 0 ) 154 if( verboseLevel > 0 )
155 { 155 {
156 G4cout<<"G4DiffuseElastic::Initialise() 156 G4cout<<"G4DiffuseElastic::Initialise() the element: "
157 <<(*theElementTable)[jEl]->GetName()<<G4 157 <<(*theElementTable)[jEl]->GetName()<<G4endl;
158 } 158 }
159 fElementNumberVector.push_back(fAtomicNumb 159 fElementNumberVector.push_back(fAtomicNumber);
160 fElementNameVector.push_back((*theElementT 160 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
161 161
162 BuildAngleTable(); 162 BuildAngleTable();
163 fAngleBank.push_back(fAngleTable); 163 fAngleBank.push_back(fAngleTable);
164 } 164 }
165 return; 165 return;
166 } 166 }
167 167
168 ////////////////////////////////////////////// 168 ////////////////////////////////////////////////////////////////////////////
169 // 169 //
170 // return differential elastic cross section d 170 // return differential elastic cross section d(sigma)/d(omega)
171 171
172 G4double 172 G4double
173 G4DiffuseElastic::GetDiffuseElasticXsc( const 173 G4DiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
174 G4doub 174 G4double theta,
175 G4double momentum, 175 G4double momentum,
176 G4doub 176 G4double A )
177 { 177 {
178 fParticle = particle; 178 fParticle = particle;
179 fWaveVector = momentum/hbarc; 179 fWaveVector = momentum/hbarc;
180 fAtomicWeight = A; 180 fAtomicWeight = A;
181 fAddCoulomb = false; 181 fAddCoulomb = false;
182 fNuclearRadius = CalculateNuclearRad(A); 182 fNuclearRadius = CalculateNuclearRad(A);
183 183
184 G4double sigma = fNuclearRadius*fNuclearRadi 184 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
185 185
186 return sigma; 186 return sigma;
187 } 187 }
188 188
189 ////////////////////////////////////////////// 189 ////////////////////////////////////////////////////////////////////////////
190 // 190 //
191 // return invariant differential elastic cross 191 // return invariant differential elastic cross section d(sigma)/d(tMand)
192 192
193 G4double 193 G4double
194 G4DiffuseElastic::GetInvElasticXsc( const G4Pa 194 G4DiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle,
195 G4doub 195 G4double tMand,
196 G4double plab, 196 G4double plab,
197 G4doub 197 G4double A, G4double Z )
198 { 198 {
199 G4double m1 = particle->GetPDGMass(); 199 G4double m1 = particle->GetPDGMass();
200 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 200 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
201 201
202 G4int iZ = static_cast<G4int>(Z+0.5); 202 G4int iZ = static_cast<G4int>(Z+0.5);
203 G4int iA = static_cast<G4int>(A+0.5); 203 G4int iA = static_cast<G4int>(A+0.5);
204 G4ParticleDefinition * theDef = 0; 204 G4ParticleDefinition * theDef = 0;
205 205
206 if (iZ == 1 && iA == 1) theDef = thePro 206 if (iZ == 1 && iA == 1) theDef = theProton;
207 else if (iZ == 1 && iA == 2) theDef = theDeu 207 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
208 else if (iZ == 1 && iA == 3) theDef = G4Trit 208 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
209 else if (iZ == 2 && iA == 3) theDef = G4He3: 209 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
210 else if (iZ == 2 && iA == 4) theDef = theAlp 210 else if (iZ == 2 && iA == 4) theDef = theAlpha;
211 else theDef = G4ParticleTable::GetParticleTa 211 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
212 212
213 G4double tmass = theDef->GetPDGMass(); 213 G4double tmass = theDef->GetPDGMass();
214 214
215 G4LorentzVector lv(0.0,0.0,0.0,tmass); 215 G4LorentzVector lv(0.0,0.0,0.0,tmass);
216 lv += lv1; 216 lv += lv1;
217 217
218 G4ThreeVector bst = lv.boostVector(); 218 G4ThreeVector bst = lv.boostVector();
219 lv1.boost(-bst); 219 lv1.boost(-bst);
220 220
221 G4ThreeVector p1 = lv1.vect(); 221 G4ThreeVector p1 = lv1.vect();
222 G4double ptot = p1.mag(); 222 G4double ptot = p1.mag();
223 G4double ptot2 = ptot*ptot; 223 G4double ptot2 = ptot*ptot;
224 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 224 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
225 225
226 if( cost >= 1.0 ) cost = 1.0; 226 if( cost >= 1.0 ) cost = 1.0;
227 else if( cost <= -1.0) cost = -1.0; 227 else if( cost <= -1.0) cost = -1.0;
228 228
229 G4double thetaCMS = std::acos(cost); 229 G4double thetaCMS = std::acos(cost);
230 230
231 G4double sigma = GetDiffuseElasticXsc( parti 231 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
232 232
233 sigma *= pi/ptot2; 233 sigma *= pi/ptot2;
234 234
235 return sigma; 235 return sigma;
236 } 236 }
237 237
238 ////////////////////////////////////////////// 238 ////////////////////////////////////////////////////////////////////////////
239 // 239 //
240 // return differential elastic cross section d 240 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
241 // correction 241 // correction
242 242
243 G4double 243 G4double
244 G4DiffuseElastic::GetDiffuseElasticSumXsc( con 244 G4DiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
245 G4doub 245 G4double theta,
246 G4double momentum, 246 G4double momentum,
247 G4doub 247 G4double A, G4double Z )
248 { 248 {
249 fParticle = particle; 249 fParticle = particle;
250 fWaveVector = momentum/hbarc; 250 fWaveVector = momentum/hbarc;
251 fAtomicWeight = A; 251 fAtomicWeight = A;
252 fAtomicNumber = Z; 252 fAtomicNumber = Z;
253 fNuclearRadius = CalculateNuclearRad(A); 253 fNuclearRadius = CalculateNuclearRad(A);
254 fAddCoulomb = false; 254 fAddCoulomb = false;
255 255
256 G4double z = particle->GetPDGCharge(); 256 G4double z = particle->GetPDGCharge();
257 257
258 G4double kRt = fWaveVector*fNuclearRadius* 258 G4double kRt = fWaveVector*fNuclearRadius*theta;
259 G4double kRtC = 1.9; 259 G4double kRtC = 1.9;
260 260
261 if( z && (kRt > kRtC) ) 261 if( z && (kRt > kRtC) )
262 { 262 {
263 fAddCoulomb = true; 263 fAddCoulomb = true;
264 fBeta = CalculateParticleBeta( parti 264 fBeta = CalculateParticleBeta( particle, momentum);
265 fZommerfeld = CalculateZommerfeld( fBeta, 265 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
266 fAm = CalculateAm( momentum, fZomm 266 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
267 } 267 }
268 G4double sigma = fNuclearRadius*fNuclearRadi 268 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
269 269
270 return sigma; 270 return sigma;
271 } 271 }
272 272
273 ////////////////////////////////////////////// 273 ////////////////////////////////////////////////////////////////////////////
274 // 274 //
275 // return invariant differential elastic cross 275 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
276 // correction 276 // correction
277 277
278 G4double 278 G4double
279 G4DiffuseElastic::GetInvElasticSumXsc( const G 279 G4DiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle,
280 G4doub 280 G4double tMand,
281 G4double plab, 281 G4double plab,
282 G4doub 282 G4double A, G4double Z )
283 { 283 {
284 G4double m1 = particle->GetPDGMass(); 284 G4double m1 = particle->GetPDGMass();
285 285
286 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 286 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
287 287
288 G4int iZ = static_cast<G4int>(Z+0.5); 288 G4int iZ = static_cast<G4int>(Z+0.5);
289 G4int iA = static_cast<G4int>(A+0.5); 289 G4int iA = static_cast<G4int>(A+0.5);
290 290
291 G4ParticleDefinition* theDef = 0; 291 G4ParticleDefinition* theDef = 0;
292 292
293 if (iZ == 1 && iA == 1) theDef = thePro 293 if (iZ == 1 && iA == 1) theDef = theProton;
294 else if (iZ == 1 && iA == 2) theDef = theDeu 294 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
295 else if (iZ == 1 && iA == 3) theDef = G4Trit 295 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
296 else if (iZ == 2 && iA == 3) theDef = G4He3: 296 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
297 else if (iZ == 2 && iA == 4) theDef = theAlp 297 else if (iZ == 2 && iA == 4) theDef = theAlpha;
298 else theDef = G4ParticleTable::GetParticleTa 298 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
299 299
300 G4double tmass = theDef->GetPDGMass(); 300 G4double tmass = theDef->GetPDGMass();
301 301
302 G4LorentzVector lv(0.0,0.0,0.0,tmass); 302 G4LorentzVector lv(0.0,0.0,0.0,tmass);
303 lv += lv1; 303 lv += lv1;
304 304
305 G4ThreeVector bst = lv.boostVector(); 305 G4ThreeVector bst = lv.boostVector();
306 lv1.boost(-bst); 306 lv1.boost(-bst);
307 307
308 G4ThreeVector p1 = lv1.vect(); 308 G4ThreeVector p1 = lv1.vect();
309 G4double ptot = p1.mag(); 309 G4double ptot = p1.mag();
310 G4double ptot2 = ptot*ptot; 310 G4double ptot2 = ptot*ptot;
311 G4double cost = 1 - 0.5*std::fabs(tMand)/ 311 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
312 312
313 if( cost >= 1.0 ) cost = 1.0; 313 if( cost >= 1.0 ) cost = 1.0;
314 else if( cost <= -1.0) cost = -1.0; 314 else if( cost <= -1.0) cost = -1.0;
315 315
316 G4double thetaCMS = std::acos(cost); 316 G4double thetaCMS = std::acos(cost);
317 317
318 G4double sigma = GetDiffuseElasticSumXsc( pa 318 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
319 319
320 sigma *= pi/ptot2; 320 sigma *= pi/ptot2;
321 321
322 return sigma; 322 return sigma;
323 } 323 }
324 324
325 ////////////////////////////////////////////// 325 ////////////////////////////////////////////////////////////////////////////
326 // 326 //
327 // return invariant differential elastic cross 327 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
328 // correction 328 // correction
329 329
330 G4double 330 G4double
331 G4DiffuseElastic::GetInvCoulombElasticXsc( con 331 G4DiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
332 G4doub 332 G4double tMand,
333 G4double plab, 333 G4double plab,
334 G4doub 334 G4double A, G4double Z )
335 { 335 {
336 G4double m1 = particle->GetPDGMass(); 336 G4double m1 = particle->GetPDGMass();
337 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 337 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
338 338
339 G4int iZ = static_cast<G4int>(Z+0.5); 339 G4int iZ = static_cast<G4int>(Z+0.5);
340 G4int iA = static_cast<G4int>(A+0.5); 340 G4int iA = static_cast<G4int>(A+0.5);
341 G4ParticleDefinition * theDef = 0; 341 G4ParticleDefinition * theDef = 0;
342 342
343 if (iZ == 1 && iA == 1) theDef = thePro 343 if (iZ == 1 && iA == 1) theDef = theProton;
344 else if (iZ == 1 && iA == 2) theDef = theDeu 344 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
345 else if (iZ == 1 && iA == 3) theDef = G4Trit 345 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
346 else if (iZ == 2 && iA == 3) theDef = G4He3: 346 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
347 else if (iZ == 2 && iA == 4) theDef = theAlp 347 else if (iZ == 2 && iA == 4) theDef = theAlpha;
348 else theDef = G4ParticleTable::GetParticleTa 348 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
349 349
350 G4double tmass = theDef->GetPDGMass(); 350 G4double tmass = theDef->GetPDGMass();
351 351
352 G4LorentzVector lv(0.0,0.0,0.0,tmass); 352 G4LorentzVector lv(0.0,0.0,0.0,tmass);
353 lv += lv1; 353 lv += lv1;
354 354
355 G4ThreeVector bst = lv.boostVector(); 355 G4ThreeVector bst = lv.boostVector();
356 lv1.boost(-bst); 356 lv1.boost(-bst);
357 357
358 G4ThreeVector p1 = lv1.vect(); 358 G4ThreeVector p1 = lv1.vect();
359 G4double ptot = p1.mag(); 359 G4double ptot = p1.mag();
360 G4double ptot2 = ptot*ptot; 360 G4double ptot2 = ptot*ptot;
361 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 361 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
362 362
363 if( cost >= 1.0 ) cost = 1.0; 363 if( cost >= 1.0 ) cost = 1.0;
364 else if( cost <= -1.0) cost = -1.0; 364 else if( cost <= -1.0) cost = -1.0;
365 365
366 G4double thetaCMS = std::acos(cost); 366 G4double thetaCMS = std::acos(cost);
367 367
368 G4double sigma = GetCoulombElasticXsc( parti 368 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
369 369
370 sigma *= pi/ptot2; 370 sigma *= pi/ptot2;
371 371
372 return sigma; 372 return sigma;
373 } 373 }
374 374
375 ////////////////////////////////////////////// 375 ////////////////////////////////////////////////////////////////////////////
376 // 376 //
377 // return differential elastic probability d(p 377 // return differential elastic probability d(probability)/d(omega)
378 378
379 G4double 379 G4double
380 G4DiffuseElastic::GetDiffElasticProb( // G4Par 380 G4DiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
381 G4doub 381 G4double theta
382 // G4double momentum, 382 // G4double momentum,
383 // G4double A 383 // G4double A
384 ) 384 )
385 { 385 {
386 G4double sigma, bzero, bzero2, bonebyarg, bo 386 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
387 G4double delta, diffuse, gamma; 387 G4double delta, diffuse, gamma;
388 G4double e1, e2, bone, bone2; 388 G4double e1, e2, bone, bone2;
389 389
390 // G4double wavek = momentum/hbarc; // wave 390 // G4double wavek = momentum/hbarc; // wave vector
391 // G4double r0 = 1.08*fermi; 391 // G4double r0 = 1.08*fermi;
392 // G4double rad = r0*G4Pow::GetInstance()- 392 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
393 393
394 if (fParticle == theProton) 394 if (fParticle == theProton)
395 { 395 {
396 diffuse = 0.63*fermi; 396 diffuse = 0.63*fermi;
397 gamma = 0.3*fermi; 397 gamma = 0.3*fermi;
398 delta = 0.1*fermi*fermi; 398 delta = 0.1*fermi*fermi;
399 e1 = 0.3*fermi; 399 e1 = 0.3*fermi;
400 e2 = 0.35*fermi; 400 e2 = 0.35*fermi;
401 } 401 }
402 else if (fParticle == theNeutron) 402 else if (fParticle == theNeutron)
403 { 403 {
404 diffuse = 0.63*fermi; // 1.63*fermi; // 404 diffuse = 0.63*fermi; // 1.63*fermi; //
405 G4double k0 = 1*GeV/hbarc; 405 G4double k0 = 1*GeV/hbarc;
406 diffuse *= k0/fWaveVector; 406 diffuse *= k0/fWaveVector;
407 407
408 gamma = 0.3*fermi; 408 gamma = 0.3*fermi;
409 delta = 0.1*fermi*fermi; 409 delta = 0.1*fermi*fermi;
410 e1 = 0.3*fermi; 410 e1 = 0.3*fermi;
411 e2 = 0.35*fermi; 411 e2 = 0.35*fermi;
412 } 412 }
413 else // as proton, if were not defined 413 else // as proton, if were not defined
414 { 414 {
415 diffuse = 0.63*fermi; 415 diffuse = 0.63*fermi;
416 gamma = 0.3*fermi; 416 gamma = 0.3*fermi;
417 delta = 0.1*fermi*fermi; 417 delta = 0.1*fermi*fermi;
418 e1 = 0.3*fermi; 418 e1 = 0.3*fermi;
419 e2 = 0.35*fermi; 419 e2 = 0.35*fermi;
420 } 420 }
421 G4double kr = fWaveVector*fNuclearRadius; 421 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
422 G4double kr2 = kr*kr; 422 G4double kr2 = kr*kr;
423 G4double krt = kr*theta; 423 G4double krt = kr*theta;
424 424
425 bzero = BesselJzero(krt); 425 bzero = BesselJzero(krt);
426 bzero2 = bzero*bzero; 426 bzero2 = bzero*bzero;
427 bone = BesselJone(krt); 427 bone = BesselJone(krt);
428 bone2 = bone*bone; 428 bone2 = bone*bone;
429 bonebyarg = BesselOneByArg(krt); 429 bonebyarg = BesselOneByArg(krt);
430 bonebyarg2 = bonebyarg*bonebyarg; 430 bonebyarg2 = bonebyarg*bonebyarg;
431 431
432 G4double lambda = 15.; // 15 ok 432 G4double lambda = 15.; // 15 ok
433 433
434 // G4double kgamma = fWaveVector*gamma; 434 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
435 435
436 G4double kgamma = lambda*(1.-G4Exp(-fWave 436 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
437 G4double kgamma2 = kgamma*kgamma; 437 G4double kgamma2 = kgamma*kgamma;
438 438
439 // G4double dk2t = delta*fWaveVector*fWaveV 439 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
440 // G4double dk2t2 = dk2t*dk2t; 440 // G4double dk2t2 = dk2t*dk2t;
441 // G4double pikdt = pi*fWaveVector*diffuse*t 441 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
442 442
443 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 443 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
444 444
445 damp = DampFactor(pikdt); 445 damp = DampFactor(pikdt);
446 damp2 = damp*damp; 446 damp2 = damp*damp;
447 447
448 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 448 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
449 G4double e2dk3t = -2.*e2*delta*fWaveVector* 449 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
450 450
451 451
452 sigma = kgamma2; 452 sigma = kgamma2;
453 // sigma += dk2t2; 453 // sigma += dk2t2;
454 sigma *= bzero2; 454 sigma *= bzero2;
455 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 455 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
456 sigma += kr2*bonebyarg2; 456 sigma += kr2*bonebyarg2;
457 sigma *= damp2; // *rad*rad; 457 sigma *= damp2; // *rad*rad;
458 458
459 return sigma; 459 return sigma;
460 } 460 }
461 461
462 ////////////////////////////////////////////// 462 ////////////////////////////////////////////////////////////////////////////
463 // 463 //
464 // return differential elastic probability d(p 464 // return differential elastic probability d(probability)/d(omega) with
465 // Coulomb correction 465 // Coulomb correction
466 466
467 G4double 467 G4double
468 G4DiffuseElastic::GetDiffElasticSumProb( // G4 468 G4DiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
469 G4doub 469 G4double theta
470 // G4double momentum, 470 // G4double momentum,
471 // G4double A 471 // G4double A
472 ) 472 )
473 { 473 {
474 G4double sigma, bzero, bzero2, bonebyarg, bo 474 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
475 G4double delta, diffuse, gamma; 475 G4double delta, diffuse, gamma;
476 G4double e1, e2, bone, bone2; 476 G4double e1, e2, bone, bone2;
477 477
478 // G4double wavek = momentum/hbarc; // wave 478 // G4double wavek = momentum/hbarc; // wave vector
479 // G4double r0 = 1.08*fermi; 479 // G4double r0 = 1.08*fermi;
480 // G4double rad = r0*G4Pow::GetInstance()- 480 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
481 481
482 G4double kr = fWaveVector*fNuclearRadius; 482 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
483 G4double kr2 = kr*kr; 483 G4double kr2 = kr*kr;
484 G4double krt = kr*theta; 484 G4double krt = kr*theta;
485 485
486 bzero = BesselJzero(krt); 486 bzero = BesselJzero(krt);
487 bzero2 = bzero*bzero; 487 bzero2 = bzero*bzero;
488 bone = BesselJone(krt); 488 bone = BesselJone(krt);
489 bone2 = bone*bone; 489 bone2 = bone*bone;
490 bonebyarg = BesselOneByArg(krt); 490 bonebyarg = BesselOneByArg(krt);
491 bonebyarg2 = bonebyarg*bonebyarg; 491 bonebyarg2 = bonebyarg*bonebyarg;
492 492
493 if (fParticle == theProton) 493 if (fParticle == theProton)
494 { 494 {
495 diffuse = 0.63*fermi; 495 diffuse = 0.63*fermi;
496 // diffuse = 0.6*fermi; 496 // diffuse = 0.6*fermi;
497 gamma = 0.3*fermi; 497 gamma = 0.3*fermi;
498 delta = 0.1*fermi*fermi; 498 delta = 0.1*fermi*fermi;
499 e1 = 0.3*fermi; 499 e1 = 0.3*fermi;
500 e2 = 0.35*fermi; 500 e2 = 0.35*fermi;
501 } 501 }
502 else if (fParticle == theNeutron) 502 else if (fParticle == theNeutron)
503 { 503 {
504 diffuse = 0.63*fermi; 504 diffuse = 0.63*fermi;
505 // diffuse = 0.6*fermi; 505 // diffuse = 0.6*fermi;
506 G4double k0 = 1*GeV/hbarc; 506 G4double k0 = 1*GeV/hbarc;
507 diffuse *= k0/fWaveVector; 507 diffuse *= k0/fWaveVector;
508 gamma = 0.3*fermi; 508 gamma = 0.3*fermi;
509 delta = 0.1*fermi*fermi; 509 delta = 0.1*fermi*fermi;
510 e1 = 0.3*fermi; 510 e1 = 0.3*fermi;
511 e2 = 0.35*fermi; 511 e2 = 0.35*fermi;
512 } 512 }
513 else // as proton, if were not defined 513 else // as proton, if were not defined
514 { 514 {
515 diffuse = 0.63*fermi; 515 diffuse = 0.63*fermi;
516 gamma = 0.3*fermi; 516 gamma = 0.3*fermi;
517 delta = 0.1*fermi*fermi; 517 delta = 0.1*fermi*fermi;
518 e1 = 0.3*fermi; 518 e1 = 0.3*fermi;
519 e2 = 0.35*fermi; 519 e2 = 0.35*fermi;
520 } 520 }
521 G4double lambda = 15.; // 15 ok 521 G4double lambda = 15.; // 15 ok
522 // G4double kgamma = fWaveVector*gamma; 522 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
523 G4double kgamma = lambda*(1.-G4Exp(-fWave 523 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
524 524
525 // G4cout<<"kgamma = "<<kgamma<<G4endl; 525 // G4cout<<"kgamma = "<<kgamma<<G4endl;
526 526
527 if(fAddCoulomb) // add Coulomb correction 527 if(fAddCoulomb) // add Coulomb correction
528 { 528 {
529 G4double sinHalfTheta = std::sin(0.5*thet 529 G4double sinHalfTheta = std::sin(0.5*theta);
530 G4double sinHalfTheta2 = sinHalfTheta*sinH 530 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
531 531
532 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 532 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
533 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 533 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
534 } 534 }
535 535
536 G4double kgamma2 = kgamma*kgamma; 536 G4double kgamma2 = kgamma*kgamma;
537 537
538 // G4double dk2t = delta*fWaveVector*fWaveV 538 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
539 // G4cout<<"dk2t = "<<dk2t<<G4endl; 539 // G4cout<<"dk2t = "<<dk2t<<G4endl;
540 // G4double dk2t2 = dk2t*dk2t; 540 // G4double dk2t2 = dk2t*dk2t;
541 // G4double pikdt = pi*fWaveVector*diffuse*t 541 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
542 542
543 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 543 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
544 544
545 // G4cout<<"pikdt = "<<pikdt<<G4endl; 545 // G4cout<<"pikdt = "<<pikdt<<G4endl;
546 546
547 damp = DampFactor(pikdt); 547 damp = DampFactor(pikdt);
548 damp2 = damp*damp; 548 damp2 = damp*damp;
549 549
550 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 550 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
551 G4double e2dk3t = -2.*e2*delta*fWaveVector* 551 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
552 552
553 sigma = kgamma2; 553 sigma = kgamma2;
554 // sigma += dk2t2; 554 // sigma += dk2t2;
555 sigma *= bzero2; 555 sigma *= bzero2;
556 sigma += mode2k2*bone2; 556 sigma += mode2k2*bone2;
557 sigma += e2dk3t*bzero*bone; 557 sigma += e2dk3t*bzero*bone;
558 558
559 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 559 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
560 sigma += kr2*bonebyarg2; // correction at J 560 sigma += kr2*bonebyarg2; // correction at J1()/()
561 561
562 sigma *= damp2; // *rad*rad; 562 sigma *= damp2; // *rad*rad;
563 563
564 return sigma; 564 return sigma;
565 } 565 }
566 566
567 567
568 ////////////////////////////////////////////// 568 ////////////////////////////////////////////////////////////////////////////
569 // 569 //
570 // return differential elastic probability d(p 570 // return differential elastic probability d(probability)/d(t) with
571 // Coulomb correction. It is called from Build 571 // Coulomb correction. It is called from BuildAngleTable()
572 572
573 G4double 573 G4double
574 G4DiffuseElastic::GetDiffElasticSumProbA( G4do 574 G4DiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
575 { 575 {
576 G4double theta; 576 G4double theta;
577 577
578 theta = std::sqrt(alpha); 578 theta = std::sqrt(alpha);
579 579
580 // theta = std::acos( 1 - alpha/2. ); 580 // theta = std::acos( 1 - alpha/2. );
581 581
582 G4double sigma, bzero, bzero2, bonebyarg, bo 582 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
583 G4double delta, diffuse, gamma; 583 G4double delta, diffuse, gamma;
584 G4double e1, e2, bone, bone2; 584 G4double e1, e2, bone, bone2;
585 585
586 // G4double wavek = momentum/hbarc; // wave 586 // G4double wavek = momentum/hbarc; // wave vector
587 // G4double r0 = 1.08*fermi; 587 // G4double r0 = 1.08*fermi;
588 // G4double rad = r0*G4Pow::GetInstance()- 588 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
589 589
590 G4double kr = fWaveVector*fNuclearRadius; 590 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
591 G4double kr2 = kr*kr; 591 G4double kr2 = kr*kr;
592 G4double krt = kr*theta; 592 G4double krt = kr*theta;
593 593
594 bzero = BesselJzero(krt); 594 bzero = BesselJzero(krt);
595 bzero2 = bzero*bzero; 595 bzero2 = bzero*bzero;
596 bone = BesselJone(krt); 596 bone = BesselJone(krt);
597 bone2 = bone*bone; 597 bone2 = bone*bone;
598 bonebyarg = BesselOneByArg(krt); 598 bonebyarg = BesselOneByArg(krt);
599 bonebyarg2 = bonebyarg*bonebyarg; 599 bonebyarg2 = bonebyarg*bonebyarg;
600 600
601 if ( fParticle == theProton ) 601 if ( fParticle == theProton )
602 { 602 {
603 diffuse = 0.63*fermi; 603 diffuse = 0.63*fermi;
604 // diffuse = 0.6*fermi; 604 // diffuse = 0.6*fermi;
605 gamma = 0.3*fermi; 605 gamma = 0.3*fermi;
606 delta = 0.1*fermi*fermi; 606 delta = 0.1*fermi*fermi;
607 e1 = 0.3*fermi; 607 e1 = 0.3*fermi;
608 e2 = 0.35*fermi; 608 e2 = 0.35*fermi;
609 } 609 }
610 else if ( fParticle == theNeutron ) 610 else if ( fParticle == theNeutron )
611 { 611 {
612 diffuse = 0.63*fermi; 612 diffuse = 0.63*fermi;
613 // diffuse = 0.6*fermi; 613 // diffuse = 0.6*fermi;
614 // G4double k0 = 0.8*GeV/hbarc; 614 // G4double k0 = 0.8*GeV/hbarc;
615 // diffuse *= k0/fWaveVector; 615 // diffuse *= k0/fWaveVector;
616 gamma = 0.3*fermi; 616 gamma = 0.3*fermi;
617 delta = 0.1*fermi*fermi; 617 delta = 0.1*fermi*fermi;
618 e1 = 0.3*fermi; 618 e1 = 0.3*fermi;
619 e2 = 0.35*fermi; 619 e2 = 0.35*fermi;
620 } 620 }
621 else // as proton, if were not defined 621 else // as proton, if were not defined
622 { 622 {
623 diffuse = 0.63*fermi; 623 diffuse = 0.63*fermi;
624 gamma = 0.3*fermi; 624 gamma = 0.3*fermi;
625 delta = 0.1*fermi*fermi; 625 delta = 0.1*fermi*fermi;
626 e1 = 0.3*fermi; 626 e1 = 0.3*fermi;
627 e2 = 0.35*fermi; 627 e2 = 0.35*fermi;
628 } 628 }
629 G4double lambda = 15.; // 15 ok 629 G4double lambda = 15.; // 15 ok
630 // G4double kgamma = fWaveVector*gamma; 630 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
631 G4double kgamma = lambda*(1.-G4Exp(-fWave 631 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
632 632
633 // G4cout<<"kgamma = "<<kgamma<<G4endl; 633 // G4cout<<"kgamma = "<<kgamma<<G4endl;
634 634
635 if( fAddCoulomb ) // add Coulomb correction 635 if( fAddCoulomb ) // add Coulomb correction
636 { 636 {
637 G4double sinHalfTheta = theta*0.5; // std 637 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
638 G4double sinHalfTheta2 = sinHalfTheta*sinH 638 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
639 639
640 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 640 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
641 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 641 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
642 } 642 }
643 G4double kgamma2 = kgamma*kgamma; 643 G4double kgamma2 = kgamma*kgamma;
644 644
645 // G4double dk2t = delta*fWaveVector*fWaveV 645 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
646 // G4cout<<"dk2t = "<<dk2t<<G4endl; 646 // G4cout<<"dk2t = "<<dk2t<<G4endl;
647 // G4double dk2t2 = dk2t*dk2t; 647 // G4double dk2t2 = dk2t*dk2t;
648 // G4double pikdt = pi*fWaveVector*diffuse*t 648 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
649 649
650 G4double pikdt = lambda*(1. - G4Exp( -pi* 650 G4double pikdt = lambda*(1. - G4Exp( -pi*fWaveVector*diffuse*theta/lambda ) ); // wavek*delta;
651 651
652 // G4cout<<"pikdt = "<<pikdt<<G4endl; 652 // G4cout<<"pikdt = "<<pikdt<<G4endl;
653 653
654 damp = DampFactor( pikdt ); 654 damp = DampFactor( pikdt );
655 damp2 = damp*damp; 655 damp2 = damp*damp;
656 656
657 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVe 657 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVector*fWaveVector;
658 G4double e2dk3t = -2.*e2*delta*fWaveVector* 658 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
659 659
660 sigma = kgamma2; 660 sigma = kgamma2;
661 // sigma += dk2t2; 661 // sigma += dk2t2;
662 sigma *= bzero2; 662 sigma *= bzero2;
663 sigma += mode2k2*bone2; 663 sigma += mode2k2*bone2;
664 sigma += e2dk3t*bzero*bone; 664 sigma += e2dk3t*bzero*bone;
665 665
666 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 666 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
667 sigma += kr2*bonebyarg2; // correction at J 667 sigma += kr2*bonebyarg2; // correction at J1()/()
668 668
669 sigma *= damp2; // *rad*rad; 669 sigma *= damp2; // *rad*rad;
670 670
671 return sigma; 671 return sigma;
672 } 672 }
673 673
674 674
675 ////////////////////////////////////////////// 675 ////////////////////////////////////////////////////////////////////////////
676 // 676 //
677 // return differential elastic probability 2*p 677 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
678 678
679 G4double 679 G4double
680 G4DiffuseElastic::GetIntegrandFunction( G4doub 680 G4DiffuseElastic::GetIntegrandFunction( G4double alpha )
681 { 681 {
682 G4double result; 682 G4double result;
683 683
684 result = GetDiffElasticSumProbA(alpha); 684 result = GetDiffElasticSumProbA(alpha);
685 685
686 // result *= 2*pi*std::sin(theta); 686 // result *= 2*pi*std::sin(theta);
687 687
688 return result; 688 return result;
689 } 689 }
690 690
691 ////////////////////////////////////////////// 691 ////////////////////////////////////////////////////////////////////////////
692 // 692 //
693 // return integral elastic cross section d(sig 693 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
694 694
695 G4double 695 G4double
696 G4DiffuseElastic::IntegralElasticProb( const 696 G4DiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle,
697 G4doub 697 G4double theta,
698 G4double momentum, 698 G4double momentum,
699 G4doub 699 G4double A )
700 { 700 {
701 G4double result; 701 G4double result;
702 fParticle = particle; 702 fParticle = particle;
703 fWaveVector = momentum/hbarc; 703 fWaveVector = momentum/hbarc;
704 fAtomicWeight = A; 704 fAtomicWeight = A;
705 705
706 fNuclearRadius = CalculateNuclearRad(A); 706 fNuclearRadius = CalculateNuclearRad(A);
707 707
708 708
709 G4Integrator<G4DiffuseElastic,G4double(G4Dif 709 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
710 710
711 // result = integral.Legendre10(this,&G4Diff 711 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
712 result = integral.Legendre96(this,&G4Diffuse 712 result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
713 713
714 return result; 714 return result;
715 } 715 }
716 716
717 ////////////////////////////////////////////// 717 ////////////////////////////////////////////////////////////////////////////
718 // 718 //
719 // Return inv momentum transfer -t > 0 719 // Return inv momentum transfer -t > 0
720 720
721 G4double G4DiffuseElastic::SampleT( const G4Pa 721 G4double G4DiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
722 { 722 {
723 G4double theta = SampleThetaCMS( aParticle, 723 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
724 G4double t = 2*p*p*( 1 - std::cos(theta) 724 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
725 return t; 725 return t;
726 } 726 }
727 727
728 ////////////////////////////////////////////// 728 ////////////////////////////////////////////////////////////////////////////
729 // 729 //
730 // Return scattering angle sampled in cms 730 // Return scattering angle sampled in cms
731 731
732 732
733 G4double 733 G4double
734 G4DiffuseElastic::SampleThetaCMS(const G4Parti 734 G4DiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle,
735 G4doubl 735 G4double momentum, G4double A)
736 { 736 {
737 G4int i, iMax = 100; 737 G4int i, iMax = 100;
738 G4double norm, theta1, theta2, thetaMax; << 738 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
739 G4double result = 0., sum = 0.; <<
740 739
741 fParticle = particle; 740 fParticle = particle;
742 fWaveVector = momentum/hbarc; 741 fWaveVector = momentum/hbarc;
743 fAtomicWeight = A; 742 fAtomicWeight = A;
744 743
745 fNuclearRadius = CalculateNuclearRad(A); 744 fNuclearRadius = CalculateNuclearRad(A);
746 745
747 thetaMax = 10.174/fWaveVector/fNuclearRadius 746 thetaMax = 10.174/fWaveVector/fNuclearRadius;
748 747
749 if (thetaMax > pi) thetaMax = pi; 748 if (thetaMax > pi) thetaMax = pi;
750 749
751 G4Integrator<G4DiffuseElastic,G4double(G4Dif 750 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
752 751
753 // result = integral.Legendre10(this,&G4Diff 752 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
754 norm = integral.Legendre96(this,&G4DiffuseEl 753 norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
755 754
756 norm *= G4UniformRand(); 755 norm *= G4UniformRand();
757 756
758 for(i = 1; i <= iMax; i++) 757 for(i = 1; i <= iMax; i++)
759 { 758 {
760 theta1 = (i-1)*thetaMax/iMax; 759 theta1 = (i-1)*thetaMax/iMax;
761 theta2 = i*thetaMax/iMax; 760 theta2 = i*thetaMax/iMax;
762 sum += integral.Legendre10(this,&G4Diffu 761 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
763 762
764 if ( sum >= norm ) 763 if ( sum >= norm )
765 { 764 {
766 result = 0.5*(theta1 + theta2); 765 result = 0.5*(theta1 + theta2);
767 break; 766 break;
768 } 767 }
769 } 768 }
770 if (i > iMax ) result = 0.5*(theta1 + theta2 769 if (i > iMax ) result = 0.5*(theta1 + theta2);
771 770
772 G4double sigma = pi*thetaMax/iMax; 771 G4double sigma = pi*thetaMax/iMax;
773 772
774 result += G4RandGauss::shoot(0.,sigma); 773 result += G4RandGauss::shoot(0.,sigma);
775 774
776 if(result < 0.) result = 0.; 775 if(result < 0.) result = 0.;
777 if(result > thetaMax) result = thetaMax; 776 if(result > thetaMax) result = thetaMax;
778 777
779 return result; 778 return result;
780 } 779 }
781 780
782 ////////////////////////////////////////////// 781 /////////////////////////////////////////////////////////////////////////////
783 ///////////////////// Table preparation and r 782 ///////////////////// Table preparation and reading ////////////////////////
784 ////////////////////////////////////////////// 783 ////////////////////////////////////////////////////////////////////////////
785 // 784 //
786 // Return inv momentum transfer -t > 0 from in 785 // Return inv momentum transfer -t > 0 from initialisation table
787 786
788 G4double G4DiffuseElastic::SampleInvariantT( c 787 G4double G4DiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
789 788 G4int Z, G4int A)
790 { 789 {
791 fParticle = aParticle; 790 fParticle = aParticle;
792 G4double m1 = fParticle->GetPDGMass(), t; 791 G4double m1 = fParticle->GetPDGMass(), t;
793 G4double totElab = std::sqrt(m1*m1+p*p); 792 G4double totElab = std::sqrt(m1*m1+p*p);
794 G4double mass2 = G4NucleiProperties::GetNucl 793 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
795 G4LorentzVector lv1(p,0.0,0.0,totElab); 794 G4LorentzVector lv1(p,0.0,0.0,totElab);
796 G4LorentzVector lv(0.0,0.0,0.0,mass2); 795 G4LorentzVector lv(0.0,0.0,0.0,mass2);
797 lv += lv1; 796 lv += lv1;
798 797
799 G4ThreeVector bst = lv.boostVector(); 798 G4ThreeVector bst = lv.boostVector();
800 lv1.boost(-bst); 799 lv1.boost(-bst);
801 800
802 G4ThreeVector p1 = lv1.vect(); 801 G4ThreeVector p1 = lv1.vect();
803 G4double momentumCMS = p1.mag(); 802 G4double momentumCMS = p1.mag();
804 803
805 if( aParticle == theNeutron) 804 if( aParticle == theNeutron)
806 { 805 {
807 G4double Tmax = NeutronTuniform( Z ); 806 G4double Tmax = NeutronTuniform( Z );
808 G4double pCMS2 = momentumCMS*momentumCMS; 807 G4double pCMS2 = momentumCMS*momentumCMS;
809 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1; 808 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
810 809
811 if( Tkin <= Tmax ) 810 if( Tkin <= Tmax )
812 { 811 {
813 t = 4.*pCMS2*G4UniformRand(); 812 t = 4.*pCMS2*G4UniformRand();
814 // G4cout<<Tkin<<", "<<Tmax<<", "<<std:: 813 // G4cout<<Tkin<<", "<<Tmax<<", "<<std::sqrt(t)<<"; ";
815 return t; 814 return t;
816 } 815 }
817 } 816 }
818 817
819 t = SampleTableT( aParticle, momentumCMS, G 818 t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
820 819
821 return t; 820 return t;
822 } 821 }
823 822
824 ////////////////////////////////////////////// 823 ///////////////////////////////////////////////////////
825 824
826 G4double G4DiffuseElastic::NeutronTuniform(G4i 825 G4double G4DiffuseElastic::NeutronTuniform(G4int Z)
827 { 826 {
828 G4double elZ = G4double(Z); 827 G4double elZ = G4double(Z);
829 elZ -= 1.; 828 elZ -= 1.;
830 // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.; 829 // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.;
831 G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.; 830 G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.;
832 return Tkin; 831 return Tkin;
833 } 832 }
834 833
835 834
836 ////////////////////////////////////////////// 835 ////////////////////////////////////////////////////////////////////////////
837 // 836 //
838 // Return inv momentum transfer -t > 0 from in 837 // Return inv momentum transfer -t > 0 from initialisation table
839 838
840 G4double G4DiffuseElastic::SampleTableT( const 839 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
841 840 G4double Z, G4double A)
842 { 841 {
843 G4double alpha = SampleTableThetaCMS( aParti 842 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
844 G4double t = 2*p*p*( 1 - std::cos(std::s 843 G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
845 // G4double t = p*p*alpha; / 844 // G4double t = p*p*alpha; // -t !!!
846 return t; 845 return t;
847 } 846 }
848 847
849 ////////////////////////////////////////////// 848 ////////////////////////////////////////////////////////////////////////////
850 // 849 //
851 // Return scattering angle2 sampled in cms acc 850 // Return scattering angle2 sampled in cms according to precalculated table.
852 851
853 852
854 G4double 853 G4double
855 G4DiffuseElastic::SampleTableThetaCMS(const G4 854 G4DiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
856 G4doubl 855 G4double momentum, G4double Z, G4double A)
857 { 856 {
858 std::size_t iElement; << 857 size_t iElement;
859 G4int iMomentum, iAngle; 858 G4int iMomentum, iAngle;
860 G4double randAngle, position, theta1, theta2 859 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
861 G4double m1 = particle->GetPDGMass(); 860 G4double m1 = particle->GetPDGMass();
862 861
863 for(iElement = 0; iElement < fElementNumberV 862 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
864 { 863 {
865 if( std::fabs(Z - fElementNumberVector[iEl 864 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
866 } 865 }
867 if ( iElement == fElementNumberVector.size() 866 if ( iElement == fElementNumberVector.size() )
868 { 867 {
869 InitialiseOnFly(Z,A); // table preparation 868 InitialiseOnFly(Z,A); // table preparation, if needed
870 869
871 // iElement--; 870 // iElement--;
872 871
873 // G4cout << "G4DiffuseElastic: Element wi 872 // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
874 // << " is not found, return zero angle" < 873 // << " is not found, return zero angle" << G4endl;
875 // return 0.; // no table for this element 874 // return 0.; // no table for this element
876 } 875 }
877 // G4cout<<"iElement = "<<iElement<<G4endl; 876 // G4cout<<"iElement = "<<iElement<<G4endl;
878 877
879 fAngleTable = fAngleBank[iElement]; 878 fAngleTable = fAngleBank[iElement];
880 879
881 G4double kinE = std::sqrt(momentum*momentum 880 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882 881
883 for( iMomentum = 0; iMomentum < fEnergyBin; 882 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
884 { 883 {
885 if( kinE < fEnergyVector->GetLowEdgeEnergy 884 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
886 } 885 }
887 if ( iMomentum >= fEnergyBin ) iMomentum = f 886 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
888 if ( iMomentum < 0 ) iMomentum = 0 887 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
889 888
890 // G4cout<<"iMomentum = "<<iMomentum<<G4endl 889 // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
891 890
892 if (iMomentum == fEnergyBin -1 || iMomentum 891 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
893 { 892 {
894 position = (*(*fAngleTable)(iMomentum))(fA 893 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
895 894
896 // G4cout<<"position = "<<position<<G4endl 895 // G4cout<<"position = "<<position<<G4endl;
897 896
898 for(iAngle = 0; iAngle < fAngleBin-1; iAng 897 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
899 { 898 {
900 if( position > (*(*fAngleTable)(iMomentu 899 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
901 } 900 }
902 if (iAngle >= fAngleBin-1) iAngle = fAngle 901 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
903 902
904 // G4cout<<"iAngle = "<<iAngle<<G4endl; 903 // G4cout<<"iAngle = "<<iAngle<<G4endl;
905 904
906 randAngle = GetScatteringAngle(iMomentum, 905 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
907 906
908 // G4cout<<"randAngle = "<<randAngle<<G4en 907 // G4cout<<"randAngle = "<<randAngle<<G4endl;
909 } 908 }
910 else // kinE inside between energy table ed 909 else // kinE inside between energy table edges
911 { 910 {
912 // position = (*(*fAngleTable)(iMomentum)) 911 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
913 position = (*(*fAngleTable)(iMomentum))(0) 912 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
914 913
915 // G4cout<<"position = "<<position<<G4endl 914 // G4cout<<"position = "<<position<<G4endl;
916 915
917 for(iAngle = 0; iAngle < fAngleBin-1; iAng 916 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
918 { 917 {
919 // if( position < (*(*fAngleTable)(iMome 918 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
920 if( position > (*(*fAngleTable)(iMomentu 919 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
921 } 920 }
922 if (iAngle >= fAngleBin-1) iAngle = fAngle 921 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
923 922
924 // G4cout<<"iAngle = "<<iAngle<<G4endl; 923 // G4cout<<"iAngle = "<<iAngle<<G4endl;
925 924
926 theta2 = GetScatteringAngle(iMomentum, iA 925 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
927 926
928 // G4cout<<"theta2 = "<<theta2<<G4endl; 927 // G4cout<<"theta2 = "<<theta2<<G4endl;
929 E2 = fEnergyVector->GetLowEdgeEnergy(iMome 928 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
930 929
931 // G4cout<<"E2 = "<<E2<<G4endl; 930 // G4cout<<"E2 = "<<E2<<G4endl;
932 931
933 iMomentum--; 932 iMomentum--;
934 933
935 // position = (*(*fAngleTable)(iMomentum)) 934 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
936 935
937 // G4cout<<"position = "<<position<<G4endl 936 // G4cout<<"position = "<<position<<G4endl;
938 937
939 for(iAngle = 0; iAngle < fAngleBin-1; iAng 938 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
940 { 939 {
941 // if( position < (*(*fAngleTable)(iMome 940 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
942 if( position > (*(*fAngleTable)(iMomentu 941 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
943 } 942 }
944 if (iAngle >= fAngleBin-1) iAngle = fAngle 943 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
945 944
946 theta1 = GetScatteringAngle(iMomentum, iA 945 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
947 946
948 // G4cout<<"theta1 = "<<theta1<<G4endl; 947 // G4cout<<"theta1 = "<<theta1<<G4endl;
949 948
950 E1 = fEnergyVector->GetLowEdgeEnergy(iMome 949 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
951 950
952 // G4cout<<"E1 = "<<E1<<G4endl; 951 // G4cout<<"E1 = "<<E1<<G4endl;
953 952
954 W = 1.0/(E2 - E1); 953 W = 1.0/(E2 - E1);
955 W1 = (E2 - kinE)*W; 954 W1 = (E2 - kinE)*W;
956 W2 = (kinE - E1)*W; 955 W2 = (kinE - E1)*W;
957 956
958 randAngle = W1*theta1 + W2*theta2; 957 randAngle = W1*theta1 + W2*theta2;
959 958
960 // randAngle = theta2; 959 // randAngle = theta2;
961 // G4cout<<"randAngle = "<<randAngle<<G4en 960 // G4cout<<"randAngle = "<<randAngle<<G4endl;
962 } 961 }
963 // G4double angle = randAngle; 962 // G4double angle = randAngle;
964 // if (randAngle > 0.) randAngle /= 2*pi*std 963 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
965 964
966 if(randAngle < 0.) randAngle = 0.; 965 if(randAngle < 0.) randAngle = 0.;
967 966
968 return randAngle; 967 return randAngle;
969 } 968 }
970 969
971 ////////////////////////////////////////////// 970 //////////////////////////////////////////////////////////////////////////////
972 // 971 //
973 // Initialisation for given particle on fly us 972 // Initialisation for given particle on fly using new element number
974 973
975 void G4DiffuseElastic::InitialiseOnFly(G4doubl 974 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
976 { 975 {
977 fAtomicNumber = Z; // atomic number 976 fAtomicNumber = Z; // atomic number
978 fAtomicWeight = G4NistManager::Instance()-> 977 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
979 978
980 fNuclearRadius = CalculateNuclearRad(fAtomic 979 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
981 980
982 if( verboseLevel > 0 ) 981 if( verboseLevel > 0 )
983 { 982 {
984 G4cout<<"G4DiffuseElastic::InitialiseOnFly 983 G4cout<<"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
985 <<Z<<"; and A = "<<A<<G4endl; 984 <<Z<<"; and A = "<<A<<G4endl;
986 } 985 }
987 fElementNumberVector.push_back(fAtomicNumber 986 fElementNumberVector.push_back(fAtomicNumber);
988 987
989 BuildAngleTable(); 988 BuildAngleTable();
990 989
991 fAngleBank.push_back(fAngleTable); 990 fAngleBank.push_back(fAngleTable);
992 991
993 return; 992 return;
994 } 993 }
995 994
996 ////////////////////////////////////////////// 995 ///////////////////////////////////////////////////////////////////////////////
997 // 996 //
998 // Build for given particle and element table 997 // Build for given particle and element table of momentum, angle probability.
999 // For the moment in lab system. 998 // For the moment in lab system.
1000 999
1001 void G4DiffuseElastic::BuildAngleTable() 1000 void G4DiffuseElastic::BuildAngleTable()
1002 { 1001 {
1003 G4int i, j; 1002 G4int i, j;
1004 G4double partMom, kinE, a = 0., z = fPartic 1003 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1005 G4double alpha1, alpha2, alphaMax, alphaCou 1004 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1006 1005
1007 G4Integrator<G4DiffuseElastic,G4double(G4Di 1006 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
1008 1007
1009 fAngleTable = new G4PhysicsTable( fEnergyBi 1008 fAngleTable = new G4PhysicsTable( fEnergyBin );
1010 1009
1011 for( i = 0; i < fEnergyBin; i++) 1010 for( i = 0; i < fEnergyBin; i++)
1012 { 1011 {
1013 kinE = fEnergyVector->GetLowEdgeEn 1012 kinE = fEnergyVector->GetLowEdgeEnergy(i);
1014 partMom = std::sqrt( kinE*(kinE + 2*m 1013 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1015 1014
1016 fWaveVector = partMom/hbarc; 1015 fWaveVector = partMom/hbarc;
1017 1016
1018 G4double kR = fWaveVector*fNuclearRad 1017 G4double kR = fWaveVector*fNuclearRadius;
1019 G4double kR2 = kR*kR; 1018 G4double kR2 = kR*kR;
1020 G4double kRmax = 18.6; // 10.6; 10.6, 18 1019 G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1021 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; 1020 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
1022 // G4double kRlim = 1.2; 1021 // G4double kRlim = 1.2;
1023 // G4double kRlim2 = kRlim*kRlim/kR2; 1022 // G4double kRlim2 = kRlim*kRlim/kR2;
1024 1023
1025 alphaMax = kRmax*kRmax/kR2; 1024 alphaMax = kRmax*kRmax/kR2;
1026 1025
1027 1026
1028 // if (alphaMax > 4.) alphaMax = 4.; // 1027 // if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1029 // if ( alphaMax > 4. || alphaMax < 1. ) 1028 // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = 15.; // vmg27.11.14
1030 1029
1031 // if ( alphaMax > 4. || alphaMax < 1. ) 1030 // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg06.01.15
1032 1031
1033 // G4cout<<"alphaMax = "<<alphaMax<<", "; 1032 // G4cout<<"alphaMax = "<<alphaMax<<", ";
1034 1033
1035 if ( alphaMax >= CLHEP::pi*CLHEP::pi ) al 1034 if ( alphaMax >= CLHEP::pi*CLHEP::pi ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg21.10.15
1036 1035
1037 alphaCoulomb = kRcoul*kRcoul/kR2; 1036 alphaCoulomb = kRcoul*kRcoul/kR2;
1038 1037
1039 if( z ) 1038 if( z )
1040 { 1039 {
1041 a = partMom/m1; // be 1040 a = partMom/m1; // beta*gamma for m1
1042 fBeta = a/std::sqrt(1+a*a); 1041 fBeta = a/std::sqrt(1+a*a);
1043 fZommerfeld = CalculateZommerfeld( fBet 1042 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1044 fAm = CalculateAm( partMom, fZo 1043 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1045 } 1044 }
1046 G4PhysicsFreeVector* angleVector = new G4 1045 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1047 1046
1048 // G4PhysicsLogVector* angleBins = new G 1047 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1049 1048
1050 G4double delth = alphaMax/fAngleBin; 1049 G4double delth = alphaMax/fAngleBin;
1051 1050
1052 sum = 0.; 1051 sum = 0.;
1053 1052
1054 // fAddCoulomb = false; 1053 // fAddCoulomb = false;
1055 fAddCoulomb = true; 1054 fAddCoulomb = true;
1056 1055
1057 // for(j = 1; j < fAngleBin; j++) 1056 // for(j = 1; j < fAngleBin; j++)
1058 for(j = fAngleBin-1; j >= 1; j--) 1057 for(j = fAngleBin-1; j >= 1; j--)
1059 { 1058 {
1060 // alpha1 = angleBins->GetLowEdgeEnergy 1059 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1061 // alpha2 = angleBins->GetLowEdgeEnergy 1060 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1062 1061
1063 // alpha1 = alphaMax*(j-1)/fAngleBin; 1062 // alpha1 = alphaMax*(j-1)/fAngleBin;
1064 // alpha2 = alphaMax*( j )/fAngleBin; 1063 // alpha2 = alphaMax*( j )/fAngleBin;
1065 1064
1066 alpha1 = delth*(j-1); 1065 alpha1 = delth*(j-1);
1067 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1066 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1068 alpha2 = alpha1 + delth; 1067 alpha2 = alpha1 + delth;
1069 1068
1070 // if( ( alpha2 > alphaCoulomb ) && z ) 1069 // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1071 if( ( alpha1 < alphaCoulomb ) && z ) fA 1070 if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
1072 1071
1073 delta = integral.Legendre10(this, &G4Di 1072 delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1074 // delta = integral.Legendre96(this, &G 1073 // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1075 1074
1076 sum += delta; 1075 sum += delta;
1077 1076
1078 angleVector->PutValue( j-1 , alpha1, su 1077 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1079 // G4cout<<"j-1 = "<<j-1<<"; alpha 1078 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2 << " delta "<< delta <<"; sum = "<<sum<<G4endl;
1080 } 1079 }
1081 fAngleTable->insertAt(i, angleVector); 1080 fAngleTable->insertAt(i, angleVector);
1082 1081
1083 // delete[] angleVector; 1082 // delete[] angleVector;
1084 // delete[] angleBins; 1083 // delete[] angleBins;
1085 } 1084 }
1086 return; 1085 return;
1087 } 1086 }
1088 1087
1089 ///////////////////////////////////////////// 1088 /////////////////////////////////////////////////////////////////////////////////
1090 // 1089 //
1091 // 1090 //
1092 1091
1093 G4double 1092 G4double
1094 G4DiffuseElastic:: GetScatteringAngle( G4int 1093 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1095 { 1094 {
1096 G4double x1, x2, y1, y2, randAngle; 1095 G4double x1, x2, y1, y2, randAngle;
1097 1096
1098 if( iAngle == 0 ) 1097 if( iAngle == 0 )
1099 { 1098 {
1100 randAngle = (*fAngleTable)(iMomentum)->Ge 1099 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1101 // iAngle++; 1100 // iAngle++;
1102 } 1101 }
1103 else 1102 else
1104 { 1103 {
1105 if ( iAngle >= G4int((*fAngleTable)(iMome 1104 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1106 { 1105 {
1107 iAngle = G4int((*fAngleTable)(iMomentum << 1106 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1108 } 1107 }
1109 y1 = (*(*fAngleTable)(iMomentum))(iAngle- 1108 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1110 y2 = (*(*fAngleTable)(iMomentum))(iAngle) 1109 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1111 1110
1112 x1 = (*fAngleTable)(iMomentum)->GetLowEdg 1111 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1113 x2 = (*fAngleTable)(iMomentum)->GetLowEdg 1112 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1114 1113
1115 if ( x1 == x2 ) randAngle = x2; 1114 if ( x1 == x2 ) randAngle = x2;
1116 else 1115 else
1117 { 1116 {
1118 if ( y1 == y2 ) randAngle = x1 + ( x2 - 1117 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1119 else 1118 else
1120 { 1119 {
1121 randAngle = x1 + ( position - y1 )*( 1120 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1122 } 1121 }
1123 } 1122 }
1124 } 1123 }
1125 return randAngle; 1124 return randAngle;
1126 } 1125 }
1127 1126
1128 1127
1129 1128
1130 ///////////////////////////////////////////// 1129 ////////////////////////////////////////////////////////////////////////////
1131 // 1130 //
1132 // Return scattering angle sampled in lab sys 1131 // Return scattering angle sampled in lab system (target at rest)
1133 1132
1134 1133
1135 1134
1136 G4double 1135 G4double
1137 G4DiffuseElastic::SampleThetaLab( const G4Had 1136 G4DiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
1138 G4dou 1137 G4double tmass, G4double A)
1139 { 1138 {
1140 const G4ParticleDefinition* theParticle = a 1139 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1141 G4double m1 = theParticle->GetPDGMass(); 1140 G4double m1 = theParticle->GetPDGMass();
1142 G4double plab = aParticle->GetTotalMomentum 1141 G4double plab = aParticle->GetTotalMomentum();
1143 G4LorentzVector lv1 = aParticle->Get4Moment 1142 G4LorentzVector lv1 = aParticle->Get4Momentum();
1144 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1143 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1145 lv += lv1; 1144 lv += lv1;
1146 1145
1147 G4ThreeVector bst = lv.boostVector(); 1146 G4ThreeVector bst = lv.boostVector();
1148 lv1.boost(-bst); 1147 lv1.boost(-bst);
1149 1148
1150 G4ThreeVector p1 = lv1.vect(); 1149 G4ThreeVector p1 = lv1.vect();
1151 G4double ptot = p1.mag(); 1150 G4double ptot = p1.mag();
1152 G4double tmax = 4.0*ptot*ptot; 1151 G4double tmax = 4.0*ptot*ptot;
1153 G4double t = 0.0; 1152 G4double t = 0.0;
1154 1153
1155 1154
1156 // 1155 //
1157 // Sample t 1156 // Sample t
1158 // 1157 //
1159 1158
1160 t = SampleT( theParticle, ptot, A); 1159 t = SampleT( theParticle, ptot, A);
1161 1160
1162 // NaN finder 1161 // NaN finder
1163 if(!(t < 0.0 || t >= 0.0)) 1162 if(!(t < 0.0 || t >= 0.0))
1164 { 1163 {
1165 if (verboseLevel > 0) 1164 if (verboseLevel > 0)
1166 { 1165 {
1167 G4cout << "G4DiffuseElastic:WARNING: A 1166 G4cout << "G4DiffuseElastic:WARNING: A = " << A
1168 << " mom(GeV)= " << plab/GeV 1167 << " mom(GeV)= " << plab/GeV
1169 << " S-wave will be sampled" 1168 << " S-wave will be sampled"
1170 << G4endl; 1169 << G4endl;
1171 } 1170 }
1172 t = G4UniformRand()*tmax; 1171 t = G4UniformRand()*tmax;
1173 } 1172 }
1174 if(verboseLevel>1) 1173 if(verboseLevel>1)
1175 { 1174 {
1176 G4cout <<" t= " << t << " tmax= " << tmax 1175 G4cout <<" t= " << t << " tmax= " << tmax
1177 << " ptot= " << ptot << G4endl; 1176 << " ptot= " << ptot << G4endl;
1178 } 1177 }
1179 // Sampling of angles in CM system 1178 // Sampling of angles in CM system
1180 1179
1181 G4double phi = G4UniformRand()*twopi; 1180 G4double phi = G4UniformRand()*twopi;
1182 G4double cost = 1. - 2.0*t/tmax; 1181 G4double cost = 1. - 2.0*t/tmax;
1183 G4double sint; 1182 G4double sint;
1184 1183
1185 if( cost >= 1.0 ) 1184 if( cost >= 1.0 )
1186 { 1185 {
1187 cost = 1.0; 1186 cost = 1.0;
1188 sint = 0.0; 1187 sint = 0.0;
1189 } 1188 }
1190 else if( cost <= -1.0) 1189 else if( cost <= -1.0)
1191 { 1190 {
1192 cost = -1.0; 1191 cost = -1.0;
1193 sint = 0.0; 1192 sint = 0.0;
1194 } 1193 }
1195 else 1194 else
1196 { 1195 {
1197 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1196 sint = std::sqrt((1.0-cost)*(1.0+cost));
1198 } 1197 }
1199 if (verboseLevel>1) 1198 if (verboseLevel>1)
1200 { 1199 {
1201 G4cout << "cos(t)=" << cost << " std::sin 1200 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1202 } 1201 }
1203 G4ThreeVector v1(sint*std::cos(phi),sint*st 1202 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1204 v1 *= ptot; 1203 v1 *= ptot;
1205 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1204 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1206 1205
1207 nlv1.boost(bst); 1206 nlv1.boost(bst);
1208 1207
1209 G4ThreeVector np1 = nlv1.vect(); 1208 G4ThreeVector np1 = nlv1.vect();
1210 1209
1211 // G4double theta = std::acos( np1.z()/np 1210 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1212 1211
1213 G4double theta = np1.theta(); 1212 G4double theta = np1.theta();
1214 1213
1215 return theta; 1214 return theta;
1216 } 1215 }
1217 1216
1218 ///////////////////////////////////////////// 1217 ////////////////////////////////////////////////////////////////////////////
1219 // 1218 //
1220 // Return scattering angle in lab system (tar 1219 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1221 1220
1222 1221
1223 1222
1224 G4double 1223 G4double
1225 G4DiffuseElastic::ThetaCMStoThetaLab( const G 1224 G4DiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
1226 G4dou 1225 G4double tmass, G4double thetaCMS)
1227 { 1226 {
1228 const G4ParticleDefinition* theParticle = a 1227 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1229 G4double m1 = theParticle->GetPDGMass(); 1228 G4double m1 = theParticle->GetPDGMass();
1230 // G4double plab = aParticle->GetTotalMomen 1229 // G4double plab = aParticle->GetTotalMomentum();
1231 G4LorentzVector lv1 = aParticle->Get4Moment 1230 G4LorentzVector lv1 = aParticle->Get4Momentum();
1232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1231 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1233 1232
1234 lv += lv1; 1233 lv += lv1;
1235 1234
1236 G4ThreeVector bst = lv.boostVector(); 1235 G4ThreeVector bst = lv.boostVector();
1237 1236
1238 lv1.boost(-bst); 1237 lv1.boost(-bst);
1239 1238
1240 G4ThreeVector p1 = lv1.vect(); 1239 G4ThreeVector p1 = lv1.vect();
1241 G4double ptot = p1.mag(); 1240 G4double ptot = p1.mag();
1242 1241
1243 G4double phi = G4UniformRand()*twopi; 1242 G4double phi = G4UniformRand()*twopi;
1244 G4double cost = std::cos(thetaCMS); 1243 G4double cost = std::cos(thetaCMS);
1245 G4double sint; 1244 G4double sint;
1246 1245
1247 if( cost >= 1.0 ) 1246 if( cost >= 1.0 )
1248 { 1247 {
1249 cost = 1.0; 1248 cost = 1.0;
1250 sint = 0.0; 1249 sint = 0.0;
1251 } 1250 }
1252 else if( cost <= -1.0) 1251 else if( cost <= -1.0)
1253 { 1252 {
1254 cost = -1.0; 1253 cost = -1.0;
1255 sint = 0.0; 1254 sint = 0.0;
1256 } 1255 }
1257 else 1256 else
1258 { 1257 {
1259 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1258 sint = std::sqrt((1.0-cost)*(1.0+cost));
1260 } 1259 }
1261 if (verboseLevel>1) 1260 if (verboseLevel>1)
1262 { 1261 {
1263 G4cout << "cos(tcms)=" << cost << " std:: 1262 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1264 } 1263 }
1265 G4ThreeVector v1(sint*std::cos(phi),sint*st 1264 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1266 v1 *= ptot; 1265 v1 *= ptot;
1267 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1266 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1268 1267
1269 nlv1.boost(bst); 1268 nlv1.boost(bst);
1270 1269
1271 G4ThreeVector np1 = nlv1.vect(); 1270 G4ThreeVector np1 = nlv1.vect();
1272 1271
1273 1272
1274 G4double thetaLab = np1.theta(); 1273 G4double thetaLab = np1.theta();
1275 1274
1276 return thetaLab; 1275 return thetaLab;
1277 } 1276 }
1278 ///////////////////////////////////////////// 1277 ////////////////////////////////////////////////////////////////////////////
1279 // 1278 //
1280 // Return scattering angle in CMS system (tar 1279 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1281 1280
1282 1281
1283 1282
1284 G4double 1283 G4double
1285 G4DiffuseElastic::ThetaLabToThetaCMS( const G 1284 G4DiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
1286 G4dou 1285 G4double tmass, G4double thetaLab)
1287 { 1286 {
1288 const G4ParticleDefinition* theParticle = a 1287 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1289 G4double m1 = theParticle->GetPDGMass(); 1288 G4double m1 = theParticle->GetPDGMass();
1290 G4double plab = aParticle->GetTotalMomentum 1289 G4double plab = aParticle->GetTotalMomentum();
1291 G4LorentzVector lv1 = aParticle->Get4Moment 1290 G4LorentzVector lv1 = aParticle->Get4Momentum();
1292 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1291 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1293 1292
1294 lv += lv1; 1293 lv += lv1;
1295 1294
1296 G4ThreeVector bst = lv.boostVector(); 1295 G4ThreeVector bst = lv.boostVector();
1297 1296
1298 // lv1.boost(-bst); 1297 // lv1.boost(-bst);
1299 1298
1300 // G4ThreeVector p1 = lv1.vect(); 1299 // G4ThreeVector p1 = lv1.vect();
1301 // G4double ptot = p1.mag(); 1300 // G4double ptot = p1.mag();
1302 1301
1303 G4double phi = G4UniformRand()*twopi; 1302 G4double phi = G4UniformRand()*twopi;
1304 G4double cost = std::cos(thetaLab); 1303 G4double cost = std::cos(thetaLab);
1305 G4double sint; 1304 G4double sint;
1306 1305
1307 if( cost >= 1.0 ) 1306 if( cost >= 1.0 )
1308 { 1307 {
1309 cost = 1.0; 1308 cost = 1.0;
1310 sint = 0.0; 1309 sint = 0.0;
1311 } 1310 }
1312 else if( cost <= -1.0) 1311 else if( cost <= -1.0)
1313 { 1312 {
1314 cost = -1.0; 1313 cost = -1.0;
1315 sint = 0.0; 1314 sint = 0.0;
1316 } 1315 }
1317 else 1316 else
1318 { 1317 {
1319 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1318 sint = std::sqrt((1.0-cost)*(1.0+cost));
1320 } 1319 }
1321 if (verboseLevel>1) 1320 if (verboseLevel>1)
1322 { 1321 {
1323 G4cout << "cos(tlab)=" << cost << " std:: 1322 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1324 } 1323 }
1325 G4ThreeVector v1(sint*std::cos(phi),sint*st 1324 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1326 v1 *= plab; 1325 v1 *= plab;
1327 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1326 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1328 1327
1329 nlv1.boost(-bst); 1328 nlv1.boost(-bst);
1330 1329
1331 G4ThreeVector np1 = nlv1.vect(); 1330 G4ThreeVector np1 = nlv1.vect();
1332 1331
1333 1332
1334 G4double thetaCMS = np1.theta(); 1333 G4double thetaCMS = np1.theta();
1335 1334
1336 return thetaCMS; 1335 return thetaCMS;
1337 } 1336 }
1338 1337
1339 ///////////////////////////////////////////// 1338 ///////////////////////////////////////////////////////////////////////////////
1340 // 1339 //
1341 // Test for given particle and element table 1340 // Test for given particle and element table of momentum, angle probability.
1342 // For the moment in lab system. 1341 // For the moment in lab system.
1343 1342
1344 void G4DiffuseElastic::TestAngleTable(const G 1343 void G4DiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1345 G 1344 G4double Z, G4double A)
1346 { 1345 {
1347 fAtomicNumber = Z; // atomic number 1346 fAtomicNumber = Z; // atomic number
1348 fAtomicWeight = A; // number of nucleo 1347 fAtomicWeight = A; // number of nucleons
1349 fNuclearRadius = CalculateNuclearRad(fAtomi 1348 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1350 1349
1351 1350
1352 1351
1353 G4cout<<"G4DiffuseElastic::TestAngleTable() 1352 G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1354 <<Z<<"; and A = "<<A<<G4endl; 1353 <<Z<<"; and A = "<<A<<G4endl;
1355 1354
1356 fElementNumberVector.push_back(fAtomicNumbe 1355 fElementNumberVector.push_back(fAtomicNumber);
1357 1356
1358 1357
1359 1358
1360 1359
1361 G4int i=0, j; 1360 G4int i=0, j;
1362 G4double a = 0., z = theParticle->GetPDGCha 1361 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1363 G4double alpha1=0., alpha2=0., alphaMax=0., 1362 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1364 G4double deltaL10 = 0., deltaL96 = 0., delt 1363 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1365 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0. 1364 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1366 G4double epsilon = 0.001; 1365 G4double epsilon = 0.001;
1367 1366
1368 G4Integrator<G4DiffuseElastic,G4double(G4Di 1367 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
1369 1368
1370 fAngleTable = new G4PhysicsTable(fEnergyBin 1369 fAngleTable = new G4PhysicsTable(fEnergyBin);
1371 1370
1372 fWaveVector = partMom/hbarc; 1371 fWaveVector = partMom/hbarc;
1373 1372
1374 G4double kR = fWaveVector*fNuclearRadiu 1373 G4double kR = fWaveVector*fNuclearRadius;
1375 G4double kR2 = kR*kR; 1374 G4double kR2 = kR*kR;
1376 G4double kRmax = 10.6; // 10.6, 18, 10.174 1375 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1377 G4double kRcoul = 1.2; // 1.4, 2.5; // on t 1376 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1378 1377
1379 alphaMax = kRmax*kRmax/kR2; 1378 alphaMax = kRmax*kRmax/kR2;
1380 1379
1381 if (alphaMax > 4.) alphaMax = 4.; // vmg05 1380 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1382 1381
1383 alphaCoulomb = kRcoul*kRcoul/kR2; 1382 alphaCoulomb = kRcoul*kRcoul/kR2;
1384 1383
1385 if( z ) 1384 if( z )
1386 { 1385 {
1387 a = partMom/m1; // beta*gamma 1386 a = partMom/m1; // beta*gamma for m1
1388 fBeta = a/std::sqrt(1+a*a); 1387 fBeta = a/std::sqrt(1+a*a);
1389 fZommerfeld = CalculateZommerfeld( fBet 1388 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1390 fAm = CalculateAm( partMom, fZo 1389 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1391 } 1390 }
1392 G4PhysicsFreeVector* angleVector = new G4Ph 1391 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1393 1392
1394 // G4PhysicsLogVector* angleBins = new G4P 1393 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1395 1394
1396 1395
1397 fAddCoulomb = false; 1396 fAddCoulomb = false;
1398 1397
1399 for(j = 1; j < fAngleBin; j++) 1398 for(j = 1; j < fAngleBin; j++)
1400 { 1399 {
1401 // alpha1 = angleBins->GetLowEdgeEnergy 1400 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1402 // alpha2 = angleBins->GetLowEdgeEnergy 1401 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1403 1402
1404 alpha1 = alphaMax*(j-1)/fAngleBin; 1403 alpha1 = alphaMax*(j-1)/fAngleBin;
1405 alpha2 = alphaMax*( j )/fAngleBin; 1404 alpha2 = alphaMax*( j )/fAngleBin;
1406 1405
1407 if( ( alpha2 > alphaCoulomb ) && z ) fAdd 1406 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1408 1407
1409 deltaL10 = integral.Legendre10(this, &G4D 1408 deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1410 deltaL96 = integral.Legendre96(this, &G4D 1409 deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1411 deltaAG = integral.AdaptiveGauss(this, & 1410 deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1412 alpha1 1411 alpha1, alpha2,epsilon);
1413 1412
1414 // G4cout<<alpha1<<"\t"<<std::sqrt(alph 1413 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1415 // <<deltaL10<<"\t"<<deltaL96<<"\t" 1414 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1416 1415
1417 sumL10 += deltaL10; 1416 sumL10 += deltaL10;
1418 sumL96 += deltaL96; 1417 sumL96 += deltaL96;
1419 sumAG += deltaAG; 1418 sumAG += deltaAG;
1420 1419
1421 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d 1420 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1422 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1421 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1423 1422
1424 angleVector->PutValue( j-1 , alpha1, sumL 1423 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1425 } 1424 }
1426 fAngleTable->insertAt(i,angleVector); 1425 fAngleTable->insertAt(i,angleVector);
1427 fAngleBank.push_back(fAngleTable); 1426 fAngleBank.push_back(fAngleTable);
1428 1427
1429 /* 1428 /*
1430 // Integral over all angle range - Bad accu 1429 // Integral over all angle range - Bad accuracy !!!
1431 1430
1432 sumL10 = integral.Legendre10(this, &G4Diffu 1431 sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1433 sumL96 = integral.Legendre96(this, &G4Diffu 1432 sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1434 sumAG = integral.AdaptiveGauss(this, &G4Di 1433 sumAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1435 0., al 1434 0., alpha2,epsilon);
1436 G4cout<<G4endl; 1435 G4cout<<G4endl;
1437 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg 1436 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1438 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1437 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1439 */ 1438 */
1440 return; 1439 return;
1441 } 1440 }
1442 1441
1443 // 1442 //
1444 // 1443 //
1445 ///////////////////////////////////////////// 1444 /////////////////////////////////////////////////////////////////////////////////
1446 1445