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