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 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	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 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
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.

  : G4HadronElastic("HadrHadrElastic")
{
  SetMinEnergy( 1.*GeV );
  SetMaxEnergy( 10000.*TeV );
  verboseLevel = 0;
  lowEnergyRecoilLimit = 100.*keV;  
  lowEnergyLimitQ  = 0.0*GeV;  
  lowEnergyLimitHE = 0.0*GeV;  
  lowestEnergyLimit= 0.0*keV;  
  plabLowLimit     = 20.0*MeV;

  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
  fInTkin=0;
  theProton   = G4Proton::Proton();
  theNeutron  = G4Neutron::Neutron();
  thePionPlus = G4PionPlus::PionPlus();
  thePionMinus= G4PionMinus::PionMinus();

  fTarget  = G4Proton::Proton();
  fProjectile  = 0;
  fHadrNuclXsc = new G4HadronNucleonXsc();

  fEnergyBin = 200;
  fBinT  = 514; // 514; // 500; // 200;

  fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );

  fTableT = 0;
  fOldTkin = 0.;
  SetParameters();

  Initialise();
}

Here is the call graph for this function:

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

Definition at line 114 of file G4hhElastic.cc.

  : G4HadronElastic("HadrHadrElastic")
{
  SetMinEnergy( 1.*GeV );
  SetMaxEnergy( 10000.*TeV );
  verboseLevel         = 0;
  lowEnergyRecoilLimit = 100.*keV;  
  lowEnergyLimitQ      = 0.0*GeV;  
  lowEnergyLimitHE     = 0.0*GeV;  
  lowestEnergyLimit    = 0.0*keV;  
  plabLowLimit         = 20.0*MeV;

  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
  fInTkin=0;
  theProton   = G4Proton::Proton();
  theNeutron  = G4Neutron::Neutron();
  thePionPlus = G4PionPlus::PionPlus();
  thePionMinus= G4PionMinus::PionMinus();

  fTarget      = target;
  fProjectile  = projectile;
  fMassTarg   = fTarget->GetPDGMass();
  fMassProj   = fProjectile->GetPDGMass();
  fMassSum2   = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
  fMassDif2   = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
  fHadrNuclXsc = new G4HadronNucleonXsc();

  fEnergyBin = 200;
  fBinT      = 514; // 200;

  fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
  fTableT       = 0;
  fOldTkin = 0.;


  SetParameters();
  SetParametersCMS( plab);
}

Here is the call graph for this function:

G4hhElastic::G4hhElastic	(	G4ParticleDefinition *	target,
		G4ParticleDefinition *	projectile
	)

Definition at line 159 of file G4hhElastic.cc.

  : G4HadronElastic("HadrHadrElastic")
{
  SetMinEnergy( 1.*GeV );
  SetMaxEnergy( 10000.*TeV );
  verboseLevel = 0;
  lowEnergyRecoilLimit = 100.*keV;  
  lowEnergyLimitQ  = 0.0*GeV;  
  lowEnergyLimitHE = 0.0*GeV;  
  lowestEnergyLimit= 0.0*keV;  
  plabLowLimit     = 20.0*MeV;

  fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
  fInTkin=0;

  fTarget = target; // later vmg
  fProjectile = projectile;
  theProton   = G4Proton::Proton();
  theNeutron  = G4Neutron::Neutron();
  thePionPlus = G4PionPlus::PionPlus();
  thePionMinus= G4PionMinus::PionMinus();

  fTarget  = G4Proton::Proton(); // later vmg
  fMassTarg   = fTarget->GetPDGMass();
  fMassProj   = fProjectile->GetPDGMass();
  fMassSum2   = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
  fMassDif2   = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
  fHadrNuclXsc = new G4HadronNucleonXsc();

  fEnergyBin = 200;
  fBinT  = 514; // 514; // 500; // 200;

  fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );

  fTableT = 0;
  fOldTkin = 0.;

  SetParameters();
}

Here is the call graph for this function:

G4hhElastic::~G4hhElastic ( )

virtual

Definition at line 205 of file G4hhElastic.cc.

{
  if ( fEnergyVector ) {
    delete fEnergyVector;
    fEnergyVector = 0;
  }

  for ( std::vector<G4PhysicsTable*>::iterator it = fBankT.begin();
        it != fBankT.end(); ++it ) {
    if ( (*it) ) (*it)->clearAndDestroy();
    delete *it;
    *it = 0;
  }
  fTableT = 0;
  if(fHadrNuclXsc) delete fHadrNuclXsc;
}

Member Function Documentation

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

Definition at line 254 of file G4hhElastic.cc.

{
  G4int iTkin, jTransfer;
  G4double plab, Tkin, tMax;
  G4double t1, t2, dt, delta = 0., sum = 0.;

  fTarget     = target;
  fProjectile = projectile;
  fMassTarg   = fTarget->GetPDGMass();
  fMassProj   = fProjectile->GetPDGMass();
  fMassSum2   = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
  fMassDif2   = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);

  G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
  // G4HadronNucleonXsc*          hnXsc = new G4HadronNucleonXsc();
  fTableT = new G4PhysicsTable(fEnergyBin);

  for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
  {
    Tkin  = fEnergyVector->GetLowEdgeEnergy(iTkin);
    plab  = std::sqrt( Tkin*( Tkin + 2*fMassProj ) );
    // G4DynamicParticle*  theDynamicParticle = new G4DynamicParticle(projectile,
    //                                           G4ParticleMomentum(0.,0.,1.),
    //                                           Tkin);
    // fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, target );

    SetParametersCMS( plab );

    tMax  = 4.*fPcms*fPcms;
    if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???

    G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
    sum = 0.;
    dt  = tMax/fBinT;

    // for(j = 1; j < fBinT; j++)

    for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer--)
    {
      t1 = dt*(jTransfer-1);
      t2 = t1 + dt;

      if( fMassProj > 900.*MeV ) // pp, pn
      {
        delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
        // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
      }
      else // pi+-p, K+-p
      {
        delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
        // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
      }
      sum += delta;
      vectorT->PutValue( jTransfer-1, t1, sum ); // t2
    }
    // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
    fTableT->insertAt( iTkin, vectorT );
    // delete theDynamicParticle;
  }
  // delete hnXsc;

  return;
}

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.

{
  G4int jTransfer;
  G4double tMax; // , sQq, sQG;
  G4double t1, t2, dt, delta = 0., sum = 0. ; // , threshold;

  fTarget     = target;
  fProjectile = projectile;
  fMassTarg =  fTarget->GetPDGMass();
  fMassProj =  fProjectile->GetPDGMass();
  fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
  fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
  fSpp  = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
  fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);

  G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fMassProj = "<<fMassProj<<" MeV"<<G4endl;
  tMax = 4.*fPcms*fPcms;
  if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???

 
  G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
  fTableT = new G4PhysicsTable(1);
  G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);

  sum = 0.;
  dt = tMax/G4double(fBinT);
  G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; fPcms = "<<fPcms/GeV
    <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt = "<<dt/GeV/GeV<<" GeV2"<<G4endl;

  // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fRB*GeV<<G4endl;

    // for(jTransfer = 1; jTransfer < fBinT; jTransfer++)
  for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer-- )
  {
    t1 = dt*(jTransfer-1);
    t2 = t1 + dt;

    if( fMassProj > 900.*MeV ) // pp, pn
    {
      delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
      // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, tMax);
    }
    else // pi+-p, K+-p
    {
      delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
      // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, tMax);
      // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
    }
    sum += delta;
    // G4cout<<delta<<"\t"<<sum<<"\t"<<threshold<<G4endl;   
  
    // sQq = GetdsdtF123(q1);
    // sQG = GetdsdtF123qQgG(q1);
    // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/millibarn<<"\t"<<sQq*GeV*GeV/millibarn<<G4endl;
    // G4cout<<"sum = "<<sum<<", "; 
    
    vectorT->PutValue( jTransfer-1, t1, sum );  // t2
  }
  // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
  fTableT->insertAt( 0, vectorT );
  fBankT.push_back( fTableT );  // 0

  // for(jTransfer = 0; jTransfer < fBinT-1; jTransfer++) 
  //   G4cout<<(*(*fTableT)(0))(jTransfer)/sum<<"\t\t"<<G4Pow::GetInstance()->powN(2.,-jTransfer)<<G4endl;
  
  return;
}

Here is the call graph for this function:

void G4hhElastic::CalculateBQ ( G4double  b )

inline

Definition at line 959 of file G4hhElastic.hh.

{
  fBq = b1;
  G4double s1 = GetCofS1();
  G4double s2 = GetCofS2();
  G4double s3 = GetCofS3();
  G4double sqrtBq = std::sqrt(fBq);

  // cofs of the fBQ 3rd equation

  G4double a = s3*sqrtBq;
  G4double b = s1*fBq - 1.;
  G4double c = (s2*fBq - 2.)*sqrtBq;
  G4double d = 1. - fBq;

  // cofs of the incomplete 3rd equation

  G4double p  = c/a;
           p -= b*b/a/a/3.;
  G4double q  = d/a;
           q -= b*c/a/a/3.;
           q += 2*b*b*b/a/a/a/27.;

  // cofs for the incomplete colutions

  G4double D  = p*p*p/3./3./3.;
           D += q*q/2./2.;
  G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
  G4complex A  = std::pow(A1,1./3.);

  G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
  G4complex B  = std::pow(B1,1./3.);

  // roots of the incomplete 3rd equation

  G4complex y1 = A + B;
  G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
  G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
 
  G4complex x1 = y1 - b/a/3.;
  G4complex x2 = y2 - b/a/3.;
  G4complex x3 = y3 - b/a/3.;

  G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = "<<real(x2)<<"; re_x3 = "<<real(x3)<<G4endl;
  G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = "<<imag(x2)<<"; im_x3 = "<<imag(x3)<<G4endl;

  G4double r1 = real(x1)*real(x1);
  G4double r2 = real(x2)*real(x2);
  G4double r3 = real(x3)*real(x3);

  if( r1 <= 1. && r1 >= 0. )      fBQ = r1;
  else if( r2 <= 1. && r2 >= 0. ) fBQ = r2;
  else if( r3 <= 1. && r3 >= 0. )      fBQ = r3;
  else fBQ = 1.;
  // fBQ = real(x3)*real(x3);
  G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
  fOptRatio  = fBq + fBQ + 2.*sqrtBqBQ - 1.;
  fOptRatio /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
  G4cout<<"F123, fOptRatio = "<<fOptRatio<<G4endl;
  
  return ;
}

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ12 ( G4double  b )

inline

Definition at line 788 of file G4hhElastic.hh.

{
  fBq = b1;

  G4complex z1324 = G4complex(1./8./CLHEP::pi,0.);
            z1324 /= Phi13() + Phi24() + fLambda + fEta;
        G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z1423 = G4complex(1./8./CLHEP::pi,0.);
            z1423 /= Phi14() + Phi23() + fLambda + fEta;
        G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);

  fBQ  = 1. - 2.*fBq;
  fBQ /= 2. - fSigmaTot*fBq*(c1324+1423); 

  G4double ratio  = 2.*(fBq + fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423); 
  G4cout<<"ratio = "<<ratio<<G4endl;

  return ;
}

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ123 ( G4double  b )

inline

Definition at line 814 of file G4hhElastic.hh.

{
  fBq = b1;

  G4complex z1324 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
            z1324 /= Phi13() + Phi24() + fLambda + fEta;
        G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z1423 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
            z1423 /= Phi14() + Phi23() + fLambda + fEta;
        G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z1314 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
            z1314 /= Phi13() + Phi14() + fEta;
        G4double c1314 = real(z1314)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z2324 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
            z2324 /= Phi23() + Phi24() + fEta;
        G4double c2324 = real(z2324)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z1323 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
            z1323 /= Phi13() + Phi23() + fLambda;
        G4double c1323 = real(z1323)/(CLHEP::hbarc*CLHEP::hbarc);

  G4complex z1424 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
            z1424 /= Phi14() + Phi24() + fLambda;
        G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);

  G4double A = fSigmaTot*c2324;
  G4double B = fSigmaTot*fBq*(c1324 + c1423 + c1323 + c1424) - 2.;
  G4double C = 1. + fSigmaTot*fBq*fBq*c1314 - 2*fBq;
  G4cout<<"A = "<<A<<"; B = "<<B<<"; C = "<<C<<G4endl;
  G4cout<<"determinant = "<<B*B-4.*A*C<<G4endl;

  G4double x1 = ( -B - std::sqrt(B*B-4.*A*C) )/2./A;
  G4double x2 = ( -B + std::sqrt(B*B-4.*A*C) )/2./A;
  G4cout<<"x1 = "<<x1<<"; x2 = "<<x2<<G4endl;

  if( B*B-4.*A*C < 1.e-6 ) fBQ = std::abs(-B/2./A);
  else if ( B < 0.)        fBQ = std::abs( ( -B - std::sqrt(B*B-4.*A*C) )/2./A);
  else                     fBQ  = std::abs( ( -B + std::sqrt(B*B-4.*A*C) )/2./A);
 
  fOptRatio  = 2*(fBq+fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423 + c1323 + c1424);
  fOptRatio -= fSigmaTot*fBq*fBq*c1314 + fSigmaTot*c2324*fBQ*fBQ;
  G4cout<<"BqQ123, fOptRatio = "<<fOptRatio<<G4endl;

  return ;
}

Here is the call graph for this function:

void G4hhElastic::CalculateBqQ13 ( G4double  b )

inline

Definition at line 765 of file G4hhElastic.hh.

{
  fBQ = b2;

  G4complex z1424 = G4complex(1./8./CLHEP::pi,0.);
            z1424 /= Phi14() + Phi24() + fAlpha;
        G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);

  fBqQ  = 1. - fBQ;
  fBQ /= 1. - fSigmaTot*fBQ*c1424; 

  G4cout<<"fSigmaTot*fBQ*c1424 = "<<fSigmaTot*fBQ*c1424<<G4endl;

  G4double ratio  = fBqQ + fBQ - fSigmaTot*fBqQ*fBQ*c1424; 
  G4cout<<"ratio = "<<ratio<<G4endl;

  return ;
}

Here is the call graph for this function:

G4complex G4hhElastic::GetAqq ( )

inline

Definition at line 869 of file G4hhElastic.hh.

{
  G4double re, im;

  re = fRq*fRq/8. + fAlphaP*G4Log(fSpp/fSo) + 8.*fLambda/9.;
  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
  return G4complex(re,im);
}

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.

{
  G4double re, im;

  re = fRQ*fRQ/8. + fAlphaP*G4Log(fSpp/fSo) + 2.*fLambda/9.;
  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
  return G4complex(re,im);
}

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.

{
  G4complex z = 0.5*( GetAqq() + GetAQQ() );
  return z;
}

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.

{ return fBq;};

G4double G4hhElastic::GetBQ ( )

inline

Definition at line 172 of file G4hhElastic.hh.

{ return fBQ;};

G4double G4hhElastic::GetBqQ ( )

inline

Definition at line 173 of file G4hhElastic.hh.

{ return fBqQ;};

G4double G4hhElastic::GetCofF2 ( )

inline

Definition at line 193 of file G4hhElastic.hh.

{return fCofF2;};

G4double G4hhElastic::GetCofF3 ( )

inline

Definition at line 194 of file G4hhElastic.hh.

{return fCofF3;};

G4double G4hhElastic::GetCofS1 ( )

inline

Definition at line 905 of file G4hhElastic.hh.

{
  G4complex z = 1./( GetAqQ() + 4.*fLambda/9. );
  G4double result = real(z);
  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
  result *= fSigmaTot*fCofF2;
  return result;
}

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.

{
  G4complex z = 1./( GetAqq() + GetAqQ() - 4.*fLambda/9. );
  G4double result = real(z);
  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
  result *= fSigmaTot*fCofF3;
  return result;
}

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.

{
  G4complex z = 1./( GetAQQ() + GetAqQ() + 2.*fLambda/9. );
  G4double result = real(z);
  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
  result *= fSigmaTot*fCofF3;
  return result;
}

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.

{
  fSpp = spp;
  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4complex F1 = GetF1(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
  return dsdt;
}

Here is the call graph for this function:

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

inline

Definition at line 1087 of file G4hhElastic.hh.

{
  fSpp = spp;
  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4complex F1 = GetF1(t) - GetF2(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
  return dsdt;
}

Here is the call graph for this function:

G4double G4hhElastic::GetdsdtF123 ( G4double  q )

inline

Definition at line 1130 of file G4hhElastic.hh.

{
  G4double  p   = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4complex F1  = GetF1(t);
            F1 -= fCofF2*GetF2(t);
            F1 -= fCofF3*GetF3(t);
  G4double dsdt = CLHEP::pi/p/p;
  dsdt         *= real(F1)*real(F1) + imag(F1)*imag(F1);  
  return dsdt;
}

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.

{
  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp );

  G4complex F123  = GetF1qQgG(t); // - fCofF2*GetF2qQgG(t) - fCofF3*GetF3qQgG(t);
            F123 -= fCofF2*GetF2qQgG(t);
            F123 -= fCofF3*GetF3qQgG(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F123)*real(F123) + imag(F123)*imag(F123);  
  return dsdt;
}

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.

{
  fSpp = spp;
  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);

  G4complex F12 = GetF1qQgG(t) - GetF2qQgG(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F12)*real(F12) + imag(F12)*imag(F12);  
  return dsdt;
}

Here is the call graph for this function:

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

inline

Definition at line 586 of file G4hhElastic.hh.

{
  fSpp = spp;
  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
  G4double k = p/CLHEP::hbarc;

  G4complex exp14 = fBqQ*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);

  G4complex F1  = exp14 + exp24;

  F1 *=  0.25*k*fSigmaTot/CLHEP::pi;
  F1 *= G4complex(0.,1.);

  // 1424

  G4complex z1424  = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
            z1424 /= Phi14() + Phi24() + fLambda;
            z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;

 
  G4complex exp1424  = std::exp(-z1424*t);
            exp1424 /= Phi14() + Phi24() + fLambda;

  G4complex F3  = fBqQ*fBQ*exp1424;


  F3 *=  0.25*k/CLHEP::pi;
  F3 *= G4complex(0.,1.);
  F3 *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);

  G4complex F13  = F1 - F3;

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F13)*real(F13) + imag(F13)*imag(F13);  

  return dsdt;
}

Here is the call graph for this function:

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

inline

Definition at line 629 of file G4hhElastic.hh.

{
  fSpp = spp;
  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);

  G4complex F1 = GetF1qQgG(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
  return dsdt;
}

Here is the call graph for this function:

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

inline

Definition at line 1145 of file G4hhElastic.hh.

{
  fSpp = spp;
  G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);

  // qQ-ds/dt

  G4complex F1 = GetF1(t) - fCofF2*GetF2(t) - fCofF3*GetF3(t);

  G4double dsdt = CLHEP::pi/p/p;
  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1); 

  // exponent ds/dt

  G4complex F10 = GetF1(0.) - fCofF2*GetF2(0.) - fCofF3*GetF3(0.);

  fRhoReIm = real(F10)/imag(F10);

  G4double dsdt0 = CLHEP::pi/p/p;
  dsdt0 *= real(F10)*real(F10) + imag(F10)*imag(F10); 

  dsdt0 *= G4Exp(-fExpSlope*t);

  G4double ratio = dsdt/dsdt0;
 
  return ratio;
}

Here is the call graph for this function:

G4complex G4hhElastic::GetF1 ( G4double  qp )

inline

Definition at line 1026 of file G4hhElastic.hh.

{
  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4double  k    = p/CLHEP::hbarc;
  G4complex exp1 = fBq*std::exp(-GetAqq()*t);
  G4complex exp2 = fBQ*std::exp(-GetAQQ()*t);
  G4complex exp3 = 2.*std::sqrt(fBq*fBQ)*std::exp(-GetAqQ()*t);

  G4complex res  = exp1 + exp2 + exp3;
  res *=  0.25*k*fSigmaTot/CLHEP::pi;
  res *= G4complex(0.,1.);
  return res;
}

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.

{
  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
  G4double k = p/CLHEP::hbarc;

  G4complex exp13 = fBq*std::exp(-(Phi13() + fBeta*fBeta*fLambda + fDelta*fDelta*fEta)*t);
  G4complex exp14 = fBq*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
  G4complex exp23 = fBQ*std::exp(-(Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta)*t);
  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);

  G4complex res  = exp13 + exp14 + exp23 + exp24;

  res *=  0.25*k*fSigmaTot/CLHEP::pi;
  res *= G4complex(0.,1.);

  return res;
}

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.

{
  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4double  k   = p/CLHEP::hbarc;
  G4complex z1  = GetAqq()*GetAQQ() - 16.*fLambda*fLambda/81.;
            z1 /= 2.*(GetAqQ() + 4.*fLambda/9.);
  G4complex exp1 = std::exp(-z1*t);

  G4complex z2 = 0.5*( GetAqQ() - 4.*fLambda/9.);

  G4complex exp2 = std::exp(-z2*t);

  G4complex res  = exp1 + exp2;

  G4complex z3 = GetAqQ() + 4.*fLambda/9.;

  res *=  0.25*k/CLHEP::pi;
  res *= G4complex(0.,1.);
  res /= z3;
  res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);

  return res;
}

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.

{
  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
  G4double k = p/CLHEP::hbarc;

  G4complex z1324  = -(Phi24() + fAlpha*fLambda + fGamma*fEta)*(Phi24() + fAlpha*fLambda + fGamma*fEta);
            z1324 /= Phi13() + Phi24() + fLambda + fEta;
            z1324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;

  G4complex exp1324  = std::exp(-z1324*t);
            exp1324 /= Phi13() + Phi24() + fLambda + fEta;

  G4complex z1423 = -(Phi23() + fAlpha*fLambda + fDelta*fEta)*(Phi24() + fAlpha*fLambda + fDelta*fEta);;
            z1423 /= Phi14() + Phi23() + fLambda + fEta;
            z1423 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;

  G4complex exp1423 = std::exp(-z1423*t);
            exp1423 /= Phi14() + Phi23() + fLambda + fEta; 

  G4complex res  = exp1324 + exp1423;


  res *=  0.25*k/CLHEP::pi;
  res *= G4complex(0.,1.);
  res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc); // or 4. ???

  return res;
}

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.

{
  G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
  G4double  k   = p/CLHEP::hbarc;
  G4complex z1  = GetAqq()*GetAqQ() - 4.*fLambda*fLambda/81.;
            z1 /= GetAqq() + GetAqQ() - 4.*fLambda/9.;

  G4complex exp1 = std::exp(-z1*t)*fBq/(GetAqq() + GetAqQ() - 4.*fLambda/9.);

  G4complex z2 = GetAqQ()*GetAQQ() - 1.*fLambda*fLambda/81.;
            z2 /= GetAQQ() + GetAqQ() + 2.*fLambda/9.;

  G4complex exp2 = std::exp(-z2*t)*fBQ/(GetAQQ() + GetAqQ() + 2.*fLambda/9.);

  G4complex res  = exp1 + exp2;


  res *=  0.25*k/CLHEP::pi;
  res *= G4complex(0.,1.);
  res *= std::sqrt(fBq*fBQ)*fSigmaTot*fSigmaTot/(4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);

  return res;
}

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.

{
  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
  G4double k = p/CLHEP::hbarc;

            // 1314

  G4complex z1314  = -(Phi14() + fGamma*fEta)*(Phi14() + fGamma*fEta);
            z1314 /= Phi13() + Phi14() + fEta;
            z1314 += Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta;

  G4complex exp1314  = std::exp(-z1314*t);
            exp1314 /= Phi13() + Phi14() + fEta;

        // 2324

  G4complex z2324 = -(Phi24() + fGamma*fEta)*(Phi24() + fGamma*fEta);;
            z2324 /= Phi24() + Phi23() + fEta;
            z2324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;

  G4complex exp2324 = std::exp(-z2324*t);
            exp2324 /= Phi24() + Phi23() + fEta; 

        // 1323

  G4complex z1323  = -(Phi23() + fAlpha*fLambda)*(Phi23() + fAlpha*fLambda);
            z1323 /= Phi13() + Phi23() + fLambda;
            z1323 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;

  G4complex exp1323  = std::exp(-z1323*t);
            exp1323 /= Phi13() + Phi23() + fLambda;

        // 1424

  G4complex z1424  = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
            z1424 /= Phi14() + Phi24() + fLambda;
            z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;

  G4complex exp1424  = std::exp(-z1424*t);
            exp1424 /= Phi14() + Phi24() + fLambda;

        G4complex res  = fBq*fBq*exp1314 + fBQ*fBQ*exp2324 + fBq*fBQ*exp1323 + fBq*fBQ*exp1424;

  res *=  0.25*k/CLHEP::pi;
  res *= G4complex(0.,1.);
  res *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);

  return res;
}

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.

{return fImCof;};

G4double G4hhElastic::GetOpticalRatio ( )

inline

Definition at line 944 of file G4hhElastic.hh.

{
  return fOptRatio;
  // G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
  // G4double result = fBq + fBQ + 2.*sqrtBqBQ - 1.;
  // result /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
  // return result;
}

G4double G4hhElastic::GetRA ( )

inline

Definition at line 196 of file G4hhElastic.hh.

{ return fRA;};

Here is the caller graph for this function:

G4double G4hhElastic::GetRB ( )

inline

Definition at line 200 of file G4hhElastic.hh.

{ return fRB;};

Here is the caller graph for this function:

G4double G4hhElastic::GetRg ( )

inline

Definition at line 201 of file G4hhElastic.hh.

{ return fRg;};

G4double G4hhElastic::GetRG ( )

inline

Definition at line 202 of file G4hhElastic.hh.

{ return fRG;};

G4double G4hhElastic::GetRhoReIm ( )

inline

Definition at line 177 of file G4hhElastic.hh.

{ return fRhoReIm;};

G4double G4hhElastic::GetRq ( )

inline

Definition at line 197 of file G4hhElastic.hh.

{ return fRq;};

G4double G4hhElastic::GetRQ ( )

inline

Definition at line 198 of file G4hhElastic.hh.

{ return fRQ;};

G4double G4hhElastic::GetSpp ( )

inline

Definition at line 168 of file G4hhElastic.hh.

{return fSpp;};

G4double G4hhElastic::GetTransfer	(	G4int	iMomentum,
		G4int	iTransfer,
		G4double	position
	)

Definition at line 629 of file G4hhElastic.cc.

{
  G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;

  if( iTransfer == 0 )
  {
    randTransfer = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
    // iTransfer++;
  }
  else
  {
    if ( iTransfer >= G4int((*fTableT)(iTkin)->GetVectorLength()) )
    {
      iTransfer = (*fTableT)(iTkin)->GetVectorLength() - 1;
    }
    y1 = (*(*fTableT)(iTkin))(iTransfer-1);
    y2 = (*(*fTableT)(iTkin))(iTransfer);

    x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
    x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);

    delta = y2 - y1;
    mean  = y2 + y1;

    if ( x1 == x2 ) randTransfer = x2;
    else
    {
      // if ( y1 == y2 ) 
      if ( delta < epsilon*mean ) 
           randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
      else randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
    }
  }
  return randTransfer;
}

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.

{
  // pp,pn

  fProjectile = G4Proton::Proton();
  BuildTableT(fTarget, fProjectile);
  fBankT.push_back(fTableT); // 0

  // pi+-p
  
  fProjectile = G4PionPlus::PionPlus();
  BuildTableT(fTarget, fProjectile);
  fBankT.push_back(fTableT);  // 1
  //K+-p
  fProjectile = G4KaonPlus::KaonPlus();
  BuildTableT(fTarget, fProjectile);
  fBankT.push_back(fTableT);  // 2
  
}

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.

{
  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
        projectile.GetDefinition() == G4Neutron::Neutron() ||
        projectile.GetDefinition() == G4PionPlus::PionPlus() ||
        projectile.GetDefinition() == G4PionMinus::PionMinus() ||
        projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
        projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&

        nucleus.GetZ_asInt() < 2 ) return true;
  else                              return false;
}

Here is the call graph for this function:

G4complex G4hhElastic::Phi13 ( )

inline

Definition at line 522 of file G4hhElastic.hh.

{
  G4double re = (fRq*fRq + fRg*fRg)/16.;
  G4complex result(re,0.);
  result += Pomeron();
  return result;
}

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.

{
  G4double re = (fRq*fRq + fRG*fRG)/16.;
  G4complex result(re,0.);
  result += Pomeron();
  return result;
}

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.

{
  G4double re = (fRQ*fRQ + fRg*fRg)/16.;
  G4complex result(re,0.);
  result += Pomeron();
  return result;
}

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.

{
  G4double re = (fRQ*fRQ + fRG*fRG)/16.;
  G4complex result(re,0.);
  result += Pomeron();
  return result;
}

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.

{
  G4double re, im;

  re = fAlphaP*G4Log(fSpp/fSo);
  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
  return G4complex(re,im);
} 

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.

{
  G4int iTkin, iTransfer;
  G4double t,  position, m1 = aParticle->GetPDGMass();
  G4double Tkin = std::sqrt(m1*m1+p*p) - m1;

  if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
  {
    fTableT = fBankT[0];
  }
  if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
  {
    fTableT = fBankT[1];
  }
  if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
  {
    fTableT = fBankT[2];
  }
  G4double delta    = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
  G4double deltaMax = 1.e-2;

  if ( delta < deltaMax ) iTkin = fInTkin; 
  else
  {  
    for( iTkin = 0; iTkin < fEnergyBin; iTkin++ )
    {
      if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
    }
  }
  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1;   // Tkin is more then theMaxEnergy
  if ( iTkin < 0 )           iTkin = 0; // against negative index, Tkin < theMinEnergy

  fOldTkin = Tkin;
  fInTkin = iTkin;

  if (iTkin == fEnergyBin -1 || iTkin == 0 )   // the table edges
  {
    position = (*(*fTableT)(iTkin))(0)*G4UniformRand();

    for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
    {
      if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
    }
    if (iTransfer >= fBinT-1) iTransfer = fBinT-2;

    t = GetTransfer(iTkin, iTransfer, position);


  }
  else  // Tkin inside between energy table edges
  {
    G4double rand = G4UniformRand();
    position = (*(*fTableT)(iTkin))(0)*rand;

    //
    // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBinT-2);
    G4int sTransfer = 0, fTransfer =  fBinT - 2, dTransfer = fTransfer - sTransfer;
    G4double y2;

    for( iTransfer = 0; iTransfer < fBinT - 1; iTransfer++ )
    {
      // dTransfer %= 2;
      dTransfer /= 2;
      // dTransfer *= 0.5;
      y2 = (*(*fTableT)(iTkin))( sTransfer + dTransfer );

      if( y2 > position ) sTransfer += dTransfer;

      // if( dTransfer <= 1 ) break;
      if( dTransfer < 1 ) break;
    }
    t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sTransfer); // +(-0.5+rand)*(*fTableT)(iTkin)->GetLowEdgeEnergy(3);
  }
  return t;
}

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.

{
  G4int iTkin, iTransfer;
  G4double t, t2, position, m1 = aParticle->GetPDGMass();
  G4double Tkin = std::sqrt(m1*m1+p*p) - m1;

  if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
  {
    fTableT = fBankT[0];
  }
  if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
  {
    fTableT = fBankT[1];
  }
  if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
  {
    fTableT = fBankT[2];
  }

  G4double delta    = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
  G4double deltaMax = 1.e-2;

  if ( delta < deltaMax ) iTkin = fInTkin; 
  else
  {  
    for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
    {
      if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
    }
  }
  if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1;   // Tkin is more then theMaxEnergy
  if ( iTkin < 0 )           iTkin = 0; // against negative index, Tkin < theMinEnergy

  fOldTkin = Tkin;
  fInTkin = iTkin;

  if (iTkin == fEnergyBin -1 || iTkin == 0 )   // the table edges
  {
    position = (*(*fTableT)(iTkin))(0)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
    {
      if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
    }
    if (iTransfer >= fBinT-1) iTransfer = fBinT-2;

    // G4cout<<"iTransfer = "<<iTransfer<<G4endl;

    t = GetTransfer(iTkin, iTransfer, position);

    // G4cout<<"t = "<<t<<G4endl;
  }
  else  // Tkin inside between energy table edges
  {
    // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
    position = (*(*fTableT)(iTkin))(0)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
    {
      // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
      if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
    }
    if (iTransfer >= fBinT-1) iTransfer = fBinT-2;

    // G4cout<<"iTransfer = "<<iTransfer<<G4endl;

    t2  = GetTransfer(iTkin, iTransfer, position);
    return t2;
    /*
    G4double t1, E1, E2, W, W1, W2;
    // G4cout<<"t2 = "<<t2<<G4endl;

    E2 = fEnergyVector->GetLowEdgeEnergy(iTkin);

    // G4cout<<"E2 = "<<E2<<G4endl;
    
    iTkin--;
    
    // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();

    // G4cout<<"position = "<<position<<G4endl;

    for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
    {
      // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
      if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
    }
    if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
    
    t1  = GetTransfer(iTkin, iTransfer, position);

    // G4cout<<"t1 = "<<t1<<G4endl;

    E1 = fEnergyVector->GetLowEdgeEnergy(iTkin);

    // G4cout<<"E1 = "<<E1<<G4endl;

    W  = 1.0/(E2 - E1);
    W1 = (E2 - Tkin)*W;
    W2 = (Tkin - E1)*W;

    t = W1*t1 + W2*t2;
    */
  }
  return t;
}

Here is the call graph for this function:

G4double G4hhElastic::SampleTest

(

G4double

tMin

)

Definition at line 593 of file G4hhElastic.cc.

{
  G4int iTkin, iTransfer, iTmin;
  G4double t, position;
  // G4double qMin = std::sqrt(tMin);

  fTableT = fBankT[0];
  iTkin = 0;

  for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
  {
    // if( qMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
    if( tMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
  }
  iTmin = iTransfer-1;
  if(iTmin < 0 ) iTmin = 0;

  position = (*(*fTableT)(iTkin))(iTmin)*G4UniformRand();

  for( iTmin = 0; iTransfer < fBinT-1; iTransfer++)
  {
    if( position > (*(*fTableT)(iTkin))(iTransfer) ) break;
  }
  if (iTransfer >= fBinT-1) iTransfer = fBinT-2;

  t = GetTransfer(iTkin, iTransfer, position);

  return t;
}

Here is the call graph for this function:

void G4hhElastic::SetAlphaP ( G4double  a )

inline

Definition at line 186 of file G4hhElastic.hh.

{fAlphaP = a;};

void G4hhElastic::SetBq ( G4double  b )

inline

Definition at line 174 of file G4hhElastic.hh.

{fBq = b;};

Here is the caller graph for this function:

void G4hhElastic::SetBQ ( G4double  b )

inline

Definition at line 175 of file G4hhElastic.hh.

{fBQ = b;};

Here is the caller graph for this function:

void G4hhElastic::SetBqQ ( G4double  b )

inline

Definition at line 176 of file G4hhElastic.hh.

{fBqQ = b;};

void G4hhElastic::SetCofF2 ( G4double  f )

inline

Definition at line 191 of file G4hhElastic.hh.

{fCofF2 = f;};

Here is the caller graph for this function:

void G4hhElastic::SetCofF3 ( G4double  f )

inline

Definition at line 192 of file G4hhElastic.hh.

{fCofF3 = f;};

Here is the caller graph for this function:

void G4hhElastic::SetEta ( G4double  E )

inline

Definition at line 190 of file G4hhElastic.hh.

{fEta = E;};

Here is the caller graph for this function:

void G4hhElastic::SetImCof ( G4double  a )

inline

Definition at line 187 of file G4hhElastic.hh.

{fImCof = a;};

Here is the caller graph for this function:

void G4hhElastic::SetLambda ( G4double  L )

inline

Definition at line 189 of file G4hhElastic.hh.

{fLambda = L;};

Here is the caller graph for this function:

void G4hhElastic::SetParameters ( )

inline

Definition at line 273 of file G4hhElastic.hh.

{
  // masses

  fMq     = 0.36*CLHEP::GeV; // 0.441*GeV; // 0.36*GeV;
  fMQ     = 0.441*CLHEP::GeV;
  fMff2   = 0.26*CLHEP::GeV*CLHEP::GeV; // 0.25*GeV*GeV; // 0.5*GeV*GeV;

  fAlpha = 1./3.;
  fBeta  = 1. - fAlpha;

  fGamma = 1./2.; // 1./3.; // 
  fDelta = 1. - fGamma; // 1./2.;

  // radii and exp cof

  fRA       = 6.5/CLHEP::GeV; // 7.3/GeV;  //  3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
  fRq       = 0.173*fRA; // 2.4/GeV;
  fRQ       = 0.316*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
  fRB       = 6.5/CLHEP::GeV; // 7.3/GeV;  //  3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
  fRg       = 0.173*fRA; // 2.4/GeV;
  fRG       = 0.173*fRA; // 1./GeV; // 2./GeV; // 1./GeV;

  fAlphaP   = 0.15/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
  fLambda   = 0.25*fRA*fRA; // 0.25
  fEta      = 0.25*fRB*fRB; // 0.25
  fImCof    = 6.5;
  fCofF2 = 1.;
  fCofF3 = 1.;

  fBq = 0.02; // 0.21; // 1./3.;
  fBQ = 1. + fBq - 2*std::sqrt(fBq); // 1 - fBq; // 2./3.;
  fBqQ = std::sqrt(fBq*fBQ);

  fLambdaFF = 1.5/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
  fSo     = 1.*CLHEP::GeV*CLHEP::GeV;
  fQcof   = 0.009*CLHEP::GeV;
  fExpSlope = 19.9/CLHEP::GeV/CLHEP::GeV;
}

Here is the caller graph for this function:

void G4hhElastic::SetParametersCMS ( G4double  plab )

inline

Definition at line 318 of file G4hhElastic.hh.

{
  G4int i;
  G4double trMass = 900.*CLHEP::MeV, Tkin;
  G4double sl, sh, ds, rAl, rAh, drA, rBl, rBh, drB, bql, bqh, dbq, bQl, bQh, dbQ, cIl, cIh, dcI;

  Tkin = std::sqrt(fMassProj*fMassProj + plab*plab) - fMassProj;

  G4DynamicParticle*  theDynamicParticle = new G4DynamicParticle(fProjectile,
                                              G4ParticleMomentum(0.,0.,1.),
                                              Tkin);
  fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, fTarget );

  delete theDynamicParticle;

  fSpp  = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
  fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);

  G4double sCMS = std::sqrt(fSpp);

  if( fMassProj > trMass ) // p,n,pb on p
  {
    this->SetCofF2(1.);
    this->SetCofF3(1.);
    fGamma = 1./3.; // 1./3.; // 
    fDelta = 1. - fGamma; // 1./2.;

    if( sCMS <= theNuclNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
    {
      this->SetRA(theNuclNuclData[0][1]/CLHEP::GeV,0.173,0.316);
      this->SetRB(theNuclNuclData[0][2]/CLHEP::GeV,0.173,0.316); 

      this->SetBq(theNuclNuclData[0][3]);
      this->SetBQ(theNuclNuclData[0][4]);
      this->SetImCof(theNuclNuclData[0][5]);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(0.25*this->GetRB()*this->GetRB());
    }
    else if( sCMS >= theNuclNuclData[17][0]*CLHEP::GeV ) // high edge, as s=7000 ???
    {
      this->SetRA(theNuclNuclData[17][1]/CLHEP::GeV,0.173,0.316);
      this->SetRB(theNuclNuclData[17][2]/CLHEP::GeV,0.173,0.316); 

      this->SetBq(theNuclNuclData[17][3]);
      this->SetBQ(theNuclNuclData[17][4]);
      this->SetImCof(theNuclNuclData[17][5]);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(0.25*this->GetRB()*this->GetRB());
    }
    else // in approximation between array points
    {
      for( i = 0; i < 18; i++ ) if( sCMS <= theNuclNuclData[i][0]*CLHEP::GeV ) break;
      if( i == 0 ) i++;
      if( i == 18 ) i--;

      sl = theNuclNuclData[i-1][0]*CLHEP::GeV;
      sh = theNuclNuclData[i][0]*CLHEP::GeV;
      ds = (sCMS - sl)/(sh - sl);

      rAl = theNuclNuclData[i-1][1]/CLHEP::GeV;
      rAh = theNuclNuclData[i][1]/CLHEP::GeV;
      drA = rAh - rAl;

      rBl = theNuclNuclData[i-1][2]/CLHEP::GeV;
      rBh = theNuclNuclData[i][2]/CLHEP::GeV;
      drB = rBh - rBl;

      bql = theNuclNuclData[i-1][3];
      bqh = theNuclNuclData[i][3];
      dbq = bqh - bql;

      bQl = theNuclNuclData[i-1][4];
      bQh = theNuclNuclData[i][4];
      dbQ = bQh - bQl;

      cIl = theNuclNuclData[i-1][5];
      cIh = theNuclNuclData[i][5];
      dcI = cIh - cIl;

      this->SetRA(rAl+drA*ds,0.173,0.316);
      this->SetRB(rBl+drB*ds,0.173,0.316); 

      this->SetBq(bql+dbq*ds);
      this->SetBQ(bQl+dbQ*ds);
      this->SetImCof(cIl+dcI*ds);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(0.25*this->GetRB()*this->GetRB());
    }
  }
  else // pi, K
  {
    this->SetCofF2(1.);
    this->SetCofF3(-1.);
    fGamma = 1./2.; // 1./3.; // 
    fDelta = 1. - fGamma; // 1./2.;
  
    if( sCMS <= thePiKaNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
    {
      this->SetRA(thePiKaNuclData[0][1]/CLHEP::GeV,0.173,0.316);
      this->SetRB(thePiKaNuclData[0][2]/CLHEP::GeV,0.173,0.173); 

      this->SetBq(thePiKaNuclData[0][3]);
      this->SetBQ(thePiKaNuclData[0][4]);
      this->SetImCof(thePiKaNuclData[0][5]);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(this->GetRB()*this->GetRB()/6.);
    }
    else if( sCMS >= thePiKaNuclData[7][0]*CLHEP::GeV ) // high edge, as s=7000 ???
    {
      this->SetRA(thePiKaNuclData[7][1]/CLHEP::GeV,0.173,0.316);
      this->SetRB(thePiKaNuclData[7][2]/CLHEP::GeV,0.173,0.173); 

      this->SetBq(thePiKaNuclData[7][3]);
      this->SetBQ(thePiKaNuclData[7][4]);
      this->SetImCof(thePiKaNuclData[7][5]);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(this->GetRB()*this->GetRB()/6.);
    }
    else // in approximation between array points
    {
      for( i = 0; i < 8; i++ ) if( sCMS <= thePiKaNuclData[i][0]*CLHEP::GeV ) break;
      if( i == 0 ) i++;
      if( i == 8 ) i--;

      sl = thePiKaNuclData[i-1][0]*CLHEP::GeV;
      sh = thePiKaNuclData[i][0]*CLHEP::GeV;
      ds = (sCMS - sl)/(sh - sl);

      rAl = thePiKaNuclData[i-1][1]/CLHEP::GeV;
      rAh = thePiKaNuclData[i][1]/CLHEP::GeV;
      drA = rAh - rAl;

      rBl = thePiKaNuclData[i-1][2]/CLHEP::GeV;
      rBh = thePiKaNuclData[i][2]/CLHEP::GeV;
      drB = rBh - rBl;

      bql = thePiKaNuclData[i-1][3];
      bqh = thePiKaNuclData[i][3];
      dbq = bqh - bql;

      bQl = thePiKaNuclData[i-1][4];
      bQh = thePiKaNuclData[i][4];
      dbQ = bQh - bQl;

      cIl = thePiKaNuclData[i-1][5];
      cIh = thePiKaNuclData[i][5];
      dcI = cIh - cIl;

      this->SetRA(rAl+drA*ds,0.173,0.316);
      this->SetRB(rBl+drB*ds,0.173,0.173); 

      this->SetBq(bql+dbq*ds);
      this->SetBQ(bQl+dbQ*ds);
      this->SetImCof(cIl+dcI*ds);

      this->SetLambda(0.25*this->GetRA()*this->GetRA());
      this->SetEta(this->GetRB()*this->GetRB()/6.);
    }
   }
  return;
}

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.

{
  fRA = rA;
  fRq = fRA*pq;
  fRQ = fRA*pQ;
}

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.

{
  fRB = rB;
  fRg = fRB*pg;
  fRG = fRB*pG;
}

Here is the caller graph for this function:

void G4hhElastic::SetSigmaTot ( G4double  stot )

inline

Definition at line 166 of file G4hhElastic.hh.

{fSigmaTot = stot;};

void G4hhElastic::SetSpp ( G4double  spp )

inline

Definition at line 167 of file G4hhElastic.hh.

{fSpp = spp;};

---

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

• geant4.10.03.p01/source/processes/hadronic/models/coherent_elastic/include/G4hhElastic.hh
• geant4.10.03.p01/source/processes/hadronic/models/coherent_elastic/src/G4hhElastic.cc