class QweakMonteCarlo
{
public:
QweakMonteCarlo();
double beamEnergy;
double thetaMin;
double thetaMax;
double phiMin;
double phiMax;
double SummedCrossSection[10000];
int SummedCrossSection_LookupTable_Flag;
double thetaE, thetaP, phiE, phiP;
double E0;
double pi;
double ElectronMomentumX;
double ElectronMomentumY;
double ElectronMomentumZ;
double ElectronEnergy_AfterVertex;
double ProtonMomentumX;
double ProtonMomentumY;
double ProtonMomentumZ;
double ProtonEnergy_AfterVertex;
int nStepTheta;
double dSigma_dOmega_CrossSection_Rutherford (double e_i, double theta);
double dSigma_dOmega_CrossSection_Mott (double e_i, double theta);
double dSigma_dOmega_CrossSection_epElastic (double e_i, double theta);
double Asymmetry (double e_i, double theta);
void InitializeEventGenerator(double myBeamEnergy, double myThetaMin, double myThetaMax, double myPhiMin, double myPhiMax);
void EventGenerator();
//private:
// Sine squared of the weak mixing angle
double Sin2W;
// The fine structure constant
double Alpha;
// the electron mass in GeV
double M_e;
// the proton mass in GeV
double M_p;
// The Fermi Constant in inverse GEV squared
double G_F;
// the magnetic moment of the proton
double Mu_p;
// the magnetic moment of the neutron
double Mu_n;
// the vector part of the proton's weak radiative correction
double RV_p;
// the vector part of the neutron's weak radiative correction
double RV_n;
// the stange vector part of the proton's weak radiative correction
double RV_s;
// the isocalar part of the wealk axial radiative correction
double RA_T0;
// The isovector part of the wealk axial radiative correction
double RA_T1;
// the singlet contribution to the weal axial radiative correction
double RA_0;
// The radiative corrected value in the standard model of the weak charge of the proton
double QWeak;
// Dipole model parameter for the vector form factor of the neuleon
double Lambda_vector;
// Dipole model parameter for the axialvector form factor of the neuleon
double Lambda_axial;
// Dipole model parameter for the strange vector form factor of the neuleon
double Lambda_strange;
// Galster parameter for the electric form factor of the neutron
double Gaster;
// galster parameter for the strange electric form factor
double Gaster_strange;
// the charge radius parameter slope of neutron's charge form factor as tau goes to zero
double Rho_n;
// the charge radius parameter for the strange quark currents
double Rho_s;
// the contibution of strange quarks to the magnetic moment of the the proton
double Mu_s;
// an axial form factor parameter
double g_A;
// an axial form factor parameter
double F_A;
// mikes code converted from C#
double Beta(double energy, double mass);
double Gamma(double energy, double mass);
double S2(double theta);
double E_f(double e_i, double theta);
double E_i(double e_f, double theta);
double Proton_KineticEnergy_AfterVertex(double e_i, double theta);
double Proton_ThetaAngle_AfterVertex(double e_i, double theta);
double Q2(double e_i, double e_f, double theta);
double Q2(double e_i, double theta);
double Tau_p(double q2, double m_A);
double Tau_p(double q2);
double Epsilon(double e_i, double theta);
double GDipole(double q2, double lambda);
double GE_p(double q2);
double GE_n(double q2);
double GM_p(double q2);
double GM_n(double q2);
double FF2(double e_i, double theta);
double GE_s(double q2);
double GM_s(double q2);
double GE_Z(double q2);
double GM_Z(double q2);
double GAxial(double q2);
double GA_3(double q2);
double GA_8(double q2);
double GA_s(double q2);
double GA(double q2);
double Schwinger(double q2);
double T_virtual(double q2);
double T_virtual_ave(double angle);
double Nu(double random,double t);
double GammaRadProbability(double tv, double tr, double nu);
};
QweakMonteCarlo::QweakMonteCarlo()
{
E0 = beamEnergy = 1.165; // GeV
thetaMin = 0.0; // deg
thetaMax = 0.0;
phiMin = 0.0;
phiMax = 0.0;
SummedCrossSection_LookupTable_Flag = 0; // not yet generated
ElectronMomentumX = 0.0;
ElectronMomentumY = 0.0;
ElectronMomentumZ = 0.0;
ElectronEnergy_AfterVertex = 0.0;
ProtonMomentumX = 0.0;
ProtonMomentumY = 0.0;
ProtonMomentumZ = 0.0;
ProtonEnergy_AfterVertex = 0.0;
pi = 4.0*atan(1.0);
// Sine squared of the weak mixing angle
Sin2W = 0.23124;
// The fine structure constant
Alpha = 1/137.03599;
// the electron mass in GeV
M_e = 0.000511; // to be fixed
// the proton mass in GeV
M_p = 0.938;
// The Fermi Constant in inverse GEV squared
G_F = 1.16639E-5;
// the magnetic moment of the proton
Mu_p = 2.79;
// the magnetic moment of the neutron
Mu_n = -1.91;
// the vector part of the proton's weak radiative correction
RV_p = -0.048;
// the vector part of the neutron's weak radiative correction
RV_n = -0.021;
// the stange vector part of the proton's weak radiative correction
RV_s = -0.0113;
// the isocalar part of the wealk axial radiative correction
RA_T0 = -0.0358;
// The isovector part of the wealk axial radiative correction
RA_T1 = -0.42;
// the singlet contribution to the weal axial radiative correction
RA_0 = 0.0;
// The radiative corrected value in the standard model of the weak charge of the proton
QWeak = (1 - 4*Sin2W)*(1 - RV_p)/4.0; // weak charge of proton
// Dipole model parameter for the vector form factor of the neuleon
Lambda_vector = 4.960;
// Dipole model parameter for the axialvector form factor of the neuleon
Lambda_axial = 3.32;
// Dipole model parameter for the strange vector form factor of the neuleon
Lambda_strange = Lambda_vector;
// Galster parameter for the electric form factor of the neutron
Gaster = 5.6;
// galster parameter for the strange electric form factor
Gaster_strange = Gaster;
// the charge radius parameter slope of neutron's charge form factor as tau goes to zero
Rho_n = -Mu_n;
// the charge radius parameter for the strange quark currents
Rho_s = 0.0;
// the contibution of strange quarks to the magnetic moment of the the proton
Mu_s = 0.0;
// an axial form factor parameter
g_A = 1.262;
// an axial form factor parameter
F_A = 0.64;
return;
}
double QweakMonteCarlo::dSigma_dOmega_CrossSection_Rutherford(double e_i, double theta)
{
// Rutherford differential cross section taken from Povh: "Particles and Nuclei", formula (5.36), page 59
// (setting hbar = c = 1)
//
// / dSigma \ Z^2 * Alpha^2
// | ------ | = ---------------------
// \ dOmega / Rutherford 4*E^2*sin(theta/2)^4)
//
// (differential cross section: in lab frame for electron detection)
double Z_Target = 1; // up now we only scatter on hydrogen (proton)
double dSigma_Rutherford = (Z_Target*Z_Target)*(Alpha*Alpha) / ( 4*e_i*e_i*TMath::Power(TMath::Sin(theta/2),4) );
return dSigma_Rutherford;
}
double QweakMonteCarlo::dSigma_dOmega_CrossSection_Mott(double e_i, double theta)
{
// This Mott cross section *does* include the proton recoil factor E_f/E_i
//
// differentail Mott c.s. taken from Thomas, Weise: "The Structure of the Nucleon", formula (2.7), page 9
//
// / dSigma \ / dSigma \ e_f
// | -------- | = | -------- | * ---- * cos(theta/2)^2
// \ dOmega_e / Mott \ dOmega_e / Rutherford e_i
//
// (differential cross section: in lab frame for electron detection)
double e_f = E_f(e_i,theta);
double dSigma_Mott = dSigma_dOmega_CrossSection_Rutherford(e_i, theta) * (e_i/e_f) * TMath::Power(TMath::Cos(theta/2),2);
return dSigma_Mott;
}
double QweakMonteCarlo::dSigma_dOmega_CrossSection_epElastic(double e_i, double theta)
{
// differential cross section of the elastic electron-proton scattering in one-photon approximation
// (non-radiative corrected)
//
// expression taken from arXiv:nucl-ex/0410010 v1 8 Oct 2004
// J. Arrington et al. : "Precision Rosenbluth measurement of the proton elastic form factors"
//
// in the paper Sigma_Mott stands for dSigma/dOmega_Mott
//
//
// / dSigma \ / dSigma \ [tau*GM_p^2 + epsilon*GE_p^2]
// | -------- | = | -------- | * -----------------------------
// \ dOmega_e / ep_elas \ dOmega_e / Mott epsilon*(1+tau)
//
//
// (differential cross section: in lab frame for electron detection)
double q2 = Q2(e_i,theta);
double dSigma_ep = dSigma_dOmega_CrossSection_Mott(e_i,theta)/(Epsilon(e_i,theta)*(1 + Tau_p(q2))) *FF2(e_i,theta);
return dSigma_ep;
}
void QweakMonteCarlo::InitializeEventGenerator(double myBeamEnergy,
double myThetaMin, double myThetaMax,
double myPhiMin, double myPhiMax)
{
beamEnergy = myBeamEnergy;
thetaMin = myThetaMin*pi/180.;
thetaMax = myThetaMax*pi/180.;
phiMin = myPhiMin*pi/180.;
phiMax = myPhiMax*pi/180.;
EventGenerator();
return;
}
void QweakMonteCarlo::EventGenerator(){
int nStepTheta;
double stepTheta;
double theta;
int j, jLow, jHigh, iflag;
double r, fstep, Eprime;
double pi = 4.0*atan(1.0);
E0 = beamEnergy;
nStepTheta = 10000;
stepTheta = (thetaMax-thetaMin)/(nStepTheta-1);
if (thetaMin == thetaMax || phiMin == phiMax) return;
// fill lookuptable only once, check status flag
if (SummedCrossSection_LookupTable_Flag == 0 ){
//==================================================
// summing (integrating) the diff. c.s. over theta
//==================================================
//
// ==> This needs only be performed the first time the Event Generator is called
//
// with the diff cross section (not phi dependend)
//
// dSigma/dOmega = f(theta)
//
// the integration over theta will look like (see next section too)
//
// dOmega => delta_Omega = cos(theta)*delta_theta*delta_phi
//
//
// theta_max
// ------
// \
// Sigma = f(theta)*cos(theta)
// /
// ------
// theta = theta_min
//
// (leaving out the (scaling) factor delta_theta*delta_phi)
//==================================================================================
//
// It is important how the electron scattering angle is defined in respect to
// the decomposition into the coordinate axes and for the definition of dOmega !!!
//
// In our case the electron scattering angle is defined by theta and is therefore
// defined in respect to the beam axis along the z axis.
// Therefore the c.s. weighting factor is here cos(theta) and *not* sin(theta)
//
// | / scattered electron
// |delta/
// | / ==> dOmega = sin(delta)*ddelta*dphi ,
// | /
// | / ==> dOmega = cos(theta)*dtheta*dphi
// | /
// |/ theta
// |---------------------------> beam axis / z axis
//
//===================================================================================
for (int i = 0; i< nStepTheta ; i++){
theta = thetaMin + i*stepTheta;
if (i == 0) SummedCrossSection[i] = dSigma_dOmega_CrossSection_epElastic(E0,theta) * TMath::Cos(theta);
if (i != 0) SummedCrossSection[i] = SummedCrossSection[i-1]
+ (dSigma_dOmega_CrossSection_epElastic(E0,theta) * TMath::Cos(theta));
}
SummedCrossSection_LookupTable_Flag = 1; // set flag for lookup table has been filled
}
// This is the uniform distribution.
// The GenSigmaSum[nStepTheta-1] is to normalize it
r = gRandom->Rndm()*SummedCrossSection[nStepTheta-1];
// This is searching table to find a F(x) which matches the
// uniform random number r
jLow = 0;
jHigh = nStepTheta;
iflag = j = 0;
while((jHigh-jLow)>1 && iflag == 0)
{
j=(jHigh+jLow)/2;
if(r<SummedCrossSection[j-1] ) jHigh = j;
if(r>SummedCrossSection[j-1] ) jLow = j;
if(SummedCrossSection[j-1] <= r && r <= SummedCrossSection[j]) iflag = 1;
}
// In general r will fall between two table values.
// Here we make a linear interpolation (assuming small step size in theta)
fstep = (r-SummedCrossSection[j-1])/(SummedCrossSection[j] - SummedCrossSection[j-1]);
theta = thetaMin + (j+fstep)*stepTheta;
// Once we have theta of the electron, we can extract all the kinematics
// electron theta scattering angle
thetaE = theta;
// Generate the phi angle for the scattered electron with an uniform(flat) angle distribution
// within the defined detectors acceptance
phiE = phiMin + (phiMax-phiMin)*gRandom->Rndm();
// corresponding decomposition
ElectronMomentumX = ElectronEnergy_AfterVertex*sin(thetaE)*cos(phiE); //GeV
ElectronMomentumY = ElectronEnergy_AfterVertex*sin(thetaE)*sin(phiE); //GeV
ElectronMomentumZ = ElectronEnergy_AfterVertex*cos(thetaE) ; //GeV
//double q2, Eq;
// energy of scattered electron
// Eprime = (M_p*beamEnergy)/(2.0*beamEnergy*pow(sin(theta/2.0),2) + M_p);
ElectronEnergy_AfterVertex = E_f(E0,theta);
// Transfered Momentum (only valid for elastic scattering)
// see Povh, expression (7.10), page 87
//q2 = 2*M_p*(beamEnergy-Eprime);
// Transferred energy to proton thru virtual photon (in lab frame)
// see Povh, expression (7.3), page 84
//Eq = beamEnergy-Eprime;
// Proton energy (GeV)
ProtonEnergy_AfterVertex = Proton_KineticEnergy_AfterVertex(E0, theta);
// proton scattering angle (rad)
thetaP = Proton_ThetaAngle_AfterVertex(E0, theta);
// the proton will be scattered to the opposite electron phi angle
// incident electron, scattered electron and scattered proton are all
// in the scattering plane
phiP = phiE + pi;
}
// the square of the half scattering angle, a commonly used kinematic quantity in electron scattering
// input parameter: "theta" , scattering angle
// output : sin^2(theta/2)
double QweakMonteCarlo::S2(double theta)
{
double s = sin(theta/2);
return s*s;
}
// calculates the reativistic gamma of the electron
// <param name="energy">electron energy</param>
// <param name="mass">electron mass</param>
// <returns>a double</returns>
double QweakMonteCarlo::Gamma(double energy, double mass)
{
return energy/mass;
}
/// <summary>
/// calculates the dimensionless velocity of the electron
/// </summary>
/// <param name="energy">electron energy</param>
/// <param name="mass">electron mass</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Beta(double energy, double mass)
{
double gamma = Gamma(energy,mass);
double beta = sqrt(1 - 1/gamma*gamma);
return beta;
}
/// <summary>
/// Invariant four-momentum-transfer-Squared for *inelastic* electron scattering =-(e-e')^2
/// in a time like basis
/// </summary>
/// <param name="e_i">energy of incident electron </param>
/// <param name="e_f">final electron energy</param>
/// <param name="theta">scattering angle</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Q2(double e_i, double e_f, double theta)
{
double m_e = M_e;
double beta_i = Beta(e_i,m_e);
double beta_f = Beta(e_f,m_e);
double Q2 = -2.0*(m_e*m_e - e_i*e_f*(1 - beta_i*beta_f*cos(theta)));
return Q2;
}
/// <summary>
/// Invariant four-momentum-transfer-Squared for *elastic* electron scattering =-(e-e')^2
/// in a time like basis
/// </summary>
/// <param name="e_i">energy of incident electron </param>
/// <param name="theta">scattering angle</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Q2(double e_i, double theta)
{
double e_f = E_f(e_i,theta);
double q2 = Q2(e_i,e_f,theta);
return q2;
}
double QweakMonteCarlo::Proton_KineticEnergy_AfterVertex(double e_i, double theta)
{
// expression taken from PVA4 Geant : junk
// double E_f_Proton = e_i / ( 1 + M_p / (e_i*( 1-cos(theta))) );
double q2 = Q2(e_i,theta);
double E_f_Proton = q2/(2*M_p*M_p);
return E_f_Proton;
}
double QweakMonteCarlo::Proton_ThetaAngle_AfterVertex(double e_i, double theta)
{
// based on Blast code (Tim Smith)
double e_f = E_f(e_i,theta);
double Eq = e_i - e_f;
// Transfered Momentum (only valid for elastic scattering)
// see Povh, expression (7.10), page 87
double q2 = 2*M_p*Eq;
double pq = sqrt(Eq*Eq + q2);
// expression taken from PVA4 Geant : junk
//double ThetaP = asin(sqrt( e_f*e_f/q2/(1+tau)*sin(theta)*sin(theta) ));
double ThetaP = asin( e_f/pq * sin(theta));
return ThetaP;
}
/// <summary>
/// Final electron energy in elastic electron scattering;
/// </summary>
/// <param name="e_i"></param>
/// <param name="theta"></param>
/// <returns></returns>
/// <remarks> evaluated in the ERL</remarks>
double QweakMonteCarlo::E_f(double e_i, double theta)
{
// standard expression (see e.g. Povh, page 78)
double E_f = e_i/(1 + e_i/M_p*(1-cos(theta)));
return E_f;
}
/// <summary>
/// initial electron energy in elastic electron scattering
/// </summary>
/// <param name="e_f">scattered electron energy</param>
/// <param name="theta">scattering angle</param>
/// <returns></returns>
double QweakMonteCarlo::E_i(double e_f, double theta)
{
double E_i = e_f/(1 - (2*e_f/M_p)*S2(theta));
return E_i;
}
/// <summary>
/// a kinematic quantity representing the dimensionless Q2 for elastic scattering from a mass A
/// </summary>
/// <param name="q2">four momentum squared</param>
/// <param name="m_A">the mass of target A</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Tau_p(double q2, double m_A)
{
double Tau_p = q2/4.0/m_A/m_A;
return Tau_p;
}
/// <summary>
/// dimensionless Q2 for the proton
/// </summary>
/// <param name="q2">four momentum squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Tau_p(double q2)
{
double m_p = M_p;
double Tau_p = q2/4.0/m_p/m_p;
return Tau_p;
}
/// <summary>
/// Longitudinal polarization of the virtual photon
/// </summary>
/// <param name="e_i">electron energy</param>
/// <param name="theta">scattering angle</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Epsilon(double e_i, double theta)
{
double m_p = M_p;
double q2 = Q2(e_i,theta);
double tau = Tau_p(q2);
double w = tau/m_p/2;
double e_f = e_i - w;
double a = e_i*e_f/m_p/m_p; //the invariant (s*P_T)(p*P_T)/(P_T*P_T)) in the fixed target lab frame
double Epsilon = 1.0/(1.0 + 2.0*(1 + tau)*tau/(a - tau));
return Epsilon;
}
/// <summary>
/// Dipole parameterization
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <param name="lambda">dipole scale parameter</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GDipole(double q2, double lambda)
{
double GDipole = 1.0/TMath::Power((1.0 + lambda*Tau_p(q2)),4.0);
return GDipole;
}
/// <summary>
/// Electric form factor of the proton
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GE_p(double q2)
{
double GE_p = GDipole(q2, Lambda_vector);
return GE_p;
}
/// <summary>
/// Electric form factor of the neutron
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GE_n(double q2)
{
double GE_n = Rho_n*GDipole(q2, Lambda_vector)/(1 + Gaster*Tau_p(q2));
return GE_n;
}
/// <summary>
/// Magnetic form factor of the proton
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GM_p(double q2)
{
double GM_p = Mu_p*GDipole(q2,Lambda_vector);
return GM_p;
}
/// <summary>
/// Magnetic form factor of the neutron
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GM_n(double q2)
{
double GM_n = Mu_n*GDipole(q2, Lambda_vector);
return GM_n;
}
/// <summary>
/// Total form factor (Electric + Magnetic)
/// </summary>
/// <param name="e_i">electron energy</param>
/// <param name="theta">scattering angle</param>
/// <returns>a double</returns>
double QweakMonteCarlo::FF2(double e_i, double theta)
{
double q2 = Q2(e_i,theta);
double FF2 = Epsilon(e_i,theta)*GE_p(q2)*GE_p(q2) + Tau_p(q2)*GM_p(q2)*GM_p(q2);
return FF2;
}
/// <summary>
/// Strange electric form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GE_s(double q2)
{
double GE_s = Rho_s*Tau_p(q2)*GDipole(q2, Lambda_strange)/(1 + Gaster_strange*Tau_p(q2));
return GE_s;
}
/// <summary>
/// Strange mafnetic form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GM_s(double q2)
{
double GM_s = Mu_s*GDipole(q2, Lambda_strange);
return GM_s;
}
/// <summary>
/// weak electric form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GE_Z(double q2)
{
double GE_Z = QWeak*GE_p(q2)-(1.0 + RV_n)*GE_n(q2)/4.0 - (1.0 + RV_s)*GE_s(q2)/4.0;
return GE_Z;
}
/// <summary>
/// weak magnetic form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GM_Z(double q2)
{
double GM_Z = QWeak*GM_p(q2) - (1.0 + RV_n)*GM_n(q2)/4.0 - (1.0 + RV_s)*GM_s(q2)/4.0;
return GM_Z;
}
/// <summary>
/// generic dipole type parameterization of the shape of the axial form factor
/// </summary>
/// <param name="q2">momentum transfered squared</param>
/// <returns>a double</returns>
// Axial form factors
double QweakMonteCarlo::GAxial(double q2)
{
double GAxial = GDipole(q2,Lambda_axial);
return GAxial;
}
/// <summary>
/// the triplet part of the weak axial form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GA_3(double q2)
{
double GA_3 = g_A*GAxial(q2)/2.0;
return GA_3;
}
/// <summary>
/// the octet part of the axial form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GA_8(double q2)
{
double GA_8 = g_A*GAxial(q2)*(3*F_A - 1.0)/(F_A + 1)/2.0/TMath::Sqrt(3.0);
return GA_8;
}
/// <summary>
/// strange quark contribution to the axial form factor
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GA_s(double q2)
{
GA_s = 0.0;
return GA_s;
}
///<summary>
///The total axial form factor contribution to the asymmetry
///</summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::GA(double q2)
{
double GA = -2.0*GA_3(q2)*(1.0 + RA_T1) + Math.Sqrt(3.0)*RA_T0*GA_8(q2) + (1 + RA_0)*GA_s(q2);
return GA;
}
/// <summary>
/// the parity violating assymetry
/// </summary>
/// <param name="e_i">electron energy</param>
/// <param name="theta">scatterning angle</param>
/// <returns></returns>
double QweakMonteCarlo::Asymmetry(double e_i,double theta)
{
double q2 = Q2(e_i,theta);
double e_f = E_f(e_i,theta);
double ep = Epsilon(e_i,theta);
double tau = Tau_p(q2);
double k = G_F*q2/Math.Sqrt(2.0)/Math.PI/Alpha/FF2(e_i,theta);
double A_V = k*(ep*GE_p(q2)*GE_Z(q2) + tau*GM_p(q2)*GM_Z(q2));
double A_A = k*(Math.Sqrt(ep)*(1.0/4.0 - Sin2W)*Math.Sqrt(tau*(1 + tau))*GM_p(q2)*GA(q2));
double Asymmetry = A_V + A_A;
return Asymmetry;
}
/// <summary>
/// the finite renormalization part of the Swinger radiative correction
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::Schwinger(double q2)
{
double m_e2 = M_e*M_e;
double Schwinger = -(2.0*Alpha/Math.PI)*((28.0/18.0) - (13.0/12.0)*Math.Log(q2/m_e2));
return Schwinger;
}
/// <summary>
/// the virtual photon's effective radiator thickness
/// </summary>
/// <param name="q2">momentum transfer squared</param>
/// <returns>a double</returns>
double QweakMonteCarlo::T_virtual(double q2)
{
double m_e2 = M_e*M_e;
double T_virtual = (Alpha/Math.PI)*(Math.Log(q2/m_e2) - 1.0);
return T_virtual;
}
/// <summary>
/// an estimate of t_virtual made by using the point scattering result at the target center.
/// </summary>
/// <parameter name="angle">the nominal central scattering angle</parameter>
/// <returns>virtual radiator thickness in radiation lengths</returns>
double QweakMonteCarlo::T_virtual_ave(double angle)
{
// Target t = Target.Estimate(angle);
double m_e2 = M_e*M_e;
//double e0 =Setup.Parameters.BeamLine.IncidentEnergy*Units.GEV;
double e0 = E0;
// double eVertex = e0 - t.EnergyLossBefore;
double eVertex = e0;
double q2 = Q2(eVertex, angle);
double T_virtual_ave = (Alpha/TMath::Pi())*(TMath::Log(q2/m_e2) - 1.0);
return T_virtual_ave;
}
/// <summary>
/// the factional energy loss due to radiation using a randomn probability weighting function
/// </summary>
/// <param name="random">a random number between zero and 1</param>
/// <param name="t"> the local value of the effective radiator</param>
/// <returns>the fradctional energy loss along the direction of radiation</returns>
double QweakMonteCarlo::Nu(double random,double t)
{
double Nu = TMath::Exp(TMath::Log(random)/t);
return Nu;
}
/// <summary>
/// a simple estimate of the fraction of events expected to be lost due to target and internal bremstralung
/// </summary>
/// <param name="tv">virtual radiator</param>
/// <param name="tr">thick target radiator</param>
/// <param name="nu">fraction energy loss</param>
/// <returns></returns>
double QweakMonteCarlo::GammaRadProbability(double tv, double tr, double nu)
{
double b0 = 4.0/3.0;
double a = ((1 - nu) + (b0*tr/tv)*(1 - nu - 3.0*nu*nu/4.0));
double b = (tv*(1 - nu)) + b0*tr*((5.0/8.0) - nu - (3.0/8.0)*nu*nu);
double GammaRadProbability = a*Math.Exp(b);
return GammaRadProbability;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
QweakMonteCarlo()
{
int Nevents = 100000;
... 162 more lines ...
|