G4QNuNuNuclearCrossSection Class Reference

#include <G4QNuNuNuclearCrossSection.hh>

Inheritance diagram for G4QNuNuNuclearCrossSection:

G4VQCrossSection

Public Member Functions

 ~G4QNuNuNuclearCrossSection ()
G4double ThresholdEnergy (G4int Z, G4int N, G4int PDG=14)
virtual G4double GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0)
G4double CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
G4int GetExchangePDGCode ()
G4double GetDirectPart (G4double Q2)
G4double GetNPartons (G4double Q2)
G4double GetQEL_ExchangeQ2 ()
G4double GetNQE_ExchangeQ2 ()
G4double GetLastTOTCS ()
G4double GetLastQELCS ()

Static Public Member Functions

static G4VQCrossSectionGetPointer ()

Protected Member Functions

 G4QNuNuNuclearCrossSection ()

Detailed Description

Definition at line 49 of file G4QNuNuNuclearCrossSection.hh.


Constructor & Destructor Documentation

G4QNuNuNuclearCrossSection::G4QNuNuNuclearCrossSection (  )  [inline, protected]

Definition at line 53 of file G4QNuNuNuclearCrossSection.hh.

00053 {};

G4QNuNuNuclearCrossSection::~G4QNuNuNuclearCrossSection (  ) 

Definition at line 78 of file G4QNuNuNuclearCrossSection.cc.

00079 {
00080   G4int lens=TX->size();
00081   for(G4int i=0; i<lens; ++i) delete[] (*TX)[i];
00082   delete TX;
00083   G4int hens=QE->size();
00084   for(G4int i=0; i<hens; ++i) delete[] (*QE)[i];
00085   delete QE;
00086 }


Member Function Documentation

G4double G4QNuNuNuclearCrossSection::CalculateCrossSection ( G4bool  CS,
G4int  F,
G4int  I,
G4int  PDG,
G4int  Z,
G4int  N,
G4double  Momentum 
) [virtual]

Implements G4VQCrossSection.

Definition at line 262 of file G4QNuNuNuclearCrossSection.cc.

References G4cerr, G4cout, and G4endl.

Referenced by GetCrossSection().

00264 {
00265   static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
00266   static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
00267   static const G4int mL=nE-1;
00268   static const G4double EMi=0.;      // Universal threshold of the reaction in GeV
00269   static const G4double EMa=300.;    // Maximum tabulated Energy of nu_mu in GeV 
00270   // *** Begin of the Associative memory for acceleration of the cross section calculations
00271   static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
00272   static G4bool first=true;          // Flag of initialization of the energy axis
00273   // *** End of Static Definitions (Associative Memory) ***
00274   //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
00275   //G4double TotEnergy2=Momentum;
00276   onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
00277   lastE=Momentum/GeV;                // Kinetic energy of the muon neutrino (in GeV!)
00278   if (lastE<=EMi)                    // Energy is below the minimum energy in the table
00279   {
00280     lastE=0.;
00281     lastSig=0.;
00282     return 0.;
00283   }
00284   G4int Z=targZ;                     // New Z, which can change the sign
00285   if(F<=0)                           // This isotope was not the last used isotop
00286   {
00287     if(F<0)                          // This isotope was found in DAMDB =-------=> RETRIEVE
00288     {
00289       lastTX =(*TX)[I];              // Pointer to the prepared TX function (same isotope)
00290       lastQE =(*QE)[I];              // Pointer to the prepared QE function (same isotope)
00291    }
00292    else                              // This isotope wasn't calculated previously => CREATE
00293    {
00294       if(first)
00295       {
00296         lastEN = new G4double[nE];   // This must be done only once!
00297         Z=-Z;                        // To explain GetFunctions that E-axis must be filled
00298         first=false;                 // To make it only once
00299       }
00300       lastTX = new G4double[nE];     // Allocate memory for the new TX function
00301       lastQE = new G4double[nE];     // Allocate memory for the new QE function
00302       G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
00303       if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
00304       // *** The synchronization check ***
00305       G4int sync=TX->size();
00306       if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
00307       TX->push_back(lastTX);
00308       QE->push_back(lastQE);
00309     } // End of creation of the new set of parameters
00310   } // End of parameters udate
00311   // =-------------= NOW Calculate the Cross Section =---------------------=
00312   if (lastE<=EMi)                   // Check that the neutrinoEnergy is higher than ThreshE
00313   {
00314     lastE=0.;
00315     lastSig=0.;
00316     return 0.;
00317   }
00318   if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
00319   {
00320     G4int chk=1;
00321     G4int ran=mL/2;
00322     G4int sep=ran;  // as a result = an index of the left edge of the interval
00323     while(ran>=2)
00324     {
00325       G4int newran=ran/2;
00326       if(lastE<=lastEN[sep]) sep-=newran;
00327       else                   sep+=newran;
00328       ran=newran;
00329       chk=chk+chk; 
00330     }
00331     if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
00332     G4double lowE=lastEN[sep];
00333     G4double highE=lastEN[sep+1];
00334     G4double lowTX=lastTX[sep];
00335     if(lastE<lowE||sep>=mL||lastE>highE)
00336       G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
00337             <<", sep="<<sep<<", mL="<<mL<<G4endl;
00338     lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
00339     if(!onlyCS)                       // Skip the differential cross-section parameters
00340     {
00341       G4double lowQE=lastQE[sep];
00342       lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
00343 #ifdef pdebug
00344       G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
00345 #endif
00346     }
00347   }
00348   else
00349   {
00350     lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
00351     lastQEL=lastQE[mL];
00352   }
00353   if(lastQEL<0.) lastQEL = 0.;
00354   if(lastSig<0.) lastSig = 0.;
00355   // The cross-sections are expected to be in mb
00356   lastSig*=mb38;
00357   if(!onlyCS) lastQEL*=mb38;
00358   return lastSig;
00359 }

G4double G4QNuNuNuclearCrossSection::GetCrossSection ( G4bool  fCS,
G4double  pMom,
G4int  tgZ,
G4int  tgN,
G4int  pPDG = 0 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 90 of file G4QNuNuNuclearCrossSection.cc.

References CalculateCrossSection(), G4cout, G4endl, ThresholdEnergy(), and G4VQCrossSection::tolerance.

00092 {
00093   static G4int j;                      // A#0f records found in DB for this projectile
00094   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00095   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00096   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00097   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00098   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00099   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00100   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00101   G4double pEn=pMom;
00102 #ifdef debug
00103   G4cout<<"G4QNMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00104         <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00105         <<colN.size()<<G4endl;
00106   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00107 #endif
00108   if(pPDG!=14)
00109   {
00110 #ifdef pdebug
00111     G4cout<<"G4QNMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00112     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00113 #endif
00114     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00115   }
00116   G4bool in=false;                     // By default the isotope must be found in the AMDB
00117   if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00118   {
00119     in = false;                        // By default the isotope haven't be found in AMDB  
00120     lastP   = 0.;                      // New momentum history (nothing to compare with)
00121     lastPDG = pPDG;                    // The last PDG of the projectile
00122     lastN   = tgN;                     // The last N of the calculated nucleus
00123     lastZ   = tgZ;                     // The last Z of the calculated nucleus
00124     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00125     j  = 0;                            // A#0f records found in DB for this projectile
00126     if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
00127     {                                  // The nucleus with projPDG is found in AMDB
00128       if(colN[i]==tgN && colZ[i]==tgZ)
00129       {
00130         lastI=i;
00131         lastTH =colTH[i];                // Last THreshold (A-dependent)
00132 #ifdef pdebug
00133         G4cout<<"G4QNMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00134         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00135 #endif
00136         if(pEn<=lastTH)
00137         {
00138 #ifdef pdebug
00139           G4cout<<"G4QNMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
00140           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00141 #endif
00142           return 0.;                     // Energy is below the Threshold value
00143         }
00144         lastP  =colP [i];                // Last Momentum  (A-dependent)
00145         lastCS =colCS[i];                // Last CrossSect (A-dependent)
00146         if(std::fabs(lastP/pMom-1.)<tolerance)
00147         {
00148 #ifdef pdebug
00149           G4cout<<"G4QNMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00150           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00151 #endif
00152           return lastCS*millibarn;     // Use theLastCS
00153         }
00154         in = true;                       // This is the case when the isotop is found in DB
00155         // Momentum pMom is in IU ! @@ Units
00156 #ifdef pdebug
00157         G4cout<<"G4QNMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00158 #endif
00159         lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
00160 #ifdef pdebug
00161         G4cout<<"G4QNMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00162         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00163 #endif
00164         if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
00165         {
00166 #ifdef pdebug
00167           G4cout<<"G4QNMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00168 #endif
00169           lastTH=pEn;
00170         }
00171         break;                           // Go out of the LOOP
00172       }
00173 #ifdef pdebug
00174       G4cout<<"---G4QNMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00175             <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
00176       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00177 #endif
00178       j++;                             // Increment a#0f records found in DB for this pPDG
00179     }
00180     if(!in)                            // This nucleus has not been calculated previously
00181     {
00182 #ifdef pdebug
00183       G4cout<<"G4QNMNCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00184 #endif
00186       lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
00187       if(lastCS<=0.)
00188       {
00189         lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
00190 #ifdef pdebug
00191         G4cout<<"G4QNMNCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00192 #endif
00193         if(pEn>lastTH)
00194         {
00195 #ifdef pdebug
00196           G4cout<<"G4QNMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00197 #endif
00198           lastTH=pEn;
00199         }
00200       }
00201 #ifdef pdebug
00202       G4cout<<"G4QNMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00203       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00204 #endif
00205       colN.push_back(tgN);
00206       colZ.push_back(tgZ);
00207       colPDG.push_back(pPDG);
00208       colP.push_back(pMom);
00209       colTH.push_back(lastTH);
00210       colCS.push_back(lastCS);
00211 #ifdef pdebug
00212       G4cout<<"G4QNMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
00213       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00214 #endif
00215       return lastCS*millibarn;
00216     } // End of creation of the new set of parameters
00217     else
00218     {
00219 #ifdef pdebug
00220       G4cout<<"G4QNMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00221 #endif
00222       colP[lastI]=pMom;
00223       colPDG[lastI]=pPDG;
00224       colCS[lastI]=lastCS;
00225     }
00226   } // End of parameters udate
00227   else if(pEn<=lastTH)
00228   {
00229 #ifdef pdebug
00230     G4cout<<"G4QNMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00231     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00232 #endif
00233     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00234   }
00235   else if(std::fabs(lastP/pMom-1.)<tolerance)
00236   {
00237 #ifdef pdebug
00238     G4cout<<"G4QNMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00239     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00240 #endif
00241     return lastCS*millibarn;     // Use theLastCS
00242   }
00243   else
00244   {
00245 #ifdef pdebug
00246     G4cout<<"G4QNMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00247 #endif
00248     lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00249     lastP=pMom;
00250   }
00251 #ifdef pdebug
00252   G4cout<<"G4QNMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00253   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00254 #endif
00255   return lastCS*millibarn;
00256 }

G4double G4QNuNuNuclearCrossSection::GetDirectPart ( G4double  Q2  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 725 of file G4QNuNuNuclearCrossSection.cc.

00726 {
00727   G4double f=Q2/4.62;
00728   G4double ff=f*f;
00729   G4double r=ff*ff;
00730   G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
00731   //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
00732   return 1.-s_value*(1.-s_value/2);
00733 }

G4int G4QNuNuNuclearCrossSection::GetExchangePDGCode (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 742 of file G4QNuNuNuclearCrossSection.cc.

00742 {return 22;}

G4double G4QNuNuNuclearCrossSection::GetLastQELCS (  )  [inline, virtual]

Reimplemented from G4VQCrossSection.

Definition at line 82 of file G4QNuNuNuclearCrossSection.hh.

00082 {return lastQEL;}

G4double G4QNuNuNuclearCrossSection::GetLastTOTCS (  )  [inline, virtual]

Reimplemented from G4VQCrossSection.

Definition at line 81 of file G4QNuNuNuclearCrossSection.hh.

00081 {return lastSig;}

G4double G4QNuNuNuclearCrossSection::GetNPartons ( G4double  Q2  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 736 of file G4QNuNuNuclearCrossSection.cc.

00737 {
00738   return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
00739 }

G4double G4QNuNuNuclearCrossSection::GetNQE_ExchangeQ2 (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 501 of file G4QNuNuNuclearCrossSection.cc.

References G4UniformRand.

00502 {
00503   static const double mpi=.13957018;    // charged pi meson mass in GeV
00504   static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
00505   static const double dMN=MN+MN;        // 2*M_N in GeV
00506   static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
00507   static const G4int power=7;           // direct power for the magic variable
00508   static const G4double pconv=1./power; // conversion power for the magic variable
00509   static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
00510   static const G4int lX=nX-1;           // index of the last in the Xl table
00511   static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
00512   static const G4int nE=20;             // #Of point in the El table (in GeV^2)
00513   static const G4int bE=nE-1;           // index of the last in the El table
00514   static const G4int pE=bE-1;           // index of the before last in the El table
00515   // Reversed table
00516   static const G4double X0[nX]={6.14081e-05,
00517  .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
00518  2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
00519   static const G4double X1[nX]={.00125268,
00520  .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
00521  4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
00522   static const G4double X2[nX]={.015694,
00523  1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
00524  10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
00525   static const G4double X3[nX]={.0866877,
00526  4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
00527  21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
00528   static const G4double X4[nX]={.160483,
00529  5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
00530  29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
00531   static const G4double X5[nX]={.0999307,
00532  5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
00533  26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
00534   static const G4double X6[nX]={.0276367,
00535  3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
00536  16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
00537   static const G4double X7[nX]={.00472383,
00538  2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
00539  9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
00540   static const G4double X8[nX]={.000630783,
00541  1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
00542  5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
00543   static const G4double X9[nX]={7.49179e-05,
00544  .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
00545  3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
00546   static const G4double XA[nX]={8.43437e-06,
00547  .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
00548  2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
00549   static const G4double XB[nX]={9.27028e-07,
00550  .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
00551  1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
00552   static const G4double XC[nX]={1.00807e-07,
00553  .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
00554  1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
00555   static const G4double XD[nX]={1.09102e-08,
00556  .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
00557  1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
00558   static const G4double XE[nX]={1.17831e-09,
00559  .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
00560  .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
00561   static const G4double XF[nX]={1.27141e-10,
00562  .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
00563  .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
00564   static const G4double XG[nX]={1.3713e-11,
00565  .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
00566  .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
00567   static const G4double XH[nX]={1.47877e-12,
00568  .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
00569  .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
00570   static const G4double XI[nX]={1.59454e-13,
00571  .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
00572  .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
00573   static const G4double XJ[nX]={1.71931e-14,
00574  .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
00575  .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
00576   // Direct table
00577   static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
00578                         X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
00579   static const G4double dX[nE]={
00580     (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
00581     (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
00582     (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
00583     (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
00584     (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
00585   static const G4double* Xl[nE]=
00586                              {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
00587   static const G4double I0[nX]={0,
00588  .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
00589  13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
00590   static const G4double I1[nX]={0,
00591  .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
00592  13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
00593   static const G4double I2[nX]={0,
00594  .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
00595  12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
00596   static const G4double I3[nX]={0,
00597  .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
00598  12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
00599   static const G4double I4[nX]={0,
00600  .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
00601  11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
00602   static const G4double I5[nX]={0,
00603  .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
00604  10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
00605   static const G4double I6[nX]={0,
00606  .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
00607  9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
00608   static const G4double I7[nX]={0,
00609  .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
00610  8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
00611   static const G4double I8[nX]={0,
00612  .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
00613  8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
00614   static const G4double I9[nX]={0,
00615  .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
00616  8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
00617   static const G4double IA[nX]={0,
00618  .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
00619  7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
00620   static const G4double IB[nX]={0,
00621  .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
00622  7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
00623   static const G4double IC[nX]={0,
00624  .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
00625  7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
00626   static const G4double ID[nX]={0,
00627  .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
00628  7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
00629   static const G4double IE[nX]={0,
00630  .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
00631  7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
00632   static const G4double IF[nX]={0,
00633  .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
00634  7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
00635   static const G4double IG[nX]={0,
00636  .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
00637  7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
00638   static const G4double IH[nX]={0,
00639  .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
00640  7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
00641   static const G4double II[nX]={0,
00642  .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
00643  7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
00644   static const G4double IJ[nX]={0,
00645  .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
00646  6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
00647   static const G4double* Il[nE]=
00648                              {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
00649   static const G4double lE[nE]={
00650 -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
00651  2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
00652   static const G4double lEmi=lE[0];
00653   static const G4double lEma=lE[nE-1];
00654   static const G4double dlE=(lEma-lEmi)/bE;
00655   //***************************************************************************************
00656   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00657   G4double lEn=std::log(Enu);         // log(E) for interpolation
00658   G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
00659   G4int fE=static_cast<int>(rE);      // Left bin for interpolation
00660   if(fE<0) fE=0;
00661   if(fE>pE)fE=pE;
00662   G4int    sE=fE+1;                   // Right bin for interpolation
00663   G4double dE=rE-fE;                  // relative log shift from the left bin
00664   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00665   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00666   G4double Emu=Enu;                   // Free Energy of neutrino/anti-neutrino
00667   G4double ME=Enu*MN;                 // M*E
00668   G4double dME=ME+ME;                 // 2*M*E
00669   G4double dEMN=(dEnu+MN)*ME;
00670   G4double sqE=Enu*ME;
00671   G4double E2M=MN*Enu2;
00672   G4double ymax=(E2M+sqE)/dEMN;
00673   G4double Q2mi=0.;                   // Q2_min(E_nu)
00674   G4double Q2ma=dME*ymax;             // Q2_max(E_nu)
00675   G4double Q2nq=Emu*dMN-mcV;
00676   if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
00677   // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
00678   G4double Rmi=Q2mi/Q2ma;
00679   G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
00680   // --- E-interpolation must be done in a log scale ---
00681   G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
00682   G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
00683   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00684   G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
00685   G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
00686   G4double rXi=(Xmi-iXmi)/idX;
00687   G4int    iXi=static_cast<int>(rXi);
00688   if(iXi<0) iXi=0;
00689   if(iXi>bX) iXi=bX;
00690   G4double dXi=rXi-iXi;
00691   G4double bntil=Il[fE][iXi];
00692   G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
00693   G4double bntir=Il[sE][iXi];
00694   G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
00695   G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
00696   //
00697   G4double rXa=(Xma-iXmi)/idX;
00698   G4int    iXa=static_cast<int>(rXa);
00699   if(iXa<0) iXa=0;
00700   if(iXa>bX) iXa=bX;
00701   G4double dXa=rXa-iXa;
00702   G4double bntal=Il[fE][iXa];
00703   G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
00704   G4double bntar=Il[sE][iXa];
00705   G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
00706   G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
00707   //
00708   // *** Find X using the reversed table ***
00709   G4double intx=inti+(inta-inti)*G4UniformRand(); 
00710   G4int    intc=static_cast<int>(intx);
00711   if(intc<0) intc=0;
00712   if(intc>bX) intc=bX;
00713   G4double dint=intx-intc;
00714   G4double mXl=Xl[fE][intc];
00715   G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
00716   G4double mXr=Xl[sE][intc];
00717   G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
00718   G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
00719   G4double R=shift-std::pow(X,pconv);
00720   G4double Q2=R*Q2ma;
00721   return Q2*GeV*GeV;
00722 }

G4VQCrossSection * G4QNuNuNuclearCrossSection::GetPointer (  )  [static]

Definition at line 72 of file G4QNuNuNuclearCrossSection.cc.

Referenced by G4QInelastic::GetMeanFreePath(), and G4QInelastic::PostStepDoIt().

00073 {
00074   static G4QNuNuNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
00075   return &theCrossSection;
00076 }

G4double G4QNuNuNuclearCrossSection::GetQEL_ExchangeQ2 (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 424 of file G4QNuNuNuclearCrossSection.cc.

References G4UniformRand.

00425 {
00426   static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
00427   static const G4double power=-3.5;    // direct power for the magic variable
00428   static const G4double pconv=1./power;// conversion power for the magic variable
00429   static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
00430   static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
00431   static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
00432   // Reversed table
00433   static const G4double Xl[nQ2]={1.87905e-10,
00434  .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
00435  .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
00436  .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
00437  .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
00438  .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
00439  .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
00440  .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
00441  .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
00442  .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
00443  .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
00444   // Direct table
00445   static const G4double Xmax=Xl[lQ2];
00446   static const G4double Xmin=Xl[0];
00447   static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
00448   static const G4double inl[nQ2]={0,
00449  1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
00450  15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
00451  26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
00452  36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
00453  46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
00454  55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
00455  65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
00456  74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
00457  83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
00458  92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
00459   G4double Enu=lastE;                 // Get energy of the last calculated cross-section
00460   G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
00461   G4double Enu2=Enu*Enu;              // squared energy of nu/anu
00462   G4double ME=Enu*MN;                 // M*E
00463   G4double dME=ME+ME;                 // 2*M*E
00464   G4double dEMN=(dEnu+MN)*ME;
00465   G4double sqE=Enu*ME;
00466   G4double E2M=MN*Enu2;
00467   G4double ymax=(E2M+sqE)/dEMN;
00468   G4double Q2mi=0; // Q2_min(E_nu)
00469   G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
00470   G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
00471   G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
00472   // Find the integral values integ(Xmi) & integ(Xma) using the direct table
00473   G4double rXi=(Xmi-Xmin)/dX;
00474   G4int    iXi=static_cast<int>(rXi);
00475   if(iXi<0) iXi=0;
00476   if(iXi>bQ2) iXi=bQ2;
00477   G4double dXi=rXi-iXi;
00478   G4double bnti=inl[iXi];
00479   G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
00480   //
00481   G4double rXa=(Xma-Xmin)/dX;
00482   G4int    iXa=static_cast<int>(rXa);
00483   if(iXa<0) iXa=0;
00484   if(iXa>bQ2) iXa=bQ2;
00485   G4double dXa=rXa-iXa;
00486   G4double bnta=inl[iXa];
00487   G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
00488   // *** Find X using the reversed table ***
00489   G4double intx=inti+(inta-inti)*G4UniformRand();
00490   G4int    intc=static_cast<int>(intx);
00491   if(intc<0) intc=0;
00492   if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
00493   G4double dint=intx-intc;
00494   G4double mX=Xl[intc];
00495   G4double X=mX+dint*(Xl[intc+1]-mX);
00496   G4double Q2=std::pow(X,pconv)-1.;
00497   return Q2*GeV*GeV;
00498 }

G4double G4QNuNuNuclearCrossSection::ThresholdEnergy ( G4int  Z,
G4int  N,
G4int  PDG = 14 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 259 of file G4QNuNuNuclearCrossSection.cc.

Referenced by GetCrossSection().

00259 {return 0.;}


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:13 2013 for Geant4 by  doxygen 1.4.7