Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QIonIonCrossSection.cc
Go to the documentation of this file.
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 // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 19-Aug-07
28 // GEANT4 tag $Name: not supported by cvs2svn $
29 //
30 //
31 // G4 Physics class: G4QIonIonCrossSection for gamma+A cross sections
32 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03
33 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04
34 // --------------------------------------------------------------------------------
35 // ****************************************************************************************
36 // This Header is a part of the CHIPS physics package (author: M. Kosov)
37 // ****************************************************************************************
38 // Short description: CHIPS cross-sectons for Ion-Ion interactions
39 // ---------------------------------------------------------------
40 //
41 //#define debug
42 //#define pdebug
43 //#define debug3
44 //#define debugn
45 //#define debugs
46 
47 #include "G4QIonIonCrossSection.hh"
48 #include "G4SystemOfUnits.hh"
49 
50 // Initialization of the
51 G4double* G4QIonIonCrossSection::lastLENI=0;// Pointer to the lastArray of LowEn Inelast CS
52 G4double* G4QIonIonCrossSection::lastHENI=0;// Pointer to the lastArray of HighEn InelastCS
53 G4double* G4QIonIonCrossSection::lastLENE=0;// Pointer to the lastArray of LowEn Elastic CS
54 G4double* G4QIonIonCrossSection::lastHENE=0;// Pointer to the lastArray of HighEn ElasticCS
55 G4int G4QIonIonCrossSection::lastPDG=0; // The last PDG code of the projectile
56 G4int G4QIonIonCrossSection::lastN=0; // The last N of calculated nucleus
57 G4int G4QIonIonCrossSection::lastZ=0; // The last Z of calculated nucleus
58 G4double G4QIonIonCrossSection::lastP=0.; // Last used in cross section Momentum
59 G4double G4QIonIonCrossSection::lastTH=0.; // Last threshold momentum
60 G4double G4QIonIonCrossSection::lastICS=0.;// Last value of the Inelastic Cross Section
61 G4double G4QIonIonCrossSection::lastECS=0.;// Last value of the Elastic Cross Section
62 G4int G4QIonIonCrossSection::lastI=0; // The last position in the DAMDB
63 
64 // Returns Pointer to the G4VQCrossSection class
65 G4VQCrossSection* G4QIonIonCrossSection::GetPointer()
66 {
67  static G4QIonIonCrossSection theCrossSection; //**Static body of Cross Section**
68  return &theCrossSection;
69 }
70 
71 // The main member function giving the collision cross section (P is in IU, CS is in mb)
72 // Make pMom in independent units !(Now it is MeV): fCS=true->Inelastic, fCS=false->Elastic
73 G4double G4QIonIonCrossSection::GetCrossSection(G4bool fCS, G4double pMom, G4int tZ,
74  G4int tN, G4int pPDG)
75 {
76  static G4int j; // A#0f records found in DB for this projectile
77  static std::vector <G4int> colPDG;// Vector of the projectile PDG code
78  static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops)
79  static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops)
80  static std::vector <G4double> colP; // Vector of last momenta for the reaction
81  static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
82  static std::vector <G4double> colICS;// Vector of last inelastic cross-sections
83  static std::vector <G4double> colECS;// Vector of last elastic cross-sections
84  // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
85 #ifdef pdebug
86  G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN
87  <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz="
88  <<colN.size()<<G4endl;
89 #endif
90  if(!pPDG)
91  {
92 #ifdef pdebug
93  G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
94 #endif
95  return 0.; // projectile PDG=0 is a mistake (?!) @@
96  }
97  G4bool in=false; // By default the isotope must be found in the AMDB
98  if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
99  {
100  in = false; // By default the isotope haven't be found in AMDB
101  lastP = 0.; // New momentum history (nothing to compare with)
102  lastPDG = pPDG; // The last PDG of the projectile
103  lastN = tN; // The last N of the calculated nucleus
104  lastZ = tZ; // The last Z of the calculated nucleus
105  lastI = colN.size(); // Size of the Associative Memory DB in the heap
106  j = 0; // A#0f records found in DB for this projectile
107 #ifdef pdebug
108  G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl;
109 #endif
110  if(lastI) for(G4int i=0; i<lastI; i++) // Loop over all DB
111  { // The nucleus with projPDG is found in AMDB
112 #ifdef pdebug
113  G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl;
114 #endif
115  if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ)
116  {
117  lastI=i;
118  lastTH =colTH[i]; // Last THreshold (A-dependent)
119 #ifdef pdebug
120  G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
121 #endif
122  if(pMom<=lastTH)
123  {
124 #ifdef pdebug
125  G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl;
126 #endif
127  return 0.; // Energy is below the Threshold value
128  }
129  lastP =colP [i]; // Last Momentum (A-dependent)
130  lastICS=colICS[i]; // Last Inelastic Cross-Section (A-dependent)
131  lastECS=colECS[i]; // Last Elastic Cross-Section (A-dependent)
132  if(std::fabs(lastP/pMom-1.)<tolerance)
133  {
134 #ifdef pdebug
135  G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS="
136  <<lastECS*millibarn<<G4endl;
137 #endif
138  CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
139  if(fCS) return lastICS*millibarn; // Use theLastInelasticCS
140  return lastECS*millibarn; // Use theLastElasticCS
141  }
142  in = true; // This is the case when the isotop is found in DB
143  // Momentum pMom is in IU ! @@ Units
144 #ifdef pdebug
145  G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
146 #endif
147  lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
148  lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
149 #ifdef pdebug
150  G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl;
151 #endif
152  if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold
153  {
154 #ifdef pdebug
155  G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
156 #endif
157  lastTH=pMom;
158  }
159  break; // Go out of the LOOP
160  }
161 #ifdef pdebug
162  G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
163  <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl;
164 #endif
165  j++; // Increment a#0f records found in DB for this pPDG
166  }
167  if(!in) // This nucleus has not been calculated previously
168  {
169 #ifdef pdebug
170  G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
171 #endif
172  lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
174  lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
175  if(lastICS<=0. || lastECS<=0.)
176  {
177  lastTH = ThresholdEnergy(tZ, tN); // Threshold Energy=Mom=0 which is now the last
178 #ifdef pdebug
179  G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
180 #endif
181  if(pMom>lastTH)
182  {
183 #ifdef pdebug
184  G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl;
185 #endif
186  lastTH=pMom;
187  }
188  }
189 #ifdef pdebug
190  G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z="
191  <<lastZ<<G4endl;
192 #endif
193  colN.push_back(tN);
194  colZ.push_back(tZ);
195  colPDG.push_back(pPDG);
196  colP.push_back(pMom);
197  colTH.push_back(lastTH);
198  colICS.push_back(lastICS);
199  colECS.push_back(lastECS);
200 #ifdef pdebug
201  G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn
202  <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
203 #endif
204  if(fCS) return lastICS*millibarn; // Use theLastInelasticCS
205  return lastECS*millibarn; // Use theLastElasticCS
206  } // End of creation of the new set of parameters
207  else
208  {
209 #ifdef pdebug
210  G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
211 #endif
212  colP[lastI]=pMom;
213  colPDG[lastI]=pPDG;
214  colICS[lastI]=lastICS;
215  colECS[lastI]=lastECS;
216  }
217  } // End of parameters udate
218  else if(pMom<=lastTH)
219  {
220 #ifdef pdebug
221  G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
222 #endif
223  return 0.; // Momentum is below the Threshold Value -> CS=0
224  }
225  else if(std::fabs(lastP/pMom-1.)<tolerance)
226  {
227 #ifdef pdebug
228  G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn
229  <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
230 #endif
231  if(fCS) return lastICS*millibarn; // Use theLastInelasticCS
232  return lastECS*millibarn; // Use theLastElasticCS
233  }
234  else
235  {
236 #ifdef pdebug
237  G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
238 #endif
239  lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
240  lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
241  lastP=pMom;
242  }
243 #ifdef pdebug
244  G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS="
245  <<lastECS*millibarn<<"(mb)"<<G4endl;
246 #endif
247  if(fCS) return lastICS*millibarn; // Use theLastInelasticCS
248  return lastECS*millibarn; // Use theLastElasticCS
249 }
250 
251 // The main member function giving the A-A cross section (Momentum in MeV, CS in mb)
252 G4double G4QIonIonCrossSection::CalculateCrossSection(G4bool XS,G4int F,G4int I,G4int pPDG,
253  G4int tZ,G4int tN, G4double TotMom)
254 {
255  //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)]
256  //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21
257  // If change the following, please change in ::GetFunctions:
258  static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold
259  static const G4double dP=10.; // step for the LEN table
260  static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106)
261  static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part
262  static const G4double Pmax=300000.; // maxE for the HighE part
263  static const G4int nH=100; // A#of HResonance points in lnE
264  static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part
265  static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent)
266  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
267  //
268  // Associative memory for acceleration
269  static std::vector <G4double*> LENI; // Vector of pointers: LowEnIneIonIonCrossSection
270  static std::vector <G4double*> HENI; // Vector of pointers: HighEnIneIonIonCrossSection
271  static std::vector <G4double*> LENE; // Vector of pointers: LowEnElaIonIonCrossSection
272  static std::vector <G4double*> HENE; // Vector of pointers: HighEnElaIonIonCrossSection
273 #ifdef debug
274  G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl;
275 #endif
276  G4int dPDG=pPDG/10; // 10SZZZAAA
277  G4int zPDG=dPDG/1000; // 10SZZZ (?)
278  G4int zA=dPDG%1000; // proj A
279  G4int pZ=zPDG%1000; // proj Z (?)
280  G4int pN=zA-pZ; // proj N (?)
281  G4double Momentum=TotMom/zA; // Momentum per nucleon
282  if (Momentum<THmin) return 0.; // @@ This can be dangerouse for the heaviest nuc.!
283  G4double sigma=0.;
284  if(F&&I) sigma=0.; // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
285  G4double tA=tN+tZ; // Target weight
286  G4double pA=zA; // Projectile weight
287  if(F<=0) // This isotope was not the last used isotop
288  {
289  if(F<0 || !XS) // This isotope was found in DAMDB or Elast =>RETRIEVE
290  {
291  lastLENI=LENI[I]; // Pointer to Low Energy inelastic cross sections
292  lastHENI=HENI[I]; // Pointer to High Energy inelastic cross sections
293  lastLENE=LENE[I]; // Pointer to Low Energy inelastic cross sections
294  lastHENE=HENE[I]; // Pointer to High Energy inelastic cross sections
295  }
296  else // This isotope wasn't calculated previously => CREATE
297  {
298  lastLENI = new G4double[nL]; // Allocate memory for the new LEN cross sections
299  lastHENI = new G4double[nH]; // Allocate memory for the new HEN cross sections
300  lastLENE = new G4double[nL]; // Allocate memory for the new LEN cross sections
301  lastHENE = new G4double[nH]; // Allocate memory for the new HEN cross sections
302  G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE);
303  if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl;
304 #ifdef debug
305  G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA
306  <<G4endl;
307 #endif
308  // *** The synchronization check ***
309  G4int sync=LENI.size();
310  if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl;
311  LENI.push_back(lastLENI); // added LEN Inelastic
312  HENI.push_back(lastHENI); // added HEN Inelastic
313  LENE.push_back(lastLENE); // added LEN Elastic
314  HENE.push_back(lastHENE); // added HEN Elastic
315  } // End of creation of the new set of parameters
316  } // End of parameters udate
317  // =------------= NOW the Magic Formula =--------------------------=
318  if (Momentum<lastTH) return 0.; // It must be already checked in the interface class
319  else if (Momentum<Pmin) // LEN region (approximated in E, not in lnE)
320  {
321 #ifdef debug
322  G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl;
323 #endif
324  if(tA<1. || pA<1.)
325  {
326  G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl;
327  sigma=0.;
328  }
329  else
330  {
331  G4double dPp=dP*pA;
332  if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI);
333  else sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE);
334  }
335 #ifdef debugn
336  if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P="
337  <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl;
338 #endif
339  }
340  else if (Momentum<Pmax*pA) // High Energy region
341  {
342  G4double lP=std::log(Momentum);
343 #ifdef debug
344  G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
345 #endif
346  if(tA<1. || pA<1.)
347  {
348  G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl;
349  sigma=0.;
350  }
351  else
352  {
353  G4double milPp=milP+std::log(pA);
354  if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI);
355  else sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE);
356  }
357  }
358  else // UltraHighE region (not frequent)
359  {
360  std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum);
361  if(XS) sigma=inelel.first;
362  else sigma=inelel.second;
363  }
364 #ifdef debug
365  G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl;
366 #endif
367  if(sigma<0.) return 0.;
368  return sigma;
369 }
370 
371 // Linear fit for YN[N] tabulated (from X0 with fixed step DX) function to X point
372 
373 // Calculate the functions for the log(A)
374 G4int G4QIonIonCrossSection::GetFunctions(G4int pZ,G4int pN,G4int tZ,G4int tN,G4double* li,
375  G4double* hi, G4double* le, G4double* he)
376 {
377  // If change the following, please change in ::CalculateCrossSection:
378  static const G4double THmin=0.; // @@ start from threshold (?) minimum Energy Threshold
379  static const G4double dP=10.; // step for the LEN table
380  static const G4int nL=100; // A#of LENesonance points in E (each MeV from 2 to 106)
381  static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part
382  static const G4double Pmax=300000.; // maxE for the HighE part
383  static const G4int nH=100; // A#of HResonance points in lnE
384  static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part
385  static const G4double malP=std::log(Pmax); // High logarithm energy
386  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
387  static const G4double lP=std::exp(dlP); // Multiplication factor in the HighE part
388  // If the cross section approximation formula is changed - replace from file.
389  if(pZ<1 || pN<0 || tZ<1 || tN<0)
390  {
391  G4cout<<"-W-G4QIonIonCS::GetFunct:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
392  return -1;
393  }
394  G4int pA=pN+pZ;
395  G4double dPp=dP*pA;
396  G4double Mom=THmin;
397  for(G4int k=0; k<nL; k++)
398  {
399  std::pair<G4double,G4double> len = CalculateXS(pZ, pN, tZ, tN, Mom);
400  li[k]=len.first;
401  le[k]=len.second;
402  Mom+=dPp;
403  }
404  G4double lMom=Pmin*pA;
405  for(G4int j=0; j<nH; j++)
406  {
407  std::pair<G4double,G4double> len = CalculateXS(pZ, pN, pZ, pN, lMom);
408  hi[j]=len.first;
409  he[j]=len.second;
410  lMom*=lP;
411  }
412 #ifdef debug
413  G4cout<<"G4QIonIonCS::GetFunctions: pZ="<<pZ<<", tZ="<<tZ<<" pair is calculated"<<G4endl;
414 #endif
415  return 1;
416 }
417 
418 // Momentum (Mom=p/A) is in MeV/c, first=InelasticXS, second=ElasticXS (mb)
419 std::pair<G4double,G4double> G4QIonIonCrossSection::CalculateXS(G4int pZ,G4int pN,G4int tZ,
420  G4int tN, G4double Mom)
421 {
422  static G4VQCrossSection* PElCSman = G4QProtonElasticCrossSection::GetPointer();
423  static G4VQCrossSection* NElCSman = G4QNeutronElasticCrossSection::GetPointer();
424  static G4VQCrossSection* InelPCSman = G4QProtonNuclearCrossSection::GetPointer();
425  static G4VQCrossSection* InelNCSman = G4QNeutronNuclearCrossSection::GetPointer();
426  G4double pA=pZ+pN;
427  G4double tA=tZ+tN;
428  if(pA<.9 || tA<.9 ||pA>239. || tA>239 || Mom < 0.) return std::make_pair(0.,0.);
429  G4double inCS=0.;
430  G4double elCS=0.;
431  if(pA<1.1 ) // nucleon-ion interaction use NA(in,el)
432  {
433  if (pZ == 1 && !pN) // proton-nuclear
434  {
435  inCS=InelPCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
436  elCS=PElCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
437  }
438  else if(pN == 1 && !pZ) // neutron-nuclear
439  {
440  inCS=InelNCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
441  elCS=NElCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
442  }
443  else G4cerr<<"-Warn-G4QIICS::CaCS:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
444  }
445  else
446  {
447  G4double T=ThresholdMomentum(pZ, pN, tZ, tN); // @@ Can be cashed as lastTH (?)
448  if(Mom<=T)
449  {
450  elCS=0.;
451  inCS=0.;
452  }
453  else
454  {
455  G4double P2=Mom*Mom;
456  G4double T2=T*T;
457  G4double R=1.-T2/P2; // @@ Very rough threshold effect
458  //G4double P4=P2*P2;
459  //G4double P8=P4*P4;
460  //G4double T4=T2*T2;
461  //G4double tot=CalculateTotal(pA, tA, Mom)*P8/(P8+T4*T4); // @@ convert to IndepUnits
462  G4double tot=R*CalculateTotal(pA, tA, Mom); // @@ convert to IndepUnits
463  G4double rat=CalculateElTot(pA, tA, Mom);
464  elCS=tot*rat;
465  inCS=tot-elCS;
466  }
467  }
468  return std::make_pair(inCS,elCS);
469 }
470 
471 // Total Ion-ion cross-section (mb), Momentum (Mom) here is p/A (MeV/c=IU) (No Threshold)
472 G4double G4QIonIonCrossSection::CalculateTotal(G4double pA, G4double tA, G4double Mom)
473 {
474  G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c
475  G4double ab=pA+tA;
476  G4double al=std::log(ab);
477  G4double ap=std::log(pA*tA);
478  G4double e=std::pow(pA,0.1)+std::pow(tA,0.1);
479  G4double d=e-1.55/std::pow(al,0.2);
480  G4double f=4.;
481  if(pA>4. && tA>4.) f=3.3+130./ab/ab+2.25/e;
482  G4double c=(385.+11.*ab/(1.+.02*ab*al)+(.5*ab+500./al/al)/al)*std::pow(d,f);
483  G4double r=y-3.-4./ap/ap;
484 #ifdef pdebug
485  G4cout<<"G4QIonIonCS::CalcTot:P="<<Mom<<", stot="<<c+d*r*r<<", d="<<d<<", r="<<r<<G4endl;
486 #endif
487  return c+d*r*r;
488 }
489 
490 // Ratio elastic/Total, Momentum (Mom) here is p/A (MeV/c=IU)
491 G4double G4QIonIonCrossSection::CalculateElTot(G4double pA, G4double tA, G4double Mom)
492 {
493  G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c
494  G4double ab=pA*tA;
495  G4double ap=std::log(ab);
496  G4double r=y-3.92-1.73/ap/ap;
497  G4double d=.00166/(1.+.002*std::pow(ab,1.33333));
498  G4double al=std::log(pA+tA);
499  G4double e=1.+.235*(std::fabs(ap-5.78));
500  if (std::fabs(pA-2.)<.1 && std::fabs(tA-2.)<.1) e=2.18;
501  else if(std::fabs(pA-2.)<.1 && std::fabs(tA-3.)<.1) e=1.90;
502  else if(std::fabs(pA-3.)<.1 && std::fabs(tA-2.)<.1) e=1.90;
503  else if(std::fabs(pA-2.)<.1 && std::fabs(tA-4.)<.1) e=1.65;
504  else if(std::fabs(pA-4.)<.1 && std::fabs(tA-2.)<.1) e=1.65;
505  else if(std::fabs(pA-3.)<.1 && std::fabs(tA-4.)<.1) e=1.32;
506  else if(std::fabs(pA-4.)<.1 && std::fabs(tA-3.)<.1) e=1.32;
507  else if(std::fabs(pA-4.)<.1 && std::fabs(tA-4.)<.1) e=1.;
508  G4double f=.37+.0188*al;
509  G4double g_value=std::log(std::pow(pA,0.35)+std::pow(tA,0.35));
510  G4double h=g_value*g_value;
511  G4double c=f/(1.+e/h/h);
512 #ifdef pdebug
513  G4cout<<"G4QIonIonCS::CalcElT:P="<<Mom<<",el/tot="<<c+d*r*r<<",d="<<d<<", r="<<r<<G4endl;
514 #endif
515  return c+d*r*r;
516 }
517 
518 // Electromagnetic momentum/A-threshold (in MeV/c)
519 G4double G4QIonIonCrossSection::ThresholdMomentum(G4int pZ, G4int pN, G4int tZ, G4int tN)
520 {
521  static const G4double third=1./3.;
522  static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
523  static const G4double tpM= pM+pM; // Doubled proton mass (MeV)
524  if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.;
525  G4double tA=tZ+tN;
526  G4double pA=pZ+pN;
527  //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
528  G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA (per projNucleon)
529  return std::sqrt(dE*(tpM+dE));
530 }
Generated on Sat May 25 2013 14:33:48 for Geant4 by   doxygen 1.8.4