ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChipsPionPlusInelasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ChipsPionPlusInelasticXS.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: G4ChipsPionPlusInelasticXS 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 // Short description: Cross-sections extracted (by W.Pokorski) from the CHIPS package for
36 // pion interactions. Original author: M. Kossov
37 // -------------------------------------------------------------------------------------
38 //
39 
40 #include "G4ChipsPionPlusInelasticXS.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4DynamicParticle.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4PionPlus.hh"
45 
46 #include "G4Log.hh"
47 #include "G4Exp.hh"
48 #include "G4Pow.hh"
49 
50 // factory
51 #include "G4CrossSectionFactory.hh"
52 //
53 G4_DECLARE_XS_FACTORY(G4ChipsPionPlusInelasticXS);
54 
55 G4ChipsPionPlusInelasticXS::G4ChipsPionPlusInelasticXS():G4VCrossSectionDataSet(Default_Name())
56 {
57  // Initialization of the
58  lastLEN=0; // Pointer to lastArray of LowEn CS
59  lastHEN=0; // Pointer to lastArray of HighEn CS
60  lastN=0; // The last N of calculated nucleus
61  lastZ=0; // The last Z of calculated nucleus
62  lastP=0.; // Last used in cross section Momentum
63  lastTH=0.; // Last threshold momentum
64  lastCS=0.; // Last value of the Cross Section
65  lastI=0; // The last position in the DAMDB
66  LEN = new std::vector<G4double*>;
67  HEN = new std::vector<G4double*>;
68 }
69 
70 
71 G4ChipsPionPlusInelasticXS::~G4ChipsPionPlusInelasticXS()
72 {
73  G4int lens=LEN->size();
74  for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
75  delete LEN;
76  G4int hens=HEN->size();
77  for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
78  delete HEN;
79 }
80 
81 void
82 G4ChipsPionPlusInelasticXS::CrossSectionDescription(std::ostream& outFile) const
83 {
84  outFile << "G4ChipsPionPlusInelasticXS provides the inelastic cross\n"
85  << "section for pion+ nucleus scattering as a function of incident\n"
86  << "momentum. The cross section is calculated using M. Kossov's\n"
87  << "CHIPS parameterization of cross section data.\n";
88 }
89 
90 G4bool G4ChipsPionPlusInelasticXS::IsIsoApplicable(const G4DynamicParticle*, G4int, G4int,
91  const G4Element*,
92  const G4Material*)
93 {
94  return true;
95 }
96 
97 // The main member function giving the collision cross section (P is in IU, CS is in mb)
98 // Make pMom in independent units ! (Now it is MeV)
99 G4double G4ChipsPionPlusInelasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
100  const G4Isotope*,
101  const G4Element*,
102  const G4Material*)
103 {
104  G4double pMom=Pt->GetTotalMomentum();
105  G4int tgN = A - tgZ;
106 
107  return GetChipsCrossSection(pMom, tgZ, tgN, 211);
108 }
109 
110 
111 G4double G4ChipsPionPlusInelasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
112 {
113 
114  G4bool in=false; // By default the isotope must be found in the AMDB
115  if(tgN!=lastN || tgZ!=lastZ) // The nucleus was not the last used isotope
116  {
117  in = false; // By default the isotope haven't be found in AMDB
118  lastP = 0.; // New momentum history (nothing to compare with)
119  lastN = tgN; // The last N of the calculated nucleus
120  lastZ = tgZ; // The last Z of the calculated nucleus
121  lastI = colN.size(); // Size of the Associative Memory DB in the heap
122  j = 0; // A#0f records found in DB for this projectile
123  if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope
124  {
125  if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
126  {
127  lastI=i; // Remember the index for future fast/last use
128  lastTH =colTH[i]; // The last THreshold (A-dependent)
129  if(pMom<=lastTH)
130  {
131  return 0.; // Energy is below the Threshold value
132  }
133  lastP =colP [i]; // Last Momentum (A-dependent)
134  lastCS =colCS[i]; // Last CrossSect (A-dependent)
135  in = true; // This is the case when the isotop is found in DB
136  // Momentum pMom is in IU ! @@ Units
137  lastCS=CalculateCrossSection(-1,j,211,lastZ,lastN,pMom); // read & update
138  if(lastCS<=0. && pMom>lastTH) // Correct the threshold (@@ No intermediate Zeros)
139  {
140  lastCS=0.;
141  lastTH=pMom;
142  }
143  break; // Go out of the LOOP
144  }
145  j++; // Increment a#0f records found in DB
146  }
147  if(!in) // This isotope has not been calculated previously
148  {
150  lastCS=CalculateCrossSection(0,j,211,lastZ,lastN,pMom); //calculate & create
151  //if(lastCS>0.) // It means that the AMBD was initialized
152  //{
153 
154  lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
155  colN.push_back(tgN);
156  colZ.push_back(tgZ);
157  colP.push_back(pMom);
158  colTH.push_back(lastTH);
159  colCS.push_back(lastCS);
160  //} // M.K. Presence of H1 with high threshold breaks the syncronization
161  return lastCS*millibarn;
162  } // End of creation of the new set of parameters
163  else
164  {
165  colP[lastI]=pMom;
166  colCS[lastI]=lastCS;
167  }
168  } // End of parameters udate
169  else if(pMom<=lastTH)
170  {
171  return 0.; // Momentum is below the Threshold Value -> CS=0
172  }
173  else // It is the last used -> use the current tables
174  {
175  lastCS=CalculateCrossSection(1,j,211,lastZ,lastN,pMom); // Only read and UpdateDB
176  lastP=pMom;
177  }
178  return lastCS*millibarn;
179 }
180 
181 // The main member function giving the gamma-A cross section (E in GeV, CS in mb)
182 G4double G4ChipsPionPlusInelasticXS::CalculateCrossSection(G4int F, G4int I,
183  G4int, G4int targZ, G4int targN, G4double Momentum)
184 {
185  static const G4double THmin=27.; // default minimum Momentum (MeV/c) Threshold
186  static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
187  static const G4double dP=10.; // step for the LEN (Low ENergy) table MeV/c
188  static const G4double dPG=dP*.001; // step for the LEN (Low ENergy) table GeV/c
189  static const G4int nL=105; // A#of LEN points in E (step 10 MeV/c)
190  static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
191  static const G4double Pmax=227000.; // maxP for the HEN (High ENergy) part 227 GeV
192  static const G4int nH=224; // A#of HEN points in lnE
193  static const G4double milP=G4Log(Pmin);// Low logarithm energy for the HEN part
194  static const G4double malP=G4Log(Pmax);// High logarithm energy (each 2.75 percent)
195  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
196  static const G4double milPG=G4Log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
197  G4double sigma=0.;
198  if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
199  //G4double A=targN+targZ; // A of the target
200  if(F<=0) // This isotope was not the last used isotop
201  {
202  if(F<0) // This isotope was found in DAMDB =-----=> RETRIEVE
203  {
204  G4int sync=LEN->size();
205  if(sync<=I) G4cerr<<"*!*G4ChipsPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
206  lastLEN=(*LEN)[I]; // Pointer to prepared LowEnergy cross sections
207  lastHEN=(*HEN)[I]; // Pointer to prepared High Energy cross sections
208  }
209  else // This isotope wasn't calculated before => CREATE
210  {
211  lastLEN = new G4double[nL]; // Allocate memory for the new LEN cross sections
212  lastHEN = new G4double[nH]; // Allocate memory for the new HEN cross sections
213  // --- Instead of making a separate function ---
214  G4double P=THmiG; // Table threshold in GeV/c
215  for(G4int k=0; k<nL; k++)
216  {
217  lastLEN[k] = CrossSectionLin(targZ, targN, P);
218  P+=dPG;
219  }
220  G4double lP=milPG;
221  for(G4int n=0; n<nH; n++)
222  {
223  lastHEN[n] = CrossSectionLog(targZ, targN, lP);
224  lP+=dlP;
225  }
226  // --- End of possible separate function
227  // *** The synchronization check ***
228  G4int sync=LEN->size();
229  if(sync!=I)
230  {
231  G4cerr<<"***G4ChipsPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
232  <<", N="<<targN<<", F="<<F<<G4endl;
233  //G4Exception("G4PiMinusNuclearCS::CalculateCS:","39",FatalException,"DBoverflow");
234  }
235  LEN->push_back(lastLEN); // remember the Low Energy Table
236  HEN->push_back(lastHEN); // remember the High Energy Table
237  } // End of creation of the new set of parameters
238  } // End of parameters udate
239  // =-----------------= NOW the Magic Formula =-------------------------=
240  if (Momentum<lastTH) return 0.; // It must be already checked in the interface class
241  else if (Momentum<Pmin) // High Energy region
242  {
243  sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
244  }
245  else if (Momentum<Pmax) // High Energy region
246  {
247  G4double lP=G4Log(Momentum);
248  sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
249  }
250  else // UHE region (calculation, not frequent)
251  {
252  G4double P=0.001*Momentum; // Approximation formula is for P in GeV/c
253  sigma=CrossSectionFormula(targZ, targN, P, G4Log(P));
254  }
255  if(sigma<0.) return 0.;
256  return sigma;
257 }
258 
259 // Electromagnetic momentum-threshold (in MeV/c)
260 G4double G4ChipsPionPlusInelasticXS::ThresholdMomentum(G4int tZ, G4int tN)
261 {
262  static const G4double third=1./3.;
263  static const G4double pM = G4PionPlus::PionPlus()->Definition()->GetPDGMass(); // Projectile mass in MeV
264  static const G4double tpM= pM+pM; // Doubled projectile mass (MeV)
265  G4double tA=tZ+tN;
266  if(tZ<.99 || tN<0.) return 0.;
267  else if(tZ==1 && tN==0) return 300.; // A threshold on the free proton
268  //G4double dE=1.263*tZ/(1.+G4Pow::GetInstance()->powA(tA,third));
269  G4double dE=tZ/(1.+G4Pow::GetInstance()->powA(tA,third)); // Safety for diffused edge of the nucleus (QE)
270  G4double tM=931.5*tA;
271  G4double T=dE+dE*(dE/2+pM)/tM;
272  return std::sqrt(T*(tpM+T));
273 }
274 
275 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) (P in GeV/c)
276 G4double G4ChipsPionPlusInelasticXS::CrossSectionLin(G4int tZ, G4int tN, G4double P)
277 {
278  G4double lP=G4Log(P);
279  return CrossSectionFormula(tZ, tN, P, lP);
280 }
281 
282 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
283 G4double G4ChipsPionPlusInelasticXS::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
284 {
285  G4double P=G4Exp(lP);
286  return CrossSectionFormula(tZ, tN, P, lP);
287 }
288 // Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
289 G4double G4ChipsPionPlusInelasticXS::CrossSectionFormula(G4int tZ, G4int tN,
290  G4double P, G4double lP)
291 {
292  G4double sigma=0.;
293  if(tZ==1 && !tN) // PiPlus-Proton interaction from G4QuasiElRatios
294  {
295  G4double ld=lP-3.5;
296  G4double ld2=ld*ld;
297  G4double p2=P*P;
298  G4double p4=p2*p2;
299  G4double sp=std::sqrt(P);
300  G4double lm=lP-.32;
301  G4double md=lm*lm+.04;
302  G4double El=(.0557*ld2+2.4+6./sp)/(1.+3./p4);
303  G4double To=(.3*ld2+22.3+5./sp)/(1.+1./p4);
304  sigma=(To-El)+.1/md;
305  }
306  else if(tZ==1 && tN==1) // pimp_tot
307  {
308  G4double p2=P*P;
309  G4double d=lP-2.7;
310  G4double f=lP+1.25;
311  G4double gg=lP-.017;
312  sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+.3/p2/p2)+18./(f*f+.1089)+.02/(gg*gg+.0025);
313  }
314  else if(tZ<97 && tN<152) // General solution
315  {
316  G4double d=lP-4.2;
317  G4double p2=P*P;
318  G4double p4=p2*p2;
319  G4double a=tN+tZ; // A of the target
320  G4double al=G4Log(a);
321  G4double sa=std::sqrt(a);
322  G4double ssa=std::sqrt(sa);
323  G4double a2=a*a;
324  G4double c=41.*G4Exp(al*.68)*(1.+44./a2)/(1.+8./a)/(1.+200./a2/a2);
325  G4double f=290.*ssa/(1.+34./a/ssa);
326  G4double gg=-1.32-al*.043;
327  G4double u=lP-gg;
328  G4double h=al*(.4-.055*al);
329  G4double r=.01+a2*5.E-8;
330  sigma=(c+d*d)/(1.+(.2-.009*sa)/p4)+f/(u*u+h*h)/(1.+r/p2);
331  }
332  else
333  {
334  G4cerr<<"-Warning-G4ChipsPiPlusNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
335  sigma=0.;
336  }
337  if(sigma<0.) return 0.;
338  return sigma;
339 }
340 
341 G4double G4ChipsPionPlusInelasticXS::EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double* Y)
342 {
343  if(DX<=0. || N<2)
344  {
345  G4cerr<<"***G4ChipsPionPlusInelasticXS::EquLinearFit: DX="<<DX<<", N="<<N<<G4endl;
346  return Y[0];
347  }
348 
349  G4int N2=N-2;
350  G4double d=(X-X0)/DX;
351  G4int jj=static_cast<int>(d);
352  if (jj<0) jj=0;
353  else if(jj>N2) jj=N2;
354  d-=jj; // excess
355  G4double yi=Y[jj];
356  G4double sigma=yi+(Y[jj+1]-yi)*d;
357 
358  return sigma;
359 }