Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
calculate_g1.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file calculate_g1.C
1 
3 double
4 get_parametrization_F2(double x, double Q2)
5 {
6  /* based on reference from Ernst's thesis */
7  double a1 = -0.24997;
8  double a2 = 2.3963;
9  double a3 = 0.23643;
10  double a4 = 0.08498;
11  double a5 = 3.8608;
12  double a6 = -7.4143;
13  double a7 = 3.4342;
14  double b1 = 0.11411;
15  double b2 = -2.2356;
16  double b3 = 0.03115;
17  double b4 = 0.02135;
18  double c1 = -1.4517;
19  double c2 = 8.4745;
20  double c3 = -34.379;
21  double c4 = 45.888;
22 
23  double A = pow(x, a1) * pow((1 - x), a2) * (a3 * pow((1 - x), 0) +
24  a4 * pow((1 - x), 1) +
25  a5 * pow((1 - x), 2) +
26  a6 * pow((1 - x), 3) +
27  a7 * pow((1 - x), 4));
28  double B = b1 + b2 * x + (b3 / (x + b4));
29  double C = ( c1 * pow(x, 1) +
30  c2 * pow(x, 2) +
31  c3 * pow(x, 3) +
32  c4 * pow(x, 4) );
33  double Q02 = 20;
34  double lambda2 = 1;
35 
36  double F2_fit = A * pow( ( log(Q2 / lambda2) / log(Q02 / lambda2) ), B) * (1 + (C / Q2));
37 
38  return F2_fit;
39 }
40 
43 double
44 get_parametrization_R(double x, double Q2)
45 {
46  double R_tmp = 1;
47 
48  if ( x < 0.12 )
49  {
50  R_tmp = 1.509 * pow(x, -0.0458) * pow((1-x), -0.0084) / log((Q2 / 0.04)) - 0.203;
51  }
52  else
53  {
54  /* below parametrization returns NAN- debug and use instead of this line for x > 0.12 */
55  R_tmp = 1.509 * pow(x, -0.0458) * pow((1-x), -0.0084) / log((Q2 / 0.04)) - 0.203;
56  }
57  // {
58  // // R from fit to SLAC data
59  // double b1 = 0.635;
60  // double b2 = 0.5747;
61  // double b3 = -0.3534;
62  // double lambda = 0.2;
63  // double Theta = 1 + 12 * (Q2 / (Q2 + 1)) * (pow(0.125, 2) / (pow(0.125, 2) + pow(x, 2)));
64  // R_tmp = b1 / log( Q2 / pow(lambda, 2)) * Theta + b2 / Q2 + b3 / (pow(Q2, 2) + pow(0.3, 2)));
65  // }
66 
67  return R_tmp;
68 }
69 
72 double
73 get_depolarization_factor(double x, double Q2, double y, double mlepton=0.000511)
74 {
75  /* Parameters needed */
76  double M = 0.938; // nucleon mass in GeV
77  double ml = mlepton; // lepton mass in GeV
78 
79  double gamma2 = pow(( 2 * M * x ), 2) / Q2; // ?
80 
81  double R = get_parametrization_R(x, Q2); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
82 
83  //double D_num = (
84  // y * ( 2 - y ) * ( 1 + gamma2 * y / 2. )
85  // );
86  double D_num = (
87  y * ( 2 - y ) * ( 1 + gamma2 * y / 2. ) - 2 * pow(y, 3) * pow(ml, 2) / Q2
88  ); // Ernst's thesis has the extra y^3 term
89  double D_denom = (
90  pow(y, 2) * ( 1 + gamma2 ) * (1 - 2 * pow(ml, 2) / Q2 )
91  + 2 * (1 - y - gamma2 * pow(y, 2) / 4 ) * (1 + R )
92  );
93  double D = D_num / D_denom; // Depolarization factor
94 
95  return D;
96 }
97 
100 double
101 get_beam_polarization()
102 {
103  /* assumed beam polarizarion factors */
104  double pol_electron = 0.7;
105  double pol_proton = 0.7;
106  double pol_prod = pol_electron * pol_proton;
107 
108  return pol_prod;
109 }
110 
113 double
114 get_g1_value( double x_val, double Q2_val, double y_val, double mlepton=0.000511 )
115 {
116  /* F2 and R parametrizations (evaluated at kinmatics of given data point) */
117  double F2 = get_parametrization_F2(x_val, Q2_val);
118  double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
119 
120  /* Depolarization factor */
121  double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
122 
123  /* Asymmetry measured in parallel spin configuration */
124  // double A_parralel = D * A1; // assuming A2 = 0 // Ernst's thesis Eq. 3.41
125 
126  /* A_raw = (N++ - N+-) / (N++ + N+-) */
127  //double A_raw = ?
128  //double A_parallel = 1. / pol_prod * A_raw;
129 
130  double A1_Ernst = 0.033;
131  double A_parallel = D * A1_Ernst;
132 
133  double g1_val = 1. / (2 * x_val * (1 + R)) * F2 * A_parallel / D;
134 
135  // return g1_val;
136  return 0;
137 }
138 
141 double
142 get_g1_sigma( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
143 {
144  /* F2 and R parametrizations (evaluated at kinmatics of given data point) */
145  double F2 = get_parametrization_F2(x_val, Q2_val);
146  double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
147 
148  /* Depolarization factor */
149  double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
150 
151  /* Beam polarization */
152  double pol_prod = get_beam_polarization();
153 
154  double A_parallel_sig = (sqrt(N_val) / N_val) * (1. / pol_prod);
155  double g1_sig = 1. / (2 * x_val * (1 + R)) * F2 * A_parallel_sig / D;
156 
157  return g1_sig;
158 }
159 
162 double
163 get_A1_sigma( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
164 {
165  /* Calculate F1 */
166  double M = 0.938; // nucleon mass in GeV
167 
168  double gamma2 = pow(( 2 * M * x_val ), 2) / Q2_val; // ?
169 
170  double F2 = get_parametrization_F2(x_val, Q2_val);
171  double R = get_parametrization_R(x_val, Q2_val); // ratio of longitudinal and transverse photoabsorption cross section, R = simga_L / sigma_T
172 
173  double F1 = ( ( 1 + gamma2 ) / ( 2 * x_val * ( 1 + R ) ) ) * F2;
174 
175  /* Calculate A1 sigma (note: A1 = g1 / F1) */
176  double g1_sig = get_g1_sigma( x_val, Q2_val, y_val, N_val, mlepton );
177  double A1_sig = g1_sig / F1;
178 
179  return A1_sig;
180 }
181 
184 void
185 calculate_g1( double x_val, double Q2_val, double y_val, double N_val, double mlepton=0.000511 )
186 {
187  double g1_val = get_g1_value(x_val, Q2_val, y_val, mlepton);
188  double g1_sig = get_g1_sigma(x_val, Q2_val, y_val, N_val, mlepton);
189  double F2 = get_parametrization_F2(x_val, Q2_val);
190  double R = get_parametrization_R(x_val, Q2_val);
191  double D = get_depolarization_factor( x_val, Q2_val, y_val, mlepton );
192 
193  /* print parameters */
194  cout << "-----------------------------------------------------------------------" << endl;
195  cout << "x = " << x_val << " , Q2 = " << Q2_val << " , y = " << y_val << " , N = " << N_val << endl;
196  cout << "F2 = " << F2 << ", R = " << R << ", D = " << D << endl;
197  cout << "==> g1 = " << g1_val << " +/- " << g1_sig << endl;
198  cout << "-----------------------------------------------------------------------" << endl;
199 
200  return;
201 }