source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QNeutronNuclearCrossSection.cc @ 1196

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

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