Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4QKaonPlusNuclearCrossSection.cc
Go to the documentation of this file.
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// GEANT4 tag $Name: not supported by cvs2svn $
29//
30//
31// G4 Physics class: G4QKaonPlusNuclearCrossSection for gamma+A cross sections
32// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03
33// The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04
34// --------------------------------------------------------------------------------
35// ****************************************************************************************
36// This Header is a part of the CHIPS physics package (author: M. Kosov)
37// ****************************************************************************************
38// Short description: CHIPS cross-sections for kaon(plus)-nuclear interactions
39// -------------------------------------------------------------------------------------
40//
41//#define debug
42//#define pdebug
43//#define debug3
44//#define debugn
45//#define debugs
46
48#include "G4SystemOfUnits.hh"
49
50// Initialization of the
51G4double* G4QKaonPlusNuclearCrossSection::lastLEN=0; // Pointer to the lastArray of LowEn CS
52G4double* G4QKaonPlusNuclearCrossSection::lastHEN=0; // Pointer to the lastArray of HighEn CS
53G4int G4QKaonPlusNuclearCrossSection::lastN=0; // The last N of calculated nucleus
54G4int G4QKaonPlusNuclearCrossSection::lastZ=0; // The last Z of calculated nucleus
55G4double G4QKaonPlusNuclearCrossSection::lastP=0.; // Last used in cross section Momentum
56G4double G4QKaonPlusNuclearCrossSection::lastTH=0.; // Last threshold momentum
57G4double G4QKaonPlusNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
58G4int G4QKaonPlusNuclearCrossSection::lastI=0; // The last position in the DAMDB
59std::vector<G4double*>* G4QKaonPlusNuclearCrossSection::LEN = new std::vector<G4double*>;
60std::vector<G4double*>* G4QKaonPlusNuclearCrossSection::HEN = new std::vector<G4double*>;
61
62// Returns Pointer to the G4VQCrossSection class
64{
65 static G4QKaonPlusNuclearCrossSection theCrossSection; //**Static body of Cross Section**
66 return &theCrossSection;
67}
68
70{
71 G4int lens=LEN->size();
72 for(G4int i=0; i<lens; ++i) delete[] (*LEN)[i];
73 delete LEN;
74 G4int hens=HEN->size();
75 for(G4int i=0; i<hens; ++i) delete[] (*HEN)[i];
76 delete HEN;
77}
78
79// The main member function giving the collision cross section (P is in IU, CS is in mb)
80// Make pMom in independent units ! (Now it is MeV)
82 G4int tgZ, G4int tgN, G4int PDG)
83{
84 //A.R.23-Oct-2012 Shadowed variable static G4double tolerance=0.001; // Tolerance (0.1%) to consider as "the same mom"
85 static G4int j; // A#0f Z/N-records already tested in AMDB
86 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops)
87 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops)
88 static std::vector <G4double> colP; // Vector of last momenta for the reaction
89 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
90 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
91 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
92#ifdef debug
93 G4cout<<"G4QKpCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
94 <<"("<<lastN<<"),PDG="<<PDG<<", thresh="<<lastTH<<",Sz="<<colN.size()<<G4endl;
95#endif
96 if(PDG!=321) G4cout<<"-Warning-G4QKaonPlusCS::GetCS:***Not a KPlus***,PDG="<<PDG<<G4endl;
97 G4bool in=false; // By default the isotope must be found in the AMDB
98 if(tgN!=lastN || tgZ!=lastZ) // The nucleus was not the last used isotope
99 {
100 in = false; // By default the isotope haven't be found in AMDB
101 lastP = 0.; // New momentum history (nothing to compare with)
102 lastN = tgN; // The last N of the calculated nucleus
103 lastZ = tgZ; // The last Z of the calculated nucleus
104 lastI = colN.size(); // Size of the Associative Memory DB in the heap
105 j = 0; // A#0f records found in DB for this projectile
106#ifdef debug
107 G4cout<<"G4QKpCS::GetCS: the amount of records in the AMDB lastI="<<lastI<<G4endl;
108#endif
109 if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope
110 {
111 if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
112 {
113 lastI=i; // Remember the index for future fast/last use
114 lastTH =colTH[i]; // The last THreshold (A-dependent)
115#ifdef debug
116 G4cout<<"G4QKpCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
117#endif
118 if(pMom<=lastTH)
119 {
120#ifdef debug
121 G4cout<<"G4QPCS::GetCS:Found,P="<<pMom<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
122#endif
123 return 0.; // Energy is below the Threshold value
124 }
125 lastP =colP [i]; // Last Momentum (A-dependent)
126 lastCS =colCS[i]; // Last CrossSect (A-dependent)
127 if(std::fabs(lastP-pMom)<tolerance*pMom)
128 //if(lastP==pMom) // VI do not use tolerance
129 {
130#ifdef debug
131 G4cout<<"..G4QKpCS::GetCS:.DoNothing.P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
132#endif
133 //CalculateCrossSection(fCS,-1,j,321,lastZ,lastN,pMom); // Update param's only
134 return lastCS*millibarn; // Use theLastCS
135 }
136 in = true; // This is the case when the isotop is found in DB
137 // Momentum pMom is in IU ! @@ Units
138#ifdef debug
139 G4cout<<"G4QKpCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
140#endif
141 lastCS=CalculateCrossSection(fCS,-1,j,321,lastZ,lastN,pMom); // read & update
142#ifdef debug
143 G4cout<<"G4QKpCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
144#endif
145 if(lastCS<=0. && pMom>lastTH) // Correct the threshold (@@ No intermediate Zeros)
146 {
147#ifdef debug
148 G4cout<<"G4QKpCS::GetCS: New P="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
149#endif
150 lastCS=0.;
151 lastTH=pMom;
152 }
153 break; // Go out of the LOOP
154 }
155#ifdef debug
156 G4cout<<"-->G4QKpCrossSec::GetCrosSec: pPDG=321, j="<<j<<", N="<<colN[i]
157 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
158#endif
159 j++; // Increment a#0f records found in DB
160 }
161#ifdef debug
162 G4cout<<"-?-G4QKpCS::GetCS:RC Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl;
163#endif
164 if(!in) // This isotope has not been calculated previously
165 {
166#ifdef debug
167 G4cout<<"^^^G4QKpCS::GetCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
168#endif
169 //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
170 lastCS=CalculateCrossSection(fCS,0,j,321,lastZ,lastN,pMom); //calculate & create
171 //if(lastCS>0.) // It means that the AMBD was initialized
172 //{
173
174 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
175#ifdef debug
176 G4cout<<"G4QKpCrossSection::GetCrossSect: NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
177#endif
178 colN.push_back(tgN);
179 colZ.push_back(tgZ);
180 colP.push_back(pMom);
181 colTH.push_back(lastTH);
182 colCS.push_back(lastCS);
183#ifdef debug
184 G4cout<<"G4QKpCS::GetCrosSec:recCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
185#endif
186 //} // M.K. Presence of H1 with high threshold breaks the syncronization
187#ifdef pdebug
188 G4cout<<"G4QKpCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
189#endif
190 return lastCS*millibarn;
191 } // End of creation of the new set of parameters
192 else
193 {
194#ifdef debug
195 G4cout<<"G4QKpCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
196#endif
197 colP[lastI]=pMom;
198 colCS[lastI]=lastCS;
199 }
200 } // End of parameters udate
201 else if(pMom<=lastTH)
202 {
203#ifdef debug
204 G4cout<<"G4QKpCS::GetCS: Current P="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
205#endif
206 return 0.; // Momentum is below the Threshold Value -> CS=0
207 }
208 else if(std::fabs(lastP-pMom)<tolerance*pMom)
209 //else if(lastP==pMom) // VI do not use tolerance
210 {
211#ifdef debug
212 G4cout<<"..G4QPCS::GetCS:OldNZ&P="<<lastP<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
213#endif
214 return lastCS*millibarn; // Use theLastCS
215 }
216 else // It is the last used -> use the current tables
217 {
218#ifdef debug
219 G4cout<<"-!-G4QPCS::GetCS:UseCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
220#endif
221 lastCS=CalculateCrossSection(fCS,1,j,321,lastZ,lastN,pMom); // Only read and UpdateDB
222 lastP=pMom;
223 }
224#ifdef debug
225 G4cout<<"==>G4QKpCS::GetCroSec: P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
226#endif
227 return lastCS*millibarn;
228}
229
230// The main member function giving the gamma-A cross section (E in GeV, CS in mb)
232 G4int, G4int targZ, G4int targN, G4double Momentum)
233{
234 static const G4double THmin=27.; // default minimum Momentum (MeV/c) Threshold
235 static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
236 static const G4double dP=10.; // step for the LEN (Low ENergy) table MeV/c
237 static const G4double dPG=dP*.001; // step for the LEN (Low ENergy) table GeV/c
238 static const G4int nL=105; // A#of LEN points in E (step 10 MeV/c)
239 static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
240 static const G4double Pmax=227000.; // maxP for the HEN (High ENergy) part 227 GeV
241 static const G4int nH=224; // A#of HEN points in lnE
242 static const G4double milP=std::log(Pmin);// Low logarithm energy for the HEN part
243 static const G4double malP=std::log(Pmax);// High logarithm energy (each 2.75 percent)
244 static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
245 static const G4double milPG=std::log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
246#ifdef debug
247 G4cout<<"G4QKpNucCS::CalCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<">"<<THmin<<G4endl;
248#endif
249 G4double sigma=0.;
250 if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
251 G4double A=targN+targZ; // A of the target
252#ifdef debug
253 G4cout<<"G4QKpNucCS::CalCS: A="<<A<<",F="<<F<<",I="<<I<<",nL="<<nL<<",nH="<<nH<<G4endl;
254#endif
255 if(F<=0) // This isotope was not the last used isotop
256 {
257 if(F<0) // This isotope was found in DAMDB =-----=> RETRIEVE
258 {
259 G4int sync=LEN->size();
260 if(sync<=I) G4cerr<<"*!*G4QPiMinusNuclCS::CalcCrosSect:Sync="<<sync<<"<="<<I<<G4endl;
261 lastLEN=(*LEN)[I]; // Pointer to prepared LowEnergy cross sections
262 lastHEN=(*HEN)[I]; // Pointer to prepared High Energy cross sections
263 }
264 else // This isotope wasn't calculated before => CREATE
265 {
266 lastLEN = new G4double[nL]; // Allocate memory for the new LEN cross sections
267 lastHEN = new G4double[nH]; // Allocate memory for the new HEN cross sections
268 // --- Instead of making a separate function ---
269 G4double P=THmiG; // Table threshold in GeV/c
270 for(G4int n=0; n<nL; n++)
271 {
272 lastLEN[n] = CrossSectionLin(targZ, targN, P);
273 P+=dPG;
274 }
275 G4double lP=milPG;
276 for(G4int n=0; n<nH; n++)
277 {
278 lastHEN[n] = CrossSectionLog(targZ, targN, lP);
279 lP+=dlP;
280 }
281#ifdef debug
282 G4cout<<"-*->G4QKpNucCS::CalcCS:Tab for Z="<<targZ<<",N="<<targN<<",I="<<I<<G4endl;
283#endif
284 // --- End of possible separate function
285 // *** The synchronization check ***
286 G4int sync=LEN->size();
287 if(sync!=I)
288 {
289 G4cerr<<"***G4QPiMinusNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
290 <<", N="<<targN<<", F="<<F<<G4endl;
291 //G4Exception("G4PiMinusNuclearCS::CalculateCS:","39",FatalException,"DBoverflow");
292 }
293 LEN->push_back(lastLEN); // remember the Low Energy Table
294 HEN->push_back(lastHEN); // remember the High Energy Table
295 } // End of creation of the new set of parameters
296 } // End of parameters udate
297 // =--------------------------= NOW the Magic Formula =---------------------------------=
298#ifdef debug
299 G4cout<<"G4QKpNCS::CalcCS:lTH="<<lastTH<<",Pmi="<<Pmin<<",dP="<<dP<<",dlP="<<dlP<<G4endl;
300#endif
301 if (Momentum<lastTH) return 0.; // It must be already checked in the interface class
302 else if (Momentum<Pmin) // Low Energy region
303 {
304#ifdef debug
305 G4cout<<"G4QKpNCS::CalcCS:bLEN A="<<A<<", nL="<<nL<<", P="<<Momentum<<G4endl;
306#endif
307 if(A<=1. && Momentum < 600.) sigma=0.; // Approximation tot/el uncertainty
308 else sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
309#ifdef debugn
310 if(sigma<0.)
311 G4cout<<"G4QKpNuCS::CalcCS: E="<<Momentum<<",T="<<THmin<<",dP="<<dP<<G4endl;
312#endif
313 }
314 else if (Momentum<Pmax) // High Energy region
315 {
316 G4double lP=std::log(Momentum);
317#ifdef debug
318 G4cout<<"G4QKpNucCS::CalcCS: before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
319#endif
320 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
321 }
322 else // UHE region (calculation, not frequent)
323 {
324 G4double P=0.001*Momentum; // Approximation formula is for P in GeV/c
325 sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
326 }
327#ifdef debug
328 G4cout<<"G4QKaonPlusNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
329#endif
330 if(sigma<0.) return 0.;
331 return sigma;
332}
333
334// Electromagnetic momentum-threshold (in MeV/c)
335G4double G4QKaonPlusNuclearCrossSection::ThresholdMomentum(G4int tZ, G4int tN)
336{
337 static const G4double third=1./3.;
338 static const G4double prM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
339 static const G4double piM = G4QPDGCode(111).GetMass()+.1; // Pion mass in MeV+Safety
340 static const G4double pM = G4QPDGCode(321).GetMass(); // Projectile mass in MeV
341 static const G4double tpM= pM+pM; // Doubled projectile mass (MeV)
342 G4double tA=tZ+tN;
343 if(tZ<.99 || tN<0.) return 0.;
344 G4double tM=931.5*tA;
345 G4double dE=piM; // At least one Pi0 must be created
346 if(tZ==1 && tN==0) tM=prM; // A threshold on the free proton
347 else dE=tZ/(1.+std::pow(tA,third)); // Safety for diffused edge of the nucleus (QE)
348 //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
349 G4double T=dE+dE*(dE/2+pM)/tM;
350 return std::sqrt(T*(tpM+T));
351}
352
353// Calculation formula for piMinus-nuclear inelastic cross-section (mb) (P in GeV/c)
354G4double G4QKaonPlusNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
355{
356 G4double lP=std::log(P);
357 return CrossSectionFormula(tZ, tN, P, lP);
358}
359
360// Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
361G4double G4QKaonPlusNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
362{
363 G4double P=std::exp(lP);
364 return CrossSectionFormula(tZ, tN, P, lP);
365}
366// Calculation formula for piMinus-nuclear inelastic cross-section (mb) log(P in GeV/c)
367G4double G4QKaonPlusNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
368 G4double P, G4double lP)
369{
370 G4double sigma=0.;
371 if(tZ==1 && !tN) // KPlus-Proton interaction from G4QuasiElRatios
372 {
373 G4double ld=lP-3.5;
374 G4double ld2=ld*ld;
375 G4double sp=std::sqrt(P);
376 G4double p2=P*P;
377 G4double p4=p2*p2;
378 G4double lm=P-1.;
379 G4double md=lm*lm+.372;
380 G4double El=(.0557*ld2+2.23)/(1.-.7/sp+.1/p4);
381 G4double To=(.3*ld2+19.5)/(1.+.46/sp+1.6/p4);
382 sigma=(To-El)+.6/md;
383 }
384 else if(tZ<97 && tN<152) // General solution
385 {
386 G4double p2=P*P;
387 G4double p4=p2*p2;
388 G4double a=tN+tZ; // A of the target
389 G4double al=std::log(a);
390 G4double sa=std::sqrt(a);
391 G4double asa=a*sa;
392 G4double a2=a*a;
393 G4double a3=a2*a;
394 G4double a4=a2*a2;
395 G4double a8=a4*a4;
396 G4double a12=a8*a4;
397 G4double f=.6; // Default values for deutrons
398 G4double r=.5;
399 G4double g_value=3.7;
400 G4double c=36.;
401 G4double s_value=3.5;
402 G4double t=3.;
403 G4double u=.44;
404 G4double v=5.E-9;
405 if(tZ>1 && tN>1) // More than deuteron
406 {
407 f=1.;
408 r=1./(1.+.007*a2);
409 g_value=4.2;
410 c=52.*std::exp(al*.6)*(1.+95./a2)/(1.+9./a)/(1.+46./a2);
411 s_value=(40.+.14*a)/(1.+12./a);
412 G4double y=std::exp(al*1.7);
413 t=.185*y/(1.+.00012*y);
414 u=(1.+80./asa)/(1.+200./asa);
415 v=(1.+3.E-6*a4*(1.+6.E-7*a3+4.E10/a12))/a3/20000.;
416 }
417 G4double d=lP-g_value;
418 G4double w=P-1.;
419 G4double rD=s_value/(w*w+.36);
420 G4double h=P-.44;
421 G4double rR=t/(h*h+u*u);
422 sigma=(f*d*d+c)/(1.+r/std::sqrt(P)+1./p4)+(rD+rR)/(1+v/p4/p4);
423#ifdef pdebug
424 G4cout<<"G4QKaonPlusNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g_value
425 <<",d="<<d<<",r="<<r<<",h="<<h<<G4endl;
426#endif
427 }
428 else
429 {
430 G4cerr<<"-Warning-G4QKaonPlusNuclearCroSect::CSForm:Bad A, Z="<<tZ<<", N="<<tN<<G4endl;
431 sigma=0.;
432 }
433 if(sigma<0.) return 0.;
434 return sigma;
435}
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
virtual G4double GetCrossSection(G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=321)
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4double EquLinearFit(G4double X, G4int N, G4double X0, G4double DX, G4double *Y)
virtual G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=0)
static G4double tolerance