QwGeant4
QweakSimEPEvent Class Reference

#include <QweakSimEPEvent.hh>

+ Collaboration diagram for QweakSimEPEvent:

Public Member Functions

 QweakSimEPEvent (QweakSimUserInformation *myUI)
 
virtual ~QweakSimEPEvent ()
 
void SetActiveOctantNumber (G4int kaot)
 
G4int GetActiveOctantNumber ()
 
void SetIsotropy (G4int isot)
 
G4int GetIsotropy ()
 
void SetEPrime_Min (G4double energy)
 
G4double GetEPrime_Min ()
 
void SetEPrime_Max (G4double energy)
 
G4double GetEPrime_Max ()
 
void SetThetaAngle_Min (G4double ang)
 
G4double GetThetaAngle_Min ()
 
void SetThetaAngle_Max (G4double ang)
 
G4double GetThetaAngle_Max ()
 
void SetPhiAngle_Min (G4double ang)
 
G4double GetPhiAngle_Min ()
 
void SetPhiAngle_Max (G4double ang)
 
G4double GetPhiAngle_Max ()
 
void SetLuminosity (G4double lum)
 
G4double GetLuminosity ()
 
void SetPhaseSpace (G4double ps)
 
G4double GetPhaseSpace ()
 
void SetBeamEnergy (G4double energy=1.160 *GeV)
 
G4double GetBeamEnergy ()
 
void SetElasticPeakDeltaE (G4double energy=15 *MeV)
 
G4double GetElasticPeakDeltaE ()
 
void SetSchwingerDeltaE (G4double energy=15 *MeV)
 
G4double GetSchwingerDeltaE ()
 
void SetReactionType (G4int rt)
 
G4int GetReactionType ()
 
void SetReactionRegion (G4int rr)
 
G4int GetReactionRegion ()
 
G4double GetVertexZ ()
 
void GetanEvent (G4double E_in, std::vector< G4double > &CrossSection, G4double &weight_n, G4double &Q2, G4double &E_out, G4ThreeVector &MomentumDirection, G4double &theta, G4double &phi, G4double &Asymmetry)
 
G4double Elastic_Cross_Section_Proton (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Elastic_Cross_Section_Aluminum (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Quasi_Elastic_Neutron (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Delta_Resonance (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Moller_Scattering (G4double E_in, G4double theta1, G4double &E_out1, G4double &E_out2, G4double &theta2, G4double &q2, G4double &fWeightN, G4double &asymmetry)
 
const std::vector< G4double > Radiative_Cross_Section_Lookup (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
void CreateLookupTable ()
 
void CheckLookupTableBounds ()
 
G4double AlNuclInel (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double AlGDR (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Elastic_Cross_Section_Carbon (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double CGDR (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double CNuclInel (G4double E_in, G4double Theta, G4int nState, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Pion_PhotoProduction (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Pion_PhotoProductionAl (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Pion_PhotoProductionCarbon (G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Quasi_Elastic_Bosted (G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double NuclearInelastic_Bosted (G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Aluminum_Excited_State (G4double E_in, G4double Theta, G4double E_lvl, G4double fit_c, G4double fit_mean, G4double fit_sigma, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Horowitz_DW_Xsect (G4double E_in, G4double theta, const G4String path, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Horowitz_DW_Asym (G4double theta, const G4String path)
 
void F1F2QE09 (G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
 
void F1F2IN09 (G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
 
void christy507 (G4double wsq, G4double Q2, G4double &F1, G4double &R, G4double &sigt, G4double &sigl)
 
void resmodd (G4double w2, G4double q2, G4double xval[50], G4double &sig)
 
G4double resmod507_v2 (G4double sf, G4double w2, G4double q2, G4double xval[50])
 
G4double MEC2009 (G4int a, G4double q2, G4double w2)
 
G4double fitemc (G4double X, G4int A)
 
G4double AlloyScattering (G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
 
G4double Fshell (G4int Z, G4int A, G4double q2)
 
G4double Fgauss (G4int Z, G4int A, G4double q2)
 
G4double FF_BESSEL (G4int Z, G4int A, G4double q2, G4bool &ofr)
 
G4double GetAsymmetry_EP (G4double theta, G4double energy)
 
G4double GetAsymmetry_EN (G4double theta, G4double energy)
 
G4double GetAsymmetry_AL (G4double theta, G4double energy)
 
G4double GetAsymmetry_Be (G4double theta, G4double energy)
 
G4double GetAsymmetry_Pi (G4double Q2_pi)
 
G4double Born_Approx_Nuclei_Asym (G4double Z, G4double A, G4double energy, G4double eprime, G4double qq)
 

Private Types

typedef float value_t
 

Private Member Functions

G4ThreeVector GetMomentumDirection ()
 
G4int SuperElasticCheck (G4double E_in, G4double E_out, G4double theta, G4double &xsec)
 
G4double ResMod507 (G4int sf, G4double w2, G4double q2, G4double *xval)
 
G4double Sigma_EEPrime (G4double eni, G4double eprime, G4double theta, G4double &q2)
 

Private Attributes

G4int TypeSetting
 
G4int ReactionType
 reaction type used for event generation More...
 
G4int ReactionRegion
 reaction region used for event generation More...
 
G4int kActiveOctantNumber
 Active octant number in the simulation, 0 will enable all octants. More...
 
G4int Isotropy
 isotropy used for event generation More...
 
G4double ElasticPeakDeltaE
 
G4double SchwingerDeltaE
 
G4double BeamEnergy
 
G4double myPositionZ
 
QweakSimFieldMap< value_t,
value_n > * 
fLookupTable
 
QweakSimEPEventMessengerEventGen_Messenger
 
QweakSimUserInformationmyUserInfo
 

Static Private Attributes

static const G4double TargetLength
 
static const G4double mil
 
static const G4double TargetWindowThickness
 
static const G4double M_n = 939.5656*MeV
 
static const G4double M_p = 938.2796*MeV
 
static const G4double Theta_Min = 1.745329E-4
 
static const G4int value_n = 15
 
static const G4int value_d = 4
 

Detailed Description

Definition at line 35 of file QweakSimEPEvent.hh.

Member Typedef Documentation

typedef float QweakSimEPEvent::value_t
private

Definition at line 71 of file QweakSimEPEvent.hh.

Constructor & Destructor Documentation

QweakSimEPEvent::QweakSimEPEvent ( QweakSimUserInformation myUI)

Definition at line 42 of file QweakSimEPEvent.cc.

References ElasticPeakDeltaE, EventGen_Messenger, Isotropy, kActiveOctantNumber, myUserInfo, ReactionRegion, ReactionType, SchwingerDeltaE, QweakSimUserInformation::SetEPEvent(), and TypeSetting.

43 : fLookupTable(0)
44 {
45  G4cout << "###### Calling QweakSimEPEvent::QweakSimEPEvent () " << G4endl;
46 
47  Isotropy = 1;
48 
49  //SetPhiAngle_Min(-16.0*degree);
50  //SetPhiAngle_Max(16.0*degree);
51 
52  //SetThetaAngle_Min(4.0*degree);
53  //SetThetaAngle_Max(13.5*degree);
54 
55  //SetEPrime_Max(1.159*GeV);
56  //SetEPrime_Min(0.059*GeV);
57 
58  //SetBeamEnergy(1.16*GeV);
59 
60  ElasticPeakDeltaE = 15*MeV;
61  SchwingerDeltaE = 15*MeV;
62 
63  TypeSetting = 1;
64  ReactionType = 1;
65  ReactionRegion = 1;
66  kActiveOctantNumber = 0; //default octant all, choose from [0,8]
67 
68  myUserInfo = myUI;
69  myUserInfo->SetEPEvent(this);
70 
72 
73  G4cout << "###### Leaving QweakSimEPEvent::QweakSimEPEvent () " << G4endl;
74 }
G4int ReactionRegion
reaction region used for event generation
QweakSimEPEventMessenger * EventGen_Messenger
G4double SchwingerDeltaE
G4int ReactionType
reaction type used for event generation
QweakSimFieldMap< value_t, value_n > * fLookupTable
UI control for the event generator.
G4double ElasticPeakDeltaE
void SetEPEvent(QweakSimEPEvent *EP)
G4int Isotropy
isotropy used for event generation
G4int kActiveOctantNumber
Active octant number in the simulation, 0 will enable all octants.
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

QweakSimEPEvent::~QweakSimEPEvent ( )
virtual

Definition at line 78 of file QweakSimEPEvent.cc.

References EventGen_Messenger, QweakSimUserInformation::GetNumberOfEventToBeProcessed(), and myUserInfo.

79 {
80 
81  G4cout << "###### Calling QweakSimEPEvent::~QweakSimEPEvent () " << G4endl;
82 
84 
85  G4cout << "Created a great many events: " << myUserInfo->GetNumberOfEventToBeProcessed() << G4endl;
86 
87  G4cout << "###### Leaving QweakSimEPEvent::~QweakSimEPEvent () " << G4endl;
88 }
QweakSimEPEventMessenger * EventGen_Messenger
G4int GetNumberOfEventToBeProcessed() const
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

Member Function Documentation

G4double QweakSimEPEvent::AlGDR ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 1934 of file QweakSimEPEvent.cc.

References Horowitz_DW_Xsect(), and M_p.

1938  {
1939 
1940  G4double M_p = 938.272081; // Proton mass [MeV] - 2016 PDG
1941  G4double Z_Al = 13.0;
1942  G4double N_Al = 14.0;
1943  G4double A_Al = 27.0;
1944  G4double M_Al = M_p*A_Al; // [MeV]
1945 
1946  G4double CTH = cos(Theta/2.0);
1947  G4double STH = sin(Theta/2.0);
1948 
1949  // now get qsq
1950  G4double ETA = 1.0+2.0*E_in*STH*STH/M_Al;
1951  E_out = E_in/ETA;
1952  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
1953 
1954  // GDR peak energy
1955  G4double E_peak = 20.8; // [MeV]
1956  G4double fac = 9.0; // FWHM [MeV]
1957  //Use fac to scale xsection to reflect energy integration.
1958 
1959  G4double mu = ((N_Al*Z_Al)/A_Al)*M_Al; //Reduced mass [MeV]
1960 
1961  G4double long_ratio = pow(N_Al/A_Al,2.0)/(2.0*mu)*Q2/E_peak;
1962  G4double tran_ratio = 0.5*pow(N_Al/A_Al,2.0)*(E_peak/mu)*((1.0+STH*STH)/(CTH*CTH));
1963 
1964  G4double elasticAl = Horowitz_DW_Xsect(E_in,
1965  Theta,
1966  "horowitz_alloy_tables/elasticAl27.csv",
1967  fWeightN,
1968  Q2,
1969  E_out);
1970 
1971  G4double xsect = elasticAl*fac*(long_ratio + tran_ratio); // [ub/sr]
1972 
1973  fWeightN = xsect*sin(Theta);
1974 
1975  //G4cerr << "-----CrossSection: " << xsect << G4endl;
1976 
1977  return xsect;
1978 }
G4double Horowitz_DW_Xsect(G4double E_in, G4double theta, const G4String path, G4double &fWeightN, G4double &Q2, G4double &E_out)
static const G4double M_p

+ Here is the call graph for this function:

G4double QweakSimEPEvent::AlloyScattering ( G4double  E_in,
G4double  Theta,
G4int  Zin,
G4int  Ain,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

use if FF_BESSEL out of range

Definition at line 3706 of file QweakSimEPEvent.cc.

References FF_BESSEL(), Fgauss(), Fshell(), M_p, and SchwingerDeltaE.

Referenced by CGDR().

3713 {
3714  G4double M_p = 938.2796 * MeV;
3715  G4double M = M_p * Ain;
3716  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
3717  G4double m_e = 0.511 * MeV; // electron mass in MeV
3718 
3719 
3720  G4double STH = sin(Theta/2.);
3721  G4double CTH = cos(Theta/2.);
3722 
3723  E_out = E_in/(1+2*E_in/M*STH*STH);
3724  Q2 = 4.*E_in*E_out*STH*STH; // MeV^2
3725 
3726 
3727  // Mott cross section
3728  G4double sigMott = pow(Zin*0.719982/E_in*CTH/(STH*STH),2);
3729  sigMott = sigMott/(1+2*E_in/M*STH*STH)*10000; // 10000 is conversion from fm^2 to ub
3730 
3731  // Form factor
3732  G4int NUC_MODEL = 1; // USE Fourier-Bessel if available
3733  G4bool OUT_OF_RANGE = 0;
3734  //G4double eng = E_in/1e-3; // GeV // unused
3735  G4double qsq = Q2/1e6; // GeV^2
3736 
3737  G4double FF = 1;
3738 
3739  if (Ain < 3) {
3740  G4cout << "** No generator fo A < 3 ** " << G4endl;
3741  return 0;
3742  }
3743  else if (Ain >= 3) {
3744  OUT_OF_RANGE = 0;
3745  if (NUC_MODEL == 1) {
3746  FF = FF_BESSEL(Zin,Ain,qsq,OUT_OF_RANGE);// !only for some Nuclei,for limited Q2
3747  }
3748  if (OUT_OF_RANGE || (NUC_MODEL == 0)) { //!use if FF_BESSEL out of range
3749  if (Ain == 3) { // !3HE ! added 5/30/96 SER
3750  G4cout << "** No generator for He3 **" << G4endl;
3751  return 0;
3752  }
3753  else if (Ain <= 20) {
3754  FF = Fshell(Zin,Ain,qsq);
3755  }
3756  else { // !ia >20
3757  FF = Fgauss(Zin,Ain,qsq);
3758  }
3759  }
3760  }
3761 
3762  G4double xsec = sigMott*FF*FF;
3763  // Schwinger correction from Mo and Tsai Rev.Mod.Phys. 41 205 (1969) eq. (II.2)
3764  G4double delta_Schwinger = 0.0;
3765  if(SchwingerDeltaE!=0) {
3766  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
3767  delta_Schwinger = (-2.0*alpha_FS/pi) * ((log(E_in/SchwingerDeltaE) - 13.0/12.0)
3768  * (log(Q2/(m_e*m_e)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
3769  }
3770  xsec = xsec*(1.0 + delta_Schwinger);
3771  fWeightN = xsec*sin(Theta);
3772 
3773  //G4cout << E_in << "\t" << Theta << "\t" << Q2 << "\t"
3774  // << SchwingerDeltaE << "\t" << delta_Schwinger << G4endl;
3775 
3776  return xsec;
3777 }
G4double SchwingerDeltaE
G4double Fshell(G4int Z, G4int A, G4double q2)
static const G4double M_p
G4double FF_BESSEL(G4int Z, G4int A, G4double q2, G4bool &ofr)
G4double Fgauss(G4int Z, G4int A, G4double q2)

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::AlNuclInel ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 1894 of file QweakSimEPEvent.cc.

References Alinel_crsec_gen(), and M_p.

1899 {
1900  ///
1901  G4double M_p = 938.2796; // proton mass in MeV
1902 //G4double Z_Al = 13.0;
1903  G4double A_Al = 27.0;
1904  G4double M_Al = M_p*A_Al;
1905 
1906 //G4double CTH = cos(Theta/2.0);
1907  G4double STH = sin(Theta/2.0);
1908 
1909  // now get qsq
1910  G4double ETA = 1.0+2.0*E_in*STH*STH/M_Al;
1911  E_out = E_in/ETA;
1912 
1913  Q2 = 4*E_in*E_out*STH*STH; //unit: MeV^2
1914 
1915  Theta *= 180/3.14159; //deg
1916  // in ub/sr
1917  G4double xsect = Alinel_crsec_gen(E_in,Theta);
1918 
1919  Theta *= 3.14159/180; // bk to radians
1920 
1921  fWeightN = xsect*sin(Theta);
1922 
1923  return xsect;
1924 }
static const G4double M_p
G4double Alinel_crsec_gen(G4double E_in=1160, G4double th=8.0)

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Aluminum_Excited_State ( G4double  E_in,
G4double  Theta,
G4double  E_lvl,
G4double  fit_c,
G4double  fit_mean,
G4double  fit_sigma,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2377 of file QweakSimEPEvent.cc.

References M_n, and M_p.

2385  {
2386  G4double xsect = 0.0;
2387  G4double FF2 = 0.0;
2388  const G4double alpha = 1/137.035999139; // Fine Structure Constant 2016-PDG
2389  const G4double M_p = 938.272081; // Proton Mass [MeV] 2016-PDG
2390  const G4double M_n = 939.565413; // Neutron Mass [MeV] 2016-PDG
2391  const G4double m_e = 0.5109989461*MeV; // Electron mass 2016 PDG
2392  const G4double Z_Al = 13.0; // # of Protons
2393  const G4double A_Al = 27.0; // Atomic Mass [amu]
2394  const G4double N_Al = A_Al-Z_Al; // # of Neutrons
2395  const G4double M_Al = M_p*Z_Al + M_n*N_Al; // Atomic Mass [MeV]
2396  const G4double hbar_c = 197.3269788; // Conversion Constant [MeV fm]
2397  const G4double R = 3.94;// Al Uniform-density charge radius [fm]
2398  const G4double deltaE = 20.0*MeV; // Schwinger correction delta E parameter
2399 
2400  G4double STH = sin(Theta/2.0);
2401  G4double CTH = cos(Theta/2.0);
2402 
2403  //Calculate Recoil Factor Eta
2404  G4double ETA = 1.0 + 2.0*(E_in/M_Al)*STH*STH;
2405 
2406  //Calculate Eprime
2407  //See H. Uberall, "Electron Scattering from Complex Nuclei",
2408  //(Academic, NY, 1971) eq.2-10a and eq.2-10d
2409  E_out = (E_in - E_lvl - ((E_lvl*E_lvl)/(2.0*M_Al)))/ETA; //[MeV]
2410 
2411  //Calculate/Return Q^2 and convert to q [fm^-1]
2412  Q2 = 4.0*E_in*E_out*STH*STH; //[MeV^2]
2413  G4double q = sqrt(Q2/(hbar_c*hbar_c)); //[fm^-1]
2414 
2415  //Calculate Q_eff (Coulomb Corrections)
2416  //See Ryan et al. Phys.Rev.C27, 2515 (1983)
2417  G4double q_eff = q*(1.0 + (3.0*Z_Al*alpha*hbar_c)/(2.0*E_in*R)); //[fm^-1]
2418 
2419  //Calculate Mott Cross-section [fm^2/sr]
2420  G4double sigma_mott = pow((Z_Al*alpha*hbar_c*CTH)/(2.0*E_in*STH*STH),2.0)/ETA;
2421  //Convert from [fm^2/sr] to [ubarn/sr]
2422  sigma_mott = 10000.0*sigma_mott;
2423 
2424  //Calculate Gaussian Parameterization of Form Factor Squared
2425  //using q_eff in range 0.3-1.8 [fm^-1]
2426  if((q_eff >= 0.3) && (q_eff <= 1.8)){
2427  FF2 = fit_c*exp(-0.5*pow((q_eff-fit_mean)/(fit_sigma),2.0));
2428  } else{
2429  FF2 = 0.0;
2430  }
2431 
2432  //Calculate Cross-section
2433  xsect = sigma_mott*FF2;
2434 
2435  // Calculate Schwinger Radiative Correction
2436  // Using eq.6-8d pg.488 from "Electron Scattering from Complex Nuclei"
2437  // Part B Herber Uberall - 2017/8/24 K. Bartlett
2438  G4double delta_Schwinger = 0.0;
2439  delta_Schwinger = (TMath::Log(TMath::Sqrt(E_in*E_out)/deltaE)-(13.0/12.0))*(TMath::Log(Q2/(m_e*m_e))-1.0);
2440  delta_Schwinger += (17.0/36.0) + 0.25*TMath::Log(E_in/E_out)*TMath::Log(E_in/E_out) + 0.5*((1.0/6.0)*TMath::Pi()*TMath::Pi()-TMath::DiLog(TMath::Cos(Theta/2.0)*TMath::Cos(Theta/2.0)));
2441  delta_Schwinger *= (2.0*alpha/TMath::Pi());
2442 
2443  xsect *= TMath::Exp(-delta_Schwinger);
2444 
2445 
2446  //Return weighted Cross-Section
2447  fWeightN = xsect*sin(Theta);
2448 
2449  //Debug statement
2450  //G4cout << "E_out: " << E_out << " Q2: " << Q2 << " q_eff: " << q_eff << " FF2: " << FF2 << " Xsect: " << xsect << G4endl;
2451 
2452  return xsect;
2453 }
static const G4double M_p
static const G4double M_n
G4double QweakSimEPEvent::Born_Approx_Nuclei_Asym ( G4double  Z,
G4double  A,
G4double  energy,
G4double  eprime,
G4double  qq 
)

Definition at line 4411 of file QweakSimEPEvent.cc.

4411  {
4412 
4413  //Define constants
4414  const G4double G_F = 1.1663787e-5; //Fermi coupling constant [GeV^-2] PDG-2014
4415  const G4double alpha = 1/137.035999074; //Fine structure constant PDG-2014
4416  const G4double Qw_p = 0.0721; //Horowitz paper
4417  const G4double Qw_n = -0.9878; //Horowitz paper
4418 
4419  //Convert beam energy/outgoing energy from MeV to GeV
4420  energy = energy/1000.0;
4421  eprime = eprime/1000.0;
4422  //Convert Q^2 from MeV^2 to GeV^2
4423  qq = qq/1000000.0;
4424 
4425  G4double A_0 = -G_F/(4*pi*alpha*sqrt(2.0));//[GeV^2]
4426  G4double born_asym = 1.0e6*A_0*qq*(Qw_p+Qw_n*((A-Z)/Z));//[ppm]
4427 
4428  return born_asym;
4429 }
G4double QweakSimEPEvent::CGDR ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2038 of file QweakSimEPEvent.cc.

References AlloyScattering().

2042  {
2043 
2044  // Nuclear Properties (Carbon-12)
2045  G4double M_A = 931.494 * MeV; // a.m.u. in [MeV]
2046  G4double Z = 6.0; // # of Protons
2047  G4double N = 6.0; // # of Neutrons
2048  G4double A = 12.0; // Atomic Weight
2049  G4double M = M_A*A; // Atomic Mass [MeV]
2050  G4double E_peak = 24.0*MeV; // GDR Peak Energy [MeV]
2051  G4double Energy_Integration_Factor = 6.5; // [MeV] FWHM of the GDR
2052  // G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2053  // G4double m_e = 0.511 * MeV; // electron mass in MeV
2054 
2055  // Kinematic Calculations
2056  G4double CTH = cos(Theta/2.0);
2057  G4double STH = sin(Theta/2.0);
2058  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2059  E_out = E_in/ETA-E_peak; // Outgoing electron energy [MeV]
2060  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2061 
2062  // Dynamics
2063  G4double LongRatio = (N*Q2)/(2*A*Z*M*E_peak); // Longitudinal component
2064  G4double TranRatio = (N*E_peak)/(2*A*Z*M)*(1+STH*STH)/(CTH*CTH); // Transverse component
2065  G4double sigma_el = AlloyScattering(E_in, Theta, Z, A, fWeightN, Q2, E_out);
2066  G4double sigma = (LongRatio+TranRatio)*sigma_el*Energy_Integration_Factor; // [ub/sr]
2067 
2068  // Commented out because this is possibly double counting.
2069  // There is allready a Schwinger correction in AlloyScattering.
2070  /* "Schwinger"-like correction for soft photon radiation: from Meister and Griffy PhysRev.133.B1032 (1964)
2071  G4double delta_Schwinger = 0.0;
2072  if(SchwingerDeltaE!=0) {
2073  delta_Schwinger = (alpha_FS/pi)*( log(SchwingerDeltaE*SchwingerDeltaE/(E_in*E_out) ) * ( log(Q2/(m_e*m_e)) - 1.0 )
2074  -0.5*log(E_in/E_out)*log(E_in/E_out) + 13.0/6.0*log(Q2/(m_e*m_e))
2075  - 28.0/9.0);
2076  }
2077  sigma = sigma*(1.0 + delta_Schwinger);
2078  */
2079  fWeightN = sigma*sin(Theta);
2080 
2081  return sigma;
2082 }
G4double AlloyScattering(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)

+ Here is the call graph for this function:

void QweakSimEPEvent::CheckLookupTableBounds ( )

Definition at line 1845 of file QweakSimEPEvent.cc.

References CreateLookupTable(), fLookupTable, GetBeamEnergy(), GetEPrime_Max(), GetEPrime_Min(), GetThetaAngle_Max(), GetThetaAngle_Min(), SetBeamEnergy(), SetEPrime_Max(), SetEPrime_Min(), SetThetaAngle_Max(), and SetThetaAngle_Min().

Referenced by CreateLookupTable(), SetBeamEnergy(), SetEPrime_Max(), SetEPrime_Min(), SetThetaAngle_Max(), and SetThetaAngle_Min().

1846 {
1847  G4cout << "---> Calling CheckLookupTableBounds() <---" << G4endl;
1848  if (fLookupTable == 0) return;
1849  if (GetBeamEnergy() != (G4double)fLookupTable->GetMinimum(1) ) {
1850  G4cout << "!!!! Beam Energy is out of bounds" << G4endl;
1851  G4cout << "---> Changing Beam Energy to default value (" << fLookupTable->GetMinimum(1)/1000 << " GeV)" << G4endl;
1852  SetBeamEnergy( (G4double)fLookupTable->GetMinimum(1) );
1854  return;
1855  }
1856 
1857  if (GetEPrime_Min() < fLookupTable->GetMinimum(2) || GetEPrime_Min() > fLookupTable->GetMaximum(2)) {
1858  SetEPrime_Min(fLookupTable->GetMinimum(2));
1859  G4cout << "!!!! Minimum E\' is out of bounds" << G4endl;
1860  G4cout << "---> Changing Minimum E\' to lower bound (" << fLookupTable->GetMinimum(2)/1000 << " GeV)" << G4endl;
1861  }
1862  if (GetEPrime_Max() > fLookupTable->GetMaximum(2) || GetEPrime_Max() < fLookupTable->GetMinimum(2)) {
1863  SetEPrime_Max(fLookupTable->GetMaximum(2));
1864  G4cout << "!!!! Maximum E\' is out of bounds" << G4endl;
1865  G4cout << "---> Changing Maximum E\' to upper bound (" << fLookupTable->GetMaximum(2)/1000 << " GeV)" << G4endl;
1866  }
1867  if (GetEPrime_Min() > GetEPrime_Max()) {
1869  G4cout << "!!!! Minimum E\' greater than maximum E\'" << G4endl;
1870  G4cout << "---> Changing Minimum E\' to " << GetEPrime_Max()/1000 << " GeV" << G4endl;
1871  }
1872 
1873  if (GetThetaAngle_Min() < fLookupTable->GetMinimum(3) || GetThetaAngle_Min() > fLookupTable->GetMaximum(3)) {
1874  SetThetaAngle_Min(fLookupTable->GetMinimum(3));
1875  G4cout << "!!!! Minimum Theta is out of bounds" << G4endl;
1876  G4cout << "---> Changing Minimum Theta to lower bound (" << fLookupTable->GetMinimum(3)/degree << " degrees)" << G4endl;
1877  }
1878  if (GetThetaAngle_Max() > fLookupTable->GetMaximum(3) || GetThetaAngle_Max() < fLookupTable->GetMinimum(3)) {
1879  SetThetaAngle_Max(fLookupTable->GetMaximum(3));
1880  G4cout << "!!!! Maximum Theta is out of bounds" << G4endl;
1881  G4cout << "---> Changing Maximum Theta to upper bound (" << fLookupTable->GetMaximum(3)/degree << " degrees)" << G4endl;
1882  }
1885  G4cout << "!!!! Minimum Theta greater than maximum Theta" << G4endl;
1886  G4cout << "---> Changing Minimum Theta to " << GetThetaAngle_Max() << " degrees" << G4endl;
1887  }
1888  G4cout << "---> Leaving CheckLookupTableBounds() <---" << G4endl;
1889 }
void SetEPrime_Max(G4double energy)
G4double GetThetaAngle_Max()
QweakSimFieldMap< value_t, value_n > * fLookupTable
void SetEPrime_Min(G4double energy)
G4double GetThetaAngle_Min()
void SetThetaAngle_Min(G4double ang)
G4double GetBeamEnergy()
void SetBeamEnergy(G4double energy=1.160 *GeV)
void SetThetaAngle_Max(G4double ang)
G4double GetEPrime_Max()
G4double GetEPrime_Min()

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::christy507 ( G4double  wsq,
G4double  Q2,
G4double &  F1,
G4double &  R,
G4double &  sigt,
G4double &  sigl 
)

Definition at line 3137 of file QweakSimEPEvent.cc.

References resmod507_v2().

Referenced by F1F2IN09().

3139 {
3140 // M.E. Christy and P.E. Bosted, ``Empirical Fit to Precision
3141 // Inclusive Electron-Proton Cross Sections in the Resonance Region'',
3142 // (arXiv:0712.3731). To be submitted to Phys. Rev. C.
3143 
3144  G4double xval1[50],xvalL[50];
3145 
3146  G4double mp = .9382727;
3147  G4double mp2 = mp*mp;
3148  G4double alpha = 1./137.036;
3149 
3150  G4double xval[100] = {
3151  0.12298E+01,0.15304E+01,0.15057E+01,0.16980E+01,0.16650E+01,
3152  0.14333E+01,0.13573E+00,0.22000E+00,0.82956E-01,0.95782E-01,
3153  0.10936E+00,0.37944E+00,0.77805E+01,0.42291E+01,0.12598E+01,
3154  0.21242E+01,0.63351E+01,0.68232E+04,0.33521E+05,0.25686E+01,
3155  0.60347E+00,0.21240E+02,0.55746E-01,0.24886E+01,0.23305E+01,
3156  -.28789E+00,0.18607E+00,0.63534E-01,0.19790E+01,-.56175E+00,
3157  0.38964E+00,0.54883E+00,0.22506E-01,0.46213E+03,0.19221E+00,
3158  0.19141E+01,0.24606E+03,0.67469E-01,0.13501E+01,0.12054E+00,
3159  -.89360E+02,0.20977E+00,0.15715E+01,0.90736E-01,-.38495E-02,
3160  0.10362E-01,0.19341E+01,0.38000E+00,0.34187E+01,0.14462E+00,
3161  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
3162  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
3163  0.00000E+00,0.00000E+00,0.29414E+02,0.19910E+02,0.22587E+00,
3164  0.00000E+00,0.00000E+00,0.38565E+04,0.65717E+00,0.00000E+00,
3165  0.15792E+03,0.97046E+02,0.31042E+00,0.00000E+00,0.42160E+01,
3166  0.38200E-01,0.12182E+01,0.00000E+00,0.13764E+02,0.31393E+00,
3167  0.29997E+01,0.00000E+00,0.55124E+01,0.53743E-01,0.13091E+01,
3168  0.00000E+00,0.86746E+02,0.40864E-04,0.40294E+01,0.31285E+01,
3169  0.33403E+00,0.49623E+01,0.00000E+00,0.00000E+00,0.11000E+02,
3170  0.18951E+01,0.51376E+00,0.00000E+00,0.42802E+01,0.00000E+00};
3171 
3172  for (G4int i = 0; i<50; i++) {
3173  xval1[i] = xval[i];
3174  xvalL[i] = xval[50+i] ;
3175  if(i < 12) xvalL[i] = xval1[i];
3176  }
3177  xvalL[42] = xval1[46];
3178  xvalL[43] = xval1[47];
3179  xvalL[49] = xval1[49];
3180 
3181  //G4double xb = q2/(w2+q2-mp2);
3182  sigT = resmod507_v2(1,w2,q2,xval1);
3183  sigL = resmod507_v2(2,w2,q2,xvalL);
3184 
3185  F1 = sigT*(w2-mp2)/8./pi/pi/alpha/0.3894e3;
3186  //G4double FL = sigL*2.*xb*(w2-mp2)/8./pi/pi/alpha/0.3894e3;
3187  R = sigL/sigT;
3188 
3189  return;
3190 }
G4double resmod507_v2(G4double sf, G4double w2, G4double q2, G4double xval[50])

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::CNuclInel ( G4double  E_in,
G4double  Theta,
G4int  nState,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2092 of file QweakSimEPEvent.cc.

References SchwingerDeltaE.

2098 {
2099  // Physical Constants
2100  G4double M_A = 931.494 * MeV; // a.m.u. in MeV
2101  G4double Z = 6.0; // # of protons
2102  G4double A = 12.0; // Atomic Weight
2103  G4double M = M_A*A; // Nuclear Mass
2104  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2105  G4double myhbarc = hbarc / MeV / fermi; // 197.3269631 MeV fm
2106  G4double m_e = 0.511 * MeV; // electron mass in MeV
2107  G4double EState[3] = { 4.439*MeV, // Energy Levels
2108  7.543*MeV,
2109  9.641*MeV};
2110 
2111  // Kinematic Calculations
2112  G4double CTH = cos(Theta/2.0);
2113  G4double STH = sin(Theta/2.0);
2114  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2115  E_out = E_in/ETA-EState[nState-1]; // Outgoing electron energy [MeV] subtracted excitation energy
2116  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2117  G4double q = sqrt(Q2/(myhbarc*myhbarc)); // [fm^-1]
2118 
2119  // Mott Cross Section
2120  G4double SigmaMott = pow(alpha_FS*myhbarc*Z/(2.0*E_in),2); // [fm^2]
2121  SigmaMott *= (CTH*CTH)/pow(STH,4)*(E_out/E_in)*10000; // [ub] 10000 is the conversion factor
2122 
2123  // FF^2 only from first three excited states (others are at least 1 order of magnitude smaller)
2124  G4double FF2;
2125  // Fits to published data. Descrbed in Ancillary Elog 257
2126  if(nState==1) { // 4.44 MeV, 2+: PhysRevLett27.745 and PhysRev.136.B1580
2127  TF1* ff2 = new TF1("4439","0.0147*exp(-3.25*(x-1.15)**2)",0.5,1.8);
2128  FF2 = ff2->Eval(q);
2129  delete ff2;
2130  } else if(nState==2) { // 7.65 Mev, 0+: PhysRev.136.B1580
2131  TF1* ff2 = new TF1("7654","0.0358*exp(-5.00*(x-0.96)**2)",0.5,1.8);
2132  FF2 = ff2->Eval(q);
2133  delete ff2;
2134  } else if(nState==3) { // 9.64 MeV, 3-: PhysRev.136.B1580
2135  TF1* ff2 = new TF1("9641","0.0526*exp(-2.45*(x-1.25)**2)",0.5,1.8);
2136  FF2 = ff2->Eval(q);
2137  delete ff2;
2138  } else {
2139  G4cout << "Please select from the first three states. nState = {1,2,3} only" << G4endl;
2140  return 0;
2141  }
2142 
2143  G4double sigma = SigmaMott*FF2;
2144 
2145  // "Schwinger"-like correction for soft photon radiation: from Meister and Griffy PhysRev.133.B1032 (1964)
2146  G4double delta_Schwinger = 0.0;
2147  if(SchwingerDeltaE!=0) {
2148  delta_Schwinger = (alpha_FS/pi)*( log(SchwingerDeltaE*SchwingerDeltaE/(E_in*E_out) ) * ( log(Q2/(m_e*m_e)) - 1.0 )
2149  -0.5*log(E_in/E_out)*log(E_in/E_out) + 13.0/6.0*log(Q2/(m_e*m_e))
2150  - 28.0/9.0);
2151  }
2152  sigma = sigma*(1.0 + delta_Schwinger);
2153 
2154  fWeightN = sigma*sin(Theta);
2155  return sigma;
2156 }
G4double SchwingerDeltaE
void QweakSimEPEvent::CreateLookupTable ( )

Definition at line 1548 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), fLookupTable, GetBeamEnergy(), ReactionRegion, SetBeamEnergy(), and value_d.

Referenced by CheckLookupTableBounds(), QweakSimEPEventMessenger::SetNewValue(), and SetReactionRegion().

1549 {
1550 
1551  G4cout << "---> Calling CreateLookupTable() <---" << G4endl;
1552 
1553  G4double energy = GetBeamEnergy();
1554  if (energy != 3350*MeV && energy != 1160*MeV && energy != 877*MeV) {
1555  G4cout << "#### Current beam energy is not a valid choice for lookup table" << G4endl;
1556  G4cout << "#### Setting beam energy to default value for lookup table (3.35 GeV)" << G4endl;
1557  energy = 3350*MeV;
1558  SetBeamEnergy(energy);
1559  }
1560 
1561  TString target;
1562  if (ReactionRegion == 1) target = "lh2";
1563  else if (ReactionRegion == 2) target = "USALWindow";
1564  else if (ReactionRegion == 3) target = "DSALWindow";
1565  else if (ReactionRegion == 4) target = "USALDummy1";
1566  else if (ReactionRegion == 5) target = "USALDummy2";
1567  else if (ReactionRegion == 6) target = "USALDummy4";
1568  else if (ReactionRegion == 7) target = "DSALDummy2";
1569  else if (ReactionRegion == 8) target = "DSALDummy4";
1570  else if (ReactionRegion == 9) target = "DSALDummy8";
1571  else if (ReactionRegion == 10) target = "USCDummy";
1572  else if (ReactionRegion == 11) target = "DSCDummy";
1573 
1574  TString filename = "radiative_cross_section_";
1575  filename += target;
1576  filename += "_";
1577  filename += (Int_t)energy;
1578  filename += "MeV";
1579  //if (ReactionRegion == 1 && energy != 3350*MeV) filename += "_v2";
1580  if (energy != 3350*MeV) filename += "_v2";
1581  filename += ".dat";
1582 
1583  TString filepath = "./radiative_lookup_tables/";
1584 
1585  std::ifstream in;
1586  in.open(filepath + filename);
1587  if (!in.is_open()) {
1588  G4cout << "#### Failed to open lookup table data file \"" << filepath << filename << "\"" << G4endl;
1589  return;
1590  }
1591  else
1592  {
1593  in.close();
1594  G4cout << "#### Found lookup table data file \"" << filepath << filename << "\"" << G4endl;
1595  if (fLookupTable != 0) delete fLookupTable;
1597  fLookupTable->ReadTextFile((filepath +filename).Data());
1598  }
1599  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1600  G4cout << "++ Lookup Table Parameters ++" << G4endl;
1601  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1602  G4cout << "Beam Energy: " << fLookupTable->GetMinimum(1)/1000 << " GeV" << G4endl;
1603  G4cout << "Target: " << target << G4endl;
1604  G4cout << "E\' Lower Bound: " << fLookupTable->GetMinimum(2)/1000 << " GeV" << G4endl;
1605  G4cout << "E\' Upper Bound: " << fLookupTable->GetMaximum(2)/1000 << " GeV" << G4endl;
1606  G4cout << "Theta Lower Bound: " << fLookupTable->GetMinimum(3)/degree << " degrees" << G4endl;
1607  G4cout << "Theta Upper Bound: " << fLookupTable->GetMaximum(3)/degree << " degrees" << G4endl;
1608  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1609 
1610  if (energy != GetBeamEnergy()) SetBeamEnergy(energy);
1611  else CheckLookupTableBounds();
1612 
1613 
1614  G4cout << "---> Leaving CreateLookupTable() <---" << G4endl;
1615 
1616 
1617  //********************************************************************
1618 
1619  /*
1620  // **********************************************************************
1621  // **********************************************************************
1622  // The code below is used to convert the raw output of the Bosted code
1623  // into a format the GEANT4 can read.
1624  // **********************************************************************
1625 
1626 
1627  //G4int entries = 93610;
1628  G4int entries = 1;
1629  //G4int entries[3] = { 1295, 1665, 1073};
1630 
1631  coord_t coord[value_d] = {0.0};
1632  value_t field[value_n] = {0.0};
1633 
1634  std::vector< G4double > fMin;
1635  std::vector< G4double > fMax;
1636  std::vector< G4double > fStep;
1637 
1638  //std::vector< G4double > fMin[3];
1639  //std::vector< G4double > fMax[3];
1640  //std::vector< G4double > fStep[3];
1641 
1642  std::string target[10] = { "DSALDummy2",
1643  "DSALDummy4",
1644  "DSALDummy8",
1645  "USALDummy1",
1646  "USALDummy2",
1647  "USALDummy4",
1648  "DSALWindow",
1649  "USALWindow",
1650  "DSCDummy",
1651  "USCDummy" };
1652 
1653  std::string energy_in[3] = { "877mev",
1654  "1gev",
1655  "3gev" };
1656 
1657  std::string energy_out[3] = { "877MeV",
1658  "1160MeV",
1659  "3350MeV" };
1660 
1661  std::string filepath = "/group/qweak/jdowd/externals_qweak/OUT/";
1662  std::string name, name1, name2, name3;
1663 
1664  // Set min, max, and step size for each coordinate in the lookup table.
1665  // Z position in radiation lengths
1666  fMin.push_back(myUserInfo->TargetCenterPositionZ - 0.5*myUserInfo->TargetLength);
1667  fMax.push_back(myUserInfo->TargetCenterPositionZ + 0.5*myUserInfo->TargetLength);
1668  fStep.push_back(myUserInfo->TargetLength/10.0);
1669 
1670  //fMin.push_back(0.0);
1671  //fMax.push_back(0.0);
1672  //fStep.push_back(1.0);
1673 
1674  //fMin[0].push_back(0.0);
1675  //fMax[0].push_back(0.0);
1676  //fStep[0].push_back(1.0);
1677 
1678  //fMin[1].push_back(0.0);
1679  //fMax[1].push_back(0.0);
1680  //fStep[1].push_back(1.0);
1681 
1682  //fMin[2].push_back(0.0);
1683  //fMax[2].push_back(0.0);
1684  //fStep[2].push_back(1.0);
1685 
1686 
1687  // Beam Energy in GeV
1688  fMin.push_back(1.16*GeV);
1689  fMax.push_back(1.16*GeV);
1690  fStep.push_back(0.005*GeV);
1691 
1692  //fMin[0].push_back(0.877*GeV);
1693  //fMax[0].push_back(0.877*GeV);
1694  //fStep[0].push_back(0.05*GeV);
1695 
1696  //fMin[1].push_back(1.16*GeV);
1697  //fMax[1].push_back(1.16*GeV);
1698  //fStep[1].push_back(0.05*GeV);
1699 
1700  //fMin[2].push_back(3.35*GeV);
1701  //fMax[2].push_back(3.35*GeV);
1702  //fStep[2].push_back(0.05*GeV);
1703 
1704  // E prime in GeV
1705  fMin.push_back(0.015*GeV);
1706  fMax.push_back(1.16*GeV);
1707  fStep.push_back(0.005*GeV);
1708 
1709  //fMin[0].push_back(0.026*GeV);
1710  //fMax[0].push_back(0.876*GeV);
1711  //fStep[0].push_back(0.025*GeV);
1712 
1713  //fMin[1].push_back(0.059*GeV);
1714  //fMax[1].push_back(1.159*GeV);
1715  //fStep[1].push_back(0.025*GeV);
1716 
1717  //fMin[2].push_back(0.150*GeV);
1718  //fMax[2].push_back(1.550*GeV);
1719  //fStep[2].push_back(0.050*GeV);
1720 
1721  // Theta angle
1722  fMin.push_back(2.00*degree);
1723  fMax.push_back(20.0*degree);
1724  fStep.push_back(0.5*degree);
1725 
1726  //fMin[0].push_back(2.00*degree);
1727  //fMax[0].push_back(20.0*degree);
1728  //fStep[0].push_back(0.5*degree);
1729 
1730  //fMin[1].push_back(2.00*degree);
1731  //fMax[1].push_back(20.0*degree);
1732  //fStep[1].push_back(0.5*degree);
1733 
1734  //fMin[2].push_back(2.00*degree);
1735  //fMax[2].push_back(20.0*degree);
1736  //fStep[2].push_back(0.5*degree);
1737 
1738  //G4cout << "Target Min: " << fMin[0] << G4endl;
1739  //G4cout << "Target Max: " << fMax[0] << G4endl;
1740 
1741  for (Int_t n = 0; n < (Int_t)fStep.size(); n++) {
1742  entries *= (G4int)( (fMax[n]-fMin[n])/fStep[n] + 1.5 );
1743  }
1744  std::ifstream in;
1745 
1746  G4cout << "!!!!!!!!!!!!!!!!!!!!!! Debug !!!!!!!!!!!!!!!!!!!!!!" << G4endl;
1747  G4cout << " fMin: " << fMin[0] << " " << fMin[1] << " " << fMin[2] << " " << fMin[3] << G4endl;
1748  G4cout << " fMax: " << fMax[0] << " " << fMax[1] << " " << fMax[2] << " " << fMax[3] << G4endl;
1749  G4cout << " fStep: " << fStep[0]<< " " << fStep[1]<< " " << fStep[2]<< " " << fStep[3]<< G4endl;
1750  G4cout << "===================================================================" << G4endl;
1751 
1752  //for (Int_t target_num = 0; target_num < 10; target_num++) {
1753  for (Int_t target_num = 0; target_num < 1; target_num++) {
1754  //for (Int_t energy_num = 0; energy_num < 3; energy_num++) {
1755  for (Int_t energy_num = 0; energy_num < 1; energy_num++) {
1756  //in.open("./radiative_lookup.dat");
1757 
1758  name = filepath + "lookup_" + target[target_num] + "_" + energy_out[energy_num] + ".out";
1759 
1760  in.open("./test_v2.out");
1761  //in.open(name.c_str());
1762  if (!in.is_open())
1763  G4cout << "#### Failed to open data file for lookup table" << G4endl;
1764  else
1765  {
1766  G4cout << "#### Found lookup table data file" << G4endl;
1767  if (fLookupTable != 0) delete fLookupTable;
1768  fLookupTable = new QweakSimFieldMap<value_t,value_n>(value_d);
1769  //fLookupTable->SetMinimumMaximumStep(fMin[energy_num],fMax[energy_num],fStep[energy_num]);
1770  fLookupTable->SetMinimumMaximumStep(fMin,fMax,fStep);
1771 
1772 
1773  //for (Int_t n = 0; n < (Int_t)fStep[energy_num].size(); n++) {
1774  // entries *= (G4int)( (fMax[energy_num][n]-fMin[energy_num][n])/fStep[energy_num][n] + 1.5 );
1775  //}
1776 
1777  // Filling Lookup Table
1778  //for (Int_t line = 0; line < entries[energy_num]; line++) {
1779  for (Int_t line = 0; line < entries; line++) {
1780  if (in.peek() == EOF) {
1781  G4cout << "#### Error reading \'elastic_lookup.dat\': File contains only "
1782  << line << " of " << entries << " expected lines. ####" << G4endl;
1783  break;
1784  }
1785 
1786  for (G4int i = 0; i < value_d; i++) in >> coord[i];
1787  for (G4int j = 0; j < value_n; j++) in >> field[j];
1788 
1789  //G4cout << "Initial: " << coord[0] << " " << coord[1] << " " << coord[2] << " " << coord[3] << G4endl;
1790 
1791  // Add units
1792  // Z position converted from rad lengths to GEANT coords
1793  coord[0] = myUserInfo->TargetCenterPositionZ - 0.5*myUserInfo->TargetLength
1794  + coord[0]*myUserInfo->TargetLength/0.0396;
1795  coord[1] *= GeV; // Beam Energy
1796  coord[2] *= GeV; // E prime
1797  coord[3] *= degree; // Theta
1798 
1799  //G4cout << "Final: " << coord[0] << " " << coord[1] << " " << coord[2] << " " << coord[3] << G4endl;
1800  //G4cout << "-----------------------------------------------------------------" << G4endl;
1801 
1802  //field[0] // Bjorken x (unitless) // This value has too few sig figs
1803  field[1] *= GeV*GeV; // Q2 // This value has too few sig figs
1804  //field[2] // total born cross section (ub/Sr)
1805  //field[3] // inelastic born cross section(ub/Sr)
1806  //field[4] // quasi-elastic born cross section (ub/Sr)
1807  //field[5] // quasi-elastic factor
1808  //field[6] // total radiated cross section (ub/Sr)
1809  //field[7] // elastic radiated cross section (ub/Sr)
1810  //field[8] // quasi-elastic radiated cross section (ub/Sr)
1811  //field[9] // deep-inelastic radiated cross section (ub/Sr)
1812  //field[10] // charge correction
1813  //field[11] // total radiated cross section internal only (ub/Sr)
1814  //field[12] // elastic radiated cross section internal only (ub/Sr)
1815  //field[13] // deep-inelastic radiated cross section internal only (ub/Sr)
1816  //field[14] // quasi-elastic radiated cross section internal only (ub/Sr)
1817 
1818  fLookupTable->Set(coord,field);
1819  }
1820  G4cout << "===== Filling of Lookup Table complete! =====" << G4endl;
1821  }
1822  in.close();
1823  fLookupTable->WriteTextFile("./radiative_lookup_1160_v2.out");
1824 
1825  name = filepath + "final/radiative_lookup_" + target[target_num] + "_" + energy_out[energy_num] + ".dat";
1826 
1827  //fLookupTable->WriteTextFile("./radiative_lookup_test.out");
1828  //fLookupTable->WriteTextFile(name.c_str());
1829 
1830  G4cout << "Wrote: " << name << G4endl;
1831 
1832  }
1833  }
1834  */
1835 }
G4int ReactionRegion
reaction region used for event generation
QweakSimFieldMap< value_t, value_n > * fLookupTable
A multi-dimensional grid of values with interpolation methods.
G4double GetBeamEnergy()
void SetBeamEnergy(G4double energy=1.160 *GeV)
static const G4int value_d

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::Delta_Resonance ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 1008 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetBeamEnergy(), myUserInfo, Sigma_EEPrime(), and Theta_Min.

1013 {
1014  G4double E_beam = myUserInfo->GetBeamEnergy(); // maximum energy (beam energy) in MeV
1015  G4double M_electron = 0.511; // minimum energy (electron mass) in MeV
1016 
1017  // G4double Theta_Min = 1.745329E-4;
1018  if (Theta<Theta_Min)
1019  Theta = Theta_Min;
1020 
1021  // Generate flat energy distribution of outgoing electron
1022  E_out = M_electron + G4UniformRand()*(E_beam - M_electron);
1023 
1024  // TODO: total energy phase space should be reduced to improve the efficiency.
1025  G4double xsect = Sigma_EEPrime(E_in/1000.0, E_out/1000.0, Theta, Q2); // ub/sr/GeV
1026  Q2 = Q2*1e6; // convert to MeV^2 for Q2
1027 
1028  //std::cout<<"Q2="<<Q2<<" xsect="<<xsect<<std::endl;
1029  xsect = xsect*(E_beam - M_electron)/1000.0;
1030  //std::cout<<"xsect*1.165 = "<<xsect<<std::endl;
1031  fWeightN = xsect*sin(Theta);
1032 
1033  if(xsect == 0) // if E > E_max, reject the event
1034  {
1035  E_out = M_electron;
1036  Q2 = 0.0;
1037  }
1038 
1039  return xsect;
1040 }
static const G4double Theta_Min
QweakSimUserInformation * myUserInfo
G4double Sigma_EEPrime(G4double eni, G4double eprime, G4double theta, G4double &q2)

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Elastic_Cross_Section_Aluminum ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 858 of file QweakSimEPEvent.cc.

References M_p, and Theta_Min.

863 {
864  G4double Theta_Min = 1.745329E-4;
865  G4double M_p = 938.2796; // proton mass in MeV
866  G4double Z = 13.0;
867  G4double A = 27.0;
868  G4double M = M_p*A;
869 
870 // harmonic oscillator well parameter a0 ~1.76 fm
871  G4double a = 2.98; //unit: fm
872  G4double ap = sqrt(0.427); //unit: fm
873  G4double a0 = sqrt((a*a-1.5*ap*ap)/(3.5-10/Z-1.5/A));
874 
875  G4double Q = 14.6; //unit fm^(-2)
876 // G4double J = 5.0/2.0;
877 
878 // G4double gamma = 2.792847351; //This is the magnetic moment
879 // G4double Omega = -3/2*(1+2*gamma)*a0*a0;
880 // G4double Gamma = 15/4*gamma*pow(a0,4);
881 // G4double mu = 3.69; // This is NOT the magnetic moment
882 // G4double Omega_mu = -6.36; // Omega/mu
883 // G4double Gamma_mu = 20.7; // Gamma/mu
884 
885  if (Theta<Theta_Min)
886  Theta = Theta_Min;
887 
888 // E_in unit is MeV, q2 unit is fm^(-2)
889  G4double CTH = cos(Theta/2.0);
890  G4double STH = sin(Theta/2.0);
891 // G4double T2THE = STH*STH/CTH/CTH;
892  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
893  E_out = E_in/ETA;
894 
895  Q2 = 4*E_in*E_out*STH*STH; //unit: MeV^2
896  G4double q2 = Q2/1000000*(1.0/0.197)*(1.0/0.197); //convert MeV^2 into fm^(-2)
897  G4double x = (1.0/4.0)*q2*a0*a0;
898 
899 // std::ofstream EventDataFile("Event.dat", std::ios::app);
900 // EventDataFile<<"E_f="<<E_f<<" q2="<<q2<<" "<<"x="<<x<<G4endl;
901 
902 // Electric form factor (data fit)
903  G4double F0 = (1.0/Z)*( Z-4.0/3.0*(Z-5.0)*x+4.0/15.0*(Z-8.0)*x*x)*exp(-x);
904  G4double F2 = (1.0-2.0/7.0*x)*exp(-x);
905  G4double Fe = sqrt( F0*F0+(7.0/450.0)*q2*q2*(Q*Q/Z/Z)*F2*F2 );
906  Fe=Fe*exp(-0.25*q2*ap*ap); //correction for finite proton size
907  Fe=Fe*exp(x/A); //correction for center-of-well motion
908  G4double Fe_2 = pow(Fe,2);
909 
910 
911 // EventDataFile<<"F0="<<F0<<" "<<"F2="<<F2<<" "<<"Fe_2="<<Fe_2<<G4endl;
912 
913 // Magnetic form factor (theoretical)
914 // G4double Fm1 = (1.0-2.0/5.0*(1.0+2.0*gamma)/(1.0+gamma/2.0)*x+6.0/35.0*gamma*x*x/(1.0+gamma/2.0))*exp(-x);
915 // G4double Fm3 = (1.0-2.0/3.0*gamma*x/(1.0+2.0*gamma))*exp(-x);
916 // G4double Fm5 = exp(-x);
917 // G4double Fm_2 = Fm1*Fm1+(4.0/525.0)*pow(q2*Omega_mu*Fm3,2)+(2.0/33075.0)*pow(q2,4)*pow(Gamma_mu*Fm5,2);
918 
919 // EventDataFile<<"Fm1="<<Fm1<<" "<<"Fm3="<<Fm3<<" "<<"Fm5="<<Fm5<<" "<<"Fm_2="<<Fm_2<<G4endl;
920 
921 // form factor square
922 // G4double F_2 = Z*Z*Fe_2+(1.0+2.0*T2THE)*mu*mu/3.0*(J+1.0)/J*q2/4.0/pow((M_p/0.197),2)*Fm_2;
923  G4double F_2 = Fe_2;//+(1.0+2.0*T2THE)*mu*mu/3.0*(J+1.0)/J*q2/4.0/pow((M_p/0.197),2)*Fm_2;
924 
925 // EventDataFile<<"F_2="<<F_2<<" "<<G4endl;
926 
927  //The next lines are to add schwinger
928  G4double Electron_Mass = 0.511 * MeV; //electorn mass in MeV
929  G4double alpha = 1.0/137.035999074;
930 
931  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
932  G4double delta_Schwinger = (-2.0*alpha/pi) * ((log(E_in/15.0) - 13.0/12.0) * (log(Q2/(Electron_Mass*Electron_Mass)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
933 
934 // cross section units: ub/sr
935 // G4double SigmaMott = ((0.72/E_in)*cos(Theta/2)/(sin(Theta/2)^2))^2 /( 1+2*E_in/M*(sin(Theta/2)^2))*10000 ;
936  G4double SigmaMott = pow(((0.72/E_in)*CTH/(STH*STH)),2)/(1+2*E_in/M*STH*STH)*10000 ;
937  SigmaMott *= (Z*Z);
938 
939 // EventDataFile<<"SigmaMott="<<SigmaMott<<" "<<"SigmaMott*F_2="<<SigmaMott*F_2<<" "<<"weight_n="<<SigmaMott*F_2*sin(Theta)<<G4endl;
940 
941  G4double Xsect = SigmaMott*F_2;
942  Xsect = Xsect * (1.0 + delta_Schwinger);
943  fWeightN = Xsect*sin(Theta);
944 
945  // G4cout<< "q(1/fm)"<<sqrt(q2) <<" SigmaMott="<<SigmaMott<<" F_2="<<F_2<<" SigmaMott*F_2="<<Xsect<<G4endl;
946  return Xsect;
947 }
static const G4double M_p
static const G4double Theta_Min
G4double QweakSimEPEvent::Elastic_Cross_Section_Carbon ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 1986 of file QweakSimEPEvent.cc.

References SchwingerDeltaE, and Theta_Min.

1990  {
1991 
1992  // Minumum Scattering angle allowed
1993  G4double Theta_Min = 0.5*pi/180;
1994  if (Theta<Theta_Min) Theta = Theta_Min;
1995 
1996  // Physical Constants
1997  G4double M_A = 931.494 * MeV; // a.m.u. in MeV
1998  G4double Z = 6.0; // # of protons
1999  G4double A = 12.0; // Atomic Weight
2000  G4double M = M_A*A; // Nuclear Mass
2001  G4double a = 1.65; // well parameter [fm]
2002  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2003  G4double m_e = 0.511 * MeV; // electron mass in MeV
2004 
2005  // Kinematic variables
2006  G4double CTH = cos(Theta/2.0);
2007  G4double STH = sin(Theta/2.0);
2008  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2009  E_out = E_in/ETA;
2010  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2011  G4double myhbarc = hbarc / MeV / fermi; // 197.3269631 MeV fm
2012  G4double q2 = Q2/(myhbarc*myhbarc); //convert MeV^2 into fm^(-2)
2013 
2014  // Dynamics
2015  // Mott Cross Section
2016  G4double SigmaMott = pow(alpha_FS*myhbarc*Z/(2.0*E_in),2); // [fm^2]
2017  SigmaMott *= (CTH*CTH)/pow(STH,4)*(E_out/E_in)*10000; // [ub] 10000 is the conversion factor
2018  G4double FF = (1.0 - (q2*a*a)/9.0)*exp(-(q2*a*a)/4.0); // Form Factor
2019  G4double sigma = SigmaMott*FF*FF;
2020  // Schwinger correction from Mo and Tsai Rev.Mod.Phys. 41 205 (1969) eq. (II.2)
2021  G4double delta_Schwinger = 0.0;
2022  if(SchwingerDeltaE!=0) {
2023  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
2024  delta_Schwinger = (-2.0*alpha_FS/pi) * ((log(E_in/SchwingerDeltaE) - 13.0/12.0)
2025  * (log(Q2/(m_e*m_e)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
2026  }
2027  sigma = sigma*(1.0 + delta_Schwinger);
2028  fWeightN = sigma*sin(Theta);
2029 
2030  return sigma;
2031 }
G4double SchwingerDeltaE
static const G4double Theta_Min
G4double QweakSimEPEvent::Elastic_Cross_Section_Proton ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 801 of file QweakSimEPEvent.cc.

References M_p.

Referenced by QweakSimEventAction::EndOfEventAction(), and SuperElasticCheck().

806 {
807  G4double Lamda_2 = 0.710;
808  G4double M_p = 938.2796 * MeV; // proton mass in MeV
809  G4double mu = 2.793;
810  G4double Z = 1.0;
811  G4double A = 1.0;
812  G4double M = M_p*A;
813  G4double myhbarc = hbarc / MeV / fermi; // 197.3269718 in MeV fm
814  G4double alpha = 1.0/137.035999074;
815  G4double CC = myhbarc*alpha/2.0; // 0.719982242379
816  G4double Electron_Mass = 0.511 * MeV; //electron mass in MeV
817 
818 // E_in units is MeV
819 
820  const G4double theta_min = 0.01 * degree;
821  if (Theta < theta_min) {
822  Theta = theta_min;
823  G4cout << "Warning: Elastic_Cross_Section_Proton: theta less than " << theta_min << G4endl;
824  G4cout << "Warning: Elastic_Cross_Section_Proton: theta was set to " << theta_min << G4endl;
825  }
826 
827  G4double CTH = cos(Theta/2.);
828  G4double STH = sin(Theta/2.);
829  G4double T2THE = STH*STH/CTH/CTH;
830  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
831  E_out = E_in/ETA;
832  Q2 = 4.0*E_in*E_out*STH*STH; // MeV^2
833  G4double tau = Q2/4.0/M/M;
834 // Mott scatering cross-section, including recoil correction
835  G4double CrossSection = (Z*CC/E_in*CTH/STH/STH)*(Z*CC/E_in*CTH/STH/STH)/ETA;
836 // Units: ub/sr
837  G4double Mott = CrossSection*10000.0;
838 // Cross section
839  G4double GEP_DIPOLE = 1.0/(1.0+Q2/1.E6/Lamda_2)/(1.0+Q2/1.E6/Lamda_2);
840  G4double GMP_DIPOLE = GEP_DIPOLE*mu;
841  G4double FAC = 1.0/(1.0+tau);
842 
843  //The next two line is to add schwinger
844  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
845  G4double delta_Schwinger = (-2.0*alpha/pi) * ((log(E_in/15.0) - 13.0/12.0) * (log(Q2/(Electron_Mass*Electron_Mass)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
846 // G4double omega_Sch = E_in - 15;
847  G4double Sigma_Dipole;
848  Sigma_Dipole = Mott*(GEP_DIPOLE*GEP_DIPOLE*FAC+tau*GMP_DIPOLE*GMP_DIPOLE*(FAC+2.*T2THE));
849  Sigma_Dipole *= (1.0 + delta_Schwinger);
850  fWeightN = Sigma_Dipole*sin(Theta);
851  return Sigma_Dipole;
852 }
static const G4double M_p

+ Here is the caller graph for this function:

void QweakSimEPEvent::F1F2IN09 ( G4int  Z,
G4int  A,
G4double  QSQ,
G4double  wsq,
G4double &  F1,
G4double &  F2 
)

Definition at line 2974 of file QweakSimEPEvent.cc.

References christy507(), fitemc(), MEC2009(), and resmodd().

Referenced by NuclearInelastic_Bosted().

2976 {
2977 /*--------------------------------------------------------------------
2978  Fit to inelastic cross sections for A(e,e')X
2979  valid for all W<3 GeV and all Q2<10 GeV2
2980 
2981  Inputs: Z, A (real*8) are Z and A of nucleus
2982  (use Z=0., A=1. to get free neutron)
2983  Qsq (real*8) is 4-vector momentum transfer squared (positive in
2984  chosen metric)
2985  Wsq (real*8) is invarinat mass squared of final state calculated
2986  assuming electron scattered from a free proton
2987 
2988  outputs: F1, F2 (real*8) are structure functions per nucleus
2989  Version of 10/20/2006 P. Bosted
2990 --------------------------------------------------------------------*/
2991 
2992  G4double W, x;
2993  G4double nu,siglp,sigtp,F1pp,F1dp;
2994  G4double W1,W2,sigt,Rc,sigl,F1d, F1p,qv;
2995  G4double DW2DPF,Wsqp,pf,kf,Es,dw2des,Fyuse;
2996  G4double x4, emcfac;
2997  G4int ism;
2998  G4double PM = 0.93828;
2999 
3000  // This is for exp(-xx**2/2.), from teste.f
3001  G4double XXp[99] = {
3002  -3.000,-2.939,-2.878,-2.816,-2.755,-2.694,-2.633,-2.571,-2.510,
3003  -2.449,-2.388,-2.327,-2.265,-2.204,-2.143,-2.082,-2.020,-1.959,
3004  -1.898,-1.837,-1.776,-1.714,-1.653,-1.592,-1.531,-1.469,-1.408,
3005  -1.347,-1.286,-1.224,-1.163,-1.102,-1.041,-0.980,-0.918,-0.857,
3006  -0.796,-0.735,-0.673,-0.612,-0.551,-0.490,-0.429,-0.367,-0.306,
3007  -0.245,-0.184,-0.122,-0.061, 0.000, 0.061, 0.122, 0.184, 0.245,
3008  0.306, 0.367, 0.429, 0.490, 0.551, 0.612, 0.673, 0.735, 0.796,
3009  0.857, 0.918, 0.980, 1.041, 1.102, 1.163, 1.224, 1.286, 1.347,
3010  1.408, 1.469, 1.531, 1.592, 1.653, 1.714, 1.776, 1.837, 1.898,
3011  1.959, 2.020, 2.082, 2.143, 2.204, 2.265, 2.327, 2.388, 2.449,
3012  2.510, 2.571, 2.633, 2.694, 2.755, 2.816, 2.878, 2.939, 3.000};
3013  // these are 100x bigger for convenience
3014  G4double fyp[99] = {
3015  0.0272,0.0326,0.0390,0.0464,0.0551,0.0651,0.0766,0.0898,0.1049,
3016  0.1221,0.1416,0.1636,0.1883,0.2159,0.2466,0.2807,0.3182,0.3595,
3017  0.4045,0.4535,0.5066,0.5637,0.6249,0.6901,0.7593,0.8324,0.9090,
3018  0.9890,1.0720,1.1577,1.2454,1.3349,1.4254,1.5163,1.6070,1.6968,
3019  1.7849,1.8705,1.9529,2.0313,2.1049,2.1731,2.2350,2.2901,2.3379,
3020  2.3776,2.4090,2.4317,2.4454,2.4500,2.4454,2.4317,2.4090,2.3776,
3021  2.3379,2.2901,2.2350,2.1731,2.1049,2.0313,1.9529,1.8705,1.7849,
3022  1.6968,1.6070,1.5163,1.4254,1.3349,1.2454,1.1577,1.0720,0.9890,
3023  0.9090,0.8324,0.7593,0.6901,0.6249,0.5637,0.5066,0.4535,0.4045,
3024  0.3595,0.3182,0.2807,0.2466,0.2159,0.1883,0.1636,0.1416,0.1221,
3025  0.1049,0.0898,0.0766,0.0651,0.0551,0.0464,0.0390,0.0326,0.0272};
3026 
3027  // deuteron fit parameters
3028  G4double xvald0[50] = {
3029  0.1964E+01, 0.1086E+01, 0.5313E-02, 0.1265E+01, 0.8000E+01,
3030  0.2979E+00, 0.1354E+00, 0.2200E+00, 0.8296E-01, 0.9578E-01,
3031  0.1094E+00, 0.3794E+00, 0.8122E+01, 0.5189E+01, 0.3290E+01,
3032  0.1870E+01, 0.6110E+01,-0.3464E+02, 0.9000E+03, 0.1717E+01,
3033  0.4335E-01, 0.1915E+03, 0.2232E+00, 0.2119E+01, 0.2088E+01,
3034  -0.3029E+00, 0.2012E+00, 0.1104E-02, 0.2276E-01,-0.4562E+00,
3035  0.2397E+00, 0.1204E+01, 0.2321E-01, 0.5419E+03, 0.2247E+00,
3036  0.2168E+01, 0.2266E+03, 0.7649E-01, 0.1457E+01, 0.1318E+00,
3037  -0.7534E+02, 0.1776E+00, 0.1636E+01, 0.1350E+00,-0.5596E-02,
3038  0.5883E-02, 0.1934E+01, 0.3800E+00, 0.3319E+01, 0.1446E+00};
3039 
3040  G4double P[24] = {
3041  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
3042  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
3043  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
3044  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
3045  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
3046  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
3047 
3048  G4double F1M;
3049  nu = (wsq - PM*PM + qsq) / 2. / PM;
3050  qv = sqrt(nu*nu + qsq);
3051 
3052  if (wsq <= 0.0) {
3053  W = 0.0;
3054  } else {
3055  W = sqrt(wsq);
3056  x = qsq / (2.0 * PM * nu);
3057  }
3058  if (wsq <= 0.0) x = 0.0;
3059 
3060  if (IA <= 2) return;
3061 
3062  // For nuclei
3063  if (IA > 2) {
3064  sigt = 0.;
3065  sigl = 0.;
3066  F1d = 0.;
3067  F1p = 0.;
3068  }
3069 
3070  // Modifed to use Superscaling from Sick, Donnelly, Maieron,
3071  // nucl-th/0109032
3072  if(IA==2) kf=0.085;
3073  if(IA==2) Es=0.0022;
3074  // changed 4/09
3075  if(IA==3) kf=0.115;
3076  if(IA==3) Es=0.001 ;
3077  // changed 4/09
3078  if(IA>3) kf=0.19;
3079  if(IA>3) Es=0.017;
3080  if(IA>7) kf=0.228;
3081  if(IA>7) Es=0.020;
3082  // changed 5/09
3083  if(IA>7) Es=0.0165;
3084  if(IA>16) kf=0.230;
3085  if(IA>16) Es=0.025;
3086  if(IA>25) kf=0.236;
3087  if(IA>25) Es=0.018;
3088  if(IA>38) kf=0.241;
3089  if(IA>38) Es=0.028;
3090  if(IA>55) kf=0.241;
3091  if(IA>55) Es=0.023;
3092  if(IA>60) kf=0.245;
3093  if(IA>60) Es=0.028;
3094  // changed 5/09
3095  if(IA>55) Es=0.018;
3096 
3097  // adjust pf to give right width based on kf
3098  pf = 0.5 * kf;
3099  // assume this is 2 * pf * qv
3100  DW2DPF = 2. * qv;
3101  dw2des = 2. * (nu + PM) ;
3102  // switched to using 99 bins!
3103  for (ism = 0; ism < 99; ism++) {
3104  Fyuse = fyp[ism]/100.;
3105  Wsqp = wsq + XXp[ism] * pf * DW2DPF - Es * dw2des;
3106  if (Wsqp > 1.159) {
3107  christy507(Wsqp,qsq,F1pp,Rc,sigtp,siglp);
3108  resmodd(Wsqp,qsq,xvald0,F1dp);
3109  F1d = F1d + F1dp * Fyuse;
3110  F1p = F1p + F1pp * Fyuse;
3111  sigt = sigt + sigtp * Fyuse;
3112  sigl = sigl + siglp * Fyuse;
3113  }
3114  }
3115 
3116  Rc = 0.;
3117  if(sigt > 0.) Rc = sigl / sigt;
3118  W1 = (2. * Z * F1d + (IA - 2. * Z) * (2. * F1d - F1p)) / PM;
3119  W1= W1*(1.0+P[13]*x+P[14]*pow(x,2)+P[15]*pow(x,3)+P[16]*pow(x,4)+P[17]*pow(x,5));
3120 
3121  if(W > 0.0) W1=W1*pow((1.0+(P[20]*W+P[21]*pow(W,2))/(1.0+P[22]*qsq)),2);
3122  F1M = MEC2009(IA,qsq,wsq);
3123  W1 = W1 + F1M;
3124  if(wsq > 0.0 ) Rc = Rc * ( 1.0 + P[6] + P[23]*IA );
3125  W2 = W1 * (1. + Rc) / (1. + nu*nu / qsq);
3126 
3127  x4 = qsq / 2. / PM / nu;
3128  emcfac = fitemc(x4, IA);
3129  F1 = PM * W1 * emcfac;
3130  F2 = nu * W2 * emcfac;
3131 
3132  return;
3133 
3134 }
G4double fitemc(G4double X, G4int A)
void christy507(G4double wsq, G4double Q2, G4double &F1, G4double &R, G4double &sigt, G4double &sigl)
G4double MEC2009(G4int a, G4double q2, G4double w2)
void resmodd(G4double w2, G4double q2, G4double xval[50], G4double &sig)

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::F1F2QE09 ( G4int  Z,
G4int  A,
G4double  QSQ,
G4double  wsq,
G4double &  F1,
G4double &  F2 
)

for avp2 term

Definition at line 2659 of file QweakSimEPEvent.cc.

Referenced by Quasi_Elastic_Bosted().

2661 {
2662 //=======================================================================
2663 // SUBROUTINE F1F2QE09(Z, A, QSQ, wsq, F1, F2)
2664 //
2665 // Calculates quasielastic A(e,e')X structure functions F1 and F2 PER NUCLEUS
2666 // for A>2 uses superscaling from Sick, Donnelly, Maieron, nucl-th/0109032
2667 // for A=2 uses pre-integrated Paris wave function (see ~bosted/smear.f)
2668 // coded by P. Bosted August to October, 2006
2669 //
2670 // input: Z, A (real*8) Z and A of nucleus (shoud be 2.0D0 for deueron)
2671 // Qsq (real*8) is 4-vector momentum transfer squared (positive in
2672 // chosen metric)
2673 // Wsq (real*8) is invarinat mass squared of final state calculated
2674 // assuming electron scattered from a free proton
2675 //
2676 // outputs: F1, F2 (real*8) are structure functions per nucleus
2677 //
2678 // Note: Deuteron agrees well with Laget (see ~bosted/eg1b/laget.f) for
2679 // a) Q2<1 gev**2 and dsig > 1% of peak: doesnt describe tail at high W
2680 // b) Q2>1 gev**2 on wings of q.e. peak. But, this model is up
2681 // to 50% too big at top of q.e. peak. BUT, F2 DOES agree very
2682 // nicely with Osipenko et al data from CLAS, up to 5 GeV**2
2683 
2684 // G4cout << "Z, A, Q2, W2: " << Z << ", " << IA << ", " << QSQ << ", " << wsq << G4endl;
2685 
2686  G4double avgN, Pauli_sup1, Pauli_sup2, GEP, GEN, GMP, GMN, Q, Q3, Q4;
2687  G4double amp = 0.93828;
2688  G4double amd = 1.8756;
2689  G4double RMUP = 2.792782;
2690  G4double RMUN = -1.913148;
2691  G4double pz, Nu, dpz, pznom, pzmin;
2692  G4double QV, TAU, FY, dwmin, w2p;
2693  G4double kappa, lam, lamp, taup, squigglef, psi, psip, nuL, nuT;
2694  G4double kf, Es, GM2bar, GE2bar, W2bar, Delta, GL, GT;
2695  G4int izz, izzmin, izp, izznom, izdif;
2696 
2697  // Look up tables for deuteron case
2698  G4double fyd[200] = {
2699  0.00001,0.00002,0.00003,0.00005,0.00006,0.00009,0.00010,0.00013,
2700  0.00015,0.00019,0.00021,0.00026,0.00029,0.00034,0.00038,0.00044,
2701  0.00049,0.00057,0.00062,0.00071,0.00078,0.00089,0.00097,0.00109,
2702  0.00119,0.00134,0.00146,0.00161,0.00176,0.00195,0.00211,0.00232,
2703  0.00252,0.00276,0.00299,0.00326,0.00352,0.00383,0.00412,0.00447,
2704  0.00482,0.00521,0.00560,0.00603,0.00648,0.00698,0.00747,0.00803,
2705  0.00859,0.00921,0.00985,0.01056,0.01126,0.01205,0.01286,0.01376,
2706  0.01467,0.01569,0.01671,0.01793,0.01912,0.02049,0.02196,0.02356,
2707  0.02525,0.02723,0.02939,0.03179,0.03453,0.03764,0.04116,0.04533,
2708  0.05004,0.05565,0.06232,0.07015,0.07965,0.09093,0.10486,0.12185,
2709  0.14268,0.16860,0.20074,0.24129,0.29201,0.35713,0.44012,0.54757,
2710  0.68665,0.86965,1.11199,1.43242,1.86532,2.44703,3.22681,4.24972,
2711  5.54382,7.04016,8.48123,9.40627,9.40627,8.48123,7.04016,5.54382,
2712  4.24972,3.22681,2.44703,1.86532,1.43242,1.11199,0.86965,0.68665,
2713  0.54757,0.44012,0.35713,0.29201,0.24129,0.20074,0.16860,0.14268,
2714  0.12185,0.10486,0.09093,0.07965,0.07015,0.06232,0.05565,0.05004,
2715  0.04533,0.04116,0.03764,0.03453,0.03179,0.02939,0.02723,0.02525,
2716  0.02356,0.02196,0.02049,0.01912,0.01793,0.01671,0.01569,0.01467,
2717  0.01376,0.01286,0.01205,0.01126,0.01056,0.00985,0.00921,0.00859,
2718  0.00803,0.00747,0.00698,0.00648,0.00603,0.00560,0.00521,0.00482,
2719  0.00447,0.00412,0.00383,0.00352,0.00326,0.00299,0.00276,0.00252,
2720  0.00232,0.00211,0.00195,0.00176,0.00161,0.00146,0.00134,0.00119,
2721  0.00109,0.00097,0.00089,0.00078,0.00071,0.00062,0.00057,0.00049,
2722  0.00044,0.00038,0.00034,0.00029,0.00026,0.00021,0.00019,0.00015,
2723  0.00013,0.00010,0.00009,0.00006,0.00005,0.00003,0.00002,0.00001};
2724 
2725  G4double avp2[200] = {
2726  1.0,0.98974,0.96975,0.96768,0.94782,0.94450,0.92494,0.92047,
2727  0.90090,0.89563,0.87644,0.87018,0.85145,0.84434,0.82593,0.81841,
2728  0.80021,0.79212,0.77444,0.76553,0.74866,0.73945,0.72264,0.71343,
2729  0.69703,0.68740,0.67149,0.66182,0.64631,0.63630,0.62125,0.61154,
2730  0.59671,0.58686,0.57241,0.56283,0.54866,0.53889,0.52528,0.51581,
2731  0.50236,0.49291,0.47997,0.47063,0.45803,0.44867,0.43665,0.42744,
2732  0.41554,0.40656,0.39511,0.38589,0.37488,0.36611,0.35516,0.34647,
2733  0.33571,0.32704,0.31656,0.30783,0.29741,0.28870,0.27820,0.26945,
2734  0.25898,0.25010,0.23945,0.23023,0.21943,0.20999,0.19891,0.18911,
2735  0.17795,0.16793,0.15669,0.14667,0.13553,0.12569,0.11504,0.10550,
2736  0.09557,0.08674,0.07774,0.06974,0.06184,0.05484,0.04802,0.04203,
2737  0.03629,0.03129,0.02654,0.02247,0.01867,0.01545,0.01251,0.01015,
2738  0.00810,0.00664,0.00541,0.00512,0.00512,0.00541,0.00664,0.00810,
2739  0.01015,0.01251,0.01545,0.01867,0.02247,0.02654,0.03129,0.03629,
2740  0.04203,0.04802,0.05484,0.06184,0.06974,0.07774,0.08674,0.09557,
2741  0.10550,0.11504,0.12569,0.13553,0.14667,0.15669,0.16793,0.17795,
2742  0.18911,0.19891,0.20999,0.21943,0.23023,0.23945,0.25010,0.25898,
2743  0.26945,0.27820,0.28870,0.29741,0.30783,0.31656,0.32704,0.33571,
2744  0.34647,0.35516,0.36611,0.37488,0.38589,0.39511,0.40656,0.41554,
2745  0.42744,0.43665,0.44867,0.45803,0.47063,0.47997,0.49291,0.50236,
2746  0.51581,0.52528,0.53889,0.54866,0.56283,0.57241,0.58686,0.59671,
2747  0.61154,0.62125,0.63630,0.64631,0.66182,0.67149,0.68740,0.69703,
2748  0.71343,0.72264,0.73945,0.74866,0.76553,0.77444,0.79212,0.80021,
2749  0.81841,0.82593,0.84434,0.85145,0.87018,0.87644,0.89563,0.90090,
2750  0.92047,0.92494,0.94450,0.94782,0.96768,0.96975,0.98974,1.0};
2751 
2752  // Peter Bosted's correction params
2753  /*
2754  G4double pb[20] = {
2755  0.1023E+02, 0.1052E+01, 0.2485E-01, 0.1455E+01,
2756  0.5650E+01,-0.2889E+00, 0.4943E-01,-0.8183E-01,
2757  -0.7495E+00, 0.8426E+00,-0.2829E+01, 0.1607E+01,
2758  0.1733E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00,
2759  0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00};
2760  */
2761  G4double y,R;
2762 
2763  G4double P[24] = {
2764  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
2765  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
2766  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
2767  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
2768  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
2769  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
2770 
2771  // return if proton: future change this to allow for
2772  // equivalent W resolution
2773  F1 = 0.;
2774  F2 = 0.;
2775  avgN = IA - Z;
2776  if (IA==1) return;
2777 
2778  // some kinematic factors. Return if Nu or QSQ is negative
2779  Nu = (wsq - amp*amp + QSQ) / 2. / amp;
2780 
2781  //G4cout << "In call... IA, Nu, QSQ = " << IA << ", " << Nu << ", " << QSQ << G4endl;
2782  if(Nu <= 0.0 || QSQ < 0.) return;
2783  TAU = QSQ / 4.0 / amp / amp;
2784  QV = sqrt(Nu*Nu + QSQ);
2785 
2786  // Bosted fit for nucleon form factors Phys. Rev. C 51, p. 409 (1995)
2787  Q = sqrt(QSQ);
2788  Q3 = QSQ * Q;
2789  Q4 = QSQ*QSQ;
2790  GEP = 1./ (1. + 0.14 * Q + 3.01 * QSQ + 0.02 * Q3 + 1.20 * Q4 + 0.32 * pow(Q,5));
2791  GMP = RMUP * GEP;
2792  GMN = RMUN / (1.- 1.74 * Q + 9.29 * QSQ - 7.63 * Q3 + 4.63 * Q4);
2793  GEN = 1.25 * RMUN * TAU / (1. + 18.3 * TAU) / pow((1. + QSQ / 0.71),2);
2794 
2795  //G4cout << "Form Factors: " << GEP << ", " << GMP << ", " << GEN << ", " << GMN << G4endl;
2796 
2797  // Get kf and Es from superscaling from Sick, Donnelly, Maieron,
2798  // nucl-th/0109032
2799  if(IA==2) kf=0.085;
2800  if(IA==2) Es=0.0022;
2801  // changed 4/09
2802  if(IA==3) kf=0.115;
2803  if(IA==3) Es=0.001 ;
2804  // changed 4/09
2805  if(IA>3) kf=0.19;
2806  if(IA>3) Es=0.017;
2807  if(IA>7) kf=0.228;
2808  if(IA>7) Es=0.020;
2809  // changed 5/09
2810  if(IA>7) Es=0.0165;
2811  if(IA>16) kf=0.230;
2812  if(IA>16) Es=0.025;
2813  if(IA>25) kf=0.236;
2814  if(IA>25) Es=0.018;
2815  if(IA>38) kf=0.241;
2816  if(IA>38) Es=0.028;
2817  if(IA>55) kf=0.241;
2818  if(IA>55) Es=0.023;
2819  if(IA>60) kf=0.245;
2820  if(IA>60) Es=0.028;
2821  // changed 5/09
2822  if(IA>55) Es=0.018;
2823 
2824  // Pauli suppression model from Tsai RMP 46,816(74) eq.B54
2825  if ((QV > 2.* kf) || (IA == 1)) {
2826  Pauli_sup2 =1.0;
2827  } else {
2828  Pauli_sup2 = 0.75 * (QV / kf) * (1.0 - (pow((QV / kf),2))/12.);
2829  }
2830  Pauli_sup1 = Pauli_sup2;
2831 
2832  //G4cout << "kf, Es, Paulisup1,2: " << kf << ", " << Es << ", " << Pauli_sup1 << ", " << Pauli_sup2 << G4endl;
2833 
2834  // structure functions with off shell factors
2835  kappa = QV / 2. / amp;
2836  lam = Nu / 2. / amp;
2837  lamp = lam - Es / 2. / amp;
2838  taup = kappa*kappa - lamp*lamp;
2839  squigglef = sqrt(1. + pow((kf/amp),2)) -1.;
2840 
2841  // Very close to treshold, could have a problem
2842  if(1.+lamp <= 0.) return;
2843  if(taup * (1. + taup) <= 0.) return;
2844 
2845  psi = (lam - TAU ) / sqrt(squigglef) / sqrt((1.+lam )* TAU + kappa * sqrt(TAU * (1. + TAU)));
2846  psip = (lamp - taup) / sqrt(squigglef) / sqrt((1.+lamp)*taup + kappa * sqrt(taup * (1. + taup)));
2847  nuL = pow((TAU / kappa / kappa),2);
2848 
2849  // changed definition of nuT from
2850  // nuT = TAU / 2. / kappa**2 + tan(thr/2.)**2
2851  // to this, in order to separate out F1 and F2 (F1 prop. to tan2 term)
2852  nuT = TAU / 2. / kappa / kappa;
2853 
2854  GM2bar = Pauli_sup1 * (Z * GMP*GMP + avgN * GMN*GMN);
2855  GE2bar = Pauli_sup2 * (Z * GEP*GEP + avgN * GEN*GEN);
2856  //G4double W1bar = TAU * GM2bar;
2857  W2bar = (GE2bar + TAU * GM2bar) / (1. + TAU);
2858 
2859  Delta = squigglef * (1. - psi*psi) * (sqrt(TAU * (1.+TAU)) / kappa + squigglef/3. *
2860  (1. - psi*psi) * TAU / kappa / kappa);
2861  GL = kappa*kappa / TAU * (GE2bar + Delta * W2bar) / 2. / kappa / (1. + squigglef *
2862  (1. + psi*psi) / 2.);
2863  GT = (2. * TAU * GM2bar + Delta * W2bar) / 2. / kappa / (1. + squigglef *
2864  (1. + psi*psi) / 2.);
2865 
2866  //G4cout << "nuL, nuT, GL, GT: " << nuL << ", " << nuT << ", " << GL << ", " << GT << G4endl;
2867 
2868  // added to prevent negative xsections:
2869  if (GT < 0) {
2870  GT = 0;
2871  //G4cout << "Reset GT to zero" << G4endl;
2872  }
2873 
2874  // from Maria Barbaro: see Amaro et al., PRC71,015501(2005).
2875  FY = 1.5576 / (1. + 1.7720*1.7720 * pow((psip + 0.3014),2)) / (1. + exp(-2.4291 * psip)) / kf;
2876 
2877 // == Note: Below is for deuteron and haven't == //
2878 // == verified if code conversion is correct == //
2879  // Use PWIA and Paris W.F. for deuteron to get better FY
2880  if (IA == 2) {
2881  // value assuming average p2=0.
2882  pz = (QSQ - 2. * amp * Nu ) / 2. / QV;
2883  izz = int((pz + 1.0) / 0.01) + 1;
2884  if (izz < 1) izz = 1;
2885  if (izz > 200) izz = 200;
2886  izznom = izz;
2887 
2888  // ignoring energy term, estimate change in pz to compensate
2889  //! for avp2 term
2890  dpz = avp2[izznom] / 2. / QV;
2891  izdif = dpz * 150.;
2892  dwmin=1.E6;
2893  izzmin=0;
2894 
2895  G4int izpmax;
2896  if ((izznom + izdif) < 1) izpmax = 1;
2897  if (izpmax > 200) izpmax = 200;
2898  for (izp = izznom; izp <= izpmax; izp++) {
2899  pz = -1. + 0.01 * (izp-0.5);
2900  // *** this version gives worse agreement with laget than
2901  // w2p = (amd + Nu - sqrt(amp**2 + avp2(izp)))**2 -
2902  // > QV**2 + 2. * QV * pz - avp2(izp)
2903  // this version!
2904  w2p = pow((amd + Nu - amp ),2) - QV*QV + 2. * QV * pz - avp2[izp];
2905 
2906  // if passed first minimum, quit looking so don't find second one
2907  if (abs(w2p - amp*amp) > dwmin) {
2908  if (izzmin < 2) izzmin = 2;
2909  if (izzmin > 199) izzmin = 199;
2910  izz = izzmin;
2911  } else if (abs(w2p - amp*amp) < dwmin) {
2912  dwmin = abs(w2p - amp*amp);
2913  izzmin = izp;
2914  }
2915  }
2916 
2917  // search for minimum in 1/10th bins locally
2918  pznom = -1. + 0.01 * (izz-0.5);
2919  dwmin=1.E6;
2920  for (izp = 1; izp <= 19; izp++) {
2921  pz = pznom - 0.01 + 0.001 * izp;
2922  // *** this version gives worse agreement with laget than
2923  // w2p = (amd + Nu - sqrt(amp**2 + avp2(izz)))**2 -
2924  // > QV**2 + 2. * QV * pz - avp2(izz)
2925  // this version!
2926  w2p = pow((amd + Nu - amp ),2) - QV*QV + 2. * QV * pz - avp2[izz];
2927  if (abs(w2p - amp*amp) < dwmin) {
2928  dwmin = abs(w2p - amp*amp);
2929  pzmin = pz;
2930  }
2931  }
2932 
2933  if (dwmin >= 1.e6 || abs(pznom-pzmin) > 0.01) {
2934  // > write(6,'(1x,''error in dwmin,pzmin'',3i4,6f7.3)')
2935  // > izznom,izzmin,izz,QSQ,wsq,w2p,dwmin/1.e6,pzmin,pznom
2936  if (pzmin < pznom) FY = fyd[izz] - (fyd[izz-1] - fyd[izz]) * (pzmin - pznom) / 0.01;
2937  else FY = fyd[izz] + (fyd[izz+1] - fyd[izz]) * (pzmin - pznom) / 0.01;
2938  }
2939  }
2940 // == To here == //
2941 
2942  // final results
2943  F2 = Nu * FY * (nuL * GL + nuT * GT);
2944  F1 = amp * FY * GT / 2.;
2945 
2946  //G4cout << "nu, Fy, nuL, GL, nuT, GT, amp: " << G4endl;
2947  //G4cout << Nu << ", " << FY << ", " << nuL << ", " << GL << ", " << nuT << ", " << GT << ", " << amp << G4endl;
2948 
2949  if (F1 < 0.0) F1 = 0.;
2950  if (Nu > 0. && F1 > 0.) R = (F2 / Nu) / (F1 / amp) * (1. + Nu*Nu / QSQ) - 1.0;
2951  else R = 0.4/QSQ;
2952 
2953 
2954  // apply correction factors
2955  if ( IA > 2 ) {
2956  y = (wsq -amp*amp) / QV;
2957  // F1 = F1 * (1. + pb(8) + pb(9) * y +
2958  // > pb(10)*y**2 + pb(11)*y**3 + pb(12)*y**4 )
2959  // R = R * (1. + pb(13))
2960  // F2 = Nu * F1/amp * (1. + R) / (1. + Nu**2/QSQ)
2961 
2962  // correction to correction Vahe
2963  if (wsq > 0.0) {
2964  F1=F1*(1.0+P[7]+P[8]*y+P[9]*y*y +P[10]*pow(y,3) +P[11]*pow(y,4));
2965  R = R * ( 1.0 + P[12] );
2966  F2 = Nu * F1/amp * (1. + R) / (1. + Nu*Nu/QSQ);
2967  if (F1 < 0.0) F1=0.0;
2968  }
2969  }
2970 
2971  return;
2972 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::FF_BESSEL ( G4int  Z,
G4int  A,
G4double  q2,
G4bool &  ofr 
)

line added 4/29/98

4*PI

2*PI

Definition at line 3835 of file QweakSimEPEvent.cc.

Referenced by AlloyScattering().

3836  {
3837 
3838  // Calculates PWBA form factor for various nuclei using
3839  // Fourier-Bessel coefficients. (Inverse Fourier transform)
3840  //
3841  // Sources:
3842  // Technique-
3843  // 1)B. Dreher et al., Nuc. Phys. A235 219-248 (1974)
3844  // Data-
3845  // 1)H. de Vries et al., Atomic Data and Nuclear Data Tables 36,
3846  // 495-536 (1987) - Table IV pg. 520
3847  // 2)G. Fricke et al., Atomic Data and Nuclear Data Tables 60,
3848  // 177-285 (1995) - Table IX pg. 264
3849  //
3850  // Note: Can be used for nuclei with coeffients listed below
3851  // for a limited Q^2 range.
3852  //----------------------------------------------------------------
3853  // 2003/ 3/28 - Corrected for divide 0/0 when qm =0 - Steve Rock
3854  // 2016/10/25 - Updated references, included addition FB coeffients,
3855  // for more nuclei and isotopes of existing nuclei. - K. Bartlett
3856  // 2016/10/26 - Included a few more FB coefficients for various
3857  // isotopes.
3858 
3859  G4double a[20];
3860 
3861  //3He2 Helium-3
3862  G4double a_3He[] = {0.20020E-01, 0.41934E-01, 0.36254E-01, 0.17941E-01,
3863  0.46608E-02, 0.46834E-02, 0.52042E-02, 0.38280E-02,
3864  0.25661E-02, 0.14182E-02, 0.61390E-03, 0.22929E-03};
3865  //12C6 Carbon-12
3866  G4double a_12C[] = {0.15721E-01, 0.38732E-01, 0.36808E-01, 0.14671E-01,
3867  -0.43277E-02,-0.97752E-02,-0.68908E-02,-0.27631E-02,
3868  -0.63568E-03, 0.71809E-04, 0.18441E-03, 0.75066E-04,
3869  0.51069E-04, 0.14308E-04, 0.23170E-05, 0.68465E-06};
3870  //15N7 Nitrogen-15
3871  G4double a_15N[] = {0.25491E-01, 0.50618E-01, 0.29822E-01,-0.55196E-02,
3872  -0.15913E-01,-0.76184E-02,-0.23992E-02,-0.47940E-03};
3873 
3874  //16O8 Oxygen-16
3875  G4double a_16O[] = {0.20238E-01, 0.44793E-01, 0.33533E-01, 0.35030E-02,
3876  -0.12293E-01,-0.10329E-01,-0.34036E-02,-0.41627E-03,
3877  -0.94435E-03,-0.25571E-03, 0.23759E-03,-0.10603E-03,
3878  0.41480E-04, 0.00000E-03};
3879  //26Mg12 Magnesium-26
3880  G4double a_26Mg[] = {2.9280e-2, 5.5520e-2, 2.3910e-2,-1.3940e-2,-1.5840e-2,
3881  -3.5680e-3, 1.9670e-3, 1.4470e-3,-8.6820e-5,-1.4000e-4,
3882  7.0000e-5,-2.0000e-5, 3.0000e-6, 1.0000e-7,-2.0000e-7};
3883 
3884  //27Al13 Aluminum-27
3885  G4double a_27Al[] = {0.43418E-01, 0.60298E-01, 0.28950E-02,-0.23522E-01,
3886  -0.79791E-02, 0.23010E-02, 0.10794E-02, 0.12574E-03,
3887  -0.13021E-03, 0.56563E-04,-0.18011E-04, 0.42869E-05};
3888 
3889  //28Si14 Silicon-28
3890  G4double a_28Si[] = {0.33495E-01, 0.59533E-01, 0.20979E-01,-0.16900E-01,
3891  -0.14998E-01,-0.93248E-03, 0.33266E-02, 0.59244E-03,
3892  -0.40013E-03, 0.12242E-03,-0.12994E-04,-0.92784E-05,
3893  0.72595E-05,-0.42096E-05};
3894  //29Si14 Silicon-29
3895  G4double a_29Si[] = {3.3521e-02, 5.9679e-02, 2.0593e-02,-1.8646e-02,
3896  -1.6550e-02,-1.1922e-03, 2.8025e-03,-6.7353e-05,
3897  -3.4619e-04, 1.7611e-04,-5.7173e-06, 1.2371e-05};
3898  //30Si14 Silicon-30
3899  G4double a_30Si[] = {2.8397e-02, 5.4163e-02, 2.5167e-02,-1.2858e-02,
3900  -1.7592e-02,-4.6722e-03, 2.4804e-03, 1.4760e-03,
3901  -3.0168e-04, 4.8346e-05, 0.0000e+00,-5.1570e-06,
3902  3.0261e-06};
3903  //48Ti22 Titanium-48
3904  G4double a_48Ti[] = {2.7850e-02, 5.5432e-02, 2.6369e-02,-1.7091e-02,
3905  -2.1798e-02,-2.4889e-03, 7.6631e-03, 3.4554e-03,
3906  -6.7477e-04, 1.0764e-04,-1.6564e-06,-5.5566e-06};
3907  //50Ti22 Titanium-50
3908  G4double a_50Ti[] = {3.1818e-02, 5.8556e-02, 1.9637e-02,-2.4309e-02,
3909  -1.8748e-02, 3.3741e-03, 8.9961e-03, 3.7954e-03,
3910  -4.1238e-04, 1.2540e-04};
3911  //50Cr24 Chromium-50
3912  G4double a_50Cr[] = {3.9174e-02, 6.1822e-02, 6.8550e-03,-3.0170e-02,
3913  -9.8745e-03, 8.7944e-03, 6.8502e-03,-9.3609e-04,
3914  -2.4962e-03,-1.5361e-03,-7.3687e-04};
3915  //52Cr24 Chromium-52
3916  G4double a_52Cr[] = {0.39287e-1, 0.62477e-1, 0.62482e-2,-0.32885e-1,
3917  -0.10648e-1, 0.10520e-1, 0.85478e-2,-0.24003e-3,
3918  -0.20499e-2,-0.12001e-2,-0.56649e-3};
3919  //54Cr24 Chromium-54
3920  G4double a_54Cr[] = {3.9002e-02, 6.0305e-02, 4.5845e-03,-3.0723e-02,
3921  -9.1355e-03, 9.3251e-03, 6.0583e-03,-1.5602e-03,
3922  -7.6809e-04, 7.6809e-04,-3.4804e-04};
3923  //54Fe26 Iron-54
3924  G4double a_54Fe[] = {4.2339e-02, 6.4428e-02, 1.5840e-03,-3.5171e-02,
3925  -1.0116e-02, 1.2069e-02, 6.2230e-03,-1.2045e-03,
3926  1.3561e-04, 1.0428e-05,-1.6980e-05, 9.1817e-06,
3927  -3.9988e-06, 1.5731e-06,-5.7862e-07, 2.0186e-07,
3928  -6.7892e-08};
3929  //56Fe26 Iron-56
3930  G4double a_56Fe[] = {0.42018E-1, 0.62337E-1, 0.23995e-3,-0.32776E-1,
3931  -0.79941E-2, 0.10844E-1, 0.49123e-2,-0.22144e-2,
3932  -0.18146E-3, 0.37261E-3,-0.23296E-3, 0.11494E-3,
3933  -0.50596E-4, 0.20652E-4,-0.79428E-5, 0.28986E-5,
3934  -0.10075E-5};
3935  //58Fe26 Iron-58
3936  G4double a_58Fe[] = {4.1791e-02, 6.0524e-02,-1.4978e-03,-3.1183e-02,
3937  -5.8013e-03, 1.0611e-02, 4.1629e-03,-2.9045e-06,
3938  5.4106e-04,-3.8689e-04, 2.0514e-04,-9.5237e-05,
3939  4.0707e-05,-1.6346e-05, 6.2233e-06,-2.2568e-06,
3940  7.8077e-07};
3941  //63Cu29 Copper-63
3942  G4double a_63Cu[] = {4.5598e-02, 6.0706e-02,-7.8616e-03,-3.1638e-02,
3943  -1.4447e-03, 1.0953e-02, 4.2578e-03,-2.4224e-04,
3944  -3.0067e-04, 2.3903e-04,-1.2910e-04, 6.0195e-05,
3945  -2.5755e-05, 1.0332e-05,-3.9330e-06, 1.4254e-06,
3946  -4.9221e-07};
3947  //65Cu29 Copper-65
3948  G4double a_65Cu[] = {0.45444E-01, 0.59544E-01,-0.94968E-02,-0.31561e-01,
3949  0.22898E-03, 0.11189E-01, 0.37360E-02,-0.64873E-03,
3950  -0.51133E-03, 0.43765E-03,-0.24276E-03, 0.11507E-03,
3951  -0.49761E-04, 0.20140E-04,-0.76945E-05, 0.28055E-05,
3952  -0.97411E-06};
3953  //64Zn30 Zinc-64
3954  G4double a_64Zn[] = {0.47038E-01, 0.61536E-01,-0.90045E-02,-0.30669E-01,
3955  -0.78705E-03, 0.10034E-01, 0.14053E-02,-0.20640E-02,
3956  0.35105E-03, 0.27303E-04,-0.63811E-04, 0.40893E-04,
3957  -0.20311E-04, 0.88986E-05,-0.35849E-05, 0.13522E-05,
3958  -0.38635E-06};
3959  //66Zn30 Zinc-66
3960  G4double a_66Zn[] = {4.6991e-02, 6.0995e-02,-9.6693e-03,-3.0457e-02,
3961  -5.3435e-04, 9.7083e-03, 1.4091e-03,-7.0813e-04,
3962  2.0809e-04,-4.8275e-05, 7.2680e-06, 9.1369e-07,
3963  -1.4874e-06, 8.8831e-07,-4.1689e-07, 1.7283e-07,
3964  -6.5968e-08};
3965  //68Zn30 Zinc-68
3966  G4double a_68Zn[] = {4.6654e-02, 5.8827e-02,-1.2283e-02,-2.9865e-02,
3967  2.5669e-03, 1.0235e-02, 3.1861e-03,-1.7351e-04,
3968  -4.2979e-04, 3.3700e-04,-1.8435e-04, 8.7043e-05,
3969  -3.7612e-05, 1.5220e-05,-5.8282e-06, 2.1230e-06,
3970  -7.3709e-07};
3971  //70Zn30 Zinc-70
3972  G4double a_70Zn[] = {4.6362e-02, 5.7130e-02,-1.3877e-02,-3.0030e-02,
3973  3.5341e-03, 1.0113e-02, 4.1029e-03, 7.6469e-04,
3974  -1.0138e-03, 6.0837e-04,-2.9929e-04, 1.3329e-04,
3975  -5.5502e-05, 2.1893e-05,-8.2286e-06, 2.9559e-06,
3976  -1.0148e-06};
3977 
3978  //196Pt78 Platinum-196
3979  G4double a_196Pt[] = {0.50218e-1, 0.53722e-1,-0.35015e-1,-0.34588e-1,
3980  0.23564e-1, 0.14340e-1,-0.13270e-1,-0.51212e-2,
3981  0.56088e-2, 0.14890e-2,-0.10928e-2, 0.55662e-3,
3982  -0.50557e-4,-0.19708e-3, 0.24016e-3};
3983 
3984 
3985  // Change to units of fm**(-1)
3986  G4double qq = sqrt ( T ) / 0.197328; // used to be q
3987  // cout << "q(1/fm):: " << qq << endl;
3988 
3989  G4double ff_bessel=0;
3990 
3991  if( qq <= 0.005 ){
3992  OUT_OF_RANGE=false; //! line added 4/29/98
3993  ff_bessel = 1.00000;
3994  return 0;
3995  }
3996 
3997  ff_bessel=0.0;
3998  OUT_OF_RANGE=false;
3999  Float_t R_max=0.;
4000  G4int n=0;
4001 
4002  if(iA == 3 && iZ == 2){ // 3He2
4003  if(qq>10.10 ) OUT_OF_RANGE=true;
4004  R_max=5.;
4005  n = 12;
4006  for(G4int i=0;i<n;i++)
4007  a[i] = a_3He[i];
4008  }
4009  else if(iA == 12 && iZ == 6){ // 12C6
4010  if(qq>4.01 ) OUT_OF_RANGE=true;
4011  R_max=8.;
4012  n = 16;
4013  for(G4int i=0;i<n;i++)
4014  a[i] = a_12C[i];
4015  }
4016  else if(iA == 15 && iZ == 7){ // 15N7
4017  if(qq>3.17) OUT_OF_RANGE=true;
4018  R_max=7.;
4019  n = 8;
4020  for(G4int i=0;i<n;i++)
4021  a[i] = a_15N[i];
4022  }
4023  else if(iA == 16 && iZ == 8){ // 16O8
4024  if(qq>2.77) OUT_OF_RANGE=true;
4025  R_max=8.;
4026  n = 13;
4027  for(G4int i=0;i<n;i++)
4028  a[i] = a_16O[i];
4029  }
4030  else if(iA == 26 && iZ == 26){ // 26Mg12
4031  if(qq>3.0) OUT_OF_RANGE=true;
4032  //if(qq<0.60) OUT_OF_RANGE=true;
4033  R_max=8.;
4034  n = 15;
4035  for(G4int i=0;i<n;i++)
4036  a[i] = a_26Mg[i];
4037  }
4038  else if(iA == 27 && iZ == 13){ // 27Al13
4039  if(qq>2.70) OUT_OF_RANGE=true;
4040  //if(qq<0.47) OUT_OF_RANGE=true;
4041  R_max=7.;
4042  n = 12;
4043  for(G4int i=0;i<n;i++)
4044  a[i] = a_27Al[i];
4045  }
4046  else if(iA == 28 && iZ == 14){ // 28Si14
4047  if(qq>2.64) OUT_OF_RANGE=true;
4048  //if(qq<0.25) OUT_OF_RANGE=true;
4049  R_max = 8.;
4050  n = 14;
4051  for(G4int i=0;i<n;i++)
4052  a[i] = a_28Si[i];
4053  }
4054  else if(iA == 29 && iZ == 14){ // 29Si14
4055  if(qq>2.64) OUT_OF_RANGE=true;
4056  //if(qq<0.25) OUT_OF_RANGE=true;
4057  R_max = 8.;
4058  n = 12;
4059  for(G4int i=0;i<n;i++)
4060  a[i] = a_29Si[i];
4061  }
4062  else if(iA == 30 && iZ == 14){ // 30Si14
4063  if(qq>2.64) OUT_OF_RANGE=true;
4064  //if(qq<0.25) OUT_OF_RANGE=true;
4065  R_max = 8.5;
4066  n = 13;
4067  for(G4int i=0;i<n;i++)
4068  a[i] = a_30Si[i];
4069  }
4070  else if (iA == 48 && iZ == 22) { // 48Ti22
4071  if (qq>2.20) OUT_OF_RANGE=true;
4072  //if (qq<0.15) OUT_OF_RANGE=true;
4073  R_max=10.0;
4074  n=12;
4075  for(G4int i=0;i<n;i++)
4076  a[i] = a_48Ti[i];
4077  }
4078  else if (iA == 50 && iZ == 22) { // 50Ti22
4079  if (qq>2.20) OUT_OF_RANGE=true;
4080  //if (qq<0.15) OUT_OF_RANGE=true;
4081  R_max=9.5;
4082  n=10;
4083  for(G4int i=0;i<n;i++)
4084  a[i] = a_50Ti[i];
4085  }
4086  else if (iA == 50 && iZ == 24) { // 50Cr24
4087  if (qq>2.59) OUT_OF_RANGE=true;
4088  //if (qq<0.15) OUT_OF_RANGE=true;
4089  R_max=9.0;
4090  n=11;
4091  for(G4int i=0;i<n;i++)
4092  a[i] = a_50Cr[i];
4093  }
4094  else if (iA == 52 && iZ == 24) { // 52Cr24
4095  if (qq>2.59) OUT_OF_RANGE=true;
4096  //if (qq<0.15) OUT_OF_RANGE=true;
4097  R_max=9.0;
4098  n=11;
4099  for(G4int i=0;i<n;i++)
4100  a[i] = a_52Cr[i];
4101  }
4102  else if (iA == 54 && iZ == 24) { // 54Cr24
4103  if (qq>2.59) OUT_OF_RANGE=true;
4104  //if (qq<0.15) OUT_OF_RANGE=true;
4105  R_max=9.0;
4106  n=11;
4107  for(G4int i=0;i<n;i++)
4108  a[i] = a_54Cr[i];
4109  }
4110  else if(iA == 54 && iZ == 26){ // 54Fe26
4111  if(qq>2.22) OUT_OF_RANGE=true;
4112  //if(qq<0.51) OUT_OF_RANGE=true;
4113  R_max=9.0;
4114  n = 17;
4115  for(G4int i=0;i<n;i++)
4116  a[i] = a_54Fe[i];
4117  }
4118  else if(iA == 56 && iZ == 26){ // 56Fe26
4119  if(qq>2.22) OUT_OF_RANGE=true;
4120  //if(qq<0.51) OUT_OF_RANGE=true;
4121  R_max=9.0;
4122  n = 17;
4123  for(G4int i=0;i<n;i++)
4124  a[i] = a_56Fe[i];
4125  }
4126  else if(iA == 58 && iZ == 26){ // 58Fe26
4127  if(qq>2.22) OUT_OF_RANGE=true;
4128  //if(qq<0.51) OUT_OF_RANGE=true;
4129  R_max=9.0;
4130  n = 17;
4131  for(G4int i=0;i<n;i++)
4132  a[i] = a_58Fe[i];
4133  }
4134  else if(iA == 63 && iZ == 29){ // 63Cu29
4135  if(qq>2.22) OUT_OF_RANGE=true;
4136  //if(qq<0.51) OUT_OF_RANGE=true;
4137  R_max=9.;
4138  n = 17;
4139  for(G4int i=0;i<n;i++)
4140  a[i] = a_63Cu[i];
4141  }
4142  else if(iA == 65 && iZ == 29){ // 65Cu29
4143  if(qq>2.22) OUT_OF_RANGE=true;
4144  //if(qq<0.51) OUT_OF_RANGE=true;
4145  R_max=9.;
4146  n = 17;
4147  for(G4int i=0;i<n;i++)
4148  a[i] = a_65Cu[i];
4149  }
4150  else if(iA == 64 && iZ == 30){ // 64Zn30
4151  if(qq>2.22) OUT_OF_RANGE=true;
4152  //if(qq<0.51) OUT_OF_RANGE=true;
4153  R_max=9.;
4154  n = 17;
4155  for(G4int i=0;i<n;i++)
4156  a[i] = a_64Zn[i];
4157  }
4158  else if(iA == 66 && iZ == 30){ // 66Zn30
4159  if(qq>2.22) OUT_OF_RANGE=true;
4160  //if(qq<0.51) OUT_OF_RANGE=true;
4161  R_max=9.;
4162  n = 17;
4163  for(G4int i=0;i<n;i++)
4164  a[i] = a_66Zn[i];
4165  }
4166  else if(iA == 68 && iZ == 30){ // 68Zn30
4167  if(qq>2.22) OUT_OF_RANGE=true;
4168  //if(qq<0.51) OUT_OF_RANGE=true;
4169  R_max=9.;
4170  n = 17;
4171  for(G4int i=0;i<n;i++)
4172  a[i] = a_68Zn[i];
4173  }
4174  else if(iA == 70 && iZ == 30){ // 70Zn30
4175  if(qq>2.22) OUT_OF_RANGE=true;
4176  //if(qq<0.51) OUT_OF_RANGE=true;
4177  R_max=9.;
4178  n = 17;
4179  for(G4int i=0;i<n;i++)
4180  a[i] = a_70Zn[i];
4181  }
4182  else if(iA == 196 && iZ == 78){ // 196Pt78
4183  if(qq>2.28) OUT_OF_RANGE=true;
4184  if(qq<0.34) OUT_OF_RANGE=true;
4185  R_max=12.;
4186  n = 15;
4187  for(G4int i=0;i<n;i++)
4188  a[i] = a_196Pt[i];
4189  }
4190  else
4191  OUT_OF_RANGE=true;
4192 
4193  if(OUT_OF_RANGE || (R_max == 0.)){
4194  ff_bessel=0.;
4195  return 0;
4196  }
4197 
4198  G4double qmu = 3.14159265 / R_max;
4199  G4double ff=0.;
4200  G4double sinqR = sin(qq*R_max);
4201  G4double qmux,qm,qp;
4202 
4203  G4double PI4=12.56637062; //! 4*PI
4204  G4double PI2=6.28318531; //! 2*PI
4205 
4206  for(G4int j=0;j<n;j++){
4207  qmux = qmu * Float_t(j+1);
4208  qm = qq - qmux;
4209  qp = qq + qmux;
4210  if(fabs(qm)>1.E-6)
4211  ff = ff + a[j]*(pow((-1.),(j+1)))*sinqR/(qm*qp);
4212  else
4213  ff = ff + a[j]*pow(R_max,2)/(PI2*(j+1));
4214  }
4215 
4216  if((qq*R_max>1.E-20) && (ff<1.E20))
4217  ff_bessel = PI4/Float_t(iZ)/qq *ff;
4218  else
4219  ff_bessel = 0.;
4220 
4221  if (ff_bessel < 0) ff_bessel*=-1;
4222 
4223  return ff_bessel;
4224 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::Fgauss ( G4int  Z,
G4int  A,
G4double  q2 
)

from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987) 8/9/96

Definition at line 3779 of file QweakSimEPEvent.cc.

Referenced by AlloyScattering().

3780 {
3781  // cout << "** Calling Fgauss **" << endl;
3782 
3783  G4double fgauss = 0.;
3784 
3785  G4double avgA = iA;
3786  G4double radius = 1.07*pow(avgA,1./3.);
3787 
3788  //! from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron
3789  //! Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987)
3790  //! 8/9/96
3791  if(iA==205)radius=5.470;
3792  if(iA==56) radius=3.729;
3793  if(iA==28) radius=3.085;
3794  if(iA==27) radius=3.035;
3795  if(iA==24) radius=3.08;// Mg-24
3796 
3797  G4double x2 = (T/pow(.197328,2))*pow(radius,2);
3798  G4double CHARR = (T/pow(.197328,2))*(2.4*2.4)/6.;
3799 
3800  //cout << " T (Q^2), CHARR: " << T << ", " << CHARR << endl;
3801  // cout << "** Fgauss CHARR:: " << CHARR << endl;
3802  if (CHARR<80) fgauss = exp(-CHARR)/(1.+x2/6.);
3803 
3804  return fgauss;
3805 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::fitemc ( G4double  X,
G4int  A 
)

! Fit to EMC effect. Steve Rock 8/3/94 ! Funciton returns value of sigma(A)/sigma(d) ! with no isoscalerity correction ! A= atomic number ! x = Bjorken x. ! ! Fit of sigma(A)/sigma(d) to form C*A**alpha where A is atomic number ! First data at each x was fit to form C*A**alpha. The point A=2(d) ! was includded with a value of 1 and an error of about 2%. ! For x>=.125 Javier Gomez fit of 7/93 to E139 data was used. ! For .09 >=x>=.0085 NMC data from Amaudruz et al Z. Phys C. 51,387(91) ! Steve did the fit for alpha and C to the He/d. C/d and Ca/d NMC data. ! Alpha(x) was fit to a 9 term polynomial a0 +a1*x +a2*x**2 + a3*x**3 .. ! C(x) was fit to a 3 term polynomial in natural logs as ! Ln(C) = c0 + c1*Ln(x) + c2*[Ln(x)]**2.

! 6/2/98 ***** Bug (which set x= .00885 if x was out of range) fixed ! also gave value at x=.0085 if x>.88 ! 11/05 PYB PEB modified to use value at x=0.7 if x>0.7, because ! beyond that assume Fermi motion is giving the rise, and we ! already are taking that into account with the y-smearing of ! the inelastic !--------------------------------------------------------------------—

Definition at line 3626 of file QweakSimEPEvent.cc.

Referenced by F1F2IN09().

3627 {
3628  /*!---------------------------------------------------------------------
3629  ! Fit to EMC effect. Steve Rock 8/3/94
3630  ! Funciton returns value of sigma(A)/sigma(d)
3631  ! with no isoscalerity correction
3632  ! A= atomic number
3633  ! x = Bjorken x.
3634  !
3635  ! Fit of sigma(A)/sigma(d) to form C*A**alpha where A is atomic number
3636  ! First data at each x was fit to form C*A**alpha. The point A=2(d)
3637  ! was includded with a value of 1 and an error of about 2%.
3638  ! For x>=.125 Javier Gomez fit of 7/93 to E139 data was used.
3639  ! For .09 >=x>=.0085 NMC data from Amaudruz et al Z. Phys C. 51,387(91)
3640  ! Steve did the fit for alpha and C to the He/d. C/d and Ca/d NMC data.
3641  ! Alpha(x) was fit to a 9 term polynomial a0 +a1*x +a2*x**2 + a3*x**3 ..
3642  ! C(x) was fit to a 3 term polynomial in natural logs as
3643  ! Ln(C) = c0 + c1*Ln(x) + c2*[Ln(x)]**2.
3644 
3645  ! 6/2/98 ***** Bug (which set x= .00885 if x was out of range) fixed
3646  ! also gave value at x=.0085 if x>.88
3647  ! 11/05 PYB PEB modified to use value at x=0.7 if x>0.7, because
3648  ! beyond that assume Fermi motion is giving the rise, and we
3649  ! already are taking that into account with the y-smearing of
3650  ! the inelastic
3651  !-----------------------------------------------------------------------
3652  */
3653  G4double ALPHA,C,LN_C,X_U;
3654  G4double fitval;
3655 
3656  // !Chisq= 19. for 30 points
3657  // !Term Coeficient
3658  G4double ALPHA_COEF[9] = {
3659  -6.98871401E-02,
3660  2.18888887E+00,
3661  -2.46673765E+01,
3662  1.45290967E+02,
3663  -4.97236711E+02,
3664  1.01312929E+03,
3665  -1.20839250E+03,
3666  7.75766802E+02,
3667  -2.05872410E+02};
3668 
3669  // !Chisq= 22. for 30 points
3670  // !Term Coeficient
3671  G4double C_COEF[3] = {
3672  1.69029097E-02,
3673  1.80889367E-02,
3674  5.04268396E-03};
3675 
3676  fitval = 1.;
3677  if (A < 2.5) return 0;
3678 
3679  if ( (X > 0.70) || (X < 0.0085) ) {
3680  // !Out of range of fit
3681  if (X < 0.0085) X_U =.0085;
3682  if (X > 0.70) X_U = 0.70;
3683  } else {
3684  X_U = X;
3685  }
3686 
3687  LN_C = C_COEF[0];
3688  for (G4int i = 1; i <= 2; i++) {
3689  LN_C = LN_C + C_COEF[i] * pow(log(X_U),i);
3690  }
3691 
3692  C = exp(LN_C);
3693 
3694  ALPHA = ALPHA_COEF[0];
3695  for (G4int i = 1; i <= 8; i++) {
3696  ALPHA = ALPHA + ALPHA_COEF[i] * pow(X_U,i);
3697  }
3698 
3699  fitval = C*pow(A,ALPHA);
3700  return fitval;
3701 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::Fshell ( G4int  Z,
G4int  A,
G4double  q2 
)

from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987) 8/9/96

Definition at line 3807 of file QweakSimEPEvent.cc.

Referenced by AlloyScattering().

3808 {
3809  // cout << "** Calling Fshell **" << endl;
3810 
3811  G4double fshell = 0.;
3812  G4double avgA = iA;
3813  G4double radius = 1.07*pow(avgA,1./3.);
3814  //! from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron
3815  //! Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987)
3816  //!8/9/96
3817  if(iA==16) radius=2.737;
3818  if(iA==15) radius=2.611;
3819  if(iA==14) radius=2.540;
3820  if(iA==12) radius=2.464;
3821  if(iA== 9) radius=2.519;
3822  if(iA== 7) radius=2.39;
3823  if(iA== 6) radius=2.56;
3824  if(iA== 4) radius=1.676;
3825  G4double x2 = (T/pow(.197328,2))*pow(radius,2);
3826  G4double alp = (iZ-2.)/3.;
3827  G4double CHARR = x2*(2.+3.*alp)/(12.+30.*alp);
3828 
3829  // cout << "** Fshell CHARR:: " << CHARR << endl;
3830  if (CHARR<80) fshell = exp(-CHARR)*(1.-alp*x2/(6.+15.*alp));
3831 
3832  return fshell;
3833 }

+ Here is the caller graph for this function:

G4int QweakSimEPEvent::GetActiveOctantNumber ( )
inline

Definition at line 97 of file QweakSimEPEvent.hh.

References kActiveOctantNumber.

97 {return kActiveOctantNumber;};
G4int kActiveOctantNumber
Active octant number in the simulation, 0 will enable all octants.
void QweakSimEPEvent::GetanEvent ( G4double  E_in,
std::vector< G4double > &  CrossSection,
G4double &  weight_n,
G4double &  Q2,
G4double &  E_out,
G4ThreeVector &  MomentumDirection,
G4double &  theta,
G4double &  phi,
G4double &  Asymmetry 
)

Referenced by QweakSimSteppingAction::UserSteppingAction().

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetAsymmetry_AL ( G4double  theta,
G4double  energy 
)

Definition at line 4319 of file QweakSimEPEvent.cc.

4320 {
4321 // Needed constants
4322  const G4double gf=0.0000116637;
4323  const G4double alpha=1./137.;
4324  const G4double qwp=0.0713;
4325  const G4double qwn=-0.988;
4326  const G4double NT=14.;
4327  const G4double ZT=13.;
4328 
4329 // The energy is in MeV after the cross section subroutine -
4330 // change back to GeV
4331  energy = energy/1000.0;
4332 
4333 // Claculate Q2
4334  G4double Q2=4*energy*energy*sin(theta/2.)*sin(theta/2.);
4335 
4336 // Correction for Coulomb distortion for Al
4337 // Q2=Q2*(1.+(3.*ZT*0.197/137./(2.*energy*2.98)))*(1.+(3.*ZT*0.197/137./(2.*energy*2.98)));
4338 
4339 // Aluminum Asymmetry
4340  G4double asym=-gf/(4.*pi*alpha*sqrt(2.))*1e6*Q2*(qwp+qwn*NT/ZT);
4341 
4342  return asym;
4343 }
G4double QweakSimEPEvent::GetAsymmetry_Be ( G4double  theta,
G4double  energy 
)

Definition at line 4357 of file QweakSimEPEvent.cc.

4358 {
4359  const G4double Mp=0.938;
4360  const G4double gf=0.0000116637;
4361  const G4double alpha=1./137.;
4362  const G4double qwp=0.0713;
4363  const G4double qwn=-0.988;
4364  const G4double NT=5.;
4365  const G4double ZT=4.;
4366 
4367 // The energy is in MeV after the cross section subroutine -
4368 // change back to GeV
4369  energy = energy/1000.;
4370 
4371 // Claculate Q2
4372  G4double Q2=4*energy*energy*sin(theta/2.)*sin(theta/2.);
4373  Q2=Q2/(1.+2.*energy/(9.*Mp)*sin(theta/2.)*sin(theta/2.));
4374 
4375 // Correction for Coulomb distortion
4376  Q2=Q2*(1.+(3./2.*sqrt(3./5.)*ZT*0.197/137./(energy*1.77)))*(1.+(3./2.*sqrt(3./5.)*ZT*0.197/137./(energy*1.77)));
4377 
4378 // Beryllium Asymmetry
4379  G4double asym=-gf/(4.*pi*alpha*sqrt(2.))*1e6*Q2*(qwp+qwn*NT/ZT);
4380 
4381  return asym;
4382 }
G4double QweakSimEPEvent::GetAsymmetry_EN ( G4double  theta,
G4double  energy 
)

Definition at line 4436 of file QweakSimEPEvent.cc.

4437 {
4438 
4439  const G4double Mn= 0.939565; // Mp=0.938;
4440  const G4double gf=0.0000116637;
4441  const G4double alpha=1./137.;
4442  //const G4double s2tw=0.231;
4443  const G4double s2tw_low=0.2387;
4444  const G4double qwp=0.0713;
4445  const G4double qwn=-0.988;
4446  //const G4double NT=14.;
4447  //const G4double ZT=13.;
4448 
4449 // Radiative correction factors [from Musolf's Phys. Rep. 1994]
4450  //const G4double rpv=-0.054;
4451  //const G4double rnv=-0.0143;
4452  //const G4double rt0a=0.0;
4453  const G4double rt1a=-0.23;
4454  //const G4double r0a=0.023;
4455 
4456 // Check for minimum theta
4457  const G4double theta_min = 1.745329E-4;
4458  if (theta<theta_min)
4459  theta = theta_min;
4460 
4461 // The energy is in MeV after the cross section subroutine -
4462 // change back to GeV
4463  energy=energy/1000.0;
4464 
4465 // Calculate Q2
4466  G4double Q2_ep=4*energy*energy*sin(theta/2.0)*sin(theta/2.0);
4467  Q2_ep=Q2_ep/(1.0+2.0*energy/Mn*sin(theta/2.0)*sin(theta/2.0));
4468 
4469 // Proton Asymmetry
4470 // Kinematic Factors
4471  G4double tau=Q2_ep/4.0/(Mn*Mn);
4472  G4double epsilon=1.0/(1.0+2.0*(1.0+tau)*tan(theta/2.0)*tan(theta/2.0));
4473  G4double epsilonp=sqrt(tau*(1.0+tau)*(1.0-epsilon*epsilon));
4474 
4475 // Form factor calculations
4476 // Form Factors: Electromagetic
4477  G4double gvd=1.0/( (1.0+Q2_ep/0.71)*(1.0+Q2_ep/0.71) );
4478  G4double gepg=gvd;
4479  G4double gmpg=2.793*gvd;
4480  G4double geng=1.91*tau*gvd/(1.+5.6*tau);
4481  G4double gmng=-1.91*gvd;
4482  //G4double gesg= 0;
4483  //G4double gmsg= 0;
4484 
4485 // Form Factors: Neutral-weak, Axial
4486 // Assume: Gs_E,M=0, G8_A,Gs_A=0
4487  //G4double gad=1/((1.0+Q2_ep/1.001/1.001)*(1.0+Q2_ep/1.001/1.001));
4488  //G4double gsa=-0.12/((1.+Q2_ep/1.06/1.06)*(1.+Q2_ep/1.06/1.06));
4489  G4double gt1a=1.2695/((1+Q2_ep/1.001/1.001)*(1+Q2_ep/1.001/1.001));
4490  //G4double g8a=0.0;
4491  //gsa=0.0;
4492 
4493 // Use SM radiative correction factors
4494 // Since we use SM values for Qw(p) and Qw(n)
4495 // we no longer need the radiative corrections
4496 // on those quantities.
4497 // gepz=(1.-4.*s2tw)*(1.+rpv)*gepg-(1.+rnv)*geng
4498 // gmpz=(1.-4.*s2tw)*(1.+rpv)*gmpg-(1.+rnv)*gmng
4499 
4500  G4double gapz=-(1.0+rt1a)*gt1a;
4501 
4502 // G4double gepz=qwp*gepg+qwn*geng;
4503 // G4double gmpz=qwp*gmpg+qwn*gmng;
4504 // G4double asym=-gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4505 // asym=asym/(epsilon*gepg*gepg+tau*gmpg*gmpg);
4506 // asym=asym*(epsilon*gepg*gepz+tau*gmpg*gmpz-(1.0-4.0*s2tw_low)*epsilonp*gmpg*gapz);
4507 // //std::cout<<"asym_proton="<<asym<<std::endl;
4508 
4509  G4double genz=qwp*geng+qwn*gepg;
4510  G4double gmnz=qwp*gmng+qwn*gmpg;
4511  G4double asym_neutron = -gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4512  asym_neutron =asym_neutron /(epsilon*geng*geng+tau*gmng*gmng);
4513  asym_neutron =asym_neutron *(epsilon*geng*genz+tau*gmng*gmnz-(1.0-4.0*s2tw_low)*epsilonp*gmng*gapz);
4514  //std::cout<<"asym_neutron="<<asym_neutron<<std::endl<<std::endl;
4515 
4516  return asym_neutron;
4517 }
G4double QweakSimEPEvent::GetAsymmetry_EP ( G4double  theta,
G4double  energy 
)

Definition at line 4233 of file QweakSimEPEvent.cc.

4234 {
4235 
4236  const G4double Mp=0.938;
4237  const G4double gf=0.0000116637;
4238  const G4double alpha=1./137.;
4239  //const G4double s2tw=0.231;
4240  const G4double s2tw_low=0.2387;
4241  const G4double qwp=0.0713;
4242  const G4double qwn=-0.988;
4243  //const G4double NT=14.;
4244  //const G4double ZT=13.;
4245 
4246 // Radiative correction factors [from Musolf's Phys. Rep. 1994]
4247  //const G4double rpv=-0.054;
4248  //const G4double rnv=-0.0143;
4249  //const G4double rt0a=0.0;
4250  const G4double rt1a=-0.23;
4251  //const G4double r0a=0.023;
4252 
4253 // Check for minimum theta
4254  const G4double theta_min = 1.745329E-4;
4255  if (theta<theta_min)
4256  theta = theta_min;
4257 
4258 // The energy is in MeV after the cross section subroutine -
4259 // change back to GeV
4260  energy=energy/1000.0;
4261 
4262 // Calculate Q2
4263  G4double Q2_ep=4*energy*energy*sin(theta/2.0)*sin(theta/2.0);
4264  Q2_ep=Q2_ep/(1.0+2.0*energy/Mp*sin(theta/2.0)*sin(theta/2.0));
4265 
4266 // Proton Asymmetry
4267 // Kinematic Factors
4268  G4double tau=Q2_ep/4.0/(Mp*Mp);
4269  G4double epsilon=1.0/(1.0+2.0*(1.0+tau)*tan(theta/2.0)*tan(theta/2.0));
4270  G4double epsilonp=sqrt(tau*(1.0+tau)*(1.0-epsilon*epsilon));
4271 
4272 // Form factor calculations
4273 // Form Factors: Electromagetic
4274  G4double gvd=1.0/( (1.0+Q2_ep/0.71)*(1.0+Q2_ep/0.71) );
4275  G4double gepg=gvd;
4276  G4double gmpg=2.793*gvd;
4277  G4double geng=1.91*tau*gvd/(1.+5.6*tau);
4278  G4double gmng=-1.91*gvd;
4279 
4280 // Form Factors: Neutral-weak, Axial
4281 // Assume: Gs_E,M=0, G8_A,Gs_A=0
4282  //G4double gad=1/((1.0+Q2_ep/1.001/1.001)*(1.0+Q2_ep/1.001/1.001));
4283  //G4double gsa=-0.12/((1.+Q2_ep/1.06/1.06)*(1.+Q2_ep/1.06/1.06));
4284  G4double gt1a=1.2695/((1+Q2_ep/1.001/1.001)*(1+Q2_ep/1.001/1.001));
4285  //G4double g8a=0.0;
4286  //gsa=0.0;
4287 
4288 // Use SM radiative correction factors
4289 // Since we use SM values for Qw(p) and Qw(n)
4290 // we no longer need the radiative corrections
4291 // on those quantities.
4292 // gepz=(1.-4.*s2tw)*(1.+rpv)*gepg-(1.+rnv)*geng
4293 // gmpz=(1.-4.*s2tw)*(1.+rpv)*gmpg-(1.+rnv)*gmng
4294 
4295  G4double gepz=qwp*gepg+qwn*geng;
4296  G4double gmpz=qwp*gmpg+qwn*gmng;
4297  G4double gapz=-(1.0+rt1a)*gt1a;
4298 
4299  G4double asym=-gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4300  asym=asym/(epsilon*gepg*gepg+tau*gmpg*gmpg);
4301  asym=asym*(epsilon*gepg*gepz+tau*gmpg*gmpz-(1.0-4.0*s2tw_low)*epsilonp*gmpg*gapz);
4302 
4303  return asym;
4304 }
G4double QweakSimEPEvent::GetAsymmetry_Pi ( G4double  Q2_pi)

Definition at line 4531 of file QweakSimEPEvent.cc.

4532 {
4533 // Inelastic Asymmetry
4534  Q2_pi = Q2_pi/1e6;
4535  G4double asym=-100*Q2_pi;
4536  return asym;
4537 }
G4double QweakSimEPEvent::GetBeamEnergy ( )

Definition at line 4582 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetBeamEnergy(), and myUserInfo.

Referenced by CheckLookupTableBounds(), and CreateLookupTable().

4582 {return myUserInfo->GetBeamEnergy();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetElasticPeakDeltaE ( )
inline

Definition at line 163 of file QweakSimEPEvent.hh.

References ElasticPeakDeltaE.

163 { return ElasticPeakDeltaE; } ;
G4double ElasticPeakDeltaE
G4double QweakSimEPEvent::GetEPrime_Max ( )

Definition at line 4553 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetEPrime_Max(), and myUserInfo.

Referenced by CheckLookupTableBounds(), and Pion_PhotoProductionCarbon().

4553 {return myUserInfo->GetEPrime_Max();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetEPrime_Min ( )

Definition at line 4550 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetEPrime_Min(), and myUserInfo.

Referenced by CheckLookupTableBounds(), and Pion_PhotoProductionCarbon().

4550 {return myUserInfo->GetEPrime_Min();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4int QweakSimEPEvent::GetIsotropy ( )
inline

Definition at line 100 of file QweakSimEPEvent.hh.

References Isotropy.

100 {return Isotropy;};
G4int Isotropy
isotropy used for event generation
G4double QweakSimEPEvent::GetLuminosity ( )

Definition at line 4568 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetLuminosity(), and myUserInfo.

4568 {return myUserInfo->GetLuminosity();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4ThreeVector QweakSimEPEvent::GetMomentumDirection ( )
private
G4double QweakSimEPEvent::GetPhaseSpace ( )

Definition at line 4571 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetPhaseSpace(), and myUserInfo.

4571 {return myUserInfo->GetPhaseSpace();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4double QweakSimEPEvent::GetPhiAngle_Max ( )

Definition at line 4565 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetPhiAngle_Max(), and myUserInfo.

4565 {return myUserInfo->GetPhiAngle_Max();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4double QweakSimEPEvent::GetPhiAngle_Min ( )

Definition at line 4562 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetPhiAngle_Min(), and myUserInfo.

4562 {return myUserInfo->GetPhiAngle_Min();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4int QweakSimEPEvent::GetReactionRegion ( )
inline

Definition at line 175 of file QweakSimEPEvent.hh.

References ReactionRegion.

175 {return ReactionRegion; };
G4int ReactionRegion
reaction region used for event generation
G4int QweakSimEPEvent::GetReactionType ( )
inline

Definition at line 172 of file QweakSimEPEvent.hh.

References ReactionType.

Referenced by QweakSimSteppingAction::UserSteppingAction().

172 {return ReactionType; };
G4int ReactionType
reaction type used for event generation

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetSchwingerDeltaE ( )
inline

Definition at line 168 of file QweakSimEPEvent.hh.

References SchwingerDeltaE.

168 { return SchwingerDeltaE; } ;
G4double SchwingerDeltaE
G4double QweakSimEPEvent::GetThetaAngle_Max ( )

Definition at line 4559 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetThetaAngle_Max(), and myUserInfo.

Referenced by CheckLookupTableBounds().

4559 {return myUserInfo->GetThetaAngle_Max();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetThetaAngle_Min ( )

Definition at line 4556 of file QweakSimEPEvent.cc.

References QweakSimUserInformation::GetThetaAngle_Min(), and myUserInfo.

Referenced by CheckLookupTableBounds(), NuclearInelastic_Bosted(), Pion_PhotoProduction(), Pion_PhotoProductionAl(), Pion_PhotoProductionCarbon(), Quasi_Elastic_Bosted(), and Radiative_Cross_Section_Lookup().

4556 {return myUserInfo->GetThetaAngle_Min();}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::GetVertexZ ( )
G4double QweakSimEPEvent::Horowitz_DW_Asym ( G4double  theta,
const G4String  path 
)

Definition at line 2506 of file QweakSimEPEvent.cc.

2508 {
2509  G4double radtodeg = (180.0/3.1415927);
2510  G4double asym = 0.0; // Asymmetry
2511  TGraph *graph = new TGraph(path, "%lg %*lg %lg", ",");
2512  graph->SetBit(graph->kIsSortedX); // Faster queires on TGraph object.
2513  if((theta*radtodeg) >= 4.0 && (theta*radtodeg) <= 13.5)
2514  {
2515  asym = graph->Eval((theta*radtodeg), 0, "S"); // Spline Interpolation
2516  }
2517  delete graph; // Delete TGraph object.
2518  return asym*1e6; // [ppm]
2519 }
G4double QweakSimEPEvent::Horowitz_DW_Xsect ( G4double  E_in,
G4double  theta,
const G4String  path,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2464 of file QweakSimEPEvent.cc.

Referenced by AlGDR().

2470 {
2471  const G4double M = 13.0*938.272081 + 14.0*939.565413; //2016 PDG - Al mass [MeV]
2472  const G4double alpha = 1.0/137.035999139; // Fine Structure Const 2016 PDG
2473  const G4double deltaE = 20.0*MeV; // Schwinger correction delta E parameter
2474  const G4double m_e = 0.5109989461*MeV; // Electron mass 2016 PDG
2475  // Calculate kinematic variables
2476  G4double STH = sin(theta/2.0);
2477  G4double ETA = 1.0+2.0*(E_in/M)*STH*STH;
2478  E_out = E_in/ETA; //E prime [MeV]
2479  Q2 = 4.0*E_in*E_out*STH*STH; //[MeV^2]
2480 
2481  // Interpolate Horowitz Xsect calculation
2482  G4double radtodeg = (180.0/TMath::Pi());
2483  G4double xsect = 0.0;
2484  TGraph *graph = new TGraph(path, "%lg %lg %*lg", ",");
2485  graph->SetBit(graph->kIsSortedX); // Faster queires on TGraph object.
2486  // Check theta value inside of table range (4.0-13.5 [deg])
2487  if((theta*radtodeg) >= 4.0 && (theta*radtodeg) <= 13.5)
2488  {
2489  xsect = graph->Eval((theta*radtodeg), 0, "S"); // Spline Interpolation
2490  }
2491  // Calculate Schwinger Radiative Correction
2492  // Using eq.6-8c pg.487 from "Electron Scattering from Complex Nuclei"
2493  // Part B Herbert Uberall - 2017/7/26 K. Bartlett
2494  G4double delta_Schwinger = 0.0;
2495  delta_Schwinger = (TMath::Log(E_in/deltaE)-(13.0/12.0))*(TMath::Log(Q2/(m_e*m_e))-1.0);
2496  delta_Schwinger += (17.0/36.0) + 0.5*((1.0/6.0)*TMath::Pi()*TMath::Pi()-TMath::DiLog(TMath::Cos(theta/2.0)*TMath::Cos(theta/2.0)));
2497  delta_Schwinger *= (2.0*alpha/TMath::Pi());
2498 
2499  xsect *= TMath::Exp(-delta_Schwinger)*1000.0; //Convert from [mb/sr] to [ub/sr]
2500 
2501  fWeightN = xsect*sin(theta);
2502  delete graph; // Delete TGraph object.
2503  return xsect;
2504 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::MEC2009 ( G4int  a,
G4double  q2,
G4double  w2 
)

fit to low q2 dip region: purefly empirical assume contribution is pure transverse

Definition at line 3581 of file QweakSimEPEvent.cc.

Referenced by F1F2IN09().

3582 {
3583  //! fit to low q2 dip region: purefly empirical
3584  //! assume contribution is pure transverse
3585  G4double f1 = 0.0;
3586  G4double am = 0.9383;
3587  G4double w,nu;
3588 
3589  G4double P[24] = {
3590  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
3591  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
3592  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
3593  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
3594  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
3595  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
3596 
3597  G4double p18, x, f1corr;
3598 
3599  if (w2 <= 0.0) return 0;
3600  w = pow(w2,0.5);
3601  nu = (w2 - am*am + q2) / 2. / am;
3602  x = q2 / (2.0 * am * nu );
3603 
3604  if (a < 2.5) return 0;
3605 
3606  p18 = P[18];
3607  // ! special case for 3He
3608  if (a > 2.5 && a < 3.5) p18 = 70;
3609  // ! special case for 4He
3610  if (a > 3.5 && a < 4.5) p18 = 170;
3611  // ! new values for C, Al, Cu
3612  if(a > 4.5) p18 = 215;
3613  if(a > 20.) p18 = 235;
3614  if(a > 50.) p18 = 230;
3615 
3616 
3617  f1corr = P[0]*exp(-pow(w-P[1],2)/P[2])/(pow(1.0 + std::max(0.3,q2)/P[3],P[4]))*pow(nu,P[5])*(1.0 + p18 * pow(a,(1.0 + P[19] * x)));
3618 
3619  f1 = f1corr;
3620 
3621  if (f1 <= 1.0E-9 ) f1 = 0.0;
3622 
3623  return f1;
3624 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::Moller_Scattering ( G4double  E_in,
G4double  theta1,
G4double &  E_out1,
G4double &  E_out2,
G4double &  theta2,
G4double &  q2,
G4double &  fWeightN,
G4double &  asymmetry 
)

Definition at line 1389 of file QweakSimEPEvent.cc.

1392 {
1393  G4double alpha = 1.0/137.036;
1394  G4double G_F = 1.166e-11; // Fermi constant in MeV^(-2)
1395  G4double M_electron = 0.511; // Incident Electron mass in MeV
1396  G4double sin2_theta_W = 0.02381; // sin2_theta_W at 0.026 GeV, taken from E158
1397 
1398  G4double theta_CM = 2.0*atan(tan(theta1)*sqrt((E_in+M_electron)/(2.*M_electron)));
1399  E_out1 = M_electron+(E_in-M_electron)*(cos(theta_CM/2.0)*cos(theta_CM/2.0));
1400  E_out2 = E_in - E_out1 + M_electron;
1401  theta2 = atan(2.0*M_electron/(E_in+M_electron)/tan(theta1));
1402  //G4double momentum = sqrt(E_out1*E_out1 - M_electron*M_electron);
1403  // compute q2 in MeV^2
1404  //q2 = 4.0*E_in*momentum*sin(theta1/2.0)*sin(theta1/2.0);
1405  q2 = 2*M_electron*(M_electron+E_in)*sin(theta_CM/2.0)*sin(theta_CM/2.0);
1406 
1407  // calculate Moller cross-section (units are ub/sr)
1408  G4double Xsect = pow((1+cos(theta_CM))*(3+cos(theta_CM)*cos(theta_CM)),2);
1409  // Xsect = pow(alpha/2.0/M_electron,2)*Xsect/pow(sin(theta_CM),4); // alpha^2/(4m^2) = 0.199 b/Sr
1410  Xsect = (1.99e4)*Xsect/pow(sin(theta_CM),4);
1411  fWeightN = Xsect*sin(theta1);
1412  asymmetry = -M_electron*E_in*(G_F/(sqrt(2.0)*alpha*3.1415927))
1413  *16*sin(theta_CM)*sin(theta_CM)/(3+cos(theta_CM)*cos(theta_CM))*(1-4*sin2_theta_W);
1414  asymmetry = asymmetry*1e6; // in [ppm]
1415 
1416  //std::cout<<"E_in E1 E2 theta1 theta2 Xsect Weight Q2 asym"<<std::endl;
1417  //std::cout<<E_in<<" "<<E_out1<<" "<<E_out2<<" "<<theta1*180/3.14<<" "
1418  // <<theta2*180/3.14<<" "<<Xsect<<" "<<fWeightN<<" "<<q2*1e-6<<" "<<asymmetry<<std::endl;
1419 
1420  return Xsect;
1421 }
G4double QweakSimEPEvent::NuclearInelastic_Bosted ( G4double  E_in,
G4double  Theta,
G4int  Zin,
G4int  Ain,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2593 of file QweakSimEPEvent.cc.

References F1F2IN09(), and GetThetaAngle_Min().

2600 {
2601 
2602  G4double F1 = 0.0;
2603  G4double F2 = 0.0;
2604  G4double W1 = 0.0;
2605  G4double W2 = 0.0;
2606  G4double Wsq = 0.0;
2607  G4double xsect = 0.0;
2608  G4double Ein = 0.0;
2609  G4double Eout = 0.0;
2610 
2611  // In this subroutine units are in GeV, so
2612  Ein = E_in/1000.; // maximum energy in GeV
2613  G4double M_electron = 0.000511; // minimum energy (electron mass) in GeV
2614  G4double Mp = 0.93828;
2615 
2616  // G4double Theta_Min = 1.745329E-4;
2617  if (Theta<GetThetaAngle_Min())
2618  Theta = GetThetaAngle_Min();
2619 
2620  // Generate flat energy distribution of outgoing electron
2621  Eout = M_electron + G4UniformRand()*(Ein - M_electron);
2622 
2623  G4double CTH = cos(Theta/2.0);
2624  G4double STH = sin(Theta/2.0);
2625  G4double T2THE = STH*STH/CTH/CTH;
2626  G4double Nu = Ein - Eout;
2627  Q2 = 4.0*Ein*Eout*STH*STH;
2628  Wsq = Mp*Mp + 2.0*Mp*Nu - Q2;
2629 
2630  // Mott scattering
2631  G4double MOTT = pow((0.00072/Ein*CTH/STH/STH),2);
2632  MOTT = MOTT*1.0E4; // Units: ub/sr/GeV
2633 
2634  //G4cout << "Ein, Eout, Nu, Theta: " << Ein << ", " << Eout << ", " << Nu << ", " << Theta/degree << G4endl;
2635  //G4cout << "Q2, W2, Mott: " << Q2 << ", " << Wsq << ", " << MOTT << G4endl;
2636 
2637  F1F2IN09(Zin, Ain, Q2, Wsq, F1, F2);
2638  W1 = F1/Mp;
2639  W2 = F2/Nu;
2640 
2641  xsect = MOTT*(W2 + 2.0*T2THE*W1)*(Ein - M_electron);
2642  //G4cout << "F1, F2, W1, W2: " << F1 << ", " << F2 << ", " << W1 << ", " << W2 << G4endl;
2643 
2644  //G4cout << E_in << "\t" << Theta << "\t" << MOTT*(Ein - M_electron)*Zin*Zin << "\t" << (W2 + 2.0*T2THE*W1)/(Zin*Zin) << G4endl;
2645 
2646  // In some cases a negative F2 is returned giving a negative cross section
2647  if (xsect <= 0) xsect = 0.0;
2648 
2649  fWeightN = xsect*sin(Theta);
2650  //G4cout << "xsect, fweight: " << xsect << ", "<< fWeightN << G4endl;
2651 
2652  E_out = Eout*1000.; // Convert back to MeV before leaving
2653  Q2 = Q2*1.0E6; // Seems like this needs to be returned in MeV**2
2654  return xsect;
2655 }
G4double GetThetaAngle_Min()
void F1F2IN09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Pion_PhotoProduction ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2169 of file QweakSimEPEvent.cc.

References GetThetaAngle_Min(), myPositionZ, myUserInfo, QweakSimUserInformation::TargetCenterPositionZ, QweakSimUserInformation::TargetEntranceWindowThickness, QweakSimUserInformation::TargetLength, and wiser_sigma().

2174 {
2175  const G4double Mpi = 139.57018*MeV;
2176 
2177  if (Theta<GetThetaAngle_Min())
2178  {
2179  Theta = GetThetaAngle_Min();
2180  }
2181 
2182  if(E_in < Mpi) E_in = Mpi;
2183  E_out = Mpi + G4UniformRand()*(E_in - Mpi);
2184  G4double pf = sqrt(pow(E_out,2.0) - pow(Mpi,2.0)); //--- final momentum in MeV
2185 
2186  // Calculate Q2
2187  G4double STH = sin(Theta/2.0);
2188  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
2189 
2190  //---
2191  //--- radiation length from page 12 of "The Qweak target design and safety document"
2192  //--- http://qweak.jlab.org/DocDB/0010/001041/002/Qweak%20Target%20PDR.pdf
2193  //---
2194  G4double lh2_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength);
2195  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2196  G4double rad_len = myUserInfo->TargetEntranceWindowThickness/(8.896*cm) + lh2_length/(871.9*cm) + intern_rad_len; //--- radiation length
2197  //G4cout << "radiation length: " << rad_len << G4endl;
2198 
2199  G4int type = 1; //--- 0: pi+, 1: pi-
2200 
2201  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2202  // rad_len in fraction (not %)
2203  // wiser_sigma returns crsec in nb/GeV/sr
2204  G4double xsect = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2205 
2206  fWeightN = xsect*sin(Theta);
2207 
2208  // G4cout << E_in/GeV <<"\t" << pf/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2209  // G4cout << "-- ** Pion crsec :: " << xsect << " ** --" << G4endl;
2210 
2211  return xsect;
2212 }
G4double GetThetaAngle_Min()
G4double wiser_sigma(G4double Ebeam, G4double pf, G4double thf, G4double rad_len, G4int type)
Definition: wiser_pion.h:12
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Pion_PhotoProductionAl ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2216 of file QweakSimEPEvent.cc.

References GetThetaAngle_Min(), myUserInfo, ReactionRegion, QweakSimUserInformation::TargetEntranceWindowThickness, QweakSimUserInformation::TargetExitWindowNippleThickness, QweakSimUserInformation::TargetLength, QweakSimUserInformation::TargetThicknessDSALDummy4, and wiser_sigma().

2221 {
2222  // Pi+- mass - PDG2014
2223  const G4double Mpi = 139.57018*MeV;
2224 
2225  if (Theta<GetThetaAngle_Min())
2226  {
2227  Theta = GetThetaAngle_Min();
2228  }
2229 
2230  //G4double tmp_EpMax = GetEPrime_Max();
2231  //if(GetEPrime_Max() > E_in) tmp_EpMax = E_in;
2232 
2233  //E_out = (G4UniformRand()*(tmp_EpMax - GetEPrime_Min()) + GetEPrime_Min()); //--- final total energy in MeV
2234 
2235  if(E_in < Mpi) E_in = Mpi;
2236  E_out = Mpi + G4UniformRand()*(E_in - Mpi);
2237  G4double pf = sqrt(pow(E_out,2.0) - pow(Mpi,2.0)); //--- final momentum in MeV
2238 
2239  // Calculate Q2
2240  G4double STH = sin(Theta/2.0);
2241  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
2242 
2243  //---
2244  //--- radiation length from page 12 of "The Qweak target design and safety document"
2245  //--- http://qweak.jlab.org/DocDB/0010/001041/002/Qweak%20Target%20PDR.pdf
2246  //---
2247 
2248  // The below was original code from Fang, which is incorrect, I believe (D. Armstrong)
2249  /*
2250  G4double Al_length = 0.0*mm;
2251  if (ReactionRegion == 2 || ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10) //-- Entrance & US Al dummy
2252  Al_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength) + (myUserInfo->TargetEntranceWindowThickness);
2253  else if (ReactionRegion == 3 || ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11) //-- Exit & DS Al dummy
2254  Al_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) - 0.5*(myUserInfo->TargetLength);
2255 
2256  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2257  G4double rad_len = Al_length/(8.896*cm) + intern_rad_len; //--- radiation length
2258  */
2259 
2260  // Need to calculate radiation length upstream of the scattering vertex, in order for Wiser code to estimate bremss spectrum for
2261  // the pion photoproduction. Since the Al windows are relatively thin, we just assume scattering takes place halfway along
2262  // the Al material, on average.
2263 
2264  G4double Al_length = 0.0*mm;
2265  G4double lh2_length = 0.0*mm;
2266  if (ReactionRegion == 2 ) //-- Entrance window
2267  {
2268  Al_length = 0.5*(myUserInfo->TargetEntranceWindowThickness);
2269  }
2270  else if (ReactionRegion == 3 ) //-- Exit window
2271  {
2273  lh2_length = myUserInfo->TargetLength;
2274  }
2275  else if (ReactionRegion == 8 ) //-- 4% DS window No entrance window, no hydrogen here...
2276  {
2277  Al_length = 0.5*(myUserInfo->TargetThicknessDSALDummy4);
2278  }
2279 
2280  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2281  G4double rad_len = Al_length/(8.896*cm) + intern_rad_len + lh2_length/(866.0*cm); //--- radiation length in fractional terms
2282 
2283  //G4cout << "ReactionRegion: " << ReactionRegion << G4endl;
2284  //G4cout << "myPositionZ: " << myPositionZ << G4endl;
2285  //G4cout << "myUserInfo->TargetCenterPositionZ: " << myUserInfo->TargetCenterPositionZ << G4endl;
2286  //G4cout << "myUserInfo->TargetLength: " << myUserInfo->TargetLength << G4endl;
2287  //G4cout << "Al_length: " << Al_length << G4endl;
2288  //G4cout << "lh2_length: " << lh2_length << G4endl;
2289  //G4cout << "rad_len: " << rad_len << G4endl;
2290  //G4cout << "Theta: " << Theta << G4endl;
2291 
2292  G4int type = 1; //--- 0: pi+ 1:pi-
2293 
2294  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2295  // rad_len in fraction (not %)
2296  // wiser_sigma returns crsec in nb/GeV/sr
2297  //
2298  // For a nucleus, need to get cross section from the neutrons and from the protons.
2299  // To generate a pi^- from the proton, need two-pion production (can generate pi^- from neutron via single pion production).
2300  // assume exact isospin symmetry, so to get pi^- cross section for neutron, use pi^+ cross section for proton, which is what
2301  // is in the Wiser code fits.
2302  // Ignore FSI or other nuclear effects (treat 27Al as a bag of 13 free protons and 14 free neutrons)
2303 
2304  G4double xsect_proton = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2305 
2306  // below: switch type to 0 to get the pi+ cross section, which we will use for the pi^- cross section from neutrons
2307  G4double xsect_neutron = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, 0); //--- nanobarns/GeV/str --> ub/GeV/sr
2308  G4double xsect = 13.0*xsect_proton + 14.0*xsect_neutron;
2309  fWeightN = xsect*sin(Theta);
2310 
2311  //G4cout << GetEPrime_Max() << "\t" << GetEPrime_Min() << "\t" << E_in/GeV <<"\t" << E_out/GeV << "\t" << pf/GeV << "\t" << Q2/GeV/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2312  //G4cout << "------ Pion crsec :: " << xsect << " ------" << G4endl;
2313 
2314  return xsect;
2315 }
G4int ReactionRegion
reaction region used for event generation
G4double GetThetaAngle_Min()
G4double wiser_sigma(G4double Ebeam, G4double pf, G4double thf, G4double rad_len, G4int type)
Definition: wiser_pion.h:12
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Pion_PhotoProductionCarbon ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

— GeV

Definition at line 2320 of file QweakSimEPEvent.cc.

References GetEPrime_Max(), GetEPrime_Min(), GetThetaAngle_Min(), myPositionZ, myUserInfo, ReactionRegion, QweakSimUserInformation::TargetCenterPositionZ, QweakSimUserInformation::TargetEntranceWindowThickness, QweakSimUserInformation::TargetLength, and wiser_sigma().

2325 {
2326  const G4double Mpi = 0.13957*GeV; ///--- GeV
2327 
2328  if (Theta<GetThetaAngle_Min())
2329  Theta = GetThetaAngle_Min();
2330 
2331  // need to evaluate Q2 properply, in order to evaluate
2332  // the internal radiation length more accurately below
2333  Q2 = 0.026;
2334 
2335  G4double tmp_EpMax = GetEPrime_Max();
2336  if(GetEPrime_Max() > E_in) tmp_EpMax = E_in;
2337 
2338  E_out = (G4UniformRand()*(tmp_EpMax - GetEPrime_Min()) + GetEPrime_Min()); //--- final total energy in GeV
2339  G4double pf = sqrt(pow(E_out,2) - pow(Mpi,2)); //--- final momentum in GeV
2340 
2341  //---
2342  //--- radiation length calculated from http://pdg.lbl.gov/2013/reviews/rpp2013-rev-atomic-nuclear-prop.pdf
2343  //---
2344  //G4double C_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength);
2345  G4double C_length = 0.0*mm;
2346  if (ReactionRegion == 2 || ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10) //-- Entrance & US C dummy
2348  else if (ReactionRegion == 3 || ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11) //-- Exit & DS C dummy
2350 
2351  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2352  G4double rad_len = C_length/(19.32*cm) + intern_rad_len; //--- radiation length
2353 
2354  G4int type = 1; //--- 0: pi+ 1:pi-
2355 
2356  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2357  // rad_len in fraction (not %)
2358  // wiser_sigma returns crsec in nb/GeV/sr
2359  G4double xsect = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2360 
2361  fWeightN = xsect*sin(Theta);
2362 
2363  //G4cout << E_in/GeV <<"\t" << pf/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2364  //G4cout << "------ Pion crsec :: " << xsect << " ------" << G4endl;
2365 
2366  return xsect;
2367 }
G4int ReactionRegion
reaction region used for event generation
G4double GetThetaAngle_Min()
G4double wiser_sigma(G4double Ebeam, G4double pf, G4double thf, G4double rad_len, G4int type)
Definition: wiser_pion.h:12
G4double GetEPrime_Max()
QweakSimUserInformation * myUserInfo
G4double GetEPrime_Min()

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Quasi_Elastic_Bosted ( G4double  E_in,
G4double  Theta,
G4int  Zin,
G4int  Ain,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 2525 of file QweakSimEPEvent.cc.

References F1F2QE09(), and GetThetaAngle_Min().

2532 {
2533 
2534  G4double F1 = 0.0;
2535  G4double F2 = 0.0;
2536  G4double W1 = 0.0;
2537  G4double W2 = 0.0;
2538  G4double Wsq = 0.0;
2539  G4double xsect = 0.0;
2540  G4double Ein = 0.0;
2541  G4double Eout = 0.0;
2542 
2543  // In this subroutine units are in GeV, so
2544  Ein = E_in/1000.; // maximum energy in GeV
2545  G4double M_electron = 0.000511; // minimum energy (electron mass) in GeV
2546  G4double Mp = 0.93828;
2547 
2548  // G4double Theta_Min = 1.745329E-4;
2549  if (Theta<GetThetaAngle_Min())
2550  Theta = GetThetaAngle_Min();
2551 
2552  // Generate flat energy distribution of outgoing electron
2553  Eout = M_electron + G4UniformRand()*(Ein - M_electron);
2554 
2555  G4double CTH = cos(Theta/2.0);
2556  G4double STH = sin(Theta/2.0);
2557  G4double T2THE = STH*STH/CTH/CTH;
2558  G4double Nu = Ein - Eout;
2559  Q2 = 4.0*Ein*Eout*STH*STH;
2560  Wsq = Mp*Mp + 2.0*Mp*Nu - Q2;
2561 
2562  // Mott scattering
2563  G4double MOTT = pow((0.00072/Ein*CTH/STH/STH),2);
2564  MOTT = MOTT*1.0E4; // Units: ub/sr/GeV
2565 
2566  //G4cout << "Ein, Eout, Nu, Theta: " << Ein << ", " << Eout << ", " << Nu << ", " << Theta/degree << G4endl;
2567  //G4cout << "Q2, W2, Mott: " << Q2 << ", " << Wsq << ", " << MOTT << G4endl;
2568 
2569  F1F2QE09(Zin, Ain, Q2, Wsq, F1, F2);
2570  W1 = F1/Mp;
2571  W2 = F2/Nu;
2572 
2573  xsect = MOTT*(W2 + 2.0*T2THE*W1)*(Ein - M_electron);
2574  // G4cout << "F1, F2, W1, W2: " << F1 << ", " << F2 << ", " << W1 << ", " << W2 << G4endl;
2575 
2576  // In some cases a negative F2 is returned giving a negative cross section
2577  if (xsect <= 0) xsect = 0.0;
2578 
2579  // G4cout << "xsect: " << xsect << G4endl;
2580 
2581  fWeightN = xsect*sin(Theta);
2582 
2583  E_out = Eout*1000.; // Convert back to MeV before leaving
2584  Q2 = Q2*1.0E6; // Seems like this needs to be returned in MeV**2
2585  return xsect;
2586 }
G4double GetThetaAngle_Min()
void F1F2QE09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)

+ Here is the call graph for this function:

G4double QweakSimEPEvent::Quasi_Elastic_Neutron ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 953 of file QweakSimEPEvent.cc.

References M_p, and Theta_Min.

958 {
959 // cout<<"===>>>>Calling Quasi_Elastic_Neutron"<<endl;
960  G4double Theta_Min = 1.745329E-4;
961  G4double Lamda_2 = 0.710;
962  G4double M_p = 938.2796; // proton mass in MeV
963 // G4double mu = 2.793;
964  G4double mu_n = -1.91;
965  G4double Z = 1.0;
966 // G4double A = 1.0;
967  G4double M = M_p;
968 
969 // E_in units is MeV
970 
971  if (Theta<Theta_Min)
972  Theta = Theta_Min;
973 
974  G4double CTH = cos(Theta/2.0);
975  G4double STH = sin(Theta/2.0);
976  G4double T2THE = STH*STH/CTH/CTH;
977  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
978  E_out = E_in/ETA;
979 // MeV^2
980  Q2 = 4.*E_in*E_out*STH*STH;
981  G4double tau = Q2/4./M/M;
982 // Mott scatering
983  G4double CrossSection = (Z*0.72/E_in*CTH/STH/STH)*(Z*0.72/E_in*CTH/STH/STH)/ETA;
984 // Units: ub/sr
985  G4double Mott = CrossSection*10000.0;
986 
987 // Cross section
988  G4double GEP_DIPOLE = 1.0/(1.0+Q2/1.E6/Lamda_2)*(1.0+Q2/1.E6/Lamda_2);
989 // G4double GMP_DIPOLE = GEP_DIPOLE*mu;
990  G4double GEn_DIPOLE = 0.0;
991  G4double GMn_DIPOLE = GEP_DIPOLE*mu_n;
992  G4double FAC = 1.0/(1.0+tau);
993 
994 // G4double Sigma_Dipole_p = Mott*(GEP_DIPOLE*GEP_DIPOLE*FAC+tau*GMP_DIPOLE*GMP_DIPOLE*(FAC+2.*T2THE));
995  G4double Sigma_Dipole_n = Mott*(GEn_DIPOLE*GEn_DIPOLE*FAC+tau*GMn_DIPOLE*GMn_DIPOLE*(FAC+2.*T2THE));
996 // G4double Sigma_Dipole = Z*Sigma_Dipole_p + (A-Z)*Sigma_Dipole_n;
997 
998 // cout<<"===>>>>leaving Quasi_Elastic_Neutron"<<endl;
999  fWeightN = Sigma_Dipole_n*sin(Theta);
1000  return Sigma_Dipole_n;
1001 }
static const G4double M_p
static const G4double Theta_Min
const std::vector< G4double > QweakSimEPEvent::Radiative_Cross_Section_Lookup ( G4double  E_in,
G4double  Theta,
G4double &  fWeightN,
G4double &  Q2,
G4double &  E_out 
)

Definition at line 1432 of file QweakSimEPEvent.cc.

References ElasticPeakDeltaE, fLookupTable, QweakSimUserInformation::GetBeamEnergy(), QweakSimUserInformation::GetEPrime(), QweakSimUserInformation::GetOriginVertexPositionZ(), GetThetaAngle_Min(), myUserInfo, ReactionRegion, SuperElasticCheck(), value_d, and value_n.

1437 {
1438  //G4int SuperElastic = 1;
1439  //G4double xsec = 0.0;
1440  Double_t value[value_n] = {0.0};
1441  std::vector< G4double > CrossSection;
1442  G4double xsec = 0.0;
1443  if (Theta<GetThetaAngle_Min())
1444  Theta = GetThetaAngle_Min();
1445 
1446  //while (SuperElastic) {
1447  //E_out = (G4UniformRand()*(GetEPrime_Max()-GetEPrime_Min()) + GetEPrime_Min());
1448  //SuperElastic = SuperElasticCheck(E_in, E_out, Theta, xsec);
1449  //G4cout << "El XSec: " << xsec << G4endl;
1450  //}
1451 
1452  E_out = myUserInfo->GetEPrime();
1453 
1454  Double_t Zpos;
1455 
1456  if ( ReactionRegion == 1 ) {
1458  }
1459  else {
1460  Zpos = 0.0;
1461  }
1462  const Double_t pos[value_d] = {Zpos,E_in,E_out,Theta};
1463  //Double_t value[value_n] = {0.0};
1464  //std::vector< G4double > CrossSection;
1465 
1466  fLookupTable->GetValue(pos,value);
1467  Q2 = value[1];
1468  fWeightN = value[6]*sin(Theta*degree);
1469  //fWeightN = value[6];
1470  //G4cout << "Rad XSec: " << value[12]*0.001 << G4endl;
1471  //}
1472 
1473  CrossSection.push_back(value[6]*0.001); // 0 total radiated cross section
1474  CrossSection.push_back(value[2]*0.001); // 1 total born cross section
1475  CrossSection.push_back(value[3]*0.001); // 2 inelastic born cross section
1476  CrossSection.push_back(value[4]*0.001); // 3 quasi-elastic born cross section
1477  CrossSection.push_back(value[6]*0.001); // 4 total radiated cross section
1478  CrossSection.push_back(value[7]*0.001); // 5 elastic radiated cross section
1479  CrossSection.push_back(value[8]*0.001); // 6 quasi-elastic radiated cross section
1480  CrossSection.push_back(value[9]*0.001); // 7 deep-inelastic radiated cross section
1481  CrossSection.push_back(value[11]*0.001); // 8 total radiated cross section internal only
1482  CrossSection.push_back(value[12]*0.001); // 9 elastic radiated cross section internal only
1483  CrossSection.push_back(value[14]*0.001); // 10 quasi-elastic radiated cross section internal only
1484  CrossSection.push_back(value[13]*0.001); // 11 deep-inelastic radiated cross section internal only
1485  CrossSection.push_back(0.0); // 12 elastic peak cross section
1486 
1487  if (SuperElasticCheck(myUserInfo->GetBeamEnergy(), E_out, Theta, xsec)) CrossSection[5] = 0.0;
1488  else if (SuperElasticCheck(myUserInfo->GetBeamEnergy(), E_out+ElasticPeakDeltaE, Theta, xsec)) {
1489  CrossSection[12] = xsec/(ElasticPeakDeltaE*0.001);
1490  CrossSection[5] = 0.0;
1491  }
1492  CrossSection[4] = CrossSection[5] + CrossSection[6] + CrossSection[7] + CrossSection[12];
1493  CrossSection[0] = CrossSection[4];
1494 
1495  return CrossSection; // units [ub/sr/GeV]
1496 }
G4int ReactionRegion
reaction region used for event generation
QweakSimFieldMap< value_t, value_n > * fLookupTable
G4double GetThetaAngle_Min()
G4double GetOriginVertexPositionZ() const
static const G4int value_n
G4int SuperElasticCheck(G4double E_in, G4double E_out, G4double theta, G4double &xsec)
G4double ElasticPeakDeltaE
static const G4int value_d
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

G4double QweakSimEPEvent::ResMod507 ( G4int  sf,
G4double  w2,
G4double  q2,
G4double *  xval 
)
private

Definition at line 1132 of file QweakSimEPEvent.cc.

Referenced by Sigma_EEPrime().

1133 {
1134  G4double xb,mass[8],width[8];
1135  G4double height[8],rescoef[7][5];
1136  G4double nr_coef[4][5],sigr[8],sig_nr;
1137  G4double intwidth[8],k,kcm,kr[8],kcmr[8],ppicm,ppi2cm;
1138  G4double petacm,ppicmr[8],ppi2cmr[8],petacmr[8],epicmr[8],epi2cmr[8];
1139  G4double eetacmr[8],epicm,epi2cm,eetacm,br[8][4],ang[8];
1140  G4double pgam[8],pwid[8][4],x0[8],q20;
1141  G4double sig_res,t,xpr[3],m0;
1142  G4int i,j,num;
1143 
1144  G4double xt,p1,p2,p3;
1145 
1146  G4double mp = 0.9382727;
1147  G4double mpi = 0.135;
1148  G4double meta = 0.547;
1149  G4double mp2 = mp*mp;
1150  G4double w = sqrt(w2);
1151  //wdif[1] = w - (mp + mpi);
1152  //wdif[2] = w - (mp + 2.*mpi);
1153 
1154  m0 = 0.125;
1155  if(sf==2)
1156  m0 = xval[49];
1157 
1158  if(sf==1)
1159  {
1160  q20 = 0.05;
1161  }
1162  else
1163  {
1164  q20 = 0.300;
1165  }
1166 
1167 // single pion branching ratios
1168  br[1][1] = 1.0; // P33(1232)
1169  br[2][1] = 0.45; // S11(1535)
1170  br[3][1] = 0.65; // D13(1520)
1171  br[4][1] = 0.65; // F15(1680)
1172  br[5][1] = 0.4; // S11(1650)
1173  br[6][1] = 0.65; // P11(1440) roper
1174  br[7][1] = 0.50; // F37(1950)
1175 
1176 // eta branching ratios
1177  br[1][3] = 0.0; // !!! P33(1232)
1178  br[2][3] = 0.45; // !!! S11(1535)
1179  br[3][3] = 0.0; // !!! D13(1520)
1180  br[4][3] = 0.0; // !!! F15(1680)
1181  br[5][3] = 0.1; // !!! S11(1650)
1182  br[6][3] = 0.0; // !!! P11(1440) roper
1183  br[7][3] = 0.0; // !!! F37(1950)
1184 
1185 // 2-pion branching ratios
1186  for (i=1;i<=7;i++)
1187  {
1188  br[i][2] = 1.-br[i][1]-br[i][3];
1189  }
1190 
1191 // Meson angular momentum
1192  ang[1] = 1.; // P33(1232)
1193  ang[2] = 0.; // S11(1535)
1194  ang[3] = 2.; // D13(1520)
1195  ang[4] = 3.; // F15(1680)
1196  ang[5] = 0.; // S15(1650)
1197  ang[6] = 1.; // P11(1440) roper
1198  ang[7] = 3.; // F37(1950)
1199 
1200  for( i=1;i<=7;i++) // resonance damping parameter
1201  {
1202  x0[i] = 0.145;
1203  }
1204 
1205  x0[1] = 0.130;
1206 
1207  for (i=1;i<=7;i++)
1208  {
1209  br[i][2] = 1.-br[i][1]-br[i][3];
1210  }
1211 
1212  // dip = 1./((1.+q2/0.71)*(1.+q2/0.71)); // Dipole parameterization
1213  // mon = 1./(1.+q2/0.71);
1214 
1215  xb = q2/(q2+w2-mp2);
1216  xpr[1] = 1.+(w2-(mp+mpi)*(mp+mpi))/(q2+q20);
1217  xpr[1] = 1./xpr[1];
1218  xpr[2] = 1.+(w2-(mp+mpi+mpi)*(mp+mpi+mpi))/(q2+q20);
1219  xpr[2] = 1./xpr[2];
1220 
1221  t = log(log((q2+m0)/0.330*0.330)/log(m0/0.330*0.330));
1222 
1223 // Calculate kinematics needed for threshold Relativistic B-W
1224 
1225  k = (w2 - mp2)/2./mp;
1226  kcm = (w2-mp2)/2./w;
1227 
1228  epicm = (w2 + mpi*mpi -mp2 )/2./w;
1229  ppicm = sqrt(fmax(0.0,(epicm*epicm - mpi*mpi)));
1230  epi2cm = (w2 + (2.*mpi)*(2.*mpi) -mp2 )/2./w;
1231  ppi2cm = sqrt(fmax(0.0,(epi2cm*epi2cm - (2.*mpi)*(2.*mpi))));
1232  eetacm = (w2 + meta*meta -mp2 )/2./w;
1233  petacm = sqrt(fmax(0.0,(eetacm*eetacm - meta*meta)));
1234 
1235  num = 0;
1236 
1237  for (i=1;i<=6;i++) // Read in resonance masses
1238  {
1239  num = num + 1;
1240  mass[i] = xval[i];
1241  }
1242 
1243  for (i=1;i<=6;i++) // Read in resonance widths
1244  {
1245  num = num + 1;
1246  intwidth[i] = xval[num];
1247  width[i] = intwidth[i];
1248  }
1249 
1250  if(sf==2) // Put in 4th resonance region
1251  {
1252  mass[7] = xval[43];
1253  intwidth[7] = xval[44];
1254  width[7] = intwidth[7];
1255  }
1256  else
1257  {
1258  mass[7] = xval[47];
1259  intwidth[7] = xval[48];
1260  width[7] = intwidth[7];
1261  }
1262 
1263  for (i=1;i<=7;i++)
1264  {
1265  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
1266  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
1267  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
1268  ppicmr[i] = sqrt(fmax(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)));
1269  epi2cmr[i] = (mass[i]*mass[i] + (2.*mpi)*(2.*mpi) -mp2 )/2./mass[i];
1270  ppi2cmr[i] = sqrt(fmax(0.0,(epi2cmr[i]*epi2cmr[i] - (2.*mpi)*(2.*mpi))));
1271  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
1272  petacmr[i] = sqrt(fmax(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)));
1273 
1274  // Calculate partial widths
1275  pwid[i][1] = intwidth[i]*pow(ppicm/ppicmr[i],(2.*ang[i]+1.))
1276  *pow((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i]),ang[i]);
1277  // 1-pion decay mode
1278 
1279  pwid[i][2] = intwidth[i]*pow(ppi2cm/ppi2cmr[i],(2.*ang[i]+4.))
1280  *pow((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i]),(ang[i]+2));
1281  // 2-pion decay mode
1282 
1283  pwid[i][2] = w/mass[i]*pwid[i][2];
1284  pwid[i][3] = 0.; // !!! eta decay mode
1285 
1286  if(i==2 || i==5)
1287  {
1288  pwid[i][3]=intwidth[i]*pow(petacm/petacmr[i],(2.*ang[i]+1.))
1289  *pow((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i]),ang[i]);
1290  // eta decay only for S11's
1291  }
1292 
1293  pgam[i] = (kcm/kcmr[i])*(kcm/kcmr[i])*(kcmr[i]*kcmr[i]+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
1294  pgam[i] = intwidth[i]*pgam[i];
1295  width[i] = br[i][1]*pwid[i][1]+br[i][2]*pwid[i][2]+br[i][3]*pwid[i][3];
1296  }
1297  // End resonance kinematics and Widths calculations
1298 
1299  // Begin resonance Q^2 dependence calculations
1300  for (i=1;i<=6;i++)
1301  {
1302  for (j=1;j<=4;j++)
1303  {
1304  num = num + 1;
1305  rescoef[i][j]=xval[num];
1306  }
1307 
1308  if(sf==1)
1309  {
1310  height[i] = rescoef[i][1]*(1.+rescoef[i][2]*q2/(1.+rescoef[i][3]*q2))/ pow(1.+q2/0.91,rescoef[i][4]);
1311  }
1312  else
1313  {
1314  height[i] = rescoef[i][1]*q2/(1.+rescoef[i][2]*q2)*exp(-1.*rescoef[i][3]*q2);
1315  }
1316 
1317  height[i] = height[i]*height[i];
1318  }
1319 
1320 
1321  if(sf==2) // 4th resonance region
1322  {
1323  height[7] = xval[45]*q2/(1.+xval[46]*q2)*exp(-1.*xval[47]*q2);
1324  }
1325  else
1326  {
1327  height[7] = xval[49]/(1.+q2/0.91);
1328  }
1329 
1330  height[7] = height[7]*height[7];
1331 
1332 // End resonance Q^2 dependence calculations
1333 
1334  for (i=1;i<=3;i++) // Non-Res coefficients
1335  {
1336  for (j=1;j<=4;j++)
1337  {
1338  num = num + 1;
1339  nr_coef[i][j]=xval[num];
1340  }
1341  }
1342 
1343 // Calculate Breit-Wigners for all resonances
1344 
1345  sig_res = 0.0;
1346  for (i=1;i<=7;i++)
1347  {
1348  sigr[i] = width[i]*pgam[i]/((w2 - mass[i]*mass[i])*(w2 - mass[i]*mass[i])
1349  + (mass[i]*width[i])*(mass[i]*width[i]));
1350  sigr[i] = height[i]*kr[i]/k*kcmr[i]/kcm*sigr[i]/intwidth[i];
1351  sig_res = sig_res + sigr[i];
1352  }
1353  sig_res = sig_res*w;
1354 
1355 // Finish resonances / start non-res background calculation
1356 
1357  sig_nr = 0.;
1358  if(sf==1)
1359  {
1360  xt = (q2+xval[37])/(w2-1.15+q2+xval[37]);
1361  p2 = xval[38]*exp(-1.*xval[39]*q2)+xval[40]*log(q2/xval[41]+1.);
1362  p3 = xval[42]*exp(-1.*xval[43]*q2)+xval[44]*log(q2/xval[45]+1.);
1363  p1 = xval[46]*(1.+xval[50]*q2/(1.+0.534*q2));
1364 
1365  sig_nr = p1*pow(xt,p2)*pow(1.-xt,p3);
1366  }
1367  else if(sf==2)
1368  {
1369 
1370  for (i=1;i<=1;i++)
1371  {
1372  sig_nr = sig_nr + nr_coef[i][1]*pow(1.-xpr[i],(nr_coef[i][3]+nr_coef[i][2]*t)) /(1.-xb)/(q2+q20)
1373  *pow(q2/(q2+q20), nr_coef[i][4]) *pow(xpr[i],(xval[41]+xval[42]*t));
1374  }
1375  }
1376 
1377  G4double sig = sig_res + sig_nr;
1378  return sig;
1379 }

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::resmod507_v2 ( G4double  sf,
G4double  w2,
G4double  q2,
G4double  xval[50] 
)

!! P33(1232)

!! S11(1535)

!! D13(1520)

!! F15(1680)

!! S11(1650)

!! P11(1440) roper

!! F37(1950)

!! P33(1232)

!! S11(1535)

!! D13(1520)

!! F15(1680)

!! S11(1650)

!! P11(1440) roper

!! F37(1950)

!! P33(1232)

!! S11(1535)

!! D13(1520)

!! F15(1680)

!! S15(1650)

!! P11(1440) roper

!! F37(1950)

!! resonance damping parameter !!!

Definition at line 3373 of file QweakSimEPEvent.cc.

Referenced by christy507().

3375 {
3376 
3377  G4double W,mp,mp2,xb;
3378  G4double mass[7] = {0,0,0,0,0,0,0};
3379  G4double width[7] = {0,0,0,0,0,0,0};
3380  G4double height[7] = {0,0,0,0,0,0,0};
3381  G4double intwidth[7] = {0,0,0,0,0,0,0};
3382  G4double rescoef[6][4];
3383  G4double nr_coef[3][4],sigr[7],wdif[2],sig_nr,h_nr[3];
3384  G4double mpi,meta,k,kcm,kr[7],kcmr[7],ppicm,ppi2cm;
3385  G4double petacm,ppicmr[7],ppi2cmr[7],petacmr[7],epicmr[7],epi2cmr[7];
3386  G4double eetacmr[7],epicm,epi2cm,eetacm,br[7][3],ang[7];
3387  G4double pgam[7],pwid[7][3],x0[7],q20;
3388  G4double sig_res,t,xpr[2],m0,sig;
3389  G4int i,j,num;
3390 
3391  mp = 0.9382727;
3392  mpi = 0.135;
3393  meta = 0.547;
3394  mp2 = mp*mp;
3395  W = sqrt(w2);
3396  wdif[0] = W - (mp + mpi);
3397  wdif[1] = W - (mp + 2.*mpi);
3398 
3399  m0 = 0.125;
3400  if(sf == 2) m0 = xval[48];
3401  if(sf == 1) {
3402  q20 = 0.05;
3403  } else {
3404  q20 = 0.125;
3405  }
3406 
3407  // single pion branching ratios CCCC
3408 
3409  br[0][0] = 1.0; //!!! P33(1232)
3410  br[1][0] = 0.45; //!!! S11(1535)
3411  br[2][0] = 0.65; //!!! D13(1520)
3412  br[3][0] = 0.65; //!!! F15(1680)
3413  br[4][0] = 0.4; //!!! S11(1650)
3414  br[5][0] = 0.65; //!!! P11(1440) roper
3415  br[6][0] = 0.50 ; //!!! F37(1950)
3416 
3417  br[0][2] = 0.0; //!!! P33(1232)
3418  br[1][2] = 0.45; //!!! S11(1535)
3419  br[2][2] = 0.0; //!!! D13(1520)
3420  br[3][2] = 0.0; //!!! F15(1680)
3421  br[4][2] = 0.1; //!!! S11(1650)
3422  br[5][2] = 0.0; //!!! P11(1440) roper
3423  br[6][2] = 0.0; //!!! F37(1950)
3424 
3425  // 2-pion branching ratios CCCC
3426 
3427  for (i = 0; i < 7; i++) {
3428  br[i][1] = 1.-br[i][0]-br[i][2];
3429  }
3430 
3431  // Meson angular momentum CCCC
3432  ang[0] = 1.; //!!! P33(1232)
3433  ang[1] = 0.; //!!! S11(1535)
3434  ang[2] = 2.; //!!! D13(1520)
3435  ang[3] = 3.; //!!! F15(1680)
3436  ang[4] = 0.; //!!! S15(1650)
3437  ang[5] = 1.; //!!! P11(1440) roper
3438  ang[6] = 3.; //!!! F37(1950)
3439 
3440  for (i = 0; i < 7; i++) { //!!! resonance damping parameter !!!
3441  x0[i] = 0.215;
3442  }
3443  x0[0] = 0.15;
3444  x0[0] = xval[49];
3445 
3446  xb = q2/(q2+w2-mp2);
3447  xpr[0] = 1.+(w2-pow((mp+mpi),2))/(q2+q20);
3448  xpr[0] = 1./xpr[0];
3449  xpr[1] = 1.+(w2-pow((mp+mpi+mpi),2))/(q2+q20);
3450  xpr[1] = 1./xpr[1];
3451 
3452  t = log(log((q2+m0)/0.330/0.330)/log(m0/0.330/0.330));
3453 
3454  // Calculate kinematics needed for threshold Relativistic B-W
3455 
3456  k = (w2 - mp2)/2./mp;
3457  kcm = (w2-mp2)/2./W;
3458 
3459  epicm = (w2 + mpi*mpi -mp2 )/2./W;
3460  ppicm = pow(std::max(0.0,(epicm*epicm - mpi*mpi)),0.5);
3461  epi2cm = (w2 + pow((2.*mpi),2) -mp2 )/2./W;
3462  ppi2cm = pow(std::max(0.0,(epi2cm*epi2cm - pow((2.*mpi),2))),0.5);
3463  eetacm = (w2 + meta*meta -mp2 )/2./W;
3464  petacm = pow(std::max(0.0,(eetacm*eetacm - meta*meta)),0.5);
3465 
3466  num = 0;
3467 
3468  for (i = 0; i < 6; i++) { // !!! Read in resonance masses !!!
3469  num = num + 1;
3470  mass[i] = xval[i];
3471  }
3472  for (i = 0; i < 6; i++) { // !!! Read in resonance widths !!!
3473  num = num + 1;
3474  intwidth[i] = xval[num-1];
3475  width[i] = intwidth[i];
3476  }
3477  if (sf == 2) { // !!! Put in 4th resonance region !!!
3478  mass[6] = xval[42];
3479  intwidth[6] = xval[43];
3480  width[6] = intwidth[6];
3481  } else {
3482  mass[6] = xval[46];
3483  intwidth[6] = xval[47];
3484  width[6] = intwidth[6];
3485  }
3486 
3487  for (i = 0; i < 7; i++) {
3488  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
3489  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
3490  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
3491  ppicmr[i] = pow(std::max(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)),0.5);
3492  epi2cmr[i] = (mass[i]*mass[i] + pow((2.*mpi),2) -mp2 )/2./mass[i];
3493  ppi2cmr[i] = pow(std::max(0.0,(epi2cmr[i]*epi2cmr[i] - pow((2.*mpi),2))),0.5);
3494  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
3495  petacmr[i] = pow(std::max(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)),0.5);
3496 
3497  // CCC Calculate partial widths CCC
3498 
3499  pwid[i][0] = intwidth[i]*pow((ppicm/ppicmr[i]),(2.*ang[i]+1.))*pow(((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i])),ang[i]); // !!! 1-pion decay mode
3500 
3501  pwid[i][1] = intwidth[i]*pow((ppi2cm/ppi2cmr[i]),(2.*ang[i]+4.))*pow(((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i])),(ang[i]+2)); // !!! 2-pion decay mode
3502 
3503  pwid[i][1] = W/mass[i]*pwid[i][1];
3504  pwid[i][2] = 0.; // !!! eta decay mode
3505 
3506  if(i == 1 || i == 4) {
3507  pwid[i][2] = intwidth[i]*pow((petacm/petacmr[i]),(2.*ang[i]+1.))*pow(((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i])),ang[i]); // !!! eta decay only for S11's
3508  }
3509 
3510  pgam[i] = pow((kcm/kcmr[i]),2)*(pow(kcmr[i],2)+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
3511  pgam[i] = intwidth[i]*pgam[i];
3512  width[i] = br[i][0]*pwid[i][0]+br[i][1]*pwid[i][1]+br[i][2]*pwid[i][2];
3513  }
3514 
3515  //CCC End resonance kinematics and Widths calculations CCC
3516 
3517  // CCC Begin resonance Q^2 dependence calculations CCC
3518 
3519  for (i = 0; i < 6; i++) {
3520  for (j = 0; j < 4; j++) {
3521  num = num + 1;
3522  rescoef[i][j]=xval[num-1];
3523  }
3524 
3525  if(sf == 1) {
3526  height[i] = rescoef[i][0]*(1.+rescoef[i][1]*q2/(1.+rescoef[i][2]*q2))/pow((1.+q2/0.91),rescoef[i][3]);
3527  } else {
3528  height[i] = rescoef[i][0]*q2/(1.+rescoef[i][1]*q2)*exp(-1.*rescoef[i][2]*q2);
3529  }
3530 
3531  height[i] = height[i]*height[i];
3532  }
3533 
3534  if(sf == 2) { // !!! 4th resonance region !!!
3535  height[6] = xval[44]*q2/(1.+xval[45]*q2)*exp(-1.*xval[46]*q2);
3536  } else {
3537  height[6] = xval[48]/(1.+q2/0.91);
3538  }
3539 
3540  height[6] = height[6]*height[6];
3541 
3542  // CCC End resonance Q^2 dependence calculations CCC
3543 
3544  for (i = 0; i < 3; i++) { // !!! Non-Res coefficients !!!
3545  for (j = 0; j < 4; j++) {
3546  num = num + 1;
3547  nr_coef[i][j]=xval[num-1];
3548  }
3549  }
3550 
3551  // CCC Calculate Breit-Wigners for all resonances CCC
3552 
3553  sig_res = 0.0;
3554  for (i = 0; i < 7; i++) {
3555  sigr[i] = width[i]*pgam[i]/(pow(w2 - mass[i]*mass[i],2.) + pow(mass[i]*width[i],2.));
3556  sigr[i] = height[i]*kr[i]/k*kcmr[i]/kcm*sigr[i]/intwidth[i];
3557  sig_res = sig_res + sigr[i];
3558  }
3559  sig_res = sig_res*W;
3560  // CCC Finish resonances / start non-res background calculation CCC
3561  sig_nr = 0.;
3562  if(sf == 1) {
3563  for (i = 0; i < 2; i++) {
3564  h_nr[i] = nr_coef[i][0]/pow(q2+nr_coef[i][1],nr_coef[i][2] + nr_coef[i][3]*q2+xval[44+i]*q2*q2);
3565  sig_nr = sig_nr +h_nr[i]*(pow(wdif[0],(2*(i+1)+1)/2.));
3566  }
3567 
3568  sig_nr = sig_nr*xpr[0];
3569 
3570  } else if (sf == 2) {
3571  for (i = 0; i < 1; i++) {
3572  sig_nr = sig_nr + nr_coef[i][0]*pow((1.-xpr[i]),(nr_coef[i][2]+nr_coef[i][1]*t))/(1.-xb)/(q2+q20);
3573  sig_nr = sig_nr*pow(q2/(q2+q20),nr_coef[i][3])*pow(xpr[i],(xval[40]+xval[41]*t));
3574  }
3575  }
3576  sig = sig_res + sig_nr;
3577 
3578  return sig;
3579 }

+ Here is the caller graph for this function:

void QweakSimEPEvent::resmodd ( G4double  w2,
G4double  q2,
G4double  xval[50],
G4double &  sig 
)

returns F1 for average of free proton and neutron for given W2, Q2

changed to allow delta width, mass to vary taken out again since xval[1] used in MEC

resonance height coefficients. xvals of 13-36

Begin non-resonant part uses xvals 45, 46, 50 Depends on both W2 and Q2 so can't easily precalculate

Definition at line 3193 of file QweakSimEPEvent.cc.

Referenced by F1F2IN09().

3195 {
3196  //! returns F1 for average of free proton and neutron
3197  //! for given W2, Q2
3198 
3199  G4double W,mp,mp2,mass[7],width[7];
3200  G4double height[7],rescoef[6][4];
3201  G4double nr_coef[3][4],wdif,sig_nr;
3202  G4double sigr[7];
3203  G4double mpi,meta,intwidth[7],k,kcm,kr[7],kcmr[7],ppicm,ppi2cm;
3204  G4double petacm,ppicmr[7],ppi2cmr[7],petacmr[7],epicmr[7],epi2cmr[7];
3205  G4double eetacmr[7],epicm,epi2cm,eetacm,br[7][2],ang[7],sigrsv[7];
3206  G4double pgam[7],pwid[7][2],x0[7],dip;
3207  G4double sig_res,xpr,alpha,F1;
3208  G4int i,j,num,iw;
3209 
3210  G4double xval0[12] = {
3211  0.12298E+01,0.15304E+01,0.15057E+01,0.16980E+01,0.16650E+01,
3212  0.14333E+01,0.13573E+00,0.22000E+00,0.82956E-01,0.95782E-01,
3213  0.10936E+00,0.37944E+00
3214  };
3215 
3216  G4double w2sv;
3217 
3218  mp = 0.9382727;
3219  mpi = 0.135;
3220  meta = 0.547;
3221  mp2 = mp*mp;
3222  alpha = 1./137.036;
3223 
3224  sig = 0.;
3225 
3226  if(w2 < 1.07327*1.07327 || w2 > 25 || q2 < 0.0 || q2 > 11.0) {
3227  G4cerr << "ERROR: " << __FILE__ << " line " << __LINE__ << G4endl;
3228  G4cerr << " W2/Q2 check failed: should be 1.07327^2<W2<25 and 0 < Q2 < 11" << G4endl;
3229  G4cerr << " but are w2 q2: " << w2 << " " << q2 << G4endl;
3230  exit(1);
3231  }
3232 
3233  // ! branching ratios
3234  br[0][0] = 1.0;
3235  br[1][0] = 0.5;
3236  br[2][0] = 0.65;
3237  br[3][0] = 0.65;
3238  br[4][0] = 0.4;
3239  br[5][0] = 0.65;
3240  br[6][0] = 0.6;
3241 
3242  // ! angular momenta
3243  ang[0] = 1.; // !!! P33(1232)
3244  ang[1] = 0.; // !!! S11(1535)
3245  ang[2] = 2.; // !!! D13(1520)
3246  ang[3] = 3.; // !!! F15(1680)
3247  ang[4] = 0.; // !!! S15(1650)
3248  ang[5] = 1.; // !!! P11(1440) roper
3249  ang[6] = 3.; // !!! ? 4th resonance region
3250 
3251  // ! x0 parameter
3252  for (i = 0; i < 7; i++) x0[i] = 0.215;
3253  x0[0] = 0.1446;
3254 
3255  // ! out branching ratio
3256  for (i = 0; i < 7; i++) br[i][1] = 1.-br[i][0];
3257 
3258  // ! remember w2
3259  w2sv = w2;
3260 
3261  // ! uses xvals of 1-12, 47, and 48
3262  // ! move masses, wdiths into local variables
3263  // ! pyb changed to be fixed
3264  num = 0;
3265  for (i = 0; i < 6; i++) {
3266  num = num + 1;
3267  mass[i] = xval0[i];
3268  }
3269  for (i = 0; i < 6; i++) {
3270  num = num + 1;
3271  intwidth[i] = xval0[num-1];
3272  }
3273 
3274  //! changed to allow delta width, mass to vary
3275  //! taken out again since xval[1] used in MEC
3276  mass[6] = 1.9341;
3277  intwidth[6] = 0.380;
3278 
3279  iw = G4int(1000.*sqrt(w2));
3280  W = 0.001 * (iw+0.5);
3281  w2 = W*W;
3282  wdif = W - (mp + mpi);
3283 
3284  // ! Calculate kinematics needed for threshold Relativistic B-W
3285  k = (w2 - mp2) / 2. / mp;
3286  kcm = (w2 - mp2) / 2. / W;
3287  epicm = (w2 + mpi*mpi -mp2 ) / 2. / W;
3288  ppicm = pow(std::max(0.0,(epicm*epicm - mpi*mpi)),0.5);
3289  epi2cm = (w2 + pow((2.*mpi),2) -mp2 ) / 2. / W;
3290  ppi2cm = pow(std::max(0.0,(epi2cm*epi2cm - pow((2.*mpi),2))),0.5);
3291  eetacm = (w2 + meta*meta -mp2 ) / 2. / W;
3292  petacm = pow(std::max(0.0,(eetacm*eetacm - meta*meta)),0.5);
3293  for (i = 0; i < 7; i++) {
3294  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
3295  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
3296  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
3297  ppicmr[i] = pow(std::max(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)),0.5);
3298  epi2cmr[i] = (mass[i]*mass[i] + pow(2.*mpi,2) -mp2 )/2./mass[i];
3299  ppi2cmr[i] = pow(std::max(0.0,(epi2cmr[i]*epi2cmr[i] - pow((2.*mpi),2))),0.5);
3300  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
3301  petacmr[i] = pow(std::max(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)),0.5);
3302 
3303  // ! Calculate partial widths
3304  pwid[i][0] = intwidth[i]*pow((ppicm/ppicmr[i]),(2.*ang[i]+1.))*pow(((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i])),ang[i]);
3305 
3306  if(i != 1)
3307  pwid[i][1] = intwidth[i]*pow((ppi2cm/ppi2cmr[i]),(2.*ang[i]+4.))*
3308  pow(((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i])),(ang[i]+2.))* W / mass[i];
3309  else
3310  pwid[i][1] = intwidth[1]*pow((petacm/petacmr[i]),(2.*ang[i]+1.))*
3311  pow(((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i])),ang[i]);
3312 
3313  pgam[i] = pow((kcm/kcmr[i]),2)*(kcmr[i]*kcmr[i]+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
3314  pgam[i] = intwidth[i]*pgam[i];
3315  width[i] = br[i][0]*pwid[i][0]+br[i][1]*pwid[i][1];
3316  sigr[i] = width[i] * pgam[i] / (pow(w2 - mass[i]*mass[i],2) + pow(mass[i]*width[i],2)) * kr[i] / k * kcmr[i] / kcm / intwidth[i];
3317  }
3318 
3319  w2 = w2sv;
3320 
3321  // ! get parameters into local variables
3322  num = 12;
3323  //! resonance height coefficients. xvals of 13-36
3324  for (i = 0; i < 6; i++) {
3325  for (j = 0; j < 4; j++) {
3326  num = num + 1;
3327  rescoef[i][j]=xval[num-1];
3328  }
3329  }
3330  // ! Non-Res coefficients xvals of 37-44
3331  for (i = 0; i < 2; i++) {
3332  for (j = 0; j < 4; j++) {
3333  num = num + 1;
3334  nr_coef[i][j]=xval[num-1];
3335  }
3336  }
3337 
3338  // ! Begin resonance Q^2 dependence calculations CCC
3339  // ! uses xvals 49
3340  for (i = 0; i < 6; i++) {
3341  height[i] = rescoef[i][0]*(1.+rescoef[i][1] * q2 / (1. + rescoef[i][2] * q2))/pow((1. + q2/0.91),rescoef[i][3]);
3342  }
3343  dip = 1./pow((1. + q2 / 0.91),2);
3344  height[6] = xval[48]*dip;
3345  sig_res = 0.;
3346  for (i = 0; i < 7; i++) {
3347  sigrsv[i] = height[i]*height[i] * sigr[i];
3348  sig_res = sig_res + sigrsv[i];
3349  }
3350 
3351  sig_res = sig_res * sqrt(w2);
3352 
3353  //! Begin non-resonant part uses xvals 45, 46, 50
3354  //! Depends on both W2 and Q2 so can't easily precalculate
3355  sig_nr = 0.;
3356  xpr = 1.+(w2-pow(mp+mpi,2))/(q2+0.05);
3357  xpr = 1./xpr;
3358  W = pow(w2,0.5);
3359  wdif = W - (mp + mpi);
3360  for (i = 0; i < 2; i++) {
3361  sig_nr = sig_nr +(nr_coef[i][0]*pow(wdif,((2*(i+1)+1)/2.)))/pow(q2+nr_coef[i][1],(nr_coef[i][2]+nr_coef[i][3]*q2+xval[44+i]*q2*q2));
3362  }
3363  sig_nr = sig_nr * xpr;
3364  sig = sig_res + sig_nr;
3365 
3366  F1 = sig * (w2-mp2)/8./pi/pi/alpha/0.3894e3;
3367  sig = F1;
3368 }

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetActiveOctantNumber ( G4int  kaot)
inline

Definition at line 96 of file QweakSimEPEvent.hh.

References kActiveOctantNumber.

Referenced by QweakSimEPEventMessenger::SetNewValue().

96 { kActiveOctantNumber = kaot; };
G4int kActiveOctantNumber
Active octant number in the simulation, 0 will enable all octants.

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetBeamEnergy ( G4double  energy = 1.160*GeV)

Definition at line 4573 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), myUserInfo, and QweakSimUserInformation::SetBeamEnergy().

Referenced by CheckLookupTableBounds(), CreateLookupTable(), QweakSimPrimaryGeneratorAction::QweakSimPrimaryGeneratorAction(), and QweakSimEPEventMessenger::SetNewValue().

4573  {
4574  if (energy>0) {
4576  myUserInfo->SetBeamEnergy(energy);
4577  }
4578  else {
4579  G4cout << G4endl << "##### Beam Energy must be greater than zero" << G4endl << G4endl;
4580  }
4581 }
void SetBeamEnergy(G4double energy)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetElasticPeakDeltaE ( G4double  energy = 15*MeV)
inline

Definition at line 162 of file QweakSimEPEvent.hh.

References ElasticPeakDeltaE.

Referenced by QweakSimEPEventMessenger::SetNewValue().

162 { if (energy < 15*MeV) { energy = 15*MeV; G4cout << "***** Value Out of Range: ElasticPeakDeltaE set to : 15 MeV *****" << G4endl; } ElasticPeakDeltaE = energy; } ;
G4double ElasticPeakDeltaE

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetEPrime_Max ( G4double  energy)

Definition at line 4552 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), myUserInfo, and QweakSimUserInformation::SetEPrime_Max().

Referenced by CheckLookupTableBounds(), and QweakSimEPEventMessenger::SetNewValue().

void SetEPrime_Max(G4double energy)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetEPrime_Min ( G4double  energy)

Definition at line 4549 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), myUserInfo, and QweakSimUserInformation::SetEPrime_Min().

Referenced by CheckLookupTableBounds(), and QweakSimEPEventMessenger::SetNewValue().

void SetEPrime_Min(G4double energy)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetIsotropy ( G4int  isot)
inline

Definition at line 99 of file QweakSimEPEvent.hh.

References Isotropy.

Referenced by QweakSimEPEventMessenger::SetNewValue().

99 { Isotropy = isot; };
G4int Isotropy
isotropy used for event generation

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetLuminosity ( G4double  lum)

Definition at line 4567 of file QweakSimEPEvent.cc.

References myUserInfo, and QweakSimUserInformation::SetLuminosity().

4567 {myUserInfo->SetLuminosity(lum);}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

void QweakSimEPEvent::SetPhaseSpace ( G4double  ps)

Definition at line 4570 of file QweakSimEPEvent.cc.

References myUserInfo, and QweakSimUserInformation::SetPhaseSpace().

4570 {myUserInfo->SetPhaseSpace(ps);}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

void QweakSimEPEvent::SetPhiAngle_Max ( G4double  ang)

Definition at line 4564 of file QweakSimEPEvent.cc.

References myUserInfo, and QweakSimUserInformation::SetPhiAngle_Max().

Referenced by QweakSimEPEventMessenger::SetNewValue().

4564 {myUserInfo->SetPhiAngle_Max(ang);}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetPhiAngle_Min ( G4double  ang)

Definition at line 4561 of file QweakSimEPEvent.cc.

References myUserInfo, and QweakSimUserInformation::SetPhiAngle_Min().

Referenced by QweakSimEPEventMessenger::SetNewValue().

4561 {myUserInfo->SetPhiAngle_Min(ang);}
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetReactionRegion ( G4int  rr)
inline

Definition at line 174 of file QweakSimEPEvent.hh.

References CreateLookupTable(), ReactionRegion, and ReactionType.

Referenced by QweakSimEPEventMessenger::SetNewValue().

G4int ReactionRegion
reaction region used for event generation
G4int ReactionType
reaction type used for event generation

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetReactionType ( G4int  rt)
inline

Definition at line 171 of file QweakSimEPEvent.hh.

References ReactionType, and TypeSetting.

Referenced by QweakSimEPEventMessenger::SetNewValue().

171 { ReactionType = rt; TypeSetting = rt;};
G4int ReactionType
reaction type used for event generation

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetSchwingerDeltaE ( G4double  energy = 15*MeV)
inline

Definition at line 167 of file QweakSimEPEvent.hh.

References SchwingerDeltaE.

Referenced by QweakSimEPEventMessenger::SetNewValue().

167 { SchwingerDeltaE = energy; } ;
G4double SchwingerDeltaE

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetThetaAngle_Max ( G4double  ang)

Definition at line 4558 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), myUserInfo, and QweakSimUserInformation::SetThetaAngle_Max().

Referenced by CheckLookupTableBounds(), and QweakSimEPEventMessenger::SetNewValue().

void SetThetaAngle_Max(G4double ang)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

void QweakSimEPEvent::SetThetaAngle_Min ( G4double  ang)

Definition at line 4555 of file QweakSimEPEvent.cc.

References CheckLookupTableBounds(), myUserInfo, and QweakSimUserInformation::SetThetaAngle_Min().

Referenced by CheckLookupTableBounds(), and QweakSimEPEventMessenger::SetNewValue().

void SetThetaAngle_Min(G4double ang)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4double QweakSimEPEvent::Sigma_EEPrime ( G4double  eni,
G4double  eprime,
G4double  theta,
G4double &  q2 
)
private

Definition at line 1044 of file QweakSimEPEvent.cc.

References ResMod507().

Referenced by Delta_Resonance().

1045 {
1046 
1047 // Peiqing, Nov 28, 2011
1048 // This is adopted from QWGEANT3 SUBROUTINE SIGMA_EEPRIME (ENI,EPRIME,THETA,WVAL,SIGMA_EEP)
1049 //
1050 // Calculates the differential cross section for inelastic
1051 // ep scattering into unit solid angle and scattered
1052 // electron energy EPRIME. The cross section is
1053 // returned in microbarns/sr/GeV.
1054 
1055  G4double w2,xval1[51],xvalL[51];
1056  G4double sigT,sigL,epsilon;
1057  G4int i;
1058 
1059  G4double t2,epmax,xneu,gamma;
1060 
1061  G4double mp = 0.9382727;
1062  G4double mp2 = mp*mp;
1063  G4double alpha = 1.0/137.036;
1064  G4double mpion = 0.135;
1065 
1066  G4double xval[101]={ 0,
1067  0.12286E+01,0.15200E+01,0.15053E+01,0.16960E+01,0.16466E+01,
1068  0.14055E+01,0.12247E+00,0.22000E+00,0.84591E-01,0.90400E-01,
1069  0.16956E+00,0.29002E+00,0.74046E+01,0.65122E+01,0.58506E+01,
1070  0.16990E+01,0.45729E+01,0.53546E+03,-.97628E+03,0.82539E+05,
1071  0.29494E+01,0.20352E+01,0.12422E+00,0.19560E+01,0.22333E+01,
1072  -.32484E+00,0.24212E+00,0.57737E-01,0.30497E+01,0.63111E+01,
1073  0.37579E+00,0.41100E+01,0.20668E-01,0.45490E+03,0.54493E-01,
1074  0.22565E+01,0.40369E+01,0.43734E+00,0.65625E+00,0.36182E+00,
1075  0.55216E-01,0.41789E+00,0.18104E+00,0.91306E+00,0.18116E+01,
1076  0.24152E+03,0.19329E+01,0.38000E+00,0.28732E+01,-.53116E+00,
1077  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1078  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1079  0.00000E+00,0.00000E+00,0.99495E+01,0.70370E-02,0.16172E+01,
1080  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1081  0.20782E+05,0.79523E+04,0.54933E+00,0.00000E+00,0.68629E+01,
1082  0.10369E+01,0.88112E+00,0.00000E+00,0.35659E+03,0.20281E-02,
1083  0.40336E+02,0.00000E+00,0.45242E+00,0.20000E-01,0.28691E+00,
1084  0.00000E+00,0.25115E+03,0.10663E+01,0.55422E+01,0.30350E+00,
1085  0.10541E+01,0.20527E+01,0.00000E+00,0.00000E+00,0.13055E+02,
1086  0.27997E+01,0.39546E+00,0.46372E+00,0.49972E+01,0.00000E+00};
1087 
1088 // Check that we are not above maximum allowed E assuming threshold W=Mp+Mpi
1089  t2 = theta/2.0;
1090  G4double sin_t2 = sin(t2);
1091  G4double tan_t2 = tan(t2);
1092  G4double sin2_t2 = sin_t2*sin_t2;
1093  G4double tan2_t2 = tan_t2*tan_t2;
1094  G4double sigma_eep = 0.0;
1095  epmax = eni-mpion*(1.0+0.5*mpion/mp);
1096  epmax = epmax/(1.0+(2.0*eni/mp)*sin2_t2);
1097  if (eprime>=epmax)
1098  return 0;
1099 
1100 // Calculate Q**2, W**2, epsilon, gamma
1101  q2 = 4.0*eni*eprime*sin2_t2;
1102  xneu = eni-eprime;
1103  w2 = mp2 + 2.0*mp*xneu - q2;
1104  epsilon = 1.0/(1.0+2.0*(1.0+xneu*xneu/q2)*tan2_t2);
1105  gamma = alpha*eprime*(w2-mp2)/((2.0*pi)*(2.0*pi));
1106  gamma = gamma/(q2*mp*eni*(1.0-epsilon));
1107  // wval = sqrt(w2);
1108 
1109  for (i=1;i<=50;i++)
1110  {
1111  xval1[i] = xval[i];
1112  xvalL[i] = xval[50+i];
1113  if(i<=12)
1114  xvalL[i] = xval1[i];
1115  }
1116 
1117  xvalL[43] = xval1[47];
1118  xvalL[44] = xval1[48];
1119  xvalL[50] = xval1[50];
1120 
1121  sigT = ResMod507(1,w2,q2,&xval1[0]);
1122  sigL = ResMod507(2,w2,q2,&xvalL[0]);
1123  sigma_eep=sigT+epsilon*sigL;
1124  sigma_eep=sigma_eep*gamma;
1125 
1126  return sigma_eep;
1127 }
G4double ResMod507(G4int sf, G4double w2, G4double q2, G4double *xval)

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

G4int QweakSimEPEvent::SuperElasticCheck ( G4double  E_in,
G4double  E_out,
G4double  theta,
G4double &  xsec 
)
private

Definition at line 1505 of file QweakSimEPEvent.cc.

References Elastic_Cross_Section_Proton(), QweakSimUserInformation::GetReactionType(), and myUserInfo.

Referenced by Radiative_Cross_Section_Lookup().

1506 {
1507  if (myUserInfo->GetReactionType()!=7) return 0; // Aborts if not using lookup table
1508 
1509  G4double Weight; // dummy
1510  G4double Q2; // dummy
1511 
1512  G4double E_elastic;
1513 
1514  xsec = Elastic_Cross_Section_Proton(E_in, theta, Weight, Q2, E_elastic);
1515 
1516  if (E_out > E_elastic) {
1517  //G4cout << "###################################################" << G4endl;
1518  //G4cout << "---------------------------------------------------" << G4endl;
1519  //G4cout << "E_in: " << E_in << G4endl;
1520  //G4cout << "E_out: " << E_out << G4endl;
1521  //G4cout << "E_elastic:" << E_elastic << G4endl;
1522  //G4cout << "Theta: " << theta*180.0/3.14159 << G4endl;
1523  //G4cout << "Event is SuperElestic! Generating new event." << G4endl;
1524  return 1;
1525  }
1526  else {
1527  //G4cout << "###################################################" << G4endl;
1528  //G4cout << "E_in: " << E_in << G4endl;
1529  //G4cout << "E_out: " << E_out << G4endl;
1530  //G4cout << "E_elastic:" << E_elastic << G4endl;
1531  //G4cout << "Theta: " << theta*180.0/3.14159 << G4endl;
1532  //G4cout << "Event is NOT SuperElastic" << G4endl;
1533  return 0;
1534  }
1535 }
G4double Elastic_Cross_Section_Proton(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
QweakSimUserInformation * myUserInfo

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

Field Documentation

G4double QweakSimEPEvent::BeamEnergy
private

Definition at line 66 of file QweakSimEPEvent.hh.

G4double QweakSimEPEvent::ElasticPeakDeltaE
private
QweakSimEPEventMessenger* QweakSimEPEvent::EventGen_Messenger
private

Definition at line 87 of file QweakSimEPEvent.hh.

Referenced by QweakSimEPEvent(), and ~QweakSimEPEvent().

QweakSimFieldMap<value_t,value_n>* QweakSimEPEvent::fLookupTable
private
G4int QweakSimEPEvent::Isotropy
private

isotropy used for event generation

Definition at line 57 of file QweakSimEPEvent.hh.

Referenced by GetIsotropy(), QweakSimEPEvent(), and SetIsotropy().

G4int QweakSimEPEvent::kActiveOctantNumber
private

Active octant number in the simulation, 0 will enable all octants.

Definition at line 55 of file QweakSimEPEvent.hh.

Referenced by GetActiveOctantNumber(), QweakSimEPEvent(), and SetActiveOctantNumber().

const G4double QweakSimEPEvent::M_n = 939.5656*MeV
staticprivate

Definition at line 47 of file QweakSimEPEvent.hh.

Referenced by Aluminum_Excited_State().

const G4double QweakSimEPEvent::M_p = 938.2796*MeV
staticprivate
const G4double QweakSimEPEvent::mil
staticprivate

Definition at line 44 of file QweakSimEPEvent.hh.

G4double QweakSimEPEvent::myPositionZ
private

Definition at line 68 of file QweakSimEPEvent.hh.

Referenced by Pion_PhotoProduction(), and Pion_PhotoProductionCarbon().

G4int QweakSimEPEvent::ReactionRegion
private
G4int QweakSimEPEvent::ReactionType
private

reaction type used for event generation

Definition at line 53 of file QweakSimEPEvent.hh.

Referenced by GetReactionType(), QweakSimEPEvent(), SetReactionRegion(), and SetReactionType().

G4double QweakSimEPEvent::SchwingerDeltaE
private
const G4double QweakSimEPEvent::TargetLength
staticprivate

Definition at line 41 of file QweakSimEPEvent.hh.

const G4double QweakSimEPEvent::TargetWindowThickness
staticprivate

Definition at line 45 of file QweakSimEPEvent.hh.

const G4double QweakSimEPEvent::Theta_Min = 1.745329E-4
staticprivate
G4int QweakSimEPEvent::TypeSetting
private

Definition at line 52 of file QweakSimEPEvent.hh.

Referenced by QweakSimEPEvent(), and SetReactionType().

const G4int QweakSimEPEvent::value_d = 4
staticprivate

Definition at line 75 of file QweakSimEPEvent.hh.

Referenced by CreateLookupTable(), and Radiative_Cross_Section_Lookup().

const G4int QweakSimEPEvent::value_n = 15
staticprivate

Definition at line 74 of file QweakSimEPEvent.hh.

Referenced by Radiative_Cross_Section_Lookup().


The documentation for this class was generated from the following files: