143{
144 const int par_cand( 5 );
145 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
146 double beta_G, beta, betterm, bethe_B;
147
148 int dedxflag = vflag[0];
149 int sigmaflag = vflag[1];
150 bool ifMC = false;
151 if(vflag[2] == 1) ifMC = true;
152
153 int Nmax_prob(0);
154 float max_prob(-0.01);
155 int ndf;
156 float chi2;
157
158 for( int it = 0; it < par_cand; it++ ) {
159 beta_G = mom/Charge_Mass[it];
160
161 if(dedxflag == 1){
162 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
163 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
164 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
165 }
166 else if(dedxflag == 2) {
169 if(x<4.5)
171 else if(x<10)
173 else
175 double partA = par[0]*pow(sqrt(x*x+1),par[2])/pow(x,par[2])*(par[1]-par[5]*log(pow(1/x,par[3])) ) - par[4]+
exp(par[6]+par[7]*x);
176 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*
x+par[11];
177 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
178 bethe_B = 550*(
A*partA+
B*partB+
C*partC);
179
180 if(beta_G> 100) bethe_B=550*1.0;
181 }
182
183 if (ifMC) {
186 if(x<4.5)
188 else if(x<10)
190 else
192 double partA = par[0]*pow(sqrt(x*x+1),par[2])/pow(x,par[2])*(par[1]-par[5]*log(pow(1/x,par[3])) ) - par[4]+
exp(par[6]+par[7]*x);
193 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*
x+par[11];
194 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
195 bethe_B = 550*(
A*partA+
B*partB+
C*partC);
196
197 if(beta_G> 100) bethe_B=550*1.0;
198 }
199
200
201 if (Nohit > 0) {
202 dedx_exp[it] = bethe_B;
203 double sig_the = std::sin((double)theta);
204 double f_betagamma, g_sinth, h_nhit, i_t0;
205
206 if (ifMC) {
207
208 if(sigmaflag >= 6) {
209
210 double nhit = (double)Nohit;
211 double sigma_bg=1.0;
212 double ndedx=bethe_B/550;
213 sigma_bg= sig_par[0]+sig_par[1]*ndedx+sig_par[2]*ndedx*ndedx;
214
215 double cor_nhit = 1.0;
216 if (nhit < 5) nhit = 5;
217 if (nhit <= 35)
218 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
219 sig_par[6]*nhit+sig_par[7];
220
221 double cor_sin= 1.0;
222 if(sig_the>0.99) sig_the=0.99;
223 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
224 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
225
226
227 double cor_t0 = 1.0;
228
229
230 if (trkalg == 1)
231 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
232 else
233 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[13];
234
235 }
236
237 else{
239 double nhit = (double)Nohit;
240 double sigma_bg=1.0;
241
242 if (x > 0.3)
243 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
244 else
245 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
246
247
248 double cor_nhit = 1.0;
249 if (nhit < 5) nhit = 5;
250 if (nhit <= 35)
251 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
252 sig_par[11]*nhit+sig_par[12];
253
254 double cor_sin= 1.0;
255 if(sig_the>0.99) sig_the=0.99;
256 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
257 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
258
259
260 double cor_t0 = 1.0;
261
262
263 if (trkalg == 1)
264 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
265 else
266 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
267
268 }
269 }
270 else {
271 if(sigmaflag == 1) {
272 f_betagamma=sig_par[0]*pow(beta_G,sig_par[1])+sig_par[2];
273 g_sinth=(sig_par[3]*sig_the*sig_the+sig_par[4])/(sig_par[3]*sig_par[5]*sig_par[5]+sig_par[4]);
274 h_nhit=(sig_par[6]*Nohit*Nohit+sig_par[7]*Nohit+sig_par[8]) /
275 (sig_par[6]*sig_par[9]*sig_par[9]+sig_par[7]*sig_par[9]+sig_par[8]);
276 if(sig_par[13] != 0)
277 i_t0 = (sig_par[10]*t0*t0+sig_par[11]*t0+sig_par[12]) /
278 (sig_par[10]*sig_par[13]*sig_par[13]+sig_par[11]*sig_par[13]+sig_par[12]);
279 else if(sig_par[13] == 0)
280 i_t0 =1;
281
282
283 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
284 }
285 else if(sigmaflag == 2) {
287 double nhit = (double)Nohit;
288 double sigma_bg=1.0;
289 if (x > 0.3)
290 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
291 else
292 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
293
294 double cor_nhit=1.0;
295 if(nhit<5) nhit=5;
296 if (nhit <= 35)
297 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) + sig_par[11]*nhit+sig_par[12];
298
299 double cor_sin= 1.0;
300 if(sig_the>0.99) sig_the = 0.99;
301 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
302 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
303
304 double cor_t0 = 1;
305 if (t0 > 1200) t0 = 1200;
306 if (t0 > 800)
307 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
308
309 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
310 }
311 else if(sigmaflag == 3) {
313 double nhit = (double)Nohit;
314 double sigma_bg=1.0;
315 if (x > 0.3)
316 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
317 else
318 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
319
320 double cor_nhit = 1.0;
321 if (nhit < 5) nhit = 5;
322 if (nhit <= 35)
323 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) + sig_par[11]*nhit+sig_par[12];
324
325 double cor_sin= 1.0;
326 if(sig_the>0.99) sig_the=0.99;
327 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
328 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
329
330 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
331 }
332 else if(sigmaflag == 4) {
334 double nhit = (double)Nohit;
335 double sigma_bg=1.0;
336 if (x > 0.3)
337 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
338 else
339 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
340
341 double cor_nhit = 1.0;
342 if (nhit < 5) nhit = 5;
343 if (nhit <= 35)
344 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) + sig_par[11]*nhit+sig_par[12];
345
346 double cor_sin= 1.0;
347 if(sig_the>0.99) sig_the=0.99;
348 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
349 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
350
351 if(trkalg==1)
352 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
353 else
354 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
355 }
356 else if(sigmaflag == 5) {
358 double nhit = (double)Nohit;
359 double sigma_bg=1.0;
360 if (x > 0.3)
361 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
362 else
363 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
364
365 double cor_nhit=1.0;
366 if(nhit<5) nhit=5;
367 if (nhit <= 35)
368 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
369 sig_par[11]*nhit+sig_par[12];
370 double cor_sin= 1.0;
371 if(sig_the>0.99) sig_the = 0.99;
372 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
373 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
374
375 double cor_t0 = 1;
376 if (t0 > 1200) t0 = 1200;
377 if (t0 > 800)
378 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
379
380 if(trkalg==1)
381 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
382 else
383 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
384 }
385 else if(sigmaflag == 6) {
386 double x= bethe_B/550;
387 double nhit = (double)Nohit;
388 double sigma_bg=1.0;
389
390 sigma_bg= sig_par[0]+sig_par[1]*
x+sig_par[2]*
x*
x;
391
392 double cor_nhit = 1.0;
393 if (nhit < 5) nhit = 5;
394 if (nhit <= 35)
395 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) + sig_par[6]*nhit+sig_par[7];
396
397 double cor_sin= 1.0;
398 if(sig_the>0.99) sig_the=0.99;
399 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
400 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
401
402 if(trkalg==1)
403 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
404 else
405 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[13];
406 }
407 else if(sigmaflag == 7) {
408
409 double x= bethe_B/550;
410 double nhit = (double)Nohit;
411 double sigma_bg=1.0;
412
413 sigma_bg= sig_par[0]+sig_par[1]*
x+sig_par[2]*
x*
x;
414
415
416 double cor_nhit=1.0;
417 if(nhit<5) nhit=5;
418 if (nhit <= 35)
419 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
420 sig_par[6]*nhit+sig_par[7];
421 double cor_sin= 1.0;
422 if(sig_the>0.99) sig_the = 0.99;
423 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
424 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
425
426 double cor_t0 = 1;
427 if (t0 > 1200) t0 = 1200;
428 if (t0 < 800) t0 = 800;
429 cor_t0 = sig_par[13]*pow(t0,2)+sig_par[14]*t0+sig_par[15];
430
431 if(trkalg==1)
432 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
433 else
434 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[16];
435
436
437 }
438 }
439
441 double dedx_correc = dedx;
442 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
443 chi2 = chi_dedx[it]*chi_dedx[it];
444 ndf=1;
445 pid_prob[it] =
prob_(&chi2,&ndf);
446
447
448 if (it == -999) {
449 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
450 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
451 << std::endl;
452 }
453 if( pid_prob[it] > max_prob ){
454 max_prob = pid_prob[it];
455 Nmax_prob = it;
456 }
457 }
458 else{
459 dedx_exp[it] = 0.0;
460 ex_sigma[it] = 1000.0;
461 pid_prob[it] = 0.0;
462 chi_dedx[it] = 999.0;
463 }
464 }
465
466}
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