00001 #include <cmath>
00002 #include "TMath.h"
00003
00004 #include "ParticleID/TofCorrPID.h"
00005
00006 #ifndef BEAN
00007 #include "EvtRecEvent/EvtRecTrack.h"
00008 #include "MdcRecEvent/RecMdcTrack.h"
00009 #include "MdcRecEvent/RecMdcKalTrack.h"
00010 #include "ExtEvent/RecExtTrack.h"
00011 #include "TofRecEvent/RecTofTrack.h"
00012 #include "DstEvent/TofHitStatus.h"
00013 #else
00014 #include "TofHitStatus.h"
00015 #endif
00016 #include <fstream>
00017
00018
00019 TofCorrPID * TofCorrPID::m_pointer = 0;
00020 TofCorrPID * TofCorrPID::instance() {
00021 if(!m_pointer) m_pointer = new TofCorrPID();
00022 return m_pointer;
00023 }
00024
00025 TofCorrPID::TofCorrPID():ParticleIDBase() {
00026 m_runBegin = 0;
00027 m_runEnd = 0;
00028 }
00029
00030 void TofCorrPID::init() {
00031
00032 for(int i = 0; i < 5; i++) {
00033 m_chi[i] = -100.0;
00034 m_prob[i] = -1.0;
00035 m_sigma[i] = 10.0;
00036 m_offset[i] = -1000.0;
00037 }
00038 m_chimin = -99.0;
00039 m_chimax = 99.0;
00040 m_pdfmin = 99.;
00041 m_ndof = 0;
00042
00043 m_ipmt = -1;
00044 for( unsigned int i=0; i<5; i++ ) {
00045 for( unsigned int j=0; j<7; j++ ) {
00046 m_dt[i][j] = -1000.0;
00047 m_dtCorr[i][j] = -1000.0;
00048 m_sigCorr[i][j] = 10.0;
00049 m_chiCorr[i][j] = 100.0;
00050 }
00051 }
00052
00053 int run = getRunNo();
00054 if( !( abs(run)>=m_runBegin && abs(run)<=m_runEnd ) ) {
00055 inputParameter( getRunNo() );
00056 }
00057
00058 return;
00059 }
00060
00061
00062 void TofCorrPID::calculate() {
00063 if(particleIDCalculation() == 0) m_ndof=1;
00064 }
00065
00066
00067 int TofCorrPID::particleIDCalculation() {
00068 int irc=-1;
00069
00070 EvtRecTrack* recTrk = PidTrk();
00071 if(!(recTrk->isMdcTrackValid())) return irc;
00072 if(!(recTrk->isMdcKalTrackValid())) return irc;
00073 if(!(recTrk->isExtTrackValid())) return irc;
00074 if(!(recTrk->isTofTrackValid())) return irc;
00075
00076 #ifndef BEAN
00077 SmartRefVector<RecTofTrack> tofTrk = recTrk->tofTrack();
00078 SmartRefVector<RecTofTrack>::iterator it;
00079 #else
00080 const std::vector<TTofTrack* >& tofTrk = recTrk->tofTrack();
00081 std::vector<TTofTrack* >::const_iterator it;
00082 #endif
00083
00084 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00085 int charge = mdcTrk->charge();
00086
00087 double p[5], betaGamma[5];
00088 RecMdcKalTrack* kalTrk = recTrk->mdcKalTrack();
00089 for( int i=0; i<5; i++ ) {
00090 if( i==0 ) {
00091 kalTrk->setPidType(RecMdcKalTrack::electron);
00092 }
00093 else if( i==1 ) {
00094 kalTrk->setPidType(RecMdcKalTrack::muon);
00095 }
00096 else if( i==2 ) {
00097 kalTrk->setPidType(RecMdcKalTrack::pion);
00098 }
00099 else if( i==3 ) {
00100 kalTrk->setPidType(RecMdcKalTrack::kaon);
00101 }
00102 else if( i==4 ) {
00103 kalTrk->setPidType(RecMdcKalTrack::proton);
00104 }
00105 #ifndef BEAN
00106 HepLorentzVector ptrk = kalTrk->p4();
00107 #else
00108 HepLorentzVector ptrk = kalTrk->p4( xmass(i) );
00109 #endif
00110 p[i] = ptrk.rho();
00111 betaGamma[i] = p[i]/xmass(i);
00112 }
00113
00114 double zrhit[2];
00115 RecExtTrack* extTrk = recTrk->extTrack();
00116 zrhit[0] = extTrk->tof1Position().z();
00117 zrhit[1] = extTrk->tof2Position().z();
00118
00119 int tofid[2] = { -9, -9 };
00120
00121 m_ipmt = -1;
00122 bool readFile = false;
00123 bool signal[2] = { false, false };
00124 TofHitStatus *hitst = new TofHitStatus;
00125 for( it = tofTrk.begin(); it!=tofTrk.end(); it++ ) {
00126 unsigned int st = (*it)->status();
00127 hitst->setStatus(st);
00128 if( hitst->is_raw() ) return irc;
00129 bool barrel = hitst->is_barrel();
00130 bool ismrpc = hitst->is_mrpc();
00131 if( !barrel && !ismrpc ) { zrhit[0] = extTrk->tof1Position().rho(); }
00132 bool readout = hitst->is_readout();
00133 bool counter = hitst->is_counter();
00134 bool cluster = hitst->is_cluster();
00135 int layer = hitst->layer();
00136 tofid[layer-1] = (*it)->tofID();
00137 bool east = hitst->is_east();
00138 double tof = (*it)->tof();
00139 unsigned int ipmt = 0;
00140 if( readout ) {
00141
00142
00143 if( barrel ) { ipmt = ( ( st & 0xC0 ) >> 5 ) + ( ( ( st ^ 0x20 ) & 0x20 ) >> 5 ) - 2; }
00144 else {
00145 if( !ismrpc ) {
00146 if( tofid[0]<=47 ) { ipmt = 7; }
00147 else { ipmt = 8; }
00148 }
00149 else {
00150 if( tofid[0]<=35 ) { ipmt = 7; }
00151 else { ipmt = 8; }
00152 }
00153 }
00154
00155 for( unsigned int i=0; i<5; i++ ) {
00156 double offset = (*it)->toffset(i);
00157 double texp = (*it)->texp(i);
00158 if( texp<0.0 ) continue;
00159 double dt = tof - offset - texp;
00160 if( barrel ) {
00161 m_dt[i][ipmt] = dt;
00162 m_dtCorr[i][ipmt] = offsetTof( i, barrel, ipmt, betaGamma[i], charge, zrhit[layer-1], dt );
00163 m_sigCorr[i][ipmt] = sigmaTof( i, charge, barrel, ipmt, &tofid[0], &zrhit[0], betaGamma[i] );
00164 m_chiCorr[i][ipmt] = m_dtCorr[i][ipmt]/m_sigCorr[i][ipmt];
00165 }
00166 else {
00167 if( !ismrpc ) {
00168 m_dt[i][0] = dt;
00169 m_dtCorr[i][0] = offsetTof( i, barrel, ipmt, betaGamma[i], charge, zrhit[layer-1], dt );
00170 m_sigCorr[i][0] = sigmaTof( i, charge, barrel, ipmt, &tofid[0], &zrhit[0], betaGamma[i] );
00171 m_chiCorr[i][0] = m_dtCorr[i][ipmt]/m_sigCorr[i][ipmt];
00172 }
00173 else {
00174 m_dt[i][0] = dt;
00175 m_dtCorr[i][0] = dt;
00176 m_sigCorr[i][0] = 0.065;
00177 m_chiCorr[i][0] = m_dtCorr[i][0]/m_sigCorr[i][0];
00178 }
00179 }
00180 }
00181 if( counter && cluster ) {
00182 m_ipmt = ipmt;
00183 for( unsigned int i=0; i<5; i++ ) {
00184 if( ((*it)->texp(i))>0.0 ) {
00185 if( barrel ) {
00186 m_offset[i] = m_dtCorr[i][ipmt];
00187 m_sigma[i] = m_sigCorr[i][ipmt];
00188 }
00189 else {
00190 m_offset[i] = m_dtCorr[i][0];
00191 m_sigma[i] = m_sigCorr[i][0];
00192 }
00193 }
00194 }
00195 }
00196 }
00197 else {
00198 if( counter ) {
00199 ipmt = layer+3;
00200 for( unsigned int i=0; i<5; i++ ) {
00201 double offset = (*it)->toffset(i);
00202 double texp = (*it)->texp(i);
00203 if( texp<0.0 ) continue;
00204 double dt = tof - offset - texp;
00205 if( barrel ) {
00206 m_dt[i][ipmt] = dt;
00207 m_dtCorr[i][ipmt] = offsetTof( i, tofid[layer-1], zrhit[layer-1], betaGamma[i], charge, dt );
00208 m_sigCorr[i][ipmt] = sigmaTof( i, charge, barrel, layer+3, &tofid[0], &zrhit[0], betaGamma[i] );
00209 m_chiCorr[i][ipmt] = m_dtCorr[i][ipmt]/m_sigCorr[i][ipmt];
00210 }
00211 else {
00212 if( ismrpc ) {
00213 m_dt[i][0] = dt;
00214 m_dtCorr[i][0] = dt;
00215 m_sigCorr[i][0] = 0.065;
00216 m_chiCorr[i][0] = m_dtCorr[i][0]/m_sigCorr[i][0];
00217 }
00218 else {
00219 cout << "ParticleID: TofCorr::particleIDCalculation: Endcap Scintillator TOF have more than one end!!!" << endl;
00220 }
00221 }
00222 }
00223 if( cluster ) {
00224 m_ipmt = ipmt;
00225 for( unsigned int i=0; i<5; i++ ) {
00226 if( ((*it)->texp(i))>0.0 ) {
00227 if( barrel ) {
00228 m_offset[i] = m_dtCorr[i][ipmt];
00229 m_sigma[i] = m_sigCorr[i][ipmt];
00230 }
00231 else {
00232 m_offset[i] = m_dtCorr[i][0];
00233 m_sigma[i] = m_sigCorr[i][0];
00234 }
00235 }
00236 }
00237 }
00238 else {
00239 signal[layer-1] = correlationCheck( ipmt );
00240 }
00241 }
00242 else {
00243 if( cluster ) {
00244 ipmt = 6;
00245 for( unsigned int i=0; i<5; i++ ) {
00246 double offset = (*it)->toffset(i);
00247 double texp = (*it)->texp(i);
00248 if( texp<0.0 ) continue;
00249 double dt = tof - offset - texp;
00250 m_dt[i][ipmt] = dt;
00251 m_dtCorr[i][ipmt] = offsetTof( i, tofid[0], tofid[1], zrhit[0], zrhit[1], betaGamma[i], charge, dt );
00252 m_sigCorr[i][ipmt] = sigmaTof( i, charge, barrel, ipmt, &tofid[0], &zrhit[0], betaGamma[i] );
00253 m_chiCorr[i][ipmt] = m_dtCorr[i][ipmt]/m_sigCorr[i][ipmt];
00254 }
00255 if( signal[0] && signal[1] ) {
00256 m_ipmt = 6;
00257 for( unsigned int i=0; i<5; i++ ) {
00258 m_offset[i] = m_dtCorr[i][ipmt];
00259 m_sigma[i] = m_sigCorr[i][ipmt];
00260 }
00261 }
00262 else if( signal[0] && !signal[1] ) {
00263 m_ipmt = 4;
00264 for( unsigned int i=0; i<5; i++ ) {
00265 m_offset[i] = m_dtCorr[i][4];
00266 m_sigma[i] = m_sigCorr[i][4];
00267 }
00268 }
00269 else if( !signal[0] && signal[1] ) {
00270 m_ipmt = 5;
00271 for( unsigned int i=0; i<5; i++ ) {
00272 m_offset[i] = m_dtCorr[i][5];
00273 m_sigma[i] = m_sigCorr[i][5];
00274 }
00275 }
00276 else return irc;
00277 }
00278 }
00279 }
00280 }
00281
00282
00283 double pdftemp = 0;
00284 for( unsigned int i=0; i<5; i++ ) {
00285 m_chi[i] = m_offset[i]/m_sigma[i];
00286 if( m_chi[i]<0. && m_chi[i]>m_chimin ) { m_chimin = m_chi[i]; }
00287 if( m_chi[i]>0. && m_chi[i]<m_chimax ) { m_chimax = m_chi[i]; }
00288 double ppp = pdfCalculate(m_chi[i],1);
00289 if( fabs(ppp) > pdftemp) { pdftemp = fabs(ppp); }
00290 }
00291
00292 m_pdfmin = pdftemp;
00293 if( pdftemp < pdfCalculate(pdfMinSigmaCut(),1) ) return irc;
00294 if( ( m_chimin > -90.0 && (0-m_chimin)>chiMinCut() ) && ( m_chimax < 90.0 && m_chimax>chiMaxCut() ) ) return irc;
00295 for(int i = 0; i < 5; i++) {
00296 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
00297 }
00298
00299 irc = 0;
00300 return irc;
00301 }
00302
00303
00304 void TofCorrPID::inputParameter( int run ) {
00305
00306 std::string filePath = path + "/share/TofCorrPID/";
00307
00308 filePath = filePath + "jpsi2012";
00309 m_runBegin = 0;
00310 m_runEnd = 80000;
00311
00312 if( run>0 ) {
00313 filePath = filePath + "/data/";
00314 }
00315 else {
00316 filePath = filePath + "/mc/";
00317 }
00318
00319
00320
00321 std::string fileWeight = filePath + "calib_barrel_sigma.txt";
00322 ifstream inputWeight( fileWeight.c_str(), std::ios_base::in );
00323 if( !inputWeight ) {
00324 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileWeight << endl;
00325 exit(1);
00326 }
00327
00328 for( unsigned int tofid=0; tofid<176; tofid++ ) {
00329 for( unsigned int readout=0; readout<3; readout++ ) {
00330 for( unsigned int p_i=0; p_i<5; p_i++ ) {
00331 inputWeight >> m_p_weight[tofid][readout][p_i];
00332 }
00333 }
00334 }
00335
00336
00337
00338 std::string fileCommon = filePath + "calib_barrel_common.txt";
00339 ifstream inputCommon( fileCommon.c_str(), std::ios_base::in );
00340 if( !inputCommon ) {
00341 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileCommon << endl;
00342 exit(1);
00343 }
00344 inputCommon >> m_p_common;
00345
00346
00347
00348 std::string fileEcSigma = filePath + "calib_endcap_sigma.txt";
00349 ifstream inputEcSigma( fileEcSigma.c_str(), std::ios_base::in );
00350 if( !inputEcSigma ) {
00351 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileEcSigma << endl;
00352 exit(1);
00353 }
00354
00355 for( unsigned int tofid=0; tofid<96; tofid++ ) {
00356 for( unsigned int p_i=0; p_i<3; p_i++ ) {
00357 inputEcSigma >> m_ec_sigma[tofid][p_i];
00358 }
00359 }
00360
00361
00362
00363 std::string fileQ0BetaGamma = filePath + "curve_Q0_BetaGamma.txt";
00364 ifstream inputQ0BetaGamma( fileQ0BetaGamma.c_str(), std::ios_base::in );
00365 if( !inputQ0BetaGamma ) {
00366 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileQ0BetaGamma << endl;
00367 exit(1);
00368 }
00369
00370 for( unsigned int layer=0; layer<3; layer++ ) {
00371 for( unsigned int ipar=0; ipar<5; ipar++ ) {
00372 inputQ0BetaGamma >> m_q0_bg[layer][ipar];
00373 }
00374 }
00375
00376
00377
00378 std::string fileParAB = filePath + "parameter_A_B.txt";
00379 ifstream inputParAB( fileParAB.c_str(), std::ios_base::in );
00380 if( !inputParAB ) {
00381 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileParAB << endl;
00382 exit(1);
00383 }
00384
00385
00386
00387
00388
00389
00390 for( unsigned int ispecies=0; ispecies<2; ispecies++ ) {
00391 for( unsigned int ipmt=0; ipmt<5; ipmt++ ) {
00392 for( unsigned int icharge=0; icharge<2; icharge++ ) {
00393 for( unsigned int iab=0; iab<2; iab++ ) {
00394 for( unsigned int ipar=0; ipar<11; ipar++ ) {
00395 inputParAB >> m_par_ab_12[ispecies][ipmt][icharge][iab][ipar];
00396 }
00397 }
00398 }
00399 }
00400 }
00401
00402
00403 std::string fileSigma = filePath + "parameter_sigma.txt";
00404 ifstream inputSigma( fileSigma.c_str(), std::ios_base::in );
00405 if( !inputSigma ) {
00406 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileSigma << endl;
00407 exit(1);
00408 }
00409
00410
00411
00412
00413
00414 for( unsigned int ispecies=0; ispecies<4; ispecies++ ) {
00415 for( unsigned int ipmt=0; ipmt<8; ipmt++ ) {
00416 for( unsigned int ipar=0; ipar<12; ipar++ ) {
00417 inputSigma >> m_par_sigma[ispecies][ipmt][ipar];
00418 }
00419 }
00420 }
00421
00422
00423 std::string fileProtonOffset = filePath + "parameter_offset_proton.txt";
00424 ifstream inputProtonOffset( fileProtonOffset.c_str(), std::ios_base::in );
00425 if( !inputProtonOffset ) {
00426 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileProtonOffset << endl;
00427 exit(1);
00428 }
00429
00430
00431
00432
00433
00434
00435 for( unsigned int ispecies=0; ispecies<2; ispecies++ ) {
00436 for( unsigned int ipmt=0; ipmt<8; ipmt++ ) {
00437 for( unsigned int ipar=0; ipar<20; ipar++ ) {
00438 if( ipmt!=7 ) {
00439 for( unsigned int jpar=0; jpar<46; jpar++ ) {
00440 inputProtonOffset >> m_p_offset_12[ispecies][ipmt][ipar][jpar];
00441 }
00442 }
00443 else {
00444 for( unsigned int jpar=0; jpar<7; jpar++ ) {
00445 inputProtonOffset >> m_p_offset_ec12[ispecies][ipar][jpar];
00446 }
00447 }
00448 }
00449 }
00450 }
00451
00452
00453
00454 std::string fileProtonSigma = filePath + "parameter_sigma_proton.txt";
00455 ifstream inputProtonSigma( fileProtonSigma.c_str(), std::ios_base::in );
00456 if( !inputProtonSigma ) {
00457 cout << "ParticleID::TofCorrPID: Can NOT open file: " << fileProtonSigma << endl;
00458 exit(1);
00459 }
00460
00461
00462
00463
00464
00465
00466 for( unsigned int ispecies=0; ispecies<2; ispecies++ ) {
00467 for( unsigned int ipmt=0; ipmt<8; ipmt++ ) {
00468 for( unsigned int ipar=0; ipar<20; ipar++ ) {
00469 if( ipmt!=7 ) {
00470 for( unsigned int jpar=0; jpar<46; jpar++ ) {
00471 inputProtonSigma >> m_p_sigma_12[ispecies][ipmt][ipar][jpar];
00472 }
00473 }
00474 else {
00475 for( unsigned int jpar=0; jpar<7; jpar++ ) {
00476 inputProtonSigma >> m_p_sigma_ec12[ispecies][ipar][jpar];
00477 }
00478 }
00479 }
00480 }
00481 }
00482
00483
00484 return;
00485 }
00486
00487
00488 double TofCorrPID::offsetTof( unsigned int ispecies, bool barrel, unsigned int ipmt, double betaGamma, int charge, double zrhit, double dt ) {
00489 if( ispecies==0 || ispecies==1 ) { return dt; }
00490
00491 double deltaT = -1000.0;
00492 if( ( ipmt>= 4 && barrel ) || ( ipmt!=7 && ipmt!=8 && !barrel ) || betaGamma<0.0 || abs(charge)!=1 || fabs(zrhit)>120.0 ) {
00493 cout << "Particle::TofCorrPID: offsetTof for single end: the input parameter are NOT correct! Please check them!" << endl;
00494 return deltaT;
00495 }
00496 unsigned int layer=0;
00497 if( barrel ) {
00498 if( ipmt==0 || ipmt==1 ) { layer=0; }
00499 else if( ipmt==2 || ipmt==3 ) { layer=1; }
00500 }
00501 else { layer=2; }
00502 double q0 = qCurveFunc( layer, betaGamma );
00503
00504 unsigned int inumber=ipmt;
00505 if( !barrel ) { inumber=4; }
00506
00507 double func[9];
00508 for( unsigned int i=0; i<9; i++ ) {
00509 func[i]=0.0;
00510 }
00511
00512 double parA = 0.0;
00513 double parB = 0.0;
00514
00515
00516 func[0] = 1.0;
00517 for( unsigned int i=1; i<8; i++ ) {
00518 func[i] = func[i-1]*zrhit;
00519 }
00520
00521 unsigned int iparticle=0;
00522 if( ispecies==2 || ispecies==3 ) { iparticle=0; }
00523 else if( ispecies==4 ) { iparticle=1; }
00524 unsigned int icharge=0;
00525 if( charge==1 ) { icharge=0; }
00526 else if( charge==-1 ) { icharge=1; }
00527
00528 if( ispecies!=4 || betaGamma*xmass(4)>0.8 ) {
00529 for( unsigned int i=0; i<8; i++ ) {
00530 if( i<5 ) {
00531 parA += m_par_ab_12[iparticle][inumber][icharge][0][i]*func[i];
00532 parB += m_par_ab_12[iparticle][inumber][icharge][1][i]*func[i];
00533 }
00534 else if( i>=5 ) {
00535 parA += m_par_ab_12[iparticle][inumber][icharge][0][i+3]*func[i];
00536 parB += m_par_ab_12[iparticle][inumber][icharge][1][i+3]*func[i];
00537 }
00538 }
00539 for( unsigned int iab=0; iab<2; iab++ ) {
00540 func[8] = m_par_ab_12[iparticle][inumber][icharge][iab][5]*TMath::Gaus(zrhit,m_par_ab_12[iparticle][inumber][icharge][iab][6], m_par_ab_12[iparticle][inumber][icharge][iab][7]);
00541 if( iab==0 ) {
00542 parA += func[8];
00543 }
00544 else if( iab==1 ) {
00545 parB += func[8];
00546 }
00547 }
00548 }
00549
00550 double tcorr = parA + parB/sqrt(q0);
00551
00552
00553
00554 if( barrel && ispecies==4 && betaGamma*xmass(4)<0.8 ) {
00555 int nzbin = 46;
00556 double zbegin = -115.0;
00557 double zend = 115.0;
00558 double zstep = (zend-zbegin)/nzbin;
00559
00560 double nbgbin = 20.0;
00561 double bgbegin = 0.3;
00562 double bgend = 0.8;
00563 double bgstep;
00564 bgstep = (bgend-bgbegin)/nbgbin;
00565
00566 int izbin = static_cast<int>((zrhit-zbegin)/zstep);
00567 if( izbin<0 ) { izbin=0; }
00568 else if( izbin>=nzbin ) { izbin=nzbin-1; }
00569 int ibgbin = static_cast<int>((betaGamma*xmass(4)-bgbegin)/bgstep);
00570 if( ibgbin<0 ) { ibgbin=0; }
00571 else if( ibgbin>=nbgbin ) { ibgbin=nbgbin-1; }
00572
00573 if( charge==1 ) {
00574 tcorr = m_p_offset_12[0][ipmt][ibgbin][izbin];
00575 }
00576 else {
00577 tcorr = m_p_offset_12[1][ipmt][ibgbin][izbin];
00578 }
00579 }
00580 else if( !barrel && ispecies==4 && betaGamma*xmass(4)<0.8 ) {
00581 int nrbin = 7;
00582 double rbegin = 55.0;
00583 double rend = 83.0;
00584 double rstep = (rend-rbegin)/nrbin;
00585
00586 double nbgbin = 20.0;
00587 double bgbegin = 0.3;
00588 double bgend = 0.8;
00589 double bgstep;
00590 bgstep = (bgend-bgbegin)/nbgbin;
00591
00592 int irbin = static_cast<int>((zrhit-rbegin)/rstep);
00593 if( irbin<0 ) { irbin=0; }
00594 else if( irbin>=nrbin ) { irbin=nrbin-1; }
00595 int ibgbin = static_cast<int>((betaGamma*xmass(4)-bgbegin)/bgstep);
00596 if( ibgbin<0 ) { ibgbin=0; }
00597 else if( ibgbin>=nbgbin ) { ibgbin=nbgbin-1; }
00598
00599 if( charge==1 ) {
00600 tcorr = m_p_offset_ec12[0][ibgbin][irbin];
00601 }
00602 else {
00603 tcorr = m_p_offset_ec12[1][ibgbin][irbin];
00604 }
00605 }
00606
00607 deltaT = dt - tcorr;
00608
00609 return deltaT;
00610 }
00611
00612
00613 double TofCorrPID::offsetTof( unsigned int ispecies, int tofid, double zrhit, double betaGamma, int charge, double dt ) {
00614 if( ispecies==0 || ispecies==1 ) { return dt; }
00615
00616 double deltaT = -1000.0;
00617 if( tofid<0 || tofid>=176 ) {
00618 cout << "Particle::TofCorrPID: offsetTof for single layer: the input parameter are NOT correct! Please check them!" << endl;
00619 exit(1);
00620 }
00621 int pmt[3] = { -9, -9, -9 };
00622 if( tofid>=0 && tofid<=87 ) {
00623 pmt[0] = 0;
00624 pmt[1] = 1;
00625 pmt[2] = 4;
00626 }
00627 else {
00628 pmt[0] = 2;
00629 pmt[1] = 3;
00630 pmt[2] = 5;
00631 }
00632
00633 double sigmaCorr2 = m_p_common*m_p_common;
00634 double sigmaLeft = bSigma( 0, tofid, zrhit );
00635 double sigmaLeft2 = sigmaLeft*sigmaLeft;
00636 double sigmaRight = bSigma( 1, tofid, zrhit );
00637 double sigmaRight2 = sigmaRight*sigmaRight;
00638
00639 double fraction = ( sigmaRight2 - sigmaCorr2 )/( sigmaLeft2 + sigmaRight2 - 2.0*sigmaCorr2);
00640 deltaT = fraction*m_dtCorr[ispecies][pmt[0]] + (1.0-fraction)*m_dtCorr[ispecies][pmt[1]];
00641
00642
00643 double tcorr = 0.0;
00644 unsigned int ipmt = 4;
00645 if( tofid>=88 && tofid<176 ) { ipmt = 5; }
00646 if( ispecies==4 && betaGamma*xmass(4)<0.8 ) {
00647 int nzbin = 46;
00648 double zbegin = -115.0;
00649 double zend = 115.0;
00650 double zstep = (zend-zbegin)/nzbin;
00651
00652 double nbgbin = 20.0;
00653 double bgbegin = 0.3;
00654 double bgend = 0.8;
00655 double bgstep;
00656 bgstep = (bgend-bgbegin)/nbgbin;
00657
00658 int izbin = static_cast<int>((zrhit-zbegin)/zstep);
00659 if( izbin<0 ) { izbin=0; }
00660 else if( izbin>=nzbin ) { izbin=nzbin-1; }
00661 int ibgbin = static_cast<int>((betaGamma*xmass(4)-bgbegin)/bgstep);
00662 if( ibgbin<0 ) { ibgbin=0; }
00663 else if( ibgbin>=nbgbin ) { ibgbin=nbgbin-1; }
00664
00665 if( charge==1 ) {
00666 tcorr = m_p_offset_12[0][ipmt][ibgbin][izbin];
00667 }
00668 else {
00669 tcorr = m_p_offset_12[1][ipmt][ibgbin][izbin];
00670 }
00671 }
00672 deltaT = deltaT - tcorr;
00673
00674 return deltaT;
00675 }
00676
00677
00678 double TofCorrPID::offsetTof( unsigned int ispecies, int tofid1, int tofid2, double zrhit1, double zrhit2, double betaGamma, int charge, double dt ) {
00679 if( ispecies==0 || ispecies==1 ) { return dt; }
00680
00681 double deltaT = -1000.0;
00682 if( tofid1<0 || tofid1>=88 || tofid2<88 || tofid2>=176 ) {
00683 cout << "Particle::TofCorrPID: offsetTof for double layer: the input parameter are NOT correct! Please check them!" << endl;
00684 exit(1);
00685 }
00686
00687 double sigmaCorr2 = m_p_common*m_p_common;
00688 double sigmaInner = bSigma( 2, tofid1, zrhit1 );
00689 double sigmaInner2 = sigmaInner*sigmaInner;
00690 double sigmaOuter = bSigma( 2, tofid2, zrhit2 );
00691 double sigmaOuter2 = sigmaOuter*sigmaOuter;
00692 double sigma = sqrt( (sigmaInner2*sigmaOuter2-sigmaCorr2*sigmaCorr2)/(sigmaInner2+sigmaOuter2-2.0*sigmaCorr2) );
00693
00694 m_sigma[0] = sigma;
00695 m_sigma[1] = sigma;
00696
00697 double fraction = ( sigmaOuter2 - sigmaCorr2 )/( sigmaInner2 + sigmaOuter2 - 2.0*sigmaCorr2);
00698 deltaT = fraction*m_dtCorr[ispecies][4] + (1.0-fraction)*m_dtCorr[ispecies][5];
00699
00700
00701 double tcorr = 0.0;
00702 if( ispecies==4 && betaGamma*xmass(4)<0.8 ) {
00703 int nzbin = 46;
00704 double zbegin = -115.0;
00705 double zend = 115.0;
00706 double zstep = (zend-zbegin)/nzbin;
00707
00708 double nbgbin = 20.0;
00709 double bgbegin = 0.3;
00710 double bgend = 0.8;
00711 double bgstep;
00712 bgstep = (bgend-bgbegin)/nbgbin;
00713
00714 int izbin = static_cast<int>((zrhit1-zbegin)/zstep);
00715 if( izbin<0 ) { izbin=0; }
00716 else if( izbin>=nzbin ) { izbin=nzbin-1; }
00717 int ibgbin = static_cast<int>((betaGamma*xmass(4)-bgbegin)/bgstep);
00718 if( ibgbin<0 ) { ibgbin=0; }
00719 else if( ibgbin>=nbgbin ) { ibgbin=nbgbin-1; }
00720
00721 if( charge==1 ) {
00722 tcorr = m_p_offset_12[0][6][ibgbin][izbin];
00723 }
00724 else {
00725 tcorr = m_p_offset_12[1][6][ibgbin][izbin];
00726 }
00727 }
00728 deltaT = deltaT - tcorr;
00729
00730 return deltaT;
00731 }
00732
00733
00734 double TofCorrPID::sigmaTof( unsigned int ispecies, int charge, bool barrel, unsigned int ipmt, int tofid[2], double zrhit[2], double betaGamma ) {
00735
00736 double sigma = 1.0e-6;
00737
00738 if( ispecies==0 || ispecies==1 ) {
00739 if( barrel ) {
00740 if( ipmt==0 ) {
00741 sigma = bSigma( 0, tofid[0], zrhit[0] );
00742 }
00743 else if( ipmt==1 ) {
00744 sigma = bSigma( 1, tofid[0], zrhit[0] );
00745 }
00746 else if( ipmt==2 ) {
00747 sigma = bSigma( 0, tofid[1], zrhit[1] );
00748 }
00749 else if( ipmt==3 ) {
00750 sigma = bSigma( 1, tofid[1], zrhit[1] );
00751 }
00752 else if( ipmt==4 ) {
00753 sigma = bSigma( 2, tofid[0], zrhit[0] );
00754 }
00755 else if( ipmt==5 ) {
00756 sigma = bSigma( 2, tofid[1], zrhit[1] );
00757 }
00758 else if( ipmt==6 ) {
00759 sigma = bSigma( &tofid[0], &zrhit[0] );
00760 }
00761 }
00762 else {
00763 sigma = eSigma( tofid[0], zrhit[0] );
00764 }
00765 }
00766 else {
00767 unsigned int iz = 0;
00768 if( ipmt==2 || ipmt==3 || ipmt==5 ) { iz=1; }
00769
00770 sigma = sigmaTof( ispecies, charge, barrel, ipmt, zrhit[iz], betaGamma );
00771 }
00772
00773 return sigma;
00774 }
00775
00776
00777 double TofCorrPID::sigmaTof( unsigned int ispecies, int charge, bool barrel, unsigned int ipmt, double zrhit, double betaGamma ) {
00778
00779 int ibgbin = -1;
00780 int izbin = 0;
00781
00782 if( ispecies==4 && (betaGamma*xmass(4))<0.8 ) {
00783 double nbgbin = 20.0;
00784 double bgbegin = 0.3;
00785 double bgend = 0.8;
00786 double bgstep;
00787 bgstep = (bgend-bgbegin)/nbgbin;
00788 ibgbin = static_cast<int>((betaGamma*xmass(4)-bgbegin)/bgstep);
00789
00790 if( ibgbin<0 ) { ibgbin=0; }
00791 else if( ibgbin>=nbgbin ) { ibgbin=nbgbin-1; }
00792
00793 if( barrel ) {
00794 int nzbin = 46;
00795 double zbegin = -115.0;
00796 double zend = 115.0;
00797 double zstep = (zend-zbegin)/nzbin;
00798
00799 izbin = static_cast<int>((zrhit-zbegin)/zstep);
00800 if( izbin<0 ) { izbin=0; }
00801 else if( izbin>=nzbin ) { izbin=nzbin-1; }
00802 }
00803 else {
00804 int nzbin = 7;
00805 double zbegin = 55.0;
00806 double zend = 83.0;
00807 double zstep = (zend-zbegin)/nzbin;
00808
00809 izbin = static_cast<int>((zrhit-zbegin)/zstep);
00810 if( izbin<0 ) { izbin=0; }
00811 else if( izbin>=nzbin ) { izbin=nzbin-1; }
00812 }
00813 }
00814
00815 unsigned int numParam = 4;
00816 double func[12];
00817 for( unsigned int i=0; i<4; i++ ) {
00818 if( i==0 ) { func[i] = 1.0; }
00819 else {
00820 func[i] = func[i-1]*zrhit;
00821 }
00822 }
00823 for( unsigned int i=4; i<12; i++ ) {
00824 func[i] = 0.0;
00825 }
00826
00827
00828 if( barrel ) {
00829 if( ispecies==2 || ispecies==3 ) {
00830 for( unsigned int i=4; i<10; i++ ) {
00831 func[i] = func[i-1]*zrhit;
00832 }
00833 func[10] = 1.0/(115.0-zrhit)/(115.0-zrhit);
00834 func[11] = 1.0/(115.0+zrhit)/(115.0+zrhit);
00835 numParam = 12;
00836 }
00837 else if( ispecies==4 ) {
00838 for( unsigned int i=4; i<12; i++ ) {
00839 func[i] = func[i-1]*zrhit;
00840 }
00841 numParam = 12;
00842 }
00843 }
00844 else {
00845 numParam = 4;
00846 }
00847
00848 unsigned int inumber = ipmt;
00849 if( !barrel ) { inumber=7; }
00850
00851 double sigma = 0.0;
00852 if( ispecies==2 || ispecies==3 ) {
00853 for( unsigned int i=0; i<numParam; i++ ) {
00854 sigma += m_par_sigma[ispecies-2][inumber][i]*func[i];
00855 }
00856 }
00857 else if( ispecies==4 ) {
00858 if( ibgbin!=-1 ) {
00859
00860 if( barrel ) {
00861 if( charge==1 ) {
00862 sigma = m_p_sigma_12[0][inumber][ibgbin][izbin];
00863 }
00864 else {
00865 sigma = m_p_sigma_12[1][inumber][ibgbin][izbin];
00866 }
00867 }
00868 else {
00869 if( charge==1 ) {
00870 sigma = m_p_sigma_ec12[0][ibgbin][izbin];
00871 }
00872 else {
00873 sigma = m_p_sigma_ec12[1][ibgbin][izbin];
00874 }
00875 }
00876 if( fabs(sigma+999.0)<1.0e-6 ) {
00877 sigma = 0.001;
00878 }
00879 }
00880 else {
00881 for( unsigned int i=0; i<numParam; i++ ) {
00882 if( charge==1 ) {
00883 sigma += m_par_sigma[2][inumber][i]*func[i];
00884 }
00885 else {
00886 sigma += m_par_sigma[3][inumber][i]*func[i];
00887 }
00888 }
00889 }
00890 }
00891
00892
00893 int run = getRunNo();
00894 if( run>0 ) {
00895 if( ispecies==2 ) {
00896 sigma = sigma*(TMath::Exp((betaGamma-0.356345)/(-0.767124))+0.994072);
00897 }
00898 }
00899
00900 return sigma;
00901 }
00902
00903
00904 double TofCorrPID::qCurveFunc( unsigned int layer, double betaGamma ) {
00905 double q0 = -100.0;
00906 if( layer>=3 || betaGamma<0.0 ) {
00907 cout << "Particle::TofCorrPID::qCurveFunc: the input parameter are NOT correct! Please check them!" << endl;
00908 return q0;
00909 }
00910
00911 double beta = betaGamma/sqrt(1.0+betaGamma*betaGamma);
00912 double logterm = log( m_q0_bg[layer][2]+pow((1.0/betaGamma), m_q0_bg[layer][4] ) );
00913 q0 = m_q0_bg[layer][0]*( m_q0_bg[layer][1]-pow( beta, m_q0_bg[layer][3] ) - logterm )/pow( beta, m_q0_bg[layer][3] );
00914
00915 return q0;
00916 }
00917
00918
00919 double TofCorrPID::bSigma( unsigned int end, int tofid, double zrhit ) {
00920
00921 double func[5];
00922 func[0] = 1.0;
00923 func[1] = zrhit;
00924 func[2] = zrhit*zrhit;
00925 func[3] = zrhit*zrhit*zrhit;
00926 func[4] = zrhit*zrhit*zrhit*zrhit;
00927
00928 double sigma = 0.0;
00929 for( unsigned int i=0; i<5; i++ ) {
00930 sigma += m_p_weight[tofid][end][i]*func[i];
00931 }
00932
00933 return sigma;
00934 }
00935
00936
00937 double TofCorrPID::bSigma( int tofid[2], double zrhit[2] ) {
00938
00939 double sigma1 = bSigma( 2, tofid[0], zrhit[0] );
00940 double sigma2 = bSigma( 2, tofid[1], zrhit[1] );
00941 double sigmaCorr2 = m_p_common*m_p_common;
00942 double sigma = ( sigma1*sigma1*sigma2*sigma2 - sigmaCorr2*sigmaCorr2 )/( sigma1*sigma1 + sigma2*sigma2 - 2.0*sigmaCorr2 );
00943 sigma = sqrt(fabs(sigma));
00944
00945 return sigma;
00946 }
00947
00948
00949 double TofCorrPID::eSigma( int tofid, double zrhit ) {
00950
00951 double func[5];
00952 func[0] = 1.0;
00953 func[1] = zrhit;
00954 func[2] = zrhit*zrhit;
00955
00956 double sigma = 0.0;
00957 for( unsigned int i=0; i<3; i++ ) {
00958 sigma += m_ec_sigma[tofid][i]*func[i];
00959 }
00960
00961 return sigma;
00962 }
00963
00964
00965 bool TofCorrPID::correlationCheck( unsigned int ipmt ) {
00966 bool chiCut = false;
00967 bool good = false;
00968 double pdftemp = 0;
00969 for( unsigned int i=0; i<5; i++ ) {
00970 if( ( m_chiCorr[i][ipmt]>(0.0-chiMinCut()) && m_chiCorr[i][ipmt]<chiMaxCut() ) && !good ) { good=true; }
00971 double ppp = pdfCalculate(m_chiCorr[i][ipmt],1);
00972 if( fabs(ppp) > pdftemp) { pdftemp = fabs(ppp); }
00973 }
00974 m_pdfmin = pdftemp;
00975 if( pdftemp < pdfCalculate(pdfMinSigmaCut(),1) ) return chiCut;
00976 if( !good ) return chiCut;
00977 chiCut = true;
00978 return chiCut;
00979 }