00001 #include <iostream>
00002 #include <cmath>
00003 #include <cstdlib>
00004
00005 #include "ParticleID/ParticleIDBase.h"
00006 #include "TMath.h"
00007
00008 const int ParticleIDBase::USE_DEDX = 1;
00009 const int ParticleIDBase::USE_TOF1 = 2;
00010 const int ParticleIDBase::USE_TOF2 = 4;
00011 const int ParticleIDBase::USE_TOFE = 8;
00012 const int ParticleIDBase::USE_TOFQ = 16;
00013 const int ParticleIDBase::USE_EMC = 32;
00014 const int ParticleIDBase::USE_MUC = 64;
00015 const int ParticleIDBase::USE_TOF = 128;
00016 const int ParticleIDBase::USE_TOFC = 256;
00017 const int ParticleIDBase::USE_TOFCorr = 512;
00018
00019 const int ParticleIDBase::IDENTIFY_ELECTRON = 1;
00020 const int ParticleIDBase::IDENTIFY_MUON = 2;
00021 const int ParticleIDBase::IDENTIFY_PION = 4;
00022 const int ParticleIDBase::IDENTIFY_KAON = 8;
00023 const int ParticleIDBase::IDENTIFY_PROTON = 16;
00024
00025 const int ParticleIDBase::PROBABILITY_PID = 1;
00026 const int ParticleIDBase::LIKELIHOOD_PID = 2;
00027 const int ParticleIDBase::NEURONNETWORK_PID = 4;
00028
00029 const int ParticleIDBase::DEDX_VALID = 1;
00030 const int ParticleIDBase::TOF1_VALID = 2;
00031 const int ParticleIDBase::TOF2_VALID = 4;
00032 const int ParticleIDBase::TOFE_VALID = 8;
00033 const int ParticleIDBase::TOFQ_VALID = 16;
00034 const int ParticleIDBase::EMC_VALID = 32;
00035 const int ParticleIDBase::MUC_VALID = 64;
00036 const int ParticleIDBase::TOF_VALID = 128;
00037 const int ParticleIDBase::TOFC_VALID = 256;
00038 const int ParticleIDBase::TOFCorr_VALID = 512;
00039
00040
00041 std::string ParticleIDBase::path = "";
00042
00043 ParticleIDBase::ParticleIDBase()
00044 {
00045 m_trk = 0;
00046 m_chimin_cut = 4;
00047 m_chimax_cut = 6;
00048 m_pdfsigmamin_cut=99;
00049
00050 #ifndef BEAN
00051 if( path.empty() ) set_path(0);
00052 #endif
00053 }
00054
00055 void ParticleIDBase::set_path(const char* s_path)
00056 {
00057 if ( s_path ) {
00058 path = string(s_path);
00059 } else {
00060 char* env_path = getenv("PARTICLEIDROOT");
00061 if ( !env_path ) {
00062 cout << " ParticleIDBase::set_path ERROR:"
00063 " the environment PARTICLEIDROOT not defined " << endl;
00064 exit(1);
00065 }
00066 path = string(env_path);
00067 }
00068 }
00069
00070
00071 double ParticleIDBase::xmass(int n) {
00072 double mass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00073 if(n < 0 || n >=5) return 0.0;
00074 return mass[n];
00075 }
00076
00077 double ParticleIDBase::velc() {
00078
00079 double vel = 299.792458;
00080 return vel;
00081 }
00082
00083 double ParticleIDBase::probCalculate(double chi2, int ndof) {
00084 double p = -1.0;
00085 if(chi2 < 0) return p;
00086 p = TMath::Prob(chi2, ndof);
00087 return p;
00088 }
00089
00090
00091 double ParticleIDBase::pdfCalculate(double offset,double sigma) {
00092
00093 const double pi = M_PI;
00094 const double twoPi = 2*pi;
00095 double chi2 = -0.5*offset*offset/(sigma*sigma);
00096 double pdf = exp(chi2)/(sigma*sqrt(twoPi));
00097 return pdf;
00098 }
00099
00100 double ParticleIDBase::p() {
00101 double val = 999;
00102 if(!m_trk) return val;
00103 if(!m_trk->isMdcTrackValid()) return val;
00104 RecMdcTrack *mdcTrk = m_trk->mdcTrack();
00105 val = mdcTrk->p();
00106 return val;
00107 }
00108 double ParticleIDBase::pt() {
00109 double val = 999;
00110 if(!m_trk) return val;
00111 if(!m_trk->isMdcTrackValid()) return val;
00112 RecMdcTrack *mdcTrk = m_trk->mdcTrack();
00113 val = mdcTrk->pxy();
00114 return val;
00115 }
00116 double ParticleIDBase::charge() {
00117 double val = 999;
00118 if(!m_trk) return val;
00119 if(!m_trk->isMdcTrackValid()) return val;
00120 RecMdcTrack *mdcTrk = m_trk->mdcTrack();
00121 val = mdcTrk->charge() + 0.0;
00122 return val;
00123 }
00124
00125 double ParticleIDBase::interpolation(double* x, double* y,double x1) {
00126 double c1 = (y[0]-y[1])*(x[1]-x[2])-(x[0]-x[1])*(y[1]-y[2]);
00127 double c2 = (x[0]*x[0]-x[1]*x[1])*(x[1]-x[2])-(x[1]*x[1]-x[2]*x[2])*(x[0]-x[1]);
00128 double c = c1/c2;
00129 double b1 = (y[0]-y[1])*(x[1]*x[1]-x[2]*x[2])-(x[0]*x[0]-x[1]*x[1])*(y[1]-y[2]);
00130 double b2 = (x[0]-x[1])*(x[1]*x[1]-x[2]*x[2])-(x[1]-x[2])*(x[0]*x[0]-x[1]*x[1]);
00131 double b = b1/b2;
00132 double a = y[0] - b*x[0]-c*x[0]*x[0];
00133 double y1 = a + b*x1 +c*x1*x1;
00134 return y1;
00135 }
00136
00137 double ParticleIDBase::pol2(double x, double* par) {
00138 double y=x;
00139
00140 return par[0] + y*(par[1] + y*(par[2]));
00141 }
00142
00143 double ParticleIDBase::pol3(double x, double* par) {
00144 double y=x;
00145
00146 return par[0] + y*(par[1] + y*(par[2] + y*(par[3])));
00147 }
00148
00149 double ParticleIDBase::pol4(double x, double* par) {
00150 double y=x;
00151
00152 return par[0] + y*(par[1] + y*(par[2] + y*(par[3] + y*(par[4]))));
00153 }