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