source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration_N.cc @ 1071

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

fichier ajoutes
File size: 11.5 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: NuMu DIS (x,Q2) approximation is integrated over x & Q2
28// .....................................................[ for a <nucleon> ]
29// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Oct-07
30//
31//=====================================================================
32#include "globals.hh"
33#include <iostream>
34#include <fstream>
35#include <vector>
36#include "G4ios.hh"
37#include <CLHEP/GenericFunctions/LogGamma.hh>
38
39void strucf(int A, double x, double Q2, double& f2, double& xf3, double& fL)
40{
41  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
42  //const double MN2=MN*MN;       // M_N^2 in GeV^2
43  //const double mpi=.13957018;   // charged pi meson mass in GeV
44  //const double Wt=MN+mpi;       // Delta threshold
45  //const double W2t=Wt*Wt;       // Squared Delta threshold
46  const Genfun::LogGamma lGam;
47  static int mA=0;
48  static double mQ2=0., mN, mD, mDel, mU2, mU3, mV;
49  //static double mUU;
50  double N=3., D=0., Del=0., U2=0., U3=0., V=0.;
51  //double UU=0.;
52  if(A==mA && Q2==mQ2)         // Associative memory for acceleration
53                {
54    N  =mN;
55    D  =mD;
56    Del=mDel;
57    U2 =mU2;
58    U3 =mU3;
59    //UU =mUU;
60    V  =mV;
61  }
62  else
63  {
64    double Q=std::sqrt(Q2);
65    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
66    Del=.077+.393/(1.+5./Q);
67    D=.68*std::pow(1.+.145/Q2,-1.-Del);
68    V=3*(N-1.);
69    double c3=.75;
70    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
71    U2=(c3+N-3.)*uu;
72    U3=c3*uu;
73    //UU=uu+uu+uu; // @@
74    mA  = A;
75    mQ2 = Q2;
76    mN  = N;
77    mD  = D;
78    mDel=Del;
79    mU2 =U2;
80    mU3 =U3;
81    //mUU =UU; // @@
82    mV  =V;
83  } 
84  // From here the Q2 coefficients are used
85  double x1=std::pow(1.-x,N-2.);
86  double pp=D*std::pow(x,-Del)*x1;
87  double dir=(1-D)*V*x*x1;
88  double per=U2*pp;
89                f2 = per + dir;
90  //double W2=MN2-MN2*x+Q2/x-Q2;
91  //if(W2<W2t)
92                //{
93  //  per=UU*pp;
94  //  xf3= per+dir;
95  //}
96  //else
97  xf3= U3*pp+dir;
98  fL = per/5.;
99  return;
100}
101
102void getFun(int A, double lx, double Q2, double* f)
103{
104  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
105  const double MN2=MN*MN;   // Squared Nucleon mass (inside nucleus, atomic mass unit, GeV)
106  double f2=0., xf3=0., fL=0.;
107  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
108  double x=std::exp(lx);
109  double x2=x*x;
110  double c=0.;
111  if(Q2>.0001) c=x2*MN2/Q2;
112  double c2=1+c+c;     // Q^2/nu^2 correction
113  strucf(A, x, Q2, f2, xf3, fL);
114  f[0]=f2;             // direct part
115  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
116  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
117  f[3]=(f2*c2-fL-xf3)/x2; // *y2 (neutrino) part
118  f[4]=(f2*c2-fL+xf3)/x2; // *y2 (anti-neutrino) part
119}
120
121int main()
122{
123  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
124  const double xeps=.0001;      // relative accuracy of the total X integral calculation
125  //           =========
126  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
127  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
128  const double MW=80.425;       // Mass of W-boson in GeV
129  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
130  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
131  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
132  const double pif=3.14159265*4;// 4pi
133  const double sik=GF2*hc2/pif; // precalculated coefficient
134  //const double mpi=.1349766;    // pi0 meson mass in GeV
135  //const double mpi=.13957018;   // charged pi meson mass in GeV
136  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
137  //const double me=.00051099892; // electron mass in GeV
138  //const double me2=me*me;       // m_e^2 in GeV^2
139  //const double hme2=me2/2;      // .5*m_e^2 in GeV^2
140  const double mmu=.105658369;  // mu meson mass in GeV
141  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
142  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
143  //const double mtau=1.777;      // tau meson mass in GeV
144  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
145  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
146  //const double mp=.93827203;    // proton mass in GeV
147  //const double mn=.93956536;    // neutron mass in GeV
148  //const double md=1.87561282;   // deuteron mass in GeV
149  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
150  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
151  //const double MD=1.232;        // proton mass in GeV
152  //const double mp2=mp*mp;       // m_p^2 in GeV^2
153  const double MN2=MN*MN;       // M_N^2 in GeV^2
154  const double dMN=MN+MN;       // 2*M_N in GeV
155  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
156  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
157  //const double EminE=me+me2/dMN;// Threshold for muon production
158  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
159  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
160  //
161  const double mc=.261;           // parameter of W2>(M_N^2+M_D^2)/2 cut for QuasiEl/Delta
162  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
163  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
164  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
165  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
166  // _____ Begin of Test Area
167  //Genfun::LogGamma logGamma;
168  //double n=4.9;
169                //double g=exp(logGamma(n));
170  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
171  // ^^^^^ End of Test Area
172  //
173  double f[5];                    // A working array
174  int    A=12;                    // Neucleus for which calculations should be done
175  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
176  double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of MaximumEnergy of neutrino
177  int    nE=63;
178  double dlE=(lEnuMax-lEnuMin)/nE;
179  lEnuMin+=dlE/10;
180  lEnuMax+=dlE/5;
181  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
182  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
183                {
184    double Enu=std::exp(std::exp(lEnu)-1.)*EminMu; // Energy of neutrino/anti-neutrino
185    double dEnu=Enu+Enu;          // doubled energy of nu/anu
186    double Enu2=Enu*Enu;          // squared energy of nu/anu
187    double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
188    double Emu2=Emu*Emu;          // squared energy of nu/anu
189    double ME=Enu*MN;             // M*E
190    double dME=ME+ME;             // 2*M*E
191    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
192    double DISmsig=1.e20;         // Total remembered DIS cross-section
193    double dEMN=(dEnu+MN)*ME;
194    double MEm=ME-hmmu2;
195    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
196    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
197    double ymax=(E2M+sqE)/dEMN;
198    double ymin=(E2M-sqE)/dEMN;
199    double rmin=1.-ymin;
200    double rhm2E=hmmu2/Enu2;
201    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
202    double Q2max=dME*ymax;
203    int    nQ2=8;
204    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
205    //      <<ymax<<G4endl;
206    while(std::fabs(DIStsig-DISmsig)/DIStsig>reps)
207    {
208      DISmsig=DIStsig;
209      DIStsig=0.;
210      nQ2*=2;
211      double dQ2=(Q2max-Q2min)/nQ2;
212      for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
213                                  {
214        double DISxint=1.;        // Curent DIS x-integral
215        double DISmint=1.e20;     // Remembered DIS x-integral
216        double Q2M=Q2+MW2;
217        double dik=MW4/Q2M/Q2M;
218        double qmc=Q2+mcV;
219        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // QuasielastBoundary
220        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
221        //double xN=Q2/dME;
222        double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
223        //double lXmin=log(xN/ymax);
224        double lXmin=std::log(xN);
225        // ****** QE ********
226        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
227        // *** End of QE^^^^^
228        double lXmax=0.;               // QES is in DIS
229        //double lXmax=lXQES;          // Cut off quasielastic
230        int    nX=8;
231        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
232        {
233          DISmint=DISxint;
234          DISxint=0.;
235          nX*=2;
236          double dlX=(lXmax-lXmin)/nX;
237          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
238                                      {
239            getFun(A, lX, Q2, f);
240            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
241            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
242            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
243          }
244          DISxint*=dlX;
245          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
246          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
247        }
248        //G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<", I="<<DISxint/dik/dik<<G4endl;
249        DIStsig+=DISxint*dik;
250      }
251      DIStsig*=dQ2;
252      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
253                                                //      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
254    }
255    //===== tot/qe choice ====
256    //DIStsig*=sik/Enu;
257    //G4cout<<"***total*** E="<<Enu<<",sig/E= "<<DIStsig<<G4endl;
258    //...................
259    DIStsig*=sik;
260    G4cout<<"***qelas*** E="<<Enu<<",sig= "<<DIStsig<<G4endl;
261    //===== End of the choice
262                } // End of the Enery LOOP
263  // int np=0;
264  //for(int m=0; m<2; m++)
265  //{
266  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
267  //  //G4cout<<"**** A_high="<<m<<G4endl;
268  //  np=0;
269  //  int nC=14;
270         //  for(G4int en=0; en<nC; en++)
271         //  {
272  //    //G4double sig=1.;
273  //    np++;
274  //    //if(np==7)           // Write by 7 number in brackets
275  //    //{
276  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
277  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
278         //    //}
279  //    //else           fileNuMuX<<sig<<",";
280  //    //if(np==7) np=0;
281  //  } // End of the point LOOP
282  //} // End of the isotop LOOP
283  return EXIT_SUCCESS;
284}
Note: See TracBrowser for help on using the repository browser.