00001
00002
00003
00004
00005 #include "MdcDedxAlg/MdcDedxTrk.h"
00006 #include "MdcDedxAlg/MdcDedxCorrection.h"
00007
00008 #include "TGraphErrors.h"
00009 #include "TGraph.h"
00010 #include <stdlib.h>
00011 #include <search.h>
00012 #include <TMath.h>
00013
00014 MdcDedxTrk::MdcDedxTrk()
00015 {
00016 m_trk = 0;
00017 m_trk_kal = 0;
00018 m_stat = -1;
00019 m_trk_id = 0;
00020 m_quality = -99;
00021 m_charge = 0;
00022 m_P = 0;
00023 m_theta = 0;
00024 m_phi = 0;
00025 m_pl_rp = 0;
00026 m_nsample = 0;
00027
00028 m_dEdx = 0;
00029 m_truncate = 1;
00030 for( unsigned a = 0; a < 5 ; a++ )
00031 {
00032 dedx_exp[a]=0.0;
00033 ex_sigma[a]=0.0;
00034 pid_prob[a]=0.0;
00035 chi_ex[a]=999.0;
00036 }
00037
00038 #ifdef DEBUG
00039 std::cout<<"MdcDedxTrk(1) constructed!"<<std::endl;
00040 #endif
00041 }
00042
00043 MdcDedxTrk::MdcDedxTrk( RecMdcTrack& trk )
00044 {
00045 m_trk = 0;
00046 m_trk_kal = 0;
00047 m_stat = -1;
00048 m_trk_id = 0;
00049 m_quality = -99;
00050 m_charge = 0;
00051 m_P = 0;
00052 m_theta = 0;
00053 m_phi = 0;
00054 m_pl_rp = 0;
00055 m_nsample = 0;
00056 m_dEdx = 0;
00057 m_truncate = 0.7;
00058
00059 set_ExTrk( trk );
00060 #ifdef DEBUG
00061 std::cout<<"MdcDedxTrk(2) constructed!"<<std::endl;
00062 #endif
00063 }
00064
00065 MdcDedxTrk::MdcDedxTrk( RecMdcKalTrack& trk_kal , int pid)
00066 {
00067 m_trk = 0;
00068 m_trk_kal = 0;
00069 m_stat = -1;
00070 m_trk_id = 0;
00071 m_quality = -99;
00072 m_charge = 0;
00073 m_P = 0;
00074 m_theta = 0;
00075 m_phi = 0;
00076 m_pl_rp = 0;
00077 m_nsample = 0;
00078 m_dEdx = 0;
00079 m_truncate = 0.7;
00080
00081 set_ExTrk_Kal( trk_kal, pid );
00082 #ifdef DEBUG
00083 std::cout<<"MdcDedxTrk(2) kal constructed!"<<std::endl;
00084 #endif
00085 }
00086
00087 MdcDedxTrk::~MdcDedxTrk()
00088 {
00089 #ifdef DEBUG
00090 std::cout<<"MdcDedxTrk destructed!!!"<<std::endl;
00091 #endif
00092 }
00093
00094 void MdcDedxTrk::set_ExTrk( RecMdcTrack& trk )
00095 {
00096 for( unsigned a = 0; a < 5 ; a++ )
00097 {
00098 dedx_exp[a]=0.0;
00099 ex_sigma[a]=0.0;
00100 pid_prob[a]=0.0;
00101 chi_ex[a]=0.0;
00102 }
00103 m_stat = 1;
00104 m_trk = &trk;
00105 m_trk_id = trk.trackId();
00106 m_quality = trk.stat();
00107
00108 m_charge = ( trk.helix(2) > 0 )? 1 : -1;
00109 float m_Pt = 1.0/fabs( trk.helix(2) );
00110 m_P = m_Pt*sqrt(1 + trk.helix(4)*trk.helix(4));
00111 m_theta = M_PI_2 - atan( trk.helix(4) );
00112 m_phi = ( trk.helix(1) < 3*M_PI_2 )? trk.helix(1)+M_PI_2 : trk.helix(1)-3*M_PI_2;
00113 m_nsample = trk.getNhits();
00114 m_pl_rp = 1.5;
00115
00116
00117 #ifdef DEBUG
00118 std::cout<<"MdcDedxTrk::set_ExTrk(&trk)!!!"<<std::endl;
00119 #endif
00120 }
00121
00122 void MdcDedxTrk::set_ExTrk_Kal( RecMdcKalTrack& trk_kal, int pid)
00123 {
00124 DstMdcKalTrack* dstTrk = &trk_kal;
00125 for( unsigned a = 0; a < 5 ; a++ )
00126 {
00127 dedx_exp[a]=0.0;
00128 ex_sigma[a]=0.0;
00129 pid_prob[a]=0.0;
00130 chi_ex[a]=0.0;
00131 }
00132 if(pid<0 || pid>4) pid = 2;
00133 m_stat = 1;
00134 m_trk_kal = &trk_kal;
00135 m_trk_id = trk_kal.trackId();
00136 m_quality = dstTrk->getStat(pid);
00137
00138 HepVector kal_helix = dstTrk->getFHelix(pid);
00139 m_charge = ( kal_helix[2] > 0 )? 1 : -1;
00140 float m_Pt = 1.0/fabs( kal_helix[2] );
00141 m_P = m_Pt*sqrt(1 + kal_helix[4]*kal_helix[4]);
00142 m_theta = M_PI_2 - atan( kal_helix[4] );
00143 m_phi = ( kal_helix[1] < 3*M_PI_2 )? kal_helix[1]+M_PI_2 : kal_helix[1]-3*M_PI_2;
00144 m_nsample = trk_kal.getNhits(pid);
00145 m_pl_rp = 1.5;
00146
00147 #ifdef DEBUG
00148 std::cout<<"MdcDedxTrk::set_ExTrk_Kal(&trk_kal)!!!"<<std::endl;
00149 #endif
00150 }
00151
00152
00153 double MdcDedxTrk::cal_dedx( float truncate )
00154 {
00155 m_truncate = truncate;
00156 sort(m_phlist.begin(),m_phlist.end());
00157 int nsampl = (int)( m_phlist.size()*truncate );
00158 double qSum = 0;
00159 unsigned int i = 0;
00160 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++)
00161 {
00162 i++;
00163 if(i<= nsampl) qSum += (*ql);
00164 }
00165
00166 float trunc=qSum/nsampl;
00167 return trunc;
00168 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00169 }
00170
00171 double MdcDedxTrk::cal_dedx_bitrunc(float truncate, int alg, int & usedhit )
00172 {
00173 m_truncate = truncate;
00174 sort(m_phlist.begin(),m_phlist.end());
00175 int nsampl = (int)( m_phlist.size()*truncate );
00176 int smpl = (int)(m_phlist.size()*(truncate+0.05));
00177 int min_cut = (int)( m_phlist.size()*0.05 + 0.5 );
00178 double qSum = 0;
00179 unsigned i = 0;
00180 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++)
00181 {
00182 i++;
00183 if(i<= smpl && i>=min_cut ) qSum += (*ql);
00184 }
00185 double trunc=-99;
00186 usedhit = smpl-min_cut+1;
00187 if(usedhit>0) trunc=qSum/usedhit;
00188
00189 return trunc;
00190 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00191 }
00192
00193 double MdcDedxTrk::cal_dedx_median( float truncate )
00194 {
00195 sort(m_phlist.begin(),m_phlist.end());
00196
00197 int nsample = m_phlist.size();
00198 double median;
00199 if( fmod(double(nsample),2.0) ) {
00200 median = m_phlist[(nsample-1)/2];
00201 } else {
00202 median= 0.5*( m_phlist[nsample/2] + m_phlist[nsample/2-1] );
00203 }
00204 return median;
00205 }
00206
00207 double MdcDedxTrk::cal_dedx_geometric( float truncate )
00208 {
00209 sort(m_phlist.begin(),m_phlist.end());
00210
00211 int nsampl = m_phlist.size();
00212 double qSum = 1.0;
00213 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00214 qSum *= (*ql);
00215 }
00216
00217 double trunc = pow(qSum,1/double(nsampl));
00218 return trunc;
00219 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00220 }
00221
00222 double MdcDedxTrk::cal_dedx_geometric_trunc( float truncate )
00223 {
00224 m_truncate = truncate;
00225 sort(m_phlist.begin(),m_phlist.end());
00226 int nsampl = (int)( m_phlist.size()*truncate );
00227 double qSum = 1.0;
00228 unsigned i = 0;
00229 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00230 i++;
00231 if( i<= nsampl )
00232 qSum *= (*ql);
00233 }
00234
00235 double trunc = pow(qSum,1/double(nsampl));
00236 return trunc;
00237
00238 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00239 }
00240
00241 double MdcDedxTrk::cal_dedx_harmonic( float truncate )
00242 {
00243 sort(m_phlist.begin(),m_phlist.end());
00244
00245 int nsampl = m_phlist.size();
00246 double qSum = 0;
00247 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00248 qSum += 1/(*ql);
00249 }
00250
00251 float trunc=nsampl/qSum;
00252 return trunc;
00253 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00254 }
00255
00256 double MdcDedxTrk::cal_dedx_harmonic_trunc( float truncate )
00257 {
00258 m_truncate = truncate;
00259 sort(m_phlist.begin(),m_phlist.end());
00260 int nsampl = (int)( m_phlist.size()*truncate );
00261 double qSum = 0;
00262 unsigned i = 0;
00263 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00264 i++;
00265 if( i<= nsampl )
00266 qSum += 1/(*ql);
00267 }
00268
00269 float trunc= nsampl/qSum;
00270 return trunc;
00271 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00272 }
00273
00274 double MdcDedxTrk::cal_dedx_transform( int alg )
00275 {
00276 sort(m_phlist.begin(),m_phlist.end());
00277
00278 int nsampl = m_phlist.size();
00279 double qSum = 0;
00280 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00281 qSum += 1/sqrt(*ql);
00282 }
00283
00284 float trunc=1/qSum;
00285
00286 return trunc;
00287 std::cout<<"MdcDedxTrk::cal_dedx()!!!"<<std::endl;
00288 }
00289
00290
00291 double MdcDedxTrk::cal_dedx_log( float truncate ,int alg)
00292 {
00293 double qSum = 0;
00294 unsigned i = 0;
00295 for(vector<double>::iterator ql= m_phlist.begin();ql!=m_phlist.end();ql++){
00296 i++;
00297 qSum += log(*ql);
00298 }
00299
00300 float trunc=qSum/m_phlist.size();
00301 return trunc;
00302 std::cout<<"MdcDedxTrk::cal_dedx_log()!!!"<<std::endl;
00303 }
00304
00305
00306
00307 void MdcDedxTrk::set_dEdx( int landau, float dEdx, int trkalg, int runflag, int vflag[3] , double t0, vector<double>& DedxCurve_Parameter, vector<double>& DedxSigma_Parameter, MdcDedxCorrection* ex_calib)
00308 {
00309 #ifdef DEBUG
00310 cout<<"in MdcDedxTrk::set_dEdx() landau: "<<landau<<" dedx: "<<dEdx<<" nsample(): "<<nsample()<<" ph size: "<<m_phlist.size()<<" m_P: "<<m_P<<" theta: "<<m_theta<<" pl-rp: "<<m_pl_rp<<endl;
00311 #endif
00312
00313 m_dEdx = dEdx;
00314 int dedxhit_use = (int)(m_phlist.size()*m_truncate);
00315
00316
00317 if(runflag ==1 || runflag ==2 )
00318 ex_calib->dedx_pid_exp_old( landau, runflag, (float)dEdx, (int)dedxhit_use,
00319 (float)m_P, (float)m_theta, (float)t0,(float)m_pl_rp,
00320 dedx_exp, ex_sigma, pid_prob, chi_ex);
00321
00322 else
00323 ex_calib->dedx_pid_exp( vflag, (float)dEdx, trkalg,(int)dedxhit_use,
00324 (float)m_P, (float)m_theta, (float)t0,(float)m_pl_rp,
00325 dedx_exp, ex_sigma, pid_prob, chi_ex, DedxCurve_Parameter, DedxSigma_Parameter);
00326
00327 #ifdef DEBUG
00328 std::cout<<"MdcDedxTrk::set_dEdx()!!!"<<std::endl;
00329 #endif
00330 }
00331
00332
00333 double MdcDedxTrk::SpaceChargeCorrec(double m_theta, double mom, int Particle, double dEdx )
00334 {
00335 const int par_cand( 5 );
00336 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
00337 double beta_G;
00338 double e_Par[5] = {143.349, 1.7315, 0.192616, 2.90437, 1.08248};
00339 double Beta_Gamma[22] ={0.373026, 0.479605, 0.586184, 0.692763, 0.799342, 782.779, 1565.56,
00340 2348.34, 17.2727, 18.1245, 1.43297, 2.14946, 12.1803, 13.6132, 6.62515, 10.4109,
00341 14.1967, 17.9825, 21.7683, 26.0274, 30.7596, 35.4919 };
00342 double K_par[22] ={4.64411e-05, 5.86544e-05, 8.05289e-05, 8.46981e-05, 8.92014e-05, 4.74517e-05,
00343 4.51684e-05, 5.32732e-05, 6.12803e-05, 6.14592e-05, 8.08608e-05, 6.73184e-05, 5.46448e-05,
00344 6.1377e-05, 6.57385e-05, 7.03053e-05, 6.61171e-05, 6.86824e-05, 6.246e-05, 7.25988e-05,
00345 7.11034e-05, 6.24924e-05 };
00346 double D_par[22] ={0.0871504, 0.0956379, 0.117193, 0.118647, 0.127203, 0.0566449, 0.0529198,
00347 0.0642525, 0.0764562, 0.081341, 0.0952263, 0.0987536, 0.0639901, 0.0845994,0.0777062,
00348 0.0823206, 0.0783874, 0.079537, 0.0815792, 0.0849875, 0.0824751,0.0776296 };
00349 double DSqr_par[22] = {0.00759519, 0.0091466, 0.0137341, 0.0140772, 0.0161807, 0.00320864,
00350 0.00280051, 0.00412839, 0.00584555, 0.00661636, 0.00906805, 0.00975227, 0.00409473,
00351 0.00715706, 0.00603826, 0.00677668, 0.00614458, 0.00632613, 0.00665516, 0.00722288,
00352 0.00680214, 0.00602635};
00353
00354 beta_G = mom/Charge_Mass[Particle];
00355 if(beta_G <0.3) beta_G =0.3;
00356 double bet=beta_G/TMath::Sqrt(beta_G*beta_G+1);
00357 double fterm=TMath::Log(e_Par[2]+1/pow(beta_G,e_Par[4]));
00358 double fitval=e_Par[0]/pow(bet,e_Par[3])*(e_Par[1]-pow(bet,e_Par[3])-fterm);
00359 TGraphErrors *gr1 = new TGraphErrors(22,Beta_Gamma, K_par,0,0);
00360 TGraphErrors *gr2 = new TGraphErrors(22,Beta_Gamma, DSqr_par,0,0);
00361
00362 double par[3];
00363 par[0] = fitval;
00364 par[1] = gr1->Eval(m_theta);
00365 par[2] = gr2->Eval(m_theta);
00366 Double_t y = fabs(cos(m_theta));
00367 double electron_par[3] ={334.032, 6.20658e-05, 0.00525673};
00368 double arg= TMath::Sqrt(y*y+ par[2]);
00369
00370 double cal_factor =TMath::Exp(-(par[1]* par[0])/arg);
00371 double arg_electron = TMath::Sqrt(y*y + electron_par[2]);
00372
00373 double electron_factor = TMath::Exp(-(electron_par[1]* electron_par[0])/arg_electron);
00374
00375 double dedx_cal = dEdx/(cal_factor/electron_factor);
00376
00377
00378 delete gr1;
00379 delete gr2;
00380 return dedx_cal;
00381 }