ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChipsProtonInelasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ChipsProtonInelasticXS.cc
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 // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
28 //
29 //
30 // G4 Physics class: G4ChipsProtonInelasticXS for gamma+A cross sections
31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03
32 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04
33 //
34 //
35 // ****************************************************************************************
36 // Short description: Cross-sections extracted (by W.Pokorski) from the CHIPS package for
37 // proton-nuclear interactions. Original author: M. Kossov
38 // -------------------------------------------------------------------------------------
39 //
40 
41 
42 #include "G4ChipsProtonInelasticXS.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4DynamicParticle.hh"
45 #include "G4ParticleDefinition.hh"
46 #include "G4Proton.hh"
47 #include "G4Log.hh"
48 #include "G4Exp.hh"
49 #include "G4Pow.hh"
50 
51 
52 // factory
53 #include "G4CrossSectionFactory.hh"
54 //
55 G4_DECLARE_XS_FACTORY(G4ChipsProtonInelasticXS);
56 
57 G4ChipsProtonInelasticXS::G4ChipsProtonInelasticXS():G4VCrossSectionDataSet(Default_Name())
58 {
59  // Initialization of the
60  lastLEN=0; // Pointer to the lastArray of LowEn CS
61  lastHEN=0; // Pointer to the lastArray of HighEn CS
62  lastN=0; // The last N of calculated nucleus
63  lastZ=0; // The last Z of calculated nucleus
64  lastP=0.; // Last used in cross section Momentum
65  lastTH=0.; // Last threshold momentum
66  lastCS=0.; // Last value of the Cross Section
67  lastI=0; // The last position in the DAMDB
68 
69  LEN = new std::vector<G4double*>;
70  HEN = new std::vector<G4double*>;
71 }
72 
73 G4ChipsProtonInelasticXS::~G4ChipsProtonInelasticXS()
74 {
75  G4int lens=LEN->size();
76  for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
77  delete LEN;
78  G4int hens=HEN->size();
79  for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
80  delete HEN;
81 }
82 
83 void
84 G4ChipsProtonInelasticXS::CrossSectionDescription(std::ostream& outFile) const
85 {
86  outFile << "G4ChipsProtonInelasticXS provides the inelastic cross\n"
87  << "section for proton nucleus scattering as a function of incident\n"
88  << "momentum. The cross section is calculated using M. Kossov's\n"
89  << "CHIPS parameterization of cross section data.\n";
90 }
91 
92 G4bool G4ChipsProtonInelasticXS::IsIsoApplicable(const G4DynamicParticle*, G4int, G4int,
93  const G4Element*,
94  const G4Material*)
95 {
96  return true;
97 }
98 
99 
100 // The main member function giving the collision cross section (P is in IU, CS is in mb)
101 // Make pMom in independent units ! (Now it is MeV)
102 G4double G4ChipsProtonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
103  const G4Isotope*,
104  const G4Element*,
105  const G4Material*)
106 {
107  G4double pMom=Pt->GetTotalMomentum();
108  G4int tgN = A - tgZ;
109 
110  return GetChipsCrossSection(pMom, tgZ, tgN, 2212);
111 }
112 
113 G4double G4ChipsProtonInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
114 {
115 
116  G4bool in=false; // By default the isotope must be found in the AMDB
117  if(tgN!=lastN || tgZ!=lastZ) // The nucleus was not the last used isotope
118  {
119  in = false; // By default the isotope haven't been found in AMDB
120  lastP = 0.; // New momentum history (nothing to compare with)
121  lastN = tgN; // The last N of the calculated nucleus
122  lastZ = tgZ; // The last Z of the calculated nucleus
123  lastI = colN.size(); // Size of the Associative Memory DB in the heap
124  j = 0; // A#0f records found in DB for this projectile
125  if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope
126  {
127  if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
128  {
129  lastI=i; // Remember the index for future fast/last use
130  lastTH =colTH[i]; // The last THreshold (A-dependent)
131  if(pMom<=lastTH)
132  {
133  return 0.; // Energy is below the Threshold value
134  }
135  lastP =colP [i]; // Last Momentum (A-dependent)
136  lastCS =colCS[i]; // Last CrossSect (A-dependent)
137  in = true; // This is the case when the isotop is found in DB
138  // Momentum pMom is in IU ! @@ Units
139  lastCS=CalculateCrossSection(-1,j,2212,lastZ,lastN,pMom); // read & update
140  if(lastCS<=0. && pMom>lastTH) // Correct the threshold (@@ No intermediate Zeros)
141  {
142  lastCS=0.;
143  lastTH=pMom;
144  }
145  break; // Go out of the LOOP
146  }
147  j++; // Increment a#0f records found in DB
148  }
149  if(!in) // This isotope has not been calculated previously
150  {
152  lastCS=CalculateCrossSection(0,j,2212,lastZ,lastN,pMom); //calculate & create
153  //if(lastCS>0.) // It means that the AMBD was initialized
154  //{
155 
156  lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
157  colN.push_back(tgN);
158  colZ.push_back(tgZ);
159  colP.push_back(pMom);
160  colTH.push_back(lastTH);
161  colCS.push_back(lastCS);
162  //} // M.K. Presence of H1 with high threshold breaks the syncronization
163  return lastCS*millibarn;
164  } // End of creation of the new set of parameters
165  else
166  {
167  colP[lastI]=pMom;
168  colCS[lastI]=lastCS;
169  }
170  } // End of parameters udate
171  else if(pMom<=lastTH)
172  {
173  return 0.; // Momentum is below the Threshold Value -> CS=0
174  }
175  else // It is the last used -> use the current tables
176  {
177  lastCS=CalculateCrossSection(1,j,2212,lastZ,lastN,pMom); // Only read and UpdateDB
178  lastP=pMom;
179  }
180  return lastCS*millibarn;
181 }
182 
183 // The main member function giving the gamma-A cross section (E in GeV, CS in mb)
184 G4double G4ChipsProtonInelasticXS::CalculateCrossSection(G4int F, G4int I,
185  G4int, G4int targZ, G4int targN, G4double Momentum)
186 {
187  static const G4double THmin=27.; // default minimum Momentum (MeV/c) Threshold
188  static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
189  static const G4double dP=10.; // step for the LEN (Low ENergy) table MeV/c
190  static const G4double dPG=dP*.001; // step for the LEN (Low ENergy) table GeV/c
191  static const G4int nL=105; // A#of LEN points in E (step 10 MeV/c)
192  static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
193  static const G4double Pmax=227000.; // maxP for the HEN (High ENergy) part 227 GeV
194  static const G4int nH=224; // A#of HEN points in lnE
195  static const G4double milP=G4Log(Pmin);// Low logarithm energy for the HEN part
196  static const G4double malP=G4Log(Pmax);// High logarithm energy (each 2.75 percent)
197  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
198  static const G4double milPG=G4Log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
199  G4double sigma=0.;
200  if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
201  //G4double A=targN+targZ; // A of the target
202  if(F<=0) // This isotope was not the last used isotop
203  {
204  if(F<0) // This isotope was found in DAMDB =-----=> RETRIEVE
205  {
206  G4int sync=LEN->size();
207  if(sync<=I) G4cout<<"*!*G4QProtonNuclCS::CalcCrossSect:Sync="<<sync<<"<="<<I<<G4endl;
208  lastLEN=(*LEN)[I]; // Pointer to prepared LowEnergy cross sections
209  lastHEN=(*HEN)[I]; // Pointer to prepared High Energy cross sections
210  }
211  else // This isotope wasn't calculated before => CREATE
212  {
213  lastLEN = new G4double[nL]; // Allocate memory for the new LEN cross sections
214  lastHEN = new G4double[nH]; // Allocate memory for the new HEN cross sections
215  // --- Instead of making a separate function ---
216  G4double P=THmiG; // Table threshold in GeV/c
217  for(G4int k=0; k<nL; k++)
218  {
219  lastLEN[k] = CrossSectionLin(targZ, targN, P);
220  P+=dPG;
221  }
222  G4double lP=milPG;
223  for(G4int n=0; n<nH; n++)
224  {
225  lastHEN[n] = CrossSectionLog(targZ, targN, lP);
226  lP+=dlP;
227  }
228  // --- End of possible separate function
229  // *** The synchronization check ***
230  G4int sync=LEN->size();
231  if(sync!=I)
232  {
233  G4cout<<"***G4ChipsProtonNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
234  <<", N="<<targN<<", F="<<F<<G4endl;
235  //G4Exception("G4ProtonNuclearCS::CalculateCS:","39",FatalException,"overflow DB");
236  }
237  LEN->push_back(lastLEN); // remember the Low Energy Table
238  HEN->push_back(lastHEN); // remember the High Energy Table
239  } // End of creation of the new set of parameters
240  } // End of parameters udate
241  // =------------------= NOW the Magic Formula =-----------------------=
242  if (Momentum<lastTH) return 0.; // It must be already checked in the interface class
243  else if (Momentum<Pmin) // High Energy region
244  {
245  sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
246  }
247  else if (Momentum<Pmax) // High Energy region
248  {
249  G4double lP=G4Log(Momentum);
250  sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
251  }
252  else // UHE region (calculation, not frequent)
253  {
254  G4double P=0.001*Momentum; // Approximation formula is for P in GeV/c
255  sigma=CrossSectionFormula(targZ, targN, P, G4Log(P));
256  }
257  if(sigma<0.) return 0.;
258  return sigma;
259 }
260 
261 // Electromagnetic momentum-threshold (in MeV/c)
262 G4double G4ChipsProtonInelasticXS::ThresholdMomentum(G4int tZ, G4int tN)
263 {
264  static const G4double third=1./3.;
265  static const G4double pM = G4Proton::Proton()->Definition()->GetPDGMass(); // Projectile mass in MeV
266  static const G4double tpM= pM+pM; // Doubled projectile mass (MeV)
267 
268  G4double tA=tZ+tN;
269  if(tZ<.99 || tN<0.) return 0.;
270  else if(tZ==1 && tN==0) return 800.; // A threshold on the free proton
271  //G4double dE=1.263*tZ/(1.+G4Pow::GetInstance()->powA(tA,third));
272  G4double dE=tZ/(1.+G4Pow::GetInstance()->powA(tA,third)); // Safety for diffused edge of the nucleus (QE)
273  G4double tM=931.5*tA;
274  G4double T=dE+dE*(dE/2+pM)/tM;
275  return std::sqrt(T*(tpM+T));
276 }
277 
278 // Calculation formula for proton-nuclear inelastic cross-section (mb) (P in GeV/c)
279 G4double G4ChipsProtonInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
280 {
281  G4double sigma=0.;
282  if(P<ThresholdMomentum(tZ,tN)*.001) return sigma;
283  G4double lP=G4Log(P);
284  if(tZ==1&&!tN){if(P>.35) sigma=CrossSectionFormula(tZ,tN,P,lP);}// s(pp)=0 below 350Mev/c
285  else if(tZ<97 && tN<152) // General solution
286  {
287  G4double pex=0.;
288  G4double pos=0.;
289  G4double wid=1.;
290  if(tZ==13 && tN==14) // Excited metastable states
291  {
292  pex=230.;
293  pos=.13;
294  wid=8.e-5;
295  }
296  else if(tZ<7)
297  {
298  if(tZ==6 && tN==6)
299  {
300  pex=320.;
301  pos=.14;
302  wid=7.e-6;
303  }
304  else if(tZ==5 && tN==6)
305  {
306  pex=270.;
307  pos=.17;
308  wid=.002;
309  }
310  else if(tZ==4 && tN==5)
311  {
312  pex=600.;
313  pos=.132;
314  wid=.005;
315  }
316  else if(tZ==3 && tN==4)
317  {
318  pex=280.;
319  pos=.19;
320  wid=.0025;
321  }
322  else if(tZ==3 && tN==3)
323  {
324  pex=370.;
325  pos=.171;
326  wid=.006;
327  }
328  else if(tZ==2 && tN==1)
329  {
330  pex=30.;
331  pos=.22;
332  wid=.0005;
333  }
334  }
335  sigma=CrossSectionFormula(tZ,tN,P,lP);
336  if(pex>0.)
337  {
338  G4double dp=P-pos;
339  sigma+=pex*G4Exp(-dp*dp/wid);
340  }
341  }
342  else
343  {
344  G4cerr<<"-Warning-G4ChipsProtonNuclearXS::CSLin:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
345  sigma=0.;
346  }
347  if(sigma<0.) return 0.;
348  return sigma;
349 }
350 
351 // Calculation formula for proton-nuclear inelastic cross-section (mb) log(P in GeV/c)
352 G4double G4ChipsProtonInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
353 {
354  G4double P=G4Exp(lP);
355  return CrossSectionFormula(tZ, tN, P, lP);
356 }
357 // Calculation formula for proton-nuclear inelastic cross-section (mb) log(P in GeV/c)
358 G4double G4ChipsProtonInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
359  G4double P, G4double lP)
360 {
361  G4double sigma=0.;
362  if(tZ==1 && !tN) // pp interaction (from G4QuasiElasticRatios)
363  {
364  G4double El(0.),To(0.); // Uzhi
365  if(P<0.1) // Copied from G4QuasiElasticRatios Uzhi / start
366  {
367  G4double p2=P*P;
368  El=1./(0.00012+p2*0.2);
369  To=El;
370  }
371  else if(P>1000.)
372  {
373  G4double lp=G4Log(P)-3.5;
374  G4double lp2=lp*lp;
375  El=0.0557*lp2+6.72;
376  To=0.3*lp2+38.2;
377  }
378  else
379  {
380  G4double p2=P*P;
381  G4double LE=1./(0.00012+p2*0.2);
382  G4double lp=G4Log(P)-3.5;
383  G4double lp2=lp*lp;
384  G4double rp2=1./p2;
385  El=LE+(0.0557*lp2+6.72+32.6/P)/(1.+rp2/P);
386  To=LE+(0.3 *lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
387  } // Copied from G4QuasiElasticRatios Uzhi / end
388 
389 /* // Uzhi 4.03.2013
390  G4double p2=P*P;
391  G4double lp=lP-3.5;
392  G4double lp2=lp*lp;
393  G4double rp2=1./p2;
394  G4double El=(.0557*lp2+6.72+30./P)/(1.+.49*rp2/P);
395  G4double To=(.3*lp2+38.2)/(1.+.54*rp2*rp2);
396 */ // Uzhi 4.03.2013
397 
398  sigma=To-El;
399  }
400  else if(tZ<97 && tN<152) // General solution
401  {
402  //G4double lP=G4Log(P); // Already calculated
403  G4double d=lP-4.2;
404  G4double p2=P*P;
405  G4double p4=p2*p2;
406  G4double a=tN+tZ; // A of the target
407  G4double al=G4Log(a);
408  G4double sa=std::sqrt(a);
409  G4double a2=a*a;
410  G4double a2s=a2*sa;
411  G4double a4=a2*a2;
412  G4double a8=a4*a4;
413  G4double a12=a8*a4;
414  G4double a16=a8*a8;
415  G4double c=(170.+3600./a2s)/(1.+65./a2s);
416  G4double dl=al-3.;
417  G4double dl2=dl*dl;
418  G4double r=.21+.62*dl2/(1.+.5*dl2);
419  G4double gg=40.*G4Exp(al*0.712)/(1.+12.2/a)/(1.+34./a2);
420  G4double e=318.+a4/(1.+.0015*a4/G4Exp(al*0.09))/(1.+4.e-28*a12)+
421  8.e-18/(1./a16+1.3e-20)/(1.+1.e-21*a12);
422  G4double ss=3.57+.009*a2/(1.+.0001*a2*a);
423  G4double h=(.01/a4+2.5e-6/a)*(1.+6.e-6*a2*a)/(1.+6.e7/a12/a2);
424  sigma=(c+d*d)/(1.+r/p4)+(gg+e*G4Exp(-ss*P))/(1.+h/p4/p4);
425  }
426  else
427  {
428  G4cerr<<"-Warning-G4QProtonNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
429  sigma=0.;
430  }
431  if(sigma<0.) return 0.;
432  return sigma;
433 }
434 
435 G4double G4ChipsProtonInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
436 {
437  if(DX<=0. || N<2)
438  {
439  G4cerr<<"***G4ChipsProtonInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
440  return Y[0];
441  }
442 
443  G4int N2=N-2;
444  G4double d=(X-X0)/DX;
445  G4int jj=static_cast<int>(d);
446  if (jj<0) jj=0;
447  else if(jj>N2) jj=N2;
448  d-=jj; // excess
449  G4double yi=Y[jj];
450  G4double sigma=yi+(Y[jj+1]-yi)*d;
451 
452  return sigma;
453 }