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