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