source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4NuclNuclDiffuseElastic.cc@ 1306

Last change on this file since 1306 was 1228, checked in by garnier, 16 years ago

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