Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4AntiNuclElastic Class Reference

#include <G4AntiNuclElastic.hh>

Inheritance diagram for G4AntiNuclElastic:
G4HadronElastic G4HadronicInteraction

Public Member Functions

 G4AntiNuclElastic ()
 
virtual ~G4AntiNuclElastic ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double SampleThetaCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double SampleThetaLab (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
 
G4double CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
 
G4double CalculateAm (G4double momentum, G4double n, G4double Z)
 
G4double DampFactor (G4double z)
 
G4double BesselJzero (G4double z)
 
G4double BesselJone (G4double z)
 
G4double BesselOneByArg (G4double z)
 
G4double GetcosTeta1 (G4double plab, G4int A)
 
G4ComponentAntiNuclNuclearXS * GetComponentCrossSection ()
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void Description () const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteraction * GetMyPointer () const
 
virtual G4int GetVerboseLevel () const
 
virtual void SetVerboseLevel (G4int value)
 
const G4String & GetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair
< G4double, G4double > 		GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double > 		GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 49 of file G4AntiNuclElastic.hh.

Constructor & Destructor Documentation

G4AntiNuclElastic::G4AntiNuclElastic ( )

Definition at line 54 of file G4AntiNuclElastic.cc.

References G4Alpha::Alpha(), G4AntiAlpha::AntiAlpha(), G4AntiDeuteron::AntiDeuteron(), G4AntiHe3::AntiHe3(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiTriton::AntiTriton(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), and G4Proton::Proton().

55  : G4HadronElastic("AntiAElastic")
56 {
57  //V.Ivanchenko commented out
58  //SetMinEnergy( 0.1*GeV );
59  //SetMaxEnergy( 10.*TeV );
60 
61 
62  theAProton = G4AntiProton::AntiProton();
63  theANeutron = G4AntiNeutron::AntiNeutron();
64  theADeuteron = G4AntiDeuteron::AntiDeuteron();
65  theATriton = G4AntiTriton::AntiTriton();
66  theAAlpha = G4AntiAlpha::AntiAlpha();
67  theAHe3 = G4AntiHe3::AntiHe3();
68 
69  theProton = G4Proton::Proton();
70  theNeutron = G4Neutron::Neutron();
71  theDeuteron = G4Deuteron::Deuteron();
72  theAlpha = G4Alpha::Alpha();
73 
74 
75  cs = new G4ComponentAntiNuclNuclearXS();
76  fParticle = 0;
77  fWaveVector = 0.;
78  fBeta = 0.;
79  fZommerfeld = 0.;
80  fAm = 0.;
81  fTetaCMS = 0.;
82  fRa = 0.;
83  fRef = 0.;
84  fceff = 0.;
85  fptot = 0.;
86  fTmax = 0.;
87  fThetaLab = 0.;
88 }
G4AntiHe3::AntiHe3
static G4AntiHe3 * AntiHe3()
Definition: G4AntiHe3.cc:94
G4HadronElastic::G4HadronElastic
G4HadronElastic(const G4String &name="hElasticLHEP")
Definition: G4HadronElastic.cc:45
G4AntiDeuteron::AntiDeuteron
static G4AntiDeuteron * AntiDeuteron()
Definition: G4AntiDeuteron.cc:94
G4AntiAlpha::AntiAlpha
static G4AntiAlpha * AntiAlpha()
Definition: G4AntiAlpha.cc:89
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4AntiTriton::AntiTriton
static G4AntiTriton * AntiTriton()
Definition: G4AntiTriton.cc:94
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4AntiNuclElastic::~G4AntiNuclElastic ( )
virtual

Definition at line 91 of file G4AntiNuclElastic.cc.

92 {
93  delete cs;
94 }

Member Function Documentation

G4double G4AntiNuclElastic::BesselJone ( G4double  z)

Definition at line 582 of file G4AntiNuclElastic.cc.

Referenced by BesselOneByArg().

583 {
584  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
585 
586  modvalue = std::fabs(value);
587 
588  if ( modvalue < 8.0 )
589  {
590  value2 = value*value;
591  fact1 = value*(72362614232.0 + value2*(-7895059235.0
592  + value2*( 242396853.1
593  + value2*(-2972611.439
594  + value2*( 15704.48260
595  + value2*(-30.16036606 ) ) ) ) ) );
596 
597  fact2 = 144725228442.0 + value2*(2300535178.0
598  + value2*(18583304.74
599  + value2*(99447.43394
600  + value2*(376.9991397
601  + value2*1.0 ) ) ) );
602  bessel = fact1/fact2;
603  }
604  else
605  {
606  arg = 8.0/modvalue;
607  value2 = arg*arg;
608 
609  shift = modvalue - 2.356194491;
610 
611  fact1 = 1.0 + value2*( 0.183105e-2
612  + value2*(-0.3516396496e-4
613  + value2*(0.2457520174e-5
614  + value2*(-0.240337019e-6 ) ) ) );
615 
616  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
617  + value2*( 0.8449199096e-5
618  + value2*(-0.88228987e-6
619  + value2*0.105787412e-6 ) ) );
620 
621  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
622  if (value < 0.0) bessel = -bessel;
623  }
624  return bessel;
625 }
value
const XML_Char int const XML_Char * value
Definition: include/expat.h:331
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::BesselJzero ( G4double  z)

Definition at line 531 of file G4AntiNuclElastic.cc.

Referenced by SampleInvariantT().

532 {
533  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
534 
535  modvalue = std::fabs(value);
536 
537  if ( value < 8.0 && value > -8.0 )
538  {
539  value2 = value*value;
540 
541  fact1 = 57568490574.0 + value2*(-13362590354.0
542  + value2*( 651619640.7
543  + value2*(-11214424.18
544  + value2*( 77392.33017
545  + value2*(-184.9052456 ) ) ) ) );
546 
547  fact2 = 57568490411.0 + value2*( 1029532985.0
548  + value2*( 9494680.718
549  + value2*(59272.64853
550  + value2*(267.8532712
551  + value2*1.0 ) ) ) );
552 
553  bessel = fact1/fact2;
554  }
555  else
556  {
557  arg = 8.0/modvalue;
558 
559  value2 = arg*arg;
560 
561  shift = modvalue-0.785398164;
562 
563  fact1 = 1.0 + value2*(-0.1098628627e-2
564  + value2*(0.2734510407e-4
565  + value2*(-0.2073370639e-5
566  + value2*0.2093887211e-6 ) ) );
567  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
568  + value2*(-0.6911147651e-5
569  + value2*(0.7621095161e-6
570  - value2*0.934945152e-7 ) ) );
571 
572  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
573  }
574  return bessel;
575 }
value
const XML_Char int const XML_Char * value
Definition: include/expat.h:331
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::BesselOneByArg ( G4double  z)

Definition at line 629 of file G4AntiNuclElastic.cc.

References BesselJone(), and test::x.

Referenced by SampleInvariantT().

630 {
631  G4double x2, result;
632 
633  if( std::fabs(x) < 0.01 )
634  {
635  x *= 0.5;
636  x2 = x*x;
637  result = (2.- x2 + x2*x2/6.)/4.;
638  }
639  else
640  {
641  result = BesselJone(x)/x;
642  }
643  return result;
644 }
G4AntiNuclElastic::BesselJone
G4double BesselJone(G4double z)
Definition: G4AntiNuclElastic.cc:582
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::CalculateAm ( G4double  momentum,
G4double  n,
G4double  Z 
)

Definition at line 515 of file G4AntiNuclElastic.cc.

References G4Pow::A13(), python.hepunit::Bohr_radius, G4Pow::GetInstance(), python.hepunit::hbarc, and n.

Referenced by SampleInvariantT().

516 {
517  G4double k = momentum/hbarc;
518  G4double ch = 1.13 + 3.76*n*n;
519  G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
520  G4double zn2 = zn*zn;
521  fAm = ch/zn2;
522 
523  return fAm;
524 }
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4Pow::A13
G4double A13(G4double A) const
Definition: G4Pow.hh:134
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::CalculateParticleBeta ( const G4ParticleDefinition *  particle,
G4double  momentum 
)

Definition at line 492 of file G4AntiNuclElastic.cc.

References test::a, and G4ParticleDefinition::GetPDGMass().

Referenced by SampleInvariantT().

494 {
495  G4double mass = particle->GetPDGMass();
496  G4double a = momentum/mass;
497  fBeta = a/std::sqrt(1+a*a);
498 
499  return fBeta;
500 }
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::CalculateZommerfeld ( G4double  beta,
G4double  Z1,
G4double  Z2 
)

Definition at line 506 of file G4AntiNuclElastic.cc.

References python.hepunit::fine_structure_const.

Referenced by SampleInvariantT().

507 {
508  fZommerfeld = fine_structure_const*Z1*Z2/beta;
509 
510  return fZommerfeld;
511 }
G4double G4AntiNuclElastic::DampFactor ( G4double  z)

Definition at line 472 of file G4AntiNuclElastic.cc.

Referenced by SampleInvariantT().

473 {
474  G4double df;
475  G4double f3 = 6.; // first factorials
476 
477  if( std::fabs(x) < 0.01 )
478  {
479  df=1./(1.+x*x/f3);
480  }
481  else
482  {
483  df = x/std::sinh(x);
484  }
485  return df;
486 }
G4double
double G4double
Definition: G4Types.hh:76
G4ComponentAntiNuclNuclearXS * G4AntiNuclElastic::GetComponentCrossSection ( )
inline

Definition at line 125 of file G4AntiNuclElastic.hh.

Referenced by DMXPhysicsList::ConstructHad(), G4HadronHElasticPhysics::ConstructProcess(), G4HadronDElasticPhysics::ConstructProcess(), and G4HadronElasticPhysics::ConstructProcess().

126 {
127  return cs;
128 }
G4double G4AntiNuclElastic::GetcosTeta1 ( G4double  plab,
G4int  A 
)

Definition at line 648 of file G4AntiNuclElastic.cc.

References python.hepunit::fermi, G4Pow::GetInstance(), python.hepunit::hbarc, and G4Pow::Z23().

Referenced by SampleInvariantT().

649 {
650 
651 // G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
652  G4double p0 = 1.*hbarc/fermi;
653 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
654  G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
655 //////////////////
656  if(cteta1 < -1.) cteta1 = -1.0;
657  return cteta1;
658 }
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4Pow::Z23
G4double Z23(G4int Z) const
Definition: G4Pow.hh:153
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 98 of file G4AntiNuclElastic.cc.

References BesselJzero(), BesselOneByArg(), CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), DampFactor(), CLHEP::HepLorentzVector::e(), energy(), G4UniformRand, G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc(), G4ParticleDefinition::GetBaryonNumber(), GetcosTeta1(), G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection(), G4Pow::GetInstance(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection(), python.hepunit::GeV, G4He3::He3(), CLHEP::Hep3Vector::mag(), CLHEP::HepLorentzVector::mag2(), python.hepunit::MeV, python.hepunit::millibarn, n, python.hepunit::pi, CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setPx(), CLHEP::HepLorentzVector::setPy(), CLHEP::HepLorentzVector::setPz(), sqr(), G4Triton::Triton(), python.hepunit::twopi, CLHEP::HepLorentzVector::vect(), test::x, and G4Pow::Z13().

Referenced by SampleThetaCMS(), and SampleThetaLab().

100 {
101  G4double T;
102  G4double Mproj = particle->GetPDGMass();
103  G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj));
104  G4double ctet1 = GetcosTeta1(Plab, A);
105 
106  G4double energy=Pproj.e()-Mproj;
107 
108  const G4ParticleDefinition* theParticle = particle;
109 
110  G4ParticleDefinition * theDef = 0;
111 
112  if(Z == 1 && A == 1) theDef = theProton;
113  else if (Z == 1 && A == 2) theDef = theDeuteron;
114  else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
115  else if (Z == 2 && A == 3) theDef = G4He3::He3();
116  else if (Z == 2 && A == 4) theDef = theAlpha;
117 
118 
119  G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z);
120 
121  //transform to CMS
122 
123  G4LorentzVector lv(0.0,0.0,0.0,TargMass);
124  lv += Pproj;
125  G4double S = lv.mag2()/GeV/GeV;
126 
127  G4ThreeVector bst = lv.boostVector();
128  Pproj.boost(-bst);
129 
130  G4ThreeVector p1 = Pproj.vect();
131  G4double ptot = p1.mag();
132 
133  fbst = bst;
134  fptot= ptot;
135  fTmax = 4.0*ptot*ptot;
136 
137  if(Plab/std::abs(particle->GetBaryonNumber()) < 100.*MeV) // Uzhi 24 Nov. 2011
138  {return fTmax*G4UniformRand();} // Uzhi 24 Nov. 2011
139 
140  G4double Z1 = particle->GetPDGCharge();
141  G4double Z2 = Z;
142 
143  G4double beta = CalculateParticleBeta(particle, ptot);
144  G4double n = CalculateZommerfeld( beta, Z1, Z2 );
145  G4double Am = CalculateAm( ptot, n, Z2 );
146  fWaveVector = ptot; // /hbarc;
147 
148  G4LorentzVector Fproj(0.,0.,0.,0.);
149  G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
150  XsCoulomb=XsCoulomb*0.38938e+6;
151 
152 
153  G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A);
154  G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A);
155 
156 
157  XsElastHad/=millibarn; XstotalHad/=millibarn;
158 
159 
160  G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad);
161 
162 // G4cout<<" XselastHadron " << XsElastHad << " XsCol "<< XsCoulomb <<G4endl;
163 // G4cout <<" XsTotal" << XstotalHad <<G4endl;
164 // G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl;
165 
166  if(G4UniformRand() < CoulombProb)
167  { // Simulation of Coulomb scattering
168 
169  G4double phi = twopi * G4UniformRand();
170  G4double Ksi = G4UniformRand();
171 
172  G4double par1 = 2.*(1.+Am)/(1.+ctet1);
173 
174 // ////sample ThetaCMS in Coulomb part
175 
176  G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
177 
178  G4double PtZ=ptot*cosThetaCMS;
179  Fproj.setPz(PtZ);
180  G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
181  G4double PtX= PtProjCMS * std::cos(phi);
182  G4double PtY= PtProjCMS * std::sin(phi);
183  Fproj.setPx(PtX);
184  Fproj.setPy(PtY);
185  Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
186  T = -(Pproj-Fproj).mag2();
187  } else
188 
189  {
190 ///////Simulation of strong interaction scattering////////////////////////////
191 
192 // G4double Qmax = 2.*ptot*197.33; // in fm^-1
193  G4double Qmax = 2.*3.0*197.33; // in fm^-1
194  G4double Amag = 70*70; // A1 in Magora funct:A1*exp(-q*A2)
195  G4double SlopeMag = 2.*3.0; // A2 in Magora funct:A1*exp(-q*A2)
196 
197  G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy);
198 
199  fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
200  0.227/G4Pow::GetInstance()->Z13(A);
201  if(A == 3) fRa=1.81;
202  if(A == 4) fRa=1.37;
203 
204  if((A>=12.) && (A<27) ) fRa=fRa*0.85;
205  if((A>=27.) && (A<48) ) fRa=fRa*0.90;
206  if((A>=48.) && (A<65) ) fRa=fRa*0.95;
207 
208  G4double Ref2 = 0;
209  G4double ceff2 =0;
210  G4double rho = 0;
211  if ((theParticle == theAProton) || (theParticle == theANeutron))
212  {
213  if(theDef == theProton)
214  {
215 // G4double Mp2=sqr(theDef->GetPDGMass()/GeV );
216 
217 // change 30 October
218 
219  if(Plab < 610.)
220  { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
221  13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
222  if((Plab < 5500.)&&(Plab >= 610.) )
223  { rho = 0.22; }
224  if((Plab >= 5500.)&&(Plab < 12300.) )
225  { rho = -0.32; }
226  if( Plab >= 12300.)
227  { rho = 0.135-2.26/(std::sqrt(S)) ;}
228 
229  Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*std::log(S) ;
230  ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ;
231 
232 /*
233  Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55;
234  if(S>1000.) Ref2=0.62+0.02*std::log(S) ;
235  ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * std::log(S) ;
236  if(S>1000.) ceff2 = 0.005 * std::log(S) + 0.29;
237 */
238 
239  Ref2=Ref2*Ref2;
240  ceff2 = ceff2*ceff2;
241 
242  SlopeMag = 0.5; // Uzhi
243  Amag= 1.; // Uzhi
244  }
245 
246  if(Z>2)
247  { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
248  ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*std::exp(-0.03*sig_pbarp);
249  }
250  if( (Z==2)&&(A==4) )
251  { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
252  ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
253  }
254  if( (Z==1)&&(A==3) )
255  { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
256  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
257  }
258  if( (Z==2)&&(A==3) )
259  { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
260  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
261  }
262  if( (Z==1)&&(A==2) )
263  {
264  Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
265  ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
266  }
267  }
268 
269  if (theParticle == theADeuteron)
270  {
271  sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.);
272  Ref2 = XstotalHad/10./2./pi ;
273  if(Z>2)
274  {
275  ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * std::exp(-0.03*sig_pbarp);
276  }
277  if(theDef == theProton)
278  {
279  ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
280  }
281  if(theDef == theDeuteron)
282  {
283  ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * std::exp(-0.03*sig_pbarp);
284  }
285  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
286  {
287  ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
288  }
289  if(theDef == theAlpha)
290  {
291  ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
292  }
293  }
294 
295  if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
296  {
297  sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
298  Ref2 = XstotalHad/10./2./pi ;
299  if(Z>2)
300  {
301  ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*std::exp(-0.03*sig_pbarp);
302  }
303  if(theDef == theProton)
304  {
305  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
306  }
307  if(theDef == theDeuteron)
308  {
309  ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
310  }
311  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
312  {
313  ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * std::exp(-0.02*sig_pbarp);
314  }
315  if(theDef == theAlpha)
316  {
317  ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
318  }
319  }
320 
321 
322  if (theParticle == theAAlpha)
323  {
324  sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
325  Ref2 = XstotalHad/10./2./pi ;
326  if(Z>2)
327  {
328  ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * std::exp(-0.03*sig_pbarp);
329  }
330  if(theDef == theProton)
331  {
332  ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
333  }
334  if(theDef == theDeuteron)
335  {
336  ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
337  }
338  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
339  {
340  ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
341  }
342  if(theDef == theAlpha)
343  {
344  ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * std::exp(-0.03*sig_pbarp);
345  }
346  }
347 
348  fRef=std::sqrt(Ref2);
349  fceff = std::sqrt(ceff2);
350 // G4cout<<" Ref "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl;
351 
352 
353  G4double Q = 0.0 ;
354  G4double BracFunct;
355  do
356  {
357  Q = -std::log(1.-(1.- std::exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag;
358  G4double x = fRef * Q;
359  BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) )
360 * sqr(DampFactor(pi*fceff*Q))) /(Amag*std::exp(-SlopeMag*Q));
361 
362  BracFunct = BracFunct * Q * sqr(sqr(fRef));
363  }
364  while (G4UniformRand()>BracFunct);
365 
366  T= sqr(Q);
367  T*=3.893913e+4; // fm -> MeV^2
368  }
369 
370  G4double cosTet=1.0-T/(2.*ptot*ptot);
371  if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
372  if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
373  fTetaCMS=std::acos(cosTet);
374 
375  return T;
376 }
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentAntiNuclNuclearXS.cc:268
G4AntiNuclElastic::BesselOneByArg
G4double BesselOneByArg(G4double z)
Definition: G4AntiNuclElastic.cc:629
python.hepunit.millibarn
int millibarn
Definition: hepunit.py:40
G4AntiNuclElastic::CalculateParticleBeta
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4AntiNuclElastic.cc:492
G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition *aParticle, G4double kinEnergy)
Definition: G4ComponentAntiNuclNuclearXS.cc:291
energy
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4Pow::Z13
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
Definition: G4ComponentAntiNuclNuclearXS.cc:88
G4AntiNuclElastic::BesselJzero
G4double BesselJzero(G4double z)
Definition: G4AntiNuclElastic.cc:531
G4AntiNuclElastic::CalculateAm
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4AntiNuclElastic.cc:515
G4AntiNuclElastic::DampFactor
G4double DampFactor(G4double z)
Definition: G4AntiNuclElastic.cc:472
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4AntiNuclElastic::CalculateZommerfeld
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4AntiNuclElastic.cc:506
G4double
double G4double
Definition: G4Types.hh:76
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
CLHEP::Hep3Vector::mag
double mag() const
G4AntiNuclElastic::GetcosTeta1
G4double GetcosTeta1(G4double plab, G4int A)
Definition: G4AntiNuclElastic.cc:648
G4double G4AntiNuclElastic::SampleThetaCMS ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 380 of file G4AntiNuclElastic.cc.

References G4cout, G4endl, G4UniformRand, python.hepunit::GeV, SampleInvariantT(), and G4HadronicInteraction::verboseLevel.

382 {
383  G4double T;
384  T = SampleInvariantT( p, plab, Z, A);
385 
386  // NaN finder
387  if(!(T < 0.0 || T >= 0.0))
388  {
389  if (verboseLevel > 0)
390  {
391  G4cout << "G4DiffuseElastic:WARNING: A = " << A
392  << " mom(GeV)= " << plab/GeV
393  << " S-wave will be sampled"
394  << G4endl;
395  }
396  T = G4UniformRand()*fTmax;
397 
398  }
399 
400  if(fptot > 0.) // Uzhi 24 Nov. 2011
401  {
402  G4double cosTet=1.0-T/(2.*fptot*fptot);
403  if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
404  if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
405  fTetaCMS=std::acos(cosTet);
406  return fTetaCMS;
407  } else // Uzhi 24 Nov. 2011
408  { // Uzhi 24 Nov. 2011
409  return 2.*G4UniformRand()-1.; // Uzhi 24 Nov. 2011
410  } // Uzhi 24 Nov. 2011
411 }
G4AntiNuclElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4AntiNuclElastic.cc:98
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AntiNuclElastic::SampleThetaLab ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 416 of file G4AntiNuclElastic.cc.

References CLHEP::HepLorentzVector::boost(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), python.hepunit::GeV, SampleInvariantT(), CLHEP::Hep3Vector::theta(), python.hepunit::twopi, test::v, G4HadronicInteraction::verboseLevel, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

418 {
419  G4double T;
420  T = SampleInvariantT( p, plab, Z, A);
421 
422  // NaN finder
423  if(!(T < 0.0 || T >= 0.0))
424  {
425  if (verboseLevel > 0)
426  {
427  G4cout << "G4DiffuseElastic:WARNING: A = " << A
428  << " mom(GeV)= " << plab/GeV
429  << " S-wave will be sampled"
430  << G4endl;
431  }
432  T = G4UniformRand()*fTmax;
433  }
434 
435  G4double phi = G4UniformRand()*twopi;
436 
437  G4double cost(1.);
438  if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} // Uzhi 24 Nov. 2011
439 
440  G4double sint;
441  if( cost >= 1.0 )
442  {
443  cost = 1.0;
444  sint = 0.0;
445  }
446  else if( cost <= -1.0)
447  {
448  cost = -1.0;
449  sint = 0.0;
450  }
451  else
452  {
453  sint = std::sqrt((1.0-cost)*(1.0+cost));
454  }
455 
456  G4double m1 = p->GetPDGMass();
457  G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
458  v *= fptot;
459  G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
460 
461  nlv.boost(fbst);
462 
463  G4ThreeVector np = nlv.vect();
464  G4double theta = np.theta();
465  fThetaLab = theta;
466 
467  return theta;
468 }
G4AntiNuclElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4AntiNuclElastic.cc:98
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::Hep3Vector::theta
double theta() const
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:11 for Geant4.10 by   doxygen 1.8.7