ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChipsKaonPlusElasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ChipsKaonPlusElasticXS.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 //
28 //
29 // G4 Physics class: G4ChipsKaonPlusElasticXS for pA elastic cross sections
30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 5-Feb-2010
31 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 5-Feb-2010
32 //
33 // -------------------------------------------------------------------------------
34 // Short description: Interaction cross-sections for the elastic process.
35 // Class extracted from CHIPS and integrated in Geant4 by W.Pokorski
36 // -------------------------------------------------------------------------------
37 //
38 
39 #include "G4ChipsKaonPlusElasticXS.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DynamicParticle.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4KaonPlus.hh"
44 #include "G4Nucleus.hh"
45 #include "G4ParticleTable.hh"
46 #include "G4NucleiProperties.hh"
47 #include "G4IonTable.hh"
48 #include "G4AutoLock.hh"
49 // factory
50 #include "G4CrossSectionFactory.hh"
51 //
52 G4_DECLARE_XS_FACTORY(G4ChipsKaonPlusElasticXS);
53 
54 namespace {
55  G4double mK;//= G4KaonPlus::KaonPlus()->GetPDGMass()*.001; // MeV to GeV //Cannot initialize here, needs particles
56  G4double mK2;//= mK*mK;
57  G4Mutex initM = G4MUTEX_INITIALIZER;
58  const G4double pwd=2727;
59  const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
60  const G4double third=1./3.;
61  const G4double fifth=1./5.;
62  const G4double sevth=1./7.;
63  const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
64  // @@ At present all nA==pA ---------> Each neucleus can have not more than 51 parameters
65  const G4int n_kppel=35; // #of parameters for pp-elastic (<nPoints=128)
66  // -0--1- -2- -3- -4- -5- -6--7--8--9- -10--11-12--13--14-
67  G4double kpp_el[n_kppel]={.7,.38,.0676,.0557,3.5,2.23,.7,.1,2.,1.,.372,5.,74.,3.,3.4,
68  .2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,
69  5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
70  // -15--16--17--18--19- -20- -21- -22- -23- -24- -25- -26-
71  // -27- -28- -29- -30- -31- -32- -33--34-
72 
73 
74 }
75 
76 G4ChipsKaonPlusElasticXS::G4ChipsKaonPlusElasticXS():G4VCrossSectionDataSet(Default_Name()), nPoints(128), nLast(nPoints-1)
77 {
78  G4AutoLock l(&initM);
79  mK = G4KaonPlus::KaonPlus()->GetPDGMass()*.001;// MeV to GeV
80  mK2 = mK*mK;
81  l.unlock();
82  lPMin=-8.; //Min tabulatedLogarithmMomentum/D
83  lPMax= 8.; //Max tabulatedLogarithmMomentum/D
84  dlnP=(lPMax-lPMin)/nLast;// LogStep inTable /D
85  onlyCS=true;//Flag toCalculOnlyCS(not Si/Bi)/L
86  lastSIG=0.; //Last calculated cross section /L
87  lastLP=-10.;//LastLog(mom_of IncidentHadron)/L
88  lastTM=0.; //Last t_maximum /L
89  theSS=0.; //TheLastSqSlope of 1st difr.Max/L
90  theS1=0.; //TheLastMantissa of 1st difrMax/L
91  theB1=0.; //TheLastSlope of 1st difructMax/L
92  theS2=0.; //TheLastMantissa of 2nd difrMax/L
93  theB2=0.; //TheLastSlope of 2nd difructMax/L
94  theS3=0.; //TheLastMantissa of 3d difr.Max/L
95  theB3=0.; //TheLastSlope of 3d difruct.Max/L
96  theS4=0.; //TheLastMantissa of 4th difrMax/L
97  theB4=0.; //TheLastSlope of 4th difructMax/L
98  lastTZ=0; // Last atomic number of theTarget
99  lastTN=0; // Last # of neutrons in theTarget
100  lastPIN=0.;// Last initialized max momentum
101  lastCST=0; // Elastic cross-section table
102  lastPAR=0; // ParametersForFunctionCalculation
103  lastSST=0; // E-dep ofSqardSlope of 1st difMax
104  lastS1T=0; // E-dep of mantissa of 1st dif.Max
105  lastB1T=0; // E-dep of the slope of 1st difMax
106  lastS2T=0; // E-dep of mantissa of 2nd difrMax
107  lastB2T=0; // E-dep of the slope of 2nd difMax
108  lastS3T=0; // E-dep of mantissa of 3d difr.Max
109  lastB3T=0; // E-dep of the slope of 3d difrMax
110  lastS4T=0; // E-dep of mantissa of 4th difrMax
111  lastB4T=0; // E-dep of the slope of 4th difMax
112  lastN=0; // The last N of calculated nucleus
113  lastZ=0; // The last Z of calculated nucleus
114  lastP=0.; // LastUsed inCrossSection Momentum
115  lastTH=0.; // Last threshold momentum
116  lastCS=0.; // Last value of the Cross Section
117  lastI=0; // The last position in the DAMDB
118 }
119 
120 G4ChipsKaonPlusElasticXS::~G4ChipsKaonPlusElasticXS()
121 {
122  std::vector<G4double*>::iterator pos;
123  for (pos=CST.begin(); pos<CST.end(); pos++)
124  { delete [] *pos; }
125  CST.clear();
126  for (pos=PAR.begin(); pos<PAR.end(); pos++)
127  { delete [] *pos; }
128  PAR.clear();
129  for (pos=SST.begin(); pos<SST.end(); pos++)
130  { delete [] *pos; }
131  SST.clear();
132  for (pos=S1T.begin(); pos<S1T.end(); pos++)
133  { delete [] *pos; }
134  S1T.clear();
135  for (pos=B1T.begin(); pos<B1T.end(); pos++)
136  { delete [] *pos; }
137  B1T.clear();
138  for (pos=S2T.begin(); pos<S2T.end(); pos++)
139  { delete [] *pos; }
140  S2T.clear();
141  for (pos=B2T.begin(); pos<B2T.end(); pos++)
142  { delete [] *pos; }
143  B2T.clear();
144  for (pos=S3T.begin(); pos<S3T.end(); pos++)
145  { delete [] *pos; }
146  S3T.clear();
147  for (pos=B3T.begin(); pos<B3T.end(); pos++)
148  { delete [] *pos; }
149  B3T.clear();
150  for (pos=S4T.begin(); pos<S4T.end(); pos++)
151  { delete [] *pos; }
152  S4T.clear();
153  for (pos=B4T.begin(); pos<B4T.end(); pos++)
154  { delete [] *pos; }
155  B4T.clear();
156 }
157 
158 void
159 G4ChipsKaonPlusElasticXS::CrossSectionDescription(std::ostream& outFile) const
160 {
161  outFile << "G4ChipsKaonPlusElasticXS provides the elastic cross\n"
162  << "section for K+ nucleus scattering as a function of incident\n"
163  << "momentum. The cross section is calculated using M. Kossov's\n"
164  << "CHIPS parameterization of cross section data.\n";
165 }
166 
167 
168 G4bool G4ChipsKaonPlusElasticXS::IsIsoApplicable(const G4DynamicParticle*, G4int, G4int,
169  const G4Element*,
170  const G4Material*)
171 {
172  return true;
173 }
174 
175 // The main member function giving the collision cross section (P is in IU, CS is in mb)
176 // Make pMom in independent units ! (Now it is MeV)
177 G4double G4ChipsKaonPlusElasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
178  const G4Isotope*,
179  const G4Element*,
180  const G4Material*)
181 {
182  G4double pMom=Pt->GetTotalMomentum();
183  G4int tgN = A - tgZ;
184 
185  return GetChipsCrossSection(pMom, tgZ, tgN, 321);
186 }
187 
188 G4double G4ChipsKaonPlusElasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
189 {
190 
191  G4bool fCS = false;
192  G4double pEn=pMom;
193  onlyCS=fCS;
194 
195  G4bool in=false; // By default the isotope must be found in the AMDB
196  lastP = 0.; // New momentum history (nothing to compare with)
197  lastN = tgN; // The last N of the calculated nucleus
198  lastZ = tgZ; // The last Z of the calculated nucleus
199  lastI = colN.size(); // Size of the Associative Memory DB in the heap
200  if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
201  { // The nucleus with projPDG is found in AMDB
202  if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
203  {
204  lastI=i;
205  lastTH =colTH[i]; // Last THreshold (A-dependent)
206  if(pEn<=lastTH)
207  {
208  return 0.; // Energy is below the Threshold value
209  }
210  lastP =colP [i]; // Last Momentum (A-dependent)
211  lastCS =colCS[i]; // Last CrossSect (A-dependent)
212  // if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance)
213  if(lastP == pMom) // Do not recalculate
214  {
215  CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom); // Update param's only
216  return lastCS*millibarn; // Use theLastCS
217  }
218  in = true; // This is the case when the isotop is found in DB
219  // Momentum pMom is in IU ! @@ Units
220  lastCS=CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom); // read & update
221  if(lastCS<=0. && pEn>lastTH) // Correct the threshold
222  {
223  lastTH=pEn;
224  }
225  break; // Go out of the LOOP with found lastI
226  }
227  } // End of attampt to find the nucleus in DB
228  if(!in) // This nucleus has not been calculated previously
229  {
231  lastCS=CalculateCrossSection(fCS,0,lastI,321,lastZ,lastN,pMom);//calculate&create
232  if(lastCS<=0.)
233  {
234  lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
235  if(pEn>lastTH)
236  {
237  lastTH=pEn;
238  }
239  }
240  colN.push_back(tgN);
241  colZ.push_back(tgZ);
242  colP.push_back(pMom);
243  colTH.push_back(lastTH);
244  colCS.push_back(lastCS);
245  return lastCS*millibarn;
246  } // End of creation of the new set of parameters
247  else
248  {
249  colP[lastI]=pMom;
250  colCS[lastI]=lastCS;
251  }
252  return lastCS*millibarn;
253 }
254 
255 // Calculation of total elastic cross section (p in IU, CS in mb) @@ Units (?)
256 // F=0 - create AMDB, F=-1 - read&update AMDB, F=1 - update AMDB (sinchro with higher AMDB)
257 G4double G4ChipsKaonPlusElasticXS::CalculateCrossSection(G4bool CS, G4int F,
258  G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
259 {
260  G4double pMom=pIU/GeV; // All calculations are in GeV
261  onlyCS=CS; // Flag to calculate only CS (not Si/Bi)
262  lastLP=std::log(pMom); // Make a logarithm of the momentum for calculation
263  if(F) // This isotope was found in AMDB =>RETRIEVE/UPDATE
264  {
265  if(F<0) // the AMDB must be loded
266  {
267  lastPIN = PIN[I]; // Max log(P) initialised for this table set
268  lastPAR = PAR[I]; // Pointer to the parameter set
269  lastCST = CST[I]; // Pointer to the total sross-section table
270  lastSST = SST[I]; // Pointer to the first squared slope
271  lastS1T = S1T[I]; // Pointer to the first mantissa
272  lastB1T = B1T[I]; // Pointer to the first slope
273  lastS2T = S2T[I]; // Pointer to the second mantissa
274  lastB2T = B2T[I]; // Pointer to the second slope
275  lastS3T = S3T[I]; // Pointer to the third mantissa
276  lastB3T = B3T[I]; // Pointer to the rhird slope
277  lastS4T = S4T[I]; // Pointer to the 4-th mantissa
278  lastB4T = B4T[I]; // Pointer to the 4-th slope
279  }
280  if(lastLP>lastPIN && lastLP<lPMax)
281  {
282  lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
283  PIN[I]=lastPIN; // Remember the new P-Limit of the tables
284  }
285  }
286  else // This isotope wasn't initialized => CREATE
287  {
288  lastPAR = new G4double[nPoints]; // Allocate memory for parameters of CS function
289  lastPAR[nLast]=0; // Initialization for VALGRIND
290  lastCST = new G4double[nPoints]; // Allocate memory for Tabulated CS function
291  lastSST = new G4double[nPoints]; // Allocate memory for Tabulated first sqaredSlope
292  lastS1T = new G4double[nPoints]; // Allocate memory for Tabulated first mantissa
293  lastB1T = new G4double[nPoints]; // Allocate memory for Tabulated first slope
294  lastS2T = new G4double[nPoints]; // Allocate memory for Tabulated second mantissa
295  lastB2T = new G4double[nPoints]; // Allocate memory for Tabulated second slope
296  lastS3T = new G4double[nPoints]; // Allocate memory for Tabulated third mantissa
297  lastB3T = new G4double[nPoints]; // Allocate memory for Tabulated third slope
298  lastS4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th mantissa
299  lastB4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th slope
300  lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
301  PIN.push_back(lastPIN); // Fill parameters of CS function to AMDB
302  PAR.push_back(lastPAR); // Fill parameters of CS function to AMDB
303  CST.push_back(lastCST); // Fill Tabulated CS function to AMDB
304  SST.push_back(lastSST); // Fill Tabulated first sq.slope to AMDB
305  S1T.push_back(lastS1T); // Fill Tabulated first mantissa to AMDB
306  B1T.push_back(lastB1T); // Fill Tabulated first slope to AMDB
307  S2T.push_back(lastS2T); // Fill Tabulated second mantissa to AMDB
308  B2T.push_back(lastB2T); // Fill Tabulated second slope to AMDB
309  S3T.push_back(lastS3T); // Fill Tabulated third mantissa to AMDB
310  B3T.push_back(lastB3T); // Fill Tabulated third slope to AMDB
311  S4T.push_back(lastS4T); // Fill Tabulated 4-th mantissa to AMDB
312  B4T.push_back(lastB4T); // Fill Tabulated 4-th slope to AMDB
313  } // End of creation/update of the new set of parameters and tables
314  // =----------= NOW Update (if necessary) and Calculate the Cross Section =----------=
315  if(lastLP>lastPIN && lastLP<lPMax)
316  {
317  lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
318  }
319  if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
320  if(lastLP>lPMin && lastLP<=lastPIN) // Linear fit is made using precalculated tables
321  {
322  if(lastLP==lastPIN)
323  {
324  G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
325  G4int blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
326  if(blast<0 || blast>=nLast) G4cout<<"G4QKPElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl;
327  lastSIG = lastCST[blast];
328  if(!onlyCS) // Skip the differential cross-section parameters
329  {
330  theSS = lastSST[blast];
331  theS1 = lastS1T[blast];
332  theB1 = lastB1T[blast];
333  theS2 = lastS2T[blast];
334  theB2 = lastB2T[blast];
335  theS3 = lastS3T[blast];
336  theB3 = lastB3T[blast];
337  theS4 = lastS4T[blast];
338  theB4 = lastB4T[blast];
339  }
340  }
341  else
342  {
343  G4double shift=(lastLP-lPMin)/dlnP; // a shift from the beginning of the table
344  G4int blast=static_cast<int>(shift); // the lower bin number
345  if(blast<0) blast=0;
346  if(blast>=nLast) blast=nLast-1; // low edge of the last bin
347  shift-=blast; // step inside the unit bin
348  G4int lastL=blast+1; // the upper bin number
349  G4double SIGL=lastCST[blast]; // the basic value of the cross-section
350  lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
351  if(!onlyCS) // Skip the differential cross-section parameters
352  {
353  G4double SSTL=lastSST[blast]; // the low bin of the first squared slope
354  theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
355  G4double S1TL=lastS1T[blast]; // the low bin of the first mantissa
356  theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
357  G4double B1TL=lastB1T[blast]; // the low bin of the first slope
358  theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
359  G4double S2TL=lastS2T[blast]; // the low bin of the second mantissa
360  theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
361  G4double B2TL=lastB2T[blast]; // the low bin of the second slope
362  theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
363  G4double S3TL=lastS3T[blast]; // the low bin of the third mantissa
364  theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
365  G4double B3TL=lastB3T[blast]; // the low bin of the third slope
366  theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
367  G4double S4TL=lastS4T[blast]; // the low bin of the 4-th mantissa
368  theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
369  G4double B4TL=lastB4T[blast]; // the low bin of the 4-th slope
370  theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
371  }
372  }
373  }
374  else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
375  if(lastSIG<0.) lastSIG = 0.; // @@ a Warning print can be added
376  return lastSIG;
377 }
378 
379 // It has parameter sets for all tZ/tN/PDG, using them the tables can be created/updated
380 G4double G4ChipsKaonPlusElasticXS::GetPTables(G4double LP,G4double ILP, G4int PDG,
381  G4int tgZ, G4int tgN)
382 {
383  if(PDG == 321)
384  {
385  // -- Total pp elastic cross section cs & s1/b1 (main), s2/b2 (tail1), s3/b3 (tail2) --
386  //p2=p*p;p3=p2*p;sp=sqrt(p);p2s=p2*sp;lp=log(p);dl1=lp-(3.=par(3));p4=p2*p2; p=|3-mom|
387  //CS=2.865/p2s/(1+.0022/p2s)+(18.9+.6461*dl1*dl1+9./p)/(1.+.425*lp)/(1.+.4276/p4);
388  // par(0) par(7) par(1) par(2) par(4) par(5) par(6)
389  //dl2=lp-5., s1=(74.+3.*dl2*dl2)/(1+3.4/p4/p)+(.2/p2+17.*p)/(p4+.001*sp),
390  // par(8) par(9) par(10) par(11) par(12)par(13) par(14)
391  // b1=8.*p**.055/(1.+3.64/p3); s2=5.e-5+4000./(p4+1500.*p); b2=.46+1.2e6/(p4+3.5e6/sp);
392  // par(15) par(16) par(17) par(18) par(19) par(20) par(21) par(22) par(23)
393  // s3=5.e-5+1.e10/(p4*p4+8.5e8*p2+1.e10); b3=1.1+3.4e6/(p4+6.8e6); ss=0.
394  // par(24) par(25) par(26) par(27) par(28) par(29) par(30) par(31)
395  //
396  if(lastPAR[nLast]!=pwd) // A unique flag to avoid the repeatable definition
397  {
398  if ( tgZ == 1 && tgN == 0 )
399  {
400  for (G4int ip=0; ip<n_kppel; ip++) lastPAR[ip]=kpp_el[ip]; // KPlus+P
401  }
402  else
403  {
404  G4double a=tgZ+tgN;
405  G4double sa=std::sqrt(a);
406  G4double ssa=std::sqrt(sa);
407  G4double asa=a*sa;
408  G4double a2=a*a;
409  G4double a3=a2*a;
410  G4double a4=a3*a;
411  G4double a5=a4*a;
412  G4double a6=a4*a2;
413  G4double a7=a6*a;
414  G4double a8=a7*a;
415  G4double a9=a8*a;
416  G4double a10=a5*a5;
417  G4double a12=a6*a6;
418  G4double a14=a7*a7;
419  G4double a16=a8*a8;
420  G4double a17=a16*a;
421  //G4double a20=a16*a4;
422  G4double a32=a16*a16;
423  // Reaction cross-section parameters (kpael_fit.f)
424  lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa)); // p1
425  lastPAR[1]=.75*asa/(1.+.009*a); // p2
426  lastPAR[2]=.9*asa*ssa/(1.+.03*a); // p3
427  lastPAR[3]=3.; // p4
428  lastPAR[4]=4.2; // p5
429  lastPAR[5]=0.; // p6 not used
430  lastPAR[6]=0.; // p7 not used
431  lastPAR[7]=0.; // p8 not used
432  lastPAR[8]=0.; // p9 not used
433  // @@ the differential cross-section is parameterized separately for A>6 & A<7
434  if(a<6.5)
435  {
436  G4double a28=a16*a12;
437  // The main pre-exponent (pel_sg)
438  lastPAR[ 9]=4000*a; // p1
439  lastPAR[10]=1.2e7*a8+380*a17; // p2
440  lastPAR[11]=.7/(1.+4.e-12*a16); // p3
441  lastPAR[12]=2.5/a8/(a4+1.e-16*a32); // p4
442  lastPAR[13]=.28*a; // p5
443  lastPAR[14]=1.2*a2+2.3; // p6
444  lastPAR[15]=3.8/a; // p7
445  // The main slope (pel_sl)
446  lastPAR[16]=.01/(1.+.0024*a5); // p1
447  lastPAR[17]=.2*a; // p2
448  lastPAR[18]=9.e-7/(1.+.035*a5); // p3
449  lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a); // p4
450  // The main quadratic (pel_sh)
451  lastPAR[20]=2.25*a3; // p1
452  lastPAR[21]=18.; // p2
453  lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7); // p3
454  lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a); // p4
455  // The 1st max pre-exponent (pel_qq)
456  lastPAR[24]=1.e5/(a8+2.5e12/a16); // p1
457  lastPAR[25]=8.e7/(a12+1.e-27*a28*a28); // p2
458  lastPAR[26]=.0006*a3; // p3
459  // The 1st max slope (pel_qs)
460  lastPAR[27]=10.+4.e-8*a12*a; // p1
461  lastPAR[28]=.114; // p2
462  lastPAR[29]=.003; // p3
463  lastPAR[30]=2.e-23; // p4
464  // The effective pre-exponent (pel_ss)
465  lastPAR[31]=1./(1.+.0001*a8); // p1
466  lastPAR[32]=1.5e-4/(1.+5.e-6*a12); // p2
467  lastPAR[33]=.03; // p3
468  // The effective slope (pel_sb)
469  lastPAR[34]=a/2; // p1
470  lastPAR[35]=2.e-7*a4; // p2
471  lastPAR[36]=4.; // p3
472  lastPAR[37]=64./a3; // p4
473  // The gloria pre-exponent (pel_us)
474  lastPAR[38]=1.e8*std::exp(.32*asa); // p1
475  lastPAR[39]=20.*std::exp(.45*asa); // p2
476  lastPAR[40]=7.e3+2.4e6/a5; // p3
477  lastPAR[41]=2.5e5*std::exp(.085*a3); // p4
478  lastPAR[42]=2.5*a; // p5
479  // The gloria slope (pel_ub)
480  lastPAR[43]=920.+.03*a8*a3; // p1
481  lastPAR[44]=93.+.0023*a12; // p2
482  }
483  else
484  {
485  G4double p1a10=2.2e-28*a10;
486  G4double r4a16=6.e14/a16;
487  G4double s4a16=r4a16*r4a16;
488  // a24
489  // a36
490  // The main pre-exponent (peh_sg)
491  lastPAR[ 9]=4.5*std::pow(a,1.15); // p1
492  lastPAR[10]=.06*std::pow(a,.6); // p2
493  lastPAR[11]=.6*a/(1.+2.e15/a16); // p3
494  lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32); // p4
495  lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5); // p5
496  lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12); // p6
497  // The main slope (peh_sl)
498  lastPAR[15]=400./a12+2.e-22*a9; // p1
499  lastPAR[16]=1.e-32*a12/(1.+5.e22/a14); // p2
500  lastPAR[17]=1000./a2+9.5*sa*ssa; // p3
501  lastPAR[18]=4.e-6*a*asa+1.e11/a16; // p4
502  lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16); // p5
503  lastPAR[20]=9.+100./a; // p6
504  // The main quadratic (peh_sh)
505  lastPAR[21]=.002*a3+3.e7/a6; // p1
506  lastPAR[22]=7.e-15*a4*asa; // p2
507  lastPAR[23]=9000./a4; // p3
508  // The 1st max pre-exponent (peh_qq)
509  lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4); // p1
510  lastPAR[25]=1.e-5*a2+2.e14/a16; // p2
511  lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12); // p3
512  lastPAR[27]=.016*asa/(1.+5.e16/a16); // p4
513  // The 1st max slope (peh_qs)
514  lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14)); // p1
515  lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11); // p2
516  lastPAR[30]=11.*a3/(1.+7.e23/a16/a8); // p3
517  lastPAR[31]=100./asa; // p4
518  // The 2nd max pre-exponent (peh_ss)
519  lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4); // p1
520  lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8); // p2
521  lastPAR[34]=1.3+3.e5/a4; // p3
522  lastPAR[35]=500./(a2+50.)+3; // p4
523  lastPAR[36]=1.e-9/a+s4a16*s4a16; // p5
524  // The 2nd max slope (peh_sb)
525  lastPAR[37]=.4*asa+3.e-9*a6; // p1
526  lastPAR[38]=.0005*a5; // p2
527  lastPAR[39]=.002*a5; // p3
528  lastPAR[40]=10.; // p4
529  // The effective pre-exponent (peh_us)
530  lastPAR[41]=.05+.005*a; // p1
531  lastPAR[42]=7.e-8/sa; // p2
532  lastPAR[43]=.8*sa; // p3
533  lastPAR[44]=.02*sa; // p4
534  lastPAR[45]=1.e8/a3; // p5
535  lastPAR[46]=3.e32/(a32+1.e32); // p6
536  // The effective slope (peh_ub)
537  lastPAR[47]=24.; // p1
538  lastPAR[48]=20./sa; // p2
539  lastPAR[49]=7.e3*a/(sa+1.); // p3
540  lastPAR[50]=900.*sa/(1.+500./a3); // p4
541  }
542  // Parameter for lowEnergyNeutrons
543  lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
544  }
545  lastPAR[nLast]=pwd;
546  // and initialize the zero element of the table
547  G4double lp=lPMin; // ln(momentum)
548  G4bool memCS=onlyCS; // ??
549  onlyCS=false;
550  lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables
551  onlyCS=memCS;
552  lastSST[0]=theSS;
553  lastS1T[0]=theS1;
554  lastB1T[0]=theB1;
555  lastS2T[0]=theS2;
556  lastB2T[0]=theB2;
557  lastS3T[0]=theS3;
558  lastB3T[0]=theB3;
559  lastS4T[0]=theS4;
560  lastB4T[0]=theB4;
561  }
562  if(LP>ILP)
563  {
564  G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1; // already inited till this
565  if(ini<0) ini=0;
566  if(ini<nPoints)
567  {
568  G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1; // final bin of initialization
569  if(fin>=nPoints) fin=nLast; // Limit of the tabular initialization
570  if(fin>=ini)
571  {
572  G4double lp=0.;
573  for(G4int ip=ini; ip<=fin; ip++) // Calculate tabular CS,S1,B1,S2,B2,S3,B3
574  {
575  lp=lPMin+ip*dlnP; // ln(momentum)
576  G4bool memCS=onlyCS;
577  onlyCS=false;
578  lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN); // Calculate AMDB tables (ret CS)
579  onlyCS=memCS;
580  lastSST[ip]=theSS;
581  lastS1T[ip]=theS1;
582  lastB1T[ip]=theB1;
583  lastS2T[ip]=theS2;
584  lastB2T[ip]=theB2;
585  lastS3T[ip]=theS3;
586  lastB3T[ip]=theB3;
587  lastS4T[ip]=theS4;
588  lastB4T[ip]=theB4;
589  }
590  return lp;
591  }
592  else G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
593  <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
594  <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
595  }
596  else G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
597  <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
598  <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
599  }
600  }
601  else
602  {
603  // G4cout<<"*Error*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG<<", Z="<<tgZ
604  // <<", N="<<tgN<<", while it is defined only for PDG=321"<<G4endl;
605  // throw G4QException("G4ChipsKaonPlusElasticXS::GetPTables:onlyK+ is implemented");
606  G4ExceptionDescription ed;
607  ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
608  << ", while it is defined only for PDG=321 (K+) " << G4endl;
609  G4Exception("G4ChipsKaonPlusElasticXS::GetPTables()", "HAD_CHPS_0000",
610  FatalException, ed);
611  }
612  return ILP;
613 }
614 
615 // Returns Q2=-t in independent units (MeV^2) (all internal calculations are in GeV)
616 G4double G4ChipsKaonPlusElasticXS::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
617 {
618  if(PDG!=321) G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetExT:PDG="<<PDG<<G4endl;
619  if(onlyCS) G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetExT: onlyCS=1"<<G4endl;
620  if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
621  G4double q2=0.;
622  if(tgZ==1 && tgN==0) // ===> p+p=p+p
623  {
624  G4double E1=lastTM*theB1;
625  G4double R1=(1.-std::exp(-E1));
626  G4double E2=lastTM*theB2;
627  G4double R2=(1.-std::exp(-E2*E2*E2));
628  G4double E3=lastTM*theB3;
629  G4double R3=(1.-std::exp(-E3));
630  G4double I1=R1*theS1/theB1;
631  G4double I2=R2*theS2;
632  G4double I3=R3*theS3;
633  G4double I12=I1+I2;
634  G4double rand=(I12+I3)*G4UniformRand();
635  if (rand<I1 )
636  {
637  G4double ran=R1*G4UniformRand();
638  if(ran>1.) ran=1.;
639  q2=-std::log(1.-ran)/theB1;
640  }
641  else if(rand<I12)
642  {
643  G4double ran=R2*G4UniformRand();
644  if(ran>1.) ran=1.;
645  q2=-std::log(1.-ran);
646  if(q2<0.) q2=0.;
647  q2=std::pow(q2,third)/theB2;
648  }
649  else
650  {
651  G4double ran=R3*G4UniformRand();
652  if(ran>1.) ran=1.;
653  q2=-std::log(1.-ran)/theB3;
654  }
655  }
656  else
657  {
658  G4double a=tgZ+tgN;
659  G4double E1=lastTM*(theB1+lastTM*theSS);
660  G4double R1=(1.-std::exp(-E1));
661  G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
662  G4double tm2=lastTM*lastTM;
663  G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st)
664  if(a>6.5)E2*=tm2; // for heavy nuclei
665  G4double R2=(1.-std::exp(-E2));
666  G4double E3=lastTM*theB3;
667  if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA
668  G4double R3=(1.-std::exp(-E3));
669  G4double E4=lastTM*theB4;
670  G4double R4=(1.-std::exp(-E4));
671  G4double I1=R1*theS1;
672  G4double I2=R2*theS2;
673  G4double I3=R3*theS3;
674  G4double I4=R4*theS4;
675  G4double I12=I1+I2;
676  G4double I13=I12+I3;
677  G4double rand=(I13+I4)*G4UniformRand();
678  if(rand<I1)
679  {
680  G4double ran=R1*G4UniformRand();
681  if(ran>1.) ran=1.;
682  q2=-std::log(1.-ran)/theB1;
683  if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
684  }
685  else if(rand<I12)
686  {
687  G4double ran=R2*G4UniformRand();
688  if(ran>1.) ran=1.;
689  q2=-std::log(1.-ran)/theB2;
690  if(q2<0.) q2=0.;
691  if(a<6.5) q2=std::pow(q2,third);
692  else q2=std::pow(q2,fifth);
693  }
694  else if(rand<I13)
695  {
696  G4double ran=R3*G4UniformRand();
697  if(ran>1.) ran=1.;
698  q2=-std::log(1.-ran)/theB3;
699  if(q2<0.) q2=0.;
700  if(a>6.5) q2=std::pow(q2,sevth);
701  }
702  else
703  {
704  G4double ran=R4*G4UniformRand();
705  if(ran>1.) ran=1.;
706  q2=-std::log(1.-ran)/theB4;
707  if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0)
708  }
709  }
710  if(q2<0.) q2=0.;
711  if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonPlusElasticCS::GetExchT: -t="<<q2<<G4endl;
712  if(q2>lastTM)
713  {
714  q2=lastTM;
715  }
716  return q2*GeVSQ;
717 }
718 
719 // Returns B in independent units (MeV^-2) (all internal calculations are in GeV) see ExT
720 G4double G4ChipsKaonPlusElasticXS::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
721 {
722  if(onlyCS)G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetSl:onlCS=true"<<G4endl;
723  if(lastLP<-4.3) return 0.; // S-wave for p<14 MeV/c (kinE<.1MeV)
724  if(PDG != 321)
725  {
726  G4ExceptionDescription ed;
727  ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
728  << ", while it is defined only for PDG=321 (K+)" << G4endl;
729  G4Exception("G4ChipsKaonPlusElasticXS::GetSlope()", "HAD_CHPS_0000",
730  FatalException, ed);
731  }
732  if(theB1<0.) theB1=0.;
733  if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonPlusElCS::GetSlope:B1="<<theB1<<G4endl;
734  return theB1/GeVSQ;
735 }
736 
737 // Returns half max(Q2=-t) in independent units (MeV^2)
738 G4double G4ChipsKaonPlusElasticXS::GetHMaxT()
739 {
740  return lastTM*HGeVSQ;
741 }
742 
743 // lastLP is used, so calculating tables, one need to remember and then recover lastLP
744 G4double G4ChipsKaonPlusElasticXS::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
745  G4int tgN)
746 {
747  if(PDG!=321)G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetTaV:PDG="<<PDG<<G4endl;
748 
749  //AR-24Apr2018 Switch to allow transuranic elements
750  const G4bool isHeavyElementAllowed = true;
751  if(tgZ<0 || ( !isHeavyElementAllowed && tgZ>92))
752  {
753  G4cout<<"*Warning*G4QKaonPlusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
754  return 0.;
755  }
756  G4int iZ=tgZ-1; // Z index
757  if(iZ<0)
758  {
759  iZ=0; // conversion of the neutron target to the proton target
760  tgZ=1;
761  tgN=0;
762  }
763  G4double p=std::exp(lp); // momentum
764  G4double sp=std::sqrt(p); // sqrt(p)
765  G4double p2=p*p;
766  G4double p3=p2*p;
767  G4double p4=p3*p;
768  if ( tgZ == 1 && tgN == 0 ) // KaonPlus+P
769  {
770  G4double dl2=lp-lastPAR[11];
771  theSS=lastPAR[34];
772  theS1=(lastPAR[12]+lastPAR[13]*dl2*dl2)/(1.+lastPAR[14]/p4/p)+
773  (lastPAR[15]/p2+lastPAR[16]*p)/(p4+lastPAR[17]*sp);
774  theB1=lastPAR[18]*std::pow(p,lastPAR[19])/(1.+lastPAR[20]/p3);
775  theS2=lastPAR[21]+lastPAR[22]/(p4+lastPAR[23]*p);
776  theB2=lastPAR[24]+lastPAR[25]/(p4+lastPAR[26]/sp);
777  theS3=lastPAR[27]+lastPAR[28]/(p4*p4+lastPAR[29]*p2+lastPAR[30]);
778  theB3=lastPAR[31]+lastPAR[32]/(p4+lastPAR[33]);
779  theS4=0.;
780  theB4=0.;
781  // Returns the total elastic pim-p cross-section (to avoid spoiling lastSIG)
782  G4double dp=lp-lastPAR[4];
783 //G4cout<<"lastPAR[8] "<<lastPAR[8]<<" lastPAR[9] "<<lastPAR[9]<<" lastPAR[10] "<<lastPAR[10]<<G4endl;
784  return lastPAR[0]/(lastPAR[2]+sqr(p-lastPAR[1]))+(lastPAR[3]*dp*dp+lastPAR[5])/
785  (1.-lastPAR[6]/sp+lastPAR[7]/p4)
786  +lastPAR[8]/(sqr(p-lastPAR[9])+lastPAR[10]); // Uzhi
787 
788  }
789  else
790  {
791  G4double p5=p4*p;
792  G4double p6=p5*p;
793  G4double p8=p6*p2;
794  G4double p10=p8*p2;
795  G4double p12=p10*p2;
796  G4double p16=p8*p8;
797  //G4double p24=p16*p8;
798  G4double dl=lp-5.;
799  G4double a=tgZ+tgN;
800  G4double pah=std::pow(p,a/2);
801  G4double pa=pah*pah;
802  G4double pa2=pa*pa;
803  if(a<6.5)
804  {
805  theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
806  (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
807  theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
808  theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
809  theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
810  theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
811  theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
812  theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
813  theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
814  lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
815  theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
816  }
817  else
818  {
819  theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
820  lastPAR[13]/(p5+lastPAR[14]/p16);
821  theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
822  lastPAR[17]/(1.+lastPAR[18]/p4);
823  theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
824  theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
825  theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
826  theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
827  lastPAR[33]/(1.+lastPAR[34]/p6);
828  theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
829  theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
830  (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
831  theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
832  }
833  // Returns the total elastic (n/p)A cross-section (to avoid spoiling lastSIG)
834  G4double dlp=lp-lastPAR[4]; // ax
835  // p1 p2 p3 p4
836  return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p2)/(1.+lastPAR[3]/p2/sp);
837  }
838  return 0.;
839 } // End of GetTableValues
840 
841 // Returns max -t=Q2 (GeV^2) for the momentum pP(GeV) and the target nucleus (tgN,tgZ)
842 G4double G4ChipsKaonPlusElasticXS::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
843  G4double pP)
844 {
845  G4double pP2=pP*pP; // squared momentum of the projectile
846  if(tgZ || tgN>-1) // ---> pipA
847  {
848  G4double mt=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(tgZ,tgZ+tgN,0)->GetPDGMass()*.001; // Target mass in GeV
849 
850  G4double dmt=mt+mt;
851  G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt; // Mondelstam mds
852  return dmt*dmt*pP2/mds;
853  }
854  else
855  {
856  G4ExceptionDescription ed;
857  ed << "PDG = " << PDG << ",Z = " << tgZ << ", N = " << tgN
858  << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
859  G4Exception("G4ChipsKaonPlusElasticXS::GetQ2max()", "HAD_CHPS_0000",
860  FatalException, ed);
861  return 0;
862  }
863 }