Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hhElastic Class Reference

#include <G4hhElastic.hh>

Inheritance diagram for G4hhElastic:
Collaboration diagram for G4hhElastic:

Public Member Functions

 G4hhElastic ()
 
 G4hhElastic (G4ParticleDefinition *target, G4ParticleDefinition *projectile, G4double plab)
 
 G4hhElastic (G4ParticleDefinition *target, G4ParticleDefinition *projectile)
 
virtual ~G4hhElastic ()
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
void Initialise ()
 
void BuildTableT (G4ParticleDefinition *target, G4ParticleDefinition *projectile)
 
void BuildTableTest (G4ParticleDefinition *target, G4ParticleDefinition *projectile, G4double plab)
 
G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int, G4int)
 
G4double SampleBisectionalT (const G4ParticleDefinition *p, G4double plab)
 
G4double SampleTest (G4double tMin)
 
G4double GetTransfer (G4int iMomentum, G4int iTransfer, G4double position)
 
void SetParameters ()
 
void SetSigmaTot (G4double stot)
 
void SetSpp (G4double spp)
 
G4double GetSpp ()
 
void SetParametersCMS (G4double plab)
 
G4double GetBq ()
 
G4double GetBQ ()
 
G4double GetBqQ ()
 
void SetBq (G4double b)
 
void SetBQ (G4double b)
 
void SetBqQ (G4double b)
 
G4double GetRhoReIm ()
 
void CalculateBQ (G4double b)
 
void CalculateBqQ13 (G4double b)
 
void CalculateBqQ12 (G4double b)
 
void CalculateBqQ123 (G4double b)
 
void SetRA (G4double rn, G4double pq, G4double pQ)
 
void SetRB (G4double rn, G4double pq, G4double pQ)
 
void SetAlphaP (G4double a)
 
void SetImCof (G4double a)
 
G4double GetImCof ()
 
void SetLambda (G4double L)
 
void SetEta (G4double E)
 
void SetCofF2 (G4double f)
 
void SetCofF3 (G4double f)
 
G4double GetCofF2 ()
 
G4double GetCofF3 ()
 
G4double GetRA ()
 
G4double GetRq ()
 
G4double GetRQ ()
 
G4double GetRB ()
 
G4double GetRg ()
 
G4double GetRG ()
 
G4complex Pomeron ()
 
G4complex Phi13 ()
 
G4complex Phi14 ()
 
G4complex Phi23 ()
 
G4complex Phi24 ()
 
G4complex GetF1qQgG (G4double qp)
 
G4double GetdsdtF1qQgG (G4double s, G4double q)
 
G4complex GetF2qQgG (G4double qp)
 
G4double GetdsdtF12qQgG (G4double s, G4double q)
 
G4complex GetF3qQgG (G4double qp)
 
G4double GetdsdtF123qQgG (G4double q)
 
G4double GetdsdtF13qQG (G4double s, G4double q)
 
G4complex GetAqq ()
 
G4complex GetAQQ ()
 
G4complex GetAqQ ()
 
G4double GetCofS1 ()
 
G4double GetCofS2 ()
 
G4double GetCofS3 ()
 
G4double GetOpticalRatio ()
 
G4complex GetF1 (G4double qp)
 
G4double GetdsdtF1 (G4double s, G4double q)
 
G4complex GetF2 (G4double qp)
 
G4double GetdsdtF12 (G4double s, G4double q)
 
G4complex GetF3 (G4double qp)
 
G4double GetdsdtF123 (G4double q)
 
G4double GetExpRatioF123 (G4double s, G4double q)
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
virtual G4HadFinalStateApplyYourself (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 ModelDescription (std::ostream &) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
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)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () 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
 
virtual const std::pair
< G4double, G4double
GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double
GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 

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 63 of file G4hhElastic.hh.

Constructor & Destructor Documentation

G4hhElastic::G4hhElastic ( )

Definition at line 73 of file G4hhElastic.cc.

74  : G4HadronElastic("HadrHadrElastic")
75 {
76  SetMinEnergy( 1.*GeV );
77  SetMaxEnergy( 10000.*TeV );
78  verboseLevel = 0;
79  lowEnergyRecoilLimit = 100.*keV;
80  lowEnergyLimitQ = 0.0*GeV;
81  lowEnergyLimitHE = 0.0*GeV;
82  lowestEnergyLimit= 0.0*keV;
83  plabLowLimit = 20.0*MeV;
84 
85  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
86  fInTkin=0;
87  theProton = G4Proton::Proton();
88  theNeutron = G4Neutron::Neutron();
89  thePionPlus = G4PionPlus::PionPlus();
90  thePionMinus= G4PionMinus::PionMinus();
91 
92  fTarget = G4Proton::Proton();
93  fProjectile = 0;
94  fHadrNuclXsc = new G4HadronNucleonXsc();
95 
96  fEnergyBin = 200;
97  fBinT = 514; // 514; // 500; // 200;
98 
99  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
100 
101  fTableT = 0;
102  fOldTkin = 0.;
103  SetParameters();
104 
105  Initialise();
106 }
void SetParameters()
Definition: G4hhElastic.hh:273
G4HadronElastic(const G4String &name="hElasticLHEP")
static constexpr double TeV
Definition: G4SIunits.hh:218
void SetMinEnergy(G4double anEnergy)
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetMaxEnergy(const G4double anEnergy)
static constexpr double MeV
Definition: G4SIunits.hh:214
void Initialise()
Definition: G4hhElastic.cc:230
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4hhElastic::G4hhElastic ( G4ParticleDefinition target,
G4ParticleDefinition projectile,
G4double  plab 
)

Definition at line 114 of file G4hhElastic.cc.

115  : G4HadronElastic("HadrHadrElastic")
116 {
117  SetMinEnergy( 1.*GeV );
118  SetMaxEnergy( 10000.*TeV );
119  verboseLevel = 0;
120  lowEnergyRecoilLimit = 100.*keV;
121  lowEnergyLimitQ = 0.0*GeV;
122  lowEnergyLimitHE = 0.0*GeV;
123  lowestEnergyLimit = 0.0*keV;
124  plabLowLimit = 20.0*MeV;
125 
126  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
127  fInTkin=0;
128  theProton = G4Proton::Proton();
129  theNeutron = G4Neutron::Neutron();
130  thePionPlus = G4PionPlus::PionPlus();
131  thePionMinus= G4PionMinus::PionMinus();
132 
133  fTarget = target;
134  fProjectile = projectile;
135  fMassTarg = fTarget->GetPDGMass();
136  fMassProj = fProjectile->GetPDGMass();
137  fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
138  fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
139  fHadrNuclXsc = new G4HadronNucleonXsc();
140 
141  fEnergyBin = 200;
142  fBinT = 514; // 200;
143 
144  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
145  fTableT = 0;
146  fOldTkin = 0.;
147 
148 
149  SetParameters();
150  SetParametersCMS( plab);
151 }
void SetParameters()
Definition: G4hhElastic.hh:273
void SetParametersCMS(G4double plab)
Definition: G4hhElastic.hh:318
const XML_Char * target
Definition: expat.h:268
G4HadronElastic(const G4String &name="hElasticLHEP")
static constexpr double TeV
Definition: G4SIunits.hh:218
void SetMinEnergy(G4double anEnergy)
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetPDGMass() const
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetMaxEnergy(const G4double anEnergy)
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4hhElastic::G4hhElastic ( G4ParticleDefinition target,
G4ParticleDefinition projectile 
)

Definition at line 159 of file G4hhElastic.cc.

160  : G4HadronElastic("HadrHadrElastic")
161 {
162  SetMinEnergy( 1.*GeV );
163  SetMaxEnergy( 10000.*TeV );
164  verboseLevel = 0;
165  lowEnergyRecoilLimit = 100.*keV;
166  lowEnergyLimitQ = 0.0*GeV;
167  lowEnergyLimitHE = 0.0*GeV;
168  lowestEnergyLimit= 0.0*keV;
169  plabLowLimit = 20.0*MeV;
170 
171  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
172  fInTkin=0;
173 
174  fTarget = target; // later vmg
175  fProjectile = projectile;
176  theProton = G4Proton::Proton();
177  theNeutron = G4Neutron::Neutron();
178  thePionPlus = G4PionPlus::PionPlus();
179  thePionMinus= G4PionMinus::PionMinus();
180 
181  fTarget = G4Proton::Proton(); // later vmg
182  fMassTarg = fTarget->GetPDGMass();
183  fMassProj = fProjectile->GetPDGMass();
184  fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
185  fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
186  fHadrNuclXsc = new G4HadronNucleonXsc();
187 
188  fEnergyBin = 200;
189  fBinT = 514; // 514; // 500; // 200;
190 
191  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
192 
193  fTableT = 0;
194  fOldTkin = 0.;
195 
196  SetParameters();
197 }
void SetParameters()
Definition: G4hhElastic.hh:273
const XML_Char * target
Definition: expat.h:268
G4HadronElastic(const G4String &name="hElasticLHEP")
static constexpr double TeV
Definition: G4SIunits.hh:218
void SetMinEnergy(G4double anEnergy)
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetPDGMass() const
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetMaxEnergy(const G4double anEnergy)
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4hhElastic::~G4hhElastic ( )
virtual

Definition at line 205 of file G4hhElastic.cc.

206 {
207  if ( fEnergyVector ) {
208  delete fEnergyVector;
209  fEnergyVector = 0;
210  }
211 
212  for ( std::vector<G4PhysicsTable*>::iterator it = fBankT.begin();
213  it != fBankT.end(); ++it ) {
214  if ( (*it) ) (*it)->clearAndDestroy();
215  delete *it;
216  *it = 0;
217  }
218  fTableT = 0;
219  if(fHadrNuclXsc) delete fHadrNuclXsc;
220 }

Member Function Documentation

void G4hhElastic::BuildTableT ( G4ParticleDefinition target,
G4ParticleDefinition projectile 
)

Definition at line 254 of file G4hhElastic.cc.

255 {
256  G4int iTkin, jTransfer;
257  G4double plab, Tkin, tMax;
258  G4double t1, t2, dt, delta = 0., sum = 0.;
259 
260  fTarget = target;
261  fProjectile = projectile;
262  fMassTarg = fTarget->GetPDGMass();
263  fMassProj = fProjectile->GetPDGMass();
264  fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
265  fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
266 
268  // G4HadronNucleonXsc* hnXsc = new G4HadronNucleonXsc();
269  fTableT = new G4PhysicsTable(fEnergyBin);
270 
271  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
272  {
273  Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
274  plab = std::sqrt( Tkin*( Tkin + 2*fMassProj ) );
275  // G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(projectile,
276  // G4ParticleMomentum(0.,0.,1.),
277  // Tkin);
278  // fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, target );
279 
280  SetParametersCMS( plab );
281 
282  tMax = 4.*fPcms*fPcms;
283  if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
284 
285  G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
286  sum = 0.;
287  dt = tMax/fBinT;
288 
289  // for(j = 1; j < fBinT; j++)
290 
291  for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer--)
292  {
293  t1 = dt*(jTransfer-1);
294  t2 = t1 + dt;
295 
296  if( fMassProj > 900.*MeV ) // pp, pn
297  {
298  delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
299  // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
300  }
301  else // pi+-p, K+-p
302  {
303  delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
304  // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
305  }
306  sum += delta;
307  vectorT->PutValue( jTransfer-1, t1, sum ); // t2
308  }
309  // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
310  fTableT->insertAt( iTkin, vectorT );
311  // delete theDynamicParticle;
312  }
313  // delete hnXsc;
314 
315  return;
316 }
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
void SetParametersCMS(G4double plab)
Definition: G4hhElastic.hh:318
const XML_Char * target
Definition: expat.h:268
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
G4double GetdsdtF123qQgG(G4double q)
Definition: G4hhElastic.hh:748
G4double GetdsdtF123(G4double q)
G4double GetPDGMass() const
tuple t1
Definition: plottest35.py:33
static constexpr double GeV
Definition: G4SIunits.hh:217
void insertAt(size_t, G4PhysicsVector *)
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4hhElastic::BuildTableTest ( G4ParticleDefinition target,
G4ParticleDefinition projectile,
G4double  plab 
)

Definition at line 520 of file G4hhElastic.cc.

521 {
522  G4int jTransfer;
523  G4double tMax; // , sQq, sQG;
524  G4double t1, t2, dt, delta = 0., sum = 0. ; // , threshold;
525 
526  fTarget = target;
527  fProjectile = projectile;
528  fMassTarg = fTarget->GetPDGMass();
529  fMassProj = fProjectile->GetPDGMass();
530  fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
531  fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
532  fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
533  fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
534 
535  G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fMassProj = "<<fMassProj<<" MeV"<<G4endl;
536  tMax = 4.*fPcms*fPcms;
537  if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
538 
539 
541  fTableT = new G4PhysicsTable(1);
542  G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
543 
544  sum = 0.;
545  dt = tMax/G4double(fBinT);
546  G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; fPcms = "<<fPcms/GeV
547  <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt = "<<dt/GeV/GeV<<" GeV2"<<G4endl;
548 
549  // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fRB*GeV<<G4endl;
550 
551  // for(jTransfer = 1; jTransfer < fBinT; jTransfer++)
552  for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer-- )
553  {
554  t1 = dt*(jTransfer-1);
555  t2 = t1 + dt;
556 
557  if( fMassProj > 900.*MeV ) // pp, pn
558  {
559  delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
560  // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, tMax);
561  }
562  else // pi+-p, K+-p
563  {
564  delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
565  // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, tMax);
566  // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
567  }
568  sum += delta;
569  // G4cout<<delta<<"\t"<<sum<<"\t"<<threshold<<G4endl;
570 
571  // sQq = GetdsdtF123(q1);
572  // sQG = GetdsdtF123qQgG(q1);
573  // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/millibarn<<"\t"<<sQq*GeV*GeV/millibarn<<G4endl;
574  // G4cout<<"sum = "<<sum<<", ";
575 
576  vectorT->PutValue( jTransfer-1, t1, sum ); // t2
577  }
578  // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
579  fTableT->insertAt( 0, vectorT );
580  fBankT.push_back( fTableT ); // 0
581 
582  // for(jTransfer = 0; jTransfer < fBinT-1; jTransfer++)
583  // G4cout<<(*(*fTableT)(0))(jTransfer)/sum<<"\t\t"<<G4Pow::GetInstance()->powN(2.,-jTransfer)<<G4endl;
584 
585  return;
586 }
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
const XML_Char * target
Definition: expat.h:268
void PutValue(size_t index, G4double energy, G4double dataValue)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double GetdsdtF123qQgG(G4double q)
Definition: G4hhElastic.hh:748
G4double GetdsdtF123(G4double q)
G4double GetPDGMass() const
tuple t1
Definition: plottest35.py:33
static constexpr double GeV
Definition: G4SIunits.hh:217
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4hhElastic::CalculateBQ ( G4double  b)
inline

Definition at line 959 of file G4hhElastic.hh.

960 {
961  fBq = b1;
962  G4double s1 = GetCofS1();
963  G4double s2 = GetCofS2();
964  G4double s3 = GetCofS3();
965  G4double sqrtBq = std::sqrt(fBq);
966 
967  // cofs of the fBQ 3rd equation
968 
969  G4double a = s3*sqrtBq;
970  G4double b = s1*fBq - 1.;
971  G4double c = (s2*fBq - 2.)*sqrtBq;
972  G4double d = 1. - fBq;
973 
974  // cofs of the incomplete 3rd equation
975 
976  G4double p = c/a;
977  p -= b*b/a/a/3.;
978  G4double q = d/a;
979  q -= b*c/a/a/3.;
980  q += 2*b*b*b/a/a/a/27.;
981 
982  // cofs for the incomplete colutions
983 
984  G4double D = p*p*p/3./3./3.;
985  D += q*q/2./2.;
986  G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
987  G4complex A = std::pow(A1,1./3.);
988 
989  G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
990  G4complex B = std::pow(B1,1./3.);
991 
992  // roots of the incomplete 3rd equation
993 
994  G4complex y1 = A + B;
995  G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
996  G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
997 
998  G4complex x1 = y1 - b/a/3.;
999  G4complex x2 = y2 - b/a/3.;
1000  G4complex x3 = y3 - b/a/3.;
1001 
1002  G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = "<<real(x2)<<"; re_x3 = "<<real(x3)<<G4endl;
1003  G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = "<<imag(x2)<<"; im_x3 = "<<imag(x3)<<G4endl;
1004 
1005  G4double r1 = real(x1)*real(x1);
1006  G4double r2 = real(x2)*real(x2);
1007  G4double r3 = real(x3)*real(x3);
1008 
1009  if( r1 <= 1. && r1 >= 0. ) fBQ = r1;
1010  else if( r2 <= 1. && r2 >= 0. ) fBQ = r2;
1011  else if( r3 <= 1. && r3 >= 0. ) fBQ = r3;
1012  else fBQ = 1.;
1013  // fBQ = real(x3)*real(x3);
1014  G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
1015  fOptRatio = fBq + fBQ + 2.*sqrtBqBQ - 1.;
1016  fOptRatio /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
1017  G4cout<<"F123, fOptRatio = "<<fOptRatio<<G4endl;
1018 
1019  return ;
1020 }
G4double GetCofS1()
Definition: G4hhElastic.hh:905
G4double GetCofS3()
Definition: G4hhElastic.hh:931
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
double B(double temperature)
tuple b
Definition: test.py:12
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4double GetCofS2()
Definition: G4hhElastic.hh:918
double D(double temp)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ12 ( G4double  b)
inline

Definition at line 788 of file G4hhElastic.hh.

789 {
790  fBq = b1;
791 
792  G4complex z1324 = G4complex(1./8./CLHEP::pi,0.);
793  z1324 /= Phi13() + Phi24() + fLambda + fEta;
794  G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
795 
796  G4complex z1423 = G4complex(1./8./CLHEP::pi,0.);
797  z1423 /= Phi14() + Phi23() + fLambda + fEta;
798  G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
799 
800  fBQ = 1. - 2.*fBq;
801  fBQ /= 2. - fSigmaTot*fBq*(c1324+1423);
802 
803  G4double ratio = 2.*(fBq + fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423);
804  G4cout<<"ratio = "<<ratio<<G4endl;
805 
806  return ;
807 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
G4complex Phi13()
Definition: G4hhElastic.hh:522
static constexpr double hbarc
G4complex Phi23()
Definition: G4hhElastic.hh:542
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ123 ( G4double  b)
inline

Definition at line 814 of file G4hhElastic.hh.

815 {
816  fBq = b1;
817 
818  G4complex z1324 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
819  z1324 /= Phi13() + Phi24() + fLambda + fEta;
820  G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
821 
822  G4complex z1423 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
823  z1423 /= Phi14() + Phi23() + fLambda + fEta;
824  G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
825 
826  G4complex z1314 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
827  z1314 /= Phi13() + Phi14() + fEta;
828  G4double c1314 = real(z1314)/(CLHEP::hbarc*CLHEP::hbarc);
829 
830  G4complex z2324 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
831  z2324 /= Phi23() + Phi24() + fEta;
832  G4double c2324 = real(z2324)/(CLHEP::hbarc*CLHEP::hbarc);
833 
834  G4complex z1323 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
835  z1323 /= Phi13() + Phi23() + fLambda;
836  G4double c1323 = real(z1323)/(CLHEP::hbarc*CLHEP::hbarc);
837 
838  G4complex z1424 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
839  z1424 /= Phi14() + Phi24() + fLambda;
840  G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
841 
842  G4double A = fSigmaTot*c2324;
843  G4double B = fSigmaTot*fBq*(c1324 + c1423 + c1323 + c1424) - 2.;
844  G4double C = 1. + fSigmaTot*fBq*fBq*c1314 - 2*fBq;
845  G4cout<<"A = "<<A<<"; B = "<<B<<"; C = "<<C<<G4endl;
846  G4cout<<"determinant = "<<B*B-4.*A*C<<G4endl;
847 
848  G4double x1 = ( -B - std::sqrt(B*B-4.*A*C) )/2./A;
849  G4double x2 = ( -B + std::sqrt(B*B-4.*A*C) )/2./A;
850  G4cout<<"x1 = "<<x1<<"; x2 = "<<x2<<G4endl;
851 
852  if( B*B-4.*A*C < 1.e-6 ) fBQ = std::abs(-B/2./A);
853  else if ( B < 0.) fBQ = std::abs( ( -B - std::sqrt(B*B-4.*A*C) )/2./A);
854  else fBQ = std::abs( ( -B + std::sqrt(B*B-4.*A*C) )/2./A);
855 
856  fOptRatio = 2*(fBq+fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423 + c1323 + c1424);
857  fOptRatio -= fSigmaTot*fBq*fBq*c1314 + fSigmaTot*c2324*fBQ*fBQ;
858  G4cout<<"BqQ123, fOptRatio = "<<fOptRatio<<G4endl;
859 
860  return ;
861 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
G4complex Phi13()
Definition: G4hhElastic.hh:522
static constexpr double hbarc
G4complex Phi23()
Definition: G4hhElastic.hh:542
double C(double temp)
double B(double temperature)
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
#define G4endl
Definition: G4ios.hh:61
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ13 ( G4double  b)
inline

Definition at line 765 of file G4hhElastic.hh.

766 {
767  fBQ = b2;
768 
769  G4complex z1424 = G4complex(1./8./CLHEP::pi,0.);
770  z1424 /= Phi14() + Phi24() + fAlpha;
771  G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
772 
773  fBqQ = 1. - fBQ;
774  fBQ /= 1. - fSigmaTot*fBQ*c1424;
775 
776  G4cout<<"fSigmaTot*fBQ*c1424 = "<<fSigmaTot*fBQ*c1424<<G4endl;
777 
778  G4double ratio = fBqQ + fBQ - fSigmaTot*fBqQ*fBQ*c1424;
779  G4cout<<"ratio = "<<ratio<<G4endl;
780 
781  return ;
782 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
static constexpr double hbarc
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4complex G4hhElastic::GetAqq ( )
inline

Definition at line 869 of file G4hhElastic.hh.

870 {
871  G4double re, im;
872 
873  re = fRq*fRq/8. + fAlphaP*G4Log(fSpp/fSo) + 8.*fLambda/9.;
874  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
875  return G4complex(re,im);
876 }
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetAQQ ( )
inline

Definition at line 882 of file G4hhElastic.hh.

883 {
884  G4double re, im;
885 
886  re = fRQ*fRQ/8. + fAlphaP*G4Log(fSpp/fSo) + 2.*fLambda/9.;
887  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
888  return G4complex(re,im);
889 }
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetAqQ ( )
inline

Definition at line 895 of file G4hhElastic.hh.

896 {
897  G4complex z = 0.5*( GetAqq() + GetAQQ() );
898  return z;
899 }
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAQQ()
Definition: G4hhElastic.hh:882
G4complex GetAqq()
Definition: G4hhElastic.hh:869
tuple z
Definition: test.py:28

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetBq ( )
inline

Definition at line 171 of file G4hhElastic.hh.

171 { return fBq;};
G4double G4hhElastic::GetBQ ( )
inline

Definition at line 172 of file G4hhElastic.hh.

172 { return fBQ;};
G4double G4hhElastic::GetBqQ ( )
inline

Definition at line 173 of file G4hhElastic.hh.

173 { return fBqQ;};
G4double G4hhElastic::GetCofF2 ( )
inline

Definition at line 193 of file G4hhElastic.hh.

193 {return fCofF2;};
G4double G4hhElastic::GetCofF3 ( )
inline

Definition at line 194 of file G4hhElastic.hh.

194 {return fCofF3;};
G4double G4hhElastic::GetCofS1 ( )
inline

Definition at line 905 of file G4hhElastic.hh.

906 {
907  G4complex z = 1./( GetAqQ() + 4.*fLambda/9. );
908  G4double result = real(z);
909  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
910  result *= fSigmaTot*fCofF2;
911  return result;
912 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
std::complex< G4double > G4complex
Definition: G4Types.hh:81
tuple z
Definition: test.py:28
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetCofS2 ( )
inline

Definition at line 918 of file G4hhElastic.hh.

919 {
920  G4complex z = 1./( GetAqq() + GetAqQ() - 4.*fLambda/9. );
921  G4double result = real(z);
922  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
923  result *= fSigmaTot*fCofF3;
924  return result;
925 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAqq()
Definition: G4hhElastic.hh:869
tuple z
Definition: test.py:28
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetCofS3 ( )
inline

Definition at line 931 of file G4hhElastic.hh.

932 {
933  G4complex z = 1./( GetAQQ() + GetAqQ() + 2.*fLambda/9. );
934  G4double result = real(z);
935  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
936  result *= fSigmaTot*fCofF3;
937  return result;
938 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAQQ()
Definition: G4hhElastic.hh:882
tuple z
Definition: test.py:28
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetdsdtF1 ( G4double  s,
G4double  q 
)
inline

Definition at line 1044 of file G4hhElastic.hh.

1045 {
1046  fSpp = spp;
1047  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1048  G4complex F1 = GetF1(t);
1049 
1050  G4double dsdt = CLHEP::pi/p/p;
1051  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1052  return dsdt;
1053 }
static constexpr double proton_mass_c2
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54
G4complex GetF1(G4double qp)

Here is the call graph for this function:

G4double G4hhElastic::GetdsdtF12 ( G4double  s,
G4double  q 
)
inline

Definition at line 1087 of file G4hhElastic.hh.

1088 {
1089  fSpp = spp;
1090  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1091  G4complex F1 = GetF1(t) - GetF2(t);
1092 
1093  G4double dsdt = CLHEP::pi/p/p;
1094  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1095  return dsdt;
1096 }
static constexpr double proton_mass_c2
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetF2(G4double qp)
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54
G4complex GetF1(G4double qp)

Here is the call graph for this function:

G4double G4hhElastic::GetdsdtF123 ( G4double  q)
inline

Definition at line 1130 of file G4hhElastic.hh.

1131 {
1132  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1133  G4complex F1 = GetF1(t);
1134  F1 -= fCofF2*GetF2(t);
1135  F1 -= fCofF3*GetF3(t);
1136  G4double dsdt = CLHEP::pi/p/p;
1137  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1138  return dsdt;
1139 }
static constexpr double proton_mass_c2
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetF2(G4double qp)
G4complex GetF3(G4double qp)
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54
G4complex GetF1(G4double qp)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetdsdtF123qQgG ( G4double  q)
inline

Definition at line 748 of file G4hhElastic.hh.

749 {
750  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp );
751 
752  G4complex F123 = GetF1qQgG(t); // - fCofF2*GetF2qQgG(t) - fCofF3*GetF3qQgG(t);
753  F123 -= fCofF2*GetF2qQgG(t);
754  F123 -= fCofF3*GetF3qQgG(t);
755 
756  G4double dsdt = CLHEP::pi/p/p;
757  dsdt *= real(F123)*real(F123) + imag(F123)*imag(F123);
758  return dsdt;
759 }
G4complex GetF3qQgG(G4double qp)
Definition: G4hhElastic.hh:694
G4complex GetF2qQgG(G4double qp)
Definition: G4hhElastic.hh:645
const char * p
Definition: xmltok.h:285
G4complex GetF1qQgG(G4double qp)
Definition: G4hhElastic.hh:564
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetdsdtF12qQgG ( G4double  s,
G4double  q 
)
inline

Definition at line 678 of file G4hhElastic.hh.

679 {
680  fSpp = spp;
681  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
682 
683  G4complex F12 = GetF1qQgG(t) - GetF2qQgG(t);
684 
685  G4double dsdt = CLHEP::pi/p/p;
686  dsdt *= real(F12)*real(F12) + imag(F12)*imag(F12);
687  return dsdt;
688 }
G4complex GetF2qQgG(G4double qp)
Definition: G4hhElastic.hh:645
const char * p
Definition: xmltok.h:285
G4complex GetF1qQgG(G4double qp)
Definition: G4hhElastic.hh:564
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define F12
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4double G4hhElastic::GetdsdtF13qQG ( G4double  s,
G4double  q 
)
inline

Definition at line 586 of file G4hhElastic.hh.

587 {
588  fSpp = spp;
589  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
590  G4double k = p/CLHEP::hbarc;
591 
592  G4complex exp14 = fBqQ*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
593  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
594 
595  G4complex F1 = exp14 + exp24;
596 
597  F1 *= 0.25*k*fSigmaTot/CLHEP::pi;
598  F1 *= G4complex(0.,1.);
599 
600  // 1424
601 
602  G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
603  z1424 /= Phi14() + Phi24() + fLambda;
604  z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
605 
606 
607  G4complex exp1424 = std::exp(-z1424*t);
608  exp1424 /= Phi14() + Phi24() + fLambda;
609 
610  G4complex F3 = fBqQ*fBQ*exp1424;
611 
612 
613  F3 *= 0.25*k/CLHEP::pi;
614  F3 *= G4complex(0.,1.);
615  F3 *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
616 
617  G4complex F13 = F1 - F3;
618 
619  G4double dsdt = CLHEP::pi/p/p;
620  dsdt *= real(F13)*real(F13) + imag(F13)*imag(F13);
621 
622  return dsdt;
623 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
#define F13
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4double G4hhElastic::GetdsdtF1qQgG ( G4double  s,
G4double  q 
)
inline

Definition at line 629 of file G4hhElastic.hh.

630 {
631  fSpp = spp;
632  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
633 
634  G4complex F1 = GetF1qQgG(t);
635 
636  G4double dsdt = CLHEP::pi/p/p;
637  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
638  return dsdt;
639 }
const char * p
Definition: xmltok.h:285
G4complex GetF1qQgG(G4double qp)
Definition: G4hhElastic.hh:564
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

G4double G4hhElastic::GetExpRatioF123 ( G4double  s,
G4double  q 
)
inline

Definition at line 1145 of file G4hhElastic.hh.

1146 {
1147  fSpp = spp;
1148  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1149 
1150  // qQ-ds/dt
1151 
1152  G4complex F1 = GetF1(t) - fCofF2*GetF2(t) - fCofF3*GetF3(t);
1153 
1154  G4double dsdt = CLHEP::pi/p/p;
1155  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1156 
1157  // exponent ds/dt
1158 
1159  G4complex F10 = GetF1(0.) - fCofF2*GetF2(0.) - fCofF3*GetF3(0.);
1160 
1161  fRhoReIm = real(F10)/imag(F10);
1162 
1163  G4double dsdt0 = CLHEP::pi/p/p;
1164  dsdt0 *= real(F10)*real(F10) + imag(F10)*imag(F10);
1165 
1166  dsdt0 *= G4Exp(-fExpSlope*t);
1167 
1168  G4double ratio = dsdt/dsdt0;
1169 
1170  return ratio;
1171 }
static constexpr double proton_mass_c2
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define F10
G4complex GetF2(G4double qp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4complex GetF3(G4double qp)
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54
G4complex GetF1(G4double qp)

Here is the call graph for this function:

G4complex G4hhElastic::GetF1 ( G4double  qp)
inline

Definition at line 1026 of file G4hhElastic.hh.

1027 {
1028  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1029  G4double k = p/CLHEP::hbarc;
1030  G4complex exp1 = fBq*std::exp(-GetAqq()*t);
1031  G4complex exp2 = fBQ*std::exp(-GetAQQ()*t);
1032  G4complex exp3 = 2.*std::sqrt(fBq*fBQ)*std::exp(-GetAqQ()*t);
1033 
1034  G4complex res = exp1 + exp2 + exp3;
1035  res *= 0.25*k*fSigmaTot/CLHEP::pi;
1036  res *= G4complex(0.,1.);
1037  return res;
1038 }
static constexpr double proton_mass_c2
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAQQ()
Definition: G4hhElastic.hh:882
G4complex GetAqq()
Definition: G4hhElastic.hh:869
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetF1qQgG ( G4double  qp)
inline

Definition at line 564 of file G4hhElastic.hh.

565 {
566  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
567  G4double k = p/CLHEP::hbarc;
568 
569  G4complex exp13 = fBq*std::exp(-(Phi13() + fBeta*fBeta*fLambda + fDelta*fDelta*fEta)*t);
570  G4complex exp14 = fBq*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
571  G4complex exp23 = fBQ*std::exp(-(Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta)*t);
572  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
573 
574  G4complex res = exp13 + exp14 + exp23 + exp24;
575 
576  res *= 0.25*k*fSigmaTot/CLHEP::pi;
577  res *= G4complex(0.,1.);
578 
579  return res;
580 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
G4complex Phi13()
Definition: G4hhElastic.hh:522
static constexpr double hbarc
G4complex Phi23()
Definition: G4hhElastic.hh:542
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetF2 ( G4double  qp)
inline

Definition at line 1059 of file G4hhElastic.hh.

1060 {
1061  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1062  G4double k = p/CLHEP::hbarc;
1063  G4complex z1 = GetAqq()*GetAQQ() - 16.*fLambda*fLambda/81.;
1064  z1 /= 2.*(GetAqQ() + 4.*fLambda/9.);
1065  G4complex exp1 = std::exp(-z1*t);
1066 
1067  G4complex z2 = 0.5*( GetAqQ() - 4.*fLambda/9.);
1068 
1069  G4complex exp2 = std::exp(-z2*t);
1070 
1071  G4complex res = exp1 + exp2;
1072 
1073  G4complex z3 = GetAqQ() + 4.*fLambda/9.;
1074 
1075  res *= 0.25*k/CLHEP::pi;
1076  res *= G4complex(0.,1.);
1077  res /= z3;
1078  res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1079 
1080  return res;
1081 }
static constexpr double proton_mass_c2
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAQQ()
Definition: G4hhElastic.hh:882
G4complex GetAqq()
Definition: G4hhElastic.hh:869
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetF2qQgG ( G4double  qp)
inline

Definition at line 645 of file G4hhElastic.hh.

646 {
647  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
648  G4double k = p/CLHEP::hbarc;
649 
650  G4complex z1324 = -(Phi24() + fAlpha*fLambda + fGamma*fEta)*(Phi24() + fAlpha*fLambda + fGamma*fEta);
651  z1324 /= Phi13() + Phi24() + fLambda + fEta;
652  z1324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
653 
654  G4complex exp1324 = std::exp(-z1324*t);
655  exp1324 /= Phi13() + Phi24() + fLambda + fEta;
656 
657  G4complex z1423 = -(Phi23() + fAlpha*fLambda + fDelta*fEta)*(Phi24() + fAlpha*fLambda + fDelta*fEta);;
658  z1423 /= Phi14() + Phi23() + fLambda + fEta;
659  z1423 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
660 
661  G4complex exp1423 = std::exp(-z1423*t);
662  exp1423 /= Phi14() + Phi23() + fLambda + fEta;
663 
664  G4complex res = exp1324 + exp1423;
665 
666 
667  res *= 0.25*k/CLHEP::pi;
668  res *= G4complex(0.,1.);
669  res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc); // or 4. ???
670 
671  return res;
672 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
G4complex Phi13()
Definition: G4hhElastic.hh:522
static constexpr double hbarc
G4complex Phi23()
Definition: G4hhElastic.hh:542
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetF3 ( G4double  qp)
inline

Definition at line 1102 of file G4hhElastic.hh.

1103 {
1104  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1105  G4double k = p/CLHEP::hbarc;
1106  G4complex z1 = GetAqq()*GetAqQ() - 4.*fLambda*fLambda/81.;
1107  z1 /= GetAqq() + GetAqQ() - 4.*fLambda/9.;
1108 
1109  G4complex exp1 = std::exp(-z1*t)*fBq/(GetAqq() + GetAqQ() - 4.*fLambda/9.);
1110 
1111  G4complex z2 = GetAqQ()*GetAQQ() - 1.*fLambda*fLambda/81.;
1112  z2 /= GetAQQ() + GetAqQ() + 2.*fLambda/9.;
1113 
1114  G4complex exp2 = std::exp(-z2*t)*fBQ/(GetAQQ() + GetAqQ() + 2.*fLambda/9.);
1115 
1116  G4complex res = exp1 + exp2;
1117 
1118 
1119  res *= 0.25*k/CLHEP::pi;
1120  res *= G4complex(0.,1.);
1121  res *= std::sqrt(fBq*fBQ)*fSigmaTot*fSigmaTot/(4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1122 
1123  return res;
1124 }
static constexpr double proton_mass_c2
G4complex GetAqQ()
Definition: G4hhElastic.hh:895
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex GetAQQ()
Definition: G4hhElastic.hh:882
G4complex GetAqq()
Definition: G4hhElastic.hh:869
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::GetF3qQgG ( G4double  qp)
inline

Definition at line 694 of file G4hhElastic.hh.

695 {
696  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
697  G4double k = p/CLHEP::hbarc;
698 
699  // 1314
700 
701  G4complex z1314 = -(Phi14() + fGamma*fEta)*(Phi14() + fGamma*fEta);
702  z1314 /= Phi13() + Phi14() + fEta;
703  z1314 += Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta;
704 
705  G4complex exp1314 = std::exp(-z1314*t);
706  exp1314 /= Phi13() + Phi14() + fEta;
707 
708  // 2324
709 
710  G4complex z2324 = -(Phi24() + fGamma*fEta)*(Phi24() + fGamma*fEta);;
711  z2324 /= Phi24() + Phi23() + fEta;
712  z2324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
713 
714  G4complex exp2324 = std::exp(-z2324*t);
715  exp2324 /= Phi24() + Phi23() + fEta;
716 
717  // 1323
718 
719  G4complex z1323 = -(Phi23() + fAlpha*fLambda)*(Phi23() + fAlpha*fLambda);
720  z1323 /= Phi13() + Phi23() + fLambda;
721  z1323 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
722 
723  G4complex exp1323 = std::exp(-z1323*t);
724  exp1323 /= Phi13() + Phi23() + fLambda;
725 
726  // 1424
727 
728  G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
729  z1424 /= Phi14() + Phi24() + fLambda;
730  z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
731 
732  G4complex exp1424 = std::exp(-z1424*t);
733  exp1424 /= Phi14() + Phi24() + fLambda;
734 
735  G4complex res = fBq*fBq*exp1314 + fBQ*fBQ*exp2324 + fBq*fBQ*exp1323 + fBq*fBQ*exp1424;
736 
737  res *= 0.25*k/CLHEP::pi;
738  res *= G4complex(0.,1.);
739  res *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
740 
741  return res;
742 }
G4complex Phi14()
Definition: G4hhElastic.hh:532
G4complex Phi13()
Definition: G4hhElastic.hh:522
static constexpr double hbarc
G4complex Phi23()
Definition: G4hhElastic.hh:542
const char * p
Definition: xmltok.h:285
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4complex Phi24()
Definition: G4hhElastic.hh:552
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::GetImCof ( )
inline

Definition at line 188 of file G4hhElastic.hh.

188 {return fImCof;};
G4double G4hhElastic::GetOpticalRatio ( )
inline

Definition at line 944 of file G4hhElastic.hh.

945 {
946  return fOptRatio;
947  // G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
948  // G4double result = fBq + fBQ + 2.*sqrtBqBQ - 1.;
949  // result /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
950  // return result;
951 }
G4double G4hhElastic::GetRA ( )
inline

Definition at line 196 of file G4hhElastic.hh.

196 { return fRA;};

Here is the caller graph for this function:

G4double G4hhElastic::GetRB ( )
inline

Definition at line 200 of file G4hhElastic.hh.

200 { return fRB;};

Here is the caller graph for this function:

G4double G4hhElastic::GetRg ( )
inline

Definition at line 201 of file G4hhElastic.hh.

201 { return fRg;};
G4double G4hhElastic::GetRG ( )
inline

Definition at line 202 of file G4hhElastic.hh.

202 { return fRG;};
G4double G4hhElastic::GetRhoReIm ( )
inline

Definition at line 177 of file G4hhElastic.hh.

177 { return fRhoReIm;};
G4double G4hhElastic::GetRq ( )
inline

Definition at line 197 of file G4hhElastic.hh.

197 { return fRq;};
G4double G4hhElastic::GetRQ ( )
inline

Definition at line 198 of file G4hhElastic.hh.

198 { return fRQ;};
G4double G4hhElastic::GetSpp ( )
inline

Definition at line 168 of file G4hhElastic.hh.

168 {return fSpp;};
G4double G4hhElastic::GetTransfer ( G4int  iMomentum,
G4int  iTransfer,
G4double  position 
)

Definition at line 629 of file G4hhElastic.cc.

630 {
631  G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
632 
633  if( iTransfer == 0 )
634  {
635  randTransfer = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
636  // iTransfer++;
637  }
638  else
639  {
640  if ( iTransfer >= G4int((*fTableT)(iTkin)->GetVectorLength()) )
641  {
642  iTransfer = (*fTableT)(iTkin)->GetVectorLength() - 1;
643  }
644  y1 = (*(*fTableT)(iTkin))(iTransfer-1);
645  y2 = (*(*fTableT)(iTkin))(iTransfer);
646 
647  x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
648  x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
649 
650  delta = y2 - y1;
651  mean = y2 + y1;
652 
653  if ( x1 == x2 ) randTransfer = x2;
654  else
655  {
656  // if ( y1 == y2 )
657  if ( delta < epsilon*mean )
658  randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
659  else randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
660  }
661  }
662  return randTransfer;
663 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
double epsilon(double density, double temperature)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4hhElastic::Initialise ( )

Definition at line 230 of file G4hhElastic.cc.

231 {
232  // pp,pn
233 
234  fProjectile = G4Proton::Proton();
235  BuildTableT(fTarget, fProjectile);
236  fBankT.push_back(fTableT); // 0
237 
238  // pi+-p
239 
240  fProjectile = G4PionPlus::PionPlus();
241  BuildTableT(fTarget, fProjectile);
242  fBankT.push_back(fTableT); // 1
243  //K+-p
244  fProjectile = G4KaonPlus::KaonPlus();
245  BuildTableT(fTarget, fProjectile);
246  fBankT.push_back(fTableT); // 2
247 
248 }
void BuildTableT(G4ParticleDefinition *target, G4ParticleDefinition *projectile)
Definition: G4hhElastic.cc:254
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4hhElastic::IsApplicable ( const G4HadProjectile projectile,
G4Nucleus nucleus 
)
inlinevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 258 of file G4hhElastic.hh.

260 {
261  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
262  projectile.GetDefinition() == G4Neutron::Neutron() ||
263  projectile.GetDefinition() == G4PionPlus::PionPlus() ||
264  projectile.GetDefinition() == G4PionMinus::PionMinus() ||
265  projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
266  projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
267 
268  nucleus.GetZ_asInt() < 2 ) return true;
269  else return false;
270 }
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
const G4ParticleDefinition * GetDefinition() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113

Here is the call graph for this function:

G4complex G4hhElastic::Phi13 ( )
inline

Definition at line 522 of file G4hhElastic.hh.

523 {
524  G4double re = (fRq*fRq + fRg*fRg)/16.;
525  G4complex result(re,0.);
526  result += Pomeron();
527  return result;
528 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex Pomeron()
Definition: G4hhElastic.hh:511
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::Phi14 ( )
inline

Definition at line 532 of file G4hhElastic.hh.

533 {
534  G4double re = (fRq*fRq + fRG*fRG)/16.;
535  G4complex result(re,0.);
536  result += Pomeron();
537  return result;
538 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex Pomeron()
Definition: G4hhElastic.hh:511
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::Phi23 ( )
inline

Definition at line 542 of file G4hhElastic.hh.

543 {
544  G4double re = (fRQ*fRQ + fRg*fRg)/16.;
545  G4complex result(re,0.);
546  result += Pomeron();
547  return result;
548 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex Pomeron()
Definition: G4hhElastic.hh:511
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::Phi24 ( )
inline

Definition at line 552 of file G4hhElastic.hh.

553 {
554  G4double re = (fRQ*fRQ + fRG*fRG)/16.;
555  G4complex result(re,0.);
556  result += Pomeron();
557  return result;
558 }
G4double G4ParticleHPJENDLHEData::G4double result
G4complex Pomeron()
Definition: G4hhElastic.hh:511
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4complex G4hhElastic::Pomeron ( )
inline

Definition at line 511 of file G4hhElastic.hh.

512 {
513  G4double re, im;
514 
515  re = fAlphaP*G4Log(fSpp/fSo);
516  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
517  return G4complex(re,im);
518 }
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4hhElastic::SampleBisectionalT ( const G4ParticleDefinition p,
G4double  plab 
)

Definition at line 439 of file G4hhElastic.cc.

440 {
441  G4int iTkin, iTransfer;
442  G4double t, position, m1 = aParticle->GetPDGMass();
443  G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
444 
445  if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
446  {
447  fTableT = fBankT[0];
448  }
449  if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
450  {
451  fTableT = fBankT[1];
452  }
453  if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
454  {
455  fTableT = fBankT[2];
456  }
457  G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
458  G4double deltaMax = 1.e-2;
459 
460  if ( delta < deltaMax ) iTkin = fInTkin;
461  else
462  {
463  for( iTkin = 0; iTkin < fEnergyBin; iTkin++ )
464  {
465  if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
466  }
467  }
468  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
469  if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
470 
471  fOldTkin = Tkin;
472  fInTkin = iTkin;
473 
474  if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
475  {
476  position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
477 
478  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
479  {
480  if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
481  }
482  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
483 
484  t = GetTransfer(iTkin, iTransfer, position);
485 
486 
487  }
488  else // Tkin inside between energy table edges
489  {
490  G4double rand = G4UniformRand();
491  position = (*(*fTableT)(iTkin))(0)*rand;
492 
493  //
494  // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBinT-2);
495  G4int sTransfer = 0, fTransfer = fBinT - 2, dTransfer = fTransfer - sTransfer;
496  G4double y2;
497 
498  for( iTransfer = 0; iTransfer < fBinT - 1; iTransfer++ )
499  {
500  // dTransfer %= 2;
501  dTransfer /= 2;
502  // dTransfer *= 0.5;
503  y2 = (*(*fTableT)(iTkin))( sTransfer + dTransfer );
504 
505  if( y2 > position ) sTransfer += dTransfer;
506 
507  // if( dTransfer <= 1 ) break;
508  if( dTransfer < 1 ) break;
509  }
510  t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sTransfer); // +(-0.5+rand)*(*fTableT)(iTkin)->GetLowEdgeEnergy(3);
511  }
512  return t;
513 }
int G4int
Definition: G4Types.hh:78
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4double GetTransfer(G4int iMomentum, G4int iTransfer, G4double position)
Definition: G4hhElastic.cc:629
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113

Here is the call graph for this function:

G4double G4hhElastic::SampleInvariantT ( const G4ParticleDefinition p,
G4double  plab,
G4int  ,
G4int   
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 322 of file G4hhElastic.cc.

324 {
325  G4int iTkin, iTransfer;
326  G4double t, t2, position, m1 = aParticle->GetPDGMass();
327  G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
328 
329  if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
330  {
331  fTableT = fBankT[0];
332  }
333  if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
334  {
335  fTableT = fBankT[1];
336  }
337  if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
338  {
339  fTableT = fBankT[2];
340  }
341 
342  G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
343  G4double deltaMax = 1.e-2;
344 
345  if ( delta < deltaMax ) iTkin = fInTkin;
346  else
347  {
348  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
349  {
350  if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
351  }
352  }
353  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
354  if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
355 
356  fOldTkin = Tkin;
357  fInTkin = iTkin;
358 
359  if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
360  {
361  position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
362 
363  // G4cout<<"position = "<<position<<G4endl;
364 
365  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
366  {
367  if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
368  }
369  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
370 
371  // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
372 
373  t = GetTransfer(iTkin, iTransfer, position);
374 
375  // G4cout<<"t = "<<t<<G4endl;
376  }
377  else // Tkin inside between energy table edges
378  {
379  // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
380  position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
381 
382  // G4cout<<"position = "<<position<<G4endl;
383 
384  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
385  {
386  // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
387  if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
388  }
389  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
390 
391  // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
392 
393  t2 = GetTransfer(iTkin, iTransfer, position);
394  return t2;
395  /*
396  G4double t1, E1, E2, W, W1, W2;
397  // G4cout<<"t2 = "<<t2<<G4endl;
398 
399  E2 = fEnergyVector->GetLowEdgeEnergy(iTkin);
400 
401  // G4cout<<"E2 = "<<E2<<G4endl;
402 
403  iTkin--;
404 
405  // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
406 
407  // G4cout<<"position = "<<position<<G4endl;
408 
409  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
410  {
411  // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
412  if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
413  }
414  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
415 
416  t1 = GetTransfer(iTkin, iTransfer, position);
417 
418  // G4cout<<"t1 = "<<t1<<G4endl;
419 
420  E1 = fEnergyVector->GetLowEdgeEnergy(iTkin);
421 
422  // G4cout<<"E1 = "<<E1<<G4endl;
423 
424  W = 1.0/(E2 - E1);
425  W1 = (E2 - Tkin)*W;
426  W2 = (Tkin - E1)*W;
427 
428  t = W1*t1 + W2*t2;
429  */
430  }
431  return t;
432 }
int G4int
Definition: G4Types.hh:78
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4double GetTransfer(G4int iMomentum, G4int iTransfer, G4double position)
Definition: G4hhElastic.cc:629
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113

Here is the call graph for this function:

G4double G4hhElastic::SampleTest ( G4double  tMin)

Definition at line 593 of file G4hhElastic.cc.

594 {
595  G4int iTkin, iTransfer, iTmin;
596  G4double t, position;
597  // G4double qMin = std::sqrt(tMin);
598 
599  fTableT = fBankT[0];
600  iTkin = 0;
601 
602  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
603  {
604  // if( qMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
605  if( tMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
606  }
607  iTmin = iTransfer-1;
608  if(iTmin < 0 ) iTmin = 0;
609 
610  position = (*(*fTableT)(iTkin))(iTmin)*G4UniformRand();
611 
612  for( iTmin = 0; iTransfer < fBinT-1; iTransfer++)
613  {
614  if( position > (*(*fTableT)(iTkin))(iTransfer) ) break;
615  }
616  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
617 
618  t = GetTransfer(iTkin, iTransfer, position);
619 
620  return t;
621 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4double GetTransfer(G4int iMomentum, G4int iTransfer, G4double position)
Definition: G4hhElastic.cc:629
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4hhElastic::SetAlphaP ( G4double  a)
inline

Definition at line 186 of file G4hhElastic.hh.

186 {fAlphaP = a;};
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void G4hhElastic::SetBq ( G4double  b)
inline

Definition at line 174 of file G4hhElastic.hh.

174 {fBq = b;};
tuple b
Definition: test.py:12

Here is the caller graph for this function:

void G4hhElastic::SetBQ ( G4double  b)
inline

Definition at line 175 of file G4hhElastic.hh.

175 {fBQ = b;};
tuple b
Definition: test.py:12

Here is the caller graph for this function:

void G4hhElastic::SetBqQ ( G4double  b)
inline

Definition at line 176 of file G4hhElastic.hh.

176 {fBqQ = b;};
tuple b
Definition: test.py:12
void G4hhElastic::SetCofF2 ( G4double  f)
inline

Definition at line 191 of file G4hhElastic.hh.

191 {fCofF2 = f;};

Here is the caller graph for this function:

void G4hhElastic::SetCofF3 ( G4double  f)
inline

Definition at line 192 of file G4hhElastic.hh.

192 {fCofF3 = f;};

Here is the caller graph for this function:

void G4hhElastic::SetEta ( G4double  E)
inline

Definition at line 190 of file G4hhElastic.hh.

190 {fEta = E;};

Here is the caller graph for this function:

void G4hhElastic::SetImCof ( G4double  a)
inline

Definition at line 187 of file G4hhElastic.hh.

187 {fImCof = a;};
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33

Here is the caller graph for this function:

void G4hhElastic::SetLambda ( G4double  L)
inline

Definition at line 189 of file G4hhElastic.hh.

189 {fLambda = L;};
static constexpr double L
Definition: G4SIunits.hh:124

Here is the caller graph for this function:

void G4hhElastic::SetParameters ( )
inline

Definition at line 273 of file G4hhElastic.hh.

274 {
275  // masses
276 
277  fMq = 0.36*CLHEP::GeV; // 0.441*GeV; // 0.36*GeV;
278  fMQ = 0.441*CLHEP::GeV;
279  fMff2 = 0.26*CLHEP::GeV*CLHEP::GeV; // 0.25*GeV*GeV; // 0.5*GeV*GeV;
280 
281  fAlpha = 1./3.;
282  fBeta = 1. - fAlpha;
283 
284  fGamma = 1./2.; // 1./3.; //
285  fDelta = 1. - fGamma; // 1./2.;
286 
287  // radii and exp cof
288 
289  fRA = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
290  fRq = 0.173*fRA; // 2.4/GeV;
291  fRQ = 0.316*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
292  fRB = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
293  fRg = 0.173*fRA; // 2.4/GeV;
294  fRG = 0.173*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
295 
296  fAlphaP = 0.15/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
297  fLambda = 0.25*fRA*fRA; // 0.25
298  fEta = 0.25*fRB*fRB; // 0.25
299  fImCof = 6.5;
300  fCofF2 = 1.;
301  fCofF3 = 1.;
302 
303  fBq = 0.02; // 0.21; // 1./3.;
304  fBQ = 1. + fBq - 2*std::sqrt(fBq); // 1 - fBq; // 2./3.;
305  fBqQ = std::sqrt(fBq*fBQ);
306 
307  fLambdaFF = 1.5/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
308  fSo = 1.*CLHEP::GeV*CLHEP::GeV;
309  fQcof = 0.009*CLHEP::GeV;
310  fExpSlope = 19.9/CLHEP::GeV/CLHEP::GeV;
311 }
static constexpr double GeV

Here is the caller graph for this function:

void G4hhElastic::SetParametersCMS ( G4double  plab)
inline

Definition at line 318 of file G4hhElastic.hh.

319 {
320  G4int i;
321  G4double trMass = 900.*CLHEP::MeV, Tkin;
322  G4double sl, sh, ds, rAl, rAh, drA, rBl, rBh, drB, bql, bqh, dbq, bQl, bQh, dbQ, cIl, cIh, dcI;
323 
324  Tkin = std::sqrt(fMassProj*fMassProj + plab*plab) - fMassProj;
325 
326  G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(fProjectile,
327  G4ParticleMomentum(0.,0.,1.),
328  Tkin);
329  fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, fTarget );
330 
331  delete theDynamicParticle;
332 
333  fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
334  fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
335 
336  G4double sCMS = std::sqrt(fSpp);
337 
338  if( fMassProj > trMass ) // p,n,pb on p
339  {
340  this->SetCofF2(1.);
341  this->SetCofF3(1.);
342  fGamma = 1./3.; // 1./3.; //
343  fDelta = 1. - fGamma; // 1./2.;
344 
345  if( sCMS <= theNuclNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
346  {
347  this->SetRA(theNuclNuclData[0][1]/CLHEP::GeV,0.173,0.316);
348  this->SetRB(theNuclNuclData[0][2]/CLHEP::GeV,0.173,0.316);
349 
350  this->SetBq(theNuclNuclData[0][3]);
351  this->SetBQ(theNuclNuclData[0][4]);
352  this->SetImCof(theNuclNuclData[0][5]);
353 
354  this->SetLambda(0.25*this->GetRA()*this->GetRA());
355  this->SetEta(0.25*this->GetRB()*this->GetRB());
356  }
357  else if( sCMS >= theNuclNuclData[17][0]*CLHEP::GeV ) // high edge, as s=7000 ???
358  {
359  this->SetRA(theNuclNuclData[17][1]/CLHEP::GeV,0.173,0.316);
360  this->SetRB(theNuclNuclData[17][2]/CLHEP::GeV,0.173,0.316);
361 
362  this->SetBq(theNuclNuclData[17][3]);
363  this->SetBQ(theNuclNuclData[17][4]);
364  this->SetImCof(theNuclNuclData[17][5]);
365 
366  this->SetLambda(0.25*this->GetRA()*this->GetRA());
367  this->SetEta(0.25*this->GetRB()*this->GetRB());
368  }
369  else // in approximation between array points
370  {
371  for( i = 0; i < 18; i++ ) if( sCMS <= theNuclNuclData[i][0]*CLHEP::GeV ) break;
372  if( i == 0 ) i++;
373  if( i == 18 ) i--;
374 
375  sl = theNuclNuclData[i-1][0]*CLHEP::GeV;
376  sh = theNuclNuclData[i][0]*CLHEP::GeV;
377  ds = (sCMS - sl)/(sh - sl);
378 
379  rAl = theNuclNuclData[i-1][1]/CLHEP::GeV;
380  rAh = theNuclNuclData[i][1]/CLHEP::GeV;
381  drA = rAh - rAl;
382 
383  rBl = theNuclNuclData[i-1][2]/CLHEP::GeV;
384  rBh = theNuclNuclData[i][2]/CLHEP::GeV;
385  drB = rBh - rBl;
386 
387  bql = theNuclNuclData[i-1][3];
388  bqh = theNuclNuclData[i][3];
389  dbq = bqh - bql;
390 
391  bQl = theNuclNuclData[i-1][4];
392  bQh = theNuclNuclData[i][4];
393  dbQ = bQh - bQl;
394 
395  cIl = theNuclNuclData[i-1][5];
396  cIh = theNuclNuclData[i][5];
397  dcI = cIh - cIl;
398 
399  this->SetRA(rAl+drA*ds,0.173,0.316);
400  this->SetRB(rBl+drB*ds,0.173,0.316);
401 
402  this->SetBq(bql+dbq*ds);
403  this->SetBQ(bQl+dbQ*ds);
404  this->SetImCof(cIl+dcI*ds);
405 
406  this->SetLambda(0.25*this->GetRA()*this->GetRA());
407  this->SetEta(0.25*this->GetRB()*this->GetRB());
408  }
409  }
410  else // pi, K
411  {
412  this->SetCofF2(1.);
413  this->SetCofF3(-1.);
414  fGamma = 1./2.; // 1./3.; //
415  fDelta = 1. - fGamma; // 1./2.;
416 
417  if( sCMS <= thePiKaNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
418  {
419  this->SetRA(thePiKaNuclData[0][1]/CLHEP::GeV,0.173,0.316);
420  this->SetRB(thePiKaNuclData[0][2]/CLHEP::GeV,0.173,0.173);
421 
422  this->SetBq(thePiKaNuclData[0][3]);
423  this->SetBQ(thePiKaNuclData[0][4]);
424  this->SetImCof(thePiKaNuclData[0][5]);
425 
426  this->SetLambda(0.25*this->GetRA()*this->GetRA());
427  this->SetEta(this->GetRB()*this->GetRB()/6.);
428  }
429  else if( sCMS >= thePiKaNuclData[7][0]*CLHEP::GeV ) // high edge, as s=7000 ???
430  {
431  this->SetRA(thePiKaNuclData[7][1]/CLHEP::GeV,0.173,0.316);
432  this->SetRB(thePiKaNuclData[7][2]/CLHEP::GeV,0.173,0.173);
433 
434  this->SetBq(thePiKaNuclData[7][3]);
435  this->SetBQ(thePiKaNuclData[7][4]);
436  this->SetImCof(thePiKaNuclData[7][5]);
437 
438  this->SetLambda(0.25*this->GetRA()*this->GetRA());
439  this->SetEta(this->GetRB()*this->GetRB()/6.);
440  }
441  else // in approximation between array points
442  {
443  for( i = 0; i < 8; i++ ) if( sCMS <= thePiKaNuclData[i][0]*CLHEP::GeV ) break;
444  if( i == 0 ) i++;
445  if( i == 8 ) i--;
446 
447  sl = thePiKaNuclData[i-1][0]*CLHEP::GeV;
448  sh = thePiKaNuclData[i][0]*CLHEP::GeV;
449  ds = (sCMS - sl)/(sh - sl);
450 
451  rAl = thePiKaNuclData[i-1][1]/CLHEP::GeV;
452  rAh = thePiKaNuclData[i][1]/CLHEP::GeV;
453  drA = rAh - rAl;
454 
455  rBl = thePiKaNuclData[i-1][2]/CLHEP::GeV;
456  rBh = thePiKaNuclData[i][2]/CLHEP::GeV;
457  drB = rBh - rBl;
458 
459  bql = thePiKaNuclData[i-1][3];
460  bqh = thePiKaNuclData[i][3];
461  dbq = bqh - bql;
462 
463  bQl = thePiKaNuclData[i-1][4];
464  bQh = thePiKaNuclData[i][4];
465  dbQ = bQh - bQl;
466 
467  cIl = thePiKaNuclData[i-1][5];
468  cIh = thePiKaNuclData[i][5];
469  dcI = cIh - cIl;
470 
471  this->SetRA(rAl+drA*ds,0.173,0.316);
472  this->SetRB(rBl+drB*ds,0.173,0.173);
473 
474  this->SetBq(bql+dbq*ds);
475  this->SetBQ(bQl+dbQ*ds);
476  this->SetImCof(cIl+dcI*ds);
477 
478  this->SetLambda(0.25*this->GetRA()*this->GetRA());
479  this->SetEta(this->GetRB()*this->GetRB()/6.);
480  }
481  }
482  return;
483 }
void SetRA(G4double rn, G4double pq, G4double pQ)
Definition: G4hhElastic.hh:489
void SetBQ(G4double b)
Definition: G4hhElastic.hh:175
void SetImCof(G4double a)
Definition: G4hhElastic.hh:187
int G4int
Definition: G4Types.hh:78
G4double GetHadronNucleonXscNS(const G4DynamicParticle *, const G4ParticleDefinition *)
void SetCofF3(G4double f)
Definition: G4hhElastic.hh:192
void SetEta(G4double E)
Definition: G4hhElastic.hh:190
static constexpr double MeV
void SetCofF2(G4double f)
Definition: G4hhElastic.hh:191
void SetBq(G4double b)
Definition: G4hhElastic.hh:174
void SetRB(G4double rn, G4double pq, G4double pQ)
Definition: G4hhElastic.hh:500
static constexpr double GeV
double G4double
Definition: G4Types.hh:76
void SetLambda(G4double L)
Definition: G4hhElastic.hh:189
G4ThreeVector G4ParticleMomentum
G4double GetRA()
Definition: G4hhElastic.hh:196
G4double GetRB()
Definition: G4hhElastic.hh:200

Here is the call graph for this function:

Here is the caller graph for this function:

void G4hhElastic::SetRA ( G4double  rn,
G4double  pq,
G4double  pQ 
)
inline

Definition at line 489 of file G4hhElastic.hh.

490 {
491  fRA = rA;
492  fRq = fRA*pq;
493  fRQ = fRA*pQ;
494 }

Here is the caller graph for this function:

void G4hhElastic::SetRB ( G4double  rn,
G4double  pq,
G4double  pQ 
)
inline

Definition at line 500 of file G4hhElastic.hh.

501 {
502  fRB = rB;
503  fRg = fRB*pg;
504  fRG = fRB*pG;
505 }
tuple pg
Definition: demo.py:37

Here is the caller graph for this function:

void G4hhElastic::SetSigmaTot ( G4double  stot)
inline

Definition at line 166 of file G4hhElastic.hh.

166 {fSigmaTot = stot;};
void G4hhElastic::SetSpp ( G4double  spp)
inline

Definition at line 167 of file G4hhElastic.hh.

167 {fSpp = spp;};

The documentation for this class was generated from the following files: