#include <HoughValidUpdate.h>
Public Types | |
typedef std::map< int, int > | M2 |
typedef std::map< int, M2 > | M1 |
Public Member Functions | |
HoughValidUpdate (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
StatusCode | beginRun () |
Private Member Functions | |
int | GetMcInfo () |
int | digiToHots (MdcDigiVec mdcDigiVec, TrkRecoTrk *newTrack, vector< bool > vec_truthHit, uint32_t getDigiFlag, vector< double > vec_flt_truth, map< int, double > hit_to_circle1, vector< MdcHit * > &) |
int | digiToHots2 (MdcDigiVec mdcDigiVec, TrkRecoTrk *newTrack, vector< bool > vec_truthHit, vector< double > vec_flt_truth, uint32_t getDigiFlag, vector< int > vec_slant, vector< double > vec_l_Calcu, vector< double > vec_l_axial, vector< double > vec_theta_ontrack, vector< int > &vec_hitbeforefit, map< int, double > hit_to_circle2, vector< MdcHit * > &, double) |
double | CFtrans (double x, double y) |
int | SlantId (int layer) |
int | LeastSquare (int n, vector< double > n_x, vector< double > n_y, vector< int > n_slant, vector< int > n_layer, vector< bool > vec_truthHit, double &d0, double &phi0, double &omega, double &circleX, double &circleY, double &circleR) |
int | Zposition (int n, vector< int > n_slant, double x_cirtemp, double y_cirtemp, double r_temp, vector< double > n_x_east, vector< double > n_x_west, vector< double > n_y_east, vector< double > n_y_west, vector< double > n_z_east, vector< double > n_z_west, vector< int > n_layer, vector< int > n_wire, vector< double > n_z, vector< double > &z_stereo_aver, vector< double > &l, vector< bool > vec_truthHit) |
void | Linefit (int z_stereoNum, vector< double > z_stereo_aver, vector< double > l, double &tanl, double &z0) |
void | Multi_array (int binX, int binY) |
void | FillCells (TH2D *h1, int n, int binX, int binY, vector< double > vec_u, vector< double > vec_v, vector< int > vec_layer, vector< int > vec_wire, vector< vector< int > > &countij, vector< vector< vector< int > > > &vec_selectNum, vector< vector< vector< int > > > &vec_selectHit) |
void | RhoTheta (int m_nHit, vector< double > vec_u, vector< double > vec_v, vector< double > &vec_rho, vector< double > &vec_theta, vector< vector< int > > &vec_hitNum, vector< int > vec_digitruth) |
void | RhoThetaAxial (vector< int > vec_slant, int nHit, vector< double > vec_u, vector< double > vec_v, vector< double > &vec_rho, vector< double > &vec_theta, vector< vector< int > > &vec_hitNum, vector< int > vec_digitruth) |
void | FillHist (TH2D *h1, vector< double > vec_u, vector< double > vec_v, int m_nHit, vector< vector< vector< int > > > &vec_selectNum) |
void | FillRhoTheta (TH2D *h1, vector< double > vec_theta, vector< double > vec_rho) |
int | zPosition (MdcHit &h, double t0, double x_cirtemp, double y_cirtemp, double r_temp, int digiId, double &delta_z, double &delta_l, double &l_Calcu_Left, double &l_Calcu_Right, double &z_Calcu_Left, double &z_Calcu_Right, double &z_Calcu, double &l_Calcu, int &i_q) |
void | AmbigChoose (vector< vector< vector< double > > > point, int n_use, vector< int > &ambigList, double &tanl, double &z0) |
void | Combine (TH2D *h1, int widthX, int widthY, vector< vector< vector< int > > > vec_selectNum, M1 &hit_combine, M2 &peak_combine_num, int &ntrack, vector< int > peakList1[]) |
void | Combine_two_cell (TH2D *h1, int xPeakCell, int yPeakCell, int xAround, int yaround, vector< vector< vector< int > > > vec_selectNum, int m, M1 &hit_combine, M2 &peak_combine_num) |
int | LeastLine (int, double x[], double y[], double &, double &, double &) |
Private Attributes | |
int | m_drawTuple |
bool operator<(const HoughValidUpdate::Mytrack& a,const HoughValidUpdate::Mytrack& b); | |
int | m_useMcInfo |
int | m_par_use |
int | m_method |
int | m_debug |
int | m_data |
int | m_binx |
int | m_biny |
int | m_peakCellNumX |
int | m_peakCellNumY |
int | m_ndev1 |
int | m_ndev2 |
double | m_fpro |
double | m_hit_pro |
std::string | m_pdtFile |
bool | m_pickHits |
double | m_fltcut |
double | m_nHitcut_low |
double | m_nHitcut_high |
int | m_eventcut |
double | m_rhocut |
int | m_setpeak_fit |
int | m_mcpar |
int | m_fillTruth |
double | m_bunchT0 |
int | t_eventNum |
int | t_runNum |
double | m_maxrho |
int | binX |
int | binY |
double | d0_mc |
double | phi0_mc |
double | omega_mc |
double | dz_mc |
double | tanl_mc |
int | track_fit |
uint32_t | m_getDigiFlag |
int | m_maxMdcDigi |
bool | m_keepBadTdc |
bool | m_dropHot |
bool | m_keepUnmatch |
int | m_minMdcDigi |
int | m_q |
double | m_dropTrkDrCut |
double | m_dropTrkDzCut |
double | m_dropTrkPtCut |
double | m_dropTrkChi2Cut |
std::string | m_configFile |
const MdcDetector * | m_gm |
HepPDT::ParticleDataTable * | m_particleTable |
TrkContextEv * | m_context |
BField * | m_bfield |
const MdcCalibFunSvc * | m_mdcCalibFunSvc |
RawDataProviderSvc * | m_rawDataProviderSvc |
MdcGeomSvc * | m_mdcGeomSvc |
MdcPrintSvc * | m_mdcPrintSvc |
int | t_eventNo |
IMagneticFieldSvc * | m_pIMF |
int | m_pid |
int | nfailure |
bool | m_combineTracking |
std::vector< float > | m_helixHitsSigma |
double | m_qualityFactor |
bool | m_removeBadTrack |
int | m_dropTrkNhitCut |
double | m_dropTrkChi2NdfCut |
NTuple::Tuple * | ntuplehit |
NTuple::Item< int > | m_eventNum |
NTuple::Item< int > | m_eventTime |
NTuple::Item< int > | m_runNum |
NTuple::Item< int > | m_nHit_Mc |
NTuple::Item< int > | m_nHitAxial_Mc |
NTuple::Array< int > | m_layer_Mc |
NTuple::Array< int > | m_cell_Mc |
NTuple::Item< int > | m_nTrkMC |
NTuple::Array< double > | m_pidTruth |
NTuple::Array< double > | m_costaTruth |
NTuple::Array< double > | m_ptTruth |
NTuple::Array< double > | m_pzTruth |
NTuple::Array< double > | m_pTruth |
NTuple::Array< double > | m_qTruth |
NTuple::Array< double > | m_drTruth |
NTuple::Array< double > | m_phi0Truth |
NTuple::Array< double > | m_omegaTruth |
NTuple::Array< double > | m_dzTruth |
NTuple::Array< double > | m_tanl_mc |
NTuple::Item< int > | m_nHit |
NTuple::Item< int > | m_nHitDigi |
NTuple::Item< int > | m_nHitAxial |
NTuple::Item< int > | m_nHitStereo |
NTuple::Array< int > | m_layer |
NTuple::Array< int > | m_cell |
NTuple::Array< double > | m_x_east |
NTuple::Array< double > | m_y_east |
NTuple::Array< double > | m_x_west |
NTuple::Array< double > | m_y_west |
NTuple::Array< double > | m_z_east |
NTuple::Array< double > | m_z_west |
NTuple::Array< double > | m_x |
NTuple::Array< double > | m_y |
NTuple::Array< double > | m_z |
NTuple::Array< double > | m_u |
NTuple::Array< double > | m_v |
NTuple::Array< double > | m_p |
NTuple::Array< int > | m_slant |
NTuple::Item< int > | m_layer_max |
NTuple::Array< int > | m_layerNhit |
NTuple::Array< double > | m_ltruth |
NTuple::Array< double > | m_ztruth |
NTuple::Array< double > | m_postruth |
NTuple::Array< int > | m_digi_truth |
NTuple::Item< int > | m_e_delta |
NTuple::Item< int > | m_nCross |
NTuple::Array< double > | m_rho |
NTuple::Array< double > | m_theta |
NTuple::Item< int > | m_xybinNum |
NTuple::Array< int > | m_xybin |
NTuple::Item< double > | m_xsigma |
NTuple::Item< double > | m_ysigma |
NTuple::Item< int > | m_npeak_1 |
NTuple::Item< int > | m_npeak_2 |
NTuple::Item< int > | m_trackFailure |
NTuple::Item< int > | m_npeak_after_tracking |
NTuple::Array< int > | m_nHitSelect |
NTuple::Array< int > | m_nHitAxialSelect |
NTuple::Array< int > | m_nHitStereoSelect |
NTuple::Array< int > | m_axiallose |
NTuple::Array< int > | m_stereolose |
NTuple::Item< int > | m_trackNum_Mc |
NTuple::Item< int > | m_trackNum |
NTuple::Item< int > | m_trackNum_fit |
NTuple::Item< int > | m_trackNum_tds |
NTuple::Array< double > | m_x_circle |
NTuple::Array< double > | m_y_circle |
NTuple::Array< double > | m_r |
NTuple::Array< double > | m_chis_pt |
NTuple::Array< double > | m_d0 |
NTuple::Array< double > | m_phi0 |
NTuple::Array< double > | m_omega |
NTuple::Array< double > | m_z0 |
NTuple::Array< double > | m_tanl |
NTuple::Array< double > | m_pt |
NTuple::Array< double > | m_pt2 |
NTuple::Array< double > | m_pz |
NTuple::Array< double > | m_pxyz |
NTuple::Array< int > | m_nFitFailure |
NTuple::Array< int > | m_2d_nFit |
NTuple::Array< int > | m_2dFit |
NTuple::Array< int > | m_3d_nFit |
NTuple::Array< int > | m_3dFit |
NTuple::Array< double > | m_pt2_rec_final |
NTuple::Array< double > | m_p_rec_final |
NTuple::Array< double > | m_d0_final |
NTuple::Array< double > | m_phi0_final |
NTuple::Array< double > | m_omega_final |
NTuple::Array< double > | m_z0_final |
NTuple::Array< double > | m_tanl_final |
NTuple::Array< double > | m_chis_3d_final |
NTuple::Array< int > | m_Q |
NTuple::Array< double > | m_pt2_rec_tds |
NTuple::Array< double > | m_p_rec_tds |
NTuple::Array< double > | m_d0_tds |
NTuple::Array< double > | m_phi0_tds |
NTuple::Array< double > | m_omega_tds |
NTuple::Array< double > | m_z0_tds |
NTuple::Array< double > | m_tanl_tds |
NTuple::Array< double > | m_chis_3d_tds |
NTuple::Array< int > | m_Q_tds |
double | m_truthPos [43][288][3] |
NTuple::Array< double > | m_dist_to_circle |
NTuple::Item< int > | m_nHitStereo_use |
NTuple::Array< double > | m_lCalcuLeft |
NTuple::Array< double > | m_lCalcuRight |
NTuple::Array< double > | m_zCalcuLeft |
NTuple::Array< double > | m_zCalcuRight |
NTuple::Array< double > | m_zCalcu |
NTuple::Array< double > | m_lCalcu |
NTuple::Array< double > | m_delta_z |
NTuple::Array< double > | m_delta_l |
NTuple::Array< int > | m_amb |
NTuple::Item< int > | m_ambig_no_match |
NTuple::Item< double > | m_tanl_Cal |
NTuple::Item< double > | m_z0_Cal |
NTuple::Item< double > | m_ambig_no_match_propotion |
NTuple::Item< int > | m_nhitdis0 |
NTuple::Item< int > | m_nhitdis1 |
NTuple::Item< int > | m_nhitdis2 |
NTuple::Array< double > | m_hit_dis0 |
NTuple::Array< double > | m_hit_dis1 |
NTuple::Array< double > | m_hit_dis2 |
NTuple::Array< double > | m_hit_dis1_3d |
NTuple::Array< double > | m_hit_dis2_3d |
NTuple::Item< bool > | m_decay |
NTuple::Item< int > | m_outputTrk |
NTuple::Array< int > | m_hitOnTrk |
NTuple::Item< int > | m_track_use_intrk1 |
NTuple::Item< int > | m_noise_intrk1 |
Classes | |
struct | Mytrack |
Definition at line 29 of file HoughValidUpdate.h.
typedef std::map<int, M2> HoughValidUpdate::M1 |
Definition at line 37 of file HoughValidUpdate.h.
typedef std::map<int, int> HoughValidUpdate::M2 |
Definition at line 36 of file HoughValidUpdate.h.
HoughValidUpdate::HoughValidUpdate | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 51 of file HoughValidUpdate.cxx.
References m_binx, m_biny, m_combineTracking, m_data, m_debug, m_drawTuple, m_dropHot, m_dropTrkChi2Cut, m_dropTrkChi2NdfCut, m_dropTrkDrCut, m_dropTrkDzCut, m_dropTrkNhitCut, m_dropTrkPtCut, m_eventcut, m_fillTruth, m_getDigiFlag, m_hit_pro, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_maxrho, m_mcpar, m_method, m_minMdcDigi, m_ndev1, m_ndev2, m_pdtFile, m_peakCellNumX, m_peakCellNumY, m_pickHits, m_pid, m_q, m_qualityFactor, m_removeBadTrack, m_rhocut, m_setpeak_fit, and m_useMcInfo.
00051 : 00052 Algorithm(name, pSvcLocator) 00053 { 00054 // Declare the properties 00055 declareProperty("drawTuple", m_drawTuple=0); 00056 declareProperty("useMcInfo", m_useMcInfo=0); 00057 declareProperty("method", m_method=0); 00058 declareProperty("debug", m_debug = 0); 00059 declareProperty("data", m_data= 0); 00060 declareProperty("binx", m_binx= 100); 00061 declareProperty("biny", m_biny= 200); 00062 declareProperty("maxrho", m_maxrho= 0.3); 00063 declareProperty("peakCellNumX", m_peakCellNumX); 00064 declareProperty("peakCellNumY", m_peakCellNumY); 00065 declareProperty("ndev1", m_ndev1= 1); 00066 declareProperty("ndev2", m_ndev2= 5); 00067 declareProperty("f_hit_pro", m_hit_pro= 0.4); 00068 declareProperty("pdtFile", m_pdtFile = "pdt.table"); 00069 declareProperty("pickHits", m_pickHits = true); // ?? 00070 declareProperty("pid", m_pid = 0); 00071 declareProperty("combineTracking",m_combineTracking = false); 00072 declareProperty("getDigiFlag", m_getDigiFlag = 0); 00073 declareProperty("maxMdcDigi", m_maxMdcDigi = 0); 00074 declareProperty("keepBadTdc", m_keepBadTdc = 0); 00075 declareProperty("dropHot", m_dropHot= 0); 00076 declareProperty("keepUnmatch", m_keepUnmatch = 0); 00077 declareProperty("minMdcDigi", m_minMdcDigi = 0); 00078 declareProperty("setnpeak_fit", m_setpeak_fit= -1); 00079 declareProperty("Q", m_q= 0); 00080 declareProperty("eventcut", m_eventcut= -1); 00081 declareProperty("rhocut", m_rhocut= -1.); 00082 declareProperty("mcpar", m_mcpar= 0); //if use par truth 00083 declareProperty("fillTruth", m_fillTruth= 0); 00084 declareProperty("removeBadTrack", m_removeBadTrack= false); 00085 declareProperty("dropTrkDrCut", m_dropTrkDrCut= 1.); 00086 declareProperty("dropTrkDzCut", m_dropTrkDzCut= 10.); 00087 declareProperty("dropTrkPtCut", m_dropTrkPtCut= 99999.); 00088 declareProperty("dropTrkChi2Cut", m_dropTrkChi2Cut = 99999.); 00089 declareProperty("dropTrkChi2NdfCut", m_dropTrkChi2NdfCut = 99999.); 00090 declareProperty("qualityFactor", m_qualityFactor= 3.); 00091 declareProperty("dropTrkNhitCut", m_dropTrkNhitCut= 9); 00092 }
void HoughValidUpdate::AmbigChoose | ( | vector< vector< vector< double > > > | point, | |
int | n_use, | |||
vector< int > & | ambigList, | |||
double & | tanl, | |||
double & | z0 | |||
) | [private] |
Definition at line 2689 of file HoughValidUpdate.cxx.
References abs, EvtCyclic3::combine(), genRecEmupikp::i, ganga-rec::j, LeastLine(), and x.
Referenced by execute().
02689 { 02690 //connect two points 02691 int n =n_use; 02692 int ambig_list[4]; 02693 // TF1 *fpol1 =new TF1("fpol1","pol1",7,17); 02694 double Chis_Least; 02695 //TGraph *tgr[4]; 02696 double Chis[4]; 02697 int ambig_combine=0; 02698 // three points chis 02699 double chi_k[4]; 02700 double chi_b[4]; 02701 double first_x[3]; 02702 double first_y[3]; 02703 first_x[0]=0; 02704 first_y[0]=0; 02705 vector< vector <int> > select(4,vector<int>() ); 02706 for(int icombine1=0;icombine1<2;icombine1++){ 02707 for(int icombine2=0;icombine2<2;icombine2++){ 02708 int combine=icombine1*2+icombine2; 02709 first_x[1]=point[icombine1][0][n-1]; 02710 first_x[2]=point[icombine2][0][n-2]; 02711 first_y[1]=point[icombine1][1][n-1]; 02712 first_y[2]=point[icombine2][1][n-2]; 02713 double k,b; 02714 double chi2=0; 02715 LeastLine(3,first_x,first_y,k,b,chi2); 02716 // TGraph *gr_first=new TGraph(3,first_x,first_y); 02717 //gr_first->Fit("fpol1","QN"); 02718 //delete gr_first; 02719 //double k=fpol1->GetParameter(1); 02720 //double b=fpol1->GetParameter(0); 02721 //cout<<"K: "<<k<<" B: "<<b<<endl; 02722 02723 // sprintf(hname,"h%d",combine); 02724 // h[combine]=new TH2D("hname","hname",100,7,17,100,-5,35); 02725 // double k=(point[icombine1][1][n-1]-point[icombine2][1][n-2])/(point[icombine1][0][n-1]-point[icombine2][0][n-2]); //left-left 02726 // double b=point[icombine1][1][n-1]-k*point[icombine1][0][n-1]; 02727 // cout<<"k: "<<k<<"b: "<<b<<endl; 02728 02729 for(int iHit=0;iHit<n-2;iHit++){ 02730 double dist_to_line[2]; 02731 for(int j=0;j<2;j++){ 02732 double x=point[j][0][iHit]; 02733 double y=point[j][1][iHit]; 02734 dist_to_line[j]=( abs(k*x-y+b) ) / sqrt( k*k +1); 02735 } 02736 if(dist_to_line[0]>dist_to_line[1]) select[combine].push_back(1); 02737 else select[combine].push_back(0); 02738 // cout<<select[combine][iHit]<<endl; 02739 } 02740 select[combine].push_back(icombine2); //left -left 02741 select[combine].push_back(icombine1); 02742 02743 //least 2* fit 02744 //add (0,0) 02745 // vector<double> l; 02746 // vector<double> z; 02747 double *l_new=new double[n+1]; 02748 double *z_new=new double[n+1]; 02749 for(int i=0;i<n;i++){ 02750 int ambig=select[combine][i]; 02751 l_new[i]=point[ambig][0][i]; 02752 z_new[i]=point[ambig][1][i]; 02753 } 02754 l_new[n]=0; 02755 z_new[n]=0; 02756 //for(int i=0;i<n;i++){ 02757 // int ambig=select[combine][i]; 02758 // cout<<"ambig : "<<ambig<<endl; 02759 // l.push_back(point[ambig][0][i]); 02760 // z.push_back(point[ambig][1][i]); 02761 // cout<<"hit: "<<i<<" :("<<l[i]<<","<<z[i]<<")"<<endl; 02762 //} 02763 //l.push_back(0); 02764 //z.push_back(0); 02765 02766 double k_least,b_least; 02767 double chi2_least=0; 02768 LeastLine(n+1,l_new,z_new,k_least,b_least,chi2_least); 02769 Chis[combine]=chi2_least; 02770 chi_k[combine]=k_least; 02771 chi_b[combine]=b_least; 02772 02773 // tgr[combine]=new TGraph(n+1,l_new,z_new); 02774 // tgr[combine]->Fit("fpol1","QN"); 02775 delete []l_new; 02776 delete []z_new; 02777 02778 02779 // TH2D *h1=new TH2D("h1","h1",100,7,17,100,-5,35); 02780 //for(int i =0 ;i<n+1; i++){ 02781 // h[combine]->Fill(l[i],z[i]); 02782 // //h->Draw("SAME"); 02783 //} 02784 // h[combine]->Fit("fpol1"); 02785 //Chis[combine]=fpol1->GetChisquare(); 02786 //chi_k[combine]=fpol1->GetParameter(1); 02787 //chi_b[combine]=fpol1->GetParameter(0); 02788 // cout<<" chis: "<<Chis[combine]<<endl; 02789 //delete tgr[combine]; 02790 // delete h[combine]; 02791 } //icombine1 02792 } //icombine2 02793 // for(int i =0;i<4;i++){ 02794 // delete h[i]; 02795 // } 02796 // delete fpol1; 02797 Chis_Least=Chis[0]; 02798 for(int i=0;i<4;i++){ 02799 if(Chis_Least>=Chis[i]) {Chis_Least=Chis[i]; ambig_combine=i;/* cout<<"least i: "<< i<<endl;*/} 02800 } 02801 for(int iHit=0;iHit<n;iHit++){ 02802 ambigList.push_back(select[ambig_combine][iHit]); 02803 // cout<<"ambig_combine: "<<ambig_combine<<endl; 02804 // cout<<" iHit: "<<select[ambig_combine][iHit]<<endl; 02805 // cout<<" iHit: "<<ambigList[iHit]<<endl; 02806 } 02807 02808 tanl=chi_k[ambig_combine]; 02809 z0=chi_b[ambig_combine]; 02810 // cout<<"Chi least: "<<Chis_Least<<" ambig: "<<ambig_combine<<endl; 02811 // cout<<"ambig_combine: "<<ambig_list[ambig_combine]<<endl; 02812 }
StatusCode HoughValidUpdate::beginRun | ( | ) |
Definition at line 94 of file HoughValidUpdate.cxx.
References MdcDetector::instance(), and m_gm.
00094 { 00095 00096 //Initailize MdcDetector 00097 m_gm = MdcDetector::instance(0); 00098 if(NULL == m_gm) return StatusCode::FAILURE; 00099 00100 return StatusCode::SUCCESS; 00101 }
double HoughValidUpdate::CFtrans | ( | double | x, | |
double | y | |||
) | [private] |
void HoughValidUpdate::Combine | ( | TH2D * | h1, | |
int | widthX, | |||
int | widthY, | |||
vector< vector< vector< int > > > | vec_selectNum, | |||
M1 & | hit_combine, | |||
M2 & | peak_combine_num, | |||
int & | ntrack, | |||
vector< int > | peakList1[] | |||
) | [private] |
Definition at line 2813 of file HoughValidUpdate.cxx.
References binX, binY, Combine_two_cell(), genRecEmupikp::i, ganga-rec::j, m_debug, and m_hit_pro.
Referenced by execute().
02813 { 02814 int m=0; 02815 int n=0; 02816 int addnum=999; 02817 int npeak1=peak_track; 02818 vector<bool> peaktrack(npeak1,true); 02819 while(addnum!=0) 02820 { 02821 // cout<<"Ntrack: "<<n<<endl; 02822 addnum=0; 02823 // int peakXtemp=peakList[0][0]-1; 02824 // int peakYtemp=peakList[1][0]-1; 02825 // cout<<"vec_selectpeak: ("<<peakXtemp<<","<<peakYtemp<<"): "<<vec_selectNum[peakXtemp][peakYtemp].size()<<endl; 02826 // for(int i =0;i<vec_selectNum[peakXtemp][peakYtemp].size();i++){cout<<" iN:vec: "<<vec_selectNum[peakXtemp][peakYtemp][i]<<endl;} 02827 02828 int peakX=peakList[0][n]; 02829 int peakY=peakList[1][n]; 02830 for (int px=peakX-widthX; px<=peakX+widthX; px++) { 02831 for (int py=peakY-widthY; py<=peakY+widthY; py++) { 02832 int ax; 02833 if (px<0) { ax=px+binX; } 02834 else if (px>=binX) { ax=px-binX; } 02835 else { ax=px; } 02836 if ( (ax!=peakX || py!=peakY) && py>=0 && py<binY) Combine_two_cell(h1,peakX,peakY,ax,py,vec_selectNum,m,peak_hit_list,peak_hit_num); 02837 } 02838 } 02839 // for(int i=0;i<peak_hit_list[m].size();i++){ 02840 // cout<<"COMBINE: "<<peak_hit_list[m][i]<<endl; 02841 // } 02842 02843 for(int i=n+1;i<npeak1;i++){ 02844 if(peaktrack[i]==false) continue; 02845 int peaktheta=peakList[0][i]; 02846 int peakrho=peakList[1][i]; 02847 int peakhitNum=peakList[2][i]; 02848 int count_same_hit=0; 02849 for(int j=0;j<peakhitNum;j++){ 02850 for(int k=0;k<peak_hit_num[m];k++){ 02851 if(vec_selectNum[peaktheta][peakrho][j]==peak_hit_list[m][k]){ 02852 count_same_hit++; 02853 } 02854 } 02855 } 02856 if(m_debug>1) cout<<" pro: "<<"peak: "<<i<<" "<<(double)count_same_hit/peakhitNum<<endl; 02857 double f_hit=m_hit_pro; 02858 if(count_same_hit>f_hit*peakhitNum){ 02859 peaktrack[i]=false; 02860 // cout<<" track_candi: "<<i<<" is false"<<endl; 02861 } 02862 } 02863 for(int i=n+1;i<npeak1;i++){ 02864 if(peaktrack[i]==true){ 02865 addnum=i-n; 02866 break; 02867 } 02868 } 02869 if(addnum!=0) m++; 02870 n=n+addnum; 02871 } 02872 int ntrack=0; 02873 for(int i=0;i<npeak1;i++){ 02874 //cout<<"track"<<i<<": "<<"("<<peakList[0][i]<<","<<peakList[1][i]<<","<<peakList[2][i]<<")"<<" "<<"truth: "<<peaktrack[i]<<endl; 02875 if(peaktrack[i]==true){ 02876 ntrack++; 02877 } 02878 } 02879 peak_track=ntrack; 02880 }
void HoughValidUpdate::Combine_two_cell | ( | TH2D * | h1, | |
int | xPeakCell, | |||
int | yPeakCell, | |||
int | xAround, | |||
int | yaround, | |||
vector< vector< vector< int > > > | vec_selectNum, | |||
int | m, | |||
M1 & | hit_combine, | |||
M2 & | peak_combine_num | |||
) | [private] |
Definition at line 2881 of file HoughValidUpdate.cxx.
References genRecEmupikp::i, ganga-rec::j, and delete_small_size::size.
Referenced by Combine().
02881 { 02882 int size_of_around=vec_selectNum[xAround][yAround].size(); 02883 if(size_of_around==0) return; 02884 // cout<<"size of around : "<<size_of_around<<endl; 02885 //for(int iaround=0;iaround<size_of_around;iaround++){ 02886 // cout<<"around: "<<vec_selectNum[xAround][yAround][iaround]<<endl; 02887 //} 02888 int size_of_peak=vec_selectNum[xPeakCell][yPeakCell].size(); 02889 // cout<<"size of peak: "<<size_of_peak<<endl; 02890 // for(int ipeak=0;ipeak<size_of_peak;ipeak++){ 02891 // cout<<"peak: "<<vec_selectNum[xPeakCell][yPeakCell][ipeak]<<endl; 02892 // } 02893 if(peak_hit_list[m].size()==0){ 02894 for(int inum=0;inum<size_of_peak;inum++){ 02895 peak_hit_list[m][inum]=vec_selectNum[xPeakCell][yPeakCell][inum]; 02896 } 02897 } 02898 vector<int> vec_more_select; 02899 for (int i =0;i<size_of_around;i++){ 02900 int same_vec_with_hitcombin=0; 02901 for(int j= 0;j<peak_hit_list[m].size();j++){ 02902 if(vec_selectNum[xAround][yAround][i]==peak_hit_list[m][j]){ 02903 same_vec_with_hitcombin=1; 02904 break; 02905 } 02906 } 02907 if(same_vec_with_hitcombin==0) { 02908 // cout<<" peak_hit_list push : "<<vec_selectNum[xAround][yAround][i]<<endl; 02909 vec_more_select.push_back(vec_selectNum[xAround][yAround][i]); 02910 } 02911 } 02912 int peak_hit_list_size=peak_hit_list[m].size(); 02913 for(int i=0;i<vec_more_select.size();i++ ){ 02914 peak_hit_list[m][i+peak_hit_list_size]=vec_more_select[i]; 02915 } 02916 int peak_hit_list_size_new=peak_hit_list[m].size(); 02917 peak_hit_num[m]=peak_hit_list_size_new; 02918 // cout<<"peak_combine: "<<peak_hit_num[m]<<endl; 02919 for(int i= 0;i<peak_hit_list_size_new;i++){ 02920 // cout<<" peak_hit_list: "<<peak_hit_list[m][i]<<endl; 02921 } 02922 02923 02924 02925 // for(int i =0;i<size_of_around;i++){ 02926 // int count_same_with_peak=0; 02927 // for(int j =0;j<vec_selectNum[xPeakCell][yPeakCell].size();j++){ 02928 // if(vec_selectNum[xAround][yAround][i]==vec_selectNum[xPeakCell][yPeakCell][j]) { 02929 // count_same_with_peak=1; 02930 // break; 02931 // } 02932 // } 02933 // if(count_same_with_peak==0) { 02934 // // cout<<" push : "<<vec_selectNum[xAround][yAround][i]<<endl; 02935 // vec_selectNum[xPeakCell][yPeakCell].push_back(vec_selectNum[xAround][yAround][i]); 02936 // } 02937 // } 02938 // 02939 // int size_of_peak_update=vec_selectNum[xPeakCell][yPeakCell].size(); 02940 }
int HoughValidUpdate::digiToHots | ( | MdcDigiVec | mdcDigiVec, | |
TrkRecoTrk * | newTrack, | |||
vector< bool > | vec_truthHit, | |||
uint32_t | getDigiFlag, | |||
vector< double > | vec_flt_truth, | |||
map< int, double > | hit_to_circle1, | |||
vector< MdcHit * > & | ||||
) | [private] |
Definition at line 1849 of file HoughValidUpdate.cxx.
References TrkHitList::appendHot(), MdcHit::driftDist(), MdcHit::driftTime(), RawData::getChargeChannel(), MdcCalibFunSvc::getT0(), MdcCalibFunSvc::getTimeWalk(), TrkRecoTrk::hits(), RawData::identify(), iter(), MdcDetector::Layer(), MdcID::layer(), m_bunchT0, m_debug, m_gm, m_mdcCalibFunSvc, TrkHitList::nHit(), MdcHit::rawTime(), MdcLayer::rMid(), MdcHit::setCalibSvc(), MdcHit::setCountPropTime(), TrkHitOnTrk::setFltLen(), subSeperate::temp, and MdcID::wire().
Referenced by execute().
01849 { 01850 TrkHitList* trkHitList = newTrack->hits(); 01851 int digiId=0; 01852 if(m_debug>0) cout<<"MdcDigiVec(in 2d fit): "<<mdcDigiVec.size()<<endl; 01853 MdcDigiVec::iterator iter = mdcDigiVec.begin(); 01854 for (;iter != mdcDigiVec.end(); iter++, digiId++) { 01855 if(vec_truthHit[digiId]==false) continue; 01856 const MdcDigi* aDigi = *iter; 01857 MdcHit *hit = new MdcHit(aDigi, m_gm); 01858 vec_for_clean.push_back(hit); 01859 // MdcHit *hit=MdcHit(aDigi, m_gm); 01860 Identifier id = aDigi->identify(); 01861 int layer = MdcID::layer(id); 01862 int wire = MdcID::wire(id); 01863 int hittemp=layer*1000+wire; 01864 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;} 01865 hit->setCalibSvc(m_mdcCalibFunSvc); 01866 hit->setCountPropTime(true); 01867 // new fltLen and ambig 01868 int newAmbig = 0; 01869 // calc. position of this point 01870 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here 01871 MdcHitOnTrack* newhot = &temp; 01872 double fltLen = m_gm->Layer(layer)->rMid(); 01873 // newhot->setFltLen(fltLen); 01874 newhot->setFltLen(vec_flt_least[digiId]); //caculate by mc 01875 //newhot->setFltLen(0); 01876 //newhot->setHitRms(0.02); 01877 //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl; 01878 double ddCut; 01879 //if(layer<8) ddCut=0.6; 01880 //else ddCut=0.8; 01881 if(layer<8) ddCut=1.0; 01882 else ddCut=1.0; 01883 int use_in2d=1; 01884 if(hit->driftDist(m_bunchT0,0)>ddCut||vec_truthHit[digiId]==false/*||vec_theta_ontrack[digiId]>M_PI*/||fabs(hit_to_circle1[hittemp])>10.) { use_in2d=0; } 01885 if(m_debug>1) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt "<<hit->driftTime(m_bunchT0,0)<<") T0 " << m_mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<< m_mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<<" m_bunchT0 "<<m_bunchT0<<" hit to circle: "<<hit_to_circle1[hittemp]<<" dist: "<<hit->driftDist(m_bunchT0,0)<<" use?: "<<use_in2d<<endl; 01886 //if(hit->driftTime(m_bunchT0,0)>800) continue; 01887 if(use_in2d==0) continue; 01888 trkHitList->appendHot(newhot); 01889 //delete hit; 01890 } 01891 if(m_debug>0) std::cout<<std::endl; 01892 return trkHitList->nHit(); 01893 }
int HoughValidUpdate::digiToHots2 | ( | MdcDigiVec | mdcDigiVec, | |
TrkRecoTrk * | newTrack, | |||
vector< bool > | vec_truthHit, | |||
vector< double > | vec_flt_truth, | |||
uint32_t | getDigiFlag, | |||
vector< int > | vec_slant, | |||
vector< double > | vec_l_Calcu, | |||
vector< double > | vec_l_axial, | |||
vector< double > | vec_theta_ontrack, | |||
vector< int > & | vec_hitbeforefit, | |||
map< int, double > | hit_to_circle2, | |||
vector< MdcHit * > & | , | |||
double | ||||
) | [private] |
Definition at line 2233 of file HoughValidUpdate.cxx.
References TrkHitList::appendHot(), TrkHitList::begin(), MdcHit::driftDist(), MdcHit::driftTime(), RawData::getChargeChannel(), MdcCalibFunSvc::getT0(), MdcCalibFunSvc::getTimeWalk(), TrkRecoTrk::hits(), RawData::identify(), MdcDetector::Layer(), MdcID::layer(), m_bunchT0, m_debug, m_gm, m_mdcCalibFunSvc, M_PI, TrkHitList::nHit(), MdcHit::rawTime(), MdcLayer::rMid(), MdcHit::setCalibSvc(), MdcHit::setCountPropTime(), TrkHitOnTrk::setFltLen(), subSeperate::temp, and MdcID::wire().
Referenced by execute().
02233 { 02234 TrkHitList* trkHitList = newTrack->hits(); 02235 //MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag); 02236 if(m_debug>0) cout<<"MdcDigiVec(in 3d fit): "<<mdcDigiVec.size()<<endl; 02237 MdcDigiVec::iterator iter1 = mdcDigiVec.begin(); 02238 int digiId=0; 02239 int nhit_beforefit_temp=0; 02240 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) { 02241 //cout<<"fltcut: "<<m_fltcut<<endl; 02242 if(vec_theta_ontrack[digiId]>M_PI) continue; 02243 // if(digiId>20&&vec_slant[digiId]!=0) continue; 02244 const MdcDigi* aDigi = *iter1; 02245 MdcHit *hit = new MdcHit(aDigi, m_gm); 02246 vec_for_clean.push_back(hit); 02247 //MdcHit *hit=MdcHit(aDigi, m_gm); 02248 Identifier id = aDigi->identify(); 02249 int layer = MdcID::layer(id); 02250 int wire = MdcID::wire(id); 02251 int hittemp=layer*1000+wire; 02252 hit->setCalibSvc(m_mdcCalibFunSvc); 02253 hit->setCountPropTime(true); //hit->setCountPropTime(m_countPropTime); 02254 02255 // new fltLen and ambig 02256 int newAmbig = 0; 02257 // calc. position of this point 02258 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here 02259 MdcHitOnTrack* newhot = &temp; 02260 double fltLen = m_gm->Layer(layer)->rMid(); 02261 // newhot->setFltLen(fltLen); 02262 newhot->setFltLen(vec_flt[digiId]); 02263 //if (vec_slant[digiId]==0) {newhot->setFltLen(vec_flt[digiId]);/* cout<<" ("<<layer<<","<<wire<<") "<<"l_truth: "<<vec_flt_truth[digiId]<<" l_axil: "<<vec_flt[digiId]<<endl;*/} //caculate by mc 02264 // else {newhot->setFltLen(vec_l_Calcu[digiId]);/* cout<<" ("<<layer<<","<<wire<<") "<<" l_truth: "<<vec_flt_truth[digiId]<<" l_Calcu: "<<vec_l_Calcu[digiId]<<endl;*/} //caculate by mc 02265 // newhot->setActivity(1); 02266 //bool active = newhot->isActive(); 02267 //newhot->setHitRms(0.02); 02268 //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl; 02269 double ddCut; 02270 double distcirCut; 02271 //if(layer<8) {ddCut=0.8; distcirCut=5;} 02272 //else {ddCut=0.6; distcirCut=2;} 02273 if(layer<8) {ddCut=1.0; distcirCut=5;} 02274 else {ddCut=1.0; distcirCut=2;} 02275 int use_in3d=1; 02276 if(hit->driftDist(m_bunchT0,0)>ddCut||vec_truthHit[digiId]==false||/*vec_theta_ontrack[digiId]>M_PI||*/fabs(hit_to_circle2[hittemp])>distcirCut) { use_in3d=0; } 02277 if(m_debug>2) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt "<<hit->driftTime(m_bunchT0,0)<<") T0 " << m_mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<< m_mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<<" vec_truth: "<<vec_truthHit[digiId]<<" theta: "<<vec_theta_ontrack[digiId]<<" ddcut?: "<<hit->driftDist(m_bunchT0,0)<<" distocir: "<<hit_to_circle2[hittemp]<<" use:? "<<use_in3d<<endl; 02278 if(use_in3d==0) continue; 02279 02280 vec_hitbeforefit.push_back(hittemp); 02281 trkHitList->appendHot(newhot); 02282 } 02283 if(m_debug>0) std::cout<<std::endl; 02284 return trkHitList->nHit(); 02285 }
StatusCode HoughValidUpdate::execute | ( | ) |
Definition at line 317 of file HoughValidUpdate.cxx.
References a0, abs, AmbigChoose(), MdcRawDataProvider::b_dropHot, MdcRawDataProvider::b_keepBadTdc, MdcRawDataProvider::b_keepUnmatch, MdcGeoWire::Backward(), TrkHotList::begin(), binX, binY, CFtrans(), HoughValidUpdate::Mytrack::charge, TrkAbsFit::chisq(), Combine(), compare(), cos(), count, TrkExchangePar::d0(), d0_mc, digiToHots(), digiToHots2(), dz_mc, deljobs::end, showlog::err, TrkErrCode::failure(), Bes_Common::FATAL, FillHist(), TrkHitList::fit(), TrkRecoTrk::fitResult(), MdcGeoWire::Forward(), GetMcInfo(), RawDataProviderSvc::getMdcDigiVec(), TrkFit::helix(), TrkRecoTrk::hits(), TrkHitList::hotList(), TrkRecoTrk::hots(), genRecEmupikp::i, Bes_Common::INFO, ganga-rec::j, MdcID::layer(), LeastSquare(), m_amb, m_ambig_no_match, m_ambig_no_match_propotion, m_axiallose, m_binx, m_biny, m_bunchT0, m_cell, m_chis_pt, m_combineTracking, m_context, m_d0_final, m_d0_tds, m_data, m_debug, m_delta_l, m_delta_z, m_digi_truth, m_drawTuple, m_dropHot, m_dropTrkChi2Cut, m_dropTrkChi2NdfCut, m_dropTrkDrCut, m_dropTrkDzCut, m_dropTrkNhitCut, m_dropTrkPtCut, m_drTruth, m_dzTruth, m_e_delta, m_eventcut, m_eventNum, m_fillTruth, m_gm, m_hitOnTrk, m_keepBadTdc, m_keepUnmatch, m_layer, m_layer_max, m_layerNhit, m_lCalcu, m_lCalcuLeft, m_lCalcuRight, m_ltruth, m_maxrho, m_mcpar, m_mdcCalibFunSvc, m_mdcGeomSvc, m_mdcPrintSvc, m_method, m_ndev1, m_ndev2, m_nHit, m_nHit_Mc, m_nHitAxial, m_nHitAxial_Mc, m_nHitAxialSelect, m_nHitDigi, m_nHitSelect, m_nHitStereo, m_nHitStereo_use, m_nHitStereoSelect, m_npeak_1, m_npeak_2, m_npeak_after_tracking, m_omega_final, m_omega_tds, m_omegaTruth, m_outputTrk, m_p, m_p_rec_final, m_p_rec_tds, m_peakCellNumX, m_peakCellNumY, m_phi0_final, m_phi0_tds, m_phi0Truth, M_PI, m_pickHits, m_postruth, m_pt, m_pt2_rec_final, m_pt2_rec_tds, m_Q, m_q, m_Q_tds, m_qualityFactor, m_r, m_rawDataProviderSvc, m_removeBadTrack, m_rhocut, m_runNum, m_setpeak_fit, m_slant, m_stereolose, m_tanl_Cal, m_tanl_final, m_tanl_mc, m_tanl_tds, m_trackFailure, m_trackNum, m_trackNum_fit, m_trackNum_tds, m_u, m_useMcInfo, m_v, m_x, m_x_circle, m_x_east, m_x_west, m_y, m_y_circle, m_y_east, m_y_west, m_z, m_z0_Cal, m_z0_final, m_z0_tds, m_z_east, m_z_west, m_zCalcu, m_zCalcuLeft, m_zCalcuRight, m_ztruth, TrkSimpleMaker< T >::makeTrack(), msgSvc(), TrkHotList::nActive(), TrkFit::nActive(), HoughValidUpdate::Mytrack::newTrack, TrkHitList::nHit(), ntuplehit, TrkExchangePar::omega(), omega_mc, HoughValidUpdate::Mytrack::phi, TrkExchangePar::phi0(), phi0, phi0_mc, TrkRecoTrk::print(), TrkRecoTrk::printAll(), MdcPrintSvc::printRecMdcTrackCol(), HoughValidUpdate::Mytrack::pt, push_back(), HoughValidUpdate::Mytrack::r, EventModel::Recon::RecMdcHitCol, EventModel::Recon::RecMdcTrackCol, MdcHit::setCalibSvc(), MdcHit::setCountPropTime(), TrkSimpleMaker< T >::setFlipAndDrop(), sin(), delete_small_size::size, SlantId(), MdcTrack::storeTrack(), swap, t_eventNum, t_runNum, TrkExchangePar::tanDip(), tanl_mc, MdcTrack::track(), track_fit, HoughValidUpdate::Mytrack::use_in_tds, v, HoughValidUpdate::Mytrack::vec_hit_on_trk, Bes_Common::WARNING, MdcGeomSvc::Wire(), MdcID::wire(), x, HoughValidUpdate::Mytrack::x_cir, HoughValidUpdate::Mytrack::y_cir, TrkExchangePar::z0(), and zPosition().
00317 { 00318 00319 cout.precision(6); 00320 MsgStream log(msgSvc(), name()); 00321 log << MSG::INFO << "in execute()" << endreq; 00322 00323 //event start time 00324 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol"); 00325 if (!evTimeCol) { 00326 log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq; 00327 m_bunchT0=0.; 00328 }else{ 00329 RecEsTimeCol::iterator iter_evt = evTimeCol->begin(); 00330 if (iter_evt != evTimeCol->end()){ 00331 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns 00332 } 00333 } 00334 00335 // Get the event header, print out event and run number 00336 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00337 if (!eventHeader) { 00338 log << MSG::FATAL << "Could not find Event Header" << endreq; 00339 return( StatusCode::FAILURE); 00340 } 00341 00342 DataObject *aTrackCol; 00343 DataObject *aRecHitCol; 00344 SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 00345 if(!m_combineTracking){ 00346 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol); 00347 if(aTrackCol != NULL) { 00348 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol"); 00349 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol"); 00350 } 00351 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol); 00352 if(aRecHitCol != NULL) { 00353 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol"); 00354 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol"); 00355 } 00356 } 00357 00358 RecMdcTrackCol* trackList; 00359 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol); 00360 if (aTrackCol) { 00361 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol); 00362 }else{ 00363 trackList = new RecMdcTrackCol; 00364 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList); 00365 if(!sc.isSuccess()) { 00366 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq; 00367 return StatusCode::FAILURE; 00368 } 00369 } 00370 00371 RecMdcHitCol* hitList; 00372 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol); 00373 if (aRecHitCol) { 00374 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol); 00375 }else{ 00376 hitList = new RecMdcHitCol; 00377 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList); 00378 if(!sc.isSuccess()) { 00379 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq; 00380 return StatusCode::FAILURE; 00381 } 00382 } 00383 00384 //remove bad quality or low pt track(s) 00385 if(m_removeBadTrack) { 00386 vector<RecMdcTrack*> vec_trackList; 00387 if( m_debug>0 ) cout<<"PATTSF collect "<<trackList->size()<<" track. "<<endl; 00388 if(trackList->size()!=0){ 00389 RecMdcTrackCol::iterator iter_pat = trackList->begin(); 00390 for(;iter_pat!=trackList->end();iter_pat++){ 00391 double d0=(*iter_pat)->helix(0); 00392 double kap = (*iter_pat)->helix(2); 00393 double pt = 0.00001; 00394 if(fabs(kap)>0.00001) pt = fabs(1./kap); 00395 double dz=(*iter_pat)->helix(4); 00396 double chi2=(*iter_pat)->chi2(); 00397 if( m_debug>0) cout<<"d0: "<<d0<<" z0: "<<dz<<" pt "<<pt<<" chi2 "<<chi2<<endl; 00398 if(fabs(d0)>m_dropTrkDrCut || fabs(dz)>m_dropTrkDzCut || chi2>m_dropTrkChi2Cut || pt>m_dropTrkPtCut) { 00399 vec_trackList.push_back(*iter_pat); 00400 //remove hits on track 00401 HitRefVec rechits = (*iter_pat)->getVecHits(); 00402 HitRefVec::iterator hotIter = rechits.begin(); 00403 while (hotIter!=rechits.end()) { 00404 if(m_debug>0) cout <<"remove hit " << (*hotIter)->getId() <<endl; 00405 hitList->remove(*hotIter); 00406 hotIter++; 00407 } 00408 trackList->remove(*iter_pat); 00409 iter_pat--; 00410 } 00411 // if( m_debug>0 ) cout<<"after delete trackcol size : "<<trackList->size()<<endl; 00412 } 00413 } else { 00414 if(m_debug>0) cout<<" PATTSF find 0 track. "<<endl; 00415 } 00416 }//end of remove bad quality high pt track 00417 00418 int nTdsTk=trackList->size(); 00419 00420 00421 t_eventNum=eventHeader->eventNumber(); 00422 t_runNum=eventHeader->runNumber(); 00423 if( m_drawTuple ){ 00424 m_eventNum=t_eventNum; 00425 m_runNum=t_runNum; 00426 } 00427 if(m_debug>0) cout<<" begin event :"<<t_eventNum<<endl; 00428 if(t_eventNum<=m_eventcut) { 00429 cout<<" eventNum cut "<<t_eventNum<<endl; 00430 return StatusCode::SUCCESS; 00431 } 00432 log << MSG::INFO << "HoughValidUpdate: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq; 00433 00434 vector<double> vec_u; 00435 vector<double> vec_v; 00436 vector<double> vec_p; 00437 vector<double> vec_x_east; 00438 vector<double> vec_y_east; 00439 vector<double> vec_z_east; 00440 vector<double> vec_x_west; 00441 vector<double> vec_y_west; 00442 vector<double> vec_z_west; 00443 vector<double> vec_z_Mc; 00444 vector<double> vec_x; 00445 vector<double> vec_y; 00446 vector<double> vec_z; 00447 vector<int> vec_layer; 00448 vector<int> vec_wire; 00449 vector<int> vec_slant; 00450 vector<double> vec_zStereo; 00451 vector<double> l; 00452 00453 vector<int> vec_layer_Mc; 00454 vector<int> vec_wire_Mc; 00455 vector<int> vec_hit_Mc; 00456 vector<double> vec_ztruth_Mc; 00457 vector<double> vec_flt_truth_mc; 00458 vector<double> vec_pos_flag; 00459 vector<double> vec_phit_MC; 00460 00461 // get mc particle digi 00462 int mc_particle=GetMcInfo(); 00463 if(m_debug>2) cout<<"MC particle : "<<mc_particle<<endl; 00464 00465 // retrieve Mdc truth hits 00466 if(m_fillTruth ==1 && m_useMcInfo) { 00467 int digiId_Mc=0; 00468 int nHitAxial_Mc=0; 00469 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol"); 00470 if (!mcMdcMcHitCol) { 00471 log << MSG::INFO << "Could not find MdcMcHit" << endreq; 00472 }else{ 00473 Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin(); 00474 00475 for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) { 00476 Identifier mdcid = (*iterMcHit)->identify(); 00477 int layerId_Mc = MdcID::layer(mdcid); 00478 int wireId_Mc = MdcID::wire(mdcid); 00479 double hittemp=layerId_Mc*1000+wireId_Mc; 00480 double mc_hit_distance =(*iterMcHit)->getDriftDistance(); 00481 double z_mc = (*iterMcHit)->getPositionZ()/10.;//mm 2 cm 00482 double flt_truth=(z_mc-dz_mc)/tanl_mc; 00483 double pos_flag=(*iterMcHit)->getPositionFlag(); 00484 00485 // cout<<"flag: "<<(*iterMcHit)->getPositionFlag()<<endl; 00486 if( (*iterMcHit)->getPositionFlag()==-999 ){ 00487 // cout<<"track Id: "<<(*iterMcHit)->getTrackIndex()<<endl; 00488 double px= (*iterMcHit)->getPositionX(); 00489 double py= (*iterMcHit)->getPositionY(); 00490 double pz= (*iterMcHit)->getPositionZ(); 00491 double pxyz = sqrt(px*px+py*py+pz*pz); 00492 vec_phit_MC.push_back(pxyz); 00493 // cout<<"PX: "<<px<<endl; 00494 // cout<<"PY: "<<py<<endl; 00495 // cout<<"PT: "<<sqrt(px*px+py*py)<<endl; 00496 // cout<<"PZ: "<<pz<<endl; 00497 // cout<<"PXYZ: "<<pxyz<<endl; 00498 // cout<<endl; 00499 } 00500 00501 // cout<<"( MC: " <<layerId_Mc<<","<<wireId_Mc<<"):"<<endl; 00502 vec_layer_Mc.push_back(layerId_Mc); 00503 vec_wire_Mc.push_back(wireId_Mc); 00504 vec_hit_Mc.push_back(hittemp); 00505 vec_ztruth_Mc.push_back(z_mc); 00506 vec_flt_truth_mc.push_back(flt_truth); 00507 vec_pos_flag.push_back(pos_flag); 00508 00509 if(m_debug>5) { 00510 cout<<"(" <<layerId_Mc<<","<<wireId_Mc<<"):"<<endl; 00511 cout<<" hit_distance : "<<mc_hit_distance<<endl; 00512 cout<<" position flag : "<<pos_flag<<endl; 00513 cout<<" hit_z_mc : "<<z_mc<<endl; 00514 cout<<" flt_truth : "<<flt_truth<<endl; 00515 } 00516 00517 if(m_data==0){ 00518 if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) { 00519 nHitAxial_Mc++;} 00520 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc); 00521 HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm 00522 HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm 00523 00524 double x = (*iterMcHit)->getPositionX()/10.;//mm 2 cm 00525 double y = (*iterMcHit)->getPositionY()/10.;//mm 2 cm 00526 double z = (*iterMcHit)->getPositionZ()/10.;//mm 2 cm 00527 double u=CFtrans(x,y); 00528 double v=CFtrans(y,x); 00529 double p=sqrt(u*u+v*v); 00530 00531 vec_x_east.push_back(eastP.x()); 00532 vec_y_east.push_back(eastP.y()); 00533 vec_z_east.push_back(eastP.z()); 00534 vec_x_west.push_back(westP.x()); 00535 vec_y_west.push_back(westP.y()); 00536 vec_z_west.push_back(westP.z()); 00537 vec_x.push_back(x); 00538 vec_y.push_back(y); 00539 vec_z.push_back(z); 00540 vec_u.push_back(u); 00541 vec_v.push_back(v); 00542 vec_p.push_back(p); 00543 vec_slant.push_back( SlantId(layerId_Mc) ); 00544 00545 m_x_east[digiId_Mc]=eastP.x(); 00546 m_y_east[digiId_Mc]=eastP.y(); 00547 m_z_east[digiId_Mc]=eastP.z(); 00548 m_x_west[digiId_Mc]=westP.x(); 00549 m_y_west[digiId_Mc]=westP.y(); 00550 m_z_west[digiId_Mc]=westP.z(); 00551 m_layer[digiId_Mc]=layerId_Mc; 00552 m_cell[digiId_Mc]=wireId_Mc; 00553 m_x[digiId_Mc]=x; 00554 m_y[digiId_Mc]=y; 00555 m_z[digiId_Mc]=z; 00556 m_u[digiId_Mc]=u; 00557 m_v[digiId_Mc]=v; 00558 m_p[digiId_Mc]=p; 00559 m_slant[digiId_Mc]=SlantId(layerId_Mc); 00560 } 00561 } 00562 } 00563 m_nHit_Mc=digiId_Mc; 00564 m_nHitAxial_Mc=nHitAxial_Mc; 00565 if(0==m_data) {m_nHit=digiId_Mc; m_nHitAxial=nHitAxial_Mc;} 00566 } 00567 00568 // retrieve Mdc digi vector form RawDataProviderSvc 00569 vector<int> vec_hit; 00570 vector<int> vec_track_index; 00571 set<int> hit_noise; //1 00572 uint32_t getDigiFlag = 0; 00573 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot; 00574 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc; 00575 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch; 00576 00577 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag); 00578 if(m_debug>0) cout<<"MdcDigiVec: "<<mdcDigiVec.size()<<endl; 00579 MdcDigiVec::iterator iter1 = mdcDigiVec.begin(); 00580 int digiId = 0; 00581 int nHitAxial = 0; 00582 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) { 00583 Identifier mdcId= (*iter1)->identify(); 00584 int layerId = MdcID::layer(mdcId); 00585 int wireId = MdcID::wire(mdcId); 00586 double hittemp=layerId*1000+wireId; 00587 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId); 00588 HepPoint3D eastP = geowir->Backward()/10.0;// cm 00589 HepPoint3D westP = geowir->Forward() /10.0;// cm 00590 if ((layerId>=8&&layerId<=19)||(layerId>=36)) nHitAxial++; 00591 00592 int track_index=(*iter1)->getTrackIndex(); 00593 if( track_index>=1000 ) track_index-=1000; 00594 if( track_index<0 ) hit_noise.insert(hittemp); 00595 int tchannal=(*iter1)->getTimeChannel(); 00596 int qchannal=(*iter1)->getChargeChannel(); 00597 if( tchannal!=qchannal && m_debug>0 ) cout<<"("<<layerId<<","<<wireId<<")"<< 0<<endl; 00598 if( m_debug>3 ) cout<<"track_index in digi: "<<"("<<layerId<<","<<wireId<<" "<< track_index<<endl; 00599 if( m_debug>3 ) cout<<"event: "<<t_eventNum<<" "<<"layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"zeast: "<<eastP.z()<<" "<<"zwest: "<<westP.z()<<" "<<endl; 00600 00601 //conformal trans: 00602 double x =(eastP.x()+westP.x())/2.; 00603 double y =(eastP.y()+westP.y())/2.; 00604 double u=CFtrans(x,y); 00605 double v=CFtrans(y,x); 00606 if(m_debug>3) cout<<"digiId: "<<digiId<<" layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"x: "<<x<<" "<<"y: "<<y<<" "<< "u: "<<u<<"v: "<<v<<endl; 00607 00608 vec_track_index.push_back(track_index); 00609 vec_layer.push_back(layerId); 00610 vec_wire.push_back(wireId); 00611 vec_hit.push_back(hittemp); 00612 vec_x_east.push_back(eastP.x()); 00613 vec_y_east.push_back(eastP.y()); 00614 vec_z_east.push_back(eastP.z()); 00615 vec_x_west.push_back(westP.x()); 00616 vec_y_west.push_back(westP.y()); 00617 vec_z_west.push_back(westP.z()); 00618 vec_x.push_back(x); 00619 vec_y.push_back(y); 00620 vec_p.push_back(sqrt(u*u+v*v)); 00621 vec_u.push_back(u); 00622 vec_v.push_back(v); 00623 vec_slant.push_back( SlantId(layerId) ); 00624 00625 if(m_drawTuple){ 00626 m_x_east[digiId]=eastP.x(); 00627 m_y_east[digiId]=eastP.y(); 00628 m_z_east[digiId]=eastP.z(); 00629 m_x_west[digiId]=westP.x(); 00630 m_y_west[digiId]=westP.y(); 00631 m_z_west[digiId]=westP.z(); 00632 m_layer[digiId]=layerId; 00633 m_cell[digiId]=wireId; 00634 m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2; 00635 m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2; 00636 m_u[digiId]=u; 00637 m_v[digiId]=v; 00638 m_p[digiId]=sqrt(u*u+v*v); 00639 m_slant[digiId]=SlantId( layerId ); 00640 m_layerNhit[layerId]++; 00641 } 00642 } 00643 if(m_drawTuple){ 00644 m_nHit=digiId; 00645 m_nHitAxial=nHitAxial; 00646 m_nHitStereo=m_nHit-m_nHitAxial; 00647 } 00648 if(m_debug>0) cout<<"Hit digi number: "<<digiId<<endl; 00649 if(digiId<5||nHitAxial<3) { 00650 if(m_drawTuple) m_trackFailure=1; 00651 if(m_drawTuple) ntuplehit->write(); 00652 if(m_debug>0) cout<<"not enough hit "<<endl; 00653 return StatusCode::SUCCESS; 00654 } 00655 00656 //correspond mc_fltLenth with digiHit 00657 vector<int> vec_digitruth(digiId,0); 00658 vector<double> vec_flt_truth(digiId,999.); 00659 vector<double> vec_ztruth(digiId,999.); 00660 vector<double> vec_posflag_truth(digiId,999.); 00661 if(m_useMcInfo){ 00662 int if_exit_delta_e=0; 00663 int nhit_digi=0; 00664 for(int iHit=0;iHit<digiId;iHit++){ 00665 vector<int>::iterator iter_ihit = find( vec_hit_Mc.begin(),vec_hit_Mc.end(),vec_hit[iHit] ); 00666 if( iter_ihit==vec_hit_Mc.end() ) { 00667 vec_digitruth[iHit]=0; 00668 if(m_drawTuple) m_digi_truth[iHit]=0; 00669 if_exit_delta_e=1; 00670 } 00671 else { 00672 vec_digitruth[iHit]=1; 00673 m_digi_truth[iHit]=1; 00674 nhit_digi++; 00675 } 00676 for(int iHit_Mc=0;iHit_Mc<vec_flt_truth_mc.size();iHit_Mc++){ 00677 if(vec_hit[iHit]==vec_hit_Mc[iHit_Mc]) { 00678 vec_flt_truth[iHit]=vec_flt_truth_mc[iHit_Mc]; 00679 vec_ztruth[iHit]=vec_ztruth_Mc[iHit_Mc]; 00680 vec_posflag_truth[iHit]=vec_pos_flag[iHit_Mc]; 00681 if(m_drawTuple){ 00682 m_ztruth[iHit]=vec_ztruth_Mc[iHit_Mc]; 00683 m_ltruth[iHit]=vec_flt_truth_mc[iHit_Mc]; 00684 m_postruth[iHit]=vec_pos_flag[iHit_Mc]; 00685 } 00686 } 00687 } 00688 if(m_debug>3) cout<<" hitPos: "<<"("<<vec_layer[iHit]<<","<<vec_wire[iHit]<<")"<<" flt_truth: "<<vec_flt_truth[iHit]<<" z_truth: "<<vec_ztruth[iHit]<<" pos_flag: "<<vec_posflag_truth[iHit]<<" digi_truth: "<<vec_digitruth[iHit]<<endl; 00689 } 00690 if( m_drawTuple ){ 00691 m_e_delta=if_exit_delta_e; 00692 m_nHitDigi=nhit_digi; 00693 } 00694 } 00695 00696 // calcu the last layer 00697 vector<int>::iterator laymax=max_element( vec_layer.begin(),vec_layer.end() ); 00698 if(m_drawTuple) m_layer_max=*laymax; 00699 00700 //------------------------- finish hit ------------------------------- 00701 int nhit; 00702 int nhitaxial; 00703 if(m_data==0 && m_useMcInfo && m_fillTruth) { 00704 nhit=m_nHit_Mc; 00705 nhitaxial=m_nHitAxial_Mc; 00706 } 00707 else{ 00708 nhit=digiId; 00709 nhitaxial=nHitAxial; 00710 } 00711 00712 binX=m_binx; 00713 binY=m_biny; 00714 TH2D *h1 = new TH2D("h1","h1",binX,0,M_PI,binY,-m_maxrho,m_maxrho); 00715 TH2D *h2 = new TH2D("h2","h2",binX,0,M_PI,binY,-m_maxrho,m_maxrho); 00716 00717 //method 0 00718 /* 00719 if(m_method==0){ 00720 vector<double> vec_rho; 00721 vector<double> vec_theta; 00722 vector< vector<int> > vec_hitNum(2,vector<int>()); //store 2 hits in a cross point 00723 vector<int> xybin; 00724 //RhoTheta(nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum,vec_digitruth); //calculate cross point 00725 RhoThetaAxial(vec_slant,nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum,vec_digitruth); //calculate cross point 00726 FillRhoTheta(h1,vec_theta,vec_rho); //fill histgram method by rhotheta 00727 if(m_nCross>125000) return StatusCode::SUCCESS; 00728 00729 //sum in x 00730 int xsum_bin; 00731 int xsum=0; 00732 TF1 *f =new TF1("f","gaus",0,3.1415926); 00733 TH1D *hx=new TH1D("hx","hx",binX,0,3.1415926); 00734 TH1D *hx_new=new TH1D("hx_new","hx_new",m_binx,0,3.1415926); 00735 for(int ibinx =0;ibinx<binX;ibinx++){ 00736 int t_xsum=0; 00737 for(int ibiny =0;ibiny<binY;ibiny++){ 00738 t_xsum+=h1->GetBinContent(ibinx+1,ibiny+1); 00739 } 00740 hx->SetBinContent(ibinx+1,t_xsum); 00741 if( xsum <= t_xsum ) { 00742 xsum=t_xsum; 00743 xsum_bin=ibinx+1; 00744 } 00745 } 00746 for(int ibinx =0;ibinx<binX;ibinx++){ 00747 int hx_content = hx->GetBinContent(ibinx+1); 00748 int binx_new=(ibinx+1)-(xsum_bin-(binX/2+1)); 00749 if(binx_new<1) { binx_new+=binX; } 00750 if(binx_new>binX) { binx_new-=binX; } 00751 hx_new->SetBinContent( binx_new , hx_content); 00752 } 00753 hx_new->Fit("f","QN"); 00754 double xsigma=f->GetParameter(2); 00755 delete hx; 00756 delete hx_new; 00757 00758 //sum in y axial 00759 TH1D *hy=new TH1D("hy","hy",binY,-0.3,0.3); 00760 or(int ibiny =0;ibiny<binY;ibiny++){ 00761 int t_ysum=0; 00762 for(int ibinx =0;ibinx<binX;ibinx++){ 00763 t_ysum+=h1->GetBinContent(ibinx+1,ibiny+1); 00764 } 00765 hy->SetBinContent(ibiny+1,t_ysum); 00766 } 00767 hy->Fit("f","QN"); 00768 double ysigma=f->GetParameter(2); 00769 delete hy; 00770 delete f; 00771 00772 for(int ibinx=0;ibinx<binX;ibinx++){ 00773 for(int ibiny=0;ibiny<binY;ibiny++){ 00774 int binNum=h1->GetBinContent(ibinx+1,ibiny+1); 00775 xybin.push_back(binNum); 00776 m_xybin[binY*ibinx+ibiny]=binNum; 00777 } 00778 } 00779 m_xybinNum=binX*binY; 00780 m_xsigma=xsigma; 00781 m_ysigma=ysigma; 00782 } 00783 */ 00784 00785 //method 1 00786 vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) ); 00787 M2 peak_hit_num; 00788 M1 peak_hit_list; 00789 int peak_track; 00790 if(m_method==1){ 00791 FillHist(h1,vec_u,vec_v,nhit,vec_selectNum); 00792 int max_count=0; 00793 int max_binx=0; 00794 int max_biny=0; 00795 for(int i=0;i<binX;i++){ 00796 for(int j=0;j<binY;j++){ 00797 int count=h1->GetBinContent(i+1,j+1); 00798 if(max_count<count) { 00799 max_count=count; 00800 max_binx=i+1; 00801 max_biny=j+1; 00802 } 00803 } 00804 } 00805 00806 if(vec_selectNum[max_binx-1][max_biny-1].size() != max_count ) cout<< "ERROR IN VEC!!! "<<endl; 00807 if(m_debug>4) { 00808 cout<<"maxBin: ["<<max_binx-1<<","<<max_biny-1<<"]: "<<vec_selectNum[max_binx-1][max_biny-1].size()<<endl; 00809 for(int i =0;i<vec_selectNum[max_binx-1][max_biny-1].size();i++) cout<<" maxBin select hit: "<<vec_selectNum[max_binx-1][max_biny-1][i]<<endl; 00810 } 00811 00812 if(m_debug>4) { 00813 for(int ibinx=0;ibinx<binX;ibinx++){ 00814 for(int ibiny=0;ibiny<binY;ibiny++){ 00815 int bincount=h1->GetBinContent(ibinx+1,ibiny+1); 00816 if(vec_selectNum[ibinx][ibiny].size() != bincount ) cout<< "ERROR IN VECTORSELECT filling "<<endl; 00817 if(vec_selectNum[ibinx][ibiny].size() == 0 ) continue; 00818 cout<<"bin:"<<"["<<ibinx<<","<<ibiny<<"]"<<" select:"<<vec_selectNum[ibinx][ibiny].size()<<endl; 00819 for(int i=0;i<vec_selectNum[ibinx][ibiny].size();i++){ 00820 int ihit=vec_selectNum[ibinx][ibiny][i]; 00821 cout<<" hit: "<<ihit<<" ("<<vec_layer[ihit]<<","<<vec_wire[ihit]<<")"<<endl; 00822 } 00823 } 00824 } 00825 } 00826 00827 //find peak 00828 // set threshold 00829 int count_use=0; 00830 double sum=0; 00831 double sum2=0; 00832 for(int i=0;i<binX;i++){ 00833 for(int j=0;j<binY;j++){ 00834 int count=h1->GetBinContent(i+1,j+1); 00835 if(count!=0){ 00836 count_use++; 00837 sum+=count; 00838 sum2+=count*count; 00839 } 00840 } 00841 } 00842 double f_n=m_ndev1; 00843 double aver=sum/count_use; 00844 double dev=sqrt(sum2/count_use-aver*aver); 00845 int min_counts=(int)(aver+f_n*dev); 00846 if (min_counts<m_ndev2) min_counts=m_ndev2; 00847 if(m_debug>2) cout<<"aver: "<<aver<<" "<<"dev: "<<dev<<" min: "<<min_counts<<endl; 00848 00849 vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) ); 00850 vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) ); 00851 for(int i=0;i<binX;i++){ 00852 for(int j=0;j<binY;j++){ 00853 hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1); 00854 } 00855 } 00856 00857 int npeak_1=0; 00858 for (int r=0; r<binY; r++) { 00859 double binHigh=2*m_maxrho/binY; 00860 double rho_peak=r*binHigh-m_maxrho; 00861 for (int a=0; a<binX; a++) { 00862 long max_around=0; 00863 for (int ar=a-1; ar<=a+1; ar++) { 00864 for (int rx=r-1; rx<=r+1; rx++) { 00865 int ax; 00866 if (ar<0) { ax=ar+binX; } 00867 else if (ar>=binX) { ax=ar-binX; } 00868 else { ax=ar; } 00869 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) { 00870 if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; } 00871 } 00872 } 00873 } 00874 00875 if (hough_trans_CS[a][r]<=min_counts||fabs(rho_peak)<m_rhocut) { hough_trans_CS_peak[a][r]=0; } 00876 else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; } 00877 else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; } 00878 else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; } 00879 if(hough_trans_CS_peak[a][r]!=0) { 00880 if(m_debug>2) cout<<" possible peak1:"<<"["<<a<<","<<r<<"]"<<": "<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<endl; 00881 npeak_1++; 00882 } 00883 } 00884 } 00885 00886 00887 //find double-point-peaks in parameter space 00888 int npeak_2=0; 00889 for (int r=0; r<binY; r++) { 00890 for (int a=0; a<binX; a++) { 00891 if (hough_trans_CS_peak[a][r]==1) { 00892 hough_trans_CS_peak[a][r]=2; 00893 for (int ar=a-1; ar<=a+1; ar++) { 00894 for (int rx=r-1; rx<=r+1; rx++) { 00895 int ax; 00896 if (ar<0) { ax=ar+binX; } 00897 else if (ar>=binX) { ax=ar-binX; } 00898 else { ax=ar; } 00899 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) { 00900 if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; } 00901 } 00902 } 00903 } 00904 } 00905 if(hough_trans_CS_peak[a][r]==2) { 00906 double binHigh=2*m_maxrho/binY; 00907 double rho_peak=r*binHigh-m_maxrho; 00908 npeak_2++; 00909 if(m_debug>2){ 00910 cout<<" possible peak2: "<<"["<<a<<","<<r<<"]"<<": "<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<" rhopeak: "<<rho_peak<<endl; 00911 for(int i=0;i<hough_trans_CS[a][r];i++){ 00912 int hitNum=vec_selectNum[a][r][i]; 00913 if(m_debug>2) cout<<" select hit: "<<hitNum<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<endl; 00914 } 00915 } 00916 } 00917 } 00918 } 00919 00920 //fill the histogram again 00921 for(int i=0;i<binX;i++){ 00922 for(int j=0;j<binY;j++){ 00923 if(hough_trans_CS_peak[i][j]==2){ 00924 h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]); 00925 } 00926 } 00927 } 00928 00929 vector<int> peakList[3]; 00930 for(int n=0;n<npeak_2;n++){ 00931 for (int r=0; r<binY; r++) { 00932 for (int a=0; a<binX; a++) { 00933 if (hough_trans_CS_peak[a][r]==2) { 00934 peakList[0].push_back(a); 00935 peakList[1].push_back(r); 00936 peakList[2].push_back(hough_trans_CS[a][r]); 00937 } 00938 } 00939 } 00940 } 00941 00942 if(m_drawTuple){ 00943 m_npeak_1=npeak_1; 00944 m_npeak_2=npeak_2; 00945 } 00946 if(m_debug>0) { 00947 cout<<"npeak_1: "<<npeak_1<<endl; 00948 cout<<"npeak_2: "<<npeak_2<<endl; 00949 } 00950 00951 //sort peaks by height 00952 int n_max; 00953 for (int na=0; na<npeak_2-1; na++) { 00954 n_max=na; 00955 for (int nb=na+1; nb<npeak_2; nb++) { 00956 if (peakList[2][n_max]<peakList[2][nb]) { n_max=nb; } 00957 } 00958 if (n_max!=na) { // swap 00959 float swap[3]; 00960 for (int i=0; i<3; i++) { 00961 swap[i]=peakList[i][na]; 00962 peakList[i][na]=peakList[i][n_max]; 00963 peakList[i][n_max]=swap[i]; 00964 } 00965 } 00966 } 00967 00968 if(npeak_2==0||npeak_2>100){ 00969 if(m_debug>0) cout<<" FAILURE IN NPEAK2 !!!!!! "<<endl; 00970 delete h1; 00971 delete h2; 00972 if(m_drawTuple){ 00973 m_trackFailure=2; 00974 m_npeak_2=-999; 00975 } 00976 if(m_drawTuple) ntuplehit->write(); 00977 return StatusCode::SUCCESS; 00978 } 00979 00980 if(m_debug>2){ 00981 for(int n=0;n<npeak_2;n++){ 00982 int bintheta=peakList[0][n]; 00983 int binrho=peakList[1][n]; 00984 int count=peakList[2][n]; 00985 cout<<"possible peakList after SORT: "<<n<<": "<<"["<<bintheta<<","<<binrho<<"]: "<<count<<endl; 00986 for(int i=0;i<count;i++){ 00987 int hitNum=vec_selectNum[bintheta][binrho][i]; 00988 if(m_debug>2) cout<<" "<<" select hit:"<<":"<<hitNum<<" ------ "<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<endl; 00989 } 00990 } 00991 } 00992 00993 // combine peak to track 00994 peak_track=npeak_2; 00995 Combine(h1,m_peakCellNumX,m_peakCellNumY,vec_selectNum,peak_hit_list,peak_hit_num,peak_track,peakList); 00996 if(m_drawTuple) m_npeak_after_tracking=peak_track; 00997 00998 if(m_debug>0) cout<<"event"<<t_eventNum<<": "<<"has found: "<<peak_track<<" track."<<endl; 00999 for( int i=0;i<peak_track;i++){ 01000 if(m_debug>0) cout<<"peak["<<i<<"] collect "<<peak_hit_num[i]<<" hits "<<endl; 01001 int nhitaxialselect=0; // temp for each peak_track 01002 for(int j =0;j<peak_hit_num[i];j++){ 01003 int hit_number=peak_hit_list[i][j]; 01004 if(vec_slant[hit_number]==0) nhitaxialselect++ ; 01005 if(m_debug>0) cout<<" "<<" collect hits: ("<<vec_layer[hit_number]<<","<<vec_wire[hit_number]<<")"<<endl; 01006 } 01007 if(m_drawTuple){ 01008 m_nHitSelect[i]=peak_hit_num[i]; 01009 m_nHitAxialSelect[i]=nhitaxialselect; 01010 m_nHitStereoSelect[i]=peak_hit_num[i]-nhitaxialselect; 01011 m_axiallose[i]=m_nHitAxial-m_nHitAxialSelect[i]; 01012 m_stereolose[i]=m_nHit-m_nHitSelect[i]-m_axiallose[i]; 01013 } 01014 } 01015 } 01016 01017 vector<int> vec_hitbeforefit; 01018 vector<bool> vec_truthHit(nhit,false); 01019 int peak_fit=peak_track; 01020 if(m_setpeak_fit!=-1) { 01021 peak_fit=m_setpeak_fit; 01022 } 01023 if(m_drawTuple) m_trackNum=peak_fit; 01024 01025 vector<int> vec_hit_use; 01026 // vector<int> vec_hit_use_in2d; 01027 // vector<int> vec_hit_use_intrack1; 01028 // vector<int> vec_Qchoise; 01029 vector<TrkRecoTrk*> vec_newTrack; 01030 vector<TrkRecoTrk*> vec_track_in_tds; 01031 vector<int> vec_hitOnTrack; 01032 vector<MdcHit*> vec_for_clean; 01033 vector<TrkRecoTrk*> vec_trk_for_clean; 01034 vector<Mytrack> vec_myTrack; 01035 int track_fit_success=0; 01036 01037 for(track_fit=0;track_fit<peak_fit;track_fit++){ 01038 double d0_before_least=-999.; 01039 double phi0_before_least=-999.; 01040 double omega_before_least=-999.; 01041 map<int,double> hit_to_circle1; 01042 map<int,double> hit_to_circle2; 01043 map<int,double> hit_to_circle3; 01044 if(m_debug>0) cout<<"Do fitting track: "<<track_fit<<endl; 01045 01046 for(int i=0;i<nhit;i++) vec_truthHit[i]=false; 01047 01048 // confirm every track has which hits 01049 if(m_debug>1) cout<<"candi track["<<track_fit<<"] collect "<<peak_hit_num[track_fit]<<" hits "<<endl; 01050 01051 for(int i=0;i<peak_hit_num[track_fit];i++){ 01052 int hitNum=peak_hit_list[track_fit][i]; 01053 int hittemp=vec_layer[hitNum]*1000+vec_wire[hitNum]; 01054 vector<int>::iterator iter_ihit = find( vec_hit_use.begin(),vec_hit_use.end(),hittemp ); 01055 if( iter_ihit==vec_hit_use.end() ) vec_truthHit[hitNum]=true; 01056 if( m_debug >1 && vec_truthHit[hitNum]==true) cout<<" "<<"collect hit :"<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<" "<<endl;; 01057 } 01058 if( m_debug >1){ 01059 for(int i=0;i<nhit;i++){ 01060 if(vec_truthHit[i]==false){ 01061 cout<<"candi track "<<"uncollect hit: "<<endl; 01062 cout<<" "<<"uncollect hit :"<<"("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<" " <<endl; 01063 } 01064 } 01065 } 01066 // begin chisq fit 01067 int t_nHitAxialSelect=0; 01068 for(int i=0;i<nhit;i++){ 01069 if(vec_truthHit[i]==true){ 01070 if(vec_slant[i]==0) t_nHitAxialSelect++; 01071 } 01072 } 01073 if(m_debug>1) cout<<"track "<<track_fit<<" with axial select: "<<t_nHitAxialSelect<<endl; 01074 if(t_nHitAxialSelect<3 && m_drawTuple){ 01075 m_x_circle[track_fit]=-999.; 01076 m_y_circle[track_fit]=-999.; 01077 m_r[track_fit]=-999.; 01078 m_chis_pt[track_fit] = -999.; 01079 m_pt[track_fit]=-999.; 01080 continue; 01081 } 01082 double circle_x=0; 01083 double circle_y=0; 01084 double circle_r=0; 01085 int leastSquare=LeastSquare(nhit,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0_before_least,phi0_before_least,omega_before_least,circle_x,circle_y,circle_r); 01086 //hitdis to circle1 01087 for(int i=0;i<nhit;i++){ 01088 int hittemp=vec_layer[i]*1000+vec_wire[i]; 01089 double dist_temp=sqrt(pow( (vec_y[i]-circle_y),2)+pow( (vec_x[i]-circle_x),2))-circle_r; 01090 if(dist_temp>10.) vec_truthHit[i]=false; 01091 if(m_debug>1&&vec_truthHit[i]==true) cout<<"IN LEAST1: "<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<dist_temp<<endl; 01092 } 01093 01094 t_nHitAxialSelect=0; 01095 for(int i=0;i<nhit;i++){ 01096 if(vec_truthHit[i]==true){ 01097 if(vec_slant[i]==0) t_nHitAxialSelect++; 01098 } 01099 } 01100 if(m_debug>1) cout<<"track "<<track_fit<<"with axial2 select: "<<t_nHitAxialSelect<<endl; 01101 if(t_nHitAxialSelect<3) continue; 01102 01103 int leastSquare2=LeastSquare(nhit,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0_before_least,phi0_before_least,omega_before_least,circle_x,circle_y,circle_r); 01104 for(int i=0;i<nhit;i++){ 01105 int hittemp=vec_layer[i]*1000+vec_wire[i]; 01106 double dist_temp=sqrt(pow( (vec_y[i]-circle_y),2)+pow( (vec_x[i]-circle_x),2))-circle_r; 01107 hit_to_circle1[hittemp]=dist_temp; 01108 if(m_debug>1&&vec_truthHit[i]==true) cout<<" IN LEAST2: "<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<dist_temp<<endl; 01109 } 01110 01111 if(m_drawTuple){ 01112 m_x_circle[track_fit]=circle_x; 01113 m_y_circle[track_fit]=circle_y; 01114 m_r[track_fit]=circle_r; 01115 m_chis_pt[track_fit]=circle_r/333.567; 01116 } 01117 01118 // q select 01119 TrkHitList* qhits[2]; 01120 const TrkFit* newFit[2]; 01121 TrkHitList* qhits2[2]; 01122 const TrkFit* newFit2[2]; 01123 TrkRecoTrk* newTrack2[2]; 01124 01125 int Q_iq[]={-1,1}; 01126 double Q_Chis_2d[]={9999,9999}; 01127 double Q_Chis_3d[]={9999,9999}; 01128 double Q_z[]={999,999}; 01129 double Q_pt2[]={999,999}; 01130 double Q_tanl[]={999,999}; 01131 int Q_ok[2][2]={0}; 01132 01133 int q_start=0; 01134 int q_end=1; 01135 if(m_q==-1) {q_start=0;q_end=0;} 01136 if(m_q==1) {q_start=1;q_end=1;} 01137 for(int i_q=q_start;i_q<=q_end;i_q++){ 01138 double d0=d0_before_least; 01139 double phi0=phi0_before_least; 01140 double omega=omega_before_least; 01141 double z0=0; 01142 double tanl=0; 01143 if(m_debug>0 ){ 01144 cout<<"BY LEASTSQUARE: "<<endl; 01145 cout<<" d0: "<<d0<<" d0_mc "<<d0_mc<<endl; 01146 cout<<" phi0: "<<phi0<<" phi0_mc "<<phi0_mc<<endl; 01147 cout<<" omega0: "<<omega<<" omega_mc "<<omega_mc<<endl; 01148 } 01149 01150 vector<double> vec_flt_least; //use circle info 01151 for(int i=0;i<nhit;i++){ 01152 double theta_temp; 01153 if(circle_x==0||vec_x[i]-circle_x==0){ 01154 theta_temp=0; 01155 } 01156 else{ 01157 theta_temp=M_PI-atan2(vec_y[i]-circle_y,vec_x[i]-circle_x)+atan2(circle_y,circle_x); 01158 if(theta_temp>2*M_PI){ 01159 theta_temp=theta_temp-2*M_PI; 01160 } 01161 if(theta_temp<0){ 01162 theta_temp=theta_temp+2*M_PI; 01163 } 01164 } 01165 if(Q_iq[i_q]==-1) { 01166 double l_temp=circle_r*theta_temp; 01167 vec_flt_least.push_back(l_temp); 01168 } 01169 if(Q_iq[i_q]==1) { 01170 theta_temp=2*M_PI-theta_temp; 01171 double l_temp=circle_r*(theta_temp); 01172 vec_flt_least.push_back(l_temp); 01173 } 01174 } 01175 if(m_debug>2) { 01176 for(int i=0;i<nhit;i++){ 01177 cout<<"By least 2d: "<< "("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<vec_flt_least[i]<<endl; 01178 } 01179 } 01180 01181 //exchange par to 2d 01182 if(Q_iq[i_q]==-1){ 01183 d0=-d0; 01184 omega=-1./circle_r; 01185 } 01186 if(Q_iq[i_q]==1){ 01187 d0=d0; 01188 omega=1./circle_r; 01189 phi0=phi0-M_PI; 01190 } 01191 01192 // 2d global fit 01193 double x_cirtemp; 01194 double y_cirtemp; 01195 double r_temp; 01196 TrkExchangePar tt(d0,phi0,omega,z0,tanl); 01197 TrkCircleMaker circlefactory; 01198 float chisum =0.; 01199 TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0); 01200 //vec_trk_for_clean.push_back(newTrack); 01201 bool permitFlips = true; 01202 bool lPickHits = m_pickHits; 01203 circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits); 01204 int nDigi = digiToHots(mdcDigiVec,newTrack,vec_truthHit,getDigiFlag,vec_flt_least,hit_to_circle1,vec_for_clean); 01205 if(m_debug>2){ 01206 newTrack->printAll(std::cout); 01207 } 01208 //fit 01209 // TrkHitList* qhits = newTrack->hits(); 01210 qhits[i_q] = newTrack->hits(); 01211 TrkErrCode err=qhits[i_q]->fit(); 01212 newFit[i_q] = newTrack->fitResult(); 01213 01214 //test fit result 01215 if (!newFit[i_q] || (newFit[i_q]->nActive()<3) || err.failure()!=0) { 01216 if(m_debug>0){ 01217 cout << "evt "<<t_eventNum<<" global 2d fit failed "; 01218 if(newFit[i_q]) cout <<" nAct "<<newFit[i_q]->nActive(); 01219 cout<<"ERR1 failure ="<<err.failure()<<endl; 01220 } 01221 // m_2dFit[track_fit]=0; 01222 Q_ok[i_q][0]=0; 01223 Q_ok[i_q][1]=0; 01224 Q_Chis_2d[i_q]=-999.; 01225 delete newTrack; 01226 continue; 01227 }else{ 01228 Q_ok[i_q][0]=1; 01229 Q_ok[i_q][1]=0; 01230 if(m_debug>0) cout <<"evt "<<t_eventNum<< " global 2d fit success"<<endl; 01231 if(m_debug>2) { 01232 newTrack->print(std::cout); 01233 } 01234 Q_Chis_2d[i_q]=newFit[i_q]->chisq(); 01235 01236 TrkExchangePar par = newFit[i_q]->helix(0.); 01237 d0=par.d0(); 01238 phi0=par.phi0(); 01239 omega=par.omega(); 01240 r_temp=Q_iq[i_q]/par.omega(); 01241 x_cirtemp = sin(par.phi0()) *(par.d0()+1/par.omega()) * -1.; //z axis verse,x_babar = -x_belle 01242 y_cirtemp = -1.*cos(par.phi0()) *(par.d0()+1/par.omega()) * -1;//??? 01243 01244 if(m_debug>1) cout<<" circle after globle 2d: "<<"("<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<")"<<endl; 01245 if(m_debug>0) cout<<"pt_rec: "<<-1*Q_iq[i_q]/omega/333.567<<endl; 01246 01247 if(m_drawTuple){ 01248 m_x_circle[track_fit]=x_cirtemp; 01249 m_y_circle[track_fit]=y_cirtemp; 01250 m_r[track_fit]=r_temp; 01251 m_chis_pt[track_fit] = newFit[i_q]->chisq(); 01252 m_pt[track_fit]=r_temp/333.567; 01253 } 01254 } 01255 01256 int nfit2d=0; 01257 if(m_debug>1) cout<<" hot list:"<<endl; 01258 TrkHotList::hot_iterator hotIter= qhits[i_q]->hotList().begin(); 01259 int lay=((MdcHit*)(hotIter->hit()))->layernumber(); 01260 int wir=((MdcHit*)(hotIter->hit()))->wirenumber(); 01261 int hittemp=lay*1000+wir; 01262 while (hotIter!=qhits[i_q]->hotList().end()) { 01263 if(m_debug>1){ cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber() 01264 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 01265 <<":"<<hotIter->isActive()<<") "; 01266 } 01267 if (hotIter->isActive()==1) nfit2d++; 01268 hotIter++; 01269 } 01270 if(m_debug>0) cout<<"charge "<<i_q<<" In 2Dfit Num Of Active Hits: "<<nfit2d<<endl; 01271 01272 //caculate z position 01273 if(m_debug>0 && m_drawTuple){ 01274 cout<<" By global 2d charge "<<i_q<<endl; 01275 cout<<" d0: " <<d0<<" "<<m_drTruth[0]<<endl; 01276 cout<<" phi0: " <<phi0<<" "<<m_phi0Truth[0]<<endl; 01277 cout<<" omega: "<<omega<<" "<<m_omegaTruth[0]<<endl; 01278 cout<<" z0: " <<z0<<" "<<m_dzTruth[0]<<endl; 01279 cout<<" tanl: "<<tanl<<" "<<m_tanl_mc[0]<<endl; 01280 } 01281 01282 if(m_debug>2) { 01283 cout<<"least:( "<<circle_x<<","<<circle_y<<","<<circle_r<<endl; 01284 cout<<"2d:( "<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<endl; 01285 } 01286 delete newTrack; 01287 01288 //calcu every hit after 2dfit 01289 vector<double> vec_flt; 01290 vector<double> vec_theta_ontrack; 01291 for(int i=0;i<nhit;i++){ 01292 double theta_temp; 01293 if(x_cirtemp==0||vec_x[i]-x_cirtemp==0){ 01294 theta_temp=0; 01295 } 01296 else{ 01297 theta_temp=M_PI-atan2(vec_y[i]-y_cirtemp,vec_x[i]-x_cirtemp)+atan2(y_cirtemp,x_cirtemp); 01298 if(theta_temp>2*M_PI){ 01299 theta_temp=theta_temp-2*M_PI; 01300 } 01301 if(theta_temp<0){ 01302 theta_temp=theta_temp+2*M_PI; 01303 } 01304 if(Q_iq[i_q]==-1) { 01305 double l_temp=r_temp*theta_temp; 01306 vec_flt.push_back(l_temp); 01307 vec_theta_ontrack.push_back(theta_temp); 01308 // if(vec_truthHit[i]==true) cout<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):l:"<<vec_flt[i]<<" theta:"<<vec_theta_ontrack[i]<<endl; 01309 } 01310 if(Q_iq[i_q]==1) { 01311 theta_temp=2*M_PI-theta_temp; 01312 double l_temp=r_temp*(theta_temp); 01313 vec_flt.push_back(l_temp); 01314 vec_theta_ontrack.push_back(theta_temp); 01315 // if(vec_truthHit[i]==true) cout<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):l:"<<vec_flt[i]<<" theta:"<<vec_theta_ontrack[i]<<endl; 01316 } 01317 } 01318 } 01319 01320 if(m_debug>3){ 01321 for(int i=0;i<nhit;i++){ 01322 cout<<"i: "<<"theta: "<<vec_theta_ontrack[i]<<" l: "<<vec_flt[i]<<endl; 01323 } 01324 } 01325 01326 for(int i=0;i<nhit;i++){ 01327 int hittemp=vec_layer[i]*1000+vec_wire[i]; 01328 double dist_temp=sqrt(pow( (vec_y[i]-y_cirtemp),2)+pow( (vec_x[i]-x_cirtemp),2))-r_temp; 01329 hit_to_circle2[hittemp]=dist_temp; 01330 } 01331 01332 vector<double> vec_l_Calcu_Left; 01333 vector<double> vec_l_Calcu_Right; 01334 vector<double> vec_z_Calcu_Left; 01335 vector<double> vec_z_Calcu_Right; 01336 vector<double> vec_z_Calcu; 01337 vector<double> vec_l_Calcu; 01338 vector<double> vec_delta_z; 01339 vector<double> vec_delta_l; 01340 01341 double l_Calcu_Left=-999.; 01342 double l_Calcu_Right=-999.; 01343 double z_Calcu_Left=-999.; 01344 double z_Calcu_Right=-999.; 01345 double z_Calcu=-999.; 01346 double l_Calcu=-999.; 01347 double delta_z=-999.; 01348 double delta_l=-999.; 01349 // z caculate 01350 MdcDigiVec::iterator iter2 = mdcDigiVec.begin(); 01351 int zPos; 01352 int digiId = 0; 01353 for (;iter2 != mdcDigiVec.end(); iter2++, digiId++) { 01354 if(vec_truthHit[digiId]==false) continue; 01355 if(vec_slant[digiId]==0) continue; 01356 if(vec_theta_ontrack[digiId]>M_PI) continue; 01357 const MdcDigi* aDigi = *iter2; 01358 MdcHit *hit = new MdcHit(aDigi, m_gm); 01359 hit->setCalibSvc(m_mdcCalibFunSvc); 01360 hit->setCountPropTime(true); 01361 zPos=zPosition(*hit,m_bunchT0,x_cirtemp,y_cirtemp,r_temp,digiId,delta_z,delta_l,l_Calcu_Left,l_Calcu_Right,z_Calcu_Left,z_Calcu_Right,z_Calcu,l_Calcu,Q_iq[i_q]); 01362 if(m_debug>2) cout<<"in ZS calcu hitPos: "<<"("<<vec_layer[digiId]<<","<<vec_wire[digiId]<<")"<<" l_truth: "<<vec_flt_truth[digiId]<<" z_truth: "<<vec_ztruth[digiId]<<" l_Cal: "<<l_Calcu<<" z_Cal: "<<z_Calcu<<endl; 01363 delete hit; 01364 if(zPos==-1||zPos==-2) continue; 01365 vec_l_Calcu_Left.push_back(l_Calcu_Left); 01366 vec_l_Calcu_Right.push_back(l_Calcu_Right); 01367 vec_z_Calcu_Left.push_back(z_Calcu_Left); 01368 vec_z_Calcu_Right.push_back(z_Calcu_Right); 01369 vec_z_Calcu.push_back(z_Calcu); 01370 vec_l_Calcu.push_back(l_Calcu_Right); 01371 vec_delta_z.push_back(delta_z); 01372 vec_delta_l.push_back(delta_l); 01373 } 01374 int nHitUseInZs=vec_delta_z.size(); 01375 01376 if(m_drawTuple){ 01377 for(int i=0;i<nHitUseInZs;i++){ 01378 m_nHitStereo_use=vec_delta_z.size(); 01379 m_lCalcuLeft[i]=vec_l_Calcu_Left[i]; 01380 m_lCalcuRight[i]=vec_l_Calcu_Right[i]; 01381 m_zCalcuLeft[i]=vec_z_Calcu_Left[i]; 01382 m_zCalcuRight[i]=vec_z_Calcu_Right[i]; 01383 m_lCalcu[i]=vec_l_Calcu[i]; 01384 m_zCalcu[i]=vec_z_Calcu[i]; 01385 m_delta_z[i]=vec_delta_z[i]; 01386 m_delta_l[i]=vec_delta_l[i]; 01387 } 01388 } 01389 01390 //ambig choose 01391 vector< vector< vector <double> > > point(2, vector< vector<double> >(2,vector<double>() )); 01392 for(int iHit=0;iHit<nHitUseInZs;iHit++){ 01393 point[0][0].push_back(vec_l_Calcu_Left[iHit]); 01394 point[0][1].push_back(vec_z_Calcu_Left[iHit]); 01395 point[1][0].push_back(vec_l_Calcu_Right[iHit]); 01396 point[1][1].push_back(vec_z_Calcu_Right[iHit]); 01397 } 01398 vector<int> ambigList; 01399 if(nHitUseInZs!=0){ 01400 AmbigChoose(point,nHitUseInZs,ambigList,tanl,z0); 01401 if(m_drawTuple){ 01402 m_tanl_Cal=tanl; 01403 m_z0_Cal=z0; 01404 } 01405 01406 if(m_useMcInfo && m_drawTuple){ 01407 int ambig_no_match_num=0; 01408 for(int iHit=0;iHit<nHitUseInZs;iHit++){ 01409 if(ambigList[iHit]==0) ambigList[iHit]=1; 01410 else ambigList[iHit]=0; 01411 if(ambigList[iHit]!=m_postruth[iHit]) { m_ambig_no_match=1; ambig_no_match_num++; } 01412 } 01413 m_ambig_no_match_propotion=(double)(ambig_no_match_num)/m_nHitStereo_use; 01414 for (int iHit=0;iHit<ambigList.size();iHit++){ 01415 m_amb[iHit]=ambigList[iHit]; 01416 } 01417 } 01418 } 01419 01420 if(m_debug>0 && m_drawTuple){ 01421 cout<<"By 3d track zs fit:"<<endl; 01422 cout<<"d0: "<<d0<<" mc : "<<m_drTruth[0]<<endl; 01423 cout<<"phi0: "<<phi0<<" mc : "<<m_phi0Truth[0]<<endl; 01424 cout<<"omega: "<<omega<<" mc : "<<m_omegaTruth[0]<<endl; 01425 cout<<"z0: "<<z0<<" mc : "<<m_dzTruth[0]<<endl; 01426 cout<<"tanl: "<<tanl<<" mc : "<<m_tanl_mc[0]<<endl; 01427 } 01428 01429 if(m_mcpar==1){ 01430 d0=m_drTruth[0]; 01431 phi0=m_phi0Truth[0]; 01432 if(Q_iq[i_q]==-1) {phi0=m_phi0Truth[0]-3./2.*M_PI;omega=m_omegaTruth[0];} 01433 else {phi0=m_phi0Truth[0]-1./2.*M_PI;omega=-1*m_omegaTruth[0];} 01434 z0=m_dzTruth[0]; 01435 tanl=m_tanl_mc[0]; 01436 } 01437 01438 //-------------------------------------------5 parameter fit----------------------- 01439 01440 if(m_debug>0) cout<< "become 3d fit "<<endl; 01441 TrkExchangePar tt2(d0,phi0,omega,z0,tanl); 01442 TrkHelixMaker helixfactory; 01443 chisum =0.; 01444 newTrack2[i_q] = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0); 01445 vec_trk_for_clean.push_back(newTrack2[i_q]); 01446 permitFlips = true; 01447 lPickHits = m_pickHits; 01448 helixfactory.setFlipAndDrop(*newTrack2[i_q], permitFlips, lPickHits); 01449 int nDigi2 = digiToHots2(mdcDigiVec,newTrack2[i_q],vec_truthHit,vec_flt_truth,getDigiFlag,vec_slant,vec_l_Calcu,vec_flt,vec_theta_ontrack,vec_hitbeforefit,hit_to_circle2,vec_for_clean,tanl); 01450 if(m_debug>2){ 01451 cout<<__FILE__<<__LINE__<<"AFTER digiToHots 3d ---------BEFORE 3d fit"<<endl; 01452 newTrack2[i_q]->printAll(std::cout); 01453 } 01454 //fit 01455 qhits2[i_q] = newTrack2[i_q]->hits(); 01456 TrkErrCode err2=qhits2[i_q]->fit(); 01457 newFit2[i_q] = newTrack2[i_q]->fitResult(); 01458 int nActiveHit = newTrack2[i_q]->hots()->nActive(); 01459 int fitSuccess_3d = 0; 01460 if (!newFit2[i_q] || (newFit2[i_q]->nActive()<5) || err2.failure()!=0) { 01461 fitSuccess_3d=0; 01462 Q_ok[i_q][1]=0; 01463 Q_Chis_3d[i_q]=-999.; 01464 if(m_debug>0){ 01465 cout << "evt "<<t_eventNum<<" global 3d fit failed "; 01466 if(!newFit2[i_q]) cout<<" newfit2 point is NULL!!!" <<endl; 01467 if(newFit2[i_q]) cout <<" nAct "<<newFit2[i_q]->nActive(); 01468 cout<<"ERR2 failure ="<<err2.failure()<<endl; 01469 } 01470 }else{ 01471 //fit success 01472 TrkExchangePar par2 = newFit2[i_q]->helix(0.); 01473 if( abs( 1/(par2.omega()) )>0.001){ 01474 fitSuccess_3d=1; 01475 Q_Chis_3d[i_q]=newFit2[i_q]->chisq(); 01476 Q_ok[i_q][1]=1; 01477 if(m_debug>0) cout <<"evt "<<t_eventNum<< " global 3d fit success "<<err2.failure()<<endl; 01478 if(m_debug>2){ 01479 cout<<__FILE__<<__LINE__<<"AFTER fit 3d "<<endl; 01480 newTrack2[i_q]->print(std::cout); 01481 } 01482 } else { 01483 fitSuccess_3d=0; 01484 Q_ok[i_q][1]=0; 01485 Q_Chis_3d[i_q]=-999.; 01486 if(m_debug>2) cout<<"3dfit failure of omega "<<endl; 01487 } 01488 bool badQuality = false; 01489 //yzhang add 2016-06-15 01490 //test quality of track 01491 if(fabs(Q_Chis_3d[i_q])>m_qualityFactor*m_dropTrkChi2Cut ){ 01492 if(m_debug>0){ 01493 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2 "<<Q_Chis_3d[i_q]<<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2Cut <<std::endl; 01494 } 01495 badQuality=1; 01496 } 01497 if( fabs(par2.d0())>m_qualityFactor*m_dropTrkDrCut) { 01498 if(m_debug>0){ 01499 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by d0 "<<par2.d0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDrCut <<std::endl; 01500 } 01501 badQuality=1; 01502 } 01503 if( fabs(par2.z0())>m_qualityFactor*m_dropTrkDzCut) { 01504 if(m_debug>0){ 01505 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by z0 "<<par2.z0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDzCut <<std::endl; 01506 } 01507 badQuality=1; 01508 } 01509 if( (fabs(Q_Chis_3d[i_q])/qhits2[i_q]->nHit()) > m_qualityFactor*m_dropTrkChi2NdfCut) { 01510 if(m_debug>0){ 01511 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2/ndf"<<(fabs(Q_Chis_3d[i_q])/qhits2[i_q]->nHit()) <<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2NdfCut<<std::endl; 01512 } 01513 badQuality=1; 01514 } 01515 if( nActiveHit <= m_dropTrkNhitCut) { 01516 if(m_debug>0){ 01517 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by nhit"<<nActiveHit <<" <= "<<m_qualityFactor<<" * "<<m_dropTrkNhitCut<<std::endl; 01518 } 01519 badQuality=1; 01520 } 01521 if(badQuality) { 01522 fitSuccess_3d=0; 01523 Q_ok[i_q][1]=0; 01524 Q_Chis_3d[i_q]=-999.; 01525 if(m_debug>2) cout<<"3dfit failure of bad quality"<<endl; 01526 } 01527 //zhangy 01528 } 01529 01530 // -------------try to out put parameters--- 01531 if(fitSuccess_3d==1){ 01532 TrkExchangePar par2 = newFit2[i_q]->helix(0.); 01533 if(m_debug>0){ 01534 cout<<"BY 3d global fit charge "<<i_q<<endl; 01535 cout<<"d0: "<<par2.d0()<<endl; 01536 cout<<"phi0: "<<par2.phi0()<<endl; 01537 cout<<"w: "<<par2.omega()<<endl; 01538 cout<<"z: "<<par2.z0()<<endl; 01539 cout<<"tanl: "<<par2.tanDip()<<endl; 01540 } 01541 r_temp=Q_iq[i_q]/par2.omega(); 01542 x_cirtemp = sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle 01543 y_cirtemp = -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//??? 01544 if(m_debug>0) cout<<" circle after globle 3d: "<<"("<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<")"<<endl; 01545 double pxy=-1*r_temp/333.567; 01546 double pz=pxy*par2.tanDip(); 01547 double pxyz=sqrt(pxy*pxy+pz*pz); 01548 01549 Q_pt2[i_q]=fabs(pxy); 01550 Q_z[i_q]=par2.z0(); 01551 Q_tanl[i_q]=par2.tanDip(); 01552 } else{ 01553 continue; 01554 } 01555 //------------------------------------clean------------------------- 01556 01557 for(int i=0;i<nhit;i++){ 01558 int hittemp=vec_layer[i]*1000+vec_wire[i]; 01559 double dist_temp=sqrt(pow( (vec_y[i]-y_cirtemp),2)+pow( (vec_x[i]-x_cirtemp),2))-r_temp; 01560 hit_to_circle3[hittemp]=dist_temp; 01561 } 01562 } 01563 01564 if(m_debug>3){ 01565 cout<<"chis 2d: "<<Q_Chis_2d[0]<<" "<<Q_Chis_2d[1]<<endl; 01566 cout<<"chis 3d: "<<Q_Chis_3d[0]<<" "<<Q_Chis_3d[1]<<endl; 01567 cout<<"Z: "<<Q_z[0]<<" "<<Q_z[1]<<endl; 01568 cout<<"chis ok-: "<<Q_ok[0][0]<<" "<<Q_ok[0][1]<<endl; 01569 cout<<"chis ok+: "<<Q_ok[1][0]<<" "<<Q_ok[1][1]<<endl; 01570 } 01571 01572 int Q; 01573 int Q_judge[2]={0}; 01574 if(Q_ok[0][1]==1) Q_judge[0]=1; 01575 if(Q_ok[1][1]==1) Q_judge[1]=1; 01576 if(Q_judge[0]==1&&Q_judge[1]==0) Q=-1; 01577 if(Q_judge[0]==0&&Q_judge[1]==1) Q=1; 01578 if(Q_judge[0]==1&&Q_judge[1]==1) { 01579 if(fabs(Q_z[0])>fabs(Q_z[1])) Q=1; 01580 else Q=-1; 01581 } 01582 if(Q_judge[0]==0&&Q_judge[1]==0) {Q=0; continue;} 01583 int q_choose=0; 01584 if(Q==-1) q_choose=0; 01585 if(Q== 1) q_choose=1; 01586 TrkExchangePar par_final = newFit2[q_choose]->helix(0.); 01587 if(m_debug>0){ 01588 cout<<"after q choose By 3d global fit: "<<Q<<endl; 01589 cout<<"d0: "<<par_final.d0()<<endl; 01590 cout<<"phi0: "<<par_final.phi0()<<endl; 01591 cout<<"w: "<<par_final.omega()<<endl; 01592 cout<<"z: "<<par_final.z0()<<endl; 01593 cout<<"tanl: "<<par_final.tanDip()<<endl; 01594 } 01595 double pxy=-1*Q_iq[q_choose]/par_final.omega()/333.567; 01596 double pz=pxy*par_final.tanDip(); 01597 double pxyz=sqrt(pxy*pxy+pz*pz); 01598 01599 if(m_debug>0) cout<<"pt2_rec: "<<pxy<<endl; 01600 if(m_drawTuple){ 01601 m_pt2_rec_final[track_fit_success]=pxy; 01602 m_p_rec_final[track_fit_success]=pxyz; 01603 m_d0_final[track_fit_success]=par_final.d0(); 01604 m_phi0_final[track_fit_success]=par_final.phi0(); 01605 m_omega_final[track_fit_success]=par_final.omega(); 01606 m_z0_final[track_fit_success]=par_final.z0(); 01607 m_tanl_final[track_fit_success]=par_final.tanDip(); 01608 m_Q[track_fit_success]=Q; 01609 } 01610 01611 int nfit3d=0; 01612 if(m_debug>1) cout<<" hot list:"<<endl; 01613 TrkHotList::hot_iterator hotIter2= qhits2[q_choose]->hotList().begin(); 01614 while (hotIter2!=qhits2[q_choose]->hotList().end()) { 01615 int lay=((MdcHit*)(hotIter2->hit()))->layernumber(); 01616 int wir=((MdcHit*)(hotIter2->hit()))->wirenumber(); 01617 int hittemp=lay*1000+wir; 01618 if(m_debug>1){ cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber() 01619 <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber() 01620 <<":"<<hotIter2->isActive()<<") "; 01621 } 01622 if( hotIter2->isActive()==1) { 01623 nfit3d++; 01624 if(track_fit_success==0) vec_hitOnTrack.push_back(hittemp); //single track 01625 } 01626 hotIter2++; 01627 } 01628 if(m_debug>0) cout<<"In 3D fit Num Of Active Hits: "<<nfit3d<<endl; 01629 01630 track_fit_success++; 01631 int choise_temp,choise_no; 01632 if(Q==-1) choise_temp=0,choise_no=1; 01633 if(Q==1) choise_temp=1,choise_no=0; 01634 vec_newTrack.push_back( newTrack2[choise_temp] ); 01635 } //end loop of all track 01636 01637 for(int iTrack=0;iTrack<vec_newTrack.size();iTrack++){ 01638 const TrkFit* newFit_combine = vec_newTrack[iTrack]->fitResult(); 01639 TrkExchangePar par_combine= newFit_combine->helix(0.); 01640 double cx_combine=(-333.567/par_combine.omega()+par_combine.d0())*cos(par_combine.phi0()); 01641 double cy_combine=(-333.567/par_combine.omega()+par_combine.d0())*sin(par_combine.phi0()); 01642 double phi_combine=par_combine.phi0()-M_PI/2; 01643 double pxy_combine=1./par_combine.omega()/333.567; 01644 if(pxy_combine<0) phi_combine+=M_PI; 01645 if(phi_combine<0) phi_combine+=2*M_PI; 01646 if(phi_combine>2*M_PI) phi_combine-=2*M_PI; 01647 01648 Mytrack mytrack; 01649 mytrack.pt=pxy_combine; 01650 mytrack.charge=(int)(fabs(pxy_combine)/pxy_combine); 01651 mytrack.phi=phi_combine; 01652 mytrack.r=-333.567/par_combine.omega(); 01653 mytrack.x_cir=cx_combine; 01654 mytrack.y_cir=cy_combine; 01655 mytrack.newTrack=vec_newTrack[iTrack]; 01656 mytrack.use_in_tds=true; 01657 mytrack.vec_hit_on_trk=vec_hitOnTrack; 01658 // cout<<"size of hitontrk: "<<mytrack.vec_hit_on_trk.size()<<endl; 01659 vec_myTrack.push_back(mytrack); 01660 01661 } 01662 01663 std::sort(vec_myTrack.begin(),vec_myTrack.end(),compare); 01664 01665 for(int i=0;i<vec_newTrack.size();i++){ 01666 Mytrack *mytrack_i=&(vec_myTrack[i]); 01667 if(mytrack_i->use_in_tds==false) continue; 01668 for(int j=i+1;j<vec_newTrack.size();j++){ 01669 Mytrack *mytrack_j=&(vec_myTrack[j]); 01670 if(mytrack_j->use_in_tds==false) continue; 01671 if( fabs(mytrack_j->phi-mytrack_i->phi)<0.1 ) mytrack_j->use_in_tds=false; 01672 if(fabs(mytrack_j->r)*(1.-0.25) <= fabs(mytrack_i->r) && fabs(mytrack_i->r) <= fabs(mytrack_j->r)*(1.+0.25)){ 01673 if(fabs(mytrack_j->x_cir-mytrack_i->x_cir) <= fabs(mytrack_j->r)*(0.25)&& fabs(mytrack_j->y_cir-mytrack_i->y_cir) <= fabs(mytrack_j->r)*(0.25) ){ 01674 if(mytrack_j->charge==mytrack_i->charge) 01675 mytrack_j->use_in_tds=false; 01676 } 01677 } 01678 } 01679 } 01680 01681 int nTrack_tds=0; 01682 vector<Mytrack> vec_mytrack_in_tds; 01683 for(int i=0;i<track_fit_success;i++){ 01684 Mytrack mytrack=vec_myTrack[i]; 01685 if(mytrack.use_in_tds==false) continue; 01686 vec_mytrack_in_tds.push_back(mytrack); 01687 nTrack_tds++; 01688 } 01689 01690 if(m_drawTuple){ 01691 m_trackNum_tds=nTrack_tds; 01692 m_trackNum_fit=track_fit_success; 01693 } 01694 01695 // RecMdcTrackCol::iterator iteritrk = trackList->begin(); 01696 // for(;iteritrk!=trackList->end();iteritrk++){ 01697 // cout<<"IN PATTSF: "<<endl; 01698 // cout<<"parameter: "<<(*iteritrk)->helix(0)<<" "<<(*iteritrk)->helix(1)<<" "<<(*iteritrk)->helix(2)<<" "<<(*iteritrk)->helix(3)<<" "<<(*iteritrk)->helix(4)<<endl; 01699 // cout<<"chi: "<<(*iteritrk)->chi2()<<" ndof: "<<(*iteritrk)->ndof()<<"nster: "<<(*iteritrk)->nster()<<endl; 01700 // } 01701 01702 //print track in pattsf with bad vertex 01703 // RecMdcTrackCol::iterator iteritrk_pattsf = vec_trackList.begin(); 01704 // for(;iteritrk_pattsf!=vec_trackList.end();iteritrk_pattsf++){ 01705 // cout<<"in PATTSF LOST: "<<(*iteritrk_pattsf)->helix(0)<<" "<<(*iteritrk_pattsf)->helix(1)<<" "<<(*iteritrk_pattsf)->helix(2)<<" "<<(*iteritrk_pattsf)->helix(3)<<" "<<(*iteritrk_pattsf)->helix(4)<<endl; 01706 // } 01707 01708 int outputTrk=0; 01709 int itrack=0; 01710 for(int i=0;i<nTrack_tds;i++){ 01711 Mytrack mytrack=vec_mytrack_in_tds[i]; 01712 const TrkFit* newFit_tds= mytrack.newTrack->fitResult(); 01713 TrkExchangePar par_tds= newFit_tds->helix(0.); 01714 // cout<<"IN HOUGH: "<<endl; 01715 // cout<<" parameter: "<<par_tds.d0()<<" "<<par_tds.phi0()<<" "<<333.567*par_tds.omega()<<" "<<par_tds.z0()<<" "<<par_tds.tanDip()<<endl; 01716 // cout<<" chis: "<<newFit_tds->chisq()<<" ndof: "<<newFit_tds->nDof()<<" nMdc: "<<newFit_tds->nMdc()<<endl; 01717 double cr= 1./ par_tds.omega(); 01718 double cx= sin(par_tds.phi0()) *(par_tds.d0()+1/par_tds.omega()) * -1.; //z axis verse,x_babar = -x_belle 01719 double cy= -1.*cos(par_tds.phi0()) *(par_tds.d0()+1/par_tds.omega()) * -1;//??? 01720 01721 unsigned bestIndex = 0; 01722 double bestDiff = 1.0e+20; 01723 double cR, cX, cY; 01724 vector<double> a0,a1,a2,a3,a4; 01725 vector<double> zdelta; 01726 RecMdcTrackCol::iterator iteritrk = trackList->begin(); 01727 int itrk=0; 01728 for(;iteritrk!=trackList->end();iteritrk++){ 01729 double pt=(*iteritrk)->pxy(); 01730 a0.push_back( (*iteritrk)->helix(0) ); 01731 a1.push_back( (*iteritrk)->helix(1) ); 01732 a2.push_back( (*iteritrk)->helix(2) ); 01733 a3.push_back( (*iteritrk)->helix(3) ); 01734 a4.push_back( (*iteritrk)->helix(4) ); 01735 zdelta.push_back( par_tds.z0()-(*iteritrk)->helix(3) ); 01736 // cout<<"IN the NEW List: "<<endl; 01737 // cout<<" parameter: "<<a0[itrk]<<" "<<a1[itrk]<<" "<<a2[itrk]<<" "<<a3[itrk]<<" "<<a4[itrk]<<endl; 01738 // cout<<" chi: "<<(*iteritrk)->chi2()<<" ndof: "<<(*iteritrk)->ndof()<<"nster: "<<(*iteritrk)->nster()<<endl; 01739 cR=(-333.567)/a2[itrk]; 01740 cX=(cR+a0[itrk])*cos(a1[itrk]); 01741 cY=(cR+a0[itrk])*sin(a1[itrk]); 01742 01743 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){ 01744 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){ 01745 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY)); 01746 if(diff < bestDiff){ 01747 bestDiff = diff; 01748 bestIndex = itrk; 01749 } 01750 } 01751 } 01752 itrk++; 01753 } 01754 01755 if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no merge "<<endl; } 01756 else continue; 01757 01758 // TrkHitList* hitlist = mytrack.newTrack->hits(); 01759 // TrkHitList::hot_iterator hotIter= hitlist->begin(); 01760 // while (hotIter!=hitlist->end()) { 01761 // cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber() 01762 // <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 01763 // <<":"<<hotIter->isActive()<<") "; 01764 // 01765 // if (hotIter->isActive()==0) hitlist->remove() ; 01766 // hotIter++; 01767 // } 01768 01769 01770 MdcTrack* mdcTrack; 01771 mdcTrack= new MdcTrack(mytrack.newTrack); 01772 vec_track_in_tds.push_back(mytrack.newTrack); 01773 int tkStat = 4;//track find by Hough set stat=4 01774 int tkId = nTdsTk + itrack; 01775 mdcTrack->storeTrack(tkId, trackList, hitList, tkStat); 01776 if(m_drawTuple) m_hitOnTrk[outputTrk]=mdcTrack->track().hots()->nActive(); 01777 outputTrk++; 01778 delete mdcTrack; 01779 01780 double pxy=1/par_tds.omega()/333.567; 01781 double pz=pxy*par_tds.tanDip(); 01782 double pxyz=sqrt(pxy*pxy+pz*pz); 01783 if(m_drawTuple){ 01784 m_pt2_rec_tds[itrack]=pxy; 01785 m_p_rec_tds[itrack]=pxyz; 01786 m_d0_tds[itrack]=par_tds.d0(); 01787 m_phi0_tds[itrack]=par_tds.phi0(); 01788 m_omega_tds[itrack]=par_tds.omega(); 01789 m_z0_tds[itrack]=par_tds.z0(); 01790 m_tanl_tds[itrack]=par_tds.tanDip(); 01791 m_Q_tds[itrack]=pxy>0?1:-1; 01792 } 01793 itrack++; 01794 } 01795 if(m_drawTuple) m_outputTrk=outputTrk; 01796 int nTdsTk_temp=trackList->size(); 01797 if(m_debug>0) m_mdcPrintSvc->printRecMdcTrackCol(); 01798 01799 /* 01800 m_track_use_intrk1=vec_hit_use_intrack1.size(); 01801 int nNoise_intrk1=0; 01802 for(int ihit=0;ihit<vec_hit_use_intrack1.size();ihit++){ 01803 int noise_yes = hit_noise.count( vec_hit_use_intrack1[ihit] ); 01804 if(noise_yes==1) {nNoise_intrk1++; cout<<"noise hit in track1: "<<vec_hit_use_intrack1[ihit]<<endl;} 01805 } 01806 m_noise_intrk1=nNoise_intrk1; 01807 */ 01808 01809 01810 // m_nhitdis1=vec_hit_use.size(); 01811 // m_nhitdis2=vec_hitbeforefit.size(); 01812 // for(int ihit=0;ihit<vec_hit_use.size();ihit++){ 01813 // int hittemp=vec_hit_use[ihit]; 01814 // m_hit_dis1[ihit]=hit_to_circle2[hittemp]; 01815 // m_hit_dis1_3d[ihit]=hit_to_circle3[hittemp]; 01816 // } 01817 // for(int ihit=0;ihit<vec_hitbeforefit.size();ihit++){ 01818 // int hittemp=vec_hitbeforefit[ihit]; 01819 // m_hit_dis2[ihit]=hit_to_circle2[hittemp]; 01820 // m_hit_dis2_3d[ihit]=hit_to_circle3[hittemp]; 01821 // } 01822 01823 for(int i=0;i<vec_for_clean.size();i++){ 01824 delete vec_for_clean[i]; 01825 } 01826 for(int i=0;i<vec_trk_for_clean.size();i++){ 01827 //cout<<"size "<<vec_trk_for_clean.size()<<endl; 01828 //cout<<"vec_clean: "<<vec_trk_for_clean[i]<<endl; 01829 vector<TrkRecoTrk*>::iterator iterTrk =find(vec_track_in_tds.begin(),vec_track_in_tds.end(),vec_trk_for_clean[i]); 01830 if(iterTrk ==vec_track_in_tds.end() ) delete vec_trk_for_clean[i]; 01831 //for(int j=0;j<vec_newTrack.size();j++) cout<<"vec_newTrack: "<<vec_newTrack[j]<<endl; 01832 } 01833 delete h1; 01834 delete h2; 01835 01836 if(m_drawTuple) ntuplehit->write(); 01837 if(m_debug>0) cout<<"end event" <<endl; 01838 return StatusCode::SUCCESS; 01839 }
void HoughValidUpdate::FillCells | ( | TH2D * | h1, | |
int | n, | |||
int | binX, | |||
int | binY, | |||
vector< double > | vec_u, | |||
vector< double > | vec_v, | |||
vector< int > | vec_layer, | |||
vector< int > | vec_wire, | |||
vector< vector< int > > & | countij, | |||
vector< vector< vector< int > > > & | vec_selectNum, | |||
vector< vector< vector< int > > > & | vec_selectHit | |||
) | [private] |
void HoughValidUpdate::FillHist | ( | TH2D * | h1, | |
vector< double > | vec_u, | |||
vector< double > | vec_v, | |||
int | m_nHit, | |||
vector< vector< vector< int > > > & | vec_selectNum | |||
) | [private] |
Definition at line 2374 of file HoughValidUpdate.cxx.
References binX, binY, cos(), count, genRecEmupikp::i, ganga-rec::j, m_maxrho, M_PI, and sin().
Referenced by execute().
02374 { 02375 for(int n=0;n<nhit;n++){ 02376 // if(vec_digitruth[n]==0) continue; 02377 for(int i=0;i<binX;i++){ 02378 double binwidth=M_PI/binX; 02379 double binhigh=2*m_maxrho/binY; 02380 double theta=(i+0.5)*binwidth; 02381 double rho=vec_u[n]*cos(theta)+vec_v[n]*sin(theta); 02382 int j=(int)((rho+m_maxrho)/binhigh); 02383 int count=h1->GetBinContent(i+1,j+1); 02384 count+=1; 02385 h1->SetBinContent(i+1,j+1,count); 02386 // h1->Fill(theta,rho); 02387 vec_selectNum[i][j].push_back(n); 02388 } 02389 } 02390 }
void HoughValidUpdate::FillRhoTheta | ( | TH2D * | h1, | |
vector< double > | vec_theta, | |||
vector< double > | vec_rho | |||
) | [private] |
Definition at line 2391 of file HoughValidUpdate.cxx.
References binX, binY, count, genRecEmupikp::i, m_maxrho, and M_PI.
02391 { 02392 for(int i=0;i<vec_theta.size();i++){ 02393 double binwidth=M_PI/binX; 02394 double binhigh=2*m_maxrho/binY; 02395 int binxNum=vec_theta[i]/binwidth+1; 02396 int binyNum=(vec_rho[i]+m_maxrho)/binhigh+1; 02397 int count =h1->GetBinContent(binxNum,binyNum); 02398 count++; 02399 h1->SetBinContent(binxNum,binyNum,count); 02400 } 02401 }
StatusCode HoughValidUpdate::finalize | ( | ) |
Definition at line 1842 of file HoughValidUpdate.cxx.
References Bes_Common::INFO, m_bfield, m_context, and msgSvc().
01842 { 01843 MsgStream log(msgSvc(), name()); 01844 delete m_bfield ; 01845 delete m_context ; 01846 log << MSG::INFO << "in finalize()" << endreq; 01847 return StatusCode::SUCCESS; 01848 }
int HoughValidUpdate::GetMcInfo | ( | ) | [private] |
Definition at line 1912 of file HoughValidUpdate.cxx.
References Helix::a(), Helix::cosTheta(), Helix::dr(), Helix::dz(), calibUtil::ERROR, Helix, m_costaTruth, m_debug, m_decay, m_drawTuple, m_drTruth, m_dzTruth, m_nTrkMC, m_omegaTruth, m_phi0Truth, m_pid, m_pidTruth, m_pTruth, m_ptTruth, m_pzTruth, m_qTruth, m_tanl_mc, Helix::momentum(), msgSvc(), Helix::phi0(), pid, Helix::pivot(), Helix::pt(), Helix::radius(), deljobs::string, t_eventNum, t_runNum, Helix::tanl(), and Bes_Common::WARNING.
Referenced by execute().
01912 { 01913 MsgStream log(msgSvc(), name()); 01914 StatusCode sc; 01915 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol"); 01916 if (!mcParticleCol) { 01917 log << MSG::ERROR << "Could not find McParticle" << endreq; 01918 return -999; 01919 } 01920 01921 int itk = 0; 01922 int t_qTruth = 0; //zhangjin 01923 int t_pidTruth = -999; 01924 int t_McTkId = -999; 01925 Helix* mchelix; 01926 vector<int> vec_decay; 01927 McParticleCol::iterator iter_mc = mcParticleCol->begin(); 01928 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) { 01929 // if(!(*iter_mc)->primaryParticle()) continue; 01930 t_pidTruth = (*iter_mc)->particleProperty(); 01931 bool t_decay= (*iter_mc)->decayInFlight(); 01932 if(m_debug>0) cout<<"Decay : "<<t_decay<<endl; 01933 vec_decay.push_back(t_decay); 01934 01935 if(m_debug>2) cout<< "tk "<<itk<< " pid="<< t_pidTruth<< endl; 01936 if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; } //zhangjin 01937 if(itk>=10) break; 01938 01939 int pid = t_pidTruth; 01940 if( pid == 0 ) { 01941 log << MSG::WARNING << "Wrong particle id" <<endreq; 01942 }else{ 01943 IPartPropSvc* p_PartPropSvc; 01944 static const bool CREATEIFNOTTHERE(true); 01945 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE); 01946 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) { 01947 std::cout<< " Could not initialize Particle Properties Service" << std::endl; 01948 } 01949 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT(); 01950 std::string name; 01951 if( p_particleTable->particle(pid) ){ 01952 name = p_particleTable->particle(pid)->name(); 01953 t_qTruth = p_particleTable->particle(pid)->charge(); 01954 }else if( p_particleTable->particle(-pid) ){ 01955 name = "anti " + p_particleTable->particle(-pid)->name(); 01956 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge(); 01957 } 01958 if(m_debug>2) std::cout << " particle "<<name<<" charge= "<<t_qTruth<<std::endl; 01959 } 01960 01961 t_McTkId = (*iter_mc)->trackIndex(); 01962 double px = (*iter_mc)->initialFourMomentum().px();//GeV 01963 double py = (*iter_mc)->initialFourMomentum().py();//GeV 01964 double pz = (*iter_mc)->initialFourMomentum().pz();//GeV 01965 01966 01967 mchelix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);//cm 01968 mchelix->pivot( HepPoint3D(0.,0.,0.) ); 01969 //cout<<"pt, p : "<<mchelix->pt()<<" "<<mchelix->momentum()<<endl; 01970 01971 if(m_debug>2){ 01972 std::cout<<"Truth tk "<<itk<<" pid "<<pid<<" charge "<<t_qTruth 01973 << " helix "<< mchelix->a()<<" p "<<mchelix->momentum(0.)<<std::endl; 01974 } 01975 if(m_drawTuple){ 01976 m_pzTruth[itk]=pz; 01977 m_pidTruth[itk] = t_pidTruth; 01978 m_costaTruth[itk] = mchelix->cosTheta(); 01979 m_ptTruth[itk] = mchelix->pt(); 01980 m_pTruth[itk] = mchelix->momentum(0.).mag(); 01981 m_qTruth[itk] = t_qTruth; 01982 m_drTruth[itk] = mchelix->dr(); 01983 m_phi0Truth[itk] = mchelix->phi0(); 01984 m_omegaTruth[itk] =1./(mchelix->radius()); //Q 01985 m_dzTruth[itk] = mchelix->dz(); 01986 m_tanl_mc[itk] =mchelix->tanl(); 01987 } 01988 itk++; 01989 if(m_debug>2){ 01990 cout<<" d0: " <<" "<<mchelix->dr()<<endl; 01991 cout<<" phi0: " <<" "<<mchelix->phi0()<<endl; 01992 cout<<" omega: "<<" "<<1./(mchelix->radius())<<endl; 01993 cout<<" z0: " <<" "<<mchelix->dz()<<endl; 01994 cout<<" tanl: "<<" "<<mchelix->tanl()<<endl; 01995 cout<<" pt: " <<" "<<mchelix->pt()<<endl; //Q 01996 cout<<" costaTruth: " <<" "<<mchelix->cosTheta()<<endl; //Q 01997 } 01998 delete mchelix; 01999 } 02000 if(m_drawTuple){ 02001 vector<int>::iterator iter_idecay=find(vec_decay.begin(),vec_decay.end(),1 ); 02002 if (iter_idecay==vec_decay.end() ) m_decay=0; 02003 else m_decay=1; 02004 m_nTrkMC = itk; 02005 } 02006 if(itk!=1) { 02007 if(m_debug>2) std::cout<<"WARNING run "<<t_runNum<<" evt "<<t_eventNum<<" not single event. nTrkMc="<<itk<<std::endl; 02008 return -1; 02009 }else{ 02010 if(m_debug>2) std::cout<<"nTrkMc=1"<<std::endl; 02011 return 1; 02012 } 02013 02014 }
StatusCode HoughValidUpdate::initialize | ( | ) |
Definition at line 104 of file HoughValidUpdate.cxx.
References BField::bFieldNominal(), calibUtil::ERROR, Bes_Common::FATAL, Bes_Common::INFO, m_amb, m_ambig_no_match, m_ambig_no_match_propotion, m_axiallose, m_bfield, m_cell, m_cell_Mc, m_chis_3d_final, m_chis_pt, m_context, m_costaTruth, m_d0_final, m_d0_tds, TrkHelixFitter::m_debug, m_debug, m_decay, m_delta_l, m_delta_z, m_digi_truth, m_drawTuple, m_drTruth, m_dzTruth, m_e_delta, m_eventNum, m_hit_dis0, m_hit_dis1, m_hit_dis1_3d, m_hit_dis2, m_hit_dis2_3d, m_hitOnTrk, m_layer, m_layer_max, m_layer_Mc, m_layerNhit, m_lCalcu, m_lCalcuLeft, m_lCalcuRight, m_ltruth, m_mdcCalibFunSvc, m_mdcGeomSvc, m_mdcPrintSvc, m_nCross, m_nHit, m_nHit_Mc, m_nHitAxial, m_nHitAxial_Mc, m_nHitAxialSelect, m_nHitDigi, m_nhitdis0, m_nhitdis1, m_nhitdis2, m_nHitSelect, m_nHitStereo, m_nHitStereo_use, m_nHitStereoSelect, m_noise_intrk1, m_npeak_1, m_npeak_2, m_npeak_after_tracking, m_nTrkMC, m_omega_final, m_omega_tds, m_omegaTruth, m_outputTrk, m_p, m_p_rec_final, m_p_rec_tds, m_particleTable, m_pdtFile, m_phi0_final, m_phi0_tds, m_phi0Truth, m_pidTruth, m_pIMF, m_postruth, m_pt, m_pt2_rec_final, m_pt2_rec_tds, m_pTruth, m_ptTruth, m_pzTruth, m_Q, m_Q_tds, m_qTruth, m_r, m_rawDataProviderSvc, m_rho, m_runNum, m_slant, m_stereolose, m_tanl_Cal, m_tanl_final, m_tanl_mc, m_tanl_tds, m_theta, m_track_use_intrk1, m_trackFailure, m_trackNum, m_trackNum_fit, m_trackNum_Mc, m_trackNum_tds, m_u, m_v, m_x, m_x_circle, m_x_east, m_x_west, m_xsigma, m_xybin, m_xybinNum, m_y, m_y_circle, m_y_east, m_y_west, m_ysigma, m_z, m_z0_Cal, m_z0_final, m_z0_tds, m_z_east, m_z_west, m_zCalcu, m_zCalcuLeft, m_zCalcuRight, m_ztruth, msgSvc(), ntuplehit, ntupleSvc(), and Pdt::readMCppTable().
00104 { 00105 00106 MsgStream log(msgSvc(), name()); 00107 log << MSG::INFO << "in initialize()" << endreq; 00108 00109 StatusCode sc; 00110 //for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i]; 00111 00112 IPartPropSvc* p_PartPropSvc; 00113 static const bool CREATEIFNOTTHERE(true); 00114 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE); 00115 if (!sc.isSuccess() || 0 == p_PartPropSvc) { 00116 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq; 00117 return sc; 00118 } 00119 m_particleTable = p_PartPropSvc->PDT(); 00120 00121 // RawData 00122 IRawDataProviderSvc* irawDataProviderSvc; 00123 sc = service ("RawDataProviderSvc", irawDataProviderSvc); 00124 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc); 00125 if ( sc.isFailure() ){ 00126 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq; 00127 return StatusCode::FAILURE; 00128 00129 } 00130 00131 // Geometry 00132 IMdcGeomSvc* imdcGeomSvc; 00133 sc = service ("MdcGeomSvc", imdcGeomSvc); 00134 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc); 00135 if ( sc.isFailure() ){ 00136 log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq; 00137 return StatusCode::FAILURE; 00138 } 00139 00140 // MdcPrintSvc 00141 IMdcPrintSvc* iMdcPrintSvc; 00142 sc = service ("MdcPrintSvc", iMdcPrintSvc); 00143 m_mdcPrintSvc= dynamic_cast<MdcPrintSvc*> (iMdcPrintSvc); 00144 if ( sc.isFailure() ){ 00145 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq; 00146 return StatusCode::FAILURE; 00147 } 00148 00149 sc = service ("MagneticFieldSvc",m_pIMF); 00150 if(sc!=StatusCode::SUCCESS) { 00151 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq; 00152 } 00153 m_bfield = new BField(m_pIMF); 00154 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq; 00155 00156 m_context = new TrkContextEv(m_bfield); 00157 00158 Pdt::readMCppTable(m_pdtFile); 00159 00160 //Get MdcCalibFunSvc 00161 IMdcCalibFunSvc* imdcCalibSvc; 00162 sc = service ("MdcCalibFunSvc", imdcCalibSvc); 00163 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc); 00164 if ( sc.isFailure() ){ 00165 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq; 00166 return StatusCode::FAILURE; 00167 } 00168 00169 //initialize ntuplehit 00170 if(m_drawTuple){ 00171 NTuplePtr nt1(ntupleSvc(), "HoughValidUpdate/hit"); 00172 if ( nt1 ){ 00173 ntuplehit = nt1; 00174 } else { 00175 ntuplehit = ntupleSvc()->book("HoughValidUpdate/hit", CLID_ColumnWiseTuple, "hit"); 00176 if(ntuplehit){ 00177 ntuplehit->addItem("eventNum", m_eventNum); 00178 ntuplehit->addItem("runNum", m_runNum); 00179 ntuplehit->addItem("nHitMc", m_nHit_Mc,0, 6796); 00180 ntuplehit->addItem("nHitAxialMc", m_nHitAxial_Mc); 00181 ntuplehit->addItem("layerMc", m_nHit_Mc, m_layer_Mc); 00182 ntuplehit->addItem("cellMc", m_nHit_Mc, m_cell_Mc); 00183 00184 ntuplehit->addItem("nTrkMC", m_nTrkMC,0,10); 00185 ntuplehit->addItem("pidTruth", m_nTrkMC,m_pidTruth); 00186 ntuplehit->addItem("costaTruth", m_nTrkMC,m_costaTruth); 00187 ntuplehit->addItem("ptTruth", m_nTrkMC,m_ptTruth); 00188 ntuplehit->addItem("pzTruth", m_nTrkMC,m_pzTruth); 00189 ntuplehit->addItem("pTruth", m_nTrkMC,m_pTruth); 00190 ntuplehit->addItem("qTruth", m_nTrkMC,m_qTruth); 00191 ntuplehit->addItem("drTruth", m_nTrkMC,m_drTruth); 00192 ntuplehit->addItem("phiTruth", m_nTrkMC,m_phi0Truth); 00193 ntuplehit->addItem("omegaTruth", m_nTrkMC,m_omegaTruth); 00194 ntuplehit->addItem("dzTruth", m_nTrkMC,m_dzTruth); 00195 ntuplehit->addItem("tanlTruth", m_nTrkMC,m_tanl_mc); 00196 00197 ntuplehit->addItem("nHit", m_nHit,0, 6796); 00198 ntuplehit->addItem("nHitDigi", m_nHitDigi); 00199 ntuplehit->addItem("nHitAxial", m_nHitAxial); 00200 ntuplehit->addItem("nHitStereo", m_nHitStereo); 00201 ntuplehit->addItem("layer", m_nHit, m_layer); 00202 ntuplehit->addItem("cell", m_nHit, m_cell); 00203 ntuplehit->addItem("x_east", m_nHit, m_x_east); 00204 ntuplehit->addItem("y_east", m_nHit, m_y_east); 00205 ntuplehit->addItem("z_east", m_nHit, m_z_east); 00206 ntuplehit->addItem("x_west", m_nHit, m_x_west); 00207 ntuplehit->addItem("y_west", m_nHit, m_y_west); 00208 ntuplehit->addItem("z_west", m_nHit, m_z_west); 00209 ntuplehit->addItem("x", m_nHit, m_x); 00210 ntuplehit->addItem("y", m_nHit, m_y); 00211 ntuplehit->addItem("z", m_nHit, m_z); 00212 ntuplehit->addItem("u", m_nHit, m_u); 00213 ntuplehit->addItem("v", m_nHit, m_v); 00214 ntuplehit->addItem("p", m_nHit, m_p); 00215 ntuplehit->addItem("slant", m_nHit, m_slant); 00216 ntuplehit->addItem("layerNhit", 43, m_layerNhit); 00217 ntuplehit->addItem("layerMax", m_layer_max); 00218 ntuplehit->addItem("l_truth", m_nHit, m_ltruth); 00219 ntuplehit->addItem("z_truth", m_nHit, m_ztruth); 00220 ntuplehit->addItem("z_postruth", m_nHit, m_postruth); 00221 ntuplehit->addItem("digi_truth", m_nHit, m_digi_truth); 00222 ntuplehit->addItem("e_delta", m_e_delta); 00223 //hough map 00224 ntuplehit->addItem("nCross", m_nCross, 0, 125000); 00225 ntuplehit->addItem("rho", m_nCross, m_rho); 00226 ntuplehit->addItem("theta", m_nCross, m_theta); 00227 ntuplehit->addItem("xybinNum", m_xybinNum, 0,100000); 00228 ntuplehit->addItem("xybin", m_xybinNum, m_xybin); 00229 ntuplehit->addItem("xsigma", m_xsigma); 00230 ntuplehit->addItem("ysigma", m_ysigma); 00231 00232 ntuplehit->addItem("npeak_1", m_npeak_1); 00233 ntuplehit->addItem("npeak_2", m_npeak_2); 00234 ntuplehit->addItem("trackFailure", m_trackFailure); 00235 ntuplehit->addItem("npeak_after_tracking", m_npeak_after_tracking, 0,1000); 00236 ntuplehit->addItem("nHitSelect", m_npeak_after_tracking, m_nHitSelect); 00237 ntuplehit->addItem("nHitAxialSelect", m_npeak_after_tracking, m_nHitAxialSelect); 00238 ntuplehit->addItem("nHitStereoSelect", m_npeak_after_tracking, m_nHitStereoSelect); 00239 ntuplehit->addItem("naxiallose", m_npeak_after_tracking, m_axiallose); 00240 ntuplehit->addItem("nstereolose", m_npeak_after_tracking, m_stereolose); 00241 00242 ntuplehit->addItem("trackNum_Mc", m_trackNum_Mc, 0,100); 00243 ntuplehit->addItem("trackNum", m_trackNum, 0,100); 00244 ntuplehit->addItem("trackNum_fit", m_trackNum_fit, 0,100); 00245 ntuplehit->addItem("trackNum_tds", m_trackNum_tds, 0,100); 00246 ntuplehit->addItem("circleCenterX", m_trackNum, m_x_circle); 00247 ntuplehit->addItem("circleCenterY", m_trackNum, m_y_circle); 00248 ntuplehit->addItem("circleR", m_trackNum, m_r); 00249 ntuplehit->addItem("chis_pt", m_trackNum, m_chis_pt); 00250 ntuplehit->addItem("pt_rec", m_trackNum, m_pt); 00251 00252 ntuplehit->addItem("nHitStereo_use",m_nHitStereo_use,0,1000); 00253 ntuplehit->addItem("l_Calcu_Left", m_nHitStereo_use, m_lCalcuLeft); 00254 ntuplehit->addItem("l_Calcu_Right", m_nHitStereo_use, m_lCalcuRight); 00255 ntuplehit->addItem("z_Calcu_Left", m_nHitStereo_use, m_zCalcuLeft); 00256 ntuplehit->addItem("z_Calcu_Right", m_nHitStereo_use, m_zCalcuRight); 00257 ntuplehit->addItem("z_Calcu", m_nHitStereo_use, m_zCalcu); 00258 ntuplehit->addItem("l_Calcu", m_nHitStereo_use, m_lCalcu); 00259 ntuplehit->addItem("z_delta", m_nHitStereo_use, m_delta_z); 00260 ntuplehit->addItem("l_delta", m_nHitStereo_use, m_delta_l); 00261 ntuplehit->addItem("amb", m_nHitStereo_use, m_amb); 00262 ntuplehit->addItem("amb_no_match", m_ambig_no_match); 00263 ntuplehit->addItem("amb_no_match_pro", m_ambig_no_match_propotion); 00264 ntuplehit->addItem("tanl_Cal", m_tanl_Cal); 00265 ntuplehit->addItem("z0_Cal", m_z0_Cal); 00266 00267 ntuplehit->addItem("pt2_rec_final", m_trackNum_fit, m_pt2_rec_final); 00268 ntuplehit->addItem("p_rec_final", m_trackNum_fit, m_p_rec_final); 00269 ntuplehit->addItem("d0_final", m_trackNum_fit, m_d0_final); 00270 ntuplehit->addItem("phi0_final", m_trackNum_fit, m_phi0_final); 00271 ntuplehit->addItem("omega_final", m_trackNum_fit, m_omega_final); 00272 ntuplehit->addItem("z0_final", m_trackNum_fit, m_z0_final); 00273 ntuplehit->addItem("tanl_final", m_trackNum_fit, m_tanl_final); 00274 ntuplehit->addItem("chis_3d_final", m_trackNum_fit, m_chis_3d_final); 00275 ntuplehit->addItem("Q_final", m_trackNum_fit, m_Q); 00276 00277 ntuplehit->addItem("pt2_rec_tds", m_trackNum_tds, m_pt2_rec_tds); 00278 ntuplehit->addItem("p_rec_tds", m_trackNum_tds, m_p_rec_tds); 00279 ntuplehit->addItem("d0_tds", m_trackNum_tds, m_d0_tds); 00280 ntuplehit->addItem("phi0_tds", m_trackNum_tds, m_phi0_tds); 00281 ntuplehit->addItem("omega_tds", m_trackNum_tds, m_omega_tds); 00282 ntuplehit->addItem("z0_tds", m_trackNum_tds, m_z0_tds); 00283 ntuplehit->addItem("tanl_tds", m_trackNum_tds, m_tanl_tds); 00284 ntuplehit->addItem("Q_tds", m_trackNum_tds, m_Q_tds); 00285 00286 00287 ntuplehit->addItem("nhitdis0", m_nhitdis0, 0,2000); 00288 ntuplehit->addItem("nhitdis1", m_nhitdis1, 0,2000); 00289 ntuplehit->addItem("nhitdis2", m_nhitdis2, 0,2000); 00290 ntuplehit->addItem("hitdis0", m_nhitdis0, m_hit_dis0); 00291 ntuplehit->addItem("hitdis1", m_nhitdis1, m_hit_dis1); 00292 ntuplehit->addItem("hitdis2", m_nhitdis2, m_hit_dis2); 00293 ntuplehit->addItem("hitdis1_3d", m_nhitdis1, m_hit_dis1_3d); 00294 ntuplehit->addItem("hitdis2_3d", m_nhitdis2, m_hit_dis2_3d); 00295 ntuplehit->addItem("track_use_intrk1", m_track_use_intrk1); 00296 ntuplehit->addItem("noise_intrk1", m_noise_intrk1); 00297 ntuplehit->addItem("decay", m_decay); 00298 ntuplehit->addItem("outputTrk", m_outputTrk, 0,100); 00299 ntuplehit->addItem("hitOnTrk", m_outputTrk, m_hitOnTrk); 00300 00301 00302 } else { log << MSG::ERROR << "Cannot book tuple HoughValidUpdate/hough" << endmsg; 00303 return StatusCode::FAILURE; 00304 } 00305 } 00306 } 00307 00308 if(m_debug>2) TrkHelixFitter::m_debug = true; 00309 00310 return StatusCode::SUCCESS; 00311 00312 00313 }
int HoughValidUpdate::LeastLine | ( | int | , | |
double | x[], | |||
double | y[], | |||
double & | , | |||
double & | , | |||
double & | ||||
) | [private] |
Definition at line 2947 of file HoughValidUpdate.cxx.
References genRecEmupikp::i.
Referenced by AmbigChoose().
02947 { 02948 double x_sum=0; 02949 double y_sum=0; 02950 double x2_sum=0; 02951 double y2_sum=0; 02952 double xy_sum=0; 02953 for(int i=0;i<nhit;i++){ 02954 x_sum=x_sum+x[i]; 02955 y_sum=y_sum+y[i]; 02956 x2_sum=x2_sum+x[i]*x[i]; 02957 y2_sum=y2_sum+y[i]*y[i]; 02958 xy_sum=xy_sum+x[i]*y[i]; 02959 } 02960 b=(x2_sum*y_sum-xy_sum*x_sum)/(nhit*x2_sum-x_sum*x_sum); 02961 k=(nhit*xy_sum-x_sum*y_sum)/(nhit*x2_sum-x_sum*x_sum); 02962 02963 double yE[3]; 02964 for(int i=0;i<nhit;i++){ 02965 yE[i]=k*x[i]+b; 02966 double chi2_temp; 02967 if(yE[i]!=0) chi2_temp=(y[i]-yE[i])*(y[i]-yE[i])/yE[i]*yE[i]; 02968 else chi2_temp=0; 02969 chi2+=chi2_temp; 02970 } 02971 02972 return 1; 02973 }
int HoughValidUpdate::LeastSquare | ( | int | n, | |
vector< double > | n_x, | |||
vector< double > | n_y, | |||
vector< int > | n_slant, | |||
vector< int > | n_layer, | |||
vector< bool > | vec_truthHit, | |||
double & | d0, | |||
double & | phi0, | |||
double & | omega, | |||
double & | circleX, | |||
double & | circleY, | |||
double & | circleR | |||
) | [private] |
Definition at line 2017 of file HoughValidUpdate.cxx.
References genRecEmupikp::i, m_debug, and M_PI.
Referenced by execute().
02017 { 02018 02019 double x_sum=0; 02020 double y_sum=0; 02021 double x2_sum=0; 02022 double y2_sum=0; 02023 double x3_sum=0; 02024 double y3_sum=0; 02025 double x2y_sum=0; 02026 double xy2_sum=0; 02027 double xy_sum=0; 02028 double a1=0; 02029 double a2=0; 02030 double b1=0; 02031 double b2=0; 02032 double c1=0; 02033 double c2=0; 02034 double x_aver,y_aver,r2; 02035 int use_in_least=0; 02036 for(int i=0;i<n;i++){ 02037 if(vec_truthHit[i]==false) continue; 02038 if(vec_slant[i]!=0) continue; 02039 x_sum=x_sum+vec_x[i]; 02040 y_sum=y_sum+vec_y[i]; 02041 x2_sum=x2_sum+vec_x[i]*vec_x[i]; 02042 y2_sum=y2_sum+vec_y[i]*vec_y[i]; 02043 x3_sum=x3_sum+vec_x[i]*vec_x[i]*vec_x[i]; 02044 y3_sum=y3_sum+vec_y[i]*vec_y[i]*vec_y[i]; 02045 x2y_sum=x2y_sum+vec_x[i]*vec_x[i]*vec_y[i]; 02046 xy2_sum=xy2_sum+vec_x[i]*vec_y[i]*vec_y[i]; 02047 xy_sum=xy_sum+vec_x[i]*vec_y[i]; 02048 use_in_least++; 02049 } 02050 a1=x2_sum-x_sum*x_sum/n; 02051 a2=xy_sum-x_sum*y_sum/n; 02052 b1=a2; 02053 b2=y2_sum-y_sum*y_sum/n; 02054 c1=(x3_sum+xy2_sum-x_sum*(x2_sum+y2_sum)/n)/2.; 02055 c2=(y3_sum+x2y_sum-y_sum*(x2_sum+y2_sum)/n)/2.; 02056 02057 x_aver=x_sum/n; 02058 y_aver=y_sum/n; 02059 02060 double x_cirtemp =(b1*c2-b2*c1)/(b1*b1-a1*b2); 02061 double y_cirtemp =(b1*c1-a1*c2)/(b1*b1-a1*b2); 02062 r2=(x2_sum+y2_sum-2*x_cirtemp*x_sum-2*y_cirtemp*y_sum)/n+x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp; 02063 double r_temp=sqrt(r2); 02064 02065 circleX = x_cirtemp; 02066 circleY = y_cirtemp; 02067 circleR = r_temp; 02068 02069 double d0_temp=sqrt(x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp)-r_temp; 02070 double pt_temp=r_temp/333.567; 02071 if(m_debug>0){ 02072 cout<<" in least fit : "<<endl; 02073 cout<<"rtemp: "<<r_temp<<" "<<endl; 02074 cout<<"xtemp: "<<x_cirtemp<<" "<<endl; 02075 cout<<"ytemp: "<<y_cirtemp<<" "<<endl; 02076 cout<<"d0temp: "<<d0_temp<<" "<<endl; 02077 cout<<"pt_temp: "<<pt_temp<<endl; 02078 } 02079 02080 d0=d0_temp; //opposite with Q 02081 //phi0=atan2(y_cirtemp,x_cirtemp)-(m_q)*M_PI/2.; //charge- + ; charge+ -; 02082 phi0=atan2(y_cirtemp,x_cirtemp)+M_PI/2.; //charge- + ; charge+ -; 02083 omega=1/r_temp; //according with Q 02084 double z0=0.; 02085 double tanl=0.; 02086 if(m_debug>0) std::cout<<" before global fit Helix PatPar :"<<d0<<","<<phi0<<","<<omega<<","<<z0<<","<<tanl<< std::endl; 02087 return 0; 02088 02089 }
void HoughValidUpdate::Linefit | ( | int | z_stereoNum, | |
vector< double > | z_stereo_aver, | |||
vector< double > | l, | |||
double & | tanl, | |||
double & | z0 | |||
) | [private] |
Definition at line 2208 of file HoughValidUpdate.cxx.
References genRecEmupikp::i.
02208 { 02209 02210 double l_sum=0; 02211 double z_sum=0; 02212 double l2_sum=0; 02213 double z2_sum=0; 02214 double lz_sum=0; 02215 for(int i=0;i<z_stereoNum;i++){ 02216 // cout<<"event: "<<t_eventNum<<" "<<"z: "<<vec_zStereo[i]<<" "<<"l: "<<l[i]<<endl; 02217 l_sum=l_sum+l[i]; 02218 z_sum=z_sum+vec_zStereo[i]; 02219 l2_sum=l2_sum+l[i]*l[i]; 02220 // cout<<l2_sum<<endl; 02221 z2_sum=z2_sum+vec_zStereo[i]*vec_zStereo[i]; 02222 lz_sum=lz_sum+l[i]*vec_zStereo[i]; 02223 } 02224 double b=(l2_sum*z_sum-lz_sum*l_sum)/(z_stereoNum*l2_sum-l_sum); 02225 double k=(z_stereoNum*lz_sum-l_sum*z_sum)/(z_stereoNum*l2_sum-l_sum*l_sum); 02226 z0=b; 02227 tanl=k; 02228 cout<<k<<" "<<b<<endl; 02229 }
void HoughValidUpdate::Multi_array | ( | int | binX, | |
int | binY | |||
) | [private] |
Definition at line 2286 of file HoughValidUpdate.cxx.
References count, genRecEmupikp::i, and ganga-rec::j.
02286 { 02287 int **count=new int*[binX]; 02288 for(int i=0;i<binX;i++){ 02289 count[i]=new int [binY]; 02290 } 02291 //vector<int> vec_selectNum[binX][binY]; 02292 int ***selectNum=new int**[binX]; 02293 for(int i=0;i<binX;i++){ 02294 selectNum[i]=new int* [binY]; 02295 for(int j=0;j<binY;j++){ 02296 //selectNum[i][j]=new int[]; 02297 selectNum[i][j]=new int [count[i][j]]; 02298 } 02299 } 02300 }
void HoughValidUpdate::RhoTheta | ( | int | m_nHit, | |
vector< double > | vec_u, | |||
vector< double > | vec_v, | |||
vector< double > & | vec_rho, | |||
vector< double > & | vec_theta, | |||
vector< vector< int > > & | vec_hitNum, | |||
vector< int > | vec_digitruth | |||
) | [private] |
Definition at line 2302 of file HoughValidUpdate.cxx.
References genRecEmupikp::i, ganga-rec::j, m_nCross, M_PI, m_rho, and m_theta.
02302 { 02303 for(int i=0;i<nhit;i++){ 02304 if(vec_digitruth[i]==0) continue; 02305 for(int j=i+1;j<nhit;j++){ 02306 if(vec_digitruth[j]==0) continue; 02307 double k,b,x_cross,y_cross; 02308 double rho_temp,theta_temp; 02309 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]); 02310 b=vec_v[i]-k*vec_u[i]; 02311 x_cross=-b/(k+1/k); 02312 y_cross=b/(k*k+1); 02313 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross); 02314 theta_temp=atan2(y_cross,x_cross); 02315 if(theta_temp<0) { 02316 theta_temp=theta_temp+M_PI; 02317 rho_temp=-rho_temp; 02318 } 02319 vec_rho.push_back(rho_temp); 02320 vec_theta.push_back(theta_temp); 02321 vec_hitNum[0].push_back(i); 02322 vec_hitNum[1].push_back(j); 02323 //cout<<"u:"<<vec_u[i]<<" "<<"v:"<<vec_v[i]<<" "<<"("<<vec_theta[i]<<","<<vec_rho[i]<<")"<<endl; 02324 } 02325 } 02326 int nCross=vec_rho.size(); 02327 m_nCross=vec_rho.size(); 02328 if(m_nCross>125000) return ; 02329 // cout<<"numCross: "<<nCross<<endl; 02330 02331 for(int i=0;i<nCross;i++){ 02332 m_rho[i]=vec_rho[i]; 02333 m_theta[i]=vec_theta[i]; 02334 //cout<<" "<<vec_theta[i]<<" "<<vec_rho[i]<<endl; 02335 } 02336 }
void HoughValidUpdate::RhoThetaAxial | ( | vector< int > | vec_slant, | |
int | nHit, | |||
vector< double > | vec_u, | |||
vector< double > | vec_v, | |||
vector< double > & | vec_rho, | |||
vector< double > & | vec_theta, | |||
vector< vector< int > > & | vec_hitNum, | |||
vector< int > | vec_digitruth | |||
) | [private] |
Definition at line 2337 of file HoughValidUpdate.cxx.
References genRecEmupikp::i, ganga-rec::j, m_drawTuple, m_nCross, M_PI, m_rho, and m_theta.
02337 { 02338 for(int i=0;i<nhit;i++){ 02339 if(vec_digitruth[i]==0) continue; 02340 if(vec_slant[i]!=0) continue; 02341 for(int j=i+1;j<nhit;j++){ 02342 if(vec_slant[j]!=0) continue; 02343 if(vec_digitruth[j]==0) continue; 02344 double k,b,x_cross,y_cross; 02345 double rho_temp,theta_temp; 02346 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]); 02347 b=vec_v[i]-k*vec_u[i]; 02348 x_cross=-b/(k+1/k); 02349 y_cross=b/(k*k+1); 02350 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross); 02351 theta_temp=atan2(y_cross,x_cross); 02352 if(theta_temp<0) { 02353 theta_temp=theta_temp+M_PI; 02354 rho_temp=-rho_temp; 02355 } 02356 vec_rho.push_back(rho_temp); 02357 vec_theta.push_back(theta_temp); 02358 vec_hitNum[0].push_back(i); 02359 vec_hitNum[1].push_back(j); 02360 //cout<<"u:"<<vec_u[i]<<" "<<"v:"<<vec_v[i]<<" "<<"("<<vec_theta[i]<<","<<vec_rho[i]<<")"<<endl; 02361 } 02362 } 02363 int nCross=vec_rho.size(); 02364 if( m_drawTuple ) m_nCross=vec_rho.size(); 02365 02366 if(m_drawTuple){ 02367 for(int i=0;i<nCross;i++){ 02368 m_rho[i]=vec_rho[i]; 02369 m_theta[i]=vec_theta[i]; 02370 } 02371 } 02372 02373 }
int HoughValidUpdate::SlantId | ( | int | layer | ) | [private] |
Definition at line 1899 of file HoughValidUpdate.cxx.
Referenced by execute().
01899 { 01900 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) { 01901 if ((layer>=0&&layer<=3)||(layer>=20&&layer<=23)||(layer>=28&&layer<=31)){ 01902 // m_slant[digiId]=-1; 01903 return -1; 01904 }else{ 01905 return 1; 01906 } 01907 } else{ 01908 return 0; 01909 } 01910 }
int HoughValidUpdate::zPosition | ( | MdcHit & | h, | |
double | t0, | |||
double | x_cirtemp, | |||
double | y_cirtemp, | |||
double | r_temp, | |||
int | digiId, | |||
double & | delta_z, | |||
double & | delta_l, | |||
double & | l_Calcu_Left, | |||
double & | l_Calcu_Right, | |||
double & | z_Calcu_Left, | |||
double & | z_Calcu_Right, | |||
double & | z_Calcu, | |||
double & | l_Calcu, | |||
int & | i_q | |||
) | [private] |
Definition at line 2404 of file HoughValidUpdate.cxx.
References MdcHit::driftDist(), MdcSWire::getEastPoint(), MdcSWire::getWestPoint(), m_debug, m_drawTuple, m_ltruth, m_postruth, m_useMcInfo, m_ztruth, ORIGIN, pi, MdcHit::print(), unit, v, w, MdcHit::wire(), x, and MdcSWire::zAxis().
Referenced by execute().
02404 { 02405 //const MdcHit & h = *(hitUse.mdcHit()); 02406 02407 //ambig not sure 02408 //int ambig = hitUse.ambig(); 02409 int return_d[2]; 02410 int return_ok[2]; 02411 return_d[0]=1; 02412 return_d[1]=1; 02413 return_ok[0]=1; 02414 return_ok[1]=1; 02415 for(int ambig=-1;ambig<=1;ambig=ambig+2){ 02416 //int ambig =m_postruth[digiId]; 02417 //if(ambig==0) ambig =-1; 02418 HepPoint3D fp ( h.wire()->getWestPoint()->x(),h.wire()->getWestPoint()->y(),h.wire()->getWestPoint()->z()); 02419 HepPoint3D rp ( h.wire()->getEastPoint()->x(),h.wire()->getEastPoint()->y(),h.wire()->getEastPoint()->z()); 02420 02421 02422 HepVector3D X = 0.5 * (fp + rp); 02423 if(m_debug>0){ 02424 std::cout<<"---------- "<<std::endl;//yzhang debug 02425 h.print(std::cout); 02426 //h.wire()->print(std::cout); 02427 std::cout<<"fp rp "<<fp<<" "<<rp<<std::endl;//yzhang debug 02428 std::cout<<"Xmid "<<X<<std::endl;//yzhang debug 02429 } 02430 HepVector3D xx = HepVector3D(X.x(), X.y(), 0.); 02431 //center of helix/circle 02432 //change to belle param 02433 02434 /* 02435 double d0 = -par.d0(); 02436 double phi0 = (par.phi0()-pi/2); 02437 double _charge = -1; 02438 if((-1)*par.omega()/Bz > 0) _charge = 1; 02439 */ 02440 02441 //d0temp,phi0,omega,x_cirtemp,y_cirtemp, r_temp,m_q 02442 //double r; 02443 //if (fabs(par.omega())< Constants::epsilon) r = 9999.; 02444 //else r = 1/par.omega(); //+ - 02445 02446 //double xc = sin(par.phi0()) *(-d0+r) * -1.; //z axis verse,x_babar = -x_belle 02447 //double yc = -1.*cos(par.phi0()) *(-d0+r) * -1;//??? 02448 //HepPoint3D center (xc, yc, 0. ); 02449 HepPoint3D center (x_cirtemp, y_cirtemp, 0. ); 02450 02451 HepVector3D yy = center - xx; 02452 HepVector3D ww = HepVector3D(yy.x(), yy.y(), 0.); 02453 double wwmag2 = ww.mag2(); 02454 double wwmag = sqrt(wwmag2); 02455 02456 02457 //dirftDist 02458 double driftdist = fabs( h.driftDist(t0,ambig) ); 02459 HepVector3D lr(driftdist/wwmag * ww.x(), 02460 driftdist/wwmag * ww.y(), 0.); 02461 if (m_debug >4){ 02462 std::cout<<"xc "<<x_cirtemp << " yc "<<y_cirtemp<< " drfitdist "<<driftdist<<std::endl; 02463 //std::cout<<"r1 "<<r <<" d0 "<<d0 <<" phi0 "<<phi0<<std::endl;// 02464 //std::cout<<"lr "<<lr<<" hit ambig "<<hitUse.ambig()<<" ambig "<<ambig 02465 std::cout<<"lr "<<lr<<" "<<" ambig "<<ambig 02466 <<" left "<<h.driftDist(0,1) 02467 <<" right "<<h.driftDist(0,-1)<<std::endl; 02468 } 02469 02470 //...Check left or right... 02471 HepPoint3D ORIGIN = HepPoint3D(0., 0., 0.); 02472 //ambig 02473 //-1 right -phi direction driftDist-=lr 02474 //, +1 left +phi directon driftDist+=lr 02475 if (ambig == 0) lr = ORIGIN; 02476 // cout<<"m_q: "<<m_q<<" i_q: "<<i_q<<endl; 02477 if (i_q> 0){ //yzhang 02478 if (ambig == -1){ 02479 lr = - lr;//right 02480 } 02481 }else{ 02482 if (ambig == 1){ 02483 lr = - lr;//left 02484 } 02485 } 02486 X += lr; 02487 02488 //...Prepare vectors... 02489 // HepPoint3D center = _helix->center(); 02490 HepPoint3D tmp(-9999., -9999., 0.); 02491 HepVector3D x = HepVector3D(X.x(), X.y(), 0.); 02492 HepVector3D w = x - center; 02493 // //modified the next sentence because the direction are different from belle. 02494 HepVector3D V = h.wire()->zAxis(); 02495 HepVector3D v = HepVector3D(V.x(), V.y(), 0.); 02496 double vmag2 = v.mag2(); 02497 //double vmag = sqrt(vmag2); 02498 02499 //double r = _helix->curv(); 02500 double wv = w.dot(v); 02501 // wv = abs(wv); 02502 double d2 = wv * wv - vmag2 * (w.mag2() - r_temp * r_temp); 02503 if (m_debug >4){ 02504 std::cout<<"X_fix "<<X<<" center "<<center<<std::endl; 02505 std::cout<<"V "<<V<<std::endl;//yzhang debug 02506 std::cout<<"w "<<w<<" v "<<v<<std::endl; 02507 std::cout<<"wv "<<wv<<endl; 02508 cout<<"d2: "<<d2<<endl; 02509 } 02510 //...No crossing in R/Phi plane... This is too tight... 02511 02512 if (d2 < 0.) { 02513 //hitUse.position(tmp); 02514 if (m_debug>4){ 02515 cout<<"d2: "<<d2<<endl; 02516 std::cout << "in MdcSegInfoSterO !!! stereo: 0. > d2 = " << d2 << " " 02517 << ambig << std::endl; 02518 } 02519 //return -1; 02520 if (ambig==-1) return_d[ambig+1]=0; 02521 else return_d[ambig]=0; 02522 continue; 02523 } 02524 double d = sqrt(d2); 02525 // cout<<"d: "<<d<<endl; 02526 02527 02528 double l[2]; 02529 l[0] =-1*( (- wv + d) / vmag2 ); //multiply -1 ......? 02530 l[1] =-1*( (- wv - d) / vmag2 ); 02531 02532 double z[2]; 02533 //z[0] = X.z() + l[0] * V.z(); 02534 z[0] = X.z() - l[0] * V.z(); //l *-1 02535 z[1] = X.z() - l[1] * V.z(); 02536 02537 02538 //...Cal. length to crossing points... 02539 // double l[2]; 02540 //if (debug >0){ 02541 //std::cout<<"wv "<<wv<<" d "<<d<<" vmag2 "<<vmag2<<std::endl;//yzhang debug 02542 //} 02543 // l[0] = (- wv + d) / vmag2; 02544 02545 // l[1] = (- wv - d) / vmag2; 02546 02547 //...Cal. z of crossing points... 02548 bool ok[2]; 02549 ok[0] = true; 02550 ok[1] = true; 02551 //double z[2]; 02552 // z[0] = X.z() + l[0] * V.z(); 02553 //z[1] = X.z() + l[1] * V.z(); 02554 if (m_debug >4){//yzhang m_debug 02555 std::cout << "X.z "<<X.z()<<" V.z "<<V.z()<<std::endl; 02556 std::cout << "l0, l1 = " << l[0] << ", " << l[1] << std::endl; 02557 std::cout << "rpz " << rp.z() << " fpz " << fp.z() << std::endl; 02558 std::cout << "z0 " << z[0] << " z1 " << z[1] << std::endl; 02559 } 02560 //...Check z position... 02561 if (ambig == 0) { 02562 if(m_debug>4)std::cout<< " ambig = 0 " <<std::endl; 02563 if (z[0] > rp.z()+20. || z[0] < fp.z()-20.) { 02564 ok[0] = false; 02565 } 02566 if (z[1] > rp.z()+20. || z[1] < fp.z()-20.) { 02567 ok[1] = false; 02568 } 02569 } else { 02570 if(m_debug>4)std::cout<< " ambig != 0 " <<std::endl; 02571 // if (z[0] > rp.z() || z[0] < fp.z() ) 02572 if (fabs(z[0]/rp.z())>1.05 ) { ok[0] = false; } 02573 //if (z[1] > rp.z() || z[1] < fp.z() ) 02574 if (fabs(z[1]/rp.z())>1.05 ) { ok[1] = false; } 02575 } 02576 if ((! ok[0]) && (! ok[1])) { 02577 if(m_debug>4&&(ambig!=0)&&fabs(z[1]/rp.z())>1.05)std::cout<< " z[1] bad " << std::endl; 02578 if(m_debug>4&&(ambig!=0)&&fabs(z[0]/rp.z())>1.05)std::cout<< " z[0] bad " << std::endl; 02579 //hitUse.position(tmp); 02580 if(m_debug>4) { 02581 std::cout<< " z[1] bad " << "rpz " << rp.z() << " fpz " << fp.z() 02582 << "z0 " << z[0] << " z1 " << z[1] << std::endl; 02583 std::cout<< " !ok[0] && !ok[1] return -2" <<std::endl; 02584 } 02585 //return -2; 02586 if (ambig==-1) return_ok[ambig+1]=0; 02587 else return_ok[ambig]=0; 02588 continue; 02589 } 02590 if (( ok[0]) && ( ok[1])) { 02591 // cout<<" two z is ok ..........?" <<endl; 02592 } 02593 unsigned best = 0; 02594 if (ok[1]) best = 1; 02595 HepVector3D p[2]; 02596 p[0] = x + l[0] * v; 02597 p[1] = x + l[1] * v; 02598 double cosdPhi = - center.dot((p[best] - center).unit()) / center.mag(); 02599 double dPhi; 02600 if(fabs(cosdPhi)<=1.0) { 02601 dPhi = acos(cosdPhi); 02602 } else if (cosdPhi>1.0) { 02603 dPhi = 0.0; 02604 } else { 02605 dPhi = pi; 02606 } 02607 double stemp=r_temp*dPhi; 02608 if( m_debug >1) cout<<" ambig : "<<ambig<<" z: "<<z[best]<<" stemp: "<<stemp<<endl; 02609 double delta_temp=0; 02610 if(m_drawTuple && m_useMcInfo) delta_temp =z[best]-m_ztruth[digiId]; 02611 if( m_debug >3) cout<<"deltatemp: "<<delta_temp<<endl; 02612 // vec_delta_z.push_back(delta_temp); wrong cant push_back in circulate 02613 //vec_delta_z[digiId]=delta_temp; 02614 // delta_z=delta_temp; 02615 // cout<<" vec_delta in position: "<<delta_z<<endl; 02616 02617 // fill z_Cacu 02618 if(ambig==1){ 02619 z_Calcu_Left=z[best]; 02620 l_Calcu_Left=stemp; 02621 //cout<<"ambig: "<<ambig<<" z0 : "<<z[0]<<" z1: "<<z[1]<<" ztruth: "<<m_ztruth[digiId]<<endl; 02622 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" zbest : "<<z[best]<<" ztruth: "<<m_ztruth[digiId]<<endl; 02623 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" lbest : "<<stemp<<" ltruth: "<<m_ltruth[digiId]<<endl; 02624 if( m_useMcInfo && m_drawTuple && m_postruth[digiId]==1) { z_Calcu=z[best]; l_Calcu=stemp; delta_z= delta_temp; delta_l =l_Calcu-m_ltruth[digiId];} 02625 } 02626 if(ambig==-1){ 02627 z_Calcu_Right=z[best]; 02628 l_Calcu_Right=stemp; 02629 //cout<<"ambig: "<<ambig<<" z0 : "<<z[0]<<" z1: "<<z[1]<<" ztruth: "<<m_ztruth[digiId]<<endl; 02630 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" zbest : "<<z[best]<<" ztruth: "<<m_ztruth[digiId]<<endl; 02631 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" lbest : "<<stemp<<" ltruth: "<<m_ltruth[digiId]<<endl; 02632 if( m_useMcInfo && m_drawTuple && m_postruth[digiId]==0) { z_Calcu=z[best]; l_Calcu=stemp; delta_z= delta_temp; delta_l =l_Calcu-m_ltruth[digiId];} 02633 } 02634 02635 02637 02638 //if(m_debug>0){ 02639 // std::cout<<__FILE__<<" "<<__LINE__<< " p0 "<<p[0].x()<<" "<<p[0].y()<< std::endl; 02640 // std::cout<<__FILE__<<" "<<__LINE__<< " p1 "<<p[1].x()<<" "<<p[1].y()<< std::endl; 02641 // std::cout<<__FILE__<<" "<<__LINE__<< " c "<<center.x()<<" "<<center.y()<<" "<< std::endl; 02642 // std::cout<<__FILE__<<" "<<__LINE__<< " p0 centerx*y "<<center.x() * p[0].y()<<" centery*x"<<center.y() * p[0].x()<< std::endl; 02643 // std::cout<<__FILE__<<" "<<__LINE__<< " p1 centerx*y "<<center.x() * p[1].y()<<" centery*x"<<center.y() * p[1].x()<< std::endl; 02644 //} 02645 02646 02647 //...Which one is the best?... Study needed... 02648 //unsigned best = 0; 02649 //if (ok[1]) best = 1; 02650 //std::cout<<"in MdcSegInfoSterO Zbelle "<<z[best]<<std::endl;//yzhang m_debug 02651 02652 /* 02653 //...Cal. arc length... 02654 double cosdPhi = - center.dot((p[best] - center).unit()) / center.mag(); 02655 double dPhi; 02656 if(fabs(cosdPhi)<=1.0) { 02657 dPhi = acos(cosdPhi); 02658 } else if (cosdPhi>1.0) { 02659 dPhi = 0.0; 02660 } else { 02661 dPhi = pi; 02662 } 02663 02664 //...Finish... 02665 tmp.setX(abs(r * dPhi)); 02666 tmp.setY(z[best]); 02667 //hitUse.position(tmp); 02668 02669 //z 02670 span->y[hitFit] = z[best]; 02671 //if (ok[0]) span->y[hitFit] = p[0].mag();//r 02672 //else if(ok[1]) span->y[hitFit] = p[1].mag(); 02673 //else span->y[hitFit] = tmp.mag(); 02674 span->x[hitFit] = abs(r * dPhi); 02675 if (hitUse.ambig()<0) driftdist *= -1.; 02676 span->sigma[hitFit] = h.sigma(driftdist,hitUse.ambig()); 02677 02678 if (m_debug >0){ 02679 std::cout<<"("<<h.layernumber()<<","<<h.wirenumber()<<") s "<<span->x[hitFit]<<" z "<<span->y[hitFit]<<std::endl;//yzhang m_debug 02680 } 02681 return 1; 02682 */ 02683 } 02684 // cout<<"ooooooooooooooooooooooooooooooooooooooooo"<<endl; 02685 if(return_d[0]==0&&return_d[1]==0) return -1; 02686 if(return_ok[0]==0&&return_ok[1]==0) return -2; 02687 return 1; 02688 }
int HoughValidUpdate::Zposition | ( | int | n, | |
vector< int > | n_slant, | |||
double | x_cirtemp, | |||
double | y_cirtemp, | |||
double | r_temp, | |||
vector< double > | n_x_east, | |||
vector< double > | n_x_west, | |||
vector< double > | n_y_east, | |||
vector< double > | n_y_west, | |||
vector< double > | n_z_east, | |||
vector< double > | n_z_west, | |||
vector< int > | n_layer, | |||
vector< int > | n_wire, | |||
vector< double > | n_z, | |||
vector< double > & | z_stereo_aver, | |||
vector< double > & | l, | |||
vector< bool > | vec_truthHit | |||
) | [private] |
Definition at line 2091 of file HoughValidUpdate.cxx.
References abs, genRecEmupikp::i, and M_PI.
02091 { 02092 double x_stereo; 02093 double y_stereo; 02094 // double theta; 02095 int stereoId=0; 02096 //vector<double> z_stereo; 02097 // vector<double> vec_zStereo; 02098 int nDropDelta=0; 02099 for(int i=0;i<n;i++){ 02100 if(vec_truthHit[i]==false) continue; 02101 double k; 02102 double b; 02103 if(vec_slant[i]!=0){ 02104 if(vec_x_east[i]!=vec_x_west[i]){ 02105 k=(vec_y_east[i]-vec_y_west[i])/(vec_x_east[i]-vec_x_west[i]); 02106 b=-k*vec_x_east[i]+vec_y_east[i]; 02107 // cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<m_cell[i]<<" "<<"k: "<<k<<" "<<"b: "<<b<<" "<<m_y_east[i]<<" "<<m_y_west[i]<<" "<<m_x_east[i]<<" "<<m_x_west[i]<<" "<<endl; 02108 double delta=4*(k*b-k*y_cirtemp-x_cirtemp)*(k*b-k*y_cirtemp-x_cirtemp)-4*(k*k+1)*(x_cirtemp*x_cirtemp+b*b+y_cirtemp*y_cirtemp-2*b*y_cirtemp-r_temp*r_temp); 02109 //cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<endl; 02110 //cout<<"x_east: "<<vec_x_east[i]<<" "<<"x_west: "<<vec_x_west[i]<<endl; 02111 //cout<<"k: "<<k<<" "<<"b: "<<b<<endl; 02112 //cout<<"x_cirtemp: "<<x_cirtemp<<" "<<"y_cirtemp: "<<y_cirtemp<<endl; 02113 if(delta<0) { 02114 nDropDelta++; 02115 continue; 02116 } 02117 //cout<<"delta: "<<delta<<endl; 02118 double x1=(2*x_cirtemp+2*k*y_cirtemp-2*k*b+sqrt(delta))/(2*(k*k+1)); 02119 double x2=(2*x_cirtemp+2*k*y_cirtemp-2*k*b-sqrt(delta))/(2*(k*k+1)); 02120 double x_middle=(vec_x_east[i]+vec_x_west[i])/2.; 02121 // if(abs(x1-vec_x[i])>abs(x2-vec_x[i])) {x_stereo=x2;} 02122 if(abs(x1-x_middle)>abs(x2-x_middle)) {x_stereo=x2;} 02123 else {x_stereo=x1;} 02124 y_stereo=k*x_stereo+b; 02125 } 02126 else{ 02127 //cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<endl; 02128 //cout<<"x_east: "<<vec_x_east[i]<<" "<<"x_west: "<<vec_x_west[i]<<endl; 02129 x_stereo=vec_x_east[i]; 02130 double y1=sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo))+y_cirtemp; 02131 double y2=y_cirtemp-sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo)); 02132 double y_middle=(vec_y_east[i]+vec_y_west[i])/2.; 02133 if(abs(y1-y_middle)>abs(y2-y_middle)) {y_stereo=y2;} 02134 else{y_stereo=y1;} 02135 } 02136 // cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<"x_stereo: "<<x_stereo<<" "<<"y_stereo: "<<y_stereo<<endl; 02137 02138 double theta_temp; 02139 if(x_cirtemp==0||x_stereo-x_cirtemp==0){ 02140 theta_temp=0; 02141 } 02142 else{ 02143 theta_temp=M_PI-atan2(y_stereo-y_cirtemp,x_stereo-x_cirtemp)+atan2(y_cirtemp,x_cirtemp); 02144 if(theta_temp>2*M_PI){ 02145 theta_temp=theta_temp-2*M_PI; 02146 } 02147 if(theta_temp<0){ 02148 theta_temp=theta_temp+2*M_PI; 02149 } 02150 } 02151 //cout<<"thetatemp: "<<theta_temp<<endl; 02152 //if(theta_temp>2*M_PI) { 02153 // theta_temp=theta_temp-2*M_PI; 02154 //} 02155 //if(theta_temp<0.) { 02156 // theta_temp=theta_temp+2*M_PI; 02157 //} 02158 //cout<<"debug: if thetatemp is cout over"<<endl; 02159 l.push_back(r_temp*theta_temp); 02160 //cout<<"debug: if l is pushback"<<endl; 02161 //l.push_back(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp)); 02162 //l.push_back(r_temp*(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp))); 02163 02164 double d1=sqrt((x_stereo-vec_x_west[i])*(x_stereo-vec_x_west[i])+(y_stereo-vec_y_west[i])*(y_stereo-vec_y_west[i])); 02165 double d2=sqrt((vec_x_east[i]-vec_x_west[i])*(vec_x_east[i]-vec_x_west[i])+(vec_y_east[i]-vec_y_west[i])*(vec_y_east[i]-vec_y_west[i])); 02166 vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2); 02167 //cout<<"debug: if vec_zStereo is pushback"<<endl; 02168 //if(stereoId==0) {vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);} 02169 //if(stereoId>0){ 02170 // if(vec_layer[i]==vec_layer.at(i-1)) { 02171 // vec_zStereo[(stereoId-1)]=((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.); 02172 // vec_zStereo.push_back((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.); 02173 // } 02174 // else{ 02175 // vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2); 02176 // } 02177 //} 02178 02179 // double delphi=atan2(y_stereo,x_stereo)-par.phi0(); 02180 //cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<" "<<m_cell[i]<<" "<<"z_stereo: "<<z_stereo<<" "<<d1<<" "<<d2<<" "<<endl; 02181 //cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<m_cell[i]<<" "<<"k: "<<k<<" "<<"b: "<<b<<" "<<m_y_east[i]<<" "<<m_y_west[i]<<" "<<m_x_east[i]<<" "<<m_x_west[i]<<" "<<"x_stereo: "<<" "<<x_stereo<<" "<<"x: "<<vec_x[i]<<"y_stereo: "<<" "<<y_stereo<<" "<<"y: "<<vec_y[i]<<" "<<"zstereo: "<<z_stereo<<endl; 02182 // cout<<"layerid: "<<" "<<vec_layer[i]<<" "<<"cellid: "<<vec_wire[i]<<" "<<"x_stereo: "<<" "<<x_stereo<<" "<<"y_stereo: "<<" "<<y_stereo<<" "<<"zstereo: "<<z_stereo.at(stereoId)<<" "<<vec_zStereo.at(stereoId)<<" "<<"ztruth: "<<vec_z[i]<<endl; 02183 02184 // cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<"z: "<<vec_zStereo[stereoId]<<" "<<"l: "<<l[stereoId]<<" "<<vec_truthHit[i]<<endl; 02185 stereoId=stereoId+1; 02186 //cout<<"theta: "<<theta<<endl; 02187 } 02188 else {k=b=0;} 02189 } 02190 cout<<"nDropDelta: "<<nDropDelta<<endl; 02191 // cout<<"if for is finished : "<<endl; 02192 //int digiId=0; 02193 //for(int i=0;i<n;i++){ 02194 // if(vec_slant[i]!=0){ 02195 // cout<<"layerid: "<<" "<<vec_layer[i]<<" "<<"cellid: "<<vec_wire[i]<<" "<<"zstereo: "<<z_stereo.at(digiId)<<" "<<vec_zStereo.at(digiId)<<" "<<"ztruth: "<<vec_z[i]<<endl; 02196 // digiId++; 02197 // } 02198 //} 02199 //double sum_zstereo=0; 02200 //for(int i=0;i<stereoId;i++){ 02201 // sum_zstereo=sum_zstereo+z_stereo[i]; 02202 //} 02203 //cout<<"sum_zstereo: "<<sum_zstereo<<endl; 02204 return stereoId; 02205 }
int HoughValidUpdate::binX [private] |
Definition at line 104 of file HoughValidUpdate.h.
Referenced by Combine(), execute(), FillHist(), and FillRhoTheta().
int HoughValidUpdate::binY [private] |
Definition at line 105 of file HoughValidUpdate.h.
Referenced by Combine(), execute(), FillHist(), and FillRhoTheta().
double HoughValidUpdate::d0_mc [private] |
double HoughValidUpdate::dz_mc [private] |
NTuple::Array<int> HoughValidUpdate::m_2d_nFit [private] |
Definition at line 240 of file HoughValidUpdate.h.
NTuple::Array<int> HoughValidUpdate::m_2dFit [private] |
Definition at line 241 of file HoughValidUpdate.h.
NTuple::Array<int> HoughValidUpdate::m_3d_nFit [private] |
Definition at line 242 of file HoughValidUpdate.h.
NTuple::Array<int> HoughValidUpdate::m_3dFit [private] |
Definition at line 243 of file HoughValidUpdate.h.
NTuple::Array<int> HoughValidUpdate::m_amb [private] |
NTuple::Item<int> HoughValidUpdate::m_ambig_no_match [private] |
NTuple::Item<double> HoughValidUpdate::m_ambig_no_match_propotion [private] |
NTuple::Array<int> HoughValidUpdate::m_axiallose [private] |
BField* HoughValidUpdate::m_bfield [private] |
int HoughValidUpdate::m_binx [private] |
int HoughValidUpdate::m_biny [private] |
double HoughValidUpdate::m_bunchT0 [private] |
Definition at line 99 of file HoughValidUpdate.h.
Referenced by digiToHots(), digiToHots2(), and execute().
NTuple::Array<int> HoughValidUpdate::m_cell [private] |
NTuple::Array<int> HoughValidUpdate::m_cell_Mc [private] |
NTuple::Array<double> HoughValidUpdate::m_chis_3d_final [private] |
NTuple::Array<double> HoughValidUpdate::m_chis_3d_tds [private] |
Definition at line 261 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_chis_pt [private] |
bool HoughValidUpdate::m_combineTracking [private] |
std::string HoughValidUpdate::m_configFile [private] |
Definition at line 126 of file HoughValidUpdate.h.
TrkContextEv* HoughValidUpdate::m_context [private] |
Definition at line 129 of file HoughValidUpdate.h.
Referenced by execute(), finalize(), and initialize().
NTuple::Array<double> HoughValidUpdate::m_costaTruth [private] |
NTuple::Array<double> HoughValidUpdate::m_d0 [private] |
Definition at line 229 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_d0_final [private] |
NTuple::Array<double> HoughValidUpdate::m_d0_tds [private] |
int HoughValidUpdate::m_data [private] |
int HoughValidUpdate::m_debug [private] |
Definition at line 77 of file HoughValidUpdate.h.
Referenced by Combine(), digiToHots(), digiToHots2(), execute(), GetMcInfo(), HoughValidUpdate(), initialize(), LeastSquare(), and zPosition().
NTuple::Item<bool> HoughValidUpdate::m_decay [private] |
NTuple::Array<double> HoughValidUpdate::m_delta_l [private] |
NTuple::Array<double> HoughValidUpdate::m_delta_z [private] |
NTuple::Array<int> HoughValidUpdate::m_digi_truth [private] |
NTuple::Array<double> HoughValidUpdate::m_dist_to_circle [private] |
Definition at line 267 of file HoughValidUpdate.h.
int HoughValidUpdate::m_drawTuple [private] |
bool operator<(const HoughValidUpdate::Mytrack& a,const HoughValidUpdate::Mytrack& b);
Definition at line 73 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), HoughValidUpdate(), initialize(), RhoThetaAxial(), and zPosition().
bool HoughValidUpdate::m_dropHot [private] |
double HoughValidUpdate::m_dropTrkChi2Cut [private] |
double HoughValidUpdate::m_dropTrkChi2NdfCut [private] |
double HoughValidUpdate::m_dropTrkDrCut [private] |
double HoughValidUpdate::m_dropTrkDzCut [private] |
int HoughValidUpdate::m_dropTrkNhitCut [private] |
double HoughValidUpdate::m_dropTrkPtCut [private] |
NTuple::Array<double> HoughValidUpdate::m_drTruth [private] |
Definition at line 166 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), and initialize().
NTuple::Array<double> HoughValidUpdate::m_dzTruth [private] |
Definition at line 169 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), and initialize().
NTuple::Item<int> HoughValidUpdate::m_e_delta [private] |
int HoughValidUpdate::m_eventcut [private] |
NTuple::Item<int> HoughValidUpdate::m_eventNum [private] |
NTuple::Item<int> HoughValidUpdate::m_eventTime [private] |
Definition at line 151 of file HoughValidUpdate.h.
int HoughValidUpdate::m_fillTruth [private] |
double HoughValidUpdate::m_fltcut [private] |
Definition at line 89 of file HoughValidUpdate.h.
double HoughValidUpdate::m_fpro [private] |
Definition at line 85 of file HoughValidUpdate.h.
uint32_t HoughValidUpdate::m_getDigiFlag [private] |
const MdcDetector* HoughValidUpdate::m_gm [private] |
Definition at line 127 of file HoughValidUpdate.h.
Referenced by beginRun(), digiToHots(), digiToHots2(), and execute().
std::vector<float> HoughValidUpdate::m_helixHitsSigma [private] |
Definition at line 141 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_hit_dis0 [private] |
NTuple::Array<double> HoughValidUpdate::m_hit_dis1 [private] |
NTuple::Array<double> HoughValidUpdate::m_hit_dis1_3d [private] |
NTuple::Array<double> HoughValidUpdate::m_hit_dis2 [private] |
NTuple::Array<double> HoughValidUpdate::m_hit_dis2_3d [private] |
double HoughValidUpdate::m_hit_pro [private] |
NTuple::Array<int> HoughValidUpdate::m_hitOnTrk [private] |
bool HoughValidUpdate::m_keepBadTdc [private] |
bool HoughValidUpdate::m_keepUnmatch [private] |
NTuple::Array<int> HoughValidUpdate::m_layer [private] |
NTuple::Item<int> HoughValidUpdate::m_layer_max [private] |
NTuple::Array<int> HoughValidUpdate::m_layer_Mc [private] |
NTuple::Array<int> HoughValidUpdate::m_layerNhit [private] |
NTuple::Array<double> HoughValidUpdate::m_lCalcu [private] |
NTuple::Array<double> HoughValidUpdate::m_lCalcuLeft [private] |
NTuple::Array<double> HoughValidUpdate::m_lCalcuRight [private] |
NTuple::Array<double> HoughValidUpdate::m_ltruth [private] |
Definition at line 193 of file HoughValidUpdate.h.
Referenced by execute(), initialize(), and zPosition().
int HoughValidUpdate::m_maxMdcDigi [private] |
double HoughValidUpdate::m_maxrho [private] |
Definition at line 102 of file HoughValidUpdate.h.
Referenced by execute(), FillHist(), FillRhoTheta(), and HoughValidUpdate().
int HoughValidUpdate::m_mcpar [private] |
const MdcCalibFunSvc* HoughValidUpdate::m_mdcCalibFunSvc [private] |
Definition at line 131 of file HoughValidUpdate.h.
Referenced by digiToHots(), digiToHots2(), execute(), and initialize().
MdcGeomSvc* HoughValidUpdate::m_mdcGeomSvc [private] |
MdcPrintSvc* HoughValidUpdate::m_mdcPrintSvc [private] |
int HoughValidUpdate::m_method [private] |
int HoughValidUpdate::m_minMdcDigi [private] |
NTuple::Item<int> HoughValidUpdate::m_nCross [private] |
Definition at line 199 of file HoughValidUpdate.h.
Referenced by initialize(), RhoTheta(), and RhoThetaAxial().
int HoughValidUpdate::m_ndev1 [private] |
int HoughValidUpdate::m_ndev2 [private] |
NTuple::Array<int> HoughValidUpdate::m_nFitFailure [private] |
Definition at line 239 of file HoughValidUpdate.h.
NTuple::Item<int> HoughValidUpdate::m_nHit [private] |
NTuple::Item<int> HoughValidUpdate::m_nHit_Mc [private] |
NTuple::Item<int> HoughValidUpdate::m_nHitAxial [private] |
NTuple::Item<int> HoughValidUpdate::m_nHitAxial_Mc [private] |
NTuple::Array<int> HoughValidUpdate::m_nHitAxialSelect [private] |
double HoughValidUpdate::m_nHitcut_high [private] |
Definition at line 91 of file HoughValidUpdate.h.
double HoughValidUpdate::m_nHitcut_low [private] |
Definition at line 90 of file HoughValidUpdate.h.
NTuple::Item<int> HoughValidUpdate::m_nHitDigi [private] |
NTuple::Item<int> HoughValidUpdate::m_nhitdis0 [private] |
NTuple::Item<int> HoughValidUpdate::m_nhitdis1 [private] |
NTuple::Item<int> HoughValidUpdate::m_nhitdis2 [private] |
NTuple::Array<int> HoughValidUpdate::m_nHitSelect [private] |
NTuple::Item<int> HoughValidUpdate::m_nHitStereo [private] |
NTuple::Item<int> HoughValidUpdate::m_nHitStereo_use [private] |
NTuple::Array<int> HoughValidUpdate::m_nHitStereoSelect [private] |
NTuple::Item<int> HoughValidUpdate::m_noise_intrk1 [private] |
NTuple::Item<int> HoughValidUpdate::m_npeak_1 [private] |
NTuple::Item<int> HoughValidUpdate::m_npeak_2 [private] |
NTuple::Item<int> HoughValidUpdate::m_npeak_after_tracking [private] |
NTuple::Item<int> HoughValidUpdate::m_nTrkMC [private] |
NTuple::Array<double> HoughValidUpdate::m_omega [private] |
Definition at line 231 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_omega_final [private] |
NTuple::Array<double> HoughValidUpdate::m_omega_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_omegaTruth [private] |
Definition at line 168 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), and initialize().
NTuple::Item<int> HoughValidUpdate::m_outputTrk [private] |
NTuple::Array<double> HoughValidUpdate::m_p [private] |
NTuple::Array<double> HoughValidUpdate::m_p_rec_final [private] |
NTuple::Array<double> HoughValidUpdate::m_p_rec_tds [private] |
int HoughValidUpdate::m_par_use [private] |
Definition at line 75 of file HoughValidUpdate.h.
HepPDT::ParticleDataTable* HoughValidUpdate::m_particleTable [private] |
std::string HoughValidUpdate::m_pdtFile [private] |
Definition at line 87 of file HoughValidUpdate.h.
Referenced by HoughValidUpdate(), and initialize().
int HoughValidUpdate::m_peakCellNumX [private] |
int HoughValidUpdate::m_peakCellNumY [private] |
NTuple::Array<double> HoughValidUpdate::m_phi0 [private] |
Definition at line 230 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_phi0_final [private] |
NTuple::Array<double> HoughValidUpdate::m_phi0_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_phi0Truth [private] |
Definition at line 167 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), and initialize().
bool HoughValidUpdate::m_pickHits [private] |
int HoughValidUpdate::m_pid [private] |
Definition at line 138 of file HoughValidUpdate.h.
Referenced by GetMcInfo(), and HoughValidUpdate().
NTuple::Array<double> HoughValidUpdate::m_pidTruth [private] |
IMagneticFieldSvc* HoughValidUpdate::m_pIMF [private] |
NTuple::Array<double> HoughValidUpdate::m_postruth [private] |
Definition at line 195 of file HoughValidUpdate.h.
Referenced by execute(), initialize(), and zPosition().
NTuple::Array<double> HoughValidUpdate::m_pt [private] |
NTuple::Array<double> HoughValidUpdate::m_pt2 [private] |
Definition at line 235 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_pt2_rec_final [private] |
NTuple::Array<double> HoughValidUpdate::m_pt2_rec_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_pTruth [private] |
NTuple::Array<double> HoughValidUpdate::m_ptTruth [private] |
NTuple::Array<double> HoughValidUpdate::m_pxyz [private] |
Definition at line 237 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_pz [private] |
Definition at line 236 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_pzTruth [private] |
NTuple::Array<int> HoughValidUpdate::m_Q [private] |
int HoughValidUpdate::m_q [private] |
NTuple::Array<int> HoughValidUpdate::m_Q_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_qTruth [private] |
double HoughValidUpdate::m_qualityFactor [private] |
NTuple::Array<double> HoughValidUpdate::m_r [private] |
bool HoughValidUpdate::m_removeBadTrack [private] |
NTuple::Array<double> HoughValidUpdate::m_rho [private] |
Definition at line 200 of file HoughValidUpdate.h.
Referenced by initialize(), RhoTheta(), and RhoThetaAxial().
double HoughValidUpdate::m_rhocut [private] |
NTuple::Item<int> HoughValidUpdate::m_runNum [private] |
int HoughValidUpdate::m_setpeak_fit [private] |
NTuple::Array<int> HoughValidUpdate::m_slant [private] |
NTuple::Array<int> HoughValidUpdate::m_stereolose [private] |
NTuple::Array<double> HoughValidUpdate::m_tanl [private] |
Definition at line 233 of file HoughValidUpdate.h.
NTuple::Item<double> HoughValidUpdate::m_tanl_Cal [private] |
NTuple::Array<double> HoughValidUpdate::m_tanl_final [private] |
NTuple::Array<double> HoughValidUpdate::m_tanl_mc [private] |
Definition at line 170 of file HoughValidUpdate.h.
Referenced by execute(), GetMcInfo(), and initialize().
NTuple::Array<double> HoughValidUpdate::m_tanl_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_theta [private] |
Definition at line 201 of file HoughValidUpdate.h.
Referenced by initialize(), RhoTheta(), and RhoThetaAxial().
NTuple::Item<int> HoughValidUpdate::m_track_use_intrk1 [private] |
NTuple::Item<int> HoughValidUpdate::m_trackFailure [private] |
NTuple::Item<int> HoughValidUpdate::m_trackNum [private] |
NTuple::Item<int> HoughValidUpdate::m_trackNum_fit [private] |
NTuple::Item<int> HoughValidUpdate::m_trackNum_Mc [private] |
NTuple::Item<int> HoughValidUpdate::m_trackNum_tds [private] |
double HoughValidUpdate::m_truthPos[43][288][3] [private] |
Definition at line 265 of file HoughValidUpdate.h.
NTuple::Array<double> HoughValidUpdate::m_u [private] |
int HoughValidUpdate::m_useMcInfo [private] |
Definition at line 74 of file HoughValidUpdate.h.
Referenced by execute(), HoughValidUpdate(), and zPosition().
NTuple::Array<double> HoughValidUpdate::m_v [private] |
NTuple::Array<double> HoughValidUpdate::m_x [private] |
NTuple::Array<double> HoughValidUpdate::m_x_circle [private] |
NTuple::Array<double> HoughValidUpdate::m_x_east [private] |
NTuple::Array<double> HoughValidUpdate::m_x_west [private] |
NTuple::Item<double> HoughValidUpdate::m_xsigma [private] |
NTuple::Array<int> HoughValidUpdate::m_xybin [private] |
NTuple::Item<int> HoughValidUpdate::m_xybinNum [private] |
NTuple::Array<double> HoughValidUpdate::m_y [private] |
NTuple::Array<double> HoughValidUpdate::m_y_circle [private] |
NTuple::Array<double> HoughValidUpdate::m_y_east [private] |
NTuple::Array<double> HoughValidUpdate::m_y_west [private] |
NTuple::Item<double> HoughValidUpdate::m_ysigma [private] |
NTuple::Array<double> HoughValidUpdate::m_z [private] |
NTuple::Array<double> HoughValidUpdate::m_z0 [private] |
Definition at line 232 of file HoughValidUpdate.h.
NTuple::Item<double> HoughValidUpdate::m_z0_Cal [private] |
NTuple::Array<double> HoughValidUpdate::m_z0_final [private] |
NTuple::Array<double> HoughValidUpdate::m_z0_tds [private] |
NTuple::Array<double> HoughValidUpdate::m_z_east [private] |
NTuple::Array<double> HoughValidUpdate::m_z_west [private] |
NTuple::Array<double> HoughValidUpdate::m_zCalcu [private] |
NTuple::Array<double> HoughValidUpdate::m_zCalcuLeft [private] |
NTuple::Array<double> HoughValidUpdate::m_zCalcuRight [private] |
NTuple::Array<double> HoughValidUpdate::m_ztruth [private] |
Definition at line 194 of file HoughValidUpdate.h.
Referenced by execute(), initialize(), and zPosition().
int HoughValidUpdate::nfailure [private] |
Definition at line 139 of file HoughValidUpdate.h.
NTuple::Tuple* HoughValidUpdate::ntuplehit [private] |
double HoughValidUpdate::omega_mc [private] |
double HoughValidUpdate::phi0_mc [private] |
int HoughValidUpdate::t_eventNo [private] |
Definition at line 136 of file HoughValidUpdate.h.
int HoughValidUpdate::t_eventNum [private] |
int HoughValidUpdate::t_runNum [private] |
double HoughValidUpdate::tanl_mc [private] |
int HoughValidUpdate::track_fit [private] |