ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hhElastic.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hhElastic.hh
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 //
29 // G4 Model: hadron diffraction elastic scattering with 4-momentum balance
30 //
31 // Class Description
32 // Final state production model for hadron-hadron elastic scattering
33 // in the framework of quark-diquark model with springy Pomeron.
34 // Projectiles are proton, neutron, pions, kaons.
35 // Targets are proton (and neutron).
36 // Class Description - End
37 //
38 // 02.05.14 V. Grichine - 1-st implementation
39 // 10.10.14 V. Grichine - change to combine with low mass diffraction
40 
41 #ifndef G4hhElastic_h
42 #define G4hhElastic_h 1
43 
44 #include "globals.hh"
45 #include <complex>
46 #include "G4Integrator.hh"
47 
48 #include "G4HadronElastic.hh"
49 #include "G4HadProjectile.hh"
50 #include "G4Nucleus.hh"
51 #include "G4HadronNucleonXsc.hh"
52 
53 #include "G4Exp.hh"
54 #include "G4Log.hh"
55 
56 
57 class G4ParticleDefinition;
58 class G4PhysicsTable;
59 class G4PhysicsLogVector;
60 
61 class G4hhElastic : public G4HadronElastic
62 {
63 public:
64  // PL constructor
65  G4hhElastic();
66  // test constructor
67  G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile,
68  G4double plab );
69  // constructor used for low mass diffraction
70  G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile);
71 
72  virtual ~G4hhElastic();
73 
74  virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/,
75  G4Nucleus & /*targetNucleus*/);
76 
77 
78  void Initialise();
79 
80  void BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile); // , G4double plab );
81  void BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile,
82  G4double plab );
83 
84  G4double SampleInvariantT( const G4ParticleDefinition* p, G4double plab, G4int, G4int);
85  G4double SampleBisectionalT( const G4ParticleDefinition* p, G4double plab);
86 
87  G4double SampleTest(G4double tMin ); // const G4ParticleDefinition* p, );
88 
89  G4double GetTransfer( G4int iMomentum, G4int iTransfer, G4double position );
90 
91 private:
92 
93 
94  G4ParticleDefinition* fTarget;
95  G4ParticleDefinition* fProjectile;
96  G4ParticleDefinition* theProton;
97  G4ParticleDefinition* theNeutron;
98  G4ParticleDefinition* thePionPlus;
99  G4ParticleDefinition* thePionMinus;
100 
101 
102  G4double lowEnergyRecoilLimit;
103  G4double lowEnergyLimitHE;
104  G4double lowEnergyLimitQ;
105  G4double lowestEnergyLimit;
106  G4double plabLowLimit;
107 
108  G4int fEnergyBin;
109  G4int fBinT;
110 
111  G4PhysicsLogVector* fEnergyVector;
112  G4PhysicsTable* fTableT;
113  std::vector<G4PhysicsTable*> fBankT;
114 
115 
116  // Gauss model parameters
117 
118 
119  G4double fMff2;
120  G4double fMQ;
121  G4double fMq;
122 
123  G4double fMassTarg; // ~ A
124  G4double fMassProj; // ~ B
125  G4double fMassSum2;
126  G4double fMassDif2;
127 
128 
129  G4double fRA; // hadron A
130  G4double fRQ;
131  G4double fRq;
132  G4double fAlpha;
133  G4double fBeta;
134 
135  G4double fRB; // hadron B
136  G4double fRG;
137  G4double fRg;
138  G4double fGamma;
139  G4double fDelta;
140 
141  G4double fAlphaP;
142  G4double fLambdaFF;
143  G4double fLambda;
144  G4double fEta;
145  G4double fImCof;
146  G4double fCofF2;
147  G4double fCofF3;
148  G4double fRhoReIm;
149  G4double fExpSlope;
150  G4double fSo;
151 
152  G4double fSigmaTot;
153  G4double fBq;
154  G4double fBQ;
155  G4double fBqQ;
156  G4double fOptRatio;
157  G4double fSpp;
158  G4double fPcms;
159  G4double fQcof; // q prime when integrate
160 
161 public: // Gauss model methods
162 
163  void SetParameters();
164  void SetSigmaTot(G4double stot){fSigmaTot = stot;};
165  void SetSpp(G4double spp){fSpp = spp;};
166  G4double GetSpp(){return fSpp;};
167  void SetParametersCMS(G4double plab);
168 
169  G4double GetBq(){ return fBq;};
170  G4double GetBQ(){ return fBQ;};
171  G4double GetBqQ(){ return fBqQ;};
172  void SetBq(G4double b){fBq = b;};
173  void SetBQ(G4double b){fBQ = b;};
174  void SetBqQ(G4double b){fBqQ = b;};
175  G4double GetRhoReIm(){ return fRhoReIm;};
176 
177  void CalculateBQ(G4double b);
178  void CalculateBqQ13(G4double b);
179  void CalculateBqQ12(G4double b);
180  void CalculateBqQ123(G4double b);
181  void SetRA(G4double rn, G4double pq, G4double pQ);
182  void SetRB(G4double rn, G4double pq, G4double pQ);
183 
184  void SetAlphaP(G4double a){fAlphaP = a;};
185  void SetImCof(G4double a){fImCof = a;};
186  G4double GetImCof(){return fImCof;};
187  void SetLambda(G4double L){fLambda = L;};
188  void SetEta(G4double E){fEta = E;};
189  void SetCofF2(G4double f){fCofF2 = f;};
190  void SetCofF3(G4double f){fCofF3 = f;};
191  G4double GetCofF2(){return fCofF2;};
192  G4double GetCofF3(){return fCofF3;};
193 
194  G4double GetRA(){ return fRA;};
195  G4double GetRq(){ return fRq;};
196  G4double GetRQ(){ return fRQ;};
197 
198  G4double GetRB(){ return fRB;};
199  G4double GetRg(){ return fRg;};
200  G4double GetRG(){ return fRG;};
201 
202  // FqQgG stuff
203 
204  G4complex Pomeron();
205 
206  G4complex Phi13();
207  G4complex Phi14();
208  G4complex Phi23();
209  G4complex Phi24();
210 
211  G4complex GetF1qQgG(G4double qp);
212  G4double GetdsdtF1qQgG(G4double s, G4double q);
213  G4complex GetF2qQgG(G4double qp);
214  G4double GetdsdtF12qQgG(G4double s, G4double q);
215  G4complex GetF3qQgG(G4double qp);
216  G4double GetdsdtF123qQgG(G4double q); // sampling ds/dt
217  G4double GetdsdtF13qQG(G4double s, G4double q);
218 
219 
220  // F123 stuff
221 
222  G4complex GetAqq();
223  G4complex GetAQQ();
224  G4complex GetAqQ();
225 
226  G4double GetCofS1();
227  G4double GetCofS2();
228  G4double GetCofS3();
229  G4double GetOpticalRatio();
230 
231  G4complex GetF1(G4double qp);
232  G4double GetdsdtF1(G4double s, G4double q);
233  G4complex GetF2(G4double qp);
234  G4double GetdsdtF12(G4double s, G4double q);
235  G4complex GetF3(G4double qp);
236  G4double GetdsdtF123(G4double q); // sampling ds/dt
237  G4double GetExpRatioF123(G4double s, G4double q);
238 
239  // parameter arrays
240 
241 private:
242 
243  G4int fInTkin;
244  G4double fOldTkin;
245  static const G4double theNuclNuclData[19][6];
246  static const G4double thePiKaNuclData[8][6];
247  G4HadronNucleonXsc* fHadrNuclXsc;
248 };
249 
253 
254 
255 
256 inline G4bool G4hhElastic::IsApplicable(const G4HadProjectile & projectile,
257  G4Nucleus & nucleus)
258 {
259  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
260  projectile.GetDefinition() == G4Neutron::Neutron() ||
261  projectile.GetDefinition() == G4PionPlus::PionPlus() ||
262  projectile.GetDefinition() == G4PionMinus::PionMinus() ||
263  projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
264  projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
265 
266  nucleus.GetZ_asInt() < 2 ) return true;
267  else return false;
268 }
269 
270 
271 inline void G4hhElastic::SetParameters()
272 {
273  // masses
274 
275  fMq = 0.36*CLHEP::GeV; // 0.441*GeV; // 0.36*GeV;
276  fMQ = 0.441*CLHEP::GeV;
277  fMff2 = 0.26*CLHEP::GeV*CLHEP::GeV; // 0.25*GeV*GeV; // 0.5*GeV*GeV;
278 
279  fAlpha = 1./3.;
280  fBeta = 1. - fAlpha;
281 
282  fGamma = 1./2.; // 1./3.; //
283  fDelta = 1. - fGamma; // 1./2.;
284 
285  // radii and exp cof
286 
287  fRA = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
288  fRq = 0.173*fRA; // 2.4/GeV;
289  fRQ = 0.316*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
290  fRB = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
291  fRg = 0.173*fRA; // 2.4/GeV;
292  fRG = 0.173*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
293 
294  fAlphaP = 0.15/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
295  fLambda = 0.25*fRA*fRA; // 0.25
296  fEta = 0.25*fRB*fRB; // 0.25
297  fImCof = 6.5;
298  fCofF2 = 1.;
299  fCofF3 = 1.;
300 
301  fBq = 0.02; // 0.21; // 1./3.;
302  fBQ = 1. + fBq - 2*std::sqrt(fBq); // 1 - fBq; // 2./3.;
303  fBqQ = std::sqrt(fBq*fBQ);
304 
305  fLambdaFF = 1.5/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
306  fSo = 1.*CLHEP::GeV*CLHEP::GeV;
307  fQcof = 0.009*CLHEP::GeV;
308  fExpSlope = 19.9/CLHEP::GeV/CLHEP::GeV;
309 }
310 
311 
313 //
314 // Set target and projectile masses and calculate mass sum and difference squared for Pcms
315 
316 inline void G4hhElastic::SetParametersCMS(G4double plab)
317 {
318  G4int i;
319  G4double trMass = 900.*CLHEP::MeV, Tkin;
320  G4double sl, sh, ds, rAl, rAh, drA, rBl, rBh, drB, bql, bqh, dbq, bQl, bQh, dbQ, cIl, cIh, dcI;
321 
322  Tkin = std::sqrt(fMassProj*fMassProj + plab*plab) - fMassProj;
323 
324  G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(fProjectile,
325  G4ParticleMomentum(0.,0.,1.),
326  Tkin);
327  fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, fTarget );
328 
329  delete theDynamicParticle;
330 
331  fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
332  fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
333 
334  G4double sCMS = std::sqrt(fSpp);
335 
336  if( fMassProj > trMass ) // p,n,pb on p
337  {
338  this->SetCofF2(1.);
339  this->SetCofF3(1.);
340  fGamma = 1./3.; // 1./3.; //
341  fDelta = 1. - fGamma; // 1./2.;
342 
343  if( sCMS <= theNuclNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
344  {
345  this->SetRA(theNuclNuclData[0][1]/CLHEP::GeV,0.173,0.316);
346  this->SetRB(theNuclNuclData[0][2]/CLHEP::GeV,0.173,0.316);
347 
348  this->SetBq(theNuclNuclData[0][3]);
349  this->SetBQ(theNuclNuclData[0][4]);
350  this->SetImCof(theNuclNuclData[0][5]);
351 
352  this->SetLambda(0.25*this->GetRA()*this->GetRA());
353  this->SetEta(0.25*this->GetRB()*this->GetRB());
354  }
355  else if( sCMS >= theNuclNuclData[17][0]*CLHEP::GeV ) // high edge, as s=7000 ???
356  {
357  this->SetRA(theNuclNuclData[17][1]/CLHEP::GeV,0.173,0.316);
358  this->SetRB(theNuclNuclData[17][2]/CLHEP::GeV,0.173,0.316);
359 
360  this->SetBq(theNuclNuclData[17][3]);
361  this->SetBQ(theNuclNuclData[17][4]);
362  this->SetImCof(theNuclNuclData[17][5]);
363 
364  this->SetLambda(0.25*this->GetRA()*this->GetRA());
365  this->SetEta(0.25*this->GetRB()*this->GetRB());
366  }
367  else // in approximation between array points
368  {
369  for( i = 0; i < 19; i++ ) if( sCMS <= theNuclNuclData[i][0]*CLHEP::GeV ) break;
370  if( i == 0 ) i++;
371  if( i == 19 ) i--;
372 
373  sl = theNuclNuclData[i-1][0]*CLHEP::GeV;
374  sh = theNuclNuclData[i][0]*CLHEP::GeV;
375  ds = (sCMS - sl)/(sh - sl);
376 
377  rAl = theNuclNuclData[i-1][1]/CLHEP::GeV;
378  rAh = theNuclNuclData[i][1]/CLHEP::GeV;
379  drA = rAh - rAl;
380 
381  rBl = theNuclNuclData[i-1][2]/CLHEP::GeV;
382  rBh = theNuclNuclData[i][2]/CLHEP::GeV;
383  drB = rBh - rBl;
384 
385  bql = theNuclNuclData[i-1][3];
386  bqh = theNuclNuclData[i][3];
387  dbq = bqh - bql;
388 
389  bQl = theNuclNuclData[i-1][4];
390  bQh = theNuclNuclData[i][4];
391  dbQ = bQh - bQl;
392 
393  cIl = theNuclNuclData[i-1][5];
394  cIh = theNuclNuclData[i][5];
395  dcI = cIh - cIl;
396 
397  this->SetRA(rAl+drA*ds,0.173,0.316);
398  this->SetRB(rBl+drB*ds,0.173,0.316);
399 
400  this->SetBq(bql+dbq*ds);
401  this->SetBQ(bQl+dbQ*ds);
402  this->SetImCof(cIl+dcI*ds);
403 
404  this->SetLambda(0.25*this->GetRA()*this->GetRA());
405  this->SetEta(0.25*this->GetRB()*this->GetRB());
406  }
407  }
408  else // pi, K
409  {
410  this->SetCofF2(1.);
411  this->SetCofF3(-1.);
412  fGamma = 1./2.; // 1./3.; //
413  fDelta = 1. - fGamma; // 1./2.;
414 
415  if( sCMS <= thePiKaNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
416  {
417  this->SetRA(thePiKaNuclData[0][1]/CLHEP::GeV,0.173,0.316);
418  this->SetRB(thePiKaNuclData[0][2]/CLHEP::GeV,0.173,0.173);
419 
420  this->SetBq(thePiKaNuclData[0][3]);
421  this->SetBQ(thePiKaNuclData[0][4]);
422  this->SetImCof(thePiKaNuclData[0][5]);
423 
424  this->SetLambda(0.25*this->GetRA()*this->GetRA());
425  this->SetEta(this->GetRB()*this->GetRB()/6.);
426  }
427  else if( sCMS >= thePiKaNuclData[7][0]*CLHEP::GeV ) // high edge, as s=7000 ???
428  {
429  this->SetRA(thePiKaNuclData[7][1]/CLHEP::GeV,0.173,0.316);
430  this->SetRB(thePiKaNuclData[7][2]/CLHEP::GeV,0.173,0.173);
431 
432  this->SetBq(thePiKaNuclData[7][3]);
433  this->SetBQ(thePiKaNuclData[7][4]);
434  this->SetImCof(thePiKaNuclData[7][5]);
435 
436  this->SetLambda(0.25*this->GetRA()*this->GetRA());
437  this->SetEta(this->GetRB()*this->GetRB()/6.);
438  }
439  else // in approximation between array points
440  {
441  for( i = 0; i < 8; i++ ) if( sCMS <= thePiKaNuclData[i][0]*CLHEP::GeV ) break;
442  if( i == 0 ) i++;
443  if( i == 8 ) i--;
444 
445  sl = thePiKaNuclData[i-1][0]*CLHEP::GeV;
446  sh = thePiKaNuclData[i][0]*CLHEP::GeV;
447  ds = (sCMS - sl)/(sh - sl);
448 
449  rAl = thePiKaNuclData[i-1][1]/CLHEP::GeV;
450  rAh = thePiKaNuclData[i][1]/CLHEP::GeV;
451  drA = rAh - rAl;
452 
453  rBl = thePiKaNuclData[i-1][2]/CLHEP::GeV;
454  rBh = thePiKaNuclData[i][2]/CLHEP::GeV;
455  drB = rBh - rBl;
456 
457  bql = thePiKaNuclData[i-1][3];
458  bqh = thePiKaNuclData[i][3];
459  dbq = bqh - bql;
460 
461  bQl = thePiKaNuclData[i-1][4];
462  bQh = thePiKaNuclData[i][4];
463  dbQ = bQh - bQl;
464 
465  cIl = thePiKaNuclData[i-1][5];
466  cIh = thePiKaNuclData[i][5];
467  dcI = cIh - cIl;
468 
469  this->SetRA(rAl+drA*ds,0.173,0.316);
470  this->SetRB(rBl+drB*ds,0.173,0.173);
471 
472  this->SetBq(bql+dbq*ds);
473  this->SetBQ(bQl+dbQ*ds);
474  this->SetImCof(cIl+dcI*ds);
475 
476  this->SetLambda(0.25*this->GetRA()*this->GetRA());
477  this->SetEta(this->GetRB()*this->GetRB()/6.);
478  }
479  }
480  return;
481 }
482 
484 //
485 // RA for qQ
486 
487 inline void G4hhElastic::SetRA(G4double rA, G4double pq, G4double pQ)
488 {
489  fRA = rA;
490  fRq = fRA*pq;
491  fRQ = fRA*pQ;
492 }
493 
495 //
496 // RB for gG
497 
498 inline void G4hhElastic::SetRB(G4double rB, G4double pg, G4double pG)
499 {
500  fRB = rB;
501  fRg = fRB*pg;
502  fRG = fRB*pG;
503 }
504 
506 //
507 // Returns Pomeron parametrization with Im part modified, *= fImCof
508 
509 inline G4complex G4hhElastic::Pomeron()
510 {
511  G4double re, im;
512 
513  re = fAlphaP*G4Log(fSpp/fSo);
514  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
515  return G4complex(re,im);
516 }
517 
519 
520 inline G4complex G4hhElastic::Phi13()
521 {
522  G4double re = (fRq*fRq + fRg*fRg)/16.;
523  G4complex result(re,0.);
524  result += Pomeron();
525  return result;
526 }
527 
529 
530 inline G4complex G4hhElastic::Phi14()
531 {
532  G4double re = (fRq*fRq + fRG*fRG)/16.;
533  G4complex result(re,0.);
534  result += Pomeron();
535  return result;
536 }
537 
539 
540 inline G4complex G4hhElastic::Phi23()
541 {
542  G4double re = (fRQ*fRQ + fRg*fRg)/16.;
543  G4complex result(re,0.);
544  result += Pomeron();
545  return result;
546 }
547 
549 
550 inline G4complex G4hhElastic::Phi24()
551 {
552  G4double re = (fRQ*fRQ + fRG*fRG)/16.;
553  G4complex result(re,0.);
554  result += Pomeron();
555  return result;
556 }
557 
559 //
560 // F1, case qQ-gG
561 
562 inline G4complex G4hhElastic::GetF1qQgG(G4double t)
563 {
564  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
565  G4double k = p/CLHEP::hbarc;
566 
567  G4complex exp13 = fBq*std::exp(-(Phi13() + fBeta*fBeta*fLambda + fDelta*fDelta*fEta)*t);
568  G4complex exp14 = fBq*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
569  G4complex exp23 = fBQ*std::exp(-(Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta)*t);
570  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
571 
572  G4complex res = exp13 + exp14 + exp23 + exp24;
573 
574  res *= 0.25*k*fSigmaTot/CLHEP::pi;
575  res *= G4complex(0.,1.);
576 
577  return res;
578 }
579 
581 //
582 //
583 
584 inline G4double G4hhElastic::GetdsdtF13qQG(G4double spp, G4double t)
585 {
586  fSpp = spp;
587  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
588  G4double k = p/CLHEP::hbarc;
589 
590  G4complex exp14 = fBqQ*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
591  G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
592 
593  G4complex F1 = exp14 + exp24;
594 
595  F1 *= 0.25*k*fSigmaTot/CLHEP::pi;
596  F1 *= G4complex(0.,1.);
597 
598  // 1424
599 
600  G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
601  z1424 /= Phi14() + Phi24() + fLambda;
602  z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
603 
604 
605  G4complex exp1424 = std::exp(-z1424*t);
606  exp1424 /= Phi14() + Phi24() + fLambda;
607 
608  G4complex F3 = fBqQ*fBQ*exp1424;
609 
610 
611  F3 *= 0.25*k/CLHEP::pi;
612  F3 *= G4complex(0.,1.);
613  F3 *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
614 
615  G4complex F13 = F1 - F3;
616 
617  G4double dsdt = CLHEP::pi/p/p;
618  dsdt *= real(F13)*real(F13) + imag(F13)*imag(F13);
619 
620  return dsdt;
621 }
622 
624 //
625 // dsigma/dt(s,t) F1qQgG
626 
627 inline G4double G4hhElastic::GetdsdtF1qQgG(G4double spp, G4double t)
628 {
629  fSpp = spp;
630  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
631 
632  G4complex F1 = GetF1qQgG(t);
633 
634  G4double dsdt = CLHEP::pi/p/p;
635  dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
636  return dsdt;
637 }
638 
640 //
641 //
642 
643 inline G4complex G4hhElastic::GetF2qQgG(G4double t)
644 {
645  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
646  G4double k = p/CLHEP::hbarc;
647 
648  G4complex z1324 = -(Phi24() + fAlpha*fLambda + fGamma*fEta)*(Phi24() + fAlpha*fLambda + fGamma*fEta);
649  z1324 /= Phi13() + Phi24() + fLambda + fEta;
650  z1324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
651 
652  G4complex exp1324 = std::exp(-z1324*t);
653  exp1324 /= Phi13() + Phi24() + fLambda + fEta;
654 
655  G4complex z1423 = -(Phi23() + fAlpha*fLambda + fDelta*fEta)*(Phi24() + fAlpha*fLambda + fDelta*fEta);;
656  z1423 /= Phi14() + Phi23() + fLambda + fEta;
657  z1423 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
658 
659  G4complex exp1423 = std::exp(-z1423*t);
660  exp1423 /= Phi14() + Phi23() + fLambda + fEta;
661 
662  G4complex res = exp1324 + exp1423;
663 
664 
665  res *= 0.25*k/CLHEP::pi;
666  res *= G4complex(0.,1.);
667  res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc); // or 4. ???
668 
669  return res;
670 }
671 
673 //
674 // dsigma/dt(s,t) F12
675 
676 inline G4double G4hhElastic::GetdsdtF12qQgG( G4double spp, G4double t)
677 {
678  fSpp = spp;
679  G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
680 
681  G4complex F12 = GetF1qQgG(t) - GetF2qQgG(t);
682 
683  G4double dsdt = CLHEP::pi/p/p;
684  dsdt *= real(F12)*real(F12) + imag(F12)*imag(F12);
685  return dsdt;
686 }
687 
689 //
690 //
691 
692 inline G4complex G4hhElastic::GetF3qQgG(G4double t)
693 {
694  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
695  G4double k = p/CLHEP::hbarc;
696 
697  // 1314
698 
699  G4complex z1314 = -(Phi14() + fGamma*fEta)*(Phi14() + fGamma*fEta);
700  z1314 /= Phi13() + Phi14() + fEta;
701  z1314 += Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta;
702 
703  G4complex exp1314 = std::exp(-z1314*t);
704  exp1314 /= Phi13() + Phi14() + fEta;
705 
706  // 2324
707 
708  G4complex z2324 = -(Phi24() + fGamma*fEta)*(Phi24() + fGamma*fEta);;
709  z2324 /= Phi24() + Phi23() + fEta;
710  z2324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
711 
712  G4complex exp2324 = std::exp(-z2324*t);
713  exp2324 /= Phi24() + Phi23() + fEta;
714 
715  // 1323
716 
717  G4complex z1323 = -(Phi23() + fAlpha*fLambda)*(Phi23() + fAlpha*fLambda);
718  z1323 /= Phi13() + Phi23() + fLambda;
719  z1323 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
720 
721  G4complex exp1323 = std::exp(-z1323*t);
722  exp1323 /= Phi13() + Phi23() + fLambda;
723 
724  // 1424
725 
726  G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
727  z1424 /= Phi14() + Phi24() + fLambda;
728  z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
729 
730  G4complex exp1424 = std::exp(-z1424*t);
731  exp1424 /= Phi14() + Phi24() + fLambda;
732 
733  G4complex res = fBq*fBq*exp1314 + fBQ*fBQ*exp2324 + fBq*fBQ*exp1323 + fBq*fBQ*exp1424;
734 
735  res *= 0.25*k/CLHEP::pi;
736  res *= G4complex(0.,1.);
737  res *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
738 
739  return res;
740 }
741 
743 //
744 // dsigma/dt(s,t) F123 sampling ds/dt
745 
746 inline G4double G4hhElastic::GetdsdtF123qQgG(G4double t)
747 {
748  G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp );
749 
750  G4complex F123 = GetF1qQgG(t); // - fCofF2*GetF2qQgG(t) - fCofF3*GetF3qQgG(t);
751  F123 -= fCofF2*GetF2qQgG(t);
752  F123 -= fCofF3*GetF3qQgG(t);
753 
754  G4double dsdt = CLHEP::pi/p/p;
755  dsdt *= real(F123)*real(F123) + imag(F123)*imag(F123);
756  return dsdt;
757 }
758 
760 //
761 // Set fBqQ at a given fBQ=b2 according to the optical theorem,qQ-G
762 
763 inline void G4hhElastic::CalculateBqQ13(G4double b2)
764 {
765  fBQ = b2;
766 
767  G4complex z1424 = G4complex(1./8./CLHEP::pi,0.);
768  z1424 /= Phi14() + Phi24() + fAlpha;
769  G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
770 
771  fBqQ = 1. - fBQ;
772  fBQ /= 1. - fSigmaTot*fBQ*c1424;
773 
774  G4cout<<"fSigmaTot*fBQ*c1424 = "<<fSigmaTot*fBQ*c1424<<G4endl;
775 
776  G4double ratio = fBqQ + fBQ - fSigmaTot*fBqQ*fBQ*c1424;
777  G4cout<<"ratio = "<<ratio<<G4endl;
778 
779  return ;
780 }
781 
783 //
784 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2
785 
786 inline void G4hhElastic::CalculateBqQ12(G4double b1)
787 {
788  fBq = b1;
789 
790  G4complex z1324 = G4complex(1./8./CLHEP::pi,0.);
791  z1324 /= Phi13() + Phi24() + fLambda + fEta;
792  G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
793 
794  G4complex z1423 = G4complex(1./8./CLHEP::pi,0.);
795  z1423 /= Phi14() + Phi23() + fLambda + fEta;
796  G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
797 
798  fBQ = 1. - 2.*fBq;
799  fBQ /= 2. - fSigmaTot*fBq*(c1324+1423);
800 
801  G4double ratio = 2.*(fBq + fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423);
802  G4cout<<"ratio = "<<ratio<<G4endl;
803 
804  return ;
805 }
806 
808 //
809 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2-F3,
810 // simplified meson-barion case g=G=q
811 
812 inline void G4hhElastic::CalculateBqQ123(G4double b1)
813 {
814  fBq = b1;
815 
816  G4complex z1324 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
817  z1324 /= Phi13() + Phi24() + fLambda + fEta;
818  G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
819 
820  G4complex z1423 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
821  z1423 /= Phi14() + Phi23() + fLambda + fEta;
822  G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
823 
824  G4complex z1314 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
825  z1314 /= Phi13() + Phi14() + fEta;
826  G4double c1314 = real(z1314)/(CLHEP::hbarc*CLHEP::hbarc);
827 
828  G4complex z2324 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
829  z2324 /= Phi23() + Phi24() + fEta;
830  G4double c2324 = real(z2324)/(CLHEP::hbarc*CLHEP::hbarc);
831 
832  G4complex z1323 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
833  z1323 /= Phi13() + Phi23() + fLambda;
834  G4double c1323 = real(z1323)/(CLHEP::hbarc*CLHEP::hbarc);
835 
836  G4complex z1424 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
837  z1424 /= Phi14() + Phi24() + fLambda;
838  G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
839 
840  G4double A = fSigmaTot*c2324;
841  G4double B = fSigmaTot*fBq*(c1324 + c1423 + c1323 + c1424) - 2.;
842  G4double C = 1. + fSigmaTot*fBq*fBq*c1314 - 2*fBq;
843  G4cout<<"A = "<<A<<"; B = "<<B<<"; C = "<<C<<G4endl;
844  G4cout<<"determinant = "<<B*B-4.*A*C<<G4endl;
845 
846  G4double x1 = ( -B - std::sqrt(B*B-4.*A*C) )/2./A;
847  G4double x2 = ( -B + std::sqrt(B*B-4.*A*C) )/2./A;
848  G4cout<<"x1 = "<<x1<<"; x2 = "<<x2<<G4endl;
849 
850  if( B*B-4.*A*C < 1.e-6 ) fBQ = std::abs(-B/2./A);
851  else if ( B < 0.) fBQ = std::abs( ( -B - std::sqrt(B*B-4.*A*C) )/2./A);
852  else fBQ = std::abs( ( -B + std::sqrt(B*B-4.*A*C) )/2./A);
853 
854  fOptRatio = 2*(fBq+fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423 + c1323 + c1424);
855  fOptRatio -= fSigmaTot*fBq*fBq*c1314 + fSigmaTot*c2324*fBQ*fBQ;
856  G4cout<<"BqQ123, fOptRatio = "<<fOptRatio<<G4endl;
857 
858  return ;
859 }
860 
861 
862 
866 
867 inline G4complex G4hhElastic::GetAqq()
868 {
869  G4double re, im;
870 
871  re = fRq*fRq/8. + fAlphaP*G4Log(fSpp/fSo) + 8.*fLambda/9.;
872  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
873  return G4complex(re,im);
874 }
875 
877 //
878 //
879 
880 inline G4complex G4hhElastic::GetAQQ()
881 {
882  G4double re, im;
883 
884  re = fRQ*fRQ/8. + fAlphaP*G4Log(fSpp/fSo) + 2.*fLambda/9.;
885  im = -0.5*fAlphaP*fImCof*CLHEP::pi;
886  return G4complex(re,im);
887 }
888 
890 //
891 //
892 
893 inline G4complex G4hhElastic::GetAqQ()
894 {
895  G4complex z = 0.5*( GetAqq() + GetAQQ() );
896  return z;
897 }
898 
900 //
901 //
902 
903 inline G4double G4hhElastic::GetCofS1()
904 {
905  G4complex z = 1./( GetAqQ() + 4.*fLambda/9. );
906  G4double result = real(z);
907  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
908  result *= fSigmaTot*fCofF2;
909  return result;
910 }
911 
913 //
914 //
915 
916 inline G4double G4hhElastic::GetCofS2()
917 {
918  G4complex z = 1./( GetAqq() + GetAqQ() - 4.*fLambda/9. );
919  G4double result = real(z);
920  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
921  result *= fSigmaTot*fCofF3;
922  return result;
923 }
924 
926 //
927 //
928 
929 inline G4double G4hhElastic::GetCofS3()
930 {
931  G4complex z = 1./( GetAQQ() + GetAqQ() + 2.*fLambda/9. );
932  G4double result = real(z);
933  result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
934  result *= fSigmaTot*fCofF3;
935  return result;
936 }
937 
939 //
940 //
941 
942 inline G4double G4hhElastic::GetOpticalRatio()
943 {
944  return fOptRatio;
945  // G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
946  // G4double result = fBq + fBQ + 2.*sqrtBqBQ - 1.;
947  // result /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
948  // return result;
949 }
950 
951 
952 
954 //
955 // Set fBQ at a given fBq=b according to the optical theorem
956 
957 inline void G4hhElastic::CalculateBQ(G4double b1)
958 {
959  fBq = b1;
960  G4double s1 = GetCofS1();
961  G4double s2 = GetCofS2();
962  G4double s3 = GetCofS3();
963  G4double sqrtBq = std::sqrt(fBq);
964 
965  // cofs of the fBQ 3rd equation
966 
967  G4double a = s3*sqrtBq;
968  G4double b = s1*fBq - 1.;
969  G4double c = (s2*fBq - 2.)*sqrtBq;
970  G4double d = 1. - fBq;
971 
972  // cofs of the incomplete 3rd equation
973 
974  G4double p = c/a;
975  p -= b*b/a/a/3.;
976  G4double q = d/a;
977  q -= b*c/a/a/3.;
978  q += 2*b*b*b/a/a/a/27.;
979 
980  // cofs for the incomplete colutions
981 
982  G4double D = p*p*p/3./3./3.;
983  D += q*q/2./2.;
984  G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
985  G4complex A = std::pow(A1,1./3.);
986 
987  G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
988  G4complex B = std::pow(B1,1./3.);
989 
990  // roots of the incomplete 3rd equation
991 
992  G4complex y1 = A + B;
993  G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
994  G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
995 
996  G4complex x1 = y1 - b/a/3.;
997  G4complex x2 = y2 - b/a/3.;
998  G4complex x3 = y3 - b/a/3.;
999 
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 
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 
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 
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 
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 
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 
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 
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 //
1175 
1176 #endif
1177 
1178 
1179 
1180 
1181 
1182