148 {
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;
152
153 int dedxflag = vflag[0];
154 int sigmaflag = vflag[1];
155 bool ifMC = false;
156 if ( vflag[2] == 1 ) ifMC = true;
157
158 int Nmax_prob( 0 );
159 float max_prob( -0.01 );
160 int ndf;
161 float chi2;
162
163 for ( int it = 0; it < par_cand; it++ )
164 {
165 beta_G = mom / Charge_Mass[it];
166
167 if ( dedxflag == 1 )
168 {
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];
172 }
173 else if ( dedxflag == 2 )
174 {
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 );
186
187 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
188 }
189
190 if ( ifMC )
191 {
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 );
203
204 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
205 }
206
207 if ( Nohit > 0 )
208 {
209 dedx_exp[it] = bethe_B;
210 double sig_the = std::sin( (double)theta );
211 double f_betagamma, g_sinth, h_nhit, i_t0;
212
213 if ( ifMC )
214 {
215
216 if ( sigmaflag >= 6 )
217 {
218
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;
223
224 double cor_nhit = 1.0;
225 if ( nhit < 5 ) nhit = 5;
226 if ( nhit <= 35 )
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];
229
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];
234
235
236 double cor_t0 = 1.0;
237
238
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];
241 }
242
243 else
244 {
246 double nhit = (double)Nohit;
247 double sigma_bg = 1.0;
248
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];
253
254 double cor_nhit = 1.0;
255 if ( nhit < 5 ) nhit = 5;
256 if ( nhit <= 35 )
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];
259
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];
264
265
266 double cor_t0 = 1.0;
267
268
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];
271 }
272 }
273 else
274 {
275 if ( sigmaflag == 1 )
276 {
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] );
280 h_nhit =
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] +
286 sig_par[12] );
287 else if ( sig_par[13] == 0 ) i_t0 = 1;
288
289
290
291 ex_sigma[it] = f_betagamma * g_sinth * h_nhit * i_t0;
292 }
293 else if ( sigmaflag == 2 )
294 {
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];
302
303 double cor_nhit = 1.0;
304 if ( nhit < 5 ) nhit = 5;
305 if ( nhit <= 35 )
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];
308
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];
313
314 double cor_t0 = 1;
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];
317
318 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
319 }
320 else if ( sigmaflag == 3 )
321 {
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];
329
330 double cor_nhit = 1.0;
331 if ( nhit < 5 ) nhit = 5;
332 if ( nhit <= 35 )
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];
335
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];
340
341 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
342 }
343 else if ( sigmaflag == 4 )
344 {
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];
352
353 double cor_nhit = 1.0;
354 if ( nhit < 5 ) nhit = 5;
355 if ( nhit <= 35 )
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];
358
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];
363
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];
366 }
367 else if ( sigmaflag == 5 )
368 {
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];
376
377 double cor_nhit = 1.0;
378 if ( nhit < 5 ) nhit = 5;
379 if ( nhit <= 35 )
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];
386
387 double cor_t0 = 1;
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];
390
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];
393 }
394 else if ( sigmaflag == 6 )
395 {
396 double x = bethe_B / 550;
397 double nhit = (double)Nohit;
398 double sigma_bg = 1.0;
399
400 sigma_bg = sig_par[0] + sig_par[1] *
x + sig_par[2] *
x *
x;
401
402 double cor_nhit = 1.0;
403 if ( nhit < 5 ) nhit = 5;
404 if ( nhit <= 35 )
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];
407
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];
412
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];
415 }
416 else if ( sigmaflag == 7 )
417 {
418
419 double x = bethe_B / 550;
420 double nhit = (double)Nohit;
421 double sigma_bg = 1.0;
422
423 sigma_bg = sig_par[0] + sig_par[1] *
x + sig_par[2] *
x *
x;
424
425 double cor_nhit = 1.0;
426 if ( nhit < 5 ) nhit = 5;
427 if ( nhit <= 35 )
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];
434
435 double cor_t0 = 1;
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];
439
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];
442
443
444
445
446 }
447 }
448
449 MdcDedxTrk dedxtrk;
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];
453 ndf = 1;
454 pid_prob[it] =
prob_( &chi2, &ndf );
455
456
457 if ( it == -999 )
458 {
459 std::cout << " mom = " << mom << "exp" << dedx_exp[it] << " sigma " << ex_sigma[it]
460 << " prob " << pid_prob[it] << std::endl;
461 }
462 if ( pid_prob[it] > max_prob )
463 {
464 max_prob = pid_prob[it];
465 Nmax_prob = it;
466 }
467 }
468 else
469 {
470 dedx_exp[it] = 0.0;
471 ex_sigma[it] = 1000.0;
472 pid_prob[it] = 0.0;
473 chi_dedx[it] = 999.0;
474 }
475 }
476
477}
float prob_(float *, int *)
EvtComplex exp(const EvtComplex &c)
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C