00001 #include <fstream>
00002 #include <cmath>
00003 #include <cstdlib>
00004
00005 #include "ParticleID/DedxPID.h"
00006 #ifndef BEAN
00007 #include "EvtRecEvent/EvtRecTrack.h"
00008 #include "MdcRecEvent/RecMdcTrack.h"
00009 #include "MdcRecEvent/RecMdcDedx.h"
00010 #else
00011 #include "RootEventDataVer.h"
00012 #include "DstEvtRecTracks.h"
00013 #endif
00014
00015 DedxPID * DedxPID::m_pointer = 0;
00016
00017 DedxPID * DedxPID::instance() {
00018 if(!m_pointer) m_pointer = new DedxPID();
00019 return m_pointer;
00020 }
00021
00022 DedxPID::DedxPID():ParticleIDBase() {
00023 m_readstate=0;
00024 }
00025
00026 void DedxPID::init() {
00027 for(int i = 0; i < 5; i++) {
00028 m_chi[i] = 99.0;
00029 m_prob[i] = -1.0;
00030 m_offset[i] = 99.0;
00031 m_sigma[i] = 1.0;
00032 }
00033 m_chimin = 99.;
00034 m_pdfmin =99;
00035 m_ndof = 0;
00036 m_goodHits = -99;
00037 m_normPH = -99;
00038 m_probPH = -99;
00039 m_nhitcutdx=5;
00040 }
00041
00042 void DedxPID::calculate() {
00043
00044 if(!m_readstate) {
00045 inputpar();
00046 m_readstate=1;
00047 }
00048 if(particleIDCalculation() == 0) m_ndof=1;
00049 }
00050
00051 int DedxPID::particleIDCalculation() {
00052
00053 int nhitcutdedx=getNhitCutDx();
00054 int irc = -1;
00055 EvtRecTrack* recTrk = PidTrk();
00056 if(!(recTrk->isMdcTrackValid())) return irc;
00057 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00058
00059 double ptrk = mdcTrk->p();
00060 int charge = mdcTrk->charge();
00061 if(ptrk>5) return irc;
00062 double cost = cos(mdcTrk->theta());
00063
00064
00065 if(!(recTrk->isMdcDedxValid())) return irc;
00066 RecMdcDedx* dedxTrk = recTrk->mdcDedx();
00067
00068 if((dedxTrk->normPH()>30)||(dedxTrk->normPH()<0)) return irc;
00069 m_goodHits = dedxTrk->numGoodHits();
00070 if(dedxTrk->numGoodHits()<nhitcutdedx) return irc;
00071 m_normPH = dedxTrk->normPH();
00072 m_probPH = dedxTrk->probPH();
00073
00074 double chitemp = 99.;
00075 double pdftemp = 0;
00076
00077
00078
00079 for(int i = 0; i < 5; i++) {
00080 double sep = dedxTrk->chi(i);
00081
00082 #ifndef BEAN
00083 string sftver = getenv("BES_RELEASE");
00084 string sft;
00085 sft.assign(sftver,0,5);
00086 if(sft=="6.6.0"||sft=="6.5.5") {
00087 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
00088 }
00089 else
00090 m_chi[i]=sep;
00091 #else
00092
00093 #if (ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,5,5) ||\
00094 ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,6,0) )
00095 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
00096 #else
00097 m_chi[i]=sep;
00098 #endif
00099 #endif
00100
00101 m_offset[i] = offsetDedx(i, ptrk, cost);
00102 m_sigma[i] = sigmaDedx(i, ptrk, cost);
00103 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
00104 double ppp = pdfCalculate(m_chi[i],1);
00105 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
00106
00107 }
00108 m_chimin = chitemp;
00109 m_pdfmin = pdftemp;
00110 if(m_chimin > chiMinCut()) return irc;
00111 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
00112
00113
00114
00115
00116 for(int i = 0; i < 5; i++)
00117 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
00118
00119 m_ndof = 1;
00120 irc = 0;
00121 return irc;
00122 }
00123
00124
00125
00126
00127
00128
00129
00130 double DedxPID::offsetDedx(int n, double ptrk, double cost) {
00131 return 0;
00132 }
00133
00134 double DedxPID::CorrDedx(int n, double ptrk, double cost,double chi,int charge) {
00135 int rundedx2 = getRunNo();
00136 double offset = 0.0;
00137 double offsetp = 0.0;
00138 double offsetc = 0.0;
00139 double sigcos = 1;
00140 double sigp = 1;
00141 double chicor=chi;
00142
00143
00144 switch(n) {
00145 case 0: {
00146 break;
00147 }
00148
00149 case 1: {
00150 break;
00151 }
00152
00153 case 2: {
00154
00155 double costm = cost;
00156 if(ptrk<0.1||ptrk>1) break;
00157 int index = int((ptrk-0.1)/0.05);
00158 if(index<=0) index=1;
00159 if(index>=17) index=16;
00160
00161 if(fabs(costm)>=0.8) break;
00162 int index1 = int((costm+0.8)/0.1);
00163 if(index1<=0) index1=1;
00164 if(index1>=15) index1=14;
00165
00166
00167 if(rundedx2>=11414&&rundedx2<=14604) {
00168 offsetp = cal_par(index,m_psipp_pi_ptrk_offset,ptrk,0.125,0.05);
00169 sigp = cal_par(index,m_psipp_pi_ptrk_sigma,ptrk,0.125,0.05);
00170 offsetc = cal_par(index1,m_psipp_pi_theta_offset,costm,-0.75,0.1);
00171 sigcos = cal_par(index1,m_psipp_pi_theta_sigma,costm,-0.75,0.1);
00172 }
00173
00174 if(rundedx2<=-11414&&rundedx2>=-14604) {
00175 offsetp = cal_par(index,m_psipp_mc_pi_ptrk_offset,ptrk,0.125,0.05);
00176 sigp = cal_par(index,m_psipp_mc_pi_ptrk_sigma,ptrk,0.125,0.05);
00177 offsetc = cal_par(index1,m_psipp_mc_pi_theta_offset,costm,-0.75,0.1);
00178 sigcos = cal_par(index1,m_psipp_mc_pi_theta_sigma,costm,-0.75,0.1);
00179 }
00180
00181 offset=offsetp+sigp*offsetc;
00182 chicor=(chicor-offset)/(sigcos*sigp);
00183 break;
00184 }
00185
00186 case 3: {
00187
00188 double costm = cost;
00189 if(ptrk<0.3||ptrk>0.8) break;
00190 offset=0;
00191 int index = int((ptrk-0.3)/0.1);
00192 if(index<=0) index=1;
00193 if(index>=4) index=3;
00194
00195 int index1 = int((costm+0.9)/0.1);
00196 if(index1<=0) index1=1;
00197 if(index1>=17) index1=16;
00198
00199 if(rundedx2>=9947&&rundedx2<=10878) {
00200 if(charge>0) {
00201 offsetp = cal_par(index,m_jpsi_kap_ptrk_offset,ptrk,0.35,0.1);
00202 sigp = cal_par(index,m_jpsi_kap_ptrk_sigma,ptrk,0.35,0.1);
00203 if(fabs(costm)<=0.83) {
00204 offsetc = cal_par(index1,m_jpsi_kap_theta_offset,costm,-0.85,0.1);
00205 sigcos = cal_par(index1,m_jpsi_kap_theta_sigma,costm,-0.85,0.1);
00206 }
00207 }
00208 if(charge<0) {
00209 offsetp = cal_par(index,m_jpsi_kam_ptrk_offset,ptrk,0.35,0.1);
00210 sigp = cal_par(index,m_jpsi_kam_ptrk_sigma,ptrk,0.35,0.1);
00211 if(fabs(costm)<=0.83) {
00212 offsetc = cal_par(index1,m_jpsi_kam_theta_offset,costm,-0.85,0.1);
00213 sigcos = cal_par(index1,m_jpsi_kam_theta_sigma,costm,-0.85,0.1);
00214 }
00215 }
00216 }
00217
00218
00219 if(rundedx2<=-9947&&rundedx2>=-10878) {
00220 if(charge>0) {
00221 offsetp = cal_par(index,m_jpsi_mc_kap_ptrk_offset,ptrk,0.35,0.1);
00222 sigp = cal_par(index,m_jpsi_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
00223 if(fabs(costm)<=0.83) {
00224 offsetc = cal_par(index1,m_jpsi_mc_kap_theta_offset,costm,-0.85,0.1);
00225 sigcos = cal_par(index1,m_jpsi_mc_kap_theta_sigma,costm,-0.85,0.1);
00226 }
00227 }
00228 if(charge<0) {
00229 offsetp = cal_par(index,m_jpsi_mc_kam_ptrk_offset,ptrk,0.35,0.1);
00230 sigp = cal_par(index,m_jpsi_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
00231 if(fabs(costm)<=0.83) {
00232 offsetc = cal_par(index1,m_jpsi_mc_kam_theta_offset,costm,-0.85,0.1);
00233 sigcos = cal_par(index1,m_jpsi_mc_kam_theta_sigma,costm,-0.85,0.1);
00234 }
00235 }
00236 }
00237
00238
00239 if(rundedx2>=8093&&rundedx2<=9025) {
00240 if(ptrk<0.3||ptrk>1.2) break;
00241 index = int((ptrk-0.3)/0.1);
00242 if(index<=0) index=1;
00243 if(index>=8) index=7;
00244 if(charge>0) {
00245 offsetp = cal_par(index,m_psip_kap_ptrk_offset,ptrk,0.35,0.1);
00246 sigp = cal_par(index,m_psip_kap_ptrk_sigma,ptrk,0.35,0.1);
00247 }
00248 if(charge<0) {
00249 offsetp = cal_par(index,m_psip_kam_ptrk_offset,ptrk,0.35,0.1);
00250 sigp = cal_par(index,m_psip_kam_ptrk_sigma,ptrk,0.35,0.1);
00251 }
00252 }
00253
00254
00255 if(rundedx2<=-8093&&rundedx2>=-9025) {
00256
00257 if(ptrk<0.3||ptrk>1.2) break;
00258 index = int((ptrk-0.3)/0.1);
00259 if(index<=0) index=1;
00260 if(index>=8) index=7;
00261 if(charge>0) {
00262 offsetp = cal_par(index,m_psip_mc_kap_ptrk_offset,ptrk,0.35,0.1);
00263 sigp = cal_par(index,m_psip_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
00264 }
00265 if(charge<0) {
00266 offsetp = cal_par(index,m_psip_mc_kam_ptrk_offset,ptrk,0.35,0.1);
00267 sigp = cal_par(index,m_psip_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
00268 }
00269 }
00270
00271
00272
00273 if(rundedx2>=11414&&rundedx2<=14604) {
00274 if(ptrk<0.15||ptrk>1) break;
00275 index = int((ptrk-0.15)/0.05);
00276 if(index<=0) index=1;
00277 if(index>=16) index=15;
00278 if(fabs(costm)>=0.8) break;
00279 index1 = int((costm+0.8)/0.1);
00280 if(index1<=0) index1=1;
00281 if(index1>=15) index1=14;
00282
00283 offsetp = cal_par(index,m_psipp_ka_ptrk_offset,ptrk,0.175,0.05);
00284 sigp = cal_par(index,m_psipp_ka_ptrk_sigma,ptrk,0.175,0.05);
00285 offsetc = cal_par(index1,m_psipp_ka_theta_offset,costm,-0.75,0.1);
00286 sigcos = cal_par(index1,m_psipp_ka_theta_sigma,costm,-0.75,0.1);
00287 }
00288
00289 if(rundedx2<=-11414&&rundedx2>=-14604) {
00290 if(ptrk<0.15||ptrk>1) break;
00291 index = int((ptrk-0.15)/0.05);
00292 if(index<=0) index=1;
00293 if(index>=16) index=15;
00294 if(fabs(costm)>=0.8) break;
00295 index1 = int((costm+0.8)/0.1);
00296 if(index1<=0) index1=1;
00297 if(index1>=15) index1=14;
00298 offsetp = cal_par(index,m_psipp_mc_ka_ptrk_offset,ptrk,0.175,0.05);
00299 sigp = cal_par(index,m_psipp_mc_ka_ptrk_sigma,ptrk,0.175,0.05);
00300 offsetc = cal_par(index1,m_psipp_mc_ka_theta_offset,costm,-0.75,0.1);
00301 sigcos = cal_par(index1,m_psipp_mc_ka_theta_sigma,costm,-0.75,0.1);
00302 }
00303
00304 offset=offsetp+sigp*offsetc;
00305 chicor=(chicor-offset)/(sigcos*sigp);
00306 break;
00307 }
00308
00309 case 4 : {
00310
00311 double costm = cost;
00312 if(ptrk<0.3||ptrk>1.1) break;
00313 int index = int((ptrk-0.3)/0.1);
00314 if(index<=0) index=1;
00315 if(index>=7) index=6;
00316
00317 int index1 = int((costm+0.9)/0.1);
00318 if(index1<=0) index1=1;
00319 if(index1>=17) index1=16;
00320
00321
00322 offset=0;
00323 if(rundedx2>=9947&&rundedx2<=10878) {
00324 if(charge>0) {
00325 offsetp = cal_par(index,m_jpsi_protonp_ptrk_offset,ptrk,0.35,0.1);
00326 sigp = cal_par(index,m_jpsi_protonp_ptrk_sigma,ptrk,0.35,0.1);
00327 if(fabs(costm)<=0.83) {
00328 offsetc = cal_par(index1,m_jpsi_protonp_theta_offset,costm,-0.85,0.1);
00329 sigcos = cal_par(index1,m_jpsi_protonp_theta_sigma,costm,-0.85,0.1);
00330 }
00331 }
00332 if(charge<0) {
00333 offsetp = cal_par(index,m_jpsi_protonm_ptrk_offset,ptrk,0.35,0.1);
00334 sigp = cal_par(index,m_jpsi_protonm_ptrk_sigma,ptrk,0.35,0.1);
00335 if(fabs(costm)<=0.83) {
00336 offsetc = cal_par(index1,m_jpsi_protonm_theta_offset,costm,-0.85,0.1);
00337 sigcos = cal_par(index1,m_jpsi_protonm_theta_sigma,costm,-0.85,0.1);
00338 }
00339 }
00340 }
00341
00342
00343 if(rundedx2<=-9947&&rundedx2>=-10878) {
00344 if(charge>0) {
00345 offsetp = cal_par(index,m_jpsi_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
00346 sigp = cal_par(index,m_jpsi_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
00347 if(fabs(costm)<=0.83) {
00348 offsetc = cal_par(index1,m_jpsi_mc_protonp_theta_offset,costm,-0.85,0.1);
00349 sigcos = cal_par(index1,m_jpsi_mc_protonp_theta_sigma,costm,-0.85,0.1);
00350 }
00351 }
00352 if(charge<0) {
00353 offsetp = cal_par(index,m_jpsi_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
00354 sigp = cal_par(index,m_jpsi_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
00355 if(fabs(costm)<=0.83) {
00356 offsetc = cal_par(index1,m_jpsi_mc_protonm_theta_offset,costm,-0.85,0.1);
00357 sigcos = cal_par(index1,m_jpsi_mc_protonm_theta_sigma,costm,-0.85,0.1);
00358 }
00359 }
00360 }
00361
00362
00363 if(rundedx2>=8093&&rundedx2<=9025) {
00364 if(charge>0) {
00365 offsetp = cal_par(index,m_psip_protonp_ptrk_offset,ptrk,0.35,0.1);
00366 sigp = cal_par(index,m_psip_protonp_ptrk_sigma,ptrk,0.35,0.1);
00367 }
00368 if(charge<0) {
00369 offsetp = cal_par(index,m_psip_protonm_ptrk_offset,ptrk,0.35,0.1);
00370 sigp = cal_par(index,m_psip_protonm_ptrk_sigma,ptrk,0.35,0.1);
00371 }
00372 }
00373
00374
00375 if(rundedx2<=-8093&&rundedx2>=-9025) {
00376 if(charge>0) {
00377 offsetp = cal_par(index,m_psip_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
00378 sigp = cal_par(index,m_psip_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
00379 }
00380 if(charge<0) {
00381 offsetp = cal_par(index,m_psip_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
00382 sigp = cal_par(index,m_psip_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
00383 }
00384 }
00385
00386
00387 if(rundedx2>=11414&&rundedx2<=14604) {
00388 if(ptrk<0.2||ptrk>1.1) break;
00389 index = int((ptrk-0.2)/0.05);
00390 if(index<=0) index=1;
00391 if(index>=17) index=16;
00392 if(fabs(costm)>=0.83) break;
00393 index1 = int((costm+0.9)/0.1);
00394 if(index1<=0) index1=1;
00395 if(index1>=17) index1=16;
00396
00397 offsetp = cal_par(index,m_psipp_proton_ptrk_offset,ptrk,0.225,0.05);
00398 sigp = cal_par(index,m_psipp_proton_ptrk_sigma,ptrk,0.225,0.05);
00399 offsetc = cal_par(index1,m_psipp_proton_theta_offset,costm,-0.85,0.1);
00400 sigcos = cal_par(index1,m_psipp_proton_theta_sigma,costm,-0.85,0.1);
00401 }
00402
00403 if(rundedx2<=-11414&&rundedx2>=-14604) {
00404 if(ptrk<0.2||ptrk>1.1) break;
00405 index = int((ptrk-0.2)/0.1);
00406 if(index<=0) index=1;
00407 if(index>=8) index=7;
00408 if(fabs(costm)>=0.83) break;
00409 index1 = int((costm+0.9)/0.1);
00410 if(index1<=0) index1=1;
00411 if(index1>=17) index1=16;
00412 offsetp = cal_par(index,m_psipp_mc_proton_ptrk_offset,ptrk,0.25,0.1);
00413 sigp = cal_par(index,m_psipp_mc_proton_ptrk_sigma,ptrk,0.25,0.1);
00414 offsetc = cal_par(index1,m_psipp_mc_proton_theta_offset,costm,-0.85,0.1);
00415 sigcos = cal_par(index1,m_psipp_mc_proton_theta_sigma,costm,-0.85,0.1);
00416 }
00417 offset=offsetp+sigp*offsetc;
00418 chicor=(chicor-offset)/(sigcos*sigp);
00419 break;
00420 }
00421
00422 default:
00423 offset = 0.0;
00424 break;
00425 }
00426
00427 return chicor;
00428 }
00429
00430 double DedxPID::sigmaDedx(int n, double ptrk, double cost) {
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 return 1;
00478
00479 }
00480
00481 double DedxPID::mypol3(double x, double par0, double par1, double par2, double par3)
00482 {
00483 double y = x;
00484 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y);
00485
00486 }
00487
00488 double DedxPID::mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
00489 {
00490 double y = x;
00491 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y) + (par4 * y * y * y *y)+ (par5 * y * y * y * y * y);
00492
00493 }
00494
00495 void DedxPID::inputpar() {
00496
00497
00498 std::string jpsi_kap_mom = path + "/share/JPsi/kaon/dedx_kap.txt";
00499 std::string jpsi_kap_mom_mc = path + "/share/JPsi/kaon/dedx_kap_mc.txt";
00500 ifstream inputmomdata6(jpsi_kap_mom.c_str(),std::ios_base::in);
00501 if ( !inputmomdata6 ) {
00502 cout << " can not open: " << jpsi_kap_mom << endl;
00503 exit(1);
00504 }
00505 ifstream inputmomdata6mc(jpsi_kap_mom_mc.c_str(),std::ios_base::in);
00506 if ( !inputmomdata6mc ) {
00507 cout << " can not open: " << jpsi_kap_mom_mc << endl;
00508 exit(1);
00509 }
00510 for(int i=0; i<12; i++) {
00511 inputmomdata6>>m_jpsi_kap_ptrk_offset[i];
00512 inputmomdata6>>m_jpsi_kap_ptrk_sigma[i];
00513 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_offset[i];
00514 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_sigma[i];
00515 }
00516
00517
00518 std::string jpsi_kam_mom = path + "/share/JPsi/kaon/dedx_kam.txt";
00519 std::string jpsi_kam_mom_mc = path + "/share/JPsi/kaon/dedx_kam_mc.txt";
00520 ifstream inputmomdata7(jpsi_kam_mom.c_str(),std::ios_base::in);
00521 if ( !inputmomdata7 ) {
00522 cout << " can not open: " << jpsi_kam_mom << endl;
00523 exit(1);
00524 }
00525 ifstream inputmomdata7mc(jpsi_kam_mom_mc.c_str(),std::ios_base::in);
00526 if ( !inputmomdata7mc ) {
00527 cout << " can not open: " << jpsi_kam_mom_mc << endl;
00528 exit(1);
00529 }
00530 for(int i=0; i<12; i++) {
00531 inputmomdata7>>m_jpsi_kam_ptrk_offset[i];
00532 inputmomdata7>>m_jpsi_kam_ptrk_sigma[i];
00533 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_offset[i];
00534 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_sigma[i];
00535
00536 }
00537
00538
00539
00540 std::string jpsi_kap_the = path + "/share/JPsi/kaon/dedx_kap_theta.txt";
00541 std::string jpsi_kap_the_mc = path + "/share/JPsi/kaon/dedx_kap_theta_mc.txt";
00542 ifstream inputmomdata8(jpsi_kap_the.c_str(),std::ios_base::in);
00543 if ( !inputmomdata8 ) {
00544 cout << " can not open: " << jpsi_kap_the << endl;
00545 exit(1);
00546 }
00547 ifstream inputmomdata8mc(jpsi_kap_the_mc.c_str(),std::ios_base::in);
00548 if ( !inputmomdata8mc ) {
00549 cout << " can not open: " << jpsi_kap_the_mc << endl;
00550 exit(1);
00551 }
00552 for(int i=0; i<18; i++) {
00553 inputmomdata8>>m_jpsi_kap_theta_offset[i];
00554 inputmomdata8>>m_jpsi_kap_theta_sigma[i];
00555 inputmomdata8mc>>m_jpsi_mc_kap_theta_offset[i];
00556 inputmomdata8mc>>m_jpsi_mc_kap_theta_sigma[i];
00557 }
00558
00559
00560 std::string jpsi_kam_the = path + "/share/JPsi/kaon/dedx_kam_theta.txt";
00561 std::string jpsi_kam_the_mc = path + "/share/JPsi/kaon/dedx_kam_theta_mc.txt";
00562 ifstream inputmomdata9(jpsi_kam_the.c_str(),std::ios_base::in);
00563 if ( !inputmomdata9 ) {
00564 cout << " can not open: " << jpsi_kam_the << endl;
00565 exit(1);
00566 }
00567 ifstream inputmomdata9mc(jpsi_kam_the_mc.c_str(),std::ios_base::in);
00568 if ( !inputmomdata9mc ) {
00569 cout << " can not open: " << jpsi_kam_the_mc << endl;
00570 exit(1);
00571 }
00572 for(int i=0; i<18; i++) {
00573 inputmomdata9>>m_jpsi_kam_theta_offset[i];
00574 inputmomdata9>>m_jpsi_kam_theta_sigma[i];
00575 inputmomdata9mc>>m_jpsi_mc_kam_theta_offset[i];
00576 inputmomdata9mc>>m_jpsi_mc_kam_theta_sigma[i];
00577 }
00578
00579
00580 std::string jpsi_protonp_mom = path + "/share/JPsi/proton/dedx_protonp.txt";
00581 std::string jpsi_protonp_mom_mc = path + "/share/JPsi/proton/dedx_protonp_mc.txt";
00582 ifstream inputmomdata12(jpsi_protonp_mom.c_str(),std::ios_base::in);
00583 if ( !inputmomdata12 ) {
00584 cout << " can not open: " << jpsi_protonp_mom << endl;
00585 exit(1);
00586 }
00587 ifstream inputmomdata12mc(jpsi_protonp_mom_mc.c_str(),std::ios_base::in);
00588 if ( !inputmomdata12mc ) {
00589 cout << " can not open: " << jpsi_protonp_mom_mc << endl;
00590 exit(1);
00591 }
00592 for(int i=0; i<8; i++) {
00593 inputmomdata12>>m_jpsi_protonp_ptrk_offset[i];
00594 inputmomdata12>>m_jpsi_protonp_ptrk_sigma[i];
00595 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_offset[i];
00596 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_sigma[i];
00597 }
00598
00599
00600 std::string jpsi_protonm_mom = path + "/share/JPsi/proton/dedx_protonm.txt";
00601 std::string jpsi_protonm_mom_mc = path + "/share/JPsi/proton/dedx_protonm_mc.txt";
00602 ifstream inputmomdata13(jpsi_protonm_mom.c_str(),std::ios_base::in);
00603 if ( !inputmomdata13 ) {
00604 cout << " can not open: " << jpsi_protonm_mom << endl;
00605 exit(1);
00606 }
00607 ifstream inputmomdata13mc(jpsi_protonm_mom_mc.c_str(),std::ios_base::in);
00608 if ( !inputmomdata13mc ) {
00609 cout << " can not open: " << jpsi_protonm_mom_mc << endl;
00610 exit(1);
00611 }
00612 for(int i=0; i<8; i++) {
00613 inputmomdata13>>m_jpsi_protonm_ptrk_offset[i];
00614 inputmomdata13>>m_jpsi_protonm_ptrk_sigma[i];
00615 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_offset[i];
00616 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_sigma[i];
00617 }
00618
00619
00620 std::string jpsi_protonp_the = path + "/share/JPsi/proton/dedx_protonp_theta.txt";
00621 std::string jpsi_protonp_the_mc = path + "/share/JPsi/proton/dedx_protonp_theta_mc.txt";
00622
00623 ifstream inputmomdata14(jpsi_protonp_the.c_str(),std::ios_base::in);
00624 if ( !inputmomdata14 ) {
00625 cout << " can not open: " << jpsi_protonp_the << endl;
00626 exit(1);
00627 }
00628 ifstream inputmomdata14mc(jpsi_protonp_the_mc.c_str(),std::ios_base::in);
00629 if ( !inputmomdata14mc ) {
00630 cout << " can not open: " << jpsi_protonp_the_mc << endl;
00631 exit(1);
00632 }
00633 for(int i=0; i<18; i++) {
00634 inputmomdata14>>m_jpsi_protonp_theta_offset[i];
00635 inputmomdata14>>m_jpsi_protonp_theta_sigma[i];
00636 inputmomdata14mc>>m_jpsi_mc_protonp_theta_offset[i];
00637 inputmomdata14mc>>m_jpsi_mc_protonp_theta_sigma[i];
00638 }
00639
00640
00641 std::string jpsi_protonm_the = path + "/share/JPsi/proton/dedx_protonm_theta.txt";
00642 std::string jpsi_protonm_the_mc = path + "/share/JPsi/proton/dedx_protonm_theta_mc.txt";
00643 ifstream inputmomdata15(jpsi_protonm_the.c_str(),std::ios_base::in);
00644 if ( !inputmomdata15 ) {
00645 cout << " can not open: " << jpsi_protonm_the << endl;
00646 exit(1);
00647 }
00648 ifstream inputmomdata15mc(jpsi_protonm_the_mc.c_str(),std::ios_base::in);
00649 if ( !inputmomdata15mc ) {
00650 cout << " can not open: " << jpsi_protonm_the_mc << endl;
00651 exit(1);
00652 }
00653 for(int i=0; i<18; i++) {
00654 inputmomdata15>>m_jpsi_protonm_theta_offset[i];
00655 inputmomdata15>>m_jpsi_protonm_theta_sigma[i];
00656 inputmomdata15mc>>m_jpsi_mc_protonm_theta_offset[i];
00657 inputmomdata15mc>>m_jpsi_mc_protonm_theta_sigma[i];
00658 }
00659
00660
00661
00662
00663
00664 std::string psip_kap_mom = path + "/share/Psip/kaon/dedx_kap.txt";
00665 std::string psip_kap_mom_mc = path + "/share/Psip/kaon/dedx_kap_mc.txt";
00666 ifstream inputmomdata24(psip_kap_mom.c_str(),std::ios_base::in);
00667 if ( !inputmomdata24 ) {
00668 cout << " can not open: " << psip_kap_mom << endl;
00669 exit(1);
00670 }
00671 ifstream inputmomdata24mc(psip_kap_mom_mc.c_str(),std::ios_base::in);
00672 if ( !inputmomdata24mc ) {
00673 cout << " can not open: " << psip_kap_mom_mc << endl;
00674 exit(1);
00675 }
00676 for(int i=0; i<9; i++) {
00677 inputmomdata24>>m_psip_kap_ptrk_offset[i];
00678 inputmomdata24>>m_psip_kap_ptrk_sigma[i];
00679 inputmomdata24mc>>m_psip_mc_kap_ptrk_offset[i];
00680 inputmomdata24mc>>m_psip_mc_kap_ptrk_sigma[i];
00681 }
00682
00683
00684 std::string psip_kam_mom = path + "/share/Psip/kaon/dedx_kam.txt";
00685 std::string psip_kam_mom_mc = path + "/share/Psip/kaon/dedx_kam_mc.txt";
00686 ifstream inputmomdata25(psip_kam_mom.c_str(),std::ios_base::in);
00687 if ( !inputmomdata25 ) {
00688 cout << " can not open: " << psip_kam_mom << endl;
00689 exit(1);
00690 }
00691 ifstream inputmomdata25mc(psip_kam_mom_mc.c_str(),std::ios_base::in);
00692 if ( !inputmomdata25mc ) {
00693 cout << " can not open: " << psip_kam_mom_mc << endl;
00694 exit(1);
00695 }
00696 for(int i=0; i<9; i++) {
00697 inputmomdata25>>m_psip_kam_ptrk_offset[i];
00698 inputmomdata25>>m_psip_kam_ptrk_sigma[i];
00699 inputmomdata25mc>>m_psip_mc_kam_ptrk_offset[i];
00700 inputmomdata25mc>>m_psip_mc_kam_ptrk_sigma[i];
00701 }
00702
00703
00704
00705 std::string psip_protonp_mom = path + "/share/Psip/proton/dedx_protonp.txt";
00706 std::string psip_protonp_mom_mc = path + "/share/Psip/proton/dedx_protonp_mc.txt";
00707 ifstream inputmomdata26(psip_protonp_mom.c_str(),std::ios_base::in);
00708 if ( !inputmomdata26 ) {
00709 cout << " can not open: " << psip_protonp_mom << endl;
00710 exit(1);
00711 }
00712 ifstream inputmomdata26mc(psip_protonp_mom_mc.c_str(),std::ios_base::in);
00713 if ( !inputmomdata26mc ) {
00714 cout << " can not open: " << psip_protonp_mom_mc << endl;
00715 exit(1);
00716 }
00717 for(int i=0; i<9; i++) {
00718 inputmomdata26>>m_psip_protonp_ptrk_offset[i];
00719 inputmomdata26>>m_psip_protonp_ptrk_sigma[i];
00720 inputmomdata26mc>>m_psip_mc_protonp_ptrk_offset[i];
00721 inputmomdata26mc>>m_psip_mc_protonp_ptrk_sigma[i];
00722 }
00723
00724
00725 std::string psip_protonm_mom = path + "/share/Psip/proton/dedx_protonm.txt";
00726 std::string psip_protonm_mom_mc = path + "/share/Psip/proton/dedx_protonm_mc.txt";
00727 ifstream inputmomdata27(psip_protonm_mom.c_str(),std::ios_base::in);
00728 if ( !inputmomdata27 ) {
00729 cout << " can not open: " << psip_protonm_mom << endl;
00730 exit(1);
00731 }
00732 ifstream inputmomdata27mc(psip_protonm_mom_mc.c_str(),std::ios_base::in);
00733 if ( !inputmomdata27mc ) {
00734 cout << " can not open: " << psip_protonm_mom_mc << endl;
00735 exit(1);
00736 }
00737 for(int i=0; i<9; i++) {
00738 inputmomdata27>>m_psip_protonm_ptrk_offset[i];
00739 inputmomdata27>>m_psip_protonm_ptrk_sigma[i];
00740 inputmomdata27mc>>m_psip_mc_protonm_ptrk_offset[i];
00741 inputmomdata27mc>>m_psip_mc_protonm_ptrk_sigma[i];
00742 }
00743
00744
00745 std::string psipp_pi_mom = path + "/share/Psipp/pion/dedx_pi.txt";
00746 std::string psipp_pi_mom_mc = path + "/share/Psipp/pion/dedx_pi_mc.txt";
00747 ifstream inputmomdata28(psipp_pi_mom.c_str(),std::ios_base::in);
00748 if ( !inputmomdata28 ) {
00749 cout << " can not open: " << psipp_pi_mom << endl;
00750 exit(1);
00751 }
00752 ifstream inputmomdata28mc(psipp_pi_mom_mc.c_str(),std::ios_base::in);
00753 if ( !inputmomdata28mc ) {
00754 cout << " can not open: " << psipp_pi_mom_mc << endl;
00755 exit(1);
00756 }
00757 for(int i=0; i<18; i++) {
00758 inputmomdata28>>m_psipp_pi_ptrk_offset[i];
00759 inputmomdata28>>m_psipp_pi_ptrk_sigma[i];
00760 inputmomdata28mc>>m_psipp_mc_pi_ptrk_offset[i];
00761 inputmomdata28mc>>m_psipp_mc_pi_ptrk_sigma[i];
00762 }
00763
00764
00765 std::string psipp_pi_the = path + "/share/Psipp/pion/dedx_pi_theta.txt";
00766 std::string psipp_pi_the_mc = path + "/share/Psipp/pion/dedx_pi_theta_mc.txt";
00767 ifstream inputmomdata29(psipp_pi_the.c_str(),std::ios_base::in);
00768 if ( !inputmomdata29 ) {
00769 cout << " can not open: " << psipp_pi_the << endl;
00770 exit(1);
00771 }
00772 ifstream inputmomdata29mc(psipp_pi_the_mc.c_str(),std::ios_base::in);
00773 if ( !inputmomdata29mc ) {
00774 cout << " can not open: " << psipp_pi_the_mc << endl;
00775 exit(1);
00776 }
00777 for(int i=0; i<16; i++) {
00778 inputmomdata29>>m_psipp_pi_theta_offset[i];
00779 inputmomdata29>>m_psipp_pi_theta_sigma[i];
00780 inputmomdata29mc>>m_psipp_mc_pi_theta_offset[i];
00781 inputmomdata29mc>>m_psipp_mc_pi_theta_sigma[i];
00782 }
00783
00784
00785 std::string psipp_ka_mom = path + "/share/Psipp/kaon/dedx_ka.txt";
00786 std::string psipp_ka_mom_mc = path + "/share/Psipp/kaon/dedx_ka_mc.txt";
00787 ifstream inputmomdata30(psipp_ka_mom.c_str(),std::ios_base::in);
00788 if ( !inputmomdata30 ) {
00789 cout << " can not open: " << psipp_ka_mom << endl;
00790 exit(1);
00791 }
00792 ifstream inputmomdata30mc(psipp_ka_mom_mc.c_str(),std::ios_base::in);
00793 if ( !inputmomdata30mc ) {
00794 cout << " can not open: " << psipp_ka_mom_mc << endl;
00795 exit(1);
00796 }
00797 for(int i=0; i<17; i++) {
00798 inputmomdata30>>m_psipp_ka_ptrk_offset[i];
00799 inputmomdata30>>m_psipp_ka_ptrk_sigma[i];
00800 inputmomdata30mc>>m_psipp_mc_ka_ptrk_offset[i];
00801 inputmomdata30mc>>m_psipp_mc_ka_ptrk_sigma[i];
00802 }
00803
00804
00805 std::string psipp_ka_the = path + "/share/Psipp/kaon/dedx_ka_theta.txt";
00806 std::string psipp_ka_the_mc = path + "/share/Psipp/kaon/dedx_ka_theta_mc.txt";
00807 ifstream inputmomdata31(psipp_ka_the.c_str(),std::ios_base::in);
00808 if ( !inputmomdata31 ) {
00809 cout << " can not open: " << psipp_ka_the << endl;
00810 exit(1);
00811 }
00812 ifstream inputmomdata31mc(psipp_ka_the_mc.c_str(),std::ios_base::in);
00813 if ( !inputmomdata31mc ) {
00814 cout << " can not open: " << psipp_ka_the_mc << endl;
00815 exit(1);
00816 }
00817 for(int i=0; i<16; i++) {
00818 inputmomdata31>>m_psipp_ka_theta_offset[i];
00819 inputmomdata31>>m_psipp_ka_theta_sigma[i];
00820 inputmomdata31mc>>m_psipp_mc_ka_theta_offset[i];
00821 inputmomdata31mc>>m_psipp_mc_ka_theta_sigma[i];
00822 }
00823
00824
00825
00826 std::string psipp_proton_mom = path + "/share/Psipp/proton/dedx_proton.txt";
00827 std::string psipp_proton_mom_mc = path + "/share/Psipp/proton/dedx_proton_mc.txt";
00828 ifstream inputmomdata32(psipp_proton_mom.c_str(),std::ios_base::in);
00829 if ( !inputmomdata32 ) {
00830 cout << " can not open: " << psipp_proton_mom << endl;
00831 exit(1);
00832 }
00833 ifstream inputmomdata32mc(psipp_proton_mom_mc.c_str(),std::ios_base::in);
00834 if ( !inputmomdata32mc ) {
00835 cout << " can not open: " << psipp_proton_mom_mc << endl;
00836 exit(1);
00837 }
00838 for(int i=0; i<18; i++) {
00839 inputmomdata32>>m_psipp_proton_ptrk_offset[i];
00840 inputmomdata32>>m_psipp_proton_ptrk_sigma[i];
00841 }
00842 for(int i=0; i<9; i++) {
00843 inputmomdata32mc>>m_psipp_mc_proton_ptrk_offset[i];
00844 inputmomdata32mc>>m_psipp_mc_proton_ptrk_sigma[i];
00845 }
00846
00847
00848 std::string psipp_proton_the = path + "/share/Psipp/proton/dedx_proton_theta.txt";
00849 std::string psipp_proton_the_mc = path + "/share/Psipp/proton/dedx_proton_theta_mc.txt";
00850 ifstream inputmomdata33(psipp_proton_the.c_str(),std::ios_base::in);
00851 if ( !inputmomdata33 ) {
00852 cout << " can not open: " << psipp_proton_the << endl;
00853 exit(1);
00854 }
00855 ifstream inputmomdata33mc(psipp_proton_the_mc.c_str(),std::ios_base::in);
00856 if ( !inputmomdata33mc ) {
00857 cout << " can not open: " << psipp_proton_the_mc << endl;
00858 exit(1);
00859 }
00860 for(int i=0; i<18; i++) {
00861 inputmomdata33>>m_psipp_proton_theta_offset[i];
00862 inputmomdata33>>m_psipp_proton_theta_sigma[i];
00863 inputmomdata33mc>>m_psipp_mc_proton_theta_offset[i];
00864 inputmomdata33mc>>m_psipp_mc_proton_theta_sigma[i];
00865 }
00866
00867 }
00868
00869 double DedxPID::iterate(double ptrk,double *mean,double *p) {
00870 double p1,p2,p3;
00871 p2=((mean[0]-mean[1])*(p[1]*p[1]-p[2]*p[2])-(mean[1]-mean[2])*(p[0]*p[0]-p[1]*p[1]))/((p[0]-p[1])*(p[1]*p[1]-p[2]*p[2])-(p[1]-p[2])*(p[0]*p[0]-p[1]*p[1]));
00872 p3=((p[0]-p[1])*(mean[1]-mean[2])-(p[1]-p[2])*(mean[0]-mean[1]))/((p[0]-p[1])*(p[1]*p[1]-p[2]*p[2])-(p[1]-p[2])*(p[0]*p[0]-p[1]*p[1]));
00873 p1=mean[0]-p2*p[0]-p3*p[0]*p[0];
00874 double mean1 = p1+p2*ptrk+p3*ptrk*ptrk;
00875 return mean1;
00876 }
00877
00878 double DedxPID::cal_par(int index1,double *m_jpsi_pip_ptrk_offset,double ptrk,double begin,double bin) {
00879 double mean1[3],p[3];
00880 p[0]=begin+(index1-1)*bin;
00881 p[1]=begin+index1*bin;
00882 p[2]=begin+(index1+1)*bin;
00883 mean1[0]=m_jpsi_pip_ptrk_offset[index1-1];
00884 mean1[1]=m_jpsi_pip_ptrk_offset[index1];
00885 mean1[2]=m_jpsi_pip_ptrk_offset[index1+1];
00886 double res=iterate(ptrk,mean1,p);
00887 return res;
00888 }
00889
