Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:58:24 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 //
0029 // G4 Model: hadron diffraction elastic scattering with 4-momentum balance 
0030 //
0031 // Class Description
0032 // Final state production model for hadron-hadron elastic scattering 
0033 // in the framework of quark-diquark model with springy Pomeron. 
0034 // Projectiles are proton, neutron, pions, kaons. 
0035 // Targets are proton (and neutron). 
0036 // Class Description - End
0037 //
0038 // 02.05.14 V. Grichine - 1-st implementation
0039 // 10.10.14 V. Grichine - change to combine with low mass diffraction
0040 
0041 #ifndef G4hhElastic_h
0042 #define G4hhElastic_h 1
0043  
0044 #include "globals.hh"
0045 #include <complex>
0046 #include "G4Integrator.hh"
0047 
0048 #include "G4HadronElastic.hh"
0049 #include "G4HadProjectile.hh"
0050 #include "G4Nucleus.hh"
0051 #include "G4HadronNucleonXsc.hh"
0052 
0053 #include "G4Exp.hh"
0054 #include "G4Log.hh"
0055 
0056 
0057 class G4ParticleDefinition;
0058 class G4PhysicsTable;
0059 class G4PhysicsLogVector;
0060 
0061 class G4hhElastic  : public G4HadronElastic  
0062 {
0063 public:
0064   // PL constructor
0065   G4hhElastic();
0066   // test constructor
0067   G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile, 
0068                     G4double plab );
0069   // constructor used for low mass diffraction
0070   G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile);
0071 
0072   virtual ~G4hhElastic();
0073 
0074   virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/, 
0075                   G4Nucleus & /*targetNucleus*/);
0076 
0077 
0078   void Initialise();
0079 
0080   void BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile); // , G4double plab );
0081   void BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile, 
0082                     G4double plab );
0083 
0084   G4double SampleInvariantT( const G4ParticleDefinition* p,  G4double plab, G4int, G4int);
0085   G4double SampleBisectionalT( const G4ParticleDefinition* p, G4double plab);
0086 
0087   G4double SampleTest(G4double tMin ); // const G4ParticleDefinition* p, );
0088 
0089   G4double GetTransfer( G4int iMomentum, G4int iTransfer, G4double position );
0090  
0091 private:
0092 
0093 
0094   G4ParticleDefinition* fTarget;
0095   G4ParticleDefinition* fProjectile;
0096   G4ParticleDefinition* theProton;
0097   G4ParticleDefinition* theNeutron;
0098   G4ParticleDefinition* thePionPlus;
0099   G4ParticleDefinition* thePionMinus;
0100 
0101 
0102   G4double lowEnergyRecoilLimit;  
0103   G4double lowEnergyLimitHE;  
0104   G4double lowEnergyLimitQ;  
0105   G4double lowestEnergyLimit;  
0106   G4double plabLowLimit;
0107 
0108   G4int fEnergyBin;
0109   G4int fBinT;
0110 
0111   G4PhysicsLogVector*           fEnergyVector;
0112   G4PhysicsTable*               fTableT;
0113   std::vector<G4PhysicsTable*>  fBankT;
0114 
0115 
0116   // Gauss model parameters
0117 
0118 
0119   G4double fMff2;
0120   G4double fMQ;
0121   G4double fMq;
0122 
0123   G4double fMassTarg;  // ~ A
0124   G4double fMassProj;  // ~ B
0125   G4double fMassSum2;
0126   G4double fMassDif2;
0127 
0128 
0129   G4double fRA; // hadron A
0130   G4double fRQ;
0131   G4double fRq;
0132   G4double fAlpha;
0133   G4double fBeta;
0134 
0135   G4double fRB; // hadron B
0136   G4double fRG;
0137   G4double fRg;
0138   G4double fGamma;
0139   G4double fDelta;
0140 
0141   G4double fAlphaP;
0142   G4double fLambdaFF;
0143   G4double fLambda;
0144   G4double fEta;
0145   G4double fImCof;
0146   G4double fCofF2;
0147   G4double fCofF3;
0148   G4double fRhoReIm;
0149   G4double fExpSlope;
0150   G4double fSo;
0151 
0152   G4double fSigmaTot;
0153   G4double fBq;
0154   G4double fBQ;
0155   G4double fBqQ;
0156   G4double fOptRatio;
0157   G4double fSpp;
0158   G4double fPcms;
0159   G4double fQcof;        // q prime when integrate
0160 
0161 public:  // Gauss model methods
0162 
0163   void SetParameters();
0164   void SetSigmaTot(G4double stot){fSigmaTot = stot;};
0165   void SetSpp(G4double spp){fSpp = spp;};
0166   G4double GetSpp(){return fSpp;};
0167   void SetParametersCMS(G4double plab);
0168 
0169   G4double GetBq(){ return fBq;};
0170   G4double GetBQ(){ return fBQ;};
0171   G4double GetBqQ(){ return fBqQ;};
0172   void SetBq(G4double b){fBq = b;};
0173   void SetBQ(G4double b){fBQ = b;};
0174   void SetBqQ(G4double b){fBqQ = b;};
0175   G4double GetRhoReIm(){ return fRhoReIm;};
0176 
0177   void CalculateBQ(G4double b);
0178   void CalculateBqQ13(G4double b);
0179   void CalculateBqQ12(G4double b);
0180   void CalculateBqQ123(G4double b);
0181   void SetRA(G4double rn, G4double pq, G4double pQ);
0182   void SetRB(G4double rn, G4double pq, G4double pQ);
0183 
0184   void SetAlphaP(G4double a){fAlphaP = a;};
0185   void SetImCof(G4double a){fImCof = a;};
0186   G4double GetImCof(){return fImCof;};
0187   void SetLambda(G4double L){fLambda = L;};
0188   void SetEta(G4double E){fEta = E;};
0189   void SetCofF2(G4double f){fCofF2 = f;};
0190   void SetCofF3(G4double f){fCofF3 = f;};
0191   G4double GetCofF2(){return fCofF2;};
0192   G4double GetCofF3(){return fCofF3;};
0193 
0194   G4double GetRA(){ return fRA;};
0195   G4double GetRq(){ return fRq;};
0196   G4double GetRQ(){ return fRQ;};
0197 
0198   G4double GetRB(){ return fRB;};
0199   G4double GetRg(){ return fRg;};
0200   G4double GetRG(){ return fRG;};
0201 
0202   // FqQgG stuff
0203 
0204   G4complex Pomeron();
0205 
0206   G4complex Phi13();
0207   G4complex Phi14();
0208   G4complex Phi23();
0209   G4complex Phi24();
0210 
0211   G4complex GetF1qQgG(G4double qp);
0212   G4double GetdsdtF1qQgG(G4double s, G4double q);
0213   G4complex GetF2qQgG(G4double qp);
0214   G4double GetdsdtF12qQgG(G4double s, G4double q);
0215   G4complex GetF3qQgG(G4double qp);
0216   G4double GetdsdtF123qQgG(G4double q); // sampling ds/dt
0217   G4double GetdsdtF13qQG(G4double s, G4double q);
0218 
0219 
0220   // F123 stuff
0221 
0222   G4complex GetAqq();
0223   G4complex GetAQQ();
0224   G4complex GetAqQ();
0225 
0226   G4double GetCofS1();
0227   G4double GetCofS2();
0228   G4double GetCofS3();
0229   G4double GetOpticalRatio();
0230 
0231   G4complex GetF1(G4double qp);
0232   G4double GetdsdtF1(G4double s, G4double q);
0233   G4complex GetF2(G4double qp);
0234   G4double GetdsdtF12(G4double s, G4double q);
0235   G4complex GetF3(G4double qp);
0236   G4double GetdsdtF123(G4double q); // sampling ds/dt
0237   G4double GetExpRatioF123(G4double s, G4double q);
0238 
0239   // parameter arrays
0240 
0241 private:
0242 
0243   G4int fInTkin;
0244   G4double fOldTkin;
0245   static const G4double theNuclNuclData[19][6];
0246   static const G4double thePiKaNuclData[8][6];
0247   G4HadronNucleonXsc* fHadrNuclXsc;
0248 };
0249 
0250 //////////////////////////////////////////////////////////////////////
0251 //////////////////////////////////////////////////////////////////////
0252 ////////////////////////////////////////////////////////////////////////
0253 
0254 
0255 
0256 inline G4bool G4hhElastic::IsApplicable(const G4HadProjectile & projectile, 
0257                   G4Nucleus & nucleus)
0258 {
0259   if( ( projectile.GetDefinition() == G4Proton::Proton() ||
0260         projectile.GetDefinition() == G4Neutron::Neutron() ||
0261         projectile.GetDefinition() == G4PionPlus::PionPlus() ||
0262         projectile.GetDefinition() == G4PionMinus::PionMinus() ||
0263         projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
0264         projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
0265 
0266         nucleus.GetZ_asInt() < 2 ) return true;
0267   else                              return false;
0268 }
0269 
0270 
0271 inline void G4hhElastic::SetParameters()
0272 {
0273   // masses
0274 
0275   fMq     = 0.36*CLHEP::GeV; // 0.441*GeV; // 0.36*GeV;
0276   fMQ     = 0.441*CLHEP::GeV;
0277   fMff2   = 0.26*CLHEP::GeV*CLHEP::GeV; // 0.25*GeV*GeV; // 0.5*GeV*GeV;
0278 
0279   fAlpha = 1./3.;
0280   fBeta  = 1. - fAlpha;
0281 
0282   fGamma = 1./2.; // 1./3.; // 
0283   fDelta = 1. - fGamma; // 1./2.;
0284 
0285   // radii and exp cof
0286 
0287   fRA       = 6.5/CLHEP::GeV; // 7.3/GeV;  //  3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
0288   fRq       = 0.173*fRA; // 2.4/GeV;
0289   fRQ       = 0.316*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
0290   fRB       = 6.5/CLHEP::GeV; // 7.3/GeV;  //  3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
0291   fRg       = 0.173*fRA; // 2.4/GeV;
0292   fRG       = 0.173*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
0293 
0294   fAlphaP   = 0.15/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
0295   fLambda   = 0.25*fRA*fRA; // 0.25
0296   fEta      = 0.25*fRB*fRB; // 0.25
0297   fImCof    = 6.5;
0298   fCofF2 = 1.;
0299   fCofF3 = 1.;
0300 
0301   fBq = 0.02; // 0.21; // 1./3.;
0302   fBQ = 1. + fBq - 2*std::sqrt(fBq); // 1 - fBq; // 2./3.;
0303   fBqQ = std::sqrt(fBq*fBQ);
0304 
0305   fLambdaFF = 1.5/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
0306   fSo     = 1.*CLHEP::GeV*CLHEP::GeV;
0307   fQcof   = 0.009*CLHEP::GeV;
0308   fExpSlope = 19.9/CLHEP::GeV/CLHEP::GeV;
0309 }
0310 
0311 
0312 ////////////////////////////////////////////////////////////////////////
0313 //
0314 // Set target and projectile masses and calculate mass sum and difference squared for Pcms
0315 
0316 inline void G4hhElastic::SetParametersCMS(G4double plab)
0317 {
0318   G4int i;
0319   G4double trMass = 900.*CLHEP::MeV, Tkin;
0320   G4double sl, sh, ds, rAl, rAh, drA, rBl, rBh, drB, bql, bqh, dbq, bQl, bQh, dbQ, cIl, cIh, dcI;
0321 
0322   Tkin = std::sqrt(fMassProj*fMassProj + plab*plab) - fMassProj;
0323 
0324   G4DynamicParticle*  theDynamicParticle = new G4DynamicParticle(fProjectile,
0325                                               G4ParticleMomentum(0.,0.,1.),
0326                                               Tkin);
0327   fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, fTarget );
0328 
0329   delete theDynamicParticle;
0330 
0331   fSpp  = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
0332   fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0333 
0334   G4double sCMS = std::sqrt(fSpp);
0335 
0336   if( fMassProj > trMass ) // p,n,pb on p
0337   {
0338     this->SetCofF2(1.);
0339     this->SetCofF3(1.);
0340     fGamma = 1./3.; // 1./3.; // 
0341     fDelta = 1. - fGamma; // 1./2.;
0342 
0343     if( sCMS <= theNuclNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
0344     {
0345       this->SetRA(theNuclNuclData[0][1]/CLHEP::GeV,0.173,0.316);
0346       this->SetRB(theNuclNuclData[0][2]/CLHEP::GeV,0.173,0.316); 
0347 
0348       this->SetBq(theNuclNuclData[0][3]);
0349       this->SetBQ(theNuclNuclData[0][4]);
0350       this->SetImCof(theNuclNuclData[0][5]);
0351 
0352       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0353       this->SetEta(0.25*this->GetRB()*this->GetRB());
0354     }
0355     else if( sCMS >= theNuclNuclData[17][0]*CLHEP::GeV ) // high edge, as s=7000 ???
0356     {
0357       this->SetRA(theNuclNuclData[17][1]/CLHEP::GeV,0.173,0.316);
0358       this->SetRB(theNuclNuclData[17][2]/CLHEP::GeV,0.173,0.316); 
0359 
0360       this->SetBq(theNuclNuclData[17][3]);
0361       this->SetBQ(theNuclNuclData[17][4]);
0362       this->SetImCof(theNuclNuclData[17][5]);
0363 
0364       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0365       this->SetEta(0.25*this->GetRB()*this->GetRB());
0366     }
0367     else // in approximation between array points
0368     {
0369       for( i = 0; i < 19; i++ ) if( sCMS <= theNuclNuclData[i][0]*CLHEP::GeV ) break;
0370       if( i == 0 ) i++;
0371       if( i == 19 ) i--;
0372 
0373       sl = theNuclNuclData[i-1][0]*CLHEP::GeV;
0374       sh = theNuclNuclData[i][0]*CLHEP::GeV;
0375       ds = (sCMS - sl)/(sh - sl);
0376 
0377       rAl = theNuclNuclData[i-1][1]/CLHEP::GeV;
0378       rAh = theNuclNuclData[i][1]/CLHEP::GeV;
0379       drA = rAh - rAl;
0380 
0381       rBl = theNuclNuclData[i-1][2]/CLHEP::GeV;
0382       rBh = theNuclNuclData[i][2]/CLHEP::GeV;
0383       drB = rBh - rBl;
0384 
0385       bql = theNuclNuclData[i-1][3];
0386       bqh = theNuclNuclData[i][3];
0387       dbq = bqh - bql;
0388 
0389       bQl = theNuclNuclData[i-1][4];
0390       bQh = theNuclNuclData[i][4];
0391       dbQ = bQh - bQl;
0392 
0393       cIl = theNuclNuclData[i-1][5];
0394       cIh = theNuclNuclData[i][5];
0395       dcI = cIh - cIl;
0396 
0397       this->SetRA(rAl+drA*ds,0.173,0.316);
0398       this->SetRB(rBl+drB*ds,0.173,0.316); 
0399 
0400       this->SetBq(bql+dbq*ds);
0401       this->SetBQ(bQl+dbQ*ds);
0402       this->SetImCof(cIl+dcI*ds);
0403 
0404       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0405       this->SetEta(0.25*this->GetRB()*this->GetRB());
0406     }
0407   }
0408   else // pi, K
0409   {
0410     this->SetCofF2(1.);
0411     this->SetCofF3(-1.);
0412     fGamma = 1./2.; // 1./3.; // 
0413     fDelta = 1. - fGamma; // 1./2.;
0414   
0415     if( sCMS <= thePiKaNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
0416     {
0417       this->SetRA(thePiKaNuclData[0][1]/CLHEP::GeV,0.173,0.316);
0418       this->SetRB(thePiKaNuclData[0][2]/CLHEP::GeV,0.173,0.173); 
0419 
0420       this->SetBq(thePiKaNuclData[0][3]);
0421       this->SetBQ(thePiKaNuclData[0][4]);
0422       this->SetImCof(thePiKaNuclData[0][5]);
0423 
0424       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0425       this->SetEta(this->GetRB()*this->GetRB()/6.);
0426     }
0427     else if( sCMS >= thePiKaNuclData[7][0]*CLHEP::GeV ) // high edge, as s=7000 ???
0428     {
0429       this->SetRA(thePiKaNuclData[7][1]/CLHEP::GeV,0.173,0.316);
0430       this->SetRB(thePiKaNuclData[7][2]/CLHEP::GeV,0.173,0.173); 
0431 
0432       this->SetBq(thePiKaNuclData[7][3]);
0433       this->SetBQ(thePiKaNuclData[7][4]);
0434       this->SetImCof(thePiKaNuclData[7][5]);
0435 
0436       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0437       this->SetEta(this->GetRB()*this->GetRB()/6.);
0438     }
0439     else // in approximation between array points
0440     {
0441       for( i = 0; i < 8; i++ ) if( sCMS <= thePiKaNuclData[i][0]*CLHEP::GeV ) break;
0442       if( i == 0 ) i++;
0443       if( i == 8 ) i--;
0444 
0445       sl = thePiKaNuclData[i-1][0]*CLHEP::GeV;
0446       sh = thePiKaNuclData[i][0]*CLHEP::GeV;
0447       ds = (sCMS - sl)/(sh - sl);
0448 
0449       rAl = thePiKaNuclData[i-1][1]/CLHEP::GeV;
0450       rAh = thePiKaNuclData[i][1]/CLHEP::GeV;
0451       drA = rAh - rAl;
0452 
0453       rBl = thePiKaNuclData[i-1][2]/CLHEP::GeV;
0454       rBh = thePiKaNuclData[i][2]/CLHEP::GeV;
0455       drB = rBh - rBl;
0456 
0457       bql = thePiKaNuclData[i-1][3];
0458       bqh = thePiKaNuclData[i][3];
0459       dbq = bqh - bql;
0460 
0461       bQl = thePiKaNuclData[i-1][4];
0462       bQh = thePiKaNuclData[i][4];
0463       dbQ = bQh - bQl;
0464 
0465       cIl = thePiKaNuclData[i-1][5];
0466       cIh = thePiKaNuclData[i][5];
0467       dcI = cIh - cIl;
0468 
0469       this->SetRA(rAl+drA*ds,0.173,0.316);
0470       this->SetRB(rBl+drB*ds,0.173,0.173); 
0471 
0472       this->SetBq(bql+dbq*ds);
0473       this->SetBQ(bQl+dbQ*ds);
0474       this->SetImCof(cIl+dcI*ds);
0475 
0476       this->SetLambda(0.25*this->GetRA()*this->GetRA());
0477       this->SetEta(this->GetRB()*this->GetRB()/6.);
0478     }
0479    }
0480   return;
0481 }
0482 
0483 /////////////////////////////////////////////////////
0484 //
0485 // RA for qQ
0486 
0487 inline void G4hhElastic::SetRA(G4double rA, G4double pq, G4double pQ)
0488 {
0489   fRA = rA;
0490   fRq = fRA*pq;
0491   fRQ = fRA*pQ;
0492 }
0493 
0494 /////////////////////////////////////////////////////
0495 //
0496 // RB for gG
0497 
0498 inline void G4hhElastic::SetRB(G4double rB, G4double pg, G4double pG)
0499 {
0500   fRB = rB;
0501   fRg = fRB*pg;
0502   fRG = fRB*pG;
0503 }
0504 
0505 ////////////////////////////////////////////////////
0506 //
0507 // Returns Pomeron parametrization with Im part modified, *= fImCof
0508 
0509 inline G4complex G4hhElastic::Pomeron()
0510 {
0511   G4double re, im;
0512 
0513   re = fAlphaP*G4Log(fSpp/fSo);
0514   im = -0.5*fAlphaP*fImCof*CLHEP::pi;
0515   return G4complex(re,im);
0516 } 
0517 
0518 //////////////////////////////////////////////////
0519 
0520 inline G4complex G4hhElastic::Phi13()
0521 {
0522   G4double re = (fRq*fRq + fRg*fRg)/16.;
0523   G4complex result(re,0.);
0524   result += Pomeron();
0525   return result;
0526 }
0527 
0528 //////////////////////////////////////////////////
0529 
0530 inline G4complex G4hhElastic::Phi14()
0531 {
0532   G4double re = (fRq*fRq + fRG*fRG)/16.;
0533   G4complex result(re,0.);
0534   result += Pomeron();
0535   return result;
0536 }
0537 
0538 //////////////////////////////////////////////////
0539 
0540 inline G4complex G4hhElastic::Phi23()
0541 {
0542   G4double re = (fRQ*fRQ + fRg*fRg)/16.;
0543   G4complex result(re,0.);
0544   result += Pomeron();
0545   return result;
0546 }
0547 
0548 //////////////////////////////////////////////////
0549 
0550 inline G4complex G4hhElastic::Phi24()
0551 {
0552   G4double re = (fRQ*fRQ + fRG*fRG)/16.;
0553   G4complex result(re,0.);
0554   result += Pomeron();
0555   return result;
0556 }
0557 
0558 /////////////////////////////////////////////////////
0559 //
0560 // F1, case qQ-gG
0561 
0562 inline G4complex G4hhElastic::GetF1qQgG(G4double t)
0563 {
0564   G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0565   G4double k = p/CLHEP::hbarc;
0566 
0567   G4complex exp13 = fBq*std::exp(-(Phi13() + fBeta*fBeta*fLambda + fDelta*fDelta*fEta)*t);
0568   G4complex exp14 = fBq*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
0569   G4complex exp23 = fBQ*std::exp(-(Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta)*t);
0570   G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
0571 
0572   G4complex res  = exp13 + exp14 + exp23 + exp24;
0573 
0574   res *=  0.25*k*fSigmaTot/CLHEP::pi;
0575   res *= G4complex(0.,1.);
0576 
0577   return res;
0578 }
0579 
0580 /////////////////////////////////////////////////////
0581 //
0582 // 
0583 
0584 inline G4double G4hhElastic::GetdsdtF13qQG(G4double spp, G4double t)
0585 {
0586   fSpp = spp;
0587   G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0588   G4double k = p/CLHEP::hbarc;
0589 
0590   G4complex exp14 = fBqQ*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
0591   G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
0592 
0593   G4complex F1  = exp14 + exp24;
0594 
0595   F1 *=  0.25*k*fSigmaTot/CLHEP::pi;
0596   F1 *= G4complex(0.,1.);
0597 
0598   // 1424
0599 
0600   G4complex z1424  = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
0601             z1424 /= Phi14() + Phi24() + fLambda;
0602             z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
0603 
0604  
0605   G4complex exp1424  = std::exp(-z1424*t);
0606             exp1424 /= Phi14() + Phi24() + fLambda;
0607 
0608   G4complex F3  = fBqQ*fBQ*exp1424;
0609 
0610 
0611   F3 *=  0.25*k/CLHEP::pi;
0612   F3 *= G4complex(0.,1.);
0613   F3 *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
0614 
0615   G4complex F13  = F1 - F3;
0616 
0617   G4double dsdt = CLHEP::pi/p/p;
0618   dsdt *= real(F13)*real(F13) + imag(F13)*imag(F13);  
0619 
0620   return dsdt;
0621 }
0622 
0623 //////////////////////////////////////////////////////////////
0624 //
0625 // dsigma/dt(s,t) F1qQgG
0626 
0627 inline  G4double G4hhElastic::GetdsdtF1qQgG(G4double spp, G4double t)
0628 {
0629   fSpp = spp;
0630   G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0631 
0632   G4complex F1 = GetF1qQgG(t);
0633 
0634   G4double dsdt = CLHEP::pi/p/p;
0635   dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
0636   return dsdt;
0637 }
0638 
0639 /////////////////////////////////////////////////////
0640 //
0641 // 
0642 
0643 inline G4complex G4hhElastic::GetF2qQgG(G4double t)
0644 {
0645   G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0646   G4double k = p/CLHEP::hbarc;
0647 
0648   G4complex z1324  = -(Phi24() + fAlpha*fLambda + fGamma*fEta)*(Phi24() + fAlpha*fLambda + fGamma*fEta);
0649             z1324 /= Phi13() + Phi24() + fLambda + fEta;
0650             z1324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
0651 
0652   G4complex exp1324  = std::exp(-z1324*t);
0653             exp1324 /= Phi13() + Phi24() + fLambda + fEta;
0654 
0655   G4complex z1423 = -(Phi23() + fAlpha*fLambda + fDelta*fEta)*(Phi24() + fAlpha*fLambda + fDelta*fEta);;
0656             z1423 /= Phi14() + Phi23() + fLambda + fEta;
0657             z1423 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
0658 
0659   G4complex exp1423 = std::exp(-z1423*t);
0660             exp1423 /= Phi14() + Phi23() + fLambda + fEta; 
0661 
0662   G4complex res  = exp1324 + exp1423;
0663 
0664 
0665   res *=  0.25*k/CLHEP::pi;
0666   res *= G4complex(0.,1.);
0667   res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc); // or 4. ???
0668 
0669   return res;
0670 }
0671 
0672 //////////////////////////////////////////////////////////////
0673 //
0674 // dsigma/dt(s,t) F12
0675 
0676 inline  G4double G4hhElastic::GetdsdtF12qQgG( G4double spp, G4double t)
0677 {
0678   fSpp = spp;
0679   G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0680 
0681   G4complex F12 = GetF1qQgG(t) - GetF2qQgG(t);
0682 
0683   G4double dsdt = CLHEP::pi/p/p;
0684   dsdt *= real(F12)*real(F12) + imag(F12)*imag(F12);  
0685   return dsdt;
0686 }
0687 
0688 /////////////////////////////////////////////////////
0689 //
0690 // 
0691 
0692 inline G4complex G4hhElastic::GetF3qQgG(G4double t)
0693 {
0694   G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
0695   G4double k = p/CLHEP::hbarc;
0696 
0697             // 1314
0698 
0699   G4complex z1314  = -(Phi14() + fGamma*fEta)*(Phi14() + fGamma*fEta);
0700             z1314 /= Phi13() + Phi14() + fEta;
0701             z1314 += Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta;
0702 
0703   G4complex exp1314  = std::exp(-z1314*t);
0704             exp1314 /= Phi13() + Phi14() + fEta;
0705 
0706         // 2324
0707 
0708   G4complex z2324 = -(Phi24() + fGamma*fEta)*(Phi24() + fGamma*fEta);;
0709             z2324 /= Phi24() + Phi23() + fEta;
0710             z2324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
0711 
0712   G4complex exp2324 = std::exp(-z2324*t);
0713             exp2324 /= Phi24() + Phi23() + fEta; 
0714 
0715         // 1323
0716 
0717   G4complex z1323  = -(Phi23() + fAlpha*fLambda)*(Phi23() + fAlpha*fLambda);
0718             z1323 /= Phi13() + Phi23() + fLambda;
0719             z1323 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
0720 
0721   G4complex exp1323  = std::exp(-z1323*t);
0722             exp1323 /= Phi13() + Phi23() + fLambda;
0723 
0724         // 1424
0725 
0726   G4complex z1424  = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
0727             z1424 /= Phi14() + Phi24() + fLambda;
0728             z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
0729 
0730   G4complex exp1424  = std::exp(-z1424*t);
0731             exp1424 /= Phi14() + Phi24() + fLambda;
0732 
0733         G4complex res  = fBq*fBq*exp1314 + fBQ*fBQ*exp2324 + fBq*fBQ*exp1323 + fBq*fBQ*exp1424;
0734 
0735   res *=  0.25*k/CLHEP::pi;
0736   res *= G4complex(0.,1.);
0737   res *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
0738 
0739   return res;
0740 }
0741 
0742 //////////////////////////////////////////////////////////////
0743 //
0744 // dsigma/dt(s,t) F123 sampling ds/dt
0745 
0746 inline  G4double G4hhElastic::GetdsdtF123qQgG(G4double t)
0747 {
0748   G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp );
0749 
0750   G4complex F123  = GetF1qQgG(t); // - fCofF2*GetF2qQgG(t) - fCofF3*GetF3qQgG(t);
0751             F123 -= fCofF2*GetF2qQgG(t);
0752             F123 -= fCofF3*GetF3qQgG(t);
0753 
0754   G4double dsdt = CLHEP::pi/p/p;
0755   dsdt *= real(F123)*real(F123) + imag(F123)*imag(F123);  
0756   return dsdt;
0757 }
0758 
0759 /////////////////////////////////////////////////////
0760 //
0761 // Set fBqQ at a given fBQ=b2 according to the optical theorem,qQ-G
0762 
0763 inline void G4hhElastic::CalculateBqQ13(G4double b2)
0764 {
0765   fBQ = b2;
0766 
0767   G4complex z1424 = G4complex(1./8./CLHEP::pi,0.);
0768             z1424 /= Phi14() + Phi24() + fAlpha;
0769         G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
0770 
0771   fBqQ  = 1. - fBQ;
0772   fBQ /= 1. - fSigmaTot*fBQ*c1424; 
0773 
0774   G4cout<<"fSigmaTot*fBQ*c1424 = "<<fSigmaTot*fBQ*c1424<<G4endl;
0775 
0776   G4double ratio  = fBqQ + fBQ - fSigmaTot*fBqQ*fBQ*c1424; 
0777   G4cout<<"ratio = "<<ratio<<G4endl;
0778 
0779   return ;
0780 }
0781 
0782 /////////////////////////////////////////////////////
0783 //
0784 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2
0785 
0786 inline void G4hhElastic::CalculateBqQ12(G4double b1)
0787 {
0788   fBq = b1;
0789 
0790   G4complex z1324 = G4complex(1./8./CLHEP::pi,0.);
0791             z1324 /= Phi13() + Phi24() + fLambda + fEta;
0792         G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
0793 
0794   G4complex z1423 = G4complex(1./8./CLHEP::pi,0.);
0795             z1423 /= Phi14() + Phi23() + fLambda + fEta;
0796         G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
0797 
0798   fBQ  = 1. - 2.*fBq;
0799   fBQ /= 2. - fSigmaTot*fBq*(c1324+1423); 
0800 
0801   G4double ratio  = 2.*(fBq + fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423); 
0802   G4cout<<"ratio = "<<ratio<<G4endl;
0803 
0804   return ;
0805 }
0806 
0807 /////////////////////////////////////////////////////
0808 //
0809 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2-F3, 
0810 // simplified meson-barion case g=G=q
0811 
0812 inline void G4hhElastic::CalculateBqQ123(G4double b1)
0813 {
0814   fBq = b1;
0815 
0816   G4complex z1324 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
0817             z1324 /= Phi13() + Phi24() + fLambda + fEta;
0818         G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
0819 
0820   G4complex z1423 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
0821             z1423 /= Phi14() + Phi23() + fLambda + fEta;
0822         G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
0823 
0824   G4complex z1314 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
0825             z1314 /= Phi13() + Phi14() + fEta;
0826         G4double c1314 = real(z1314)/(CLHEP::hbarc*CLHEP::hbarc);
0827 
0828   G4complex z2324 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
0829             z2324 /= Phi23() + Phi24() + fEta;
0830         G4double c2324 = real(z2324)/(CLHEP::hbarc*CLHEP::hbarc);
0831 
0832   G4complex z1323 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
0833             z1323 /= Phi13() + Phi23() + fLambda;
0834         G4double c1323 = real(z1323)/(CLHEP::hbarc*CLHEP::hbarc);
0835 
0836   G4complex z1424 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
0837             z1424 /= Phi14() + Phi24() + fLambda;
0838         G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
0839 
0840   G4double A = fSigmaTot*c2324;
0841   G4double B = fSigmaTot*fBq*(c1324 + c1423 + c1323 + c1424) - 2.;
0842   G4double C = 1. + fSigmaTot*fBq*fBq*c1314 - 2*fBq;
0843   G4cout<<"A = "<<A<<"; B = "<<B<<"; C = "<<C<<G4endl;
0844   G4cout<<"determinant = "<<B*B-4.*A*C<<G4endl;
0845 
0846   G4double x1 = ( -B - std::sqrt(B*B-4.*A*C) )/2./A;
0847   G4double x2 = ( -B + std::sqrt(B*B-4.*A*C) )/2./A;
0848   G4cout<<"x1 = "<<x1<<"; x2 = "<<x2<<G4endl;
0849 
0850   if( B*B-4.*A*C < 1.e-6 ) fBQ = std::abs(-B/2./A);
0851   else if ( B < 0.)        fBQ = std::abs( ( -B - std::sqrt(B*B-4.*A*C) )/2./A);
0852   else                     fBQ  = std::abs( ( -B + std::sqrt(B*B-4.*A*C) )/2./A);
0853  
0854   fOptRatio  = 2*(fBq+fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423 + c1323 + c1424);
0855   fOptRatio -= fSigmaTot*fBq*fBq*c1314 + fSigmaTot*c2324*fBQ*fBQ;
0856   G4cout<<"BqQ123, fOptRatio = "<<fOptRatio<<G4endl;
0857 
0858   return ;
0859 }
0860 
0861 
0862 
0863 /////////////////////  F123 stuff hh-elastic, qQ-qQ ///////////////////
0864 //////////////////////////////////////////////////////////////
0865 /////////////////////////////////////////////////////////////////
0866 
0867 inline G4complex G4hhElastic::GetAqq()
0868 {
0869   G4double re, im;
0870 
0871   re = fRq*fRq/8. + fAlphaP*G4Log(fSpp/fSo) + 8.*fLambda/9.;
0872   im = -0.5*fAlphaP*fImCof*CLHEP::pi;
0873   return G4complex(re,im);
0874 }
0875 
0876 /////////////////////////////////////////////////////
0877 //
0878 //
0879 
0880 inline G4complex G4hhElastic::GetAQQ()
0881 {
0882   G4double re, im;
0883 
0884   re = fRQ*fRQ/8. + fAlphaP*G4Log(fSpp/fSo) + 2.*fLambda/9.;
0885   im = -0.5*fAlphaP*fImCof*CLHEP::pi;
0886   return G4complex(re,im);
0887 }
0888 
0889 /////////////////////////////////////////////////////
0890 //
0891 //
0892 
0893 inline G4complex G4hhElastic::GetAqQ()
0894 {
0895   G4complex z = 0.5*( GetAqq() + GetAQQ() );
0896   return z;
0897 }
0898 
0899 /////////////////////////////////////////////////////
0900 //
0901 //
0902 
0903 inline G4double G4hhElastic::GetCofS1()
0904 {
0905   G4complex z = 1./( GetAqQ() + 4.*fLambda/9. );
0906   G4double result = real(z);
0907   result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
0908   result *= fSigmaTot*fCofF2;
0909   return result;
0910 }
0911 
0912 /////////////////////////////////////////////////////
0913 //
0914 //
0915 
0916 inline G4double G4hhElastic::GetCofS2()
0917 {
0918   G4complex z = 1./( GetAqq() + GetAqQ() - 4.*fLambda/9. );
0919   G4double result = real(z);
0920   result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
0921   result *= fSigmaTot*fCofF3;
0922   return result;
0923 }
0924 
0925 /////////////////////////////////////////////////////
0926 //
0927 //
0928 
0929 inline G4double G4hhElastic::GetCofS3()
0930 {
0931   G4complex z = 1./( GetAQQ() + GetAqQ() + 2.*fLambda/9. );
0932   G4double result = real(z);
0933   result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
0934   result *= fSigmaTot*fCofF3;
0935   return result;
0936 }
0937 
0938 /////////////////////////////////////////////////////
0939 //
0940 //
0941 
0942 inline G4double G4hhElastic::GetOpticalRatio()
0943 {
0944   return fOptRatio;
0945   // G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
0946   // G4double result = fBq + fBQ + 2.*sqrtBqBQ - 1.;
0947   // result /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
0948   // return result;
0949 }
0950 
0951 
0952 
0953 /////////////////////////////////////////////////////
0954 //
0955 // Set fBQ at a given fBq=b according to the optical theorem
0956 
0957 inline void G4hhElastic::CalculateBQ(G4double b1)
0958 {
0959   fBq = b1;
0960   G4double s1 = GetCofS1();
0961   G4double s2 = GetCofS2();
0962   G4double s3 = GetCofS3();
0963   G4double sqrtBq = std::sqrt(fBq);
0964 
0965   // cofs of the fBQ 3rd equation
0966 
0967   G4double a = s3*sqrtBq;
0968   G4double b = s1*fBq - 1.;
0969   G4double c = (s2*fBq - 2.)*sqrtBq;
0970   G4double d = 1. - fBq;
0971 
0972   // cofs of the incomplete 3rd equation
0973 
0974   G4double p  = c/a;
0975            p -= b*b/a/a/3.;
0976   G4double q  = d/a;
0977            q -= b*c/a/a/3.;
0978            q += 2*b*b*b/a/a/a/27.;
0979 
0980   // cofs for the incomplete colutions
0981 
0982   G4double D  = p*p*p/3./3./3.;
0983            D += q*q/2./2.;
0984   G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
0985   G4complex A  = std::pow(A1,1./3.);
0986 
0987   G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
0988   G4complex B  = std::pow(B1,1./3.);
0989 
0990   // roots of the incomplete 3rd equation
0991 
0992   G4complex y1 = A + B;
0993   G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
0994   G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
0995  
0996   G4complex x1 = y1 - b/a/3.;
0997   G4complex x2 = y2 - b/a/3.;
0998   G4complex x3 = y3 - b/a/3.;
0999 
1000   G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = "<<real(x2)<<"; re_x3 = "<<real(x3)<<G4endl;
1001   G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = "<<imag(x2)<<"; im_x3 = "<<imag(x3)<<G4endl;
1002 
1003   G4double r1 = real(x1)*real(x1);
1004   G4double r2 = real(x2)*real(x2);
1005   G4double r3 = real(x3)*real(x3);
1006 
1007   if( r1 <= 1. && r1 >= 0. )      fBQ = r1;
1008   else if( r2 <= 1. && r2 >= 0. ) fBQ = r2;
1009   else if( r3 <= 1. && r3 >= 0. )      fBQ = r3;
1010   else fBQ = 1.;
1011   // fBQ = real(x3)*real(x3);
1012   G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
1013   fOptRatio  = fBq + fBQ + 2.*sqrtBqBQ - 1.;
1014   fOptRatio /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
1015   G4cout<<"F123, fOptRatio = "<<fOptRatio<<G4endl;
1016   
1017   return ;
1018 }
1019 
1020 /////////////////////////////////////////////////////
1021 //
1022 //
1023 
1024 inline G4complex G4hhElastic::GetF1(G4double t)
1025 {
1026   G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1027   G4double  k    = p/CLHEP::hbarc;
1028   G4complex exp1 = fBq*std::exp(-GetAqq()*t);
1029   G4complex exp2 = fBQ*std::exp(-GetAQQ()*t);
1030   G4complex exp3 = 2.*std::sqrt(fBq*fBQ)*std::exp(-GetAqQ()*t);
1031 
1032   G4complex res  = exp1 + exp2 + exp3;
1033   res *=  0.25*k*fSigmaTot/CLHEP::pi;
1034   res *= G4complex(0.,1.);
1035   return res;
1036 }
1037 
1038 //////////////////////////////////////////////////////////////
1039 //
1040 // dsigma/dt(s,t) F1
1041 
1042 inline  G4double G4hhElastic::GetdsdtF1(G4double spp, G4double t)
1043 {
1044   fSpp = spp;
1045   G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1046   G4complex F1 = GetF1(t);
1047 
1048   G4double dsdt = CLHEP::pi/p/p;
1049   dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
1050   return dsdt;
1051 }
1052 
1053 /////////////////////////////////////////////////////
1054 //
1055 //
1056 
1057 inline G4complex G4hhElastic::GetF2(G4double t)
1058 {
1059   G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1060   G4double  k   = p/CLHEP::hbarc;
1061   G4complex z1  = GetAqq()*GetAQQ() - 16.*fLambda*fLambda/81.;
1062             z1 /= 2.*(GetAqQ() + 4.*fLambda/9.);
1063   G4complex exp1 = std::exp(-z1*t);
1064 
1065   G4complex z2 = 0.5*( GetAqQ() - 4.*fLambda/9.);
1066 
1067   G4complex exp2 = std::exp(-z2*t);
1068 
1069   G4complex res  = exp1 + exp2;
1070 
1071   G4complex z3 = GetAqQ() + 4.*fLambda/9.;
1072 
1073   res *=  0.25*k/CLHEP::pi;
1074   res *= G4complex(0.,1.);
1075   res /= z3;
1076   res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1077 
1078   return res;
1079 }
1080 
1081 //////////////////////////////////////////////////////////////
1082 //
1083 // dsigma/dt(s,t) F12
1084 
1085 inline  G4double G4hhElastic::GetdsdtF12(G4double spp, G4double t)
1086 {
1087   fSpp = spp;
1088   G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1089   G4complex F1 = GetF1(t) - GetF2(t);
1090 
1091   G4double dsdt = CLHEP::pi/p/p;
1092   dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);  
1093   return dsdt;
1094 }
1095 
1096 /////////////////////////////////////////////////////
1097 //
1098 //
1099 
1100 inline G4complex G4hhElastic::GetF3(G4double t)
1101 {
1102   G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1103   G4double  k   = p/CLHEP::hbarc;
1104   G4complex z1  = GetAqq()*GetAqQ() - 4.*fLambda*fLambda/81.;
1105             z1 /= GetAqq() + GetAqQ() - 4.*fLambda/9.;
1106 
1107   G4complex exp1 = std::exp(-z1*t)*fBq/(GetAqq() + GetAqQ() - 4.*fLambda/9.);
1108 
1109   G4complex z2 = GetAqQ()*GetAQQ() - 1.*fLambda*fLambda/81.;
1110             z2 /= GetAQQ() + GetAqQ() + 2.*fLambda/9.;
1111 
1112   G4complex exp2 = std::exp(-z2*t)*fBQ/(GetAQQ() + GetAqQ() + 2.*fLambda/9.);
1113 
1114   G4complex res  = exp1 + exp2;
1115 
1116 
1117   res *=  0.25*k/CLHEP::pi;
1118   res *= G4complex(0.,1.);
1119   res *= std::sqrt(fBq*fBQ)*fSigmaTot*fSigmaTot/(4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1120 
1121   return res;
1122 }
1123 
1124 //////////////////////////////////////////////////////////////
1125 //
1126 // dsigma/dt(s,t) F123, sampling ds/dt
1127 
1128 inline  G4double G4hhElastic::GetdsdtF123(G4double t)
1129 {
1130   G4double  p   = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1131   G4complex F1  = GetF1(t);
1132             F1 -= fCofF2*GetF2(t);
1133             F1 -= fCofF3*GetF3(t);
1134   G4double dsdt = CLHEP::pi/p/p;
1135   dsdt         *= real(F1)*real(F1) + imag(F1)*imag(F1);  
1136   return dsdt;
1137 }
1138 
1139 //////////////////////////////////////////////////////////////
1140 //
1141 // dsigma/dt(s,t) F123
1142 
1143 inline  G4double G4hhElastic::GetExpRatioF123(G4double spp, G4double t)
1144 {
1145   fSpp = spp;
1146   G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1147 
1148   // qQ-ds/dt
1149 
1150   G4complex F1 = GetF1(t) - fCofF2*GetF2(t) - fCofF3*GetF3(t);
1151 
1152   G4double dsdt = CLHEP::pi/p/p;
1153   dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1); 
1154 
1155   // exponent ds/dt
1156 
1157   G4complex F10 = GetF1(0.) - fCofF2*GetF2(0.) - fCofF3*GetF3(0.);
1158 
1159   fRhoReIm = real(F10)/imag(F10);
1160 
1161   G4double dsdt0 = CLHEP::pi/p/p;
1162   dsdt0 *= real(F10)*real(F10) + imag(F10)*imag(F10); 
1163 
1164   dsdt0 *= G4Exp(-fExpSlope*t);
1165 
1166   G4double ratio = dsdt/dsdt0;
1167  
1168   return ratio;
1169 }
1170 
1171 
1172 //
1173 //
1174 ////////////////////////////////////////////////////////////////////////
1175 
1176 #endif
1177 
1178 
1179 
1180 
1181 
1182
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!