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