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