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