00001 #include <cmath>
00002
00003 #include "ParticleID/TofEPID.h"
00004
00005 #ifndef BEAN
00006 #include "MdcRecEvent/RecMdcTrack.h"
00007 #include "TofRecEvent/RecTofTrack.h"
00008 #include "EvtRecEvent/EvtRecTrack.h"
00009 #include "DstEvent/TofHitStatus.h"
00010 #else
00011 #include "TofHitStatus.h"
00012 #endif
00013
00014 TofEPID * TofEPID::m_pointer = 0;
00015
00016 TofEPID * TofEPID::instance() {
00017 if(!m_pointer) m_pointer = new TofEPID();
00018 return m_pointer;
00019 }
00020
00021 TofEPID::TofEPID():ParticleIDBase() {
00022
00023 }
00024
00025 void TofEPID::init() {
00026 for(int i = 0; i < 5; i++) {
00027 m_chi[i] = 99.0;
00028 m_prob[i] = -1.0;
00029 m_offset[i] = 99.0;
00030 m_sigma[i] = 1.0;
00031 }
00032 m_chimin = 99.;
00033 m_pdfmin =99.;
00034 m_ndof = 0;
00035 m_mass2 = -999;
00036 m_rhit = -99;
00037 }
00038
00039 void TofEPID::calculate() {
00040 if(particleIDCalculation()==0) m_ndof=1;
00041 }
00042
00043 int TofEPID::particleIDCalculation() {
00044 int irc = -1;
00045 EvtRecTrack* recTrk = PidTrk();
00046 if(!(recTrk->isMdcTrackValid())) return irc;
00047 if(!(recTrk->isTofTrackValid())) return irc;
00048
00049 #ifndef BEAN
00050 SmartRefVector<RecTofTrack> tofTrackCol = recTrk->tofTrack();
00051 SmartRefVector<RecTofTrack>::iterator it;
00052 #else
00053 const std::vector<TTofTrack* >& tofTrackCol = recTrk->tofTrack();
00054 std::vector<TTofTrack* >::const_iterator it;
00055 #endif
00056
00057 bool isMRPC = false;
00058 if( tofTrackCol.size()!=1 && tofTrackCol.size()!=3 ) { return irc; }
00059 else {
00060 TofHitStatus *hitst = new TofHitStatus;
00061 if( tofTrackCol.size()==1 ) {
00062 it = tofTrackCol.begin();
00063 hitst->setStatus( (*it)->status() );
00064 bool isRaw = hitst->is_raw();
00065 if( isRaw ) { return irc; }
00066 isMRPC = hitst->is_mrpc();
00067 bool isReadout = hitst->is_readout();
00068 bool isCounter = hitst->is_counter();
00069 bool isCluster = hitst->is_cluster();
00070 if( !isReadout || !isCounter || !isCluster ) { return irc; }
00071 }
00072 else if( tofTrackCol.size()==3 ) {
00073 unsigned int icounter = -1;
00074 unsigned int inumber = 0;
00075 for( it = tofTrackCol.begin(); it != tofTrackCol.end(); it++, inumber++ ) {
00076 hitst->setStatus( (*it)->status() );
00077 bool isRaw = hitst->is_raw();
00078 if( isRaw ) continue;
00079 isMRPC = hitst->is_mrpc();
00080 if( !isMRPC ) { return irc; }
00081 bool isReadout = hitst->is_readout();
00082 bool isCounter = hitst->is_counter();
00083 bool isCluster = hitst->is_cluster();
00084 if( !isReadout && isCounter && isCluster ) { icounter = inumber; }
00085 }
00086 if( icounter == -1 ) { return irc; }
00087 it = tofTrackCol.begin() + icounter;
00088 }
00089 }
00090
00091 double tof = (*it)->tof();
00092 if( tof <= 0 ) { return irc; }
00093 double path = (*it)->path();
00094 m_rhit = (*it)->zrhit();
00095 double beta = path/velc()/tof;
00096 double beta2 = beta*beta;
00097 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00098 double ptrk = mdcTrk->p();
00099 m_mass2 = ptrk*ptrk*(1.0/beta2 -1.0);
00100
00101 double chitemp = 99.;
00102 double pdftemp = 0;
00103 for( unsigned int i=0; i<5; i++ ) {
00104 double texp = (*it)->texp(i);
00105 m_offset[i] = tof - texp - (*it)->toffset(i);
00106 double sigma_tmp = (*it)->sigma(0);
00107 if( !isMRPC ) {
00108 if( sigma_tmp>0 ) {
00109 if( i<2 ) { m_sigma[i] = sigma_tmp; }
00110 else { m_sigma[i] = 1.2*sigma_tmp; }
00111 }
00112 else { m_sigma[i] = 0.12; }
00113 }
00114 else {
00115 m_sigma[i] = 0.065;
00116 }
00117 m_chi[i] = m_offset[i]/m_sigma[i];
00118
00119 if( fabs(m_chi[i]) < chitemp ) { chitemp = fabs(m_chi[i]); }
00120 double ppp = pdfCalculate(m_chi[i],1);
00121 if( fabs(ppp) > pdftemp ) { pdftemp = fabs(ppp); }
00122 }
00123 m_chimin = chitemp;
00124
00125
00126
00127 for( unsigned int i=0; i<5; i++ ) {
00128 m_prob[i] = probCalculate( m_chi[i]*m_chi[i], 1 );
00129 }
00130
00131 irc = 0;
00132 return irc;
00133 }
00134
00135
00136
00137
00138
00139 double TofEPID::offsetTofE(int n, int cntr, double ptrk, double rtof, double ph,double charge) {
00140 double offset;
00141
00142 switch(n) {
00143 case 0: {
00144 double ptemp = ptrk;
00145 if(ptrk<0.2) ptemp = 0.2;
00146 if(ptrk > 2.1) ptemp = 2.1;
00147 double plog = log(ptemp);
00148 offset = 0.001*(-28.8481+138.159*plog-249.334*plog*plog);
00149 break;
00150 }
00151
00152 case 1: {
00153 double ptemp = ptrk;
00154 if(ptrk<0.2) ptemp = 0.2;
00155 if(ptrk > 2.1) ptemp = 2.1;
00156 double plog = log(ptemp);
00157 offset = 0.001*(-33.6966+1.91915*plog-0.592320*plog*plog);
00158 break;
00159 }
00160 case 2: {
00161 double ptemp = ptrk;
00162 if(ptrk<0.2) ptemp = 0.2;
00163 if(ptrk > 2.1) ptemp = 2.1;
00164 double plog = log(ptemp);
00165 offset = 0.001*(-27.9965 + 1.213 * plog - 2.02740 * plog * plog);
00166 break;
00167 }
00168 case 3: {
00169 double ptemp = ptrk;
00170 if(ptrk<0.3) ptemp = 0.3;
00171 if(ptrk > 2.1) ptemp = 2.1;
00172 double plog = log(ptemp);
00173 offset = 0.001*(-23.4842 -28.7569 * plog + 78.21* plog *plog);
00174 break;
00175 }
00176
00177 case 4: {
00178 double ptemp = ptrk;
00179 if(ptrk<0.4) ptemp = 0.4;
00180 if(ptrk > 2.1) ptemp = 2.1;
00181 double plog = log(ptemp);
00182 if(charge>0)
00183 offset = 0.001*(-4.854-110.540*plog+99.8732*plog*plog);
00184 if(charge<0)
00185 offset = 0.001*(27.047-145.120*plog+167.014*plog*plog);
00186 break;
00187 }
00188
00189 default:
00190 offset = 0.0;
00191 break;
00192 }
00193
00194 return offset;
00195 }
00196
00197 double TofEPID::sigmaTofE(int n, int cntr, double ptrk, double rtof, double ph,double charge) {
00198
00199 double sigma;
00200
00201 switch(n) {
00202
00203 case 0: {
00204 double ptemp = ptrk;
00205 if(ptrk < 0.2) ptemp = 0.2;
00206 if(ptrk > 2.1) ptemp = 2.1;
00207 double plog = log(ptemp);
00208 sigma = 0.001 * (109.974 +15.2457 * plog + 36.8139 * plog * plog);
00209
00210 break;
00211 }
00212
00213 case 1: {
00214 double ptemp = ptrk;
00215 if(ptrk < 0.2) ptemp = 0.2;
00216 if(ptrk > 2.1) ptemp = 2.1;
00217 double plog = log(ptemp);
00218 sigma = 0.001 * (96.5077 -2.96232 * plog + 3.12910 * plog * plog);
00219 break;
00220 }
00221
00222 case 2: {
00223 double ptemp = ptrk;
00224 if(ptrk < 0.2) ptemp = 0.2;
00225 if(ptrk > 2.1) ptemp = 2.1;
00226 double plog = log(ptemp);
00227 sigma = 0.001 * (105.447 - 2.08044 * plog + 3.44846 * plog * plog);
00228 break;
00229 }
00230
00231 case 3: {
00232 double ptemp = ptrk;
00233 if(ptrk < 0.3) ptemp = 0.3;
00234 if(ptrk > 2.1) ptemp = 2.1;
00235 double plog = log(ptemp);
00236 sigma = 0.001*(88.8806 - 26.8464 * plog + 113.672 * plog * plog);
00237 break;
00238 }
00239 case 4: {
00240 double ptemp = ptrk;
00241 if(ptrk < 0.5) ptemp = 0.5;
00242 if(ptrk > 2.1) ptemp = 2.1;
00243 double plog = log(ptemp);
00244 if(charge>0)
00245 sigma = 0.001 * (96.3534 -44.1139 * plog + 53.9805 * plog * plog);
00246 if(charge<0)
00247 sigma = 0.001 * (157.345 -98.7357 * plog + 55.1145 * plog * plog);
00248 break;
00249 }
00250
00251 default:
00252 sigma = 0.100;
00253
00254 break;
00255 }
00256
00257 return sigma;
00258 }
00259