6#include "MdcDedxAlg/MdcDedxCorrection.h"
7#include "MdcDedxAlg/MdcDedxParam.h"
8#include "MdcDedxAlg/MdcDedxTrk.h"
16extern "C" {
float prob_ (
float *,
int*);}
21 std::cout<<
"MdcDedxCorrection constructed!!!"<<std::endl;
28 std::cout<<
"MdcDedxCorrection destructed!!!"<<std::endl;
35 int Nohit,
float mom,
float theta,
float t0,
float lsamp,
36 double dedx_exp[5],
double ex_sigma[5],
37 double pid_prob[5],
double chi_dedx[5] )
const
39 double par[5], sigma_par[4], sigma_index_nhit, sigma_index_sin;
68 const int par_cand( 5 );
69 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
70 double beta_G, beta, betterm, bethe_B, sig_param;
73 float max_prob( -0.01 );
77 for(
int it = 0; it < par_cand; it++ ) {
78 beta_G = mom/Charge_Mass[it];
79 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
80 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
81 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
84 dedx_exp[it] = bethe_B;
85 double sig_the=std::sin( (
double)theta );
87 if(runflag <3 && runflag>0){
89 sig_param = 1.6*std::sin( (
double) theta )/(lsamp*double(Nohit));
90 ex_sigma[it] = 0.05*bethe_B*sqrt( 50.0*sig_param );
95 sig_param=sigma_par[1]+sigma_par[2]*std::pow(beta_G,sigma_par[3]);
97 sig_param= sigma_par[0];
100 sig_the=std::pow(sig_the,sigma_index_sin);
102 sig_n=35.0/double(Nohit);
103 sig_n=std::pow(sig_n,sigma_index_nhit);
104 ex_sigma[it]=sig_param*sig_the*sig_n;
111 else dedx_correc = dedx;
112 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
113 chi2 = chi_dedx[it]*chi_dedx[it];
115 pid_prob[it] =
prob_(&chi2,&ndf);
119 std::cout <<
" mom = " << mom <<
"exp"<< dedx_exp[it]
120 <<
" sigma "<<ex_sigma[it] <<
" prob "<<pid_prob[it]
123 if( pid_prob[it] > max_prob ) {
124 max_prob = pid_prob[it];
130 ex_sigma[it] = 1000.0;
132 chi_dedx[it] = 999.0;
139 int Nohit,
float mom,
float theta,
float t0,
float lsamp,
140 double dedx_exp[5],
double ex_sigma[5],
141 double pid_prob[5],
double chi_dedx[5],
142 std::vector<double> & par, std::vector<double> & sig_par)
const
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;
148 int dedxflag = vflag[0];
149 int sigmaflag = vflag[1];
151 if(vflag[2] == 1) ifMC =
true;
154 float max_prob(-0.01);
158 for(
int it = 0; it < par_cand; it++ ) {
159 beta_G = mom/Charge_Mass[it];
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];
166 else if(dedxflag == 2) {
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);
180 if(beta_G> 100) bethe_B=550*1.0;
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);
197 if(beta_G> 100) bethe_B=550*1.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;
210 double nhit = (double)Nohit;
212 double ndedx=bethe_B/550;
213 sigma_bg= sig_par[0]+sig_par[1]*ndedx+sig_par[2]*ndedx*ndedx;
215 double cor_nhit = 1.0;
216 if (nhit < 5) nhit = 5;
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];
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];
231 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
233 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[13];
239 double nhit = (double)Nohit;
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];
245 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
248 double cor_nhit = 1.0;
249 if (nhit < 5) nhit = 5;
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];
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];
264 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
266 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
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]);
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)
283 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
285 else if(sigmaflag == 2) {
287 double nhit = (double)Nohit;
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];
292 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
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];
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];
305 if (t0 > 1200) t0 = 1200;
307 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
309 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
311 else if(sigmaflag == 3) {
313 double nhit = (double)Nohit;
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];
318 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
320 double cor_nhit = 1.0;
321 if (nhit < 5) nhit = 5;
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];
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];
330 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
332 else if(sigmaflag == 4) {
334 double nhit = (double)Nohit;
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];
339 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
341 double cor_nhit = 1.0;
342 if (nhit < 5) nhit = 5;
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];
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];
352 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
354 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
356 else if(sigmaflag == 5) {
358 double nhit = (double)Nohit;
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];
363 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
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];
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];
376 if (t0 > 1200) t0 = 1200;
378 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
381 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
383 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
385 else if(sigmaflag == 6) {
386 double x= bethe_B/550;
387 double nhit = (double)Nohit;
390 sigma_bg= sig_par[0]+sig_par[1]*
x+sig_par[2]*
x*
x;
392 double cor_nhit = 1.0;
393 if (nhit < 5) nhit = 5;
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];
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];
403 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
405 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[13];
407 else if(sigmaflag == 7) {
409 double x= bethe_B/550;
410 double nhit = (double)Nohit;
413 sigma_bg= sig_par[0]+sig_par[1]*
x+sig_par[2]*
x*
x;
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];
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];
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];
432 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
434 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[16];
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];
445 pid_prob[it] =
prob_(&chi2,&ndf);
449 std::cout <<
" mom = " << mom <<
"exp"<< dedx_exp[it]
450 <<
" sigma "<<ex_sigma[it] <<
" prob "<<pid_prob[it]
453 if( pid_prob[it] > max_prob ){
454 max_prob = pid_prob[it];
460 ex_sigma[it] = 1000.0;
462 chi_dedx[it] = 999.0;
float prob_(float *, int *)
EvtComplex exp(const EvtComplex &c)
float prob_(float *, int *)
***************************************************************************************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
void dedx_pid_exp_old(int landau, int runflag, float dedx, int Nhit, float mom, float theta, float t0, float lsamp, double dedx_ex[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5]) const
void dedx_pid_exp(int vflag[3], float dedx, int trkalg, int Nhit, float mom, float theta, float t0, float lsamp, double dedx_ex[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5], std::vector< double > &DedxCurve_Parameter, std::vector< double > &DedxSigma_Parameter) const
static double HV2_index_nhit
static double HV2_curvep3
static double HV1_curvep4
static double HV2_index_sin
static double HV1_sigmap2
static double HV1_curvep0
static double HV1_curvep1
static double HV2_curvep4
static double HV2_sigmap1
static double HV1_index_sin
static double HV2_sigmap0
static double HV2_curvep1
static double HV1_index_nhit
static double HV2_curvep2
static double HV2_curvep0
static double HV2_sigmap3
static double HV1_curvep3
static double HV1_sigmap1
static double HV1_sigmap0
static double HV1_sigmap3
static double HV2_sigmap2
static double HV1_curvep2
double SpaceChargeCorrec(double, double, int, double)