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