00001 #include <cmath>
00002 #include "TMath.h"
00003
00004 #include "ParticleID/Tof1PID.h"
00005
00006 #ifndef BEAN
00007 #include "MdcRecEvent/RecMdcTrack.h"
00008 #include "TofRecEvent/RecTofTrack.h"
00009 #include "EvtRecEvent/EvtRecTrack.h"
00010 #include "DstEvent/TofHitStatus.h"
00011 #else
00012 #include "TofHitStatus.h"
00013 #endif
00014
00015 Tof1PID * Tof1PID::m_pointer = 0;
00016 Tof1PID * Tof1PID::instance() {
00017 if(!m_pointer) m_pointer = new Tof1PID();
00018 return m_pointer;
00019 }
00020
00021 Tof1PID::Tof1PID():ParticleIDBase() {
00022 m_pars[0]= -0.208289;
00023 m_pars[1]= 1.63092;
00024 m_pars[2]= -0.0733139;
00025 m_pars[3]= 1.02385;
00026 m_pars[4]= 0.00145052;
00027 m_pars[5]= -1.72471e-06;
00028 m_pars[6]= 5.92726e-10;
00029 m_pars[7]= -0.0645428;
00030 m_pars[8]= -0.00143637;
00031 m_pars[9]= -0.133817;
00032 m_pars[10]= 0.0101188;
00033 m_pars[11]= -5.07622;
00034 m_pars[12]= 5.31472;
00035 m_pars[13]= -0.443142;
00036 m_pars[14]= -12.1083;
00037
00038 }
00039
00040 void Tof1PID::init() {
00041 for(int i = 0; i < 5; i++) {
00042 m_chi[i] = 99.0;
00043 m_prob[i] = -1.0;
00044 m_sigma[i] = 1;
00045 m_offset[i] = 99.0;
00046 }
00047 m_chimin = 99.;
00048 m_pdfmin = 99.;
00049 m_ndof = 0;
00050 m_mass2 = -999;
00051 m_ph1 = -99;
00052 m_zhit1 = -99;
00053 }
00054
00055 void Tof1PID::calculate() {
00056 if(particleIDCalculation() == 0) m_ndof=1;
00057 }
00058
00059 int Tof1PID::particleIDCalculation() {
00060 int irc = -1;
00061
00062 EvtRecTrack* recTrk = PidTrk();
00063 if(!(recTrk->isMdcTrackValid())) return irc;
00064 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00065 double ptrk = mdcTrk->p();
00066 double charge = mdcTrk->charge();
00067
00068 if(!(recTrk->isTofTrackValid())) return irc;
00069
00070 #ifndef BEAN
00071 SmartRefVector<RecTofTrack> tofTrk = recTrk->tofTrack();
00072 SmartRefVector<RecTofTrack>::iterator it;
00073 #else
00074 const std::vector<TTofTrack* >& tofTrk = recTrk->tofTrack();
00075 std::vector<TTofTrack* >::const_iterator it;
00076 #endif
00077
00078 TofHitStatus *hitst = new TofHitStatus;
00079 std::vector<int> tof1count;
00080 int goodtof1trk=0;
00081 for(it = tofTrk.begin(); it!=tofTrk.end(); it++,goodtof1trk++) {
00082 unsigned int st = (*it)->status();
00083 hitst->setStatus(st);
00084 if( !(hitst->is_barrel()) ) continue;
00085 if( !(hitst->is_counter()) ) continue;
00086 if( hitst->layer()==1 ) tof1count.push_back(goodtof1trk);
00087 }
00088 delete hitst;
00089 if(tof1count.size()!=1) return irc;
00090 it = tofTrk.begin()+tof1count[0];
00091 double tof = (*it)->tof();
00092 m_tof1 = tof;
00093
00094
00095 int cntr = (*it)->tofID();
00096 double path = ((*it)->path())*10.0;
00097 m_path1 = path;
00098 m_ph1 = (*it)->ph();
00099 m_zhit1 = ((*it)->zrhit())*10;
00100 double beta2 = path*path/velc()/velc()/tof/tof;
00101 m_mass2 = ptrk * ptrk * (1/beta2 -1);
00102 if ((m_mass2>20)||(m_mass2<-1)) return irc;
00103 if(tof <=0 ) return irc;
00104 double chitemp = 99.;
00105 double pdftemp = 0;
00106 double sigma_tmp= (*it)->sigma(0);
00107 for(int i = 0; i < 5; i++) {
00108
00109 m_offset[i] = tof - (*it)->texp(i)-(*it)->toffset(i);
00110 if(sigma_tmp!=0) {
00111 m_sigma[i] = 1.2*sigma_tmp;
00112 if(i<2) m_sigma[i]=sigma_tmp;
00113 }
00114 else m_sigma[i]=sigmaTof1(i, cntr,ptrk,m_zhit1, m_ph1,charge);
00115 m_chi[i] = m_offset[i]/m_sigma[i];
00116
00117 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
00118 double ppp = pdfCalculate(m_chi[i],1);
00119 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
00120 }
00121 m_chimin = chitemp;
00122 m_pdfmin = pdftemp;
00123 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
00124 if(fabs(m_chimin) > chiMinCut()) return irc;
00125 for(int i = 0; i < 5; i++) {
00126 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
00127 }
00128 irc = 0;
00129 return irc;
00130 }
00131
00132
00133
00134
00135
00136
00137 double Tof1PID::offsetTof1(int n, int cntr, double ptrk, double z, double ph1,double charge) {
00138 double offset;
00139
00140 double betagamma;
00141 switch(n) {
00142
00143 case 0: {
00144 offset = 0.0;
00145 return offset;
00146 }
00147
00148 case 1: {
00149 offset = 0.0;
00150 return offset;
00151 }
00152
00153 case 2: {
00154 betagamma=ptrk/139.57;
00155 break;
00156 }
00157 case 3: {
00158 betagamma=ptrk/493.68;
00159 break;
00160 }
00161
00162 case 4: {
00163 betagamma=ptrk/938.27;
00164 break;
00165 }
00166
00167 default: {
00168 offset = 0.0;
00169 return offset;
00170 }
00171 break;
00172 }
00173 double Q = ph1;
00174 z=z/1000.0;
00175 double beta = betagamma / TMath::Sqrt(1 + betagamma*betagamma);
00176 double Q0 = sampleQ0(betagamma,beta);
00177 double func[15]= {0.};
00178 func[0]=1.;
00179 func[1]=1./sqrt(Q);
00180 func[2]=z/sqrt(Q);
00181 func[3]=z*z/sqrt(Q);
00182 func[4]=Q;
00183 func[5]=Q*Q;
00184 func[6]=Q*Q*Q;
00185 func[7]=1./(0.84*0.84+z*z);
00186 func[8]=z;
00187 func[9]=z*z;
00188 func[10]=z*z*z;
00189 func[11]=1./sqrt(Q0);
00190 func[12]=z/sqrt(Q0);
00191 func[13]=z*z/sqrt(Q0);
00192 func[14]=z*z*z/sqrt(Q0);
00193 offset=0.;
00194 for(int i=0; i<15; i++) {
00195 offset+= m_pars[i]*func[i];
00196 }
00197
00198 return offset;
00199 }
00200
00201 double Tof1PID::sigmaTof1(int n, int cntr, double ptrk, double ztof, double ph,double charge) {
00202 double sigma;
00203
00204 double cosz = cos(ztof/1000.0);
00205
00206 switch(n) {
00207
00208 case 0: {
00209 sigma = 0.132405-0.135638*cosz+0.105528*cosz*cosz;
00210 break;
00211 }
00212
00213 case 1: {
00214 sigma = 0.164057 -0.199648*cosz+ 0.139481*cosz*cosz;
00215 break;
00216 }
00217
00218 case 2: {
00219 sigma = 0.132943-0.0996678*cosz+0.0785578*cosz*cosz;
00220 break;
00221 }
00222 case 3: {
00223 sigma = 0.146572 -0.124672*cosz+0.0920656*cosz*cosz;
00224 break;
00225 }
00226
00227 case 4: {
00228 sigma = 0.1513 -0.0806081*cosz+0.0607409*cosz*cosz;
00229 break;
00230 }
00231
00232 default:
00233 sigma = 0.100;
00234 break;
00235 }
00236 sigma = 0.110;
00237
00238 return sigma;
00239 }
00240
00241
00242 double Tof1PID::sampleQ0(double betagamma,double beta) {
00243 double val = 0.261/ TMath::Power(beta, 1.96) *( 5.43- TMath::Power(beta, 1.96) -TMath::Log(1.12 + TMath::Power(1/betagamma, 0.386)));
00244 return val*100.;
00245 }