source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4PhotoNuclearCalculation.cc @ 962

Last change on this file since 962 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 9.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// G4 Tools program: PhotoNuclearCalculation of gamma+A cross sections
28// according to the paper M.V.Kossov "Approximation of photonuclear
29// interaction cross sections" (EPJA,2002, be published). The cross
30// sections are calculated according to the exact complicated formulas
31// and tabulated for the fast response. The total energy range is
32// subdivided in three regions:
33// a) The GDR region (from the hadron production threshold to 106 MeV)
34//    covering 46 nuclei: H2,He4,Li6,Li7,Be9,C12,N14,N15,O16,F19,Na23,
35//    Mg24,Al27,Si28,S32,S34,(Ar40,Ca40),Fe54,Mn55,Fe56,Ni58,Co59,Cu,Zn
36//    Se76,Se82,Ag,Cd,Sn,I,Sm154,Gd156,Tb159,Ho165,Er168,Yb174,Hf178,Hf180,
37//    Ta181,W184,W186,Au197,Tl,Pb,Bi209,Th232,U235,U238,Pu239. For these
38//    isotops the approximation fits the existing measurements. For
39//    them the basic functions are calculated. For other isotops the
40//    approximation is interpolated linearly with A.
41// b) The resonance region (from 106 MeV to 50 GeV)
42//    covering 14 nuclei: H1,H2,He3,He4,Li6,Li7,Be,C,O,Al,Cu,Sn,Pb,U.
43//    For these isotops the melting of resonances in nuclear matter is
44//    measured and approximated. For them the basic functions are
45//    calculated. For other isotops the approximation is interpolated
46//    linearly with A.
47// c) The high energy region is an A-dependent function for E>100 GeV,
48//    taking into account nuclear shaddowing.
49// .....................................................
50// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-May-02
51// The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02
52//
53//=====================================================================
54#include "globals.hh"
55#include <iostream>
56#include <fstream>
57#include <vector>
58#include "G4ios.hh"
59
60int main()
61{
62  const G4int mN=49;                // A#of GDR basic nuclei
63  const G4int mC=105;               // A#of Resonance points in E (each MeV from 2 to 106)
64  const G4double mA[mN]={
65        2.,4.,6.,7.,9.,12.,14.,15.,16.,19.,23.,24.,27.,28.,32.,34.,40.,54.,
66    55.,56.,58.7,58.9,63.5,65.4,76.,82.,107.9,112.4,118.7,126.9,154.,156.,159.,165.,
67    168.,174.,178.,180.,181.,184.,186.,197.,204.4,207.2,209.,232.,235.,
68    238.,239.};
69  const G4int nN=14;                // A#of Resonance basic nuclei
70  const G4int nC=224;               // A#of Resonance points in lnE
71  const G4double nA[nN]={1.,2.,3.,4.,6.,7.,9.,12.,16.,27.,63.5,118.7,207.2,238.};
72  const G4double r4[6]={0.,3.38,3.51,3.23,3.57,3.60};
73  const G4double t4[6]={0.,3.14,3.09,2.80,3.25,3.17};
74  const G4double r8[mN]={
75    0.00,0.00,0.00,2.57,3.26,3.52,3.51,3.51,3.49,3.05,3.20,3.25,3.40,3.42,3.38,3.50,3.51,
76    3.59,3.45,3.46,3.52,3.56,3.40,3.40,3.29,3.38,3.50,3.49,3.41,3.46,3.40,3.37,3.40,3.29,
77    3.37,3.21,3.45,3.45,3.39,3.36,3.36,3.39,3.47,3.40,3.36,3.31,3.30,3.37,3.36};
78  const G4double t8[mN]={
79    0.00,0.00,0.00,2.48,3.05,3.11,3.09,3.08,3.10,2.70,2.90,2.95,3.00,2.97,2.98,3.00,3.05,
80    2.94,2.93,2.86,2.99,3.01,2.85,2.84,2.72,2.76,2.81,2.77,2.72,2.75,2.78,2.74,2.72,2.72,
81    2.70,2.74,2.70,2.70,2.64,2.63,2.61,2.59,2.63,2.57,2.63,2.62,2.62,2.67,2.65};
82  const G4double iE=std::log(106.);      // Start logarithm energy
83  const G4double fE=std::log(50000.);    // Finish logarithm energy (each 2.75 percent)
84  const G4double dE=(fE-iE)/(nC-1); // Step in logarithm energy
85  std::ofstream fileGDR("GDR.out", std::ios::out);
86  fileGDR.setf( std::ios::scientific, std::ios::floatfield );
87  std::ofstream fileRes("Res.out", std::ios::out);
88  fileRes.setf( std::ios::scientific, std::ios::floatfield );
89  G4int np=0;
90  for(G4int m=0; m<mN; m++)
91  {
92    G4double A=mA[m];
93    G4double lnA=std::log(A);
94    G4double A2=A*A;
95    G4double red=1.+16/A2/A2;
96    G4double rho1=(3.2+.75*lnA)/red;
97    G4double tau1=(6.6-.5*lnA)/red;
98        if(A==2.)
99        {
100          rho1=1.86;
101          tau1=1.2;
102        }
103        G4double rho2=(4.+.125*lnA)/red;
104        G4double tau2=3.4/red;
105        if(A==6.)
106        {
107          rho2=2.9;
108          tau2=2.32;
109        }
110        else if(A==2.)
111        {
112          rho2=2.11;
113          tau2=1.5;
114        }
115        G4double rho4=3.8+.05*lnA;
116        G4double tau4=3.8-.25*lnA;
117        if(A<13.)
118        {
119          rho4=r4[m];
120          tau4=t4[m];
121        }
122        G4double rho8=r8[m];
123        G4double tau8=t8[m];
124    G4double tr=5.13-.00075*A; // Resonance threshold Position
125    //if(A==1.)tr=5.24;
126        G4double rr=11.;           // Resonance threshold Slope
127        if(A<2.5)rr=25.;
128        G4double dw=.056+lnA*(.03-.001*lnA); // Delta width
129        G4double du=5.82-.07/(1.+.003*A2);   // Delta position
130    G4double da=.39*A;                   // Delta amplitude
131    if(A==2.)da=.88;
132    //else if(A==1.) da=.55;
133    G4double ekin=1.;
134    fileGDR<<"  static const G4double SL"<<m<<"[nL]={"<<G4endl<<"    ";
135    G4cout<<"**** A_low="<<A<<G4endl;
136    np=0;
137        for(G4int em=0; em<mC; em++)
138        {
139      ekin+=1.;
140      G4double z=std::log(ekin);
141      G4double ds=z-du;
142          G4double sigd=da/(1.+ds*ds/dw)/(1.+std::exp(rr*(tr-z))); // Delta contribution
143          G4double g1=std::exp(rho1-z)/(1.+std::exp(3*(tau1-z)));
144          G4double g2=std::exp(2*(rho2-z))/(1.+std::exp(6*(tau2-z)));
145      G4double g4=0.;
146      if(A>3.5)g4=std::exp(4*(rho4-z))/(1.+std::exp(12*(tau4-z)));
147      G4double g8=0.;
148      if(A>6.5)g8=std::exp(8*(rho8-z))/(1.+std::exp(24*(tau8-z)));
149      G4double sig=sigd+g1+g2+g4+g8;
150      np++;
151      if(np==7)
152      {
153        if(em==mC-1) fileGDR<<sig<<"};"<<G4endl;
154        else         fileGDR<<sig<<","<<G4endl<<"    ";
155          }
156      else      fileGDR<<sig<<",";
157      if(np==7) np=0;
158    } // End of the point LOOP
159  } // End of the isotop LOOP
160  // ===== High Energy Calculations =======================
161  //G4double shd=1.0663-.0023*log(2.);   // HE PomShadowing(D)
162  np=0;
163  for(G4int n=0; n<nN; n++)
164  {
165    G4double A=nA[n];
166    G4double A2=A*A;
167    fileRes<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
168    G4cout<<"**** A_high="<<A<<G4endl;
169    G4double lnA=std::log(A);
170    G4double slA=std::sqrt(lnA);
171    G4double red=1.+16/A2/A2;
172    G4double rho1=(3.2+.75*lnA)/red;
173    G4double tau1=(6.6-.5*lnA)/red;
174        if(A==2.)
175        {
176          rho1=1.86;
177          tau1=1.2;
178        }
179        G4double rho2=(4.+.125*lnA)/red;
180        G4double tau2=3.4/red;
181        if(A==6.)
182        {
183          rho2=2.9;
184          tau2=2.32;
185        }
186        else if(A==2.)
187        {
188          rho2=2.11;
189          tau2=1.5;
190        }
191    G4double rho4=6.27;
192    G4double tau4=7.25;
193    G4double rho8=6.66;
194    G4double tau8=6.90;
195    if(A==2.)
196        {
197     rho4=6.2;
198     tau4=7.1;
199     rho8=6.62;
200     tau8=6.91;
201    }
202    G4double tr=5.13-.00075*A; // Resonance threshold Position
203    if(A==1.)tr=5.24;
204        G4double rr=11.;           // Resonance threshold Slope
205        if(A<2.5)rr=25.;
206        // Delta
207        G4double dw=.056+lnA*(.03-.001*lnA); // Delta width
208        G4double du=5.82-.07/(1.+.003*A2);   // Delta position
209    G4double da=.39*A;                   // Delta amplitude
210    if(A==2.)da=.88;
211    else if(A==1.) da=.55;
212        // High Resonance
213        G4double hw=.045+.04*slA*slA*slA;    // HighR width
214        G4double hu=6.496+.042*lnA;          // HighR position
215    if(A==1.)hu=6.57;
216    else if(A==2.)hu=6.575;
217    G4double ha=.223;                    // HighR amplitude
218    if(A>2.5)ha=.16*A/slA;
219    else if(A==2.) ha=.348;
220    G4double sp=A*(1.-.072*lnA);         // HE TotShadowing
221    G4double sh=1.0663-.0023*lnA;        // HE PomShadowing
222    if(A==1.)sh=1.07;
223        // -- Loop over energies ---
224    G4double z=iE-dE;
225    np=0;
226        for(G4int en=0; en<nC; en++)
227        {
228      z+=dE;
229      G4double ds=z-du;
230      G4double fr=1.+std::exp(rr*(tr-z));
231          G4double sigd=da/(1.+ds*ds/dw); // Delta contribution
232      G4double hs=z-hu;
233          G4double sigh=ha/(1.+hs*hs/hw); // HighR contribution
234          G4double g1=0.;
235      if(A>1.5)g1=std::exp(rho1-z)/(1.+std::exp(3*(tau1-z)));
236          G4double g2=0.;
237      if(A>1.5)g2=std::exp(2*(rho2-z))/(1.+std::exp(6*(tau2-z)));
238      G4double g4=0.;
239      if(A<2.5)g4=std::exp(4*(rho4-z))/(1.+std::exp(12*(tau4-z)));
240      G4double g8=0.;
241      if(A<2.5)g8=std::exp(8*(rho8-z))/(1.+std::exp(24*(tau8-z)));
242      G4double hp=.0375*(z-16.5)+sh*std::exp(-.11*z);
243      G4double fp=hp/(1.+std::exp(4*(7.-z)));
244      G4double sig=(sigd+sigh)/fr+g1+g2+g4+g8+sp*fp;
245      np++;
246      if(np==7)
247      {
248        if(en==nC-1) fileRes<<sig<<"};"<<G4endl;
249        else         fileRes<<sig<<","<<G4endl<<"    ";
250          }
251      else      fileRes<<sig<<",";
252      if(np==7) np=0;
253      //if(en==nC)
254          //{
255      //  G4double fsig=sp*(.0375*(z-16.5)+shd*exp(-.11*z));
256      //  fileRes<<">>> fun="<<fsig<<G4endl;
257          //}
258    } // End of the point LOOP
259  } // End of the isotop LOOP
260  return EXIT_SUCCESS;
261}
Note: See TracBrowser for help on using the repository browser.