5#include "TGraphErrors.h"
7#include "DedxCalibAlg/DedxCalibParameters.h"
17double mylan(
double *x,
double *par)
19 double kk=(
x[0]-par[1])/par[2];
20 double fterm=kk+
exp(-1*kk);
21 double fitval=par[0]*
exp(par[3]*fterm);
28 double kk=(
x[0]-par[1])/par[2];
29 double fterm=kk+
exp(-1*kk);
30 double fitval=par[0]*
exp(-0.5*fterm);
36 double fitval = TMath::Landau(
x[0], par[0], par[1], kFALSE);
45 vavset_(&par[0], &par[1], &mode);
52 double delta, a,
b, c, d, fitval;
53 a = TMath::Sqrt(log(4.0));
54 b = (
x[0]-par[1])/par[2];
55 c = TMath::Power(par[3],2.0);
56 delta = (1+ TMath::SinH(par[3]*a)/a)*
b;
57 d = TMath::Power(log(
delta),2.0)/c+c;
58 fitval = par[0]*TMath::Exp(-0.5*d);
64 int Nohit,
float mom,
float theta,
float t0,
float lsamp,
65 double dedx_exp[5],
double ex_sigma[5],
66 double pid_prob[5],
double chi_dedx[5] )
68 double par[5], sigma_par[4], sigma_index_nhit, sigma_index_sin;
97 const int par_cand( 5 );
98 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
99 double beta_G, beta, betterm, bethe_B, sig_param;
102 float max_prob( -0.01 );
106 for(
int it = 0; it < par_cand; it++ ) {
107 beta_G = mom/Charge_Mass[it];
108 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
109 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
110 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
113 dedx_exp[it] = bethe_B;
114 double sig_the=std::sin( (
double)theta );
116 if(runflag <3 && runflag>0){
118 sig_param = 1.6*std::sin( (
double) theta )/(lsamp*double(Nohit));
119 ex_sigma[it] = 0.05*bethe_B*sqrt( 50.0*sig_param );
124 sig_param=sigma_par[1]+sigma_par[2]*std::pow(beta_G,sigma_par[3]);
126 sig_param= sigma_par[0];
129 sig_the=std::pow(sig_the,sigma_index_sin);
131 sig_n=35.0/double(Nohit);
132 sig_n=std::pow(sig_n,sigma_index_nhit);
133 ex_sigma[it]=sig_param*sig_the*sig_n;
139 else dedx_correc = dedx;
140 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
141 chi2 = chi_dedx[it]*chi_dedx[it];
143 pid_prob[it] =
prob_(&chi2,&ndf);
147 std::cout <<
" mom = " << mom <<
"exp"<< dedx_exp[it]
148 <<
" sigma "<<ex_sigma[it] <<
" prob "<<pid_prob[it]
151 if( pid_prob[it] > max_prob ) {
152 max_prob = pid_prob[it];
158 ex_sigma[it] = 1000.0;
160 chi_dedx[it] = 999.0;
166 int Nohit,
float mom,
float theta,
float t0,
float lsamp,
167 double dedx_exp[5],
double ex_sigma[5],
168 double pid_prob[5],
double chi_dedx[5],
169 std::vector<double> & par, std::vector<double> & sig_par)
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;
175 int dedxflag =
vflag[0];
176 int sigmaflag =
vflag[1];
178 if(
vflag[2] == 1) ifMC =
true;
181 float max_prob(-0.01);
185 for(
int it = 0; it < par_cand; it++ ) {
186 beta_G = mom/Charge_Mass[it];
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];
193 else if(dedxflag == 2) {
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);
207 if(beta_G> 100) bethe_B=550*1.0;
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);
224 if(beta_G> 100) bethe_B=550*1.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;
236 double nhit = (double)Nohit;
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];
241 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
243 double cor_nhit = 1.0;
244 if (nhit < 5) nhit = 5;
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];
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];
259 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
261 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
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]);
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)
276 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
278 else if(sigmaflag == 2) {
280 double nhit = (double)Nohit;
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];
285 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
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];
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];
298 if (t0 > 1200) t0 = 1200;
300 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
302 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
304 else if(sigmaflag == 3) {
306 double nhit = (double)Nohit;
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];
311 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
313 double cor_nhit = 1.0;
314 if (nhit < 5) nhit = 5;
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];
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];
323 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
325 else if(sigmaflag == 4) {
327 double nhit = (double)Nohit;
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];
332 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
334 double cor_nhit = 1.0;
335 if (nhit < 5) nhit = 5;
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];
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];
345 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
347 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
349 else if(sigmaflag == 5) {
351 double nhit = (double)Nohit;
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];
356 sigma_bg = sig_par[5]*
exp(sig_par[6]*
x)+ sig_par[7];
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];
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];
369 if (t0 > 1200) t0 = 1200;
371 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
374 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
376 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
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];
384 pid_prob[it] =
prob_(&chi2,&ndf);
388 std::cout <<
" mom = " << mom <<
"exp"<< dedx_exp[it]
389 <<
" sigma "<<ex_sigma[it] <<
" prob "<<pid_prob[it]
392 if( pid_prob[it] > max_prob ){
393 max_prob = pid_prob[it];
399 ex_sigma[it] = 1000.0;
401 chi_dedx[it] = 999.0;
408 const int par_cand( 5 );
409 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
411 double e_Par[5] = {143.349, 1.7315, 0.192616, 2.90437, 1.08248};
412 double Beta_Gamma[22] ={0.373026, 0.479605, 0.586184, 0.692763, 0.799342, 782.779, 1565.56,
413 2348.34, 17.2727, 18.1245, 1.43297, 2.14946, 12.1803, 13.6132, 6.62515, 10.4109,
414 14.1967, 17.9825, 21.7683, 26.0274, 30.7596, 35.4919 };
415 double K_par[22] ={4.64411e-05, 5.86544e-05, 8.05289e-05, 8.46981e-05, 8.92014e-05, 4.74517e-05,
416 4.51684e-05, 5.32732e-05, 6.12803e-05, 6.14592e-05, 8.08608e-05, 6.73184e-05, 5.46448e-05,
417 6.1377e-05, 6.57385e-05, 7.03053e-05, 6.61171e-05, 6.86824e-05, 6.246e-05, 7.25988e-05,
418 7.11034e-05, 6.24924e-05 };
419 double D_par[22] ={0.0871504, 0.0956379, 0.117193, 0.118647, 0.127203, 0.0566449, 0.0529198,
420 0.0642525, 0.0764562, 0.081341, 0.0952263, 0.0987536, 0.0639901, 0.0845994,0.0777062,
421 0.0823206, 0.0783874, 0.079537, 0.0815792, 0.0849875, 0.0824751,0.0776296 };
422 double DSqr_par[22] = {0.00759519, 0.0091466, 0.0137341, 0.0140772, 0.0161807, 0.00320864,
423 0.00280051, 0.00412839, 0.00584555, 0.00661636, 0.00906805, 0.00975227, 0.00409473,
424 0.00715706, 0.00603826, 0.00677668, 0.00614458, 0.00632613, 0.00665516, 0.00722288,
425 0.00680214, 0.00602635};
428 if(beta_G <0.3) beta_G =0.3;
429 double bet=beta_G/TMath::Sqrt(beta_G*beta_G+1);
430 double fterm=TMath::Log(e_Par[2]+1/pow(beta_G,e_Par[4]));
431 double fitval=e_Par[0]/pow(bet,e_Par[3])*(e_Par[1]-pow(bet,e_Par[3])-fterm);
432 TGraphErrors *gr1 =
new TGraphErrors(22,Beta_Gamma, K_par,0,0);
433 TGraphErrors *gr2 =
new TGraphErrors(22,Beta_Gamma, DSqr_par,0,0);
437 par[1] = gr1->Eval(m_theta);
438 par[2] = gr2->Eval(m_theta);
439 Double_t
y = fabs(
cos(m_theta));
440 double electron_par[3] ={334.032, 6.20658e-05, 0.00525673};
441 double arg= TMath::Sqrt(
y*
y+ par[2]);
443 double cal_factor =TMath::Exp(-(par[1]* par[0])/
arg);
444 double arg_electron = TMath::Sqrt(
y*
y + electron_par[2]);
446 double electron_factor = TMath::Exp(-(electron_par[1]* electron_par[0])/arg_electron);
448 double dedx_cal = dEdx/(cal_factor/electron_factor);
double Landau(double *x, double *par)
double Vavilov(double *x, double *par)
void dedx_pid_exp_old(int landau, int runflag, float dedx, int Nohit, float mom, float theta, float t0, float lsamp, double dedx_exp[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5])
double landaun(double *x, double *par)
double AsymGauss(double *x, double *par)
void dedx_pid_exp(int vflag[3], float dedx, int trkalg, int Nohit, float mom, float theta, float t0, float lsamp, double dedx_exp[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5], std::vector< double > &par, std::vector< double > &sig_par)
float prob_(float *, int *)
float vavset_(double *, double *, int *)
double mylan(double *x, double *par)
double SpaceChargeCorrec(double m_theta, double mom, int Particle, double dEdx)
EvtComplex exp(const EvtComplex &c)
double arg(const EvtComplex &c)
HepMC::GenParticle Particle
const double HV2_index_nhit
const double HV1_index_nhit
const double HV1_index_sin
const double HV2_index_sin
double cos(const BesAngle a)
***************************************************************************************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
int vflag[3]
Curve parameter vars.