00001
00002
00003
00004
00005
00006 #include "MdcDedxAlg/MdcDedxCorrection.h"
00007 #include "MdcDedxAlg/MdcDedxParam.h"
00008 #include "MdcDedxAlg/MdcDedxTrk.h"
00009
00010 #include <vector>
00011 #include <cmath>
00012 #include <iostream>
00013
00014 using namespace std;
00015
00016 extern "C" {float prob_ (float *, int*);}
00017
00018 MdcDedxCorrection::MdcDedxCorrection()
00019 {
00020 #ifdef DEBUG
00021 std::cout<<"MdcDedxCorrection constructed!!!"<<std::endl;
00022 #endif
00023 }
00024
00025 MdcDedxCorrection::~MdcDedxCorrection()
00026 {
00027 #ifdef DEBUG
00028 std::cout<<"MdcDedxCorrection destructed!!!"<<std::endl;
00029 #endif
00030 }
00031
00032
00033
00034 void MdcDedxCorrection::dedx_pid_exp_old( int landau, int runflag, float dedx,
00035 int Nohit, float mom, float theta, float t0, float lsamp,
00036 double dedx_exp[5], double ex_sigma[5],
00037 double pid_prob[5], double chi_dedx[5] ) const
00038 {
00039 double par[5], sigma_par[4], sigma_index_nhit, sigma_index_sin;
00040
00041 if(runflag==1) {
00042 par[0]= MdcDedxParam::HV1_curvep0;
00043 par[1]= MdcDedxParam::HV1_curvep1;
00044 par[2]= MdcDedxParam::HV1_curvep2;
00045 par[3]= MdcDedxParam::HV1_curvep3;
00046 par[4]= MdcDedxParam::HV1_curvep4;
00047 sigma_par[0] = MdcDedxParam::HV1_sigmap0;
00048 sigma_par[1] = MdcDedxParam::HV1_sigmap1;
00049 sigma_par[2] = MdcDedxParam::HV1_sigmap2;
00050 sigma_par[3] = MdcDedxParam::HV1_sigmap3;
00051 sigma_index_nhit = MdcDedxParam::HV1_index_nhit;
00052 sigma_index_sin = MdcDedxParam::HV1_index_sin;
00053 }
00054 else if(runflag==2) {
00055 par[0]= MdcDedxParam::HV2_curvep0;
00056 par[1]= MdcDedxParam::HV2_curvep1;
00057 par[2]= MdcDedxParam::HV2_curvep2;
00058 par[3]= MdcDedxParam::HV2_curvep3;
00059 par[4]= MdcDedxParam::HV2_curvep4;
00060 sigma_par[0] = MdcDedxParam::HV2_sigmap0;
00061 sigma_par[1] = MdcDedxParam::HV2_sigmap1;
00062 sigma_par[2] = MdcDedxParam::HV2_sigmap2;
00063 sigma_par[3] = MdcDedxParam::HV2_sigmap3;
00064 sigma_index_nhit = MdcDedxParam::HV2_index_nhit;
00065 sigma_index_sin = MdcDedxParam::HV2_index_sin;
00066 }
00067
00068 const int par_cand( 5 );
00069 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
00070 double beta_G, beta, betterm, bethe_B, sig_param;
00071
00072 int Nmax_prob( 0 );
00073 float max_prob( -0.01 );
00074 int ndf;
00075 float chi2;
00076
00077 for( int it = 0; it < par_cand; it++ ) {
00078 beta_G = mom/Charge_Mass[it];
00079 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
00080 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
00081 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
00082
00083 if( Nohit >0 ) {
00084 dedx_exp[it] = bethe_B;
00085 double sig_the=std::sin( (double)theta );
00086
00087 if(runflag <3 && runflag>0){
00088 if(landau == 0) {
00089 sig_param = 1.6*std::sin( (double) theta )/(lsamp*double(Nohit));
00090 ex_sigma[it] = 0.05*bethe_B*sqrt( 50.0*sig_param );
00091 }
00092 else {
00093
00094 if(beta_G < 4) {
00095 sig_param=sigma_par[1]+sigma_par[2]*std::pow(beta_G,sigma_par[3]);
00096 } else {
00097 sig_param= sigma_par[0];
00098 }
00099
00100 sig_the=std::pow(sig_the,sigma_index_sin);
00101 double sig_n;
00102 sig_n=35.0/double(Nohit);
00103 sig_n=std::pow(sig_n,sigma_index_nhit);
00104 ex_sigma[it]=sig_param*sig_the*sig_n;
00105 }
00106 }
00107
00108 MdcDedxTrk dedxtrk;
00109 double dedx_correc;
00110 if(runflag == 2 ) dedx_correc = dedxtrk.SpaceChargeCorrec(theta, mom, it, dedx);
00111 else dedx_correc = dedx;
00112 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
00113 chi2 = chi_dedx[it]*chi_dedx[it];
00114 ndf = 1;
00115 pid_prob[it] = prob_(&chi2,&ndf);
00116
00117
00118 if( it == -999 ){
00119 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
00120 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
00121 << std::endl;
00122 }
00123 if( pid_prob[it] > max_prob ) {
00124 max_prob = pid_prob[it];
00125 Nmax_prob = it;
00126 }
00127 }
00128 else {
00129 dedx_exp[it] = 0.0;
00130 ex_sigma[it] = 1000.0;
00131 pid_prob[it] = 0.0;
00132 chi_dedx[it] = 999.0;
00133 }
00134 }
00135
00136 }
00137
00138 void MdcDedxCorrection::dedx_pid_exp(int vflag[3], float dedx, int trkalg,
00139 int Nohit, float mom, float theta, float t0, float lsamp,
00140 double dedx_exp[5], double ex_sigma[5],
00141 double pid_prob[5], double chi_dedx[5],
00142 std::vector<double> & par, std::vector<double> & sig_par) const
00143 {
00144 const int par_cand( 5 );
00145 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
00146 double beta_G, beta, betterm, bethe_B;
00147
00148 int dedxflag = vflag[0];
00149 int sigmaflag = vflag[1];
00150 bool ifMC = false;
00151 if(vflag[2] == 1) ifMC = true;
00152
00153 int Nmax_prob(0);
00154 float max_prob(-0.01);
00155 int ndf;
00156 float chi2;
00157
00158 for( int it = 0; it < par_cand; it++ ) {
00159 beta_G = mom/Charge_Mass[it];
00160
00161 if(dedxflag == 1){
00162 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
00163 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
00164 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
00165 }
00166 else if(dedxflag == 2) {
00167 double A=0, B=0,C=0;
00168 double x = beta_G;
00169 if(x<4.5)
00170 A=1;
00171 else if(x<10)
00172 B=1;
00173 else
00174 C=1;
00175 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);
00176 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
00177 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
00178 bethe_B = 550*(A*partA+B*partB+C*partC);
00179
00180 if(beta_G> 100) bethe_B=550*1.0;
00181 }
00182
00183 if (ifMC) {
00184 double A=0, B=0,C=0;
00185 double x = beta_G;
00186 if(x<4.5)
00187 A=1;
00188 else if(x<10)
00189 B=1;
00190 else
00191 C=1;
00192 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);
00193 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
00194 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
00195 bethe_B = 550*(A*partA+B*partB+C*partC);
00196
00197 if(beta_G> 100) bethe_B=550*1.0;
00198 }
00199
00200
00201 if (Nohit > 0) {
00202 dedx_exp[it] = bethe_B;
00203 double sig_the = std::sin((double)theta);
00204 double f_betagamma, g_sinth, h_nhit, i_t0;
00205
00206 if (ifMC) {
00207
00208 if(sigmaflag >= 6) {
00209
00210 double nhit = (double)Nohit;
00211 double sigma_bg=1.0;
00212 double ndedx=bethe_B/550;
00213 sigma_bg= sig_par[0]+sig_par[1]*ndedx+sig_par[2]*ndedx*ndedx;
00214
00215 double cor_nhit = 1.0;
00216 if (nhit < 5) nhit = 5;
00217 if (nhit <= 35)
00218 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
00219 sig_par[6]*nhit+sig_par[7];
00220
00221 double cor_sin= 1.0;
00222 if(sig_the>0.99) sig_the=0.99;
00223 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
00224 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
00225
00226
00227 double cor_t0 = 1.0;
00228
00229
00230 if (trkalg == 1)
00231 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
00232 else
00233 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[13];
00234
00235 }
00236
00237 else{
00238 double x= beta_G;
00239 double nhit = (double)Nohit;
00240 double sigma_bg=1.0;
00241
00242 if (x > 0.3)
00243 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00244 else
00245 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00246
00247
00248 double cor_nhit = 1.0;
00249 if (nhit < 5) nhit = 5;
00250 if (nhit <= 35)
00251 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00252 sig_par[11]*nhit+sig_par[12];
00253
00254 double cor_sin= 1.0;
00255 if(sig_the>0.99) sig_the=0.99;
00256 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00257 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00258
00259
00260 double cor_t0 = 1.0;
00261
00262
00263 if (trkalg == 1)
00264 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
00265 else
00266 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
00267
00268 }
00269 }
00270 else {
00271 if(sigmaflag == 1) {
00272 f_betagamma=sig_par[0]*pow(beta_G,sig_par[1])+sig_par[2];
00273 g_sinth=(sig_par[3]*sig_the*sig_the+sig_par[4])/(sig_par[3]*sig_par[5]*sig_par[5]+sig_par[4]);
00274 h_nhit=(sig_par[6]*Nohit*Nohit+sig_par[7]*Nohit+sig_par[8]) /
00275 (sig_par[6]*sig_par[9]*sig_par[9]+sig_par[7]*sig_par[9]+sig_par[8]);
00276 if(sig_par[13] != 0)
00277 i_t0 = (sig_par[10]*t0*t0+sig_par[11]*t0+sig_par[12]) /
00278 (sig_par[10]*sig_par[13]*sig_par[13]+sig_par[11]*sig_par[13]+sig_par[12]);
00279 else if(sig_par[13] == 0)
00280 i_t0 =1;
00281
00282
00283 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
00284 }
00285 else if(sigmaflag == 2) {
00286 double x = beta_G;
00287 double nhit = (double)Nohit;
00288 double sigma_bg=1.0;
00289 if (x > 0.3)
00290 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00291 else
00292 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00293
00294 double cor_nhit=1.0;
00295 if(nhit<5) nhit=5;
00296 if (nhit <= 35)
00297 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];
00298
00299 double cor_sin= 1.0;
00300 if(sig_the>0.99) sig_the = 0.99;
00301 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00302 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00303
00304 double cor_t0 = 1;
00305 if (t0 > 1200) t0 = 1200;
00306 if (t0 > 800)
00307 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
00308
00309 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
00310 }
00311 else if(sigmaflag == 3) {
00312 double x= beta_G;
00313 double nhit = (double)Nohit;
00314 double sigma_bg=1.0;
00315 if (x > 0.3)
00316 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00317 else
00318 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00319
00320 double cor_nhit = 1.0;
00321 if (nhit < 5) nhit = 5;
00322 if (nhit <= 35)
00323 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];
00324
00325 double cor_sin= 1.0;
00326 if(sig_the>0.99) sig_the=0.99;
00327 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00328 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00329
00330 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
00331 }
00332 else if(sigmaflag == 4) {
00333 double x= beta_G;
00334 double nhit = (double)Nohit;
00335 double sigma_bg=1.0;
00336 if (x > 0.3)
00337 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00338 else
00339 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00340
00341 double cor_nhit = 1.0;
00342 if (nhit < 5) nhit = 5;
00343 if (nhit <= 35)
00344 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];
00345
00346 double cor_sin= 1.0;
00347 if(sig_the>0.99) sig_the=0.99;
00348 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00349 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00350
00351 if(trkalg==1)
00352 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
00353 else
00354 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
00355 }
00356 else if(sigmaflag == 5) {
00357 double x = beta_G;
00358 double nhit = (double)Nohit;
00359 double sigma_bg=1.0;
00360 if (x > 0.3)
00361 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00362 else
00363 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00364
00365 double cor_nhit=1.0;
00366 if(nhit<5) nhit=5;
00367 if (nhit <= 35)
00368 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00369 sig_par[11]*nhit+sig_par[12];
00370 double cor_sin= 1.0;
00371 if(sig_the>0.99) sig_the = 0.99;
00372 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00373 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00374
00375 double cor_t0 = 1;
00376 if (t0 > 1200) t0 = 1200;
00377 if (t0 > 800)
00378 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
00379
00380 if(trkalg==1)
00381 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
00382 else
00383 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
00384 }
00385 else if(sigmaflag == 6) {
00386 double x= bethe_B/550;
00387 double nhit = (double)Nohit;
00388 double sigma_bg=1.0;
00389
00390 sigma_bg= sig_par[0]+sig_par[1]*x+sig_par[2]*x*x;
00391
00392 double cor_nhit = 1.0;
00393 if (nhit < 5) nhit = 5;
00394 if (nhit <= 35)
00395 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];
00396
00397 double cor_sin= 1.0;
00398 if(sig_the>0.99) sig_the=0.99;
00399 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
00400 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
00401
00402 if(trkalg==1)
00403 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
00404 else
00405 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[13];
00406 }
00407 else if(sigmaflag == 7) {
00408
00409 double x= bethe_B/550;
00410 double nhit = (double)Nohit;
00411 double sigma_bg=1.0;
00412
00413 sigma_bg= sig_par[0]+sig_par[1]*x+sig_par[2]*x*x;
00414
00415
00416 double cor_nhit=1.0;
00417 if(nhit<5) nhit=5;
00418 if (nhit <= 35)
00419 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
00420 sig_par[6]*nhit+sig_par[7];
00421 double cor_sin= 1.0;
00422 if(sig_the>0.99) sig_the = 0.99;
00423 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
00424 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
00425
00426 double cor_t0 = 1;
00427 if (t0 > 1200) t0 = 1200;
00428 if (t0 < 800) t0 = 800;
00429 cor_t0 = sig_par[13]*pow(t0,2)+sig_par[14]*t0+sig_par[15];
00430
00431 if(trkalg==1)
00432 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
00433 else
00434 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[16];
00435
00436
00437 }
00438 }
00439
00440 MdcDedxTrk dedxtrk;
00441 double dedx_correc = dedx;
00442 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
00443 chi2 = chi_dedx[it]*chi_dedx[it];
00444 ndf=1;
00445 pid_prob[it] = prob_(&chi2,&ndf);
00446
00447
00448 if (it == -999) {
00449 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
00450 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
00451 << std::endl;
00452 }
00453 if( pid_prob[it] > max_prob ){
00454 max_prob = pid_prob[it];
00455 Nmax_prob = it;
00456 }
00457 }
00458 else{
00459 dedx_exp[it] = 0.0;
00460 ex_sigma[it] = 1000.0;
00461 pid_prob[it] = 0.0;
00462 chi_dedx[it] = 999.0;
00463 }
00464 }
00465
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517