#include <MdcDedxCorrection.h>
Public Member Functions | |
| 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 |
| 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 |
| void | dedx_pid_exp_old (int landau, int alg, 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_old (int landau, int alg, 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 | Dump_calibdata (void) const |
| void | Dump_calibdata (void) const |
| void | getCalib (void) const |
| void | getCalib (void) const |
Static Public Member Functions | |
| MdcDedxCorrection * | getMdcDedxCorrection () |
| MdcDedxCorrection * | getMdcDedxCorrection () |
Protected Member Functions | |
| MdcDedxCorrection () | |
| constructer | |
| MdcDedxCorrection () | |
| constructer | |
| virtual | ~MdcDedxCorrection () |
| destructor | |
| virtual | ~MdcDedxCorrection () |
| destructor | |
Static Private Attributes | |
| MdcDedxCorrection * | m_excalib |
| MdcDedxCorrection * | m_excalib = 0 |
|
|
constructer
00047 {
00048 std::cout<<"MdcDedxCorrection constructed!!!"<<std::endl;
00049 }
|
|
|
destructor
00051 {
00052 std::cout<<"MdcDedxCorrection destructed!!!"<<std::endl;}
|
|
|
constructer
|
|
|
destructor
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
00189 {
00190 //cout<<"Nohit = "<<Nohit<<endl;
00191
00192 const int par_cand( 5 );
00193 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
00194 const float mass_e( 0.511 );
00195 double beta_G, beta, bet_f, Xa, deci;
00196 float betterm, bethe_B;
00197
00198 int dedxflag = vflag[0];
00199 int sigmaflag = vflag[1];
00200 bool ifMC = false;
00201 if(vflag[2] == 1)
00202 ifMC = true;
00203
00204
00205 float prob(0.0);
00206 double delta;
00207 double sig_param;
00208 // for debag
00209 int Nmax_prob(0);
00210 float max_prob(-0.01);
00211 int ndf;
00212 float prob_d[5], chi2;
00213
00214 for( int it = 0; it < par_cand; it++ ) {
00215 beta_G = mom/Charge_Mass[it];
00216
00217 if(dedxflag == 1){
00218 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
00219 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
00220 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
00221 }
00222 else if(dedxflag == 2) {
00223 double A=0, B=0,C=0;
00224 double x = beta_G;
00225 if(x<4.5)
00226 A=1;
00227 else if(x<10)
00228 B=1;
00229 else
00230 C=1;
00231 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])) ) -
00232 par[4]+exp(par[6]+par[7]*x);
00233 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
00234 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
00235 bethe_B = 550*(A*partA+B*partB+C*partC);
00236 //for fermi plateau ( the electron region) we just use 1.0
00237 if(beta_G> 100)
00238 bethe_B=550*1.0;
00239
00240 }
00241
00242 if (ifMC) {
00243 double A=0, B=0,C=0;
00244 double x = beta_G;
00245 if(x<4.5)
00246 A=1;
00247 else if(x<10)
00248 B=1;
00249 else
00250 C=1;
00251 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])) ) -
00252 par[4]+exp(par[6]+par[7]*x);
00253 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
00254 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
00255 bethe_B = 550*(A*partA+B*partB+C*partC);
00256 //for fermi plateau ( the electron region) we just use 1.0
00257 if(beta_G> 100)
00258 bethe_B=550*1.0;
00259 }
00260
00261
00262 if (Nohit > 0) {
00263 dedx_exp[it] = bethe_B;
00264 double sig_the = std::sin((double)theta);
00265 //cout<<"nhit = "<<Nohit<<endl;
00266 double f_betagamma, g_sinth, h_nhit, i_t0;
00267
00268 if (ifMC) {
00269
00270 double x= beta_G;
00271 double nhit = (double)Nohit;
00272 double sigma_bg=1.0;
00273 if (x > 0.3)
00274 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00275 else
00276 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00277
00278 double cor_nhit = 1.0;
00279 if (nhit < 5) nhit = 5;
00280 if (nhit <= 35)
00281 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00282 sig_par[11]*nhit+sig_par[12];
00283
00284 double cor_sin= 1.0;
00285 if(sig_the>0.99) sig_the=0.99;
00286 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00287 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00288
00289 //sigma vs t0
00290 double cor_t0 = 1.0;
00291
00292 //calculate sigma
00293 if (trkalg == 1)
00294 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
00295 else
00296 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
00297 }
00298 else {
00299 if(sigmaflag == 1) {
00300 f_betagamma=sig_par[0]*pow(beta_G,sig_par[1])+sig_par[2];
00301 g_sinth=(sig_par[3]*sig_the*sig_the+sig_par[4])/(sig_par[3]*sig_par[5]*sig_par[5]+sig_par[4]);
00302 h_nhit=(sig_par[6]*Nohit*Nohit+sig_par[7]*Nohit+sig_par[8]) /
00303 (sig_par[6]*sig_par[9]*sig_par[9]+sig_par[7]*sig_par[9]+sig_par[8]);
00304 if(sig_par[13] != 0)
00305 i_t0 = (sig_par[10]*t0*t0+sig_par[11]*t0+sig_par[12]) /
00306 (sig_par[10]*sig_par[13]*sig_par[13]+sig_par[11]*sig_par[13]+sig_par[12]);
00307 else if(sig_par[13] == 0)
00308 i_t0 =1;
00309 //cout<<"f_betagamma : "<<f_betagamma<<" g_sinth : "<<g_sinth<<" h_nhit : "
00310 //<<h_nhit<<" i_t0 : "<<i_t0<<endl;
00311 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
00312 }
00313 else if(sigmaflag == 2) {
00314 double x = beta_G;
00315 double nhit = (double)Nohit;
00316 double sigma_bg=1.0;
00317 if (x > 0.3)
00318 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00319 else
00320 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00321
00322 double cor_nhit=1.0;
00323 if(nhit<5) nhit=5;
00324 if (nhit <= 35)
00325 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00326 sig_par[11]*nhit+sig_par[12];
00327
00328 double cor_sin= 1.0;
00329 if(sig_the>0.99) sig_the = 0.99;
00330 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00331 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00332
00333 double cor_t0 = 1;
00334 if (t0 > 1200) t0 = 1200;
00335 if (t0 > 800)
00336 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
00337
00338 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
00339 }
00340 else if(sigmaflag == 3) {
00341 double x= beta_G;
00342 double nhit = (double)Nohit;
00343 double sigma_bg=1.0;
00344 if (x > 0.3)
00345 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00346 else
00347 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00348
00349 double cor_nhit = 1.0;
00350 if (nhit < 5) nhit = 5;
00351 if (nhit <= 35)
00352 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00353 sig_par[11]*nhit+sig_par[12];
00354
00355 double cor_sin= 1.0;
00356 if(sig_the>0.99) sig_the=0.99;
00357 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00358 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00359
00360
00361 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
00362 }
00363 else if(sigmaflag == 4) {
00364 double x= beta_G;
00365 double nhit = (double)Nohit;
00366 double sigma_bg=1.0;
00367 if (x > 0.3)
00368 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00369 else
00370 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00371
00372 double cor_nhit = 1.0;
00373 if (nhit < 5) nhit = 5;
00374 if (nhit <= 35)
00375 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00376 sig_par[11]*nhit+sig_par[12];
00377
00378 double cor_sin= 1.0;
00379 if(sig_the>0.99) sig_the=0.99;
00380 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00381 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00382
00383
00384 if(trkalg==1)
00385 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
00386 else
00387 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
00388 }
00389 else if(sigmaflag == 5) {
00390 double x = beta_G;
00391 double nhit = (double)Nohit;
00392 double sigma_bg=1.0;
00393 if (x > 0.3)
00394 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
00395 else
00396 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
00397
00398 double cor_nhit=1.0;
00399 if(nhit<5) nhit=5;
00400 if (nhit <= 35)
00401 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
00402 sig_par[11]*nhit+sig_par[12];
00403
00404 double cor_sin= 1.0;
00405 if(sig_the>0.99) sig_the = 0.99;
00406 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
00407 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
00408
00409 double cor_t0 = 1;
00410 if (t0 > 1200) t0 = 1200;
00411 if (t0 > 800)
00412 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
00413
00414 if(trkalg==1)
00415 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
00416 else
00417 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
00418 }
00419 }
00420
00421
00422
00423 MdcDedxTrk dedxtrk;
00424 double dedx_correc;
00425 //dedxtrk.clears();
00426 //chi_dedx[it] = (dedx - dedx_exp[it])/ex_sigma[it];
00427 //cout<<"dedx = "<<dedx<< " dedx_correc = "<<dedx_correc<<endl;
00428 dedx_correc = dedx;
00429 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
00430 chi2 = chi_dedx[it]*chi_dedx[it];
00431 ndf=1;
00432 pid_prob[it] = prob_(&chi2,&ndf);
00433 //if(it ==0 ) cout<<"runflag: "<<runflag<<" dedx : "<<dedx<<" chi_dedx: "
00434 //<<chi_dedx[it] <<" ptrk: "<<mom<<endl;
00435 if (it == -999) { // here a debug flag
00436 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
00437 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
00438 << std::endl;
00439 }
00440 if( pid_prob[it] > max_prob ){
00441 max_prob = pid_prob[it];
00442 Nmax_prob = it;
00443 }
00444 }
00445 else{
00446 dedx_exp[it] = 0.0;
00447 ex_sigma[it] = 1000.0;
00448 pid_prob[it] = 0.0;
00449 chi_dedx[it] = 999.0;
00450 } //if Nohit > 0
00451 } //for
00452 // std::cout<<"MdcDedxCorrection::dedx_pid_exp(blum)!!!"<<std::endl;
00453 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
00062 {
00063 //cout<<"Nohit = "<<Nohit<<endl;
00064
00065 const int par_cand( 5 );
00066 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
00067 const float mass_e( 0.511 );
00068 double beta_G, beta, bet_f, Xa, deci;
00069 float betterm, bethe_B;
00070 double par[5], sig_par[14], sigma_par[4], sigma_index_nhit, sigma_index_sin;
00071
00072 if(runflag==1) {
00073 par[0]= MdcDedxParam::HV1_curvep0;
00074 par[1]= MdcDedxParam::HV1_curvep1;
00075 par[2]= MdcDedxParam::HV1_curvep2;
00076 par[3]= MdcDedxParam::HV1_curvep3;
00077 par[4]= MdcDedxParam::HV1_curvep4;
00078 sigma_par[0] = MdcDedxParam::HV1_sigmap0;
00079 sigma_par[1] = MdcDedxParam::HV1_sigmap1;
00080 sigma_par[2] = MdcDedxParam::HV1_sigmap2;
00081 sigma_par[3] = MdcDedxParam::HV1_sigmap3;
00082 sigma_index_nhit = MdcDedxParam::HV1_index_nhit;
00083 sigma_index_sin = MdcDedxParam::HV1_index_sin;
00084 }
00085 else if(runflag==2) {
00086 par[0]= MdcDedxParam::HV2_curvep0;
00087 par[1]= MdcDedxParam::HV2_curvep1;
00088 par[2]= MdcDedxParam::HV2_curvep2;
00089 par[3]= MdcDedxParam::HV2_curvep3;
00090 par[4]= MdcDedxParam::HV2_curvep4;
00091 sigma_par[0] = MdcDedxParam::HV2_sigmap0;
00092 sigma_par[1] = MdcDedxParam::HV2_sigmap1;
00093 sigma_par[2] = MdcDedxParam::HV2_sigmap2;
00094 sigma_par[3] = MdcDedxParam::HV2_sigmap3;
00095 sigma_index_nhit = MdcDedxParam::HV2_index_nhit;
00096 sigma_index_sin = MdcDedxParam::HV2_index_sin;
00097 }
00098
00099
00100
00101 float prob( 0.0 );
00102 double delta;
00103 double sig_param;
00104 // for debag
00105 int Nmax_prob( 0 );
00106 float max_prob( -0.01 );
00107 int ndf;
00108 float prob_d[5],chi2;
00109
00110
00111 for( int it = 0; it < par_cand; it++ ) {
00112 beta_G = mom/Charge_Mass[it];
00113
00114 // ----- beta square ------
00115 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
00116 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
00117 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
00118
00119
00120 if( Nohit >0 ) {
00121 dedx_exp[it] = bethe_B;
00122 double sig_the=std::sin( (double)theta );
00123 //cout<<"nhit = "<<Nohit<<endl;
00124 double f_betagamma, g_sinth, h_nhit, i_t0;
00125
00126 if(runflag <3 && runflag>0){
00127 if(landau == 0) {
00128 float pathL_f = 1.2;
00129 double phi_f = 2*M_PI/12.0; // 30 deg.
00130 sig_param = 1.6*std::sin( (double) theta )/(lsamp*double(Nohit));
00131 ex_sigma[it] = 0.05*bethe_B*sqrt( 50.0*sig_param );
00132 }
00133 else {
00134 //currently use one sigmap0
00135 if(beta_G < 4) {
00136 sig_param=sigma_par[1]+sigma_par[2]*std::pow(beta_G,sigma_par[3]);
00137 } else {
00138 sig_param= sigma_par[0];
00139 }
00140 //double sig_the=std::sin( (double)theta );
00141 sig_the=std::pow(sig_the,sigma_index_sin);
00142 double sig_n;
00143 sig_n=35.0/double(Nohit);
00144 sig_n=std::pow(sig_n,sigma_index_nhit);
00145 ex_sigma[it]=sig_param*sig_the*sig_n;
00146 }
00147 }
00148
00149 MdcDedxTrk dedxtrk;
00150 double dedx_correc;
00151 //dedxtrk.clears();
00152 //chi_dedx[it] = (dedx - dedx_exp[it])/ex_sigma[it];
00153 //cout<<"dedx = "<<dedx<< " dedx_correc = "<<dedx_correc<<endl;
00154 if(runflag == 2 ) dedx_correc = dedxtrk.SpaceChargeCorrec(theta, mom, it+1, dedx);
00155 else if(runflag != 2) dedx_correc = dedx;
00156 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
00157 chi2 = chi_dedx[it]*chi_dedx[it];
00158 ndf = 1;
00159 pid_prob[it] = prob_(&chi2,&ndf);
00160 //if(it ==0 ) cout<<"runflag: "<<runflag<<" dedx : "<<dedx<<" chi_dedx: "
00161 //<<chi_dedx[it] <<" ptrk: "<<mom<<endl;
00162 if( it == -999 ){ // here a debug flag
00163 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
00164 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
00165 << std::endl;
00166 }
00167 if( pid_prob[it] > max_prob ) {
00168 max_prob = pid_prob[it];
00169 Nmax_prob = it;
00170 }
00171 }
00172 else {
00173 dedx_exp[it] = 0.0;
00174 ex_sigma[it] = 1000.0;
00175 pid_prob[it] = 0.0;
00176 chi_dedx[it] = 999.0;
00177 }
00178 }
00179 // std::cout<<"MdcDedxCorrection::dedx_pid_exp(blum)!!!"<<std::endl;
00180 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
00037 {
00038 if(!m_excalib) {
00039 m_excalib = new MdcDedxCorrection;
00040 }
00041 #ifdef DEBUG
00042 std::cout<<"MdcDedxCorrection::getMdcDedxCorrection()!!!"<<std::endl;
00043 #endif
00044 return m_excalib;
00045 }
|
|
|
|
|
|
MdcDedxCorrection class |
1.3.9.1