55void dedx_pid_exp_old(
int landau,
int runflag,
float dedx,
int Nohit,
float mom,
float theta,
56 float t0,
float lsamp,
double dedx_exp[5],
double ex_sigma[5],
57 double pid_prob[5],
double chi_dedx[5] ) {
58 double par[5], sigma_par[4], sigma_index_nhit, sigma_index_sin;
74 else if ( runflag == 2 )
89 const int par_cand( 5 );
90 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
91 double beta_G, beta, betterm, bethe_B, sig_param;
94 float max_prob( -0.01 );
98 for (
int it = 0; it < par_cand; it++ )
100 beta_G = mom / Charge_Mass[it];
101 beta = beta_G / sqrt( 1 + ( beta_G ) * ( beta_G ) );
102 betterm = par[1] - log( par[2] + pow( 1 / beta_G, par[4] ) );
103 bethe_B = par[0] / pow( beta, par[3] ) * betterm - par[0];
107 dedx_exp[it] = bethe_B;
108 double sig_the = std::sin( (
double)theta );
110 if ( runflag < 3 && runflag > 0 )
114 sig_param = 1.6 * std::sin( (
double)theta ) / ( lsamp * double( Nohit ) );
115 ex_sigma[it] = 0.05 * bethe_B * sqrt( 50.0 * sig_param );
121 { sig_param = sigma_par[1] + sigma_par[2] * std::pow( beta_G, sigma_par[3] ); }
122 else { sig_param = sigma_par[0]; }
124 sig_the = std::pow( sig_the, sigma_index_sin );
126 sig_n = 35.0 / double( Nohit );
127 sig_n = std::pow( sig_n, sigma_index_nhit );
128 ex_sigma[it] = sig_param * sig_the * sig_n;
134 else dedx_correc = dedx;
135 chi_dedx[it] = ( dedx_correc - dedx_exp[it] ) / ex_sigma[it];
136 chi2 = chi_dedx[it] * chi_dedx[it];
138 pid_prob[it] =
prob_( &chi2, &ndf );
143 std::cout <<
" mom = " << mom <<
"exp" << dedx_exp[it] <<
" sigma " << ex_sigma[it]
144 <<
" prob " << pid_prob[it] << std::endl;
146 if ( pid_prob[it] > max_prob )
148 max_prob = pid_prob[it];
155 ex_sigma[it] = 1000.0;
157 chi_dedx[it] = 999.0;
162void dedx_pid_exp(
int vflag[3],
float dedx,
int trkalg,
int Nohit,
float mom,
float theta,
163 float t0,
float lsamp,
double dedx_exp[5],
double ex_sigma[5],
164 double pid_prob[5],
double chi_dedx[5], std::vector<double>& par,
165 std::vector<double>& sig_par ) {
166 const int par_cand( 5 );
167 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
168 double beta_G, beta, betterm, bethe_B;
170 int dedxflag = vflag[0];
171 int sigmaflag = vflag[1];
173 if ( vflag[2] == 1 ) ifMC =
true;
176 float max_prob( -0.01 );
180 for (
int it = 0; it < par_cand; it++ )
182 beta_G = mom / Charge_Mass[it];
186 beta = beta_G / sqrt( 1 + ( beta_G ) * ( beta_G ) );
187 betterm = par[1] - log( par[2] + pow( 1 / beta_G, par[4] ) );
188 bethe_B = par[0] / pow( beta, par[3] ) * betterm - par[0];
190 else if ( dedxflag == 2 )
192 double A = 0, B = 0,
C = 0;
194 if ( x < 4.5 ) A = 1;
195 else if ( x < 10 ) B = 1;
197 double partA = par[0] * pow( sqrt( x * x + 1 ), par[2] ) / pow( x, par[2] ) *
198 ( par[1] - par[5] * log( pow( 1 / x, par[3] ) ) ) -
199 par[4] +
exp( par[6] + par[7] * x );
200 double partB = par[8] * pow( x, 3 ) + par[9] * pow( x, 2 ) + par[10] * x + par[11];
201 double partC = -par[12] * log( par[15] + pow( 1 / x, par[13] ) ) + par[14];
202 bethe_B = 550 * ( A * partA + B * partB +
C * partC );
204 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
209 double A = 0, B = 0,
C = 0;
211 if ( x < 4.5 ) A = 1;
212 else if ( x < 10 ) B = 1;
214 double partA = par[0] * pow( sqrt( x * x + 1 ), par[2] ) / pow( x, par[2] ) *
215 ( par[1] - par[5] * log( pow( 1 / x, par[3] ) ) ) -
216 par[4] +
exp( par[6] + par[7] * x );
217 double partB = par[8] * pow( x, 3 ) + par[9] * pow( x, 2 ) + par[10] * x + par[11];
218 double partC = -par[12] * log( par[15] + pow( 1 / x, par[13] ) ) + par[14];
219 bethe_B = 550 * ( A * partA + B * partB +
C * partC );
221 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
226 dedx_exp[it] = bethe_B;
227 double sig_the = std::sin( (
double)theta );
228 double f_betagamma, g_sinth, h_nhit, i_t0;
234 double nhit = (double)Nohit;
235 double sigma_bg = 1.0;
237 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
238 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
239 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
241 double cor_nhit = 1.0;
242 if ( nhit < 5 ) nhit = 5;
244 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
245 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
247 double cor_sin = 1.0;
248 if ( sig_the > 0.99 ) sig_the = 0.99;
249 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
250 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
256 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
257 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
261 if ( sigmaflag == 1 )
263 f_betagamma = sig_par[0] * pow( beta_G, sig_par[1] ) + sig_par[2];
264 g_sinth = ( sig_par[3] * sig_the * sig_the + sig_par[4] ) /
265 ( sig_par[3] * sig_par[5] * sig_par[5] + sig_par[4] );
267 ( sig_par[6] * Nohit * Nohit + sig_par[7] * Nohit + sig_par[8] ) /
268 ( sig_par[6] * sig_par[9] * sig_par[9] + sig_par[7] * sig_par[9] + sig_par[8] );
269 if ( sig_par[13] != 0 )
270 i_t0 = ( sig_par[10] * t0 * t0 + sig_par[11] * t0 + sig_par[12] ) /
271 ( sig_par[10] * sig_par[13] * sig_par[13] + sig_par[11] * sig_par[13] +
273 else if ( sig_par[13] == 0 ) i_t0 = 1;
277 ex_sigma[it] = f_betagamma * g_sinth * h_nhit * i_t0;
279 else if ( sigmaflag == 2 )
282 double nhit = (double)Nohit;
283 double sigma_bg = 1.0;
285 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
286 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
287 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
289 double cor_nhit = 1.0;
290 if ( nhit < 5 ) nhit = 5;
292 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
293 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
295 double cor_sin = 1.0;
296 if ( sig_the > 0.99 ) sig_the = 0.99;
297 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
298 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
301 if ( t0 > 1200 ) t0 = 1200;
302 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
304 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
306 else if ( sigmaflag == 3 )
309 double nhit = (double)Nohit;
310 double sigma_bg = 1.0;
312 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
313 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
314 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
316 double cor_nhit = 1.0;
317 if ( nhit < 5 ) nhit = 5;
319 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
320 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
322 double cor_sin = 1.0;
323 if ( sig_the > 0.99 ) sig_the = 0.99;
324 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
325 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
327 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
329 else if ( sigmaflag == 4 )
332 double nhit = (double)Nohit;
333 double sigma_bg = 1.0;
335 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
336 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
337 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
339 double cor_nhit = 1.0;
340 if ( nhit < 5 ) nhit = 5;
342 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
343 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
345 double cor_sin = 1.0;
346 if ( sig_the > 0.99 ) sig_the = 0.99;
347 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
348 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
350 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
351 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * sig_par[18];
353 else if ( sigmaflag == 5 )
356 double nhit = (double)Nohit;
357 double sigma_bg = 1.0;
359 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
360 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
361 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
363 double cor_nhit = 1.0;
364 if ( nhit < 5 ) nhit = 5;
366 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
367 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
368 double cor_sin = 1.0;
369 if ( sig_the > 0.99 ) sig_the = 0.99;
370 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
371 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
374 if ( t0 > 1200 ) t0 = 1200;
375 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
377 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
378 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[21];
382 double dedx_correc = dedx;
383 chi_dedx[it] = ( dedx_correc - dedx_exp[it] ) / ex_sigma[it];
384 chi2 = chi_dedx[it] * chi_dedx[it];
386 pid_prob[it] =
prob_( &chi2, &ndf );
391 std::cout <<
" mom = " << mom <<
"exp" << dedx_exp[it] <<
" sigma " << ex_sigma[it]
392 <<
" prob " << pid_prob[it] << std::endl;
394 if ( pid_prob[it] > max_prob )
396 max_prob = pid_prob[it];
403 ex_sigma[it] = 1000.0;
405 chi_dedx[it] = 999.0;
411 const int par_cand( 5 );
412 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
414 double e_Par[5] = { 143.349, 1.7315, 0.192616, 2.90437, 1.08248 };
415 double Beta_Gamma[22] = { 0.373026, 0.479605, 0.586184, 0.692763, 0.799342, 782.779,
416 1565.56, 2348.34, 17.2727, 18.1245, 1.43297, 2.14946,
417 12.1803, 13.6132, 6.62515, 10.4109, 14.1967, 17.9825,
418 21.7683, 26.0274, 30.7596, 35.4919 };
419 double K_par[22] = { 4.64411e-05, 5.86544e-05, 8.05289e-05, 8.46981e-05, 8.92014e-05,
420 4.74517e-05, 4.51684e-05, 5.32732e-05, 6.12803e-05, 6.14592e-05,
421 8.08608e-05, 6.73184e-05, 5.46448e-05, 6.1377e-05, 6.57385e-05,
422 7.03053e-05, 6.61171e-05, 6.86824e-05, 6.246e-05, 7.25988e-05,
423 7.11034e-05, 6.24924e-05 };
424 double D_par[22] = { 0.0871504, 0.0956379, 0.117193, 0.118647, 0.127203, 0.0566449,
425 0.0529198, 0.0642525, 0.0764562, 0.081341, 0.0952263, 0.0987536,
426 0.0639901, 0.0845994, 0.0777062, 0.0823206, 0.0783874, 0.079537,
427 0.0815792, 0.0849875, 0.0824751, 0.0776296 };
428 double DSqr_par[22] = { 0.00759519, 0.0091466, 0.0137341, 0.0140772, 0.0161807,
429 0.00320864, 0.00280051, 0.00412839, 0.00584555, 0.00661636,
430 0.00906805, 0.00975227, 0.00409473, 0.00715706, 0.00603826,
431 0.00677668, 0.00614458, 0.00632613, 0.00665516, 0.00722288,
432 0.00680214, 0.00602635 };
434 beta_G = mom / Charge_Mass[
Particle];
435 if ( beta_G < 0.3 ) beta_G = 0.3;
436 double bet = beta_G / TMath::Sqrt( beta_G * beta_G + 1 );
437 double fterm = TMath::Log( e_Par[2] + 1 / pow( beta_G, e_Par[4] ) );
439 e_Par[0] / pow( bet, e_Par[3] ) * ( e_Par[1] - pow( bet, e_Par[3] ) - fterm );
440 TGraphErrors* gr1 =
new TGraphErrors( 22, Beta_Gamma, K_par, 0, 0 );
441 TGraphErrors* gr2 =
new TGraphErrors( 22, Beta_Gamma, DSqr_par, 0, 0 );
445 par[1] = gr1->Eval( m_theta );
446 par[2] = gr2->Eval( m_theta );
447 Double_t y = fabs(
cos( m_theta ) );
448 double electron_par[3] = { 334.032, 6.20658e-05, 0.00525673 };
449 double arg = TMath::Sqrt( y * y + par[2] );
451 double cal_factor = TMath::Exp( -( par[1] * par[0] ) /
arg );
452 double arg_electron = TMath::Sqrt( y * y + electron_par[2] );
455 double electron_factor = TMath::Exp( -( electron_par[1] * electron_par[0] ) / arg_electron );
457 double dedx_cal = dEdx / ( cal_factor / electron_factor );