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