Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration.cc

    r819 r962  
    6262  else
    6363  {
    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;
     64                  //double r=0.;
     65    double min=.077;             // Delta(0)
     66                                double max=.47;              // Hard Pomeron
     67    //double H=1.22;
     68    if(Q2<=0.) G4cout<<"strucf: Q2="<<Q2<<" <= 0"<< G4endl;
     69    double Q=std::sqrt(Q2);
     70    //if(A==1)                   // Proton
     71                  //{
     72    //  r=std::sqrt(Q2/1.66);
     73    //  max=.5;
     74    //}
     75    //else if(A<13)              // Light nuclei
     76                  //{
     77    //  double f=Q2/4.62;
     78    //  r=f*f;
     79    //  max=.3;
     80    //  if(A>2) H=1.;
     81    //}
     82    //else if(A>0)               // Heavy nuclei
     83                  //{
     84    //  double f=Q2/3.4;
     85    //  double ff=f*f;
     86    //  r=ff*ff;
     87    //  max=.5;
     88    //  H=1.;
     89    //}
     90    //else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl;
    8891    //
    8992    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    Del=min+(max-min)/(1.+5./Q);
     94    D=0.68*std::pow(1.+.145/Q2,-1.-Del); // 0.68=1-0.34, m2=.145 GeV2
    9395    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;
     96    double c3=.75;                       // 3*0.25
    9997    double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N;
    10098    U2=(c3+N-3.)*uu;
    10199    U3=c3*uu;
    102     //UU=uu+uu+uu; // @@
    103100    mA  = A;
    104101    mQ2 = Q2;
    … …  
    125122  //else
    126123  xf3= U3*pp+dir;
    127   fL = per/4.;
     124  fL = per/5.;                             // 20%
    128125  return;
    129126}
    … …  
    157154  const double sik=GF2*hc2/pif; // precalculated coefficient
    158155  //const double mpi=.1349766;    // pi0 meson mass in GeV
    159   const double mpi=.13957018;   // charged pi meson mass in GeV
     156  //const double mpi=.13957018;   // charged pi meson mass in GeV
    160157  //const double mpi2=mpi*mpi;    // m_pi^2 in GeV^2
    161158  //const double me=.00051099892; // electron mass in GeV
    … …  
    168165  //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2
    169166  //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
     167  const double mp=.93827203;    // proton mass in GeV
     168  const double mn=.93956536;    // neutron mass in GeV
    172169  //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
     170  //const double MN=.931494043;   // Nucleon mass (inside nucleus, atomic mass unit, GeV)
     171  const double MN=(mn+mp)/2;    // Nucleon mass (mean free) in GeV
     172  const double MD=1.232;        // proton mass in GeV
    176173  //const double mp2=mp*mp;       // m_p^2 in GeV^2
    177174  const double MN2=MN*MN;       // M_N^2 in GeV^2
    … …  
    184181  //
    185182  //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
     183  const double mc=std::sqrt((MN*MN+MD*MD)/2)-MN; // Squared mean between N and \Delta
     184  //const double mc=mpi;          // parameter of W>M+mc cut for Quasi-Elastic/Delta
    187185  const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic
    188186  //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out);
    … …  
    198196  int    A=12;                    // Neucleus for which calculations should be done
    199197  double lEnuMin=0;               // LogLog of Minimum energy of neutrino
    200   double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of Maximum energy of neutrino
    201   int    nE=31;
     198  double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of MaximumEnergy of neutrino
     199  int    nE=49;                   // Number of points
    202200  double dlE=(lEnuMax-lEnuMin)/nE;
    203201  lEnuMin+=dlE/10;
    … …  
    241239        double dik=MW4/Q2M/Q2M;
    242240        double qmc=Q2+mcV;
    243         double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // Quasielastic boundary
     241        double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2);//QuasielasticBoundary
    244242        //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c)
    245243        //double xN=Q2/dME;
    … …  
    248246        double lXmin=std::log(xN);
    249247        // ****** QE ********
    250         //if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES
     248        if(lXQES>lXmin) lXmin=lXQES;   // A cut which leaves only QES >>>>>>>>>>>>>>>>
    251249        // *** End of QE^^^^^
    252250        double lXmax=0.;               // QES is in DIS
    … …  
    262260                                      {
    263261            getFun(A, lX, Q2, f);
    264             DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
    265             //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
     262            // ***** Neutrino/Antineutrino switch ******>>>>>>>>>>>>>>>>>>>>>>>>>>>>
     263            //DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino
     264            DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino
    266265            //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl;
    267266          }
    … …  
    279278    //===== tot/qe choice ====
    280279    DIStsig*=sik/Enu;
    281     G4cout<<"***total*** E="<<Enu<<",sig/E="<<DIStsig<<G4endl;
     280    //G4cout<<"***total-neutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
     281    G4cout<<"***total-antineutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
     282    //G4cout<<"***quasiel-nu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
     283    //G4cout<<"***quasiel-antinu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl;
    282284    //...................
    283285    //DIStsig*=sik;
Note: See TracChangeset for help on using the changeset viewer.