ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChipsKaonMinusElasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ChipsKaonMinusElasticXS.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: G4ChipsKaonMinusElasticXS 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 "G4ChipsKaonMinusElasticXS.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DynamicParticle.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4KaonMinus.hh"
44 #include "G4Nucleus.hh"
45 #include "G4ParticleTable.hh"
46 #include "G4NucleiProperties.hh"
47 #include "G4IonTable.hh"
48 #include "G4AutoLock.hh"
49 
50 // factory
51 #include "G4CrossSectionFactory.hh"
52 //
53 G4_DECLARE_XS_FACTORY(G4ChipsKaonMinusElasticXS);
54 
55 
56 namespace {
57  G4double mK;//= G4KaonMinus::KaonMinus()->GetPDGMass()*.001; // MeV to GeV//Cannot initialize here, needs particles
58  G4double mK2;//= mK*mK;
59  G4Mutex initM = G4MUTEX_INITIALIZER;
60  const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
61  const G4double third=1./3.;
62  const G4double fifth=1./5.;
63  const G4double sevth=1./7.;
64  const G4double pwd=2727;
65  const G4int n_kmpel=36; // #of parameters for pp-elastic ( < nPoints=128)
66  // -0- -1- -2- -3- -4- -5- -6- -7- -8- -9--10- -11--12-
67  const G4double kmp_el[n_kmpel]={5.2,.0557,3.5,2.23,.7,.075,.004,.39,.000156,.15,1.,.0156,5.,
68  74.,3.,3.4,.2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,
69  1.2e6,3.5e6,5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
70  // -13-14--15-16--17--18--19- -20- -21- -22- -23- -24- -25-
71  // -26- -27- -28- -29- -30- -31- -32- -33- -34- -35-
72  const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
73 }
74 
75 G4ChipsKaonMinusElasticXS::G4ChipsKaonMinusElasticXS():G4VCrossSectionDataSet(Default_Name()), nPoints(128), nLast(nPoints-1)
76 {
77  G4AutoLock l(&initM);
78  mK = G4KaonMinus::KaonMinus()->GetPDGMass()*.001;
79  mK2 = mK*mK;
80  l.unlock();
81  lPMin=-8.; //Min tabulatedLogarithmMomentum/D
82  lPMax= 8.; //Max tabulatedLogarithmMomentum/D
83  dlnP=(lPMax-lPMin)/nLast;// LogStep inTable /D
84  onlyCS=true;//Flag toCalculOnlyCS(not Si/Bi)/L
85  lastSIG=0.; //Last calculated cross section /L
86  lastLP=-10.;//LastLog(mom_of IncidentHadron)/L
87  lastTM=0.; //Last t_maximum /L
88  theSS=0.; //TheLastSqSlope of 1st difr.Max/L
89  theS1=0.; //TheLastMantissa of 1st difrMax/L
90  theB1=0.; //TheLastSlope of 1st difructMax/L
91  theS2=0.; //TheLastMantissa of 2nd difrMax/L
92  theB2=0.; //TheLastSlope of 2nd difructMax/L
93  theS3=0.; //TheLastMantissa of 3d difr.Max/L
94  theB3=0.; //TheLastSlope of 3d difruct.Max/L
95  theS4=0.; //TheLastMantissa of 4th difrMax/L
96  theB4=0.; //TheLastSlope of 4th difructMax/L
97  lastTZ=0; // Last atomic number of theTarget
98  lastTN=0; // Last # of neutrons in theTarget
99  lastPIN=0.;// Last initialized max momentum
100  lastCST=0; // Elastic cross-section table
101  lastPAR=0; // ParametersForFunctionCalculation
102  lastSST=0; // E-dep ofSqardSlope of 1st difMax
103  lastS1T=0; // E-dep of mantissa of 1st dif.Max
104  lastB1T=0; // E-dep of the slope of 1st difMax
105  lastS2T=0; // E-dep of mantissa of 2nd difrMax
106  lastB2T=0; // E-dep of the slope of 2nd difMax
107  lastS3T=0; // E-dep of mantissa of 3d difr.Max
108  lastB3T=0; // E-dep of the slope of 3d difrMax
109  lastS4T=0; // E-dep of mantissa of 4th difrMax
110  lastB4T=0; // E-dep of the slope of 4th difMax
111  lastN=0; // The last N of calculated nucleus
112  lastZ=0; // The last Z of calculated nucleus
113  lastP=0.; // LastUsed inCrossSection Momentum
114  lastTH=0.; // Last threshold momentum
115  lastCS=0.; // Last value of the Cross Section
116  lastI=0; // The last position in the DAMDB
117 }
118 
119 G4ChipsKaonMinusElasticXS::~G4ChipsKaonMinusElasticXS()
120 {
121  std::vector<G4double*>::iterator pos;
122  for (pos=CST.begin(); pos<CST.end(); pos++)
123  { delete [] *pos; }
124  CST.clear();
125  for (pos=PAR.begin(); pos<PAR.end(); pos++)
126  { delete [] *pos; }
127  PAR.clear();
128  for (pos=SST.begin(); pos<SST.end(); pos++)
129  { delete [] *pos; }
130  SST.clear();
131  for (pos=S1T.begin(); pos<S1T.end(); pos++)
132  { delete [] *pos; }
133  S1T.clear();
134  for (pos=B1T.begin(); pos<B1T.end(); pos++)
135  { delete [] *pos; }
136  B1T.clear();
137  for (pos=S2T.begin(); pos<S2T.end(); pos++)
138  { delete [] *pos; }
139  S2T.clear();
140  for (pos=B2T.begin(); pos<B2T.end(); pos++)
141  { delete [] *pos; }
142  B2T.clear();
143  for (pos=S3T.begin(); pos<S3T.end(); pos++)
144  { delete [] *pos; }
145  S3T.clear();
146  for (pos=B3T.begin(); pos<B3T.end(); pos++)
147  { delete [] *pos; }
148  B3T.clear();
149  for (pos=S4T.begin(); pos<S4T.end(); pos++)
150  { delete [] *pos; }
151  S4T.clear();
152  for (pos=B4T.begin(); pos<B4T.end(); pos++)
153  { delete [] *pos; }
154  B4T.clear();
155 }
156 
157 void
158 G4ChipsKaonMinusElasticXS::CrossSectionDescription(std::ostream& outFile) const
159 {
160  outFile << "G4ChipsKaonMinusElasticXS provides the elastic cross\n"
161  << "section for K- nucleus scattering as a function of incident\n"
162  << "momentum. The cross section is calculated using M. Kossov's\n"
163  << "CHIPS parameterization of cross section data.\n";
164 }
165 
166 G4bool G4ChipsKaonMinusElasticXS::IsIsoApplicable(const G4DynamicParticle*, G4int, G4int,
167  const G4Element*,
168  const G4Material*)
169 {
170  return true;
171 }
172 
173 // The main member function giving the collision cross section (P is in IU, CS is in mb)
174 // Make pMom in independent units ! (Now it is MeV)
175 G4double G4ChipsKaonMinusElasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
176  const G4Isotope*,
177  const G4Element*,
178  const G4Material*)
179 {
180  G4double pMom=Pt->GetTotalMomentum();
181  G4int tgN = A - tgZ;
182 
183  return GetChipsCrossSection(pMom, tgZ, tgN, -321);
184 }
185 
186 G4double G4ChipsKaonMinusElasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
187 {
188 
189  G4bool fCS = false;
190 
191  G4double pEn=pMom;
192  onlyCS=fCS;
193 
194  G4bool in=false; // By default the isotope must be found in the AMDB
195  lastP = 0.; // New momentum history (nothing to compare with)
196  lastN = tgN; // The last N of the calculated nucleus
197  lastZ = tgZ; // The last Z of the calculated nucleus
198  lastI = colN.size(); // Size of the Associative Memory DB in the heap
199  if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
200  { // The nucleus with projPDG is found in AMDB
201  if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
202  {
203  lastI=i;
204  lastTH =colTH[i]; // Last THreshold (A-dependent)
205  if(pEn<=lastTH)
206  {
207  return 0.; // Energy is below the Threshold value
208  }
209  lastP =colP [i]; // Last Momentum (A-dependent)
210  lastCS =colCS[i]; // Last CrossSect (A-dependent)
211  // if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance)
212  if(lastP == pMom) // Do not recalculate
213  {
214  CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom); // Update param's only
215  return lastCS*millibarn; // Use theLastCS
216  }
217  in = true; // This is the case when the isotop is found in DB
218  // Momentum pMom is in IU ! @@ Units
219  lastCS=CalculateCrossSection(fCS,-1,i,-321,lastZ,lastN,pMom); // read & update
220  if(lastCS<=0. && pEn>lastTH) // Correct the threshold
221  {
222  lastTH=pEn;
223  }
224  break; // Go out of the LOOP with found lastI
225  }
226  } // End of attampt to find the nucleus in DB
227  if(!in) // This nucleus has not been calculated previously
228  {
230  lastCS=CalculateCrossSection(fCS,0,lastI,-321,lastZ,lastN,pMom);//calculate&create
231  if(lastCS<=0.)
232  {
233  lastTH = 0; //ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
234  if(pEn>lastTH)
235  {
236  lastTH=pEn;
237  }
238  }
239  colN.push_back(tgN);
240  colZ.push_back(tgZ);
241  colP.push_back(pMom);
242  colTH.push_back(lastTH);
243  colCS.push_back(lastCS);
244  return lastCS*millibarn;
245  } // End of creation of the new set of parameters
246  else
247  {
248  colP[lastI]=pMom;
249  colCS[lastI]=lastCS;
250  }
251  return lastCS*millibarn;
252 }
253 
254 // Calculation of total elastic cross section (p in IU, CS in mb) @@ Units (?)
255 // F=0 - create AMDB, F=-1 - read&update AMDB, F=1 - update AMDB (sinchro with higher AMDB)
256 G4double G4ChipsKaonMinusElasticXS::CalculateCrossSection(G4bool CS, G4int F,
257  G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
258 {
259  G4double pMom=pIU/GeV; // All calculations are in GeV
260  onlyCS=CS; // Flag to calculate only CS (not Si/Bi)
261  lastLP=std::log(pMom); // Make a logarithm of the momentum for calculation
262  if(F) // This isotope was found in AMDB =>RETRIEVE/UPDATE
263  {
264  if(F<0) // the AMDB must be loded
265  {
266  lastPIN = PIN[I]; // Max log(P) initialised for this table set
267  lastPAR = PAR[I]; // Pointer to the parameter set
268  lastCST = CST[I]; // Pointer to the total sross-section table
269  lastSST = SST[I]; // Pointer to the first squared slope
270  lastS1T = S1T[I]; // Pointer to the first mantissa
271  lastB1T = B1T[I]; // Pointer to the first slope
272  lastS2T = S2T[I]; // Pointer to the second mantissa
273  lastB2T = B2T[I]; // Pointer to the second slope
274  lastS3T = S3T[I]; // Pointer to the third mantissa
275  lastB3T = B3T[I]; // Pointer to the rhird slope
276  lastS4T = S4T[I]; // Pointer to the 4-th mantissa
277  lastB4T = B4T[I]; // Pointer to the 4-th slope
278  }
279  if(lastLP>lastPIN && lastLP<lPMax)
280  {
281  lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
282  PIN[I]=lastPIN; // Remember the new P-Limit of the tables
283  }
284  }
285  else // This isotope wasn't initialized => CREATE
286  {
287  lastPAR = new G4double[nPoints]; // Allocate memory for parameters of CS function
288  lastPAR[nLast]=0; // Initialization for VALGRIND
289  lastCST = new G4double[nPoints]; // Allocate memory for Tabulated CS function
290  lastSST = new G4double[nPoints]; // Allocate memory for Tabulated first sqaredSlope
291  lastS1T = new G4double[nPoints]; // Allocate memory for Tabulated first mantissa
292  lastB1T = new G4double[nPoints]; // Allocate memory for Tabulated first slope
293  lastS2T = new G4double[nPoints]; // Allocate memory for Tabulated second mantissa
294  lastB2T = new G4double[nPoints]; // Allocate memory for Tabulated second slope
295  lastS3T = new G4double[nPoints]; // Allocate memory for Tabulated third mantissa
296  lastB3T = new G4double[nPoints]; // Allocate memory for Tabulated third slope
297  lastS4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th mantissa
298  lastB4T = new G4double[nPoints]; // Allocate memory for Tabulated 4-th slope
299  lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
300  PIN.push_back(lastPIN); // Fill parameters of CS function to AMDB
301  PAR.push_back(lastPAR); // Fill parameters of CS function to AMDB
302  CST.push_back(lastCST); // Fill Tabulated CS function to AMDB
303  SST.push_back(lastSST); // Fill Tabulated first sq.slope to AMDB
304  S1T.push_back(lastS1T); // Fill Tabulated first mantissa to AMDB
305  B1T.push_back(lastB1T); // Fill Tabulated first slope to AMDB
306  S2T.push_back(lastS2T); // Fill Tabulated second mantissa to AMDB
307  B2T.push_back(lastB2T); // Fill Tabulated second slope to AMDB
308  S3T.push_back(lastS3T); // Fill Tabulated third mantissa to AMDB
309  B3T.push_back(lastB3T); // Fill Tabulated third slope to AMDB
310  S4T.push_back(lastS4T); // Fill Tabulated 4-th mantissa to AMDB
311  B4T.push_back(lastB4T); // Fill Tabulated 4-th slope to AMDB
312  } // End of creation/update of the new set of parameters and tables
313  // =----------= NOW Update (if necessary) and Calculate the Cross Section =-----------=
314  if(lastLP>lastPIN && lastLP<lPMax)
315  {
316  lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
317  }
318  if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
319  if(lastLP>lPMin && lastLP<=lastPIN) // Linear fit is made using precalculated tables
320  {
321  if(lastLP==lastPIN)
322  {
323  G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
324  G4int blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
325  if(blast<0 || blast>=nLast) G4cout<<"G4QKMElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl;
326  lastSIG = lastCST[blast];
327  if(!onlyCS) // Skip the differential cross-section parameters
328  {
329  theSS = lastSST[blast];
330  theS1 = lastS1T[blast];
331  theB1 = lastB1T[blast];
332  theS2 = lastS2T[blast];
333  theB2 = lastB2T[blast];
334  theS3 = lastS3T[blast];
335  theB3 = lastB3T[blast];
336  theS4 = lastS4T[blast];
337  theB4 = lastB4T[blast];
338  }
339  }
340  else
341  {
342  G4double shift=(lastLP-lPMin)/dlnP; // a shift from the beginning of the table
343  G4int blast=static_cast<int>(shift); // the lower bin number
344  if(blast<0) blast=0;
345  if(blast>=nLast) blast=nLast-1; // low edge of the last bin
346  shift-=blast; // step inside the unit bin
347  G4int lastL=blast+1; // the upper bin number
348  G4double SIGL=lastCST[blast]; // the basic value of the cross-section
349  lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
350  if(!onlyCS) // Skip the differential cross-section parameters
351  {
352  G4double SSTL=lastSST[blast]; // the low bin of the first squared slope
353  theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
354  G4double S1TL=lastS1T[blast]; // the low bin of the first mantissa
355  theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
356  G4double B1TL=lastB1T[blast]; // the low bin of the first slope
357  theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
358  G4double S2TL=lastS2T[blast]; // the low bin of the second mantissa
359  theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
360  G4double B2TL=lastB2T[blast]; // the low bin of the second slope
361  theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
362  G4double S3TL=lastS3T[blast]; // the low bin of the third mantissa
363  theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
364  G4double B3TL=lastB3T[blast]; // the low bin of the third slope
365  theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
366  G4double S4TL=lastS4T[blast]; // the low bin of the 4-th mantissa
367  theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
368  G4double B4TL=lastB4T[blast]; // the low bin of the 4-th slope
369  theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
370  }
371  }
372  }
373  else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
374  if(lastSIG<0.) lastSIG = 0.; // @@ a Warning print can be added
375  return lastSIG;
376 }
377 
378 // It has parameter sets for all tZ/tN/PDG, using them the tables can be created/updated
379 G4double G4ChipsKaonMinusElasticXS::GetPTables(G4double LP,G4double ILP, G4int PDG,
380  G4int tgZ, G4int tgN)
381 {
382  // @@ At present all nA==pA ---------> Each neucleus can have not more than 51 parameters
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_kmpel; ip++) lastPAR[ip]=kmp_el[ip]; // PiMinus+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 (kmael_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]=.1*a2*ssa/(1.+.0015*a2/ssa); // p3
427  lastPAR[3]=1./(1.+500./a2); // 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*G4ChipsKaonMinusElasticXS::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*G4ChipsKaonMinusElasticXS::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*G4ChipsKaonMinusElasticXS::GetPTables: PDG="<<PDG<<", Z="<<tgZ
604  // <<", N="<<tgN<<", while it is defined only for PDG=-321"<<G4endl;
605  // throw G4QException("G4ChipsKaonMinusElasticXS::GetPTables:onlyK-'s 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("G4ChipsKaonMinusElasticXS::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 G4ChipsKaonMinusElasticXS::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
617 {
618  if(PDG==310 || PDG==130) PDG=-321;
619  if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetET:PDG="<<PDG<<G4endl;
620  if(onlyCS) G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetExT: onlyCS=1"<<G4endl;
621  if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
622  G4double q2=0.;
623  if(tgZ==1 && tgN==0) // ===> p+p=p+p
624  {
625  G4double E1=lastTM*theB1;
626  G4double R1=(1.-std::exp(-E1));
627  G4double E2=lastTM*theB2;
628  G4double R2=(1.-std::exp(-E2*E2*E2));
629  G4double E3=lastTM*theB3;
630  G4double R3=(1.-std::exp(-E3));
631  G4double I1=R1*theS1/theB1;
632  G4double I2=R2*theS2;
633  G4double I3=R3*theS3;
634  G4double I12=I1+I2;
635  G4double rand=(I12+I3)*G4UniformRand();
636  if (rand<I1 )
637  {
638  G4double ran=R1*G4UniformRand();
639  if(ran>1.) ran=1.;
640  q2=-std::log(1.-ran)/theB1;
641  }
642  else if(rand<I12)
643  {
644  G4double ran=R2*G4UniformRand();
645  if(ran>1.) ran=1.;
646  q2=-std::log(1.-ran);
647  if(q2<0.) q2=0.;
648  q2=std::pow(q2,third)/theB2;
649  }
650  else
651  {
652  G4double ran=R3*G4UniformRand();
653  if(ran>1.) ran=1.;
654  q2=-std::log(1.-ran)/theB3;
655  }
656  }
657  else
658  {
659  G4double a=tgZ+tgN;
660  G4double E1=lastTM*(theB1+lastTM*theSS);
661  G4double R1=(1.-std::exp(-E1));
662  G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
663  G4double tm2=lastTM*lastTM;
664  G4double E2=lastTM*tm2*theB2; // power 3 for lowA, 5 for HighA (1st)
665  if(a>6.5)E2*=tm2; // for heavy nuclei
666  G4double R2=(1.-std::exp(-E2));
667  G4double E3=lastTM*theB3;
668  if(a>6.5)E3*=tm2*tm2*tm2; // power 1 for lowA, 7 (2nd) for HighA
669  G4double R3=(1.-std::exp(-E3));
670  G4double E4=lastTM*theB4;
671  G4double R4=(1.-std::exp(-E4));
672  G4double I1=R1*theS1;
673  G4double I2=R2*theS2;
674  G4double I3=R3*theS3;
675  G4double I4=R4*theS4;
676  G4double I12=I1+I2;
677  G4double I13=I12+I3;
678  G4double rand=(I13+I4)*G4UniformRand();
679  if(rand<I1)
680  {
681  G4double ran=R1*G4UniformRand();
682  if(ran>1.) ran=1.;
683  q2=-std::log(1.-ran)/theB1;
684  if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
685  }
686  else if(rand<I12)
687  {
688  G4double ran=R2*G4UniformRand();
689  if(ran>1.) ran=1.;
690  q2=-std::log(1.-ran)/theB2;
691  if(q2<0.) q2=0.;
692  if(a<6.5) q2=std::pow(q2,third);
693  else q2=std::pow(q2,fifth);
694  }
695  else if(rand<I13)
696  {
697  G4double ran=R3*G4UniformRand();
698  if(ran>1.) ran=1.;
699  q2=-std::log(1.-ran)/theB3;
700  if(q2<0.) q2=0.;
701  if(a>6.5) q2=std::pow(q2,sevth);
702  }
703  else
704  {
705  G4double ran=R4*G4UniformRand();
706  if(ran>1.) ran=1.;
707  q2=-std::log(1.-ran)/theB4;
708  if(a<6.5) q2=lastTM-q2; // u reduced for lightA (starts from 0)
709  }
710  }
711  if(q2<0.) q2=0.;
712  if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<G4endl;
713  if(q2>lastTM)
714  {
715  q2=lastTM;
716  }
717  return q2*GeVSQ;
718 }
719 
720 // Returns B in independent units (MeV^-2) (all internal calculations are in GeV) see ExT
721 G4double G4ChipsKaonMinusElasticXS::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
722 {
723  if(onlyCS)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetSl:onlCS=true"<<G4endl;
724  if(lastLP<-4.3) return 0.; // S-wave for p<14 MeV/c (kinE<.1MeV)
725  if(PDG != -321)
726  {
727  // G4cout<<"*Error*G4ChipsKaonMinusElasticXS::GetSlope: PDG="<<PDG<<", Z="<<tgZ
728  // <<", N="<<tgN<<", while it is defined only for PDG=-321"<<G4endl;
729  // throw G4QException("G4ChipsKaonMinusElasticXS::GetSlope:Only K- is implemented");
730  G4ExceptionDescription ed;
731  ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
732  << ", while it is defined only for PDG=-321 (K-)" << G4endl;
733  }
734  if(theB1<0.) theB1=0.;
735  if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<G4endl;
736  return theB1/GeVSQ;
737 }
738 
739 // Returns half max(Q2=-t) in independent units (MeV^2)
740 G4double G4ChipsKaonMinusElasticXS::GetHMaxT()
741 {
742  return lastTM*HGeVSQ;
743 }
744 
745 // lastLP is used, so calculating tables, one need to remember and then recover lastLP
746 G4double G4ChipsKaonMinusElasticXS::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
747  G4int tgN)
748 {
749  if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetTV:PDG="<<PDG<<G4endl;
750 
751  //AR-24Apr2018 Switch to allow transuranic elements
752  const G4bool isHeavyElementAllowed = true;
753  if(tgZ<0 || ( !isHeavyElementAllowed && tgZ>92))
754  {
755  G4cout<<"*Warning*G4QKaonMinusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
756  return 0.;
757  }
758  G4int iZ=tgZ-1; // Z index
759  if(iZ<0)
760  {
761  iZ=0; // conversion of the neutron target to the proton target
762  tgZ=1;
763  tgN=0;
764  }
765  G4double p=std::exp(lp); // momentum
766  G4double sp=std::sqrt(p); // sqrt(p)
767  G4double psp=p*sp; // p*sqrt(p)
768  G4double p2=p*p;
769  G4double p3=p2*p;
770  G4double p4=p3*p;
771  if ( tgZ == 1 && tgN == 0 ) // KaonMinus+P
772  {
773  G4double dl2=lp-lastPAR[12];
774  theSS=lastPAR[35];
775  theS1=(lastPAR[13]+lastPAR[14]*dl2*dl2)/(1.+lastPAR[15]/p4/p)+
776  (lastPAR[16]/p2+lastPAR[17]*p)/(p4+lastPAR[18]*sp);
777  theB1=lastPAR[19]*std::pow(p,lastPAR[20])/(1.+lastPAR[21]/p3);
778  theS2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]*p);
779  theB2=lastPAR[25]+lastPAR[26]/(p4+lastPAR[27]/sp);
780  theS3=lastPAR[28]+lastPAR[29]/(p4*p4+lastPAR[30]*p2+lastPAR[31]);
781  theB3=lastPAR[32]+lastPAR[33]/(p4+lastPAR[34]);
782  theS4=0.;
783  theB4=0.;
784  // Returns the total elastic pim-p cross-section (to avoid spoiling lastSIG)
785  G4double dp=lp-lastPAR[2];
786  return lastPAR[0]/psp+(lastPAR[1]*dp*dp+lastPAR[3])/(1.-lastPAR[4]/sp+lastPAR[5]/p4)+
787  lastPAR[6]/(sqr(p-lastPAR[7])+lastPAR[8])+lastPAR[9]/(sqr(p-lastPAR[10])+lastPAR[11]);
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]/p3)/(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 G4ChipsKaonMinusElasticXS::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("G4ChipsKaonMinusElasticXS::GetQ2max()", "HAD_CHPS_0000",
860  FatalException, ed);
861  return 0;
862  }
863 }