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

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

import all except CVS
File size: 10.6 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(Q2) fixed step integration
28// .....................................................
29// Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-2005
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
39// All calculations have been done for C12 nucleus
40
41double nuQ2(double Q2) // (Q2 is in GeV^2)
42{
43  //  x      *     *     *     *     *     *     *     *      *     *
44  // -1-    -2-   -3-   -4-   -5-   -6-   -7-   -8-   -9-   -10-  -11-
45  // 2.068 2.634 4.333 3.000 10.59 .0011 18.74 1.484 134950 .0755 4.5
46  G4double y=Q2/3;
47  G4double z=std::pow((Q2/2.634),4.333);
48                // The first (1.+Q2)**4.5 fuctor is just normalization, which payed back by 1/(1+Q2)^3.5
49  G4double f=std::pow((1.+Q2),4.5)*
50                                (1.+z*(1.+18.74*std::exp(-Q2/1.484)-134950*std::exp(-Q2/.0755)))/std::pow((1.+y+.0011*y*y),10.59);
51  //G4cout<<"Q2="<<Q2<<", y="<<y<<", z="<<z<<", f="<<f<<G4endl;
52  return f;
53}
54
55double anuQ2(double Q2)
56{
57  //    x     *     *     *     *     *     *     *     -     *
58  //   -1-   -2-   -3-   -4-   -5-   -6-   -7-   -8-   -9-   -10-
59  // .2000 .0100 .3000 .4000 3.900 .0150 .2000 6.500 0.000 .00000001
60  G4double  y=Q2/.4;
61  G4double  y2=y*y;
62  G4double  z=std::pow((Q2/.01),-.3);
63                // The first (1.+Q2)**6.5 fuctor is just normalization, which payed back by 1/(1+Q2)^5.5
64  G4double f=(1.+z)*std::pow((1+Q2),6.5)*std::pow(Q2,-.2)/std::pow((1.+y+(.015+.00000001*y2)*y2),3.9);
65  //G4cout<<"Q2="<<Q2<<", y="<<y<<", z="<<z<<", f="<<f<<G4endl;
66  return f;
67}
68
69int main()
70{
71  const double eps=.000001;      // relative accuracy of the integral calculation
72  //           =========
73  const double mmu=.105658369;  // mu meson mass in GeV
74  const double mmu2=mmu*mmu;    // m_mu^2 in GeV^2
75  const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
76  //const double mtau=1.777;      // tau meson mass in GeV
77  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
78  //const double hmtau2=mtau2/2;  // .5*m_e^2 in GeV^2
79  //const double mp=.93827203;    // proton mass in GeV
80  //const double mn=.93956536;    // neutron mass in GeV
81  //const double md=1.87561282;   // deuteron mass in GeV
82  const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
83  //const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
84  //const double MD=1.232;        // proton mass in GeV
85  //const double mp2=mp*mp;       // m_p^2 in GeV^2
86  const double MN2=MN*MN;       // M_N^2 in GeV^2
87  const double dMN=MN+MN;       // 2*M_N in GeV
88  //const double dMN2=MN2+MN2;    // 2*M_N^2 in GeV^2
89  //const double fMN2=dMN2+dMN2;  // 4*M_N^2 in GeV^2
90  //const double EminE=me+me2/dMN;// Threshold for muon production
91  const double Emin=mmu+mmu2/dMN;  // Threshold for muon production
92  //const double EminTau=mmu+mmu2/dMN; // Threshold for muon production
93  //
94  //const double mc=.3;           // parameter of W>M+mc cut for Quasi-Elastic/Delta
95  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
96  //const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
97  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
98  //fileNuMuX.setf( std::ios::scientific, std::ios::floatfield );
99  // _____ Begin of Test Area
100  //Genfun::LogGamma logGamma;
101  //double n=4.9;
102                //double g=exp(logGamma(n));
103  //G4cout<<"Gamma("<<n<<") = "<<g<<G4endl;
104  // ^^^^^ End of Test Area
105  //
106  const int niQ2=100;           // Number of points in the Q2 integration
107  const int diQ2=niQ2+1;        // Dimention for arrays
108  double Xl[diQ2];
109  double inl[diQ2];
110  double Enu=Emin;              // Initial E=minE 
111  double Emo=Emin+Emin;         // First new
112  double Emv=0.;
113  double Q2i=1.;
114  double Q2o=1.;
115  double Q2v=0.;
116                int nit=0;
117  while (std::fabs(Q2v-Q2o)>eps && nit<10)
118  {
119    nit++;
120    double Emm=Emo;               // Candidate for the next point is predefined
121    double dEnu=Enu+Enu;          // doubled energy of nu/anu
122    double Enu2=Enu*Enu;          // squared energy of nu/anu
123    double ME=Enu*MN;             // M*E
124    double dEMN=(dEnu+MN)*ME;
125    double MEm=ME-hmmu2;
126    double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
127    double E2M=MN*Enu2-(Enu+MN)*hmmu2;
128    //double ymax=(E2M+sqE)/dEMN;
129    double ymin=(E2M-sqE)/dEMN;
130    double rmin=1.-ymin;
131    double rhm2E=hmmu2/Enu2;
132    Q2i=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // minQ2 for Enu
133    //G4cout<<"MinSearch: E="<<Enu<<": Q2="<<Q2i<<G4endl;
134    if(Emv>0.0001) // Not initialization (not first two) Use three points
135                                {
136      //double d=Enu*Emo*(Enu-Emo)+Enu*Emv*(Emv-Enu)+Emo*Emv*(Emo-Emv);
137      double a=Q2i*Emo-Enu*Q2o+Q2v*Enu-Emv*Q2i+Q2o*Emv-Emo*Q2v;
138                                                double b=Enu*Enu*(Q2o-Q2v)+Emo*Emo*(Q2v-Q2i)+Emv*Emv*(Q2i-Q2o);
139      //double a
140                                  Emm=-b/(a+a);
141      if(Q2v<Q2i)                 // swap to make Q2v the biggest
142                                  {
143        double q2=Q2v;
144        double en=Emv;
145        Q2v=Q2i;
146        Emv=Enu;
147        Q2i=q2;
148        Enu=en;
149      }
150      if(Q2v<Q2o)                 // swap to make Q2v the biggest
151                                  {
152        double q2=Q2v;
153        double en=Emv;
154        Q2v=Q2o;
155        Emv=Emo;
156        Q2o=q2;
157        Emo=en;
158      }
159    }
160    else if(Q2o>0.0000000001) Emm=Enu+Emin; // the second step
161    //G4cout<<"***Enu="<<Enu<<", Emm="<<Emm<<", Emin="<<Emin<<G4endl;
162    Emv=Emo;
163    Q2v=Q2o;
164    Emo=Enu;
165    Q2o=Q2i;
166    Enu=Emm;
167    //G4cout<<"___Q2i="<<Q2i<<": Q2o="<<Q2o<<", Q2v="<<Q2v<<G4endl;
168  }
169  double Q2min=Q2i-eps;
170  double Q2max=600.;            // covers the calculated region
171  // ----------------- nu/anu switch ------------------------------
172  //bool nu=true;
173  bool nu=false;
174  // -----------------------------**************************-------
175  // *************** Convert to the magic variable ****************
176  double Xmin=0.;
177  double Xmax=0.;
178                if(nu)
179  {
180    Xmin=std::pow(1+Q2max,-3.5);
181    Xmax=std::pow(1+Q2min,-3.5);
182                } 
183  else
184  {
185    Xmin=std::pow(1+Q2max,-5.5);
186    Xmax=std::pow(1+Q2min,-5.5);
187                }
188  G4cout<<"Ei="<<Enu<<":Q2i="<<Q2min<<",Q2a="<<Q2max<<",Xmi="<<Xmin<<",Xma="<<Xmax<<G4endl;
189  int    nQ2=8;                 // nitial #Of points for the overall integration
190  double DISmsig=0.;
191  double DIStsig=1.;
192                double Q2=0.;                 // Prototype for the integration
193  while(std::fabs(DIStsig-DISmsig)/DIStsig>eps)
194  {
195    DISmsig=DIStsig;
196    DIStsig=0.;
197    nQ2*=2;
198    double dX=(Xmax-Xmin)/nQ2;
199    double hX=dX/2;
200    for(double X=Xmin+hX; X<Xmax; X+=dX)
201                                {
202      if(nu) // neutrino
203                                                {
204        Q2=std::pow(X,-.2857142); // -1./3.5
205        DIStsig+=nuQ2(Q2);
206      }
207      else  // anti-neutrino
208                                                {
209        Q2=std::pow(X,-.1818182); // -1./5.5
210        DIStsig+=anuQ2(Q2);
211      }
212    }
213    DIStsig*=dX;
214    //G4cout<<"n="<<nQ2<<",i="<<DIStsig<<",m="<<DISmsig<<",r="<<(DIStsig-DISmsig)/DIStsig
215                                //                  <<G4endl;
216                }
217  G4cout<<"Total: nQ2="<<nQ2<<", Int="<<DIStsig<<G4endl;
218  // ***************** Calculate the reversed table *****************
219  double dInt=DIStsig/niQ2;             // Step for the integral
220  DIStsig=0.;
221  DISmsig=0.;
222  nQ2*=2;
223  double dX=(Xmax-Xmin)/nQ2;
224  dInt/=dX;                     // To avoid multiplication (nQ2 is fixed)
225  double sum=dInt;
226  int nn=1;
227  Xl[0]=0.;
228  double hX=dX/2;
229  for(double X=Xmin+hX; X<Xmax; X+=dX)
230                {
231    if(nu) // neutrino
232                                {
233      Q2=std::pow(X,-.2857142); // -1./3.5
234      DIStsig+=nuQ2(Q2);
235    }
236    else  // anti-neutrino
237                                {
238      Q2=std::pow(X,-.1818182); // -1./5.5
239      DIStsig+=anuQ2(Q2);
240    }
241    if(DIStsig>sum+eps)
242                                {
243      inl[nn]=sum*dX;
244      Xl[nn]=X-(DIStsig-sum)*dX/(DIStsig-DISmsig);
245      G4cout<<"sum="<<inl[nn]<<", Xl["<<nn<<"]="<<Xl[nn]<<G4endl;
246      nn++;
247      sum+=dInt;
248    }
249    DISmsig=DIStsig;
250  }
251  inl[nn]=sum*dX;
252  Xl[nn]=Xmax;
253  G4cout<<"sum="<<inl[nn]<<", Xl["<<nn<<"]="<<Xmax<<G4endl;
254  // ************* The following is just a test (better than .5%) ************
255  for(int i=1; i<=niQ2; i++)
256  {
257    DIStsig=0.;
258    double Xm=Xl[i];
259    double dX=(Xm-Xmin)/nQ2;
260    double hX=dX/2;
261    for(double X=Xmin+hX; X<Xm; X+=dX)
262                                {
263      if(nu) // neutrino
264                                                {
265        Q2=std::pow(X,-.2857142); // -1./3.5
266        DIStsig+=nuQ2(Q2);
267      }
268      else  // anti-neutrino
269                                                {
270        Q2=std::pow(X,-.1818182); // -1./5.5
271        DIStsig+=anuQ2(Q2);
272      }
273    }
274    DIStsig*=dX;
275    G4cout<<"i="<<i<<", v="<<DIStsig<<", d="<<std::fabs(DIStsig-inl[i])/DIStsig<<G4endl;
276  }
277  // ***************** Calculate the direct table *****************
278  dX=(Xmax-Xmin)/niQ2;
279  double nor=inl[niQ2]/niQ2;
280  for(int i=1; i<=niQ2; i++)
281  {
282    DIStsig=0.;
283    double Xm=Xmin+dX*i;
284    double rX=(Xm-Xmin)/nQ2;
285    double hX=rX/2;
286    for(double X=Xmin+hX; X<Xm; X+=rX)
287                                {
288      if(nu) // neutrino
289                                                {
290        Q2=std::pow(X,-.2857142); // -1./3.5
291        DIStsig+=nuQ2(Q2);
292      }
293      else  // anti-neutrino
294                                                {
295        Q2=std::pow(X,-.1818182); // -1./5.5
296        DIStsig+=anuQ2(Q2);
297      }
298    }
299    DIStsig*=rX/nor;
300    G4cout<<"i="<<i<<", Xm="<<Xm<<", I="<<DIStsig<<G4endl;
301  }
302  G4cout<<"End"<<G4endl;
303  DIStsig=0.;
304  return 0;
305}
Note: See TracBrowser for help on using the repository browser.