144 float mom,
float theta,
float t0,
float lsamp,
145 double dedx_exp[5],
double ex_sigma[5],
146 double pid_prob[5],
double chi_dedx[5],
147 std::vector<double>& par,
148 std::vector<double>& sig_par )
const {
149 const int par_cand( 5 );
150 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
151 double beta_G, beta, betterm, bethe_B = 0;
153 int dedxflag = vflag[0];
154 int sigmaflag = vflag[1];
156 if ( vflag[2] == 1 ) ifMC =
true;
159 float max_prob( -0.01 );
163 for (
int it = 0; it < par_cand; it++ )
165 beta_G = mom / Charge_Mass[it];
169 beta = beta_G / sqrt( 1 + ( beta_G ) * ( beta_G ) );
170 betterm = par[1] - log( par[2] + pow( 1 / beta_G, par[4] ) );
171 bethe_B = par[0] / pow( beta, par[3] ) * betterm - par[0];
173 else if ( dedxflag == 2 )
175 double A = 0, B = 0,
C = 0;
177 if ( x < 4.5 ) A = 1;
178 else if ( x < 10 ) B = 1;
180 double partA = par[0] * pow( sqrt( x * x + 1 ), par[2] ) / pow( x, par[2] ) *
181 ( par[1] - par[5] * log( pow( 1 / x, par[3] ) ) ) -
182 par[4] +
exp( par[6] + par[7] * x );
183 double partB = par[8] * pow( x, 3 ) + par[9] * pow( x, 2 ) + par[10] * x + par[11];
184 double partC = -par[12] * log( par[15] + pow( 1 / x, par[13] ) ) + par[14];
185 bethe_B = 550 * ( A * partA + B * partB +
C * partC );
187 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
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 dedx_exp[it] = bethe_B;
210 double sig_the = std::sin( (
double)theta );
211 double f_betagamma, g_sinth, h_nhit, i_t0;
216 if ( sigmaflag >= 6 )
219 double nhit = (double)Nohit;
220 double sigma_bg = 1.0;
221 double ndedx = bethe_B / 550;
222 sigma_bg = sig_par[0] + sig_par[1] * ndedx + sig_par[2] * ndedx * ndedx;
224 double cor_nhit = 1.0;
225 if ( nhit < 5 ) nhit = 5;
227 cor_nhit = sig_par[3] * pow( nhit, 4 ) + sig_par[4] * pow( nhit, 3 ) +
228 sig_par[5] * pow( nhit, 2 ) + sig_par[6] * nhit + sig_par[7];
230 double cor_sin = 1.0;
231 if ( sig_the > 0.99 ) sig_the = 0.99;
232 cor_sin = sig_par[8] * pow( sig_the, 4 ) + sig_par[9] * pow( sig_the, 3 ) +
233 sig_par[10] * pow( sig_the, 2 ) + sig_par[11] * sig_the + sig_par[12];
239 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
240 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[13];
246 double nhit = (double)Nohit;
247 double sigma_bg = 1.0;
250 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
251 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
252 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
254 double cor_nhit = 1.0;
255 if ( nhit < 5 ) nhit = 5;
257 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
258 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
260 double cor_sin = 1.0;
261 if ( sig_the > 0.99 ) sig_the = 0.99;
262 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
263 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
269 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
270 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
275 if ( sigmaflag == 1 )
277 f_betagamma = sig_par[0] * pow( beta_G, sig_par[1] ) + sig_par[2];
278 g_sinth = ( sig_par[3] * sig_the * sig_the + sig_par[4] ) /
279 ( sig_par[3] * sig_par[5] * sig_par[5] + sig_par[4] );
281 ( sig_par[6] * Nohit * Nohit + sig_par[7] * Nohit + sig_par[8] ) /
282 ( sig_par[6] * sig_par[9] * sig_par[9] + sig_par[7] * sig_par[9] + sig_par[8] );
283 if ( sig_par[13] != 0 )
284 i_t0 = ( sig_par[10] * t0 * t0 + sig_par[11] * t0 + sig_par[12] ) /
285 ( sig_par[10] * sig_par[13] * sig_par[13] + sig_par[11] * sig_par[13] +
287 else if ( sig_par[13] == 0 ) i_t0 = 1;
291 ex_sigma[it] = f_betagamma * g_sinth * h_nhit * i_t0;
293 else if ( sigmaflag == 2 )
296 double nhit = (double)Nohit;
297 double sigma_bg = 1.0;
299 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
300 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
301 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
303 double cor_nhit = 1.0;
304 if ( nhit < 5 ) nhit = 5;
306 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
307 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
309 double cor_sin = 1.0;
310 if ( sig_the > 0.99 ) sig_the = 0.99;
311 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
312 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
315 if ( t0 > 1200 ) t0 = 1200;
316 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
318 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
320 else if ( sigmaflag == 3 )
323 double nhit = (double)Nohit;
324 double sigma_bg = 1.0;
326 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
327 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
328 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
330 double cor_nhit = 1.0;
331 if ( nhit < 5 ) nhit = 5;
333 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
334 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
336 double cor_sin = 1.0;
337 if ( sig_the > 0.99 ) sig_the = 0.99;
338 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
339 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
341 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
343 else if ( sigmaflag == 4 )
346 double nhit = (double)Nohit;
347 double sigma_bg = 1.0;
349 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
350 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
351 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
353 double cor_nhit = 1.0;
354 if ( nhit < 5 ) nhit = 5;
356 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
357 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
359 double cor_sin = 1.0;
360 if ( sig_the > 0.99 ) sig_the = 0.99;
361 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
362 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
364 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
365 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * sig_par[18];
367 else if ( sigmaflag == 5 )
370 double nhit = (double)Nohit;
371 double sigma_bg = 1.0;
373 sigma_bg = sig_par[0] *
exp( sig_par[1] * x ) +
374 sig_par[2] *
exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
375 else sigma_bg = sig_par[5] *
exp( sig_par[6] * x ) + sig_par[7];
377 double cor_nhit = 1.0;
378 if ( nhit < 5 ) nhit = 5;
380 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
381 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
382 double cor_sin = 1.0;
383 if ( sig_the > 0.99 ) sig_the = 0.99;
384 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
385 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
388 if ( t0 > 1200 ) t0 = 1200;
389 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
391 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
392 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[21];
394 else if ( sigmaflag == 6 )
396 double x = bethe_B / 550;
397 double nhit = (double)Nohit;
398 double sigma_bg = 1.0;
400 sigma_bg = sig_par[0] + sig_par[1] * x + sig_par[2] * x * x;
402 double cor_nhit = 1.0;
403 if ( nhit < 5 ) nhit = 5;
405 cor_nhit = sig_par[3] * pow( nhit, 4 ) + sig_par[4] * pow( nhit, 3 ) +
406 sig_par[5] * pow( nhit, 2 ) + sig_par[6] * nhit + sig_par[7];
408 double cor_sin = 1.0;
409 if ( sig_the > 0.99 ) sig_the = 0.99;
410 cor_sin = sig_par[8] * pow( sig_the, 4 ) + sig_par[9] * pow( sig_the, 3 ) +
411 sig_par[10] * pow( sig_the, 2 ) + sig_par[11] * sig_the + sig_par[12];
413 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
414 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * sig_par[13];
416 else if ( sigmaflag == 7 )
419 double x = bethe_B / 550;
420 double nhit = (double)Nohit;
421 double sigma_bg = 1.0;
423 sigma_bg = sig_par[0] + sig_par[1] * x + sig_par[2] * x * x;
425 double cor_nhit = 1.0;
426 if ( nhit < 5 ) nhit = 5;
428 cor_nhit = sig_par[3] * pow( nhit, 4 ) + sig_par[4] * pow( nhit, 3 ) +
429 sig_par[5] * pow( nhit, 2 ) + sig_par[6] * nhit + sig_par[7];
430 double cor_sin = 1.0;
431 if ( sig_the > 0.99 ) sig_the = 0.99;
432 cor_sin = sig_par[8] * pow( sig_the, 4 ) + sig_par[9] * pow( sig_the, 3 ) +
433 sig_par[10] * pow( sig_the, 2 ) + sig_par[11] * sig_the + sig_par[12];
436 if ( t0 > 1200 ) t0 = 1200;
437 if ( t0 < 800 ) t0 = 800;
438 cor_t0 = sig_par[13] * pow( t0, 2 ) + sig_par[14] * t0 + sig_par[15];
440 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
441 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[16];
450 double dedx_correc = dedx;
451 chi_dedx[it] = ( dedx_correc - dedx_exp[it] ) / ex_sigma[it];
452 chi2 = chi_dedx[it] * chi_dedx[it];
454 pid_prob[it] =
prob_( &chi2, &ndf );
459 std::cout <<
" mom = " << mom <<
"exp" << dedx_exp[it] <<
" sigma " << ex_sigma[it]
460 <<
" prob " << pid_prob[it] << std::endl;
462 if ( pid_prob[it] > max_prob )
464 max_prob = pid_prob[it];
471 ex_sigma[it] = 1000.0;
473 chi_dedx[it] = 999.0;