170{
171 const int par_cand( 5 );
172 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
173 double beta_G, beta, betterm, bethe_B;
174
175 int dedxflag = vflag[0];
176 int sigmaflag = vflag[1];
177 bool ifMC = false;
178 if(vflag[2] == 1) ifMC = true;
179
180 int Nmax_prob(0);
181 float max_prob(-0.01);
182 int ndf;
183 float chi2;
184
185 for( int it = 0; it < par_cand; it++ ) {
186 beta_G = mom/Charge_Mass[it];
187
188 if(dedxflag == 1){
189 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
190 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
191 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
192 }
193 else if(dedxflag == 2) {
196 if(x<4.5)
198 else if(x<10)
200 else
202 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);
203 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*
x+par[11];
204 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
205 bethe_B = 550*(
A*partA+
B*partB+
C*partC);
206
207 if(beta_G> 100) bethe_B=550*1.0;
208 }
209
210 if (ifMC) {
213 if(x<4.5)
215 else if(x<10)
217 else
219 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);
220 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*
x+par[11];
221 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
222 bethe_B = 550*(
A*partA+
B*partB+
C*partC);
223
224 if(beta_G> 100) bethe_B=550*1.0;
225 }
226
227
228 if (Nohit > 0) {
229 dedx_exp[it] = bethe_B;
230 double sig_the = std::sin((double)theta);
231 double f_betagamma, g_sinth, h_nhit, i_t0;
232
233 if (ifMC) {
234
236 double nhit = (double)Nohit;
237 double sigma_bg=1.0;
238 if (x > 0.3)
239 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
240 else
241 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
242
243 double cor_nhit = 1.0;
244 if (nhit < 5) nhit = 5;
245 if (nhit <= 35)
246 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
247 sig_par[11]*nhit+sig_par[12];
248
249 double cor_sin= 1.0;
250 if(sig_the>0.99) sig_the=0.99;
251 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
252 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
253
254
255 double cor_t0 = 1.0;
256
257
258 if (trkalg == 1)
259 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
260 else
261 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
262 }
263 else {
264 if(sigmaflag == 1) {
265 f_betagamma=sig_par[0]*pow(beta_G,sig_par[1])+sig_par[2];
266 g_sinth=(sig_par[3]*sig_the*sig_the+sig_par[4])/(sig_par[3]*sig_par[5]*sig_par[5]+sig_par[4]);
267 h_nhit=(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]+sig_par[12]);
272 else if(sig_par[13] == 0)
273 i_t0 =1;
274
275
276 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
277 }
278 else if(sigmaflag == 2) {
280 double nhit = (double)Nohit;
281 double sigma_bg=1.0;
282 if (x > 0.3)
283 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
284 else
285 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
286
287 double cor_nhit=1.0;
288 if(nhit<5) nhit=5;
289 if (nhit <= 35)
290 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];
291
292 double cor_sin= 1.0;
293 if(sig_the>0.99) sig_the = 0.99;
294 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
295 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
296
297 double cor_t0 = 1;
298 if (t0 > 1200) t0 = 1200;
299 if (t0 > 800)
300 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
301
302 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
303 }
304 else if(sigmaflag == 3) {
306 double nhit = (double)Nohit;
307 double sigma_bg=1.0;
308 if (x > 0.3)
309 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
310 else
311 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
312
313 double cor_nhit = 1.0;
314 if (nhit < 5) nhit = 5;
315 if (nhit <= 35)
316 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];
317
318 double cor_sin= 1.0;
319 if(sig_the>0.99) sig_the=0.99;
320 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
321 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
322
323 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
324 }
325 else if(sigmaflag == 4) {
327 double nhit = (double)Nohit;
328 double sigma_bg=1.0;
329 if (x > 0.3)
330 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
331 else
332 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
333
334 double cor_nhit = 1.0;
335 if (nhit < 5) nhit = 5;
336 if (nhit <= 35)
337 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];
338
339 double cor_sin= 1.0;
340 if(sig_the>0.99) sig_the=0.99;
341 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
342 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
343
344 if(trkalg==1)
345 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
346 else
347 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
348 }
349 else if(sigmaflag == 5) {
351 double nhit = (double)Nohit;
352 double sigma_bg=1.0;
353 if (x > 0.3)
354 sigma_bg = sig_par[0]*
exp(sig_par[1]*x)+sig_par[2]*
exp(sig_par[3]*pow(x,0.25))+sig_par[4];
355 else
356 sigma_bg = sig_par[5]*
exp(sig_par[6]*x)+ sig_par[7];
357
358 double cor_nhit=1.0;
359 if(nhit<5) nhit=5;
360 if (nhit <= 35)
361 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
362 sig_par[11]*nhit+sig_par[12];
363 double cor_sin= 1.0;
364 if(sig_the>0.99) sig_the = 0.99;
365 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
366 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
367
368 double cor_t0 = 1;
369 if (t0 > 1200) t0 = 1200;
370 if (t0 > 800)
371 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
372
373 if(trkalg==1)
374 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
375 else
376 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
377 }
378 }
379
380 double dedx_correc = dedx;
381 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
382 chi2 = chi_dedx[it]*chi_dedx[it];
383 ndf=1;
384 pid_prob[it] =
prob_(&chi2,&ndf);
385
386
387 if (it == -999) {
388 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
389 <<
" sigma "<<ex_sigma[it] <<
" prob "<<pid_prob[it]
390 << std::endl;
391 }
392 if( pid_prob[it] > max_prob ){
393 max_prob = pid_prob[it];
394 Nmax_prob = it;
395 }
396 }
397 else{
398 dedx_exp[it] = 0.0;
399 ex_sigma[it] = 1000.0;
400 pid_prob[it] = 0.0;
401 chi_dedx[it] = 999.0;
402 }
403 }
404}
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
1.9.6