source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4GlauberCrossSections.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
File size: 54.8 KB
Line 
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// G4 Tools program: Glauber Hadron-Nuclear Cross Sections
28// .....................................................
29// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Jan-2005
30//
31//=====================================================================
32
33//#define debug
34//#define bestest
35//#define integrcs
36#define eldiffcs
37
38#include "globals.hh"
39#include <iostream>
40#include <fstream>
41#include <vector>
42#include "G4ios.hh"
43#include "G4QBesIKJY.hh"
44#include "G4QPDGCode.hh"
45#include "G4NucleiPropertiesTable.hh"
46#include "G4NucleiProperties.hh"
47//#include <CLHEP/GenericFunctions/Bessel.hh>
48
49// Integrated cross-sections
50G4double TotalCrSec, InelCrSec, ProdCrSec, ElasticCrSec, QuasyElasticCrSec; 
51// Parameters of hadron-nucleon interaction
52G4double HadrTot, HadrSlope, HadrReIm, DDSect2, DDSect3, Tot00;
53// CHIPS parameters of hadron-proton & hadron-neutron interaction (Re/Im is common)
54G4double HadrTElN, HadrTElP, HadrSlpN, HadrSlpP, HadrPexN, HadrPexP;
55
56// CHIPS Parameters (only for protons, nn is necessary for neutrons np+nn)
57void CHIPSParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV
58{
59  G4double p=HadrMoment/GeV;                   // Momentum in GeV
60  G4double p2=p*p;                             // Squared momentum in GeV^2
61  G4double sp=std::sqrt(p);
62  G4double p2s=p2*sp;
63  G4double p3=p2*p;
64  G4double p4=p2*p2;
65  G4double p5=p4*p;
66  G4double lp=std::log(p);
67  G4double d2=(lp-5.)*(lp-5.);
68  G4double p8=p4*p4;
69  if (iPDG==2212)
70  {
71    HadrTElN=12./(p2s+.05*p+.0001/sp)+.35/p+(6.75+.14*d2+19./p)/(1.+.6/p4); // pn (mb)
72    HadrTElP=2.91/(p2s+.0024)+(6.75+.14*d2+23./p)/(1.+.4/p4);               // pp (mb)
73    HadrSlpN=(7.2+4.32/(p8+.012*p3))/(1.+2.5/p4);                           // pn (GeV^-2)
74    HadrSlpP=8.*std::pow(p,.055)/(1.+3.64/p4);                              // pp (GeV^-2)
75    HadrPexN=(6.75+.14*d2+13./p)/(1.+.14/p4)+.6/(p4+.00013);                // pn (mb/sr)
76    HadrPexP=(74.+3.*d2)/(1.+3.4/p5)+(.2/p2+17.*p)/(p4+.001*sp);            // pn (mb/sr)
77    HadrReIm=-.54+.0954*lp+.975/(p2+.09865/p);           // no unit (the same for pp & pn)
78    //G4cout<<"CHIPSPar:p="<<p<<",N="<<HadrTotN<<",P="<<HadrTotN<<",RI="<<HadrReIm<<G4endl;
79  }     
80  else G4cout<<"CHIPSPar: Only proton projectile (2212) is implemented PDG="<<iPDG<<G4endl;
81}
82
83// Dubna Parameters
84G4int CalculateParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV
85{
86  static const G4double mN=.938;   // Dubna's mass of the nucleon in nuclei in GeV
87  //static const G4double mN=.9315;  // AtomicMassUnit (mass of the nucleon in nuclei) in GeV
88  static const G4double dmN=mN+mN; // Doubled mass of the nucleon in nuclei in GeV
89  static const G4double mN2=mN*mN; // Squared mass of the nucleon in nuclei in GeV^2
90  G4int iHadron=-1;
91  if(iPDG==2212||iPDG==2112||iPDG==3122||iPDG==3222||iPDG==3112||
92     iPDG==3212||iPDG==3312||iPDG==3322||iPDG==3334) iHadron = 0;
93  else if(iPDG==-2212||iPDG==-2112||iPDG==-3122||iPDG==-3222||
94          iPDG==-3112||iPDG==-3212||iPDG==-3312||iPDG==-3322||
95          iPDG==-3334) iHadron = 1;       // Anti-nucleons, Anti-hyperons
96  else if(iPDG== 211)  iHadron = 2;       // Pi+
97  else if(iPDG==-211)  iHadron = 3;       // Pi-
98  else if(iPDG== 321)  iHadron = 4;       // K+
99  else if(iPDG==-321)  iHadron = 5;       // K- @@ What about K0?
100  else
101  {
102    G4cout<<"CalculateParameters: PDG= "<<iPDG<<" is not supported"<<G4endl;
103    return iHadron;
104  }
105  G4double mHadr      = G4QPDGCode(iPDG).GetMass()/1000.; // Hadron Mass in GeV
106  G4double mHadr2     = mHadr*mHadr;               // Squared Hadron Mass in GeV
107           HadrMoment/= 1000.;                     // Momentum in GeV (input parameter)
108  G4double HadrMoment2= HadrMoment*HadrMoment;     // Squared Momentum in GeV^2
109  G4double HadrEnergy2= mHadr2+HadrMoment2;        // Tot energy in GeV
110  G4double HadrEnergy = std::sqrt(HadrEnergy2);    // Tot energy in GeV (global parameter)
111  G4double sHadr      = HadrEnergy*dmN+mN2+mHadr2; // s in GeV^2
112  G4double sqrS       = std::sqrt(sHadr);          // W in GeV
113  //G4cout<<"GHAD: E="<<HadrEnergy<<",dM="<<dmN<<",sN="<<mN2<<",sH="<<mHadr2<<G4endl;
114  switch (iHadron)
115  {
116  case 0:                                      // =====>>> proton or hyperons (the same?)
117    //G4double Delta=1;                                                   // LowEnergy corr
118    //if(HadrEnergy<40.) Delta = 0.916+0.0021*HadrEnergy;
119    HadrTot = 5.2+5.2*std::log(HadrEnergy)+51*std::pow(HadrEnergy,-0.35); // mb
120    HadrSlope = 6.44+0.88*std::log(sHadr)-1;                              // GeV^-2
121    HadrReIm  = 0.13*std::log(sHadr/350)*std::pow(sHadr,-0.18);           // no unit
122    //G4cout<<"GHAD:S="<<sHadr<<",T="<<HadrTot<<",R="<<HadrReIm<<",B="<<HadrSlope<<G4endl;
123    DDSect2   = 11.;                                                      // mb*GeV^-2
124    DDSect3   = 3.;                                                       // mb*GeV^-2
125    // Hyperon corrections
126    if(iPDG==3122||iPDG==3222||iPDG==3112||iPDG==3212) //Lambda,Sig0+-
127    {
128      HadrTot   *= 0.80;
129      HadrSlope *= 0.85;
130    }
131    else if(iPDG==3312||iPDG==3322)            // ---> Xi0-
132    {
133      HadrTot   *= 0.70;
134      HadrSlope *= 0.75;
135    }           
136    else if(iPDG==3334)                        // ---> Omega-
137    {
138      HadrTot   *= 0.60;
139      HadrSlope *= 0.65;
140    }
141    break;
142  case 1:                                      // =====>>> antiproton and antihyperons
143    HadrTot = 5.2+5.2*std::log(HadrEnergy)+123.2*std::pow(HadrEnergy,-0.5); //  mb
144    HadrSlope = 8.32+0.57*std::log(sHadr);     // GeV^-2
145    if(HadrEnergy<1000.) HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
146    else HadrReIm  = 0.6*std::log(sHadr/350)*std::pow(sHadr,-0.25);
147    DDSect2 = 11.;                             //mb*GeV-2 @@ the same as in case 0
148    DDSect3 = 3.;                              //mb*GeV-2 @@ the same as in case 0
149    // Hyperon corrections
150    if(iPDG==-3122||iPDG==-3222||iPDG==-3112||iPDG==-3212)//AntiLam,Sig
151    {
152      HadrTot   *=0.75;
153      HadrSlope *=0.85;
154    }
155    else if(iPDG==-3312||iPDG==-3322)          // =====>>> AntiXi0-
156    {
157      HadrTot   *=0.65;
158      HadrSlope *=0.75;
159    }
160    if(iPDG==-3334)                            // ======>>>AntiOmega-
161    {
162      HadrTot   *=0.55;
163      HadrSlope *=0.65;
164    }
165    break;
166  case 2:                                      // =====>>>  pi plus
167    if(HadrMoment>2.) HadrTot = 10.6+2.*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.43);// mb
168    else HadrTot = 40-50*(HadrMoment-1.5)*(HadrMoment-1.7);
169    HadrSlope = 7.28+0.245*std::log(sHadr);                                // GeV^-2
170    HadrReIm  = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15);             // no dim
171    DDSect2   = 4.6;                                                       // mb*GeV^-2
172    DDSect3   = 1.33;                                                      // mb*GeV^-2
173    break;
174  case 3:                                      // =====>>>  pi minus
175    HadrTot = 10.6+2*std::log(HadrEnergy)+30*std::pow(HadrEnergy,-0.43);   // mb @@ E->inf?
176    if(HadrMoment<1.399) HadrTot = HadrTot+21.0/0.4*(1.4-HadrMoment);      // A huge jump ?
177    HadrSlope = 7.28+0.245*std::log(sHadr);                                // GeV-2
178    HadrReIm  = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15);
179    DDSect2   = 4.6;                           //mb*GeV-2 @@ the same as in case 2
180    DDSect3   = 1.33;                          //mb*GeV-2
181    break;
182  case 4:                                      // =====>>> K plus
183    HadrTot = 10.6+1.8*std::log(HadrEnergy)+9.*std::pow(HadrEnergy,-0.55); // mb
184    if(HadrEnergy>100) HadrSlope = 15.;              // Jump at 100 GeV (?)
185    else HadrSlope = 1.0+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5);     // GeV-2
186    //else HadrSlope = 5.28+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5);    // GeV-2
187    HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
188    DDSect2   = 3.5;                           //mb*GeV-2
189    DDSect3   = 1.03;                          //mb*GeV-2
190    break;
191  case 5:                                      // ======>>> K minus
192    HadrTot = 10+1.8*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.5);      // mb
193    HadrSlope = 6.98+0.127*std::log(sHadr);           // GeV-2
194    //if(HadrEnergy<8) HadrReIm = 0.7;
195    HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
196    DDSect2   = 3.5;                           //mb*GeV-2
197    DDSect3   = 1.03;                          //mb*GeV-2
198    break;
199  }     
200  return iHadron;
201}
202
203void  CalculateIntegralCS(G4int  Anuc, G4double HadrEnergy)
204{
205  static const G4double eps = 0.0001;          // Accuracy of calculations
206  static const G4double mb2G2 = 2.568;         // Transformation from mb to GeV^-2
207  static const G4double incrf = 10;            // Increase factor of radius for intergation
208  static const G4double m2G10 = incrf*mb2G2;   // Big conversion factor
209  static const G4double dR0   = 1.06+1.06;     // Doubled nominal R0
210  G4double An      = Anuc;                     // A float
211  G4double Stot    = HadrTot*mb2G2;            // Total hadron-nucleon CS in GeV-2
212  G4double Bhad    = HadrSlope;                // Diffraction cone nH In GeV-2
213  G4double Asq     = 1+HadrReIm*HadrReIm;      // |M|^2/(ImM)^2
214
215  G4double R0      = std::sqrt(0.99);          // in fermi (strange value? 1.07?)
216
217  if (Anuc == 58) R0 = std::sqrt(0.6);         // Personal correction for Fe
218  if (Anuc >20)   R0 = std::sqrt((35.34+0.5*Anuc)/(40.97+Anuc));
219  if (Anuc == 16) R0 = std::sqrt(0.75);        // Personal correction for O16
220  if (Anuc >10)   R0 = std::sqrt(0.84);        // The same R0 for A=11-20 (except O16)
221  if (Anuc ==  4) R0 *= 1.2;                   // Personal correction for Helium
222
223  G4double Rnucl   = R0*std::pow(static_cast<double>(Anuc),.3333333); // in fermi
224  G4double Rnuc2   = Rnucl*Rnucl*m2G10;        // in GeV^-2 (10 -> to make it >> R)
225  G4double RB      = Rnuc2+Bhad;               // Effective increase of R2 (Convolution?)
226  G4double R2B     = RB+Bhad;                  // Over-convolution (?)
227  G4double Delta   = Stot/R2B/twopi;           // Integration step for Total & Inelastic CS
228  G4double Delt2   = Delta+Delta;              // Two delta
229  G4double Delt3   = Stot*Asq*R2B/RB/Bhad/16/pi; // Integration step for Production CS
230
231  G4double Tot0    = 0.;                       // =SUM[(-Delta)^(i-1)*(C^A_i)/i]
232  G4double Inel0   = 0.;                       // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i]
233  G4double N       = -1./Delta; 
234  G4double N1      = N;
235  G4double N3      = -1./Delt2;
236  G4double Prod0   = 0.;
237  G4int    Ai      = Anuc+1;                   // Working value to avoid multiplication
238  G4double RBi     = 0.;                       // Working value to avoid multiplication
239  for (G4int i=1; i<= Anuc; i++)               // ==> Loop over nucleons
240  {
241    Ai--;                                      // Working value to avoid multiplication
242    RBi   += RB;                               // Working value to avoid multiplication
243    G4double Di=Delta/i;                       // Working value to avoid repetition
244    N     *= -Di*Ai;                           // (-Delta)^(i-1)*A!/(A-i)!/i! (C^A_i)
245    N1    *= -Di*(Anuc+Ai);                    // (-Delta)^(i-1)*(2A)!/(2A-i)!/i!
246    N3    *= -Delt2*Ai/i;                      // (-2*Delt)^(i-1)*A!/(A-i)!/i! (C^A_i)
247    Tot0  += N/i;                              // =SUM[(-Delta)^(i-1)*(C^A_i)/i]
248    Inel0 += N1/i;                             // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i]
249    //G4double N2    = 1.;
250    G4double N4    = 1.;
251    //G4double Inel1 = 0.;
252    //G4double Inel2 = 0.;
253    G4double Prod1 = 0.;
254    G4double Bhl   = 0.;
255    for (G4int l=0; l<= i; l++)                // ==> Loop over nucleon-nucleon pairs (0?)
256    {
257      //Inel1 += N2/(i+l);                 // SUM[(-Delta)^(l-1)*i!(i-1)!/(l-1)!/l!/(i+l)!]
258      Prod1 += N4*RB/(RBi+Bhl);                //
259      //N2    *= -Delt*(i-l)/(l+1);              // (-Delta)^l*i!/l!/(l+1)!
260      N4    *= -Delt3*(i-l)/(l+1);
261      Bhl   += Bhad;
262    }
263    //Inel2 += Inel1*N3;
264    Prod0   += Prod1*N3;
265    if(std::fabs(N1/i/Inel0) < eps)  break;
266  }
267
268  Tot0             *= HadrTot;                 // Multiply by the hN cross section
269  Inel0            *= HadrTot/2;               // Multiply by a half of hN cross section
270  //Inel2            *= HadrTot;                 // Multiply by the hN cross section
271  Prod0            *= HadrTot;                 // Multiply by the hN cross section
272  Tot00             = Tot0;                    // Remember because it will be changed
273  G4double ak       = Rnuc2*twopi/Stot;        // cross-section ratio for integration
274  G4double Ank      = An/ak;                   // Working value to avoid repetition
275  G4double rak      = 1./ak;                   // Working value to avoid repetition
276  G4double ak1      = 1.-rak/4;                // Working value to avoid repetition
277  G4double Ank1     = Ank*ak1;                 // Working value to avoid repetition
278  G4double DDSect1  = (DDSect2+DDSect3*std::log(dR0*HadrEnergy/Rnucl/std::sqrt(m2G10)/4));
279  G4double Dtot     = 8*pi*ak/HadrTot*(1.-(1.+Ank)*std::exp(-Ank))*DDSect1/mb2G2;//TotCSCor
280  G4double bk       = (1.-rak)/Stot/ak1;       // Working value to avoid power(,2)
281  G4double bd       = bk*bk*DDSect1*(1.-(1.+Ank1)*std::exp(-Ank1))*Rnuc2; // Working
282  G4double Dprod    = bd*4*pi2*mb2G2;          // Correction for the Production & Inel CS
283  TotalCrSec        = Tot0-Dtot;               // Total Cross Section         (primary)
284  InelCrSec         = Inel0-Dprod;             // Inelastic Cross Section     (primary)
285  //InelCrSec1        = Inel2;                 // Direct calculation without correction
286  ProdCrSec         = Prod0-Dprod;             // Production Cross Section    (primary)
287  ElasticCrSec      = TotalCrSec-InelCrSec;    // Elastic Cross Section       (secondary)
288  QuasyElasticCrSec = InelCrSec-ProdCrSec;     // Quasy-Elastic Cross Section (secondary)
289}
290
291#ifdef eldiffcs
292// Functions for calculation of differential elastic cross-sections
293// Static globals
294static const G4double  rAmax = 2.5;      // Factor of maximum radius
295static const G4int     NptB  = 500;      // A#of steps in the impact parameter integration
296static const G4int     nFact = 253;      // Maximum A (for Factorials and Factors)
297
298// Globals
299G4double r0 = 1.1;                       // Wood-Saxon's factor
300G4double Mnoj[nFact], Factorials[nFact]; // Factors and factorials
301G4double R1, R2, Pnucl, Aeff;            // Nuclear Parameters
302G4double InCoh;                          // Incoherent quasi-elastic cross section
303
304G4double binom(G4int N, G4int M)
305{
306  if(N>170 && N>M && M)
307                {
308    G4double fN=N;
309    G4double fM=M;
310    G4double fD=N-M;
311    G4double man=N*std::log(fN)-M*std::log(fM)-fD*std::log(fD)-1.;
312    //G4cout<<"G4GCS::binom: N="<<N<<", M="<<M<<", C="<<man<<", E="<<std::exp(man)<<G4endl;
313    return std::exp(man);
314                }
315  else if(N>1 && N>M && M) return Factorials[N]/Factorials[M]/Factorials[N-M];
316  return 1.;
317}
318
319// Initialize factorials and combinatoric coefficients
320void Initialize()
321{
322  //static const G4double  sat = 3.03;      // Stirling saturation of Dubna
323  G4double Val = Factorials[0] = Mnoj[0] = 1.;
324  for (G4int i=1; i<nFact; i++)
325  {
326    G4double Si = 0.;
327             Val *= i;                    // Calculation of the factorial values
328    Factorials[i] = Val;
329    for (G4int l = 0; l<=i; l++)
330    {
331      G4double Sm = 1.;                   // sum of reversed binom coefficients (?)
332      for (G4int m=1; m<=i-l; m++) Sm += 1./binom(i-l,m);
333      //if(i>169)G4cout<<"G4GCS::Init:S="<<Sm<<", b("<<i<<","<<l<<")="<<binom(i,l)<<G4endl;
334      Si += Sm/binom(i,l);
335    } //  End of l LOOP
336    Mnoj[i] = Si;
337    //G4double d=Si-sat;
338    //G4double r=d/sat;
339    //if(i>100)G4cout<<"G4GCS::Init:i="<<i<<":"<<Si<<"-"<<sat<<" = "<<d<<", r="<<r<<G4endl;
340  } // End of i LOOP
341} // End of Initialization
342
343// Dubna paranmeters of nuclei
344void CalculateNuclearParameters(G4int Anuc)  // R1,R2,Pnuc,Aeff
345{//  ==============================================
346  G4double A=Anuc;
347  // Personal corrections fore some nuclei
348  if(Anuc == 4)                                      // Personal correction for He (!)
349  {
350    R1    = 5.5;
351    R2    = 3.7;
352    Pnucl = 0.4;   
353    Aeff  = 0.87;
354  }
355  else if(Anuc == 9)                                 // Personal correction for Be (!)
356  {
357    R1    = 9.0;
358    R2    = 7.0;
359    Pnucl = 0.190;
360    Aeff  = 0.9;
361  }
362  //else if(Anuc == 11)                              // Personal correction for B (!)
363  //{
364  //  R1    = 10.8;
365  //  R2    = 7.5;
366  //  Pnucl = 0.85;
367  //  Aeff  = 1.2;
368  //}
369                //else if(Anuc == 12)                              // Personal correction for C (!)
370  //{
371  //   R1    = 9.336;
372  //   R2    = 5.63;
373  //   Pnucl = 0.197;
374  //   Aeff  = 1.0;
375  //}
376  else if(Anuc == 16)                                // Personal correction for O (!)
377  {
378    R1    = 10.50;
379    R2    = 5.5;
380    Pnucl = 0.7;
381    Aeff  = 0.98;
382    //R1    = 11.3;
383    //R2    = 2.5;
384    //Pnucl = 0.75;
385    //Aeff  = 0.9;
386  }
387  else if(Anuc == 58)                                // Personal correction for Ni (!)
388  {
389    R1    = 15.0;
390    R2    = 9.9;
391    Pnucl = 0.45;
392    Aeff  = 0.85;
393  }
394  else if(Anuc == 90)                                // Personal correction for Zr (!)
395  {
396    //R1    = 16.5;
397    //R2    = 11.62;
398    //Pnucl = 0.4;
399    //Aeff  = 0.9;
400    R1    = 16.5;
401    R2    = 11.62;
402    Pnucl = 0.4;
403    Aeff  = 0.7;
404  }
405  else if(Anuc == 208)                               // Personal correction for Pb (!)
406  { 
407    //      R1 = 20.73; R2 = 15.74.
408    //R1    = 4.1408*std::pow(A,0.3018);
409    //R2    = 3.806*std::pow(A-10.068,0.2685);
410    //Pnucl = 0.9;
411    //Aeff  = 1.1;
412    R1    = 19.5;
413    R2    = 15.74;
414    Pnucl = 0.4;
415    Aeff  = 0.7;
416  }
417  else
418  {
419    R1 = 4.45*std::pow(A-1.,.309);                  // First diffractional radius
420    if(Anuc == 28) R1 *= .955;                       // Personal correction for Si (?!)
421    R2    = 2.3*std::pow(A,.36);                    // Second diffraction radius
422    Pnucl = 0.176+(.00167+.00000869*Anuc)*Anuc;      // Momentum of mean Fermi motion
423    Aeff  = 0.9;                                     // Effective screaning
424  }
425}
426
427// Independent transition method calculation. Has problems for A>95 (short cut).
428G4double DiffElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) // All in MeV
429{//      =================================================================
430  static const G4double eps = 0.000001;                     // Accuracy of calculations
431  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
432  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
433  static const G4double dR0   = 1.06+1.06;                  // Doubled nominal R0
434  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
435                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
436  G4double MassH2     = MassH*MassH;
437  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
438                //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
439  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
440  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
441  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
442  CalculateParameters(hPDG, HadrMom);                   
443  CalculateIntegralCS(A, HadrEnergy);
444  CalculateNuclearParameters(A);
445  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
446                //MassN           = A*0.938;                              // @@ This is bad!
447                G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
448  G4double MassN2   = MassN*MassN;
449  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
450  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
451  G4double CMMom    = std::sqrt(EcmH*EcmH-MassH2);          // CM momentum
452
453  G4double Stot     = HadrTot*mb2G2;                        // in GeV^-2
454  G4double Bhad     = HadrSlope;                            // in GeV^-2
455  G4double Asq      = 1+HadrReIm*HadrReIm;                  // |M|^2/(ImM)^2
456  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
457  G4double R12      = R1*R1;
458  G4double R22      = R2*R2;
459  G4double R12B     = R12+Bhad+Bhad;
460  G4double R22B     = R22+Bhad+Bhad;
461  G4double R12Ap    = R12+20;
462  G4double R22Ap    = R22+20;
463  G4double R13Ap    = R12*R1/R12Ap;
464  G4double R23Ap    = R22*R2/R22Ap*Pnucl;
465  G4double R23dR13  = R23Ap/R13Ap;
466  G4double R12Apd   = 2./R12Ap;
467  G4double R22Apd   = 2./R22Ap;
468  G4double R122Apd  = (R12Apd+R22Apd)/2;
469  G4double dR0H     = dR0*HadrEnergy/4;
470  G4double DDSec1p  = DDSect2+DDSect3*std::log(dR0H/R1);
471  G4double DDSec2p  = DDSect2+DDSect3*std::log(dR0H/std::sqrt((R12+R22)/2));
472  G4double DDSec3p  = DDSect2+DDSect3*std::log(dR0H/R2);
473  G4double Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;    // Some questionable norming
474  G4double R13      = R12*R1/R12B;
475  G4double R23      = Pnucl*R22*R2/R22B;
476  G4double norFac   = Stot/twopi/Norm;               // totCS (in GeV^-2) is used
477  G4double Unucl    = norFac*R13;
478  G4double UnuclScr = norFac*R13Ap;
479  G4double SinFi    = HadrReIm/Rho2;
480  G4double FiH      = std::asin(SinFi);
481  G4double N        = -1.;
482  G4double N2       = R23/R13;
483  G4double ImElA0   = 0.;
484  G4double ReElA0   = 0.;
485  G4double Tot1     = 0.;
486  for(G4int i=1; i<=A; i++)                          // @@ Make separately for n and p
487  {
488    N              *= -Unucl*(A-i+1)*Rho2/i;         // Includes total cross-section
489    G4double N4     = 1.;
490    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i; // Includes the slope
491    G4double medTot = R12B/i;
492    for(G4int l=1; l<=i; l++)
493    {
494      G4double exp1 = l/R22B+(i-l)/R12B;
495      N4           *= -(i-l+1)*N2/l;
496      G4double Nexp = N4/exp1;
497      Prod1        += Nexp*std::exp(-Q2/exp1/4);
498      medTot       += Nexp;
499    }
500    G4double FiHi   = FiH*i;
501    G4double nCos   = N*std::cos(FiHi);
502    ReElA0 += Prod1*N*std::sin(FiHi);
503    ImElA0 += Prod1*nCos;
504    Tot1           += medTot*nCos;
505    if(std::abs(Prod1*N/ImElA0) < eps) break;
506  }
507  ImElA0   *= piMG;                                  // The amplitude in mb
508  ReElA0   *= piMG;                                  // The amplitude in mb
509  Tot1             *= piMG+piMG;
510  G4double N1p      = 1.;
511  G4double R13Ap1   = DDSec1p*R13Ap*R13Ap/2;
512  G4double R22Ap1   = DDSec2p*R13Ap*R23Ap;
513  G4double R23Ap1   = DDSec3p*R23Ap*R23Ap*R22Ap/2;
514  G4double R13Ap2   = R13Ap1*R12Ap/4;
515  G4double R22Ap2   = R22Ap1/R122Apd/2;
516  G4double R23Ap2   = R23Ap1*R22Ap/4;
517  G4double Q28      = Q2/8;
518  G4double Q24      = Q28+Q28;
519  G4double Din1     = R13Ap2*std::exp(-Q28*R12Ap)-R22Ap2*std::exp(-(Q28+Q28)/R122Apd)+
520                                                  R23Ap2*std::exp(-Q28*R22Ap);   // i=0 start value
521  G4double DTot1    = R13Ap2-R22Ap2+R23Ap2;          // i=0 start value
522  if (A<96)
523                {
524    for(G4int i=1; i<A-1; i++)
525    {
526      N1p              *= -UnuclScr*(A-i-1)/i*Rho2;
527      G4double N2p      = 1.;
528      G4double Din2     = 0.;
529      G4double DmedTot  = 0.;
530      G4double BinCoeff = 1.;                        // Start for Binomial Coefficient
531      for(G4int l = 0; l<=i; l++) 
532      {
533        if(l) BinCoeff *= (i-l+1)/l;
534        G4double exp1   = l/R22B+(i-l)/R12B;
535        G4double exp1p  = exp1+R12Apd;
536        G4double exp2p  = exp1+R122Apd;
537        G4double exp3p  = exp1+R22Apd;
538        G4double NBinC  = N2p*BinCoeff;
539        Din2           += NBinC*(R13Ap1/exp1p*std::exp(-Q24/exp1p)-
540                                 R22Ap1/exp2p*std::exp(-Q24/exp2p)+
541                                 R23Ap1/exp3p*std::exp(-Q24/exp3p));
542        DmedTot        += NBinC * (R13Ap1/exp1p - R22Ap1/exp2p + R23Ap1/exp3p);
543        N2p            *= -R23dR13;
544      }
545      G4double comFac   = N1p*Mnoj[i]/(i+2)/(i+1)*std::cos(FiH*i);
546                                  Din1             += Din2*comFac;
547      DTot1            += DmedTot*comFac;
548      if(std::abs(Din2*N1p/Din1) < eps) break;
549    }
550    G4double cFac       = A*(A-1)*4/Norm/Norm;
551    Din1               *= -cFac;
552    DTot1              *= cFac;
553  }
554  //else Din1=0.;
555
556  G4double Corr0      = Tot00/Tot1;
557  ImElA0             *= Corr0;
558                       
559  //return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; // ds/do
560  return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*mb2G2/(twopi+twopi);
561} // End of DiffElasticCS
562
563G4double CHIPSDiffElCS(G4int hPDG, G4double HadrMom, G4int Z, G4int A, G4double aQ2) // MeV
564{//      =================================================================
565  static const G4double eps = 0.000001;                     // Accuracy of calculations
566  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
567  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
568  G4double p          = HadrMom/1000.;                      // Momentum (GeV, inputParam)
569  G4double p2         = p*p;
570                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
571  G4double MassH2     = MassH*MassH;
572  G4double HadrEnergy = std::sqrt(p2+MassH2);               // Tot energy in GeV
573                //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
574  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
575  CalculateParameters(hPDG, HadrMom);                       // ?                   
576  CalculateIntegralCS(A, HadrEnergy);                       // ?
577  CalculateNuclearParameters(A);                            // ?
578  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
579                //G4doubleMassN     = A*0.938;                              // @@ This is bad!
580                G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
581  if(G4NucleiPropertiesTable::IsInTable(Z,A))
582                                MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.;    // Geant4 NuclearMass in GeV
583  G4double MassN2   = MassN*MassN;
584  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
585  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
586  G4double CMMom    = std::sqrt(EcmH*EcmH-MassH2);          // CM momentum
587  //G4cout<<"P="<<CMMom<<",E="<<HadrEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl;
588  G4double p3       = p2*p;
589  G4double p4       = p2*p2;
590  G4double sp       = std::sqrt(p);
591  G4double p2s      = p2*sp;
592  G4double ap       = std::log(p);
593  G4double dl       = ap-3.;
594  G4double dl2      = dl*dl;
595  G4double shPPTot  = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
596  G4double shPNTot  = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
597                      +(38.+.3*dl2+8./p)/(1.+1.2/p3);       // SigPP in mb
598  G4double hNTot    = (Z*shPPTot+(A-Z)*shPNTot)/A;          // SighN in mb
599  G4double Stot     = hNTot*mb2G2;                          // in GeV^-2
600  G4double shPPSl   = 8.*std::pow(p,.055)/(1.+3.64/p4);     // PPslope in GeV^-2
601  G4double shPNSl   = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
602  G4double Bhad     = (Z*shPPSl+(A-Z)*shPNSl)/A;            // B-slope in GeV^-2
603  G4double hNReIm   = -.55+ap*(.12+ap*.0045);               // Re/Im_hN in no unit
604  G4double Asq      = 1+hNReIm*hNReIm;                      // |M|^2/(ImM)^2
605  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
606  G4double r1       = 3.9*std::pow(A-1.,.309);             // Positive diffractionalRadius
607  G4double r2       = 2.*std::pow(A,.36);                  // Negative diffraction radius
608  G4double pN       = Pnucl;                                // Dubna value
609  //G4double pN       = .4;                                   // Screaning factor
610  G4double Ae       = Aeff;                                 // Dubna value
611  //G4double Ae       = .75;                                   // Normalization
612  G4double R12      = r1*r1;
613  G4double R22      = r2*r2;
614  G4double R1C      = R12*r1;
615  G4double R2C      = R22*r2;
616  G4double R12B     = R12+Bhad+Bhad;                        // Slope is used
617  G4double R22B     = R22+Bhad+Bhad;                        // Slope is used
618  G4double Norm     = (R1C-pN*R2C)*Ae;                      // ScreanFac & NormFac are used
619  G4double R13      = R1C/R12B;                             // Slope is used
620  G4double R23      = pN*R2C/R22B;                          // Slope & ScreanFact are used
621  G4double norFac   = Stot/twopi/Norm;                      // totCS (in GeV^-2) is used
622  G4double Unucl    = norFac*R13;
623  G4double SinFi    = hNReIm/Rho2;                          // Real part
624  G4double FiH      = std::asin(SinFi);
625  G4double N        = -1.;
626  G4double N2       = R23/R13;                              // Slope is used
627  G4double ImElA0   = 0.;
628  G4double ReElA0   = 0.;
629  G4double Tot1     = 0.;
630  for(G4int i=1; i<=A; i++)
631  {
632    N              *= -Unucl*(A-i+1)*Rho2/i;                // TotCS is used
633    G4double N4     = 1.;
634    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i;        // Slope is used
635    G4double medTot = R12B/i;                               // Slope is used
636    for(G4int l=1; l<=i; l++)
637    {
638      G4double exp1 = l/R22B+(i-l)/R12B;                    // Slope is used
639      N4           *= -(i-l+1)*N2/l;                        // Slope is used
640      G4double Nexp = N4/exp1;
641      Prod1        += Nexp*std::exp(-Q2/exp1/4);
642      medTot       += Nexp;
643    }
644    G4double FiHi   = FiH*i;
645    G4double nCos   = N*std::cos(FiHi);
646    ReElA0 += Prod1*N*std::sin(FiHi);
647    ImElA0 += Prod1*nCos;
648    Tot1           += medTot*nCos;
649    if(std::abs(Prod1*N/ImElA0) < eps) break;
650  }
651  ImElA0         *= piMG;                          // The amplitude in mb
652  ReElA0         *= piMG;                          // The amplitude in mb
653  Tot1           *= piMG+piMG;
654  G4double Corr0  = Tot00/Tot1;
655  ImElA0         *= Corr0;
656                       
657  //return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; // ds/do
658  return (ReElA0*ReElA0+ImElA0*ImElA0)*mb2G2/(twopi+twopi);
659} // End of DiffElasticCS
660
661G4double CHIPS_Tb(G4int  A, G4double b)              // T(b) in fm-2
662{
663  static G4int Am=0;                                 // Associative memory value A
664  static G4double B=0.;                              // Effective edge
665  static G4double C=0.;                              // Normalization factor
666  static G4double D=0.;                              // Slope parameter
667  if(A!=Am)
668                {
669    B=.0008*A*A;                                     // no units
670    D=.42*std::pow(A,-.26);                               // fm^-2
671    C=A*D/pi/std::log(1.+B);                         // fm^-2
672  }
673  G4double E=B*std::exp(-D*b*b);
674  return C*E/(1.+E);                                 // T(b) in fm^-2
675}
676
677G4double Thickness(G4int A, G4double b, G4double R)  // T(b) in fm^-2
678{//      ==========================================  // rAmax=2.5 is a GlobStatConstPar
679  static const G4double bTh  = .545;                 // Surface thicknes
680  static const G4double bTh2 = bTh*bTh;              // SquaredSurface thicknes
681  G4double b2    = b*b;                              // Squared impact parameter
682  G4double dr    = rAmax*R/(NptB-1);                 // Step of integration in z
683  G4double R2    = R*R;                              // Working parameter
684  G4double Norm  = .75/pi/R2/R/(1.+bTh2*pi2/R2);     // Corrected Volume (?)
685  G4double r     = -dr;                              // Pre value for the radius
686  G4double SumZ  = 0.;                               // Prototype of the result
687  //G4double SumN  = 0.;                             // is not used
688  for(G4int k=0; k<NptB; k++)
689  {
690    r    += dr;
691    SumZ += 1./(1.+std::exp((std::sqrt(b2+r*r)-R)/bTh)); // Woods-Saxon density
692    //SumN += r*r/(1.+std::exp((r-rAfm)/bTh));       // is not used
693  }
694  //SumN *= (twopi+twopi)*dr*Norm;                   // is not used
695  return SumZ*(A+A)*dr*Norm;                         // T(b) in fm^-2
696} // End of Thickness
697
698G4double CoherentDifElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) //in MeV
699{//      =========================================================================
700  //static const G4double eps = 0.000001;                   // CalculationAccuracy@@NotUsed
701  //static const G4double mN = .9315;                      // Atomic Unit GeV
702  static const G4double mN = .938;                      // Atomic Unit GeV
703  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
704  static const G4double incrf = 10;            // Increase factor of radius for intergation
705  static const G4double m2G10 = incrf*mb2G2;                // Big conversion factor
706  static const G4double s2G10 = std::sqrt(m2G10);           // sqrt Big conversion factor
707  G4double ReIntegrand[NptB],ImIntegrand[NptB],Thick[NptB]; // Calculated arrays
708  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
709  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
710  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
711                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
712  G4double MassH2     = MassH*MassH;
713  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
714  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
715  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
716  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
717  G4double Re         = std::pow(A, 0.33333333);
718  G4double Re2        = Re*Re;
719  CalculateParameters(hPDG, HadrMom);                   
720  G4double Q2 = aQ2/1000/1000;                              //  GeV
721  // Individual corrections for nuclei which had data
722  if    (A==208) r0   = 1.125;
723  else if(A==90) r0   = 1.12;
724  else if(A==64) r0   = 1.1;
725  else if(A==58) r0   = 1.09;
726  else if(A==48) r0   = 1.07;
727  else if(A==40) r0   = 1.15;
728  else if(A==28) r0   = 0.93;
729  else if(A==16) r0   = 0.92;
730  else if(A==12) r0   = 0.80;
731                else           r0   = 1.16*(1.-1.16/Re2);  // For other nuclei which have not been tested
732                G4double MassN      = A*0.9315;                             // @@ Atomic Unit is bad too
733  G4double MassN2     = MassN*MassN;
734  G4double S          = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S
735  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
736  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum
737  G4double rAfm       = r0*Re;                              // Nuclear radius
738  G4double rAGeV      = rAfm*s2G10;                         // Big integration radius
739  G4double stepB      = rAmax*rAGeV/(NptB-1);               // dr step of integration
740  G4double hTotG2     = HadrTot*mb2G2;
741  G4double ReSum      = 0.;
742  G4double ImSum      = 0.;
743  G4double ValB       = -stepB;
744  for(G4int i=0; i<NptB; i++)
745  {
746    ValB             += stepB;                              // An incident parameter
747    G4double ValB2    = ValB*ValB;                          // A working value
748    G4double IPH      = ValB/HadrSlope;                     // A working value
749    G4double dHS      = HadrSlope+HadrSlope;                // A working value
750    G4double Integ    = 0.;
751    G4double ValS     = 0.;
752    for(G4int j=1; j<NptB; j++)
753    {
754      ValS           += stepB;                              // Impact parameter GeV^-1
755      if(!i) Thick[j] = Thickness(A,ValS/s2G10,rAfm)/m2G10; // Calculate only once
756      G4double FunS   = IPH*ValS;
757      if(FunS > 320.) break;
758      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Thick[j];
759    } 
760    G4double InExp    = -hTotG2*Integ*stepB/dHS;
761    G4double expB     = std::exp(InExp);
762    G4double HRIE     = HadrReIm*InExp;
763    ReIntegrand[i]    = (1.-expB*std::cos(HRIE));           // Real part of the amplitude
764    ImIntegrand[i]    = expB*std::sin(HRIE);                // Imaginary part of the amplit
765  } 
766  InCoh     = 0.;                                           // incohirent (quasi-elastic)
767  ValB       = -stepB;
768  for(G4int k=0; k<NptB; k++)                               // Third integration (?)
769  {
770    ValB              += stepB;
771    InCoh             += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]);
772    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;      // Bessel0(Q*b)
773    ReSum             += J0qb*ReIntegrand[k];
774    ImSum             += J0qb*ImIntegrand[k];
775  }
776  //G4cout<<"GHAD:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
777  //      <<Q2<<",m="<<mb2G2<<G4endl;
778  InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2;
779  //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
780  //return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; // ds/do
781  return (ReSum*ReSum+ImSum*ImSum)*m2G10*stepB*stepB/12; // ds/dt
782} 
783
784G4double CHIPSDifElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
785{//      ============================================================================
786  static const G4int Npb      = 500;                        // A#of intergation points
787  //static const G4double mN = .9315;                         // Atomic Unit GeV
788  static const G4double mN = .938;                          // Atomic Unit GeV
789  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
790  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
791  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
792  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
793  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
794  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
795  G4double Lim        = 100*Re;                             // Integration accuracy limit
796  G4double Tb[Npb];                                         // Calculated T(b) array
797  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
798  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
799  G4double p          = hMom/1000.;                         // Momentum (GeV, inputParam)
800  G4double p2         = p*p;
801                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
802  G4double MassH2     = MassH*MassH;                        // Squared mass of the hadron
803  G4double hEnergy    = std::sqrt(p2+MassH2);               // Tot energy in GeV
804  G4double Q2         = aQ2/1000000.;                       // -t in GeV
805                G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
806  if(G4NucleiPropertiesTable::IsInTable(Z,A))
807                                MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.;    // Geant4 NuclearMass in GeV
808  G4double MassN2     = MassN*MassN;                        // Squared mass of the target
809  G4double S          = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s
810  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
811  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum (to norm CS)
812  //G4cout<<"2: P="<<CMMom<<",E="<<hEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl;
813  // The mean value of the total can be used
814  G4double p3         = p2*p;
815  G4double p4         = p2*p2;
816  G4double sp         = std::sqrt(p);
817  G4double p2s        = p2*sp;
818  G4double ap         = std::log(p);
819  G4double dl         = ap-3.;
820  G4double dl2        = dl*dl;
821  G4double shPPTot    = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
822  G4double shPNTot    = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
823                        +(38.+.3*dl2+8./p)/(1.+1.2/p3);     // SigPP in mb
824  G4double shNTot     = (Z*shPPTot+(A-Z)*shPNTot)/A;        // SighN in mb
825#ifdef debug
826  G4cout<<"CHIPS:SI,p="<<p<<",n="<<shNTot<<",P="<<shPPTot<<",N="<<shPNTot
827        <<",Z="<<Z<<",A="<<A<<G4endl;
828#endif
829  G4double shPPSl     = 8.*std::pow(p,.055)/(1.+3.64/p4);   // PPslope in GeV^-2
830  G4double shPNSl     = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
831  G4double shNSl      = (Z*shPPSl+(A-Z)*shPNSl)/A;          // B-slope in GeV^-2
832#ifdef debug
833  G4cout<<"CHIPS:SL,n="<<shNSl<<",P="<<shPPSl<<",N="<<shPNSl<<G4endl;
834#endif
835  G4double shNReIm    = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit
836  //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl;
837  // @@ End of temporary ^^^^^^^
838  G4double dHS        = shNSl+shNSl;                        // Working: doubled B-slope
839  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
840  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
841  G4double ReSum      = 0.;                                 // Integration of RePart of Amp
842  G4double ImSum      = 0.;                                 // Integration of ImPart of Amp
843  G4double ValB       = -stepB;
844  for(G4int i=0; i<Npb; i++)                                // First integration over b
845  {
846    ValB             += stepB;                              // An incident parameter
847    G4double ValB2    = ValB*ValB;                          // A working value b^2
848    G4double IPH      = ValB/shNSl;                         // A working value slope
849    G4double Integ    = 0.;                                 // Integral over ImpactParam.
850    G4double ValS     = 0.;                                 // Prototype of ImpactParameter
851    for(G4int j=1; j<Npb; j++)                              // Second integration over b
852    {
853      ValS           += stepB; //  back to fm               // Impact parameter GeV^-1
854      if(!i) Tb[j]    = CHIPS_Tb(A,ValS/f2Gm1)/f22G2;       // GeV^2, calculate only once
855      //if(!i) Tb[j]    = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once
856      G4double FunS   = IPH*ValS;                           // b1*b2/slope
857      if(FunS > Lim) break;                                 // To avoid NAN
858      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; // BessI0
859    } 
860    G4double InExp    = -hTotG2*Integ*stepB/dHS;            // Integrated absorption
861    G4double expB     = std::exp(InExp);                    // Exponential absorption
862    G4double HRIE     = shNReIm*InExp;                      // Phase shift
863    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;      // Bessel0(Q*b)
864    ReSum             += J0qb*(1.-expB*std::cos(HRIE));
865    ImSum             += J0qb*expB*std::sin(HRIE);
866  } 
867  //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
868  //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi;
869  //return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; // ds/do
870  return (ReSum*ReSum+ImSum*ImSum)*f22G2*stepB*stepB/12; // ds/dt
871}
872
873// Separate quasielastic calculation
874G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
875{//      =================================================================================
876  static const G4int Npb      = 500;                        // A#of intergation points
877  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
878  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
879  //static const G4double mb2G2 = 2.568;                     // Transform from mb to GeV^-2
880  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
881  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
882  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
883  G4double p          = hMom/1000.;                         // Momentum (GeV, inputParam)
884  G4double p2         = p*p;
885  G4double Q2         = aQ2/1000000.;                       // -t in GeV
886  // @@ Temporary only for nucleons
887  G4double p3         = p2*p;
888  G4double p4         = p2*p2;
889  G4double sp         = std::sqrt(p);
890  G4double p2s        = p2*sp;
891  G4double ap         = std::log(p);
892  G4double dl         = ap-3.;
893  G4double dl2        = dl*dl;
894  G4double shPPTot    = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb
895  G4double shPNTot    = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p
896                        +(38.+.3*dl2+8./p)/(1.+1.2/p3);     // SigPP in mb
897  G4double shNTot     = (Z*shPPTot+(A-Z)*shPNTot)/A;        // SighN in mb
898  G4double shPPSl     = 8.*std::pow(p,.055)/(1.+3.64/p4);   // PPslope in GeV^-2
899  G4double shPNSl     = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2
900  G4double shNSl      = (Z*shPPSl+(A-Z)*shPNSl)/A;          // B-slope in GeV^-2
901  G4double shNReIm    = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit
902  // @@ End of temporary ^^^^^^^
903  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
904  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
905  G4double InQE       = 0.;                                 // Quasielastic integral
906  G4double ValB       = -stepB;                             // Impact parameter
907  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
908  for(G4int k=0; k<Npb; k++)
909  {
910    ValB        += stepB;
911    //G4double Tb  = Thickness(A,ValB/f2Gm1,rAfm)/f22G2;      // Calculate T(b) only once
912    G4double Tb  = CHIPS_Tb(A,ValB/f2Gm1)/f22G2;            // GeV^2, calculate only once
913    InQE        += Tb*ValB*std::exp(-hTotG2*Tb);
914  }
915  //G4cout<<"CHPS:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
916  //      <<Q2<<",m="<<mb2G2<<G4endl;
917  return InQE*stepB*hTotG2*hTotG2*(1.+shNReIm*shNReIm)*std::exp(-shNSl*Q2)/8/mb2G2;
918}
919#endif
920
921int main()
922{
923#ifdef bestest
924  // *** Test of CHIPS implementation of Bessel functions (accuracy 1.E-8 and better) ***
925  const G4int nj=21;
926  const G4int ny=16;
927  const G4int ni=20;
928  const G4int nk=20;
929  const G4double jX[nj]=
930                                        {-5.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
931  const G4double yX[ny]={.1,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
932  const G4double iX[ni]=
933                    {0.,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.};
934  const G4double kX[nk]=
935                    {.1,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.};
936  const G4double vBJ0[nj]=
937                                              {-.1775968,-.3971498,-.2600520,0.2238908,0.7651976,1.0000000,0.7651977,
938                   0.2238908,-.2600520,-.3971498,-.1775968,0.1506453,0.3000793,0.1716508,
939                   -.0903336,-.2459358,-.1711903,0.0476893,0.2069261,0.1710735,-.0142245};
940  const G4double vBJ1[nj]=
941                                        {0.3275791,0.0660433,-.3390590,-.5767248,-.4400506,0.0000000,0.4400506,
942                   0.5767248,0.3390590,-.0660433,-.3275791,-.2766839,-.0046828,0.2346364,
943                   0.2453118,0.0434728,-.1767853,-.2234471,-.0703181,0.1333752,0.2051040};
944  const G4double vBY0[ny]=
945        {-1.534239,0.0882570,.51037567,.37685001,-.0169407,-.3085176,-.2881947,-.0259497,
946         0.2235215,0.2499367,0.0556712,-.1688473,-.2252373,-.0782079,0.1271926,0.2054743};
947  const G4double vBY1[ny]=
948                                {-6.458951,-.7812128,-.1070324,0.3246744,0.3979257,0.1478631,-.1750103,-.3026672,
949         -.1580605,0.1043146,0.2490154,0.1637055,-.0570992,-.2100814,-.1666448,0.0210736};
950  const G4double vBI0[ni]={1.0000000,1.0100250,1.0404018,1.0920453,1.1665149,
951                           1.2660658,1.3937256,1.5533951,1.7499807,1.9895593,
952                           2.2795852,3.2898391,4.8807925,7.3782035,11.301922,
953                           17.481172,27.239871,67.234406,427.56411,2815.7167};
954  const G4double vBI1[ni]={.00000000,.10050083,.20402675,.31370403,.43286480,
955                           .56515912,.71467794,.88609197,1.0848107,1.3171674,
956                           1.5906369,2.5167163,3.9533700,6.2058350,9.7594652,
957                           15.389221,24.335643,61.341937,399.87313,2670.9883};
958  const G4double vBK0[nk]={2.4270690,1.7527038,1.1145291,.77752208,.56534710,
959                           .42102445,.31850821,.24365506,.18795475,.14593140,
960                           .11389387,.062347553,.0347395,.019598897,.011159676,
961                           .0063998572,.0036910983,.0012439943,.00014647071,.000017780062};
962  const G4double vBK1[nk]={9.8538451,4.7759725,2.1843544,1.3028349,.86178163,
963                           .60190724,.43459241,.32083589,.24063392,.18262309,
964                                                                                                                                                                                                                        .13986588,.073890816,.040156431,.022239393,.012483499,
965                           .0070780949,.0040446134,.0013439197,.00015536921,.000018648773};
966  //Bessel J0('J',0); // CLHEP J/Y is possible insead of 'J'
967  //Bessel J1('J',1);
968  G4QBesIKJY QI0(BessI0); // CHIPS I/K/J/Y 0/1 Bessel functions
969  G4QBesIKJY QI1(BessI1);
970  G4QBesIKJY QJ0(BessJ0);
971  G4QBesIKJY QJ1(BessJ1);
972  G4QBesIKJY QK0(BessK0);
973  G4QBesIKJY QK1(BessK1);
974  G4QBesIKJY QY0(BessY0);
975  G4QBesIKJY QY1(BessY1);
976  G4cout<<"G4GCS: ***J0*** Test: ================================================"<<G4endl;
977  for(G4int j0=0; j0<nj; j0++)
978                {
979    G4double F=QJ0(jX[j0]);
980    G4double d=F-vBJ0[j0];
981    G4double r=std::abs(d/F);
982    G4cout<<"G4GCS: x="<<jX[j0]<<", J0="<<F<<" - "<<vBJ0[j0]<<" = "<<d<<", r="<<r<<G4endl;
983  }
984  G4cout<<"G4GCS: ***J1*** Test: ================================================"<<G4endl;
985  for(G4int j1=0; j1<nj; j1++)
986                {
987    G4double F=QJ1(jX[j1]);
988    G4double d=F-vBJ1[j1];
989    G4double r=std::abs(d/F);
990    G4cout<<"G4GCS: x="<<jX[j1]<<", J1="<<F<<" - "<<vBJ1[j1]<<" = "<<d<<", r="<<r<<G4endl;
991  }
992  G4cout<<"G4GCS: ***Y0*** Test: ================================================"<<G4endl;
993  for(G4int y0=0; y0<ny; y0++)
994                {
995    G4double F=QY0(yX[y0]);
996    G4double d=F-vBY0[y0];
997    G4double r=std::abs(d/F);
998    G4cout<<"G4GCS: x="<<yX[y0]<<", Y0="<<F<<" - "<<vBY0[y0]<<" = "<<d<<", r="<<r<<G4endl;
999  }
1000  G4cout<<"G4GCS: ***Y1*** Test: ================================================"<<G4endl;
1001  for(G4int y1=0; y1<ny; y1++)
1002                {
1003    G4double F=QY1(yX[y1]);
1004    G4double d=F-vBY1[y1];
1005    G4double r=std::abs(d/F);
1006    G4cout<<"G4GCS: x="<<yX[y1]<<", Y1="<<F<<" - "<<vBY1[y1]<<" = "<<d<<", r="<<r<<G4endl;
1007  }
1008  G4cout<<"G4GCS: ***I0*** Test: ================================================"<<G4endl;
1009  for(G4int i0=0; i0<ni; i0++)
1010                {
1011    G4double F=QI0(iX[i0]);
1012    G4double d=F-vBI0[i0];
1013    G4double r=std::abs(d/F);
1014    G4cout<<"G4GCS: x="<<iX[i0]<<", I0="<<F<<" - "<<vBI0[i0]<<" = "<<d<<", r="<<r<<G4endl;
1015  }
1016  G4cout<<"G4GCS: ***I1*** Test: ================================================"<<G4endl;
1017  for(G4int i1=0; i1<ni; i1++)
1018                {
1019    G4double F=QI1(iX[i1]);
1020    G4double d=F-vBI1[i1];
1021    G4double r=std::abs(d/F);
1022    G4cout<<"G4GCS: x="<<iX[i1]<<", I1="<<F<<" - "<<vBI1[i1]<<" = "<<d<<", r="<<r<<G4endl;
1023  }
1024  G4cout<<"G4GCS: ***K0*** Test: ================================================"<<G4endl;
1025  for(G4int k0=0; k0<nk; k0++)
1026                {
1027    G4double F=QK0(kX[k0]);
1028    G4double d=F-vBK0[k0];
1029    G4double r=std::abs(d/F);
1030    G4cout<<"G4GCS: x="<<kX[k0]<<", I1="<<F<<" - "<<vBK0[k0]<<" = "<<d<<", r="<<r<<G4endl;
1031  }
1032  G4cout<<"G4GCS: ***K1*** Test: ================================================"<<G4endl;
1033  for(G4int k1=0; k1<nk; k1++)
1034                {
1035    G4double F=QK1(kX[k1]);
1036    G4double d=F-vBK1[k1];
1037    G4double r=std::abs(d/F);
1038    G4cout<<"G4GCS: x="<<kX[k1]<<", I1="<<F<<" - "<<vBK1[k1]<<" = "<<d<<", r="<<r<<G4endl;
1039  }
1040  G4cout<<"End .................................................................."<<G4endl;
1041#endif
1042
1043  // Test of integrated cross sections
1044  //const G4int na=12;
1045  const G4int na=1;
1046  //
1047  const G4int np=7;
1048  const G4int nm=1;
1049  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
1050  //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei
1051  //                 He Al Pb
1052  const G4int A[na]={119}; // A's of target nuclei
1053  //                     p    n  pi+  pi-  K+  K-  antip
1054  const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles
1055  const G4double mom[nm]={120.}; // momentum in MeV/c
1056#ifdef integrc
1057  for(G4int ip=0; ip<np; ip++)
1058                {
1059    G4int PDG=pdg[ip];
1060    for(G4int im=0; im<nm; im++)
1061                  {
1062      CalculateParameters(PDG, mom[im]);
1063      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<":A,Tot,Inel,Prod,El,QEl"<<G4endl;
1064      for(G4int ia=0; ia<na; ia++)
1065                    {
1066        CalculateIntegralCrossSections(A[ia]);
1067        G4cout<<A[ia]<<" "<<TotalCrSec<<" "<<InelCrSec<<" "<<ProdCrSec
1068              <<" "<<ElasticCrSec<<" "<<QuasyElasticCrSec<<G4endl;
1069      }
1070    }
1071  }
1072#endif
1073
1074#ifdef eldiffcs
1075  const G4double ms=.000001;
1076  const G4int nt=200;
1077  G4double t[nt]; // Q2=-t in Mev^2
1078  const G4double tMin=200.;
1079  const G4double tMax=2000000.;
1080  const G4double ltMi=std::log(tMin);
1081  const G4double ltMa=std::log(tMax);
1082  const G4double dlt=(ltMa-ltMi)/(nt-1);
1083  G4double lt=ltMi-dlt;
1084  for(G4int it=0; it<nt; it++)
1085                {
1086    lt+=dlt;
1087    t[it]=std::exp(lt);
1088  }
1089  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
1090  //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei
1091  //                He Al Pb
1092  const G4int Z[na]={50}; // Z's of target nuclei
1093  // Test of differential ellastic cross sections
1094                Initialize();
1095  //for(G4int ip=0; ip<np; ip++)
1096  for(G4int ip=0; ip<1; ip++)
1097                {
1098    G4int PDG=pdg[ip];
1099    for(G4int im=0; im<nm; im++)
1100                  {
1101      CalculateParameters(PDG, mom[im]);
1102      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<G4endl;
1103      for(G4int ia=0; ia<na; ia++)
1104                    {
1105        G4cout<<"G4GCS:-------------------A="<<A[ia]<<G4endl;
1106        for(G4int it=0; it<nt; it++)
1107                      {
1108          G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]);
1109          //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); //Doesn't work
1110          G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]);
1111          G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
1112          G4double CQEl = CHIPSDifQuasiElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
1113                                                                  G4cout<<std::setw(12)<<ms*t[it]<<std::setw(12)<<Sig1<<std::setw(12)<<Sig2
1114                <<std::setw(12)<<InCoh<<std::setw(12)<<Sig3<<std::setw(12)<<CQEl<<G4endl;
1115        }
1116      }
1117    }
1118  }
1119#endif
1120  return 0;
1121}
Note: See TracBrowser for help on using the repository browser.