source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QIonIonCrossSection.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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