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

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

import all except CVS
File size: 12.1 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
28// .....................................................
29// Created: M.V. Kossov, CERN/ITEP(Moscow), 30-Sept-05
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 r=0.;
65    double max=1.;
66    double H=1.22;
67    if(A==1)                   // Proton
68                  {
69      r=std::sqrt(Q2/1.66);
70      max=.5;
71    }
72    else if(A<13)              // Light nuclei
73                  {
74      double f=Q2/4.62;
75      r=f*f;
76      max=.3;
77      if(A>2) H=1.;
78    }
79    else if(A>0)               // Heavy nuclei
80                  {
81      double f=Q2/3.4;
82      double ff=f*f;
83      r=ff*ff;
84      max=.5;
85      H=1.;
86    }
87    else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
88    //
89    N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
90    Del=(1.+r)/(12.5+r/max);
91    double S=std::pow(1.+.6/Q2,-1.-Del);
92    D=H*S*(1.-.5*S);
93    V=3*(1.-D)*(N-1.);
94    double cc=Q2/.08;
95    double cc2=cc*cc;
96    //double C=(1.+cc2)/(1.+cc2/.24);              // Weak?
97    double C=(1.+cc2)/(1.+cc2/.24)/(1.+Q2/21.6); // EM
98    double c3=C+C+C;
99    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
100    U2=(c3+N-3.)*uu;
101    U3=c3*uu;
102    //UU=uu+uu+uu; // @@
103    mA  = A;
104    mQ2 = Q2;
105    mN  = N;
106    mD  = D;
107    mDel=Del;
108    mU2 =U2;
109    mU3 =U3;
110    //mUU =UU; // @@
111    mV  =V;
112  } 
113  // From here the Q2 coefficients are used
114  double x1=std::pow(1.-x,N-2.);
115  double pp=D*std::pow(x,-Del)*x1;
116  double dir=V*x*x1;
117  double per=U2*pp;
118                f2 = per + dir;
119  //double W2=MN2-MN2*x+Q2/x-Q2;
120  //if(W2<W2t)
121                //{
122  //  per=UU*pp;
123  //  xf3= per+dir;
124  //}
125  //else
126  xf3= U3*pp+dir;
127  fL = per/4.;
128  return;
129}
130
131void getFun(int A, double lx, double Q2, double* f)
132{
133  double f2=0., xf3=0., fL=0.;
134  if (lx>0.5) G4cerr<<"***getFun: ln(x)="<<lx<<">.5"<<G4endl;
135  double x=std::exp(lx);
136  double x2=x*x;
137  strucf(A, x, Q2, f2, xf3, fL);
138  f[0]=f2;             // direct part
139  f[1]=(-f2+xf3)/x;    // *y (neutrino) part
140  f[2]=(-f2-xf3)/x;    // *y (anti-neutrino) part
141  f[3]=(f2-fL-xf3)/x2; // *y2 (neutrino) part
142  f[4]=(f2-fL+xf3)/x2; // *y2 (anti-neutrino) part
143}
144
145int main()
146{
147  const double reps=.001;        // relative accuracy of the total Q2 integral calculation
148  const double xeps=.0001;      // relative accuracy of the total X integral calculation
149  //           =========
150  const double GF=1.16637e-5;   // Fermi constant in GeV^-2
151  const double GF2=GF*GF;       // Squared Fermi constant in GeV^-4
152  const double MW=80.425;       // Mass of W-boson in GeV
153  const double MW2=MW*MW;       // Squared mass of W-boson in GeV^2
154  const double MW4=MW2*MW2;     // Quadro mass of W-boson in GeV^4
155  const double hc2=38937932300.;// (hc)^2 in GeV^2*10^-38cm2 to convert GeV^-2 to 10^-38cm2
156  const double pif=3.14159265*4;// 4pi
157  const double sik=GF2*hc2/pif; // precalculated coefficient
158  //const double mpi=.1349766;    // pi0 meson mass in GeV
159  const double mpi=.13957018;   // charged pi meson mass in GeV
160  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
161  //const double me=.00051099892; // electron mass in GeV
162  //const double me2=me*me;       // m_e^2 in GeV^2
163  //const double hme2=me2/2;      // .5*m_e^2 in GeV^2
164  const double mmu=.105658369;  // mu meson mass in GeV
165  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
166  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
167  //const double mtau=1.777;      // tau meson mass in GeV
168  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
169  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
170  //const double mp=.93827203;    // proton mass in GeV
171  //const double mn=.93956536;    // neutron mass in GeV
172  //const double md=1.87561282;   // deuteron mass in GeV
173  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
174  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
175  //const double MD=1.232;        // proton mass in GeV
176  //const double mp2=mp*mp;       // m_p^2 in GeV^2
177  const double MN2=MN*MN;       // M_N^2 in GeV^2
178  const double dMN=MN+MN;       // 2*M_N in GeV
179  const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
180  const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
181  //const double EminE=me+me2/dMN;// Threshold for muon production
182  const double EminMu=mmu+mmu2/dMN;  // Threshold for muon production
183  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
184  //
185  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
186  const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
187  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
188  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
189  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
190  // _____ Begin of Test Area
191  //Genfun::LogGamma logGamma;
192  //double n=4.9;
193                //double g=exp(logGamma(n));
194  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
195  // ^^^^^ End of Test Area
196  //
197  double f[5];                    // A working array
198  int    A=12;                    // Neucleus for which calculations should be done
199  double lEnuMin=std::log(EminMu);     // Log of Minimum energy of neutrino
200  double lEnuMax=std::log(300.);       // Log of Maximum energy of neutrino
201  int    nE=20;
202  double dlE=(lEnuMax-lEnuMin)/nE;
203  lEnuMin+=dlE/10;
204  lEnuMax+=dlE/5;
205  G4cout<<"Emin="<<EminMu<<",lEi="<<lEnuMin<<",lEa="<<lEnuMax<<",dlE="<<dlE<<G4endl;
206  for(double lEnu=lEnuMin; lEnu<lEnuMax; lEnu+=dlE)
207                {
208    double Enu=std::exp(lEnu);         // Energy of neutrino/anti-neutrino
209    double dEnu=Enu+Enu;          // doubled energy of nu/anu
210    double Enu2=Enu*Enu;          // squared energy of nu/anu
211    double Emu=Enu-mmu;           // Free Energy of neutrino/anti-neutrino
212    double Emu2=Emu*Emu;          // squared energy of nu/anu
213    double ME=Enu*MN;             // M*E
214    double dME=ME+ME;             // 2*M*E
215    double DIStsig=1.;            // Total curent DIS cross-section to be integrated
216    double DISmsig=1.e20;         // Total remembered DIS cross-section
217    double dEMN=(dEnu+MN)*ME;
218    double MEm=ME-hmmu2;
219    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
220    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
221    double ymax=(E2M+sqE)/dEMN;
222    double ymin=(E2M-sqE)/dEMN;
223    double rmin=1.-ymin;
224    double rhm2E=hmmu2/Enu2;
225    double Q2min=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E));
226    double Q2max=dME*ymax;
227    // *** Additional notQE limit ****
228    //double Q2nqe=Emu*dMN-mcV;
229    //if(Q2max>Q2nqe) Q2max=Q2nqe; // !! only for not QE
230    // *** End of notQE ^^^^^
231    int    nQ2=8;
232    //G4cout<<"*** E="<<Enu<<", Q2i="<<Q2min<<" < Q2a="<<Q2max<<", yi="<<ymin<<" < ya="
233    //      <<ymax<<G4endl;
234    //while(fabs(DIStsig-DISmsig)/DIStsig>reps)
235    //{
236      DISmsig=DIStsig;
237      DIStsig=0.;
238      nQ2*=2;
239      double dQ2=(Q2max-Q2min)/nQ2;
240      if(dQ2>0) for(double Q2=Q2min+dQ2/2; Q2<Q2max; Q2+=dQ2)
241                                  {
242        double DISxint=1.;        // Curent DIS x-integral
243        double DISmint=1.e20;     // Remembered DIS x-integral
244        double Q2M=Q2+MW2;
245        double dik=MW4/Q2M/Q2M;
246        double qmc=Q2+mcV;
247        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // Quasielastic boundary
248        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
249        //double xN=Q2/dME;
250        double xN=Q2/MN/(Emu+std::sqrt(Emu2+Q2));
251        //double lXmin=log(xN/ymax);
252        double lXmin=std::log(xN);
253        double lXmax=0.;               // QES is in DIS
254        // ****** QE ********
255        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
256        //if(lXQES<lXmax) lXmax=lXQES;   // A cut which excludes QES (see above "notQE")
257        // *** End of QE^^^^^
258        //double lXmax=lXQES;          // Cut off quasielastic
259        int    nX=8;
260        while(std::fabs(DISxint-DISmint)/DISxint>xeps)
261        {
262          DISmint=DISxint;
263          DISxint=0.;
264          nX*=2;
265          double dlX=(lXmax-lXmin)/nX;
266          for(double lX=lXmin+dlX/2; lX<lXmax; lX+=dlX)
267                                      {
268            getFun(A, lX, Q2, f);
269            DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
270            //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
271            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
272          }
273          DISxint*=dlX;
274          //G4cout<<"--- E="<<Enu<<" --- Q2="<<Q2<<" --- nX="<<nX<<", iX="<<DISxint
275          //      <<", mX="<<DISmint<<", rX="<<(DISxint-DISmint)/DISxint<<G4endl;
276        }
277        //G4cout<<"(E="<<Enu<<"), r="<<Q2/Q2mx<<" "<<DISxint*dik*sik<<G4endl;//NoNormedByEn
278        // QE
279        G4cout<<"(E="<<Enu<<"), Q2="<<Q2<<" "<<DISxint*dik*sik<<G4endl;//Not Normed By En
280        // notQE
281        //G4cout<<"(E="<<Enu<<"), r="<<Q2/Q2max<<" "<<DISxint*dik*sik/Enu<<G4endl;//NormByE
282        DIStsig+=DISxint*dik;
283      }
284      DIStsig*=dQ2;
285      //G4cout<<"=== E="<<Enu<<" ===> nQ="<<nQ2<<", iQ="<<DIStsig<<", mQ="<<DISmsig
286                                                //      <<", rQ="<<(DIStsig-DISmsig)/DIStsig<<G4endl;
287                                //}
288    DIStsig*=sik/Enu;
289    G4cout<<"*** E="<<Enu<<",sig/E="<<DIStsig<<G4endl;
290                } // End of the Enery LOOP
291  // int np=0;
292  //for(int m=0; m<2; m++)
293  //{
294  //  //fileNuMuX<<"  static const G4double SH"<<n<<"[nH]={"<<G4endl<<"    ";
295  //  //G4cout<<"**** A_high="<<m<<G4endl;
296  //  np=0;
297  //  int nC=14;
298         //  for(G4int en=0; en<nC; en++)
299         //  {
300  //    //G4double sig=1.;
301  //    np++;
302  //    //if(np==7)           // Write by 7 number in brackets
303  //    //{
304  //    //  if(en==nC-1) fileNuMuX<<sig<<"};"<<G4endl;
305  //    //  else         fileNuMuX<<sig<<","<<G4endl<<"    ";
306         //    //}
307  //    //else           fileNuMuX<<sig<<",";
308  //    //if(np==7) np=0;
309  //  } // End of the point LOOP
310  //} // End of the isotop LOOP
311  return 0;
312}
Note: See TracBrowser for help on using the repository browser.