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