00001 #include <cmath>
00002 #include "TMath.h"
00003
00004 #include "ParticleID/Tof2PID.h"
00005
00006 #ifndef BEAN
00007 #include "EvtRecEvent/EvtRecTrack.h"
00008 #include "MdcRecEvent/RecMdcTrack.h"
00009 #include "TofRecEvent/RecTofTrack.h"
00010 #include "DstEvent/TofHitStatus.h"
00011 #else
00012 #include "TofHitStatus.h"
00013 #endif
00014
00015 Tof2PID * Tof2PID::m_pointer = 0;
00016
00017 Tof2PID * Tof2PID::instance() {
00018 if(!m_pointer) m_pointer = new Tof2PID();
00019 return m_pointer;
00020 }
00021
00022
00023 Tof2PID::Tof2PID():ParticleIDBase() {
00024 m_pars[0]=-0.237207;
00025 m_pars[1]= 1.90436;
00026 m_pars[2]= -0.210625;
00027 m_pars[3]= 0.664667;
00028 m_pars[4]= 0.00165226;
00029 m_pars[5]= -1.86503e-06;
00030 m_pars[6]= 6.07045e-10;
00031 m_pars[7]=-0.0882228;
00032 m_pars[8]= 0.0125708;
00033 m_pars[9]= -0.117157;
00034 m_pars[10]= 0.00252878;
00035 m_pars[11]= 0.254343;
00036 m_pars[12]= -5.74886;
00037 m_pars[13]= 5.20507;
00038 m_pars[14]= 1.86515;
00039
00040 }
00041
00042 void Tof2PID::init() {
00043
00044 for(int i = 0; i < 5; i++) {
00045 m_chi[i] = 99.0;
00046 m_prob[i] = -1.0;
00047 m_sigma[i] = 1.0;
00048 m_offset[i] =99.0;
00049 }
00050 m_chimin = 99.;
00051 m_pdfmin =99.;
00052 m_ndof = 0;
00053 m_mass2 = -999;
00054 m_ph2 = -99;
00055 m_zhit2 = -99;
00056 }
00057
00058 void Tof2PID::calculate() {
00059 if(particleIDCalculation() == 0) m_ndof=1;
00060 }
00061
00062 int Tof2PID::particleIDCalculation() {
00063 int irc = -1;
00064 EvtRecTrack* recTrk = PidTrk();
00065 if(!(recTrk->isMdcTrackValid())) return irc;
00066 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00067
00068 double ptrk = mdcTrk->p();
00069
00070 double charge = mdcTrk->charge();
00071
00072 if(!(recTrk->isTofTrackValid())) return irc;
00073
00074 #ifndef BEAN
00075 SmartRefVector<RecTofTrack> tofTrk = recTrk->tofTrack();
00076 SmartRefVector<RecTofTrack>::iterator it;
00077 #else
00078 const std::vector<TTofTrack* >& tofTrk = recTrk->tofTrack();
00079 std::vector<TTofTrack* >::const_iterator it;
00080 #endif
00081
00082 TofHitStatus *hitst = new TofHitStatus;
00083 std::vector<int> tof2count;
00084 int goodtof2trk=0;
00085 for(it = tofTrk.begin(); it!=tofTrk.end(); it++,goodtof2trk++) {
00086 unsigned int st = (*it)->status();
00087 hitst->setStatus(st);
00088
00089 if( !(hitst->is_barrel()) ) continue;
00090 if( !(hitst->is_counter()) ) continue;
00091 if( hitst->layer()==2 ) tof2count.push_back(goodtof2trk);
00092 }
00093 delete hitst;
00094
00095 if(tof2count.size()!=1) return irc;
00096 it = tofTrk.begin()+tof2count[0];
00097 double tof = (*it)->tof();
00098 m_tof2 = tof;
00099
00100 if(tof <=0 ) return irc;
00101
00102
00103 int cntr = (*it)->tofID();
00104 double path = ((*it)->path())*10.0;
00105
00106 m_path2 = path;
00107
00108
00109
00110
00111 m_ph2 = (*it)->ph();
00112 m_zhit2 = ((*it)->zrhit())*10;
00113
00114 double beta2 = path*path/velc()/velc()/tof/tof;
00115 m_mass2 = ptrk * ptrk * (1/beta2 -1);
00116 if ((m_mass2>20)||(m_mass2<-1)) return irc;
00117
00118 double chitemp = 99.;
00119 double pdftemp = 0;
00120 double sig_tmp = (*it)->sigma(0);
00121 for(int i = 0; i < 5; i++) {
00122 m_offset[i] = tof - (*it)->texp(i)-(*it)->toffset(i);
00123 if(sig_tmp!=0) {
00124 m_sigma[i] = 1.2*sig_tmp;
00125 if(i<2) m_sigma[i]=sig_tmp;
00126 }
00127 else m_sigma[i]= sigmaTof2(i, cntr,ptrk,m_zhit2, m_ph2,charge);
00128 m_chi[i] = m_offset[i]/m_sigma[i];
00129 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
00130 double ppp = pdfCalculate(m_chi[i],1);
00131 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
00132 }
00133
00134 m_chimin = chitemp;
00135 m_pdfmin = pdftemp;
00136 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
00137 if(fabs(m_chimin) > chiMinCut()) return irc;
00138
00139
00140
00141 for(int i = 0; i < 5; i++) {
00142 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
00143 }
00144 irc = 0;
00145 return irc;
00146 }
00147
00148
00149
00150
00151
00152
00153 double Tof2PID::offsetTof2(int n, int cntr, double ptrk, double z, double ph,double charge) {
00154 double offset;
00155
00156 double betagamma;
00157 switch(n) {
00158
00159 case 0: {
00160 offset = 0.0;
00161 return offset;
00162 }
00163
00164 case 1: {
00165 offset = 0.0;
00166 return offset;
00167 }
00168
00169 case 2: {
00170 betagamma=ptrk/139.57;
00171 break;
00172 }
00173 case 3: {
00174 betagamma=ptrk/493.68;
00175 break;
00176 }
00177
00178 case 4: {
00179 betagamma=ptrk/938.27;
00180 break;
00181 }
00182
00183 default: {
00184 offset = 0.0;
00185 return offset;
00186 }
00187 break;
00188 }
00189 z=z/1000.0;
00190 double Q = ph;
00191 double beta = betagamma / TMath::Sqrt(1 + betagamma*betagamma);
00192 double Q0 = sampleQ0(betagamma,beta);
00193 double func[15]= {0.};
00194 func[0]=1.;
00195 func[1]=1./sqrt(Q);
00196 func[2]=z/sqrt(Q);
00197 func[3]=z*z/sqrt(Q);
00198 func[4]=Q;
00199 func[5]=Q*Q;
00200 func[6]=Q*Q*Q;
00201 func[7]=1./(0.89*0.89+z*z);
00202 func[8]=z;
00203 func[9]=z*z;
00204 func[10]=z*z*z;
00205 func[11]=1./sqrt(Q0);
00206 func[12]=z/sqrt(Q0);
00207 func[13]=z*z/sqrt(Q0);
00208 func[14]=z*z*z/sqrt(Q0);
00209 offset=0.;
00210 for(int i=0; i<15; i++) {
00211 offset+= m_pars[i]*func[i];
00212 }
00213
00214 return offset;
00215 }
00216
00217 double Tof2PID::sigmaTof2(int n, int cntr, double ptrk, double ztof, double ph,double charge) {
00218 double sigma;
00219
00220 double cosz = cos(ztof/1000.0);
00221
00222 switch(n) {
00223
00224 case 0: {
00225 sigma = 0.115816 -0.0933215*cosz+ 0.0788252*cosz*cosz;
00226 break;
00227 }
00228
00229 case 1: {
00230 sigma = 0.121536 -0.0935898*cosz+0.0748771*cosz*cosz;
00231 break;
00232 }
00233
00234 case 2: {
00235 sigma = 0.106882-0.0147375*cosz+0.027491*cosz*cosz;
00236 break;
00237 }
00238 case 3: {
00239 sigma =0.106882 -0.0147375*cosz +0.027491 *cosz*cosz;
00240 break;
00241 }
00242
00243 case 4: {
00244 sigma = 0.168489 -0.131762*cosz+ 0.105708*cosz*cosz;
00245 break;
00246 }
00247
00248 default:
00249 sigma = 0.100;
00250 break;
00251 }
00252 sigma = 0.110;
00253
00254 return sigma;
00255 }
00256
00257
00258 double Tof2PID::sampleQ0(double betagamma,double beta) {
00259 double p1 = 0.261;
00260 double p2 = 5.43;
00261 double p3 = 1.12;
00262 double p4 = 1.96;
00263 double p5 = 0.386;
00264 double val = p1 / TMath::Power(beta, p4) *
00265 ( p2- TMath::Power(beta, p4) -
00266 TMath::Log(p3 + TMath::Power(1/betagamma, p5)));
00267
00268 return val*100.;
00269 }
00270