source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QIonIonCrossSection.cc @ 968

Last change on this file since 968 was 962, checked in by garnier, 15 years ago

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