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