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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 52.4 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
319void Initialize()
320{
321  //static const G4double  sat = 3.03;      // Stirling saturation of Dubna
322  G4double Val = Factorials[0] = Mnoj[0] = 1.;
323  for (G4int i=1; i<nFact; i++)
324  {
325    G4double Si = 0.;
326             Val *= i;                    // Calculation of the factorial values
327    Factorials[i] = Val;
328    for (G4int l = 0; l<=i; l++)
329    {
330      G4double Sm = 1.;                   // sum of reversed binom coefficients (?)
331      for (G4int m=1; m<=i-l; m++) Sm += 1./binom(i-l,m);
332      //if(i>169)G4cout<<"G4GCS::Init:S="<<Sm<<", b("<<i<<","<<l<<")="<<binom(i,l)<<G4endl;
333      Si += Sm/binom(i,l);
334    } //  End of l LOOP
335    Mnoj[i] = Si;
336    //G4double d=Si-sat;
337    //G4double r=d/sat;
338    //if(i>100)G4cout<<"G4GCS::Init:i="<<i<<":"<<Si<<"-"<<sat<<" = "<<d<<", r="<<r<<G4endl;
339  } // End of i LOOP
340} // End of Initialization
341
342// Dubna paranmeters of nuclei
343void CalculateNuclearParameters(G4int Anuc)  // R1,R2,Pnuc,Aeff
344{//  ==============================================
345  G4double A=Anuc;
346  // Personal corrections fore some nuclei
347  if(Anuc == 4)                                      // Personal correction for He (!)
348  {
349    R1    = 5.5;
350    R2    = 3.7;
351    Pnucl = 0.4;   
352    Aeff  = 0.87;
353  }
354  else if(Anuc == 9)                                 // Personal correction for Be (!)
355  {
356    R1    = 9.0;
357    R2    = 7.0;
358    Pnucl = 0.190;
359    Aeff  = 0.9;
360  }
361  //else if(Anuc == 11)                              // Personal correction for B (!)
362  //{
363  //  R1    = 10.8;
364  //  R2    = 7.5;
365  //  Pnucl = 0.85;
366  //  Aeff  = 1.2;
367  //}
368                //else if(Anuc == 12)                              // Personal correction for C (!)
369  //{
370  //   R1    = 9.336;
371  //   R2    = 5.63;
372  //   Pnucl = 0.197;
373  //   Aeff  = 1.0;
374  //}
375  else if(Anuc == 16)                                // Personal correction for O (!)
376  {
377    R1    = 10.50;
378    R2    = 5.5;
379    Pnucl = 0.7;
380    Aeff  = 0.98;
381    //R1    = 11.3;
382    //R2    = 2.5;
383    //Pnucl = 0.75;
384    //Aeff  = 0.9;
385  }
386  else if(Anuc == 58)                                // Personal correction for Ni (!)
387  {
388    R1    = 15.0;
389    R2    = 9.9;
390    Pnucl = 0.45;
391    Aeff  = 0.85;
392  }
393  else if(Anuc == 90)                                // Personal correction for Zr (!)
394  {
395    //R1    = 16.5;
396    //R2    = 11.62;
397    //Pnucl = 0.4;
398    //Aeff  = 0.9;
399    R1    = 16.5;
400    R2    = 11.62;
401    Pnucl = 0.4;
402    Aeff  = 0.7;
403  }
404  else if(Anuc == 208)                               // Personal correction for Pb (!)
405  { 
406    //      R1 = 20.73; R2 = 15.74.
407    //R1    = 4.1408*std::pow(A,0.3018);
408    //R2    = 3.806*std::pow(A-10.068,0.2685);
409    //Pnucl = 0.9;
410    //Aeff  = 1.1;
411    R1    = 19.5;
412    R2    = 15.74;
413    Pnucl = 0.4;
414    Aeff  = 0.7;
415  }
416  else
417  {
418    R1 = 4.45*std::pow(A-1.,.309);                  // First diffractional radius
419    if(Anuc == 28) R1 *= .955;                       // Personal correction for Si (?!)
420    R2    = 2.3*std::pow(A,.36);                    // Second diffraction radius
421    Pnucl = 0.176+(.00167+.00000869*Anuc)*Anuc;      // Momentum of mean Fermi motion
422    Aeff  = 0.9;                                     // Effective screaning
423  }
424}
425
426// Independent transition method calculation. Has problems for A>95 (short cut).
427G4double DiffElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) // All in MeV
428{//      =================================================================
429  static const G4double eps = 0.000001;                     // Accuracy of calculations
430  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
431  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
432  static const G4double dR0   = 1.06+1.06;                  // Doubled nominal R0
433  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
434                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
435  G4double MassH2     = MassH*MassH;
436  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
437                //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
438  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
439  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
440  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
441  CalculateParameters(hPDG, HadrMom);                   
442  CalculateIntegralCS(A, HadrEnergy);
443  CalculateNuclearParameters(A);
444  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
445                //MassN           = A*0.938;                              // @@ This is bad!
446                G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
447  G4double MassN2   = MassN*MassN;
448  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
449  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
450  G4double CMMom = std::sqrt(EcmH*EcmH-MassH2);             // CM momentum
451
452  G4double Stot     = HadrTot*mb2G2;                        // in GeV^-2
453  G4double Bhad     = HadrSlope;                            // in GeV^-2
454  G4double Asq      = 1+HadrReIm*HadrReIm;                  // |M|^2/(ImM)^2
455  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
456  G4double R12      = R1*R1;
457  G4double R22      = R2*R2;
458  G4double R12B     = R12+Bhad+Bhad;
459  G4double R22B     = R22+Bhad+Bhad;
460  G4double R12Ap    = R12+20;
461  G4double R22Ap    = R22+20;
462  G4double R13Ap    = R12*R1/R12Ap;
463  G4double R23Ap    = R22*R2/R22Ap*Pnucl;
464  G4double R23dR13  = R23Ap/R13Ap;
465  G4double R12Apd   = 2./R12Ap;
466  G4double R22Apd   = 2./R22Ap;
467  G4double R122Apd  = (R12Apd+R22Apd)/2;
468  G4double dR0H     = dR0*HadrEnergy/4;
469  G4double DDSec1p  = DDSect2+DDSect3*std::log(dR0H/R1);
470  G4double DDSec2p  = DDSect2+DDSect3*std::log(dR0H/std::sqrt((R12+R22)/2));
471  G4double DDSec3p  = DDSect2+DDSect3*std::log(dR0H/R2);
472  G4double Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;    // Some questionable norming
473  G4double R13      = R12*R1/R12B;
474  G4double R23      = Pnucl*R22*R2/R22B;
475  G4double norFac   = Stot/twopi/Norm;               // totCS (in GeV^-2) is used
476  G4double Unucl    = norFac*R13;
477  G4double UnuclScr = norFac*R13Ap;
478  G4double SinFi    = HadrReIm/Rho2;
479  G4double FiH      = std::asin(SinFi);
480  G4double N        = -1.;
481  G4double N2       = R23/R13;
482  G4double ImElA0   = 0.;
483  G4double ReElA0   = 0.;
484  G4double Tot1     = 0.;
485  for(G4int i=1; i<=A; i++)                          // @@ Make separately for n and p
486  {
487    N              *= -Unucl*(A-i+1)*Rho2/i;
488    G4double N4     = 1.;
489    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i;
490    G4double medTot = R12B/i;
491    for(G4int l=1; l<=i; l++)
492    {
493      G4double exp1 = l/R22B+(i-l)/R12B;
494      N4           *= -(i-l+1)*N2/l;
495      G4double Nexp = N4/exp1;
496      Prod1        += Nexp*std::exp(-Q2/exp1/4);
497      medTot       += Nexp;
498    }
499    G4double FiHi   = FiH*i;
500    G4double nCos   = N*std::cos(FiHi);
501    ReElA0 += Prod1*N*std::sin(FiHi);
502    ImElA0 += Prod1*nCos;
503    Tot1           += medTot*nCos;
504    if(std::abs(Prod1*N/ImElA0) < eps) break;
505  }
506  ImElA0   *= piMG;                                  // The amplitude in mb
507  ReElA0   *= piMG;                                  // The amplitude in mb
508  Tot1             *= piMG+piMG;
509  G4double N1p      = 1.;
510  G4double R13Ap1   = DDSec1p*R13Ap*R13Ap/2;
511  G4double R22Ap1   = DDSec2p*R13Ap*R23Ap;
512  G4double R23Ap1   = DDSec3p*R23Ap*R23Ap*R22Ap/2;
513  G4double R13Ap2   = R13Ap1*R12Ap/4;
514  G4double R22Ap2   = R22Ap1/R122Apd/2;
515  G4double R23Ap2   = R23Ap1*R22Ap/4;
516  G4double Q28      = Q2/8;
517  G4double Q24      = Q28+Q28;
518  G4double Din1     = R13Ap2*std::exp(-Q28*R12Ap)-R22Ap2*std::exp(-(Q28+Q28)/R122Apd)+
519                                                  R23Ap2*std::exp(-Q28*R22Ap);   // i=0 start value
520  G4double DTot1    = R13Ap2-R22Ap2+R23Ap2;          // i=0 start value
521  if (A<96)
522                {
523    for(G4int i=1; i<A-1; i++)
524    {
525      N1p              *= -UnuclScr*(A-i-1)/i*Rho2;
526      G4double N2p      = 1.;
527      G4double Din2     = 0.;
528      G4double DmedTot  = 0.;
529      G4double BinCoeff = 1.;                        // Start for Binomial Coefficient
530      for(G4int l = 0; l<=i; l++) 
531      {
532        if(l) BinCoeff *= (i-l+1)/l;
533        G4double exp1   = l/R22B+(i-l)/R12B;
534        G4double exp1p  = exp1+R12Apd;
535        G4double exp2p  = exp1+R122Apd;
536        G4double exp3p  = exp1+R22Apd;
537        G4double NBinC  = N2p*BinCoeff;
538        Din2           += NBinC*(R13Ap1/exp1p*std::exp(-Q24/exp1p)-
539                                 R22Ap1/exp2p*std::exp(-Q24/exp2p)+
540                                 R23Ap1/exp3p*std::exp(-Q24/exp3p));
541        DmedTot        += NBinC * (R13Ap1/exp1p - R22Ap1/exp2p + R23Ap1/exp3p);
542        N2p            *= -R23dR13;
543      }
544      G4double comFac   = N1p*Mnoj[i]/(i+2)/(i+1)*std::cos(FiH*i);
545                                  Din1             += Din2*comFac;
546      DTot1            += DmedTot*comFac;
547      if(std::abs(Din2*N1p/Din1) < eps) break;
548    }
549    G4double cFac       = A*(A-1)*4/Norm/Norm;
550    Din1               *= -cFac;
551    DTot1              *= cFac;
552  }
553  //else Din1=0.;
554
555  G4double Corr0      = Tot00/Tot1;
556  ImElA0             *= Corr0;
557                       
558  return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2;
559} // End of DiffElasticCS
560
561G4double CHIPSDiffElCS(G4int hPDG, G4double HadrMom, G4int Z, G4int A, G4double aQ2) // MeV
562{//      =================================================================
563  static const G4double eps = 0.000001;                     // Accuracy of calculations
564  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
565  static const G4double piMG  = pi/mb2G2;                   // PiTrans from GeV^-2 to mb
566  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
567                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
568  G4double MassH2     = MassH*MassH;
569  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
570                //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl;
571  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
572  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
573  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
574  CalculateParameters(hPDG, HadrMom);                   
575  CalculateIntegralCS(A, HadrEnergy);
576  CalculateNuclearParameters(A);
577  G4double Q2 = aQ2/1000/1000;                              // in GeV^2
578                //G4doubleMassN     = A*0.938;                              // @@ This is bad!
579                G4double MassN    = A*0.9315;                             // @@ Atomic Unit is bad too
580  if(G4NucleiPropertiesTable::IsInTable(Z,A))
581                                MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.;    // Geant4 NuclearMass
582  G4double MassN2   = MassN*MassN;
583  G4double S        = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s
584  G4double EcmH     = (S-MassN2+MassH2)/2/std::sqrt(S);     // CM energy of a hadron
585  G4double CMMom    = std::sqrt(EcmH*EcmH-MassH2);          // CM momentum
586  G4double Stot     = HadrTot*mb2G2;                        // in GeV^-2
587  G4double Bhad     = HadrSlope;                            // in GeV^-2
588  G4double Asq      = 1+HadrReIm*HadrReIm;                  // |M|^2/(ImM)^2
589  G4double Rho2     = std::sqrt(Asq);                       // M/ImM
590  G4double R12      = R1*R1;
591  G4double R22      = R2*R2;
592  G4double R12B     = R12+Bhad+Bhad;                        // Slope is used
593  G4double R22B     = R22+Bhad+Bhad;                        // Slope is used
594  G4double Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;           // Some questionable norming
595  G4double R13      = R12*R1/R12B;                          // Slope is used
596  G4double R23      = Pnucl*R22*R2/R22B;                    // Slope is used
597  G4double norFac   = Stot/twopi/Norm;                      // totCS (in GeV^-2) is used
598  G4double Unucl    = norFac*R13;
599  G4double SinFi    = HadrReIm/Rho2;
600  G4double FiH      = std::asin(SinFi);
601  G4double N        = -1.;
602  G4double N2       = R23/R13;                              // Slope is used
603  G4double ImElA0   = 0.;
604  G4double ReElA0   = 0.;
605  G4double Tot1     = 0.;
606  for(G4int i=1; i<=A; i++)
607  {
608    N              *= -Unucl*(A-i+1)*Rho2/i;
609    G4double N4     = 1.;
610    G4double Prod1  = std::exp(-Q2*R12B/i/4)*R12B/i;        // Slope is used
611    G4double medTot = R12B/i;                               // Slope is used
612    for(G4int l=1; l<=i; l++)
613    {
614      G4double exp1 = l/R22B+(i-l)/R12B;                    // Slope is used
615      N4           *= -(i-l+1)*N2/l;                        // Slope is used
616      G4double Nexp = N4/exp1;
617      Prod1        += Nexp*std::exp(-Q2/exp1/4);
618      medTot       += Nexp;
619    }
620    G4double FiHi   = FiH*i;
621    G4double nCos   = N*std::cos(FiHi);
622    ReElA0 += Prod1*N*std::sin(FiHi);
623    ImElA0 += Prod1*nCos;
624    Tot1           += medTot*nCos;
625    if(std::abs(Prod1*N/ImElA0) < eps) break;
626  }
627  ImElA0         *= piMG;                          // The amplitude in mb
628  ReElA0         *= piMG;                          // The amplitude in mb
629  Tot1           *= piMG+piMG;
630  G4double Corr0  = Tot00/Tot1;
631  ImElA0         *= Corr0;
632                       
633  return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2;
634} // End of DiffElasticCS
635
636G4double CHIPS_Tb(G4int  A, G4double b)              // T(b) in fm-2
637{
638  static G4int Am=0;                                 // Associative memory value A
639  static G4double B=0.;                              // Effective edge
640  static G4double C=0.;                              // Normalization factor
641  static G4double D=0.;                              // Slope parameter
642  if(A!=Am)
643                {
644    B=.0008*A*A;                                     // no units
645    D=.42*std::pow(A,-.26);                               // fm^-2
646    C=A*D/pi/std::log(1.+B);                         // fm^-2
647  }
648  G4double E=B*std::exp(-D*b*b);
649  return C*E/(1.+E);                                 // T(b) in fm^-2
650}
651
652G4double Thickness(G4int A, G4double b, G4double R)  // T(b) in fm^-2
653{//      ==========================================  // rAmax=2.5 is a GlobStatConstPar
654  static const G4double bTh  = .545;                 // Surface thicknes
655  static const G4double bTh2 = bTh*bTh;              // SquaredSurface thicknes
656  G4double b2    = b*b;                              // Squared impact parameter
657  G4double dr    = rAmax*R/(NptB-1);                 // Step of integration in z
658  G4double R2    = R*R;                              // Working parameter
659  G4double Norm  = .75/pi/R2/R/(1.+bTh2*pi2/R2);     // Corrected Volume (?)
660  G4double r     = -dr;                              // Pre value for the radius
661  G4double SumZ  = 0.;                               // Prototype of the result
662  //G4double SumN  = 0.;                             // is not used
663  for(G4int k=0; k<NptB; k++)
664  {
665    r    += dr;
666    SumZ += 1./(1.+std::exp((std::sqrt(b2+r*r)-R)/bTh)); // Woods-Saxon density
667    //SumN += r*r/(1.+std::exp((r-rAfm)/bTh));       // is not used
668  }
669  //SumN *= (twopi+twopi)*dr*Norm;                   // is not used
670  return SumZ*(A+A)*dr*Norm;                         // T(b) in fm^-2
671} // End of Thickness
672
673G4double CoherentDifElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) //in MeV
674{//      =========================================================================
675  //static const G4double eps = 0.000001;                   // CalculationAccuracy@@NotUsed
676  //static const G4double mN = .9315;                      // Atomic Unit GeV
677  static const G4double mN = .938;                      // Atomic Unit GeV
678  static const G4double mb2G2 = 2.568;                      // Transform from mb to GeV^-2
679  static const G4double incrf = 10;            // Increase factor of radius for intergation
680  static const G4double m2G10 = incrf*mb2G2;                // Big conversion factor
681  static const G4double s2G10 = std::sqrt(m2G10);           // sqrt Big conversion factor
682  G4double ReIntegrand[NptB],ImIntegrand[NptB],Thick[NptB]; // Calculated arrays
683  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
684  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
685  G4double hMom       = HadrMom/1000.;                      // Momentum (GeV, inputParam)
686                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
687  G4double MassH2     = MassH*MassH;
688  G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
689  if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");
690  if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");
691  //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)");
692  G4double Re         = std::pow(A, 0.33333333);
693  G4double Re2        = Re*Re;
694  CalculateParameters(hPDG, HadrMom);                   
695  G4double Q2 = aQ2/1000/1000;                              //  GeV
696  // Individual corrections for nuclei which had data
697  if    (A==208) r0   = 1.125;
698  else if(A==90) r0   = 1.12;
699  else if(A==64) r0   = 1.1;
700  else if(A==58) r0   = 1.09;
701  else if(A==48) r0   = 1.07;
702  else if(A==40) r0   = 1.15;
703  else if(A==28) r0   = 0.93;
704  else if(A==16) r0   = 0.92;
705  else if(A==12) r0   = 0.80;
706                else           r0   = 1.16*(1.-1.16/Re2);  // For other nuclei which have not been tested
707  G4double MassN      = mN;                              // A * AtomicUnit(GeV)
708  G4double MassN2     = MassN*MassN;
709  G4double S          = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S
710  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
711  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum
712  G4double rAfm       = r0*Re;                              // Nuclear radius
713  G4double rAGeV      = rAfm*s2G10;                         // Big integration radius
714  G4double stepB      = rAmax*rAGeV/(NptB-1);               // dr step of integration
715  G4double hTotG2     = HadrTot*mb2G2;
716  G4double ReSum      = 0.;
717  G4double ImSum      = 0.;
718  G4double ValB       = -stepB;
719  for(G4int i=0; i<NptB; i++)
720  {
721    ValB             += stepB;                              // An incident parameter
722    G4double ValB2    = ValB*ValB;                          // A working value
723    G4double IPH      = ValB/HadrSlope;                     // A working value
724    G4double dHS      = HadrSlope+HadrSlope;                // A working value
725    G4double Integ    = 0.;
726    G4double ValS     = 0.;
727    for(G4int j=1; j<NptB; j++)
728    {
729      ValS           += stepB;                              // Impact parameter GeV^-1
730      if(!i) Thick[j] = Thickness(A,ValS/s2G10,rAfm)/m2G10; // Calculate only once
731      G4double FunS   = IPH*ValS;
732      if(FunS > 320.) break;
733      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Thick[j];
734    } 
735    G4double InExp    = -hTotG2*Integ*stepB/dHS;
736    G4double expB     = std::exp(InExp);
737    G4double HRIE     = HadrReIm*InExp;
738    ReIntegrand[i]    = (1.-expB*std::cos(HRIE));
739    ImIntegrand[i]    = expB*std::sin(HRIE);
740  } 
741  InCoh     = 0.;                                           // incohirent (quasi-elastic)
742  ValB       = -stepB;
743  for(G4int k=0; k<NptB; k++)
744  {
745    ValB              += stepB;
746    InCoh             += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]);
747    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;
748    ReSum             += J0qb*ReIntegrand[k];
749    ImSum             += J0qb*ImIntegrand[k];
750  }
751  //G4cout<<"GHAD:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
752  //      <<Q2<<",m="<<mb2G2<<G4endl;
753  InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2;
754  //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
755  return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi;
756} 
757
758G4double CHIPSDifElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
759{//      ============================================================================
760  static const G4int Npb      = 500;                        // A#of intergation points
761  //static const G4double mN = .9315;                         // Atomic Unit GeV
762  static const G4double mN = .938;                          // Atomic Unit GeV
763  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
764  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
765  //static const G4double mb2G2 = 2.568;                     // Transform from mb to GeV^-2
766  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
767  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
768  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
769  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
770  G4double Lim        = 50.*Re;                             // Integration accuracy limit
771  G4double Tb[Npb];                                         // Calculated T(b) array
772  G4QBesIKJY QI0(BessI0);                                   // I0 Bessel function
773  G4QBesIKJY QJ0(BessJ0);                                   // J0 Bessel function
774           hMom       = hMom/1000.;                         // Momentum (GeV, inputParam)
775                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
776  G4double MassH2     = MassH*MassH;                        // Squared mass of the hadron
777  G4double hEnergy    = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
778  G4double Q2         = aQ2/1000000.;                       // -t in GeV
779  G4double MassN      = mN;                                 // AtomicUnit(GeV)[prototype]
780  G4double MassN2     = MassN*MassN;                        // Squared mass of the target
781  G4double S          = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s
782  G4double EcmH       = (S-MassN2+MassH2)/2/std::sqrt(S);   // Hadron CM Energy
783  G4double CMMom      = std::sqrt(EcmH*EcmH-MassH2);        // CM momentum (to norm CS)
784  // @@ Temporary only for nucleons
785  //G4cout<<"CHPS:E="<<hEnergy<<",dM="<<2*MassN<<",sN="<<MassN2<<",sH="<<MassH2<<G4endl;
786  G4double shNTot     = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb
787  G4double shNSl      = 5.44+.88*std::log(S);               // B-slope in GeV^-2
788  G4double shNReIm    = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit
789  //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl;
790  // @@ End of temporary ^^^^^^^
791  G4double dHS        = shNSl+shNSl;                        // Working: doubled B-slope
792  //G4double rAfm       = 0.;
793  //if    (A==208) rAfm = 1.125*Re;
794  //else if(A==90) rAfm = 1.12*Re;
795  //else if(A==64) rAfm = 1.1*Re;
796  //else if(A==58) rAfm = 1.09*Re;
797  //else if(A==48) rAfm = 1.07*Re;
798  //else if(A==40) rAfm = 1.15*Re;
799  //else if(A==28) rAfm = 0.93*Re;
800  //else if(A==16) rAfm = 0.92*Re;
801  //else if(A==12) rAfm = 0.80*Re;
802                //else           rAfm = 1.16*(Re-1.16/Re);                  // For other nuclei
803  //G4double stepB      = 2.5*rAfm*f2Gm1/(Npb-1);           // in GeV^-1, step of integral
804  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
805  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
806  G4double ReSum      = 0.;                                 // Integration of RePart of Amp
807  G4double ImSum      = 0.;                                 // Integration of ImPart of Amp
808  G4double ValB       = -stepB;
809  for(G4int i=0; i<Npb; i++)
810  {
811    ValB             += stepB;                              // An incident parameter
812    G4double ValB2    = ValB*ValB;                          // A working value
813    G4double IPH      = ValB/shNSl;                         // A working value
814    G4double Integ    = 0.;                                 // Integral over ImpactParam.
815    G4double ValS     = 0.;                                 // Prototype of ImpactParameter
816    for(G4int j=1; j<Npb; j++)
817    {
818      ValS           += stepB; //  back to fm               // Impact parameter GeV^-1
819      if(!i) Tb[j]    = CHIPS_Tb(A,ValS/f2Gm1)/f22G2;       // GeV^2, calculate only once
820      //if(!i) Tb[j]    = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once
821      G4double FunS   = IPH*ValS;                           // Working product
822      if(FunS > Lim) break;                                 // (?)
823      Integ          += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j];
824    } 
825    G4double InExp    = -hTotG2*Integ*stepB/dHS;            // Working product
826    G4double expB     = std::exp(InExp);                    // Workung sqrt
827    G4double HRIE     = shNReIm*InExp;                      // Phase shift
828    G4double J0qb      = QJ0(std::sqrt(Q2)*ValB)*ValB;
829    ReSum             += J0qb*(1.-expB*std::cos(HRIE));
830    ImSum             += J0qb*expB*std::sin(HRIE);
831  } 
832  //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl;
833  //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi;
834  return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi;
835} 
836G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2)
837{//      =================================================================================
838  static const G4int Npb      = 500;                        // A#of intergation points
839  //static const G4double mN = .9315;                         // Atomic Unit GeV
840  static const G4double mN = .938;                          // Mass of proton GeV
841  static const G4double hc2   = .3893793;                   // Transform from GeV^-2 to mb
842  static const G4double mb2G2 = 1./hc2;                     // Transform from mb to GeV^-2
843  //static const G4double mb2G2 = 2.568;                     // Transform from mb to GeV^-2
844  static const G4double f22mb = 10;                         // Transform from fermi^2 to mb
845  static const G4double f22G2 = f22mb*mb2G2;                // Transform from fm2 to GeV^-2
846  static const G4double f2Gm1 = std::sqrt(f22G2);           // Transform from fm to GeV^-1
847           hMom       = hMom/1000.;                         // Momentum (GeV, inputParam)
848                G4double MassH      = G4QPDGCode(hPDG).GetMass()/1000.;   // Hadron Mass in GeV
849  G4double MassH2     = MassH*MassH;                        // Squared mass of the hadron
850  G4double hEnergy    = std::sqrt(hMom*hMom+MassH2);        // Tot energy in GeV
851  G4double Q2         = aQ2/1000000.;                       // -t in GeV
852  G4double MassN      = mN;                                 // A*AtomicUnit(GeV)[prototype]
853  //if(G4NucleiPropertiesTable::IsInTable(Z,A))
854                //              MassN=G4NucleiProperties::GetNuclearMass(A,Z)/A/1000.;  // Geant4 NuclearMass/A
855  G4double MassN2     = MassN*MassN;                        // Squared mass of the target
856  G4double S          = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s
857  // @@ Temporary only for nucleons
858  G4double shNTot     = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb
859  G4double shNSl      = 6.44+.88*std::log(S);               // B-slope in GeV^-2
860  G4double shNReIm    = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit
861  // @@ End of temporary ^^^^^^^
862  G4double Re         = std::pow(A,.33333333);              // A-dep coefficient
863  G4double stepB      = (Re+Re+2.7)*f2Gm1/(Npb-1);          // in GeV^-1, step of integral
864  G4double InQE       = 0.;                                 // Quasielastic integral
865  G4double ValB       = -stepB;                             // Impact parameter
866  G4double hTotG2     = shNTot*mb2G2;                       // sigma_hN in GeV^-2
867  for(G4int k=0; k<Npb; k++)
868  {
869    ValB        += stepB;
870    //G4double Tb  = Thickness(A,ValB/f2Gm1,rAfm)/f22G2;      // Calculate T(b) only once
871    G4double Tb  = CHIPS_Tb(A,ValB/f2Gm1)/f22G2;            // GeV^2, calculate only once
872    InQE        += Tb*ValB*std::exp(-hTotG2*Tb);
873  }
874  //G4cout<<"CHPS:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2="
875  //      <<Q2<<",m="<<mb2G2<<G4endl;
876  return InQE*stepB*hTotG2*hTotG2*(1.+shNReIm*shNReIm)*std::exp(-shNSl*Q2)/8/mb2G2;
877}
878#endif
879
880int main()
881{
882#ifdef bestest
883  // *** Test of CHIPS implementation of Bessel functions (accuracy 1.E-8 and better) ***
884  const G4int nj=21;
885  const G4int ny=16;
886  const G4int ni=20;
887  const G4int nk=20;
888  const G4double jX[nj]=
889                                        {-5.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
890  const G4double yX[ny]={.1,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.};
891  const G4double iX[ni]=
892                    {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.};
893  const G4double kX[nk]=
894                    {.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.};
895  const G4double vBJ0[nj]=
896                                              {-.1775968,-.3971498,-.2600520,0.2238908,0.7651976,1.0000000,0.7651977,
897                   0.2238908,-.2600520,-.3971498,-.1775968,0.1506453,0.3000793,0.1716508,
898                   -.0903336,-.2459358,-.1711903,0.0476893,0.2069261,0.1710735,-.0142245};
899  const G4double vBJ1[nj]=
900                                        {0.3275791,0.0660433,-.3390590,-.5767248,-.4400506,0.0000000,0.4400506,
901                   0.5767248,0.3390590,-.0660433,-.3275791,-.2766839,-.0046828,0.2346364,
902                   0.2453118,0.0434728,-.1767853,-.2234471,-.0703181,0.1333752,0.2051040};
903  const G4double vBY0[ny]=
904        {-1.534239,0.0882570,.51037567,.37685001,-.0169407,-.3085176,-.2881947,-.0259497,
905         0.2235215,0.2499367,0.0556712,-.1688473,-.2252373,-.0782079,0.1271926,0.2054743};
906  const G4double vBY1[ny]=
907                                {-6.458951,-.7812128,-.1070324,0.3246744,0.3979257,0.1478631,-.1750103,-.3026672,
908         -.1580605,0.1043146,0.2490154,0.1637055,-.0570992,-.2100814,-.1666448,0.0210736};
909  const G4double vBI0[ni]={1.0000000,1.0100250,1.0404018,1.0920453,1.1665149,
910                           1.2660658,1.3937256,1.5533951,1.7499807,1.9895593,
911                           2.2795852,3.2898391,4.8807925,7.3782035,11.301922,
912                           17.481172,27.239871,67.234406,427.56411,2815.7167};
913  const G4double vBI1[ni]={.00000000,.10050083,.20402675,.31370403,.43286480,
914                           .56515912,.71467794,.88609197,1.0848107,1.3171674,
915                           1.5906369,2.5167163,3.9533700,6.2058350,9.7594652,
916                           15.389221,24.335643,61.341937,399.87313,2670.9883};
917  const G4double vBK0[nk]={2.4270690,1.7527038,1.1145291,.77752208,.56534710,
918                           .42102445,.31850821,.24365506,.18795475,.14593140,
919                           .11389387,.062347553,.0347395,.019598897,.011159676,
920                           .0063998572,.0036910983,.0012439943,.00014647071,.000017780062};
921  const G4double vBK1[nk]={9.8538451,4.7759725,2.1843544,1.3028349,.86178163,
922                           .60190724,.43459241,.32083589,.24063392,.18262309,
923                                                                                                                                                                                                                        .13986588,.073890816,.040156431,.022239393,.012483499,
924                           .0070780949,.0040446134,.0013439197,.00015536921,.000018648773};
925  //Bessel J0('J',0); // CLHEP J/Y is possible insead of 'J'
926  //Bessel J1('J',1);
927  G4QBesIKJY QI0(BessI0); // CHIPS I/K/J/Y 0/1 Bessel functions
928  G4QBesIKJY QI1(BessI1);
929  G4QBesIKJY QJ0(BessJ0);
930  G4QBesIKJY QJ1(BessJ1);
931  G4QBesIKJY QK0(BessK0);
932  G4QBesIKJY QK1(BessK1);
933  G4QBesIKJY QY0(BessY0);
934  G4QBesIKJY QY1(BessY1);
935  G4cout<<"G4GCS: ***J0*** Test: ================================================"<<G4endl;
936  for(G4int j0=0; j0<nj; j0++)
937                {
938    G4double F=QJ0(jX[j0]);
939    G4double d=F-vBJ0[j0];
940    G4double r=std::abs(d/F);
941    G4cout<<"G4GCS: x="<<jX[j0]<<", J0="<<F<<" - "<<vBJ0[j0]<<" = "<<d<<", r="<<r<<G4endl;
942  }
943  G4cout<<"G4GCS: ***J1*** Test: ================================================"<<G4endl;
944  for(G4int j1=0; j1<nj; j1++)
945                {
946    G4double F=QJ1(jX[j1]);
947    G4double d=F-vBJ1[j1];
948    G4double r=std::abs(d/F);
949    G4cout<<"G4GCS: x="<<jX[j1]<<", J1="<<F<<" - "<<vBJ1[j1]<<" = "<<d<<", r="<<r<<G4endl;
950  }
951  G4cout<<"G4GCS: ***Y0*** Test: ================================================"<<G4endl;
952  for(G4int y0=0; y0<ny; y0++)
953                {
954    G4double F=QY0(yX[y0]);
955    G4double d=F-vBY0[y0];
956    G4double r=std::abs(d/F);
957    G4cout<<"G4GCS: x="<<yX[y0]<<", Y0="<<F<<" - "<<vBY0[y0]<<" = "<<d<<", r="<<r<<G4endl;
958  }
959  G4cout<<"G4GCS: ***Y1*** Test: ================================================"<<G4endl;
960  for(G4int y1=0; y1<ny; y1++)
961                {
962    G4double F=QY1(yX[y1]);
963    G4double d=F-vBY1[y1];
964    G4double r=std::abs(d/F);
965    G4cout<<"G4GCS: x="<<yX[y1]<<", Y1="<<F<<" - "<<vBY1[y1]<<" = "<<d<<", r="<<r<<G4endl;
966  }
967  G4cout<<"G4GCS: ***I0*** Test: ================================================"<<G4endl;
968  for(G4int i0=0; i0<ni; i0++)
969                {
970    G4double F=QI0(iX[i0]);
971    G4double d=F-vBI0[i0];
972    G4double r=std::abs(d/F);
973    G4cout<<"G4GCS: x="<<iX[i0]<<", I0="<<F<<" - "<<vBI0[i0]<<" = "<<d<<", r="<<r<<G4endl;
974  }
975  G4cout<<"G4GCS: ***I1*** Test: ================================================"<<G4endl;
976  for(G4int i1=0; i1<ni; i1++)
977                {
978    G4double F=QI1(iX[i1]);
979    G4double d=F-vBI1[i1];
980    G4double r=std::abs(d/F);
981    G4cout<<"G4GCS: x="<<iX[i1]<<", I1="<<F<<" - "<<vBI1[i1]<<" = "<<d<<", r="<<r<<G4endl;
982  }
983  G4cout<<"G4GCS: ***K0*** Test: ================================================"<<G4endl;
984  for(G4int k0=0; k0<nk; k0++)
985                {
986    G4double F=QK0(kX[k0]);
987    G4double d=F-vBK0[k0];
988    G4double r=std::abs(d/F);
989    G4cout<<"G4GCS: x="<<kX[k0]<<", I1="<<F<<" - "<<vBK0[k0]<<" = "<<d<<", r="<<r<<G4endl;
990  }
991  G4cout<<"G4GCS: ***K1*** Test: ================================================"<<G4endl;
992  for(G4int k1=0; k1<nk; k1++)
993                {
994    G4double F=QK1(kX[k1]);
995    G4double d=F-vBK1[k1];
996    G4double r=std::abs(d/F);
997    G4cout<<"G4GCS: x="<<kX[k1]<<", I1="<<F<<" - "<<vBK1[k1]<<" = "<<d<<", r="<<r<<G4endl;
998  }
999  G4cout<<"End .................................................................."<<G4endl;
1000#endif
1001
1002  // Test of integrated cross sections
1003  //const G4int na=12;
1004  const G4int na=1;
1005  //
1006  const G4int np=7;
1007  const G4int nm=1;
1008  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
1009  //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei
1010  //                 He Al Pb
1011  const G4int A[na]={208}; // A's of target nuclei
1012  //                     p    n  pi+  pi-  K+  K-  antip
1013  const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles
1014  const G4double mom[nm]={1090.}; // momentum in MeV/c
1015#ifdef integrc
1016  for(G4int ip=0; ip<np; ip++)
1017                {
1018    G4int PDG=pdg[ip];
1019    for(G4int im=0; im<nm; im++)
1020                  {
1021      CalculateParameters(PDG, mom[im]);
1022      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<":A,Tot,Inel,Prod,El,QEl"<<G4endl;
1023      for(G4int ia=0; ia<na; ia++)
1024                    {
1025        CalculateIntegralCrossSections(A[ia]);
1026        G4cout<<A[ia]<<" "<<TotalCrSec<<" "<<InelCrSec<<" "<<ProdCrSec
1027              <<" "<<ElasticCrSec<<" "<<QuasyElasticCrSec<<G4endl;
1028      }
1029    }
1030  }
1031#endif
1032
1033#ifdef eldiffcs
1034  const G4double ms=.000001;
1035  const G4int nt=200;
1036  G4double t[nt]; // Q2=-t in Mev^2
1037  const G4double tMin=200.;
1038  const G4double tMax=2000000.;
1039  const G4double ltMi=std::log(tMin);
1040  const G4double ltMa=std::log(tMax);
1041  const G4double dlt=(ltMa-ltMi)/(nt-1);
1042  G4double lt=ltMi-dlt;
1043  for(G4int it=0; it<nt; it++)
1044                {
1045    lt+=dlt;
1046    t[it]=std::exp(lt);
1047  }
1048  ////                He Be C O  Al Ti Ni Cu  Sn Ta  Pb   U
1049  //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei
1050  //                He Al Pb
1051  const G4int Z[na]={82}; // Z's of target nuclei
1052  // Test of differential ellastic cross sections
1053                Initialize();
1054  //for(G4int ip=0; ip<np; ip++)
1055  for(G4int ip=0; ip<1; ip++)
1056                {
1057    G4int PDG=pdg[ip];
1058    for(G4int im=0; im<nm; im++)
1059                  {
1060      CalculateParameters(PDG, mom[im]);
1061      G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<G4endl;
1062      for(G4int ia=0; ia<na; ia++)
1063                    {
1064        G4cout<<"G4GCS:-------------------A="<<A[ia]<<G4endl;
1065        for(G4int it=0; it<nt; it++)
1066                      {
1067          G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]);
1068          //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]);
1069          G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]);
1070          G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
1071          G4double CQEl = CHIPSDifQuasiElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
1072                                                                  G4cout<<std::setw(12)<<ms*t[it]<<std::setw(12)<<Sig1<<std::setw(12)<<Sig2
1073                <<std::setw(12)<<InCoh<<std::setw(12)<<Sig3<<std::setw(12)<<CQEl<<G4endl;
1074        }
1075      }
1076    }
1077  }
1078#endif
1079  return 0;
1080}
Note: See TracBrowser for help on using the repository browser.