source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QProtonNuclearCrossSection.cc @ 1197

Last change on this file since 1197 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 19.2 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) 17-June-02
28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
29//
30//
31// G4 Physics class: G4QProtonNuclearCrossSection 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// Short description: CHIPS cross-sections for proton-nuclear interactions
40// -----------------------------------------------------------------------
41//
42//#define debug
43//#define pdebug
44//#define debug3
45//#define debugn
46//#define debugs
47
48#include "G4QProtonNuclearCrossSection.hh"
49
50// Initialization of the
51G4double* G4QProtonNuclearCrossSection::lastLEN=0; // Pointer to the lastArray of LowEn CS
52G4double* G4QProtonNuclearCrossSection::lastHEN=0; // Pointer to the lastArray of HighEn CS
53G4int     G4QProtonNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
54G4int     G4QProtonNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
55G4double  G4QProtonNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
56G4double  G4QProtonNuclearCrossSection::lastTH=0.; // Last threshold momentum
57G4double  G4QProtonNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
58G4int     G4QProtonNuclearCrossSection::lastI=0;   // The last position in the DAMDB
59
60// Returns Pointer to the G4VQCrossSection class
61G4VQCrossSection* G4QProtonNuclearCrossSection::GetPointer()
62{
63  static G4QProtonNuclearCrossSection theCrossSection; //**Static body of Cross Section**
64  return &theCrossSection;
65}
66
67// The main member function giving the collision cross section (P is in IU, CS is in mb)
68// Make pMom in independent units ! (Now it is MeV)
69G4double G4QProtonNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
70                                                       G4int tgZ, G4int tgN, G4int)
71{
72  static G4double tolerance=0.001;     // Tolerance (0.1%) to consider as "the same mom"
73  static G4int j;                      // A#0f Z/N-records already tested in AMDB
74  static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
75  static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
76  static std::vector <G4double> colP;  // Vector of last momenta for the reaction
77  static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
78  static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
79  // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
80#ifdef debug
81  G4cout<<"G4QPrCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
82        <<"("<<lastN<<"),PDG=2212, thresh="<<lastTH<<",Sz="<<colN.size()<<G4endl;
83#endif
84  G4bool in=false;                     // By default the isotope must be found in the AMDB
85  if(tgN!=lastN || tgZ!=lastZ)         // The nucleus was not the last used isotope
86  {
87    in = false;                        // By default the isotope haven't be found in AMDB 
88    lastP   = 0.;                      // New momentum history (nothing to compare with)
89    lastN   = tgN;                     // The last N of the calculated nucleus
90    lastZ   = tgZ;                     // The last Z of the calculated nucleus
91    lastI   = colN.size();             // Size of the Associative Memory DB in the heap
92    j  = 0;                            // A#0f records found in DB for this projectile
93#ifdef debug
94    G4cout<<"G4QPrCS::GetCS: the amount of records in the AMDB lastI="<<lastI<<G4endl;
95#endif
96    if(lastI) for(G4int i=0; i<lastI; i++) // AMDB exists, try to find the (Z,N) isotope
97    {
98      if(colN[i]==tgN && colZ[i]==tgZ) // Try the record "i" in the AMDB
99      {
100        lastI=i;                       // Remember the index for future fast/last use
101        lastTH =colTH[i];              // The last THreshold (A-dependent)
102#ifdef debug
103        G4cout<<"G4QPrCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
104#endif
105        if(pMom<=lastTH)
106        {
107#ifdef debug
108          G4cout<<"G4QPCS::GetCS:Found,P="<<pMom<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
109#endif
110          return 0.;                   // Energy is below the Threshold value
111        }
112        lastP  =colP [i];              // Last Momentum  (A-dependent)
113        lastCS =colCS[i];              // Last CrossSect (A-dependent)
114        if(std::fabs(lastP-pMom)<tolerance*pMom)
115        //if(lastP==pMom)              // VI do not use tolerance
116        {
117#ifdef debug
118          G4cout<<"..G4QPrCS::GetCS:.DoNothing.P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
119#endif
120          //CalculateCrossSection(fCS,-1,j,2212,lastZ,lastN,pMom); // Update param's only
121          return lastCS*millibarn;     // Use theLastCS
122        }
123        in = true;                     // This is the case when the isotop is found in DB
124        // Momentum pMom is in IU ! @@ Units
125#ifdef debug
126        G4cout<<"G4QPrCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
127#endif
128        lastCS=CalculateCrossSection(fCS,-1,j,2212,lastZ,lastN,pMom); // read & update
129#ifdef debug
130        G4cout<<"G4QPrCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
131#endif
132        if(lastCS<=0. && pMom>lastTH)  // Correct the threshold (@@ No intermediate Zeros)
133        {
134#ifdef debug
135          G4cout<<"G4QPrCS::GetCS: New P="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
136#endif
137          lastCS=0.;
138          lastTH=pMom;
139        }
140        break;                         // Go out of the LOOP
141      }
142#ifdef debug
143      G4cout<<"-->G4QPrCrossSec::GetCrosSec: pPDG=2212, j="<<j<<", N="<<colN[i]
144            <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
145#endif
146      j++;                             // Increment a#0f records found in DB
147    }
148#ifdef debug
149    G4cout<<"-?-G4QPrCS::GetCS:RC Z="<<tgZ<<",N="<<tgN<<",in="<<in<<",j="<<j<<" ?"<<G4endl;
150#endif
151    if(!in)                            // This isotope has not been calculated previously
152    {
153#ifdef debug
154      G4cout<<"^^^G4QPrCS::GetCS:CalcNew P="<<pMom<<", f="<<fCS<<", lastI="<<lastI<<G4endl;
155#endif
156      //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
157      lastCS=CalculateCrossSection(fCS,0,j,2212,lastZ,lastN,pMom); //calculate & create
158      if(lastCS>0.)                   // It means that the AMBD was initialized
159      {
160
161        lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
162#ifdef debug
163        G4cout<<"G4QPrCrossSection::GetCrossSect: NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
164#endif
165        colN.push_back(tgN);
166        colZ.push_back(tgZ);
167        colP.push_back(pMom);
168        colTH.push_back(lastTH);
169        colCS.push_back(lastCS);
170#ifdef debug
171        G4cout<<"G4QPrCS::GetCrosSec:recCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
172#endif
173      }
174#ifdef pdebug
175      G4cout<<"G4QPrCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
176#endif
177      return lastCS*millibarn;
178    } // End of creation of the new set of parameters
179    else
180    {
181#ifdef debug
182      G4cout<<"G4QPrCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
183#endif
184      colP[lastI]=pMom;
185      colCS[lastI]=lastCS;
186    }
187  } // End of parameters udate
188  else if(pMom<=lastTH)
189  {
190#ifdef debug
191    G4cout<<"G4QPrCS::GetCS: Current P="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
192#endif
193    return 0.;                         // Momentum is below the Threshold Value -> CS=0
194  }
195  else if(std::fabs(lastP-pMom)<tolerance*pMom)
196  //else if(lastP==pMom)               // VI do not use tolerance
197  {
198#ifdef debug
199    G4cout<<"..G4QPCS::GetCS:OldNZ&P="<<lastP<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
200#endif
201    return lastCS*millibarn;           // Use theLastCS
202  }
203  else                                 // It is the last used -> use the current tables
204  {
205#ifdef debug
206    G4cout<<"-!-G4QPCS::GetCS:UseCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
207#endif
208    lastCS=CalculateCrossSection(fCS,1,j,2212,lastZ,lastN,pMom); // Only read and UpdateDB
209    lastP=pMom;
210  }
211#ifdef debug
212  G4cout<<"==>G4QPrCS::GetCroSec: P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
213#endif
214  return lastCS*millibarn;
215}
216
217// The main member function giving the gamma-A cross section (E in GeV, CS in mb)
218G4double G4QProtonNuclearCrossSection::CalculateCrossSection(G4bool, G4int F, G4int I,
219                                        G4int, G4int targZ, G4int targN, G4double Momentum)
220{
221  static const G4double THmin=27.;     // default minimum Momentum (MeV/c) Threshold
222  static const G4double THmiG=THmin*.001; // minimum Momentum (GeV/c) Threshold
223  static const G4double dP=10.;        // step for the LEN (Low ENergy) table MeV/c
224  static const G4double dPG=dP*.001;   // step for the LEN (Low ENergy) table GeV/c
225  static const G4int    nL=105;        // A#of LEN points in E (step 10 MeV/c)
226  static const G4double Pmin=THmin+(nL-1)*dP; // minP for the HighE part with safety
227  static const G4double Pmax=227000.;  // maxP for the HEN (High ENergy) part 227 GeV
228  static const G4int    nH=224;        // A#of HEN points in lnE
229  static const G4double milP=std::log(Pmin);// Low logarithm energy for the HEN part
230  static const G4double malP=std::log(Pmax);// High logarithm energy (each 2.75 percent)
231  static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HEN part
232  static const G4double milPG=std::log(.001*Pmin);// Low logarithmEnergy for HEN part GeV/c
233  //
234  // Associative memory for acceleration
235  //static std::vector <G4double>  spA;  // shadowing coefficients (A-dependent)
236  static std::vector <G4double*> LEN;  // Vector of pointers to LowEnProtonCrossSection
237  static std::vector <G4double*> HEN;  // Vector of pointers to HighEnProtonCrossSection
238#ifdef debug
239  G4cout<<"G4QProtNCS::CalCS:N="<<targN<<",Z="<<targZ<<",P="<<Momentum<<">"<<THmin<<G4endl;
240#endif
241  if (Momentum<THmin) return 0.;       // @@ This can be dangerouse for the heaviest nuc.?!
242  G4double sigma=0.;
243  if(F&&I) sigma=0.;                   // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
244  G4double A=targN+targZ;              // A of the target
245#ifdef debug
246  G4cout<<"G4QProtNucCS::CalCS: A="<<A<<",F="<<F<<",I="<<I<<",nL="<<nL<<",nH="<<nH<<G4endl;
247#endif
248  if(F<=0)                             // This isotope was not the last used isotop
249  {
250    if(F<0)                            // This isotope was found in DAMDB =======> RETRIEVE
251    {
252      G4int sync=LEN.size();
253      if(sync<=I) G4cerr<<"*!*G4QProtonNuclCS::CalcCrossSect:Sync="<<sync<<"<="<<I<<G4endl;
254      lastLEN=LEN[I];                  // Pointer to prepared LowEnergy cross sections
255      lastHEN=HEN[I];                  // Pointer to prepared High Energy cross sections
256    }
257    else if(Momentum<ThresholdMomentum(targZ,targN)) return 0.; // BelowThreshold -> NotIni
258    else                               // This isotope wasn't calculated before => CREATE
259    {
260      lastLEN = new G4double[nL];      // Allocate memory for the new LEN cross sections
261      lastHEN = new G4double[nH];      // Allocate memory for the new HEN cross sections
262      // --- Instead of making a separate function ---
263      G4double P=THmiG;                // Table threshold in GeV/c
264      for(G4int m=0; m<nL; m++)
265      {
266        lastLEN[m] = CrossSectionLin(targZ, targN, P);
267        P+=dPG;
268      }
269      G4double lP=milPG;
270      for(G4int n=0; n<nH; n++)
271      {
272        lastHEN[n] = CrossSectionLog(targZ, targN, lP);
273        lP+=dlP;
274      }
275#ifdef debug
276      G4cout<<"-*->G4QPr0tNucCS::CalcCS:Tab for Z="<<targZ<<",N="<<targN<<",I="<<I<<G4endl;
277#endif
278      // --- End of possible separate function
279      // *** The synchronization check ***
280      G4int sync=LEN.size();
281      if(sync!=I)
282      {
283        G4cerr<<"***G4QProtonNuclCS::CalcCrossSect: Sinc="<<sync<<"#"<<I<<", Z=" <<targZ
284              <<", N="<<targN<<", F="<<F<<G4endl;
285        //G4Exception("G4ProtonNuclearCS::CalculateCS:","39",FatalException,"overflow DB");
286      }
287      LEN.push_back(lastLEN);          // remember the Low Energy Table
288      HEN.push_back(lastHEN);          // remember the High Energy Table
289    } // End of creation of the new set of parameters
290  } // End of parameters udate
291  // ============================== NOW the Magic Formula =================================
292#ifdef debug
293  G4cout<<"G4QPrNCS::CalcCS:lTH="<<lastTH<<",Pmi="<<Pmin<<",dP="<<dP<<",dlP="<<dlP<<G4endl;
294#endif
295  if (Momentum<lastTH) return 0.;      // It must be already checked in the interface class
296  else if (Momentum<Pmin)              // High Energy region
297  {
298#ifdef debug
299    G4cout<<"G4QPrNCS::CalcCS:bLEN A="<<A<<", nL="<<nL<<",TH="<<THmin<<",dP="<<dP<<G4endl;
300#endif
301    if(A<=1.) sigma=0.;
302    else      sigma=EquLinearFit(Momentum,nL,THmin,dP,lastLEN);
303#ifdef debugn
304    if(sigma<0.)
305      G4cout<<"G4QPrNuCS::CalcCS:A="<<A<<",E="<<Momentum<<",T="<<THmin<<",dP="<<dP<<G4endl;
306#endif
307  }
308  else if (Momentum<Pmax)              // High Energy region
309  {
310    G4double lP=std::log(Momentum);
311#ifdef debug
312    G4cout<<"G4QProtNucCS::CalcCS: before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
313#endif
314    sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN);
315  }
316  else                                 // UHE region (calculation, not frequent)
317  {
318    G4double P=0.001*Momentum;         // Approximation formula is for P in GeV/c
319    sigma=CrossSectionFormula(targZ, targN, P, std::log(P));
320  }
321#ifdef debug
322  G4cout<<"G4QProtonNuclearCrossSection::CalcCS: CS="<<sigma<<G4endl;
323#endif
324  if(sigma<0.) return 0.;
325  return sigma;
326}
327
328// Electromagnetic momentum-threshold (in MeV/c)
329G4double G4QProtonNuclearCrossSection::ThresholdMomentum(G4int tZ, G4int tN)
330{
331  static const G4double third=1./3.;
332  static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
333  static const G4double tpM= pM+pM;       // Doubled proton mass (MeV)
334  G4double tA=tZ+tN;
335  if(tZ<.99 || tN<0.) return 0.;
336  else if(tZ==1 && tN==0) return 800.;    // A threshold on the free proton
337  //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
338  G4double dE=tZ/(1.+std::pow(tA,third)); // Safety for diffused edge of the nucleus (QE)
339  return std::sqrt(dE*(tpM+dE));
340}
341
342// Calculation formula for proton-nuclear inelastic cross-section (mb) (P in GeV/c)
343G4double G4QProtonNuclearCrossSection::CrossSectionLin(G4int tZ, G4int tN, G4double P)
344{
345  G4double sigma=0.;
346  if(P<ThresholdMomentum(tZ,tN)*.001) return sigma;
347  G4double lP=std::log(P);
348  if(tZ==1&&!tN){if(P>.35) sigma=CrossSectionFormula(tZ,tN,P,lP);}// s(pp)=0 below 350Mev/c
349  else if(tZ<93 && tN<239)                // General solution
350  {
351    G4double pex=0.;
352    G4double pos=0.;
353    G4double wid=1.;
354    if(tZ==13 && tN==14)                  // Excited metastable states
355    {
356      pex=230.;
357      pos=.13;
358      wid=8.e-5;
359    }
360    else if(tZ<7)
361    {
362      if(tZ==6 && tN==6)
363      {
364        pex=320.;
365        pos=.14;
366        wid=7.e-6;
367      }
368      else if(tZ==4 && tN==5)
369      {
370        pex=600.;
371        pos=.132;
372        wid=.005;
373      }
374      else if(tZ==3 && tN==4)
375      {
376        pex=280.;
377        pos=.19;
378        wid=.0025;
379      }
380      else if(tZ==3 && tN==3)
381      {
382        pex=370.;
383        pos=.171;
384        wid=.006;
385      }
386      else if(tZ==2 && tN==1)
387      {
388        pex=30.;
389        pos=.22;
390        wid=.0005;
391      }
392    }
393    sigma=CrossSectionFormula(tZ,tN,P,lP);
394    if(pex>0.)
395    {
396      G4double dp=P-pos;
397      sigma+=pex*std::exp(-dp*dp/wid);
398    }
399  }
400  else
401  {
402    G4cerr<<"-Warning-G4QProtonNuclearCroSect::CSLin:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
403    sigma=0.;
404  }
405  if(sigma<0.) return 0.;
406  return sigma; 
407}
408
409// Calculation formula for proton-nuclear inelastic cross-section (mb) log(P in GeV/c)
410G4double G4QProtonNuclearCrossSection::CrossSectionLog(G4int tZ, G4int tN, G4double lP)
411{
412  G4double P=std::exp(lP);
413  return CrossSectionFormula(tZ, tN, P, lP);
414}
415// Calculation formula for proton-nuclear inelastic cross-section (mb) log(P in GeV/c)
416G4double G4QProtonNuclearCrossSection::CrossSectionFormula(G4int tZ, G4int tN,
417                                                           G4double P, G4double lP)
418{
419  G4double sigma=0.;
420  if(tZ==1 && !tN)                        // pp interaction
421  {
422    G4double sp=std::sqrt(P);
423    G4double ds=lP-4.2;
424    G4double dp=P-.35;
425    G4double d3=dp*dp*dp;
426    sigma=(33.+.2*ds*ds)/(1.+.4/sp)/(1.+.5/d3/d3);
427  }
428  else if(tZ<93 && tN<146)                // General solution
429  {
430    //G4double lP=std::log(P);            // Already calculated
431    G4double d=lP-4.2;
432    G4double p2=P*P;
433    G4double p4=p2*p2;
434    G4double a=tN+tZ;                       // A of the target
435    G4double al=std::log(a);
436    G4double a2=a*a;
437    G4double a4=a2*a2;
438    G4double a8=a4*a4;
439    G4double a12=a8*a4;
440    G4double a16=a8*a8;
441    G4double c=170./(1.+5./std::exp(al*1.6));
442    G4double dl=al-2.8;
443    G4double r=.3+.2*dl*dl;
444    G4double g=40.*std::exp(al*0.712)/(1.+12.2/a)/(1.+34./a2);
445    G4double e=318.+a4/(1.+.0015*a4/std::exp(al*0.09))/(1.+4.e-28*a12)+
446               8.e-18/(1./a16+1.3e-20)/(1.+1.e-21*a12);
447    G4double s=3.57+.009*a2/(1.+.0001*a2*a);
448    G4double h=(.01/a4+2.5e-6/a)*(1.+7.e-8*a4)/(1.+6.e7/a12/a2);
449    sigma=(c+d*d)/(1.+r/p4)+(g+e*std::exp(-s*P))/(1.+h/p4/p4);
450#ifdef pdebug
451  G4cout<<"G4QProtNucCS::CSForm: A="<<a<<",P="<<P<<",CS="<<sigma<<",c="<<c<<",g="<<g<<",d="
452        <<d<<",r="<<r<<",e="<<e<<",h="<<h<<G4endl;
453#endif
454  }
455  else
456  {
457    G4cerr<<"-Warning-G4QProtonNuclearCroSect::CSForm:*Bad A* Z="<<tZ<<", N="<<tN<<G4endl;
458    sigma=0.;
459  }
460  if(sigma<0.) return 0.;
461  return sigma; 
462}
Note: See TracBrowser for help on using the repository browser.