EICMCAnalysis/analysis_eicsmear/calculate_g1.C

0001 /** Parametrization for F2(x,Q2)
0002  */
0003 double
0004 get_parametrization_F2(double x, double Q2)
0005 {
0006   /* based on reference from Ernst's thesis */
0007   double a1 = -0.24997;
0008   double a2 = 2.3963;
0009   double a3 = 0.23643;
0010   double a4 = 0.08498;
0011   double a5 = 3.8608;
0012   double a6 = -7.4143;
0013   double a7 = 3.4342;
0014   double b1 = 0.11411;
0015   double b2 = -2.2356;
0016   double b3 = 0.03115;
0017   double b4 = 0.02135;
0018   double c1 = -1.4517;
0019   double c2 = 8.4745;
0020   double c3 = -34.379;
0021   double c4 = 45.888;
0022
0023   double A = pow(x, a1) * pow((1 - x), a2) * (a3 * pow((1 - x), 0) +
0024                                               a4 * pow((1 - x), 1) +
0025                                               a5 * pow((1 - x), 2) +
0026                                               a6 * pow((1 - x), 3) +
0027                                               a7 * pow((1 - x), 4));
0028   double B = b1 + b2 * x + (b3 / (x + b4));
0029   double C = ( c1 * pow(x, 1) +
0030                c2 * pow(x, 2) +
0031                c3 * pow(x, 3) +
0032                c4 * pow(x, 4) );
0033   double Q02 = 20;
0034   double lambda2 = 1;
0035
0036   double F2_fit = A * pow( ( log(Q2 / lambda2) / log(Q02 / lambda2) ), B) * (1 + (C / Q2));
0037
0038   return F2_fit;
0039 }
0040
0041 /** Parametrization for R(x,Q2)
0042  */
0043 double
0044 get_parametrization_R(double x, double Q2)
0045 {
0046   double R_tmp = 1;
0047
0048   if ( x < 0.12 )
0049     {
0050       R_tmp = 1.509 * pow(x, -0.0458) * pow((1-x), -0.0084) / log((Q2 / 0.04)) - 0.203;
0051     }
0052   else
0053     {
0054       /* below parametrization returns NAN- debug and use instead of this line for x > 0.12 */
0055       R_tmp = 1.509 * pow(x, -0.0458) * pow((1-x), -0.0084) / log((Q2 / 0.04)) - 0.203;
0056     }
0057   //    {
0058   //      // R from fit to SLAC data
0059   //      double b1 = 0.635;
0060   //      double b2 = 0.5747;
0061   //      double b3 = -0.3534;
0062   //      double lambda = 0.2;
0063   //      double Theta = 1 + 12 * (Q2 / (Q2 + 1)) * (pow(0.125, 2) / (pow(0.125, 2) + pow(x, 2)));
0064   //      R_tmp = b1 / log( Q2 / pow(lambda, 2)) * Theta + b2 / Q2 + b3 / (pow(Q2, 2) + pow(0.3, 2)));
0065   //    }
0066
0067   return R_tmp;
0068 }
0069
0070 /** Calculate depolarization factor D(x,Q2,y,m_lepton)
0071  */
0072 double
0073 get_depolarization_factor(double x, double Q2, double y, double mlepton=0.000511)
0074 {
0075   /* Parameters needed */
0076   double M = 0.938; // nucleon mass in GeV
0077   double ml = mlepton; // lepton mass in GeV
0078
0079   double gamma2 = pow(( 2 * M * x ), 2) / Q2; // ?
0080
0081   double R = get_parametrization_R(x, Q2); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
0082
0083   //double D_num = (
0084   //                y * ( 2 - y ) * ( 1 + gamma2 * y / 2. )
0085   //                );
0086   double D_num = (
0087                   y * ( 2 - y ) * ( 1 + gamma2 * y / 2. ) - 2 * pow(y, 3) * pow(ml, 2) / Q2
0088                   ); // Ernst's thesis has the extra y^3 term
0089   double D_denom = (
0090                     pow(y, 2) * ( 1 + gamma2 ) * (1 - 2 * pow(ml, 2) / Q2 )
0091                     + 2 * (1 - y - gamma2 * pow(y, 2) / 4 ) * (1 + R )
0092                     );
0093   double D = D_num / D_denom; // Depolarization factor
0094
0095   return D;
0096 }
0097
0098 /** Calculate beam polarization
0099  */
0100 double
0101 get_beam_polarization()
0102 {
0103   /* assumed beam polarizarion factors */
0104   double pol_electron = 0.7;
0105   double pol_proton = 0.7;
0106   double pol_prod = pol_electron * pol_proton;
0107
0108   return pol_prod;
0109 }
0110
0111 /** Calculate value of g1
0112  */
0113 double
0114 get_g1_value( double x_val, double Q2_val, double y_val, double mlepton=0.000511 )
0115 {
0116   /* F2 and R parametrizations (evaluated at kinmatics of given data point) */
0117   double F2 = get_parametrization_F2(x_val, Q2_val);
0118   double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
0119
0120   /* Depolarization factor */
0121   double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
0122
0123   /* Asymmetry measured in parallel spin configuration */
0124   // double A_parralel = D * A1; // assuming A2 = 0  // Ernst's thesis Eq. 3.41
0125
0126   /* A_raw = (N++ - N+-) / (N++ + N+-) */
0127   //double A_raw = ?
0128   //double A_parallel = 1. / pol_prod * A_raw;
0129
0130   double A1_Ernst = 0.033;
0131   double A_parallel = D * A1_Ernst;
0132
0133   double g1_val = 1. / (2 * x_val * (1 + R)) * F2 * A_parallel / D;
0134
0135   //  return g1_val;
0136   return 0;
0137 }
0138
0139 /** Caluclate g1 statistical uncertainty
0140  */
0141 double
0142 get_g1_sigma( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
0143 {
0144   /* F2 and R parametrizations (evaluated at kinmatics of given data point) */
0145   double F2 = get_parametrization_F2(x_val, Q2_val);
0146   double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
0147
0148   /* Depolarization factor */
0149   double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
0150
0151   /* Beam polarization */
0152   double pol_prod = get_beam_polarization();
0153
0154   double A_parallel_sig = (sqrt(N_val) / N_val) * (1. / pol_prod);
0155   double g1_sig =  1. / (2 * x_val * (1 + R)) * F2 * A_parallel_sig / D;
0156
0157   return g1_sig;
0158 }
0159
0160 /** Caluclate A1 statistical uncertainty
0161  */
0162 double
0163 get_A1_sigma( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
0164 {
0165   /* Calculate F1 */
0166   double M = 0.938; // nucleon mass in GeV
0167
0168   double gamma2 = pow(( 2 * M * x_val ), 2) / Q2_val; // ?
0169
0170   double F2 = get_parametrization_F2(x_val, Q2_val);
0171   double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
0172
0173   double F1 = ( ( 1 + gamma2 ) / ( 2 * x_val * ( 1 + R ) ) ) * F2;
0174
0175   /* Calculate A1 sigma (note: A1 = g1 / F1) */
0176   double g1_sig = get_g1_sigma( x_val, Q2_val, y_val, N_val, mlepton );
0177   double A1_sig = g1_sig / F1;
0178
0179   return A1_sig;
0180 }
0181
0182 /** Printing values on screen
0183  */
0184 void
0185 calculate_g1( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
0186 {
0187   double g1_val = get_g1_value(x_val, Q2_val, y_val, mlepton);
0188   double g1_sig = get_g1_sigma(x_val, Q2_val, y_val, N_val, mlepton);
0189   double F2 = get_parametrization_F2(x_val, Q2_val);
0190   double R = get_parametrization_R(x_val, Q2_val);
0191   double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
0192
0193   /* print parameters */
0194   cout << "-----------------------------------------------------------------------" << endl;
0195   cout << "x = " << x_val << " , Q2 = " << Q2_val << " , y = " << y_val << " , N = " << N_val << endl;
0196   cout << "F2 = " << F2 << ", R = " << R << ", D = " << D << endl;
0197   cout << "==> g1 = " << g1_val << " +/- " << g1_sig << endl;
0198   cout << "-----------------------------------------------------------------------" << endl;
0199
0200   return;
0201 }