BOSS_7.0.2: HoughValidUpdate Class Reference

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]`

Definition at line 1895 of file HoughValidUpdate.cxx.

Referenced by execute().

01895                                                  {
01896   return 2*x/(x*x+y*y);
01897 }

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.

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]

Definition at line 107 of file HoughValidUpdate.h.

Referenced by execute().

double HoughValidUpdate::dz_mc [private]

Definition at line 110 of file HoughValidUpdate.h.

Referenced by execute().

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]

Definition at line 279 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_ambig_no_match [private]

Definition at line 280 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<double> HoughValidUpdate::m_ambig_no_match_propotion [private]

Definition at line 283 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_axiallose [private]

Definition at line 218 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

BField* HoughValidUpdate::m_bfield [private]

Definition at line 130 of file HoughValidUpdate.h.

Referenced by finalize(), and initialize().

int HoughValidUpdate::m_binx [private]

Definition at line 79 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

int HoughValidUpdate::m_biny [private]

Definition at line 80 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 177 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_cell_Mc [private]

Definition at line 156 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_chis_3d_final [private]

Definition at line 251 of file HoughValidUpdate.h.

Referenced by initialize().

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]

Definition at line 228 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

bool HoughValidUpdate::m_combineTracking [private]

Definition at line 140 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 161 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_d0 [private]

Definition at line 229 of file HoughValidUpdate.h.

NTuple::Array<double> HoughValidUpdate::m_d0_final [private]

Definition at line 246 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_d0_tds [private]

Definition at line 256 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_data [private]

Definition at line 78 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 301 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_delta_l [private]

Definition at line 278 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_delta_z [private]

Definition at line 277 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_digi_truth [private]

Definition at line 196 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 117 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

double HoughValidUpdate::m_dropTrkChi2Cut [private]

Definition at line 124 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

double HoughValidUpdate::m_dropTrkChi2NdfCut [private]

Definition at line 145 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

double HoughValidUpdate::m_dropTrkDrCut [private]

Definition at line 121 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

double HoughValidUpdate::m_dropTrkDzCut [private]

Definition at line 122 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

int HoughValidUpdate::m_dropTrkNhitCut [private]

Definition at line 144 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

double HoughValidUpdate::m_dropTrkPtCut [private]

Definition at line 123 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 197 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_eventcut [private]

Definition at line 92 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Item<int> HoughValidUpdate::m_eventNum [private]

Definition at line 150 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_eventTime [private]

Definition at line 151 of file HoughValidUpdate.h.

int HoughValidUpdate::m_fillTruth [private]

Definition at line 96 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 114 of file HoughValidUpdate.h.

Referenced by HoughValidUpdate().

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]

Definition at line 288 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_hit_dis1 [private]

Definition at line 289 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_hit_dis1_3d [private]

Definition at line 291 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_hit_dis2 [private]

Definition at line 290 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_hit_dis2_3d [private]

Definition at line 292 of file HoughValidUpdate.h.

Referenced by initialize().

double HoughValidUpdate::m_hit_pro [private]

Definition at line 86 of file HoughValidUpdate.h.

Referenced by Combine(), and HoughValidUpdate().

NTuple::Array<int> HoughValidUpdate::m_hitOnTrk [private]

Definition at line 303 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

bool HoughValidUpdate::m_keepBadTdc [private]

Definition at line 116 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

bool HoughValidUpdate::m_keepUnmatch [private]

Definition at line 118 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<int> HoughValidUpdate::m_layer [private]

Definition at line 176 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_layer_max [private]

Definition at line 191 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_layer_Mc [private]

Definition at line 155 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<int> HoughValidUpdate::m_layerNhit [private]

Definition at line 192 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_lCalcu [private]

Definition at line 276 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_lCalcuLeft [private]

Definition at line 271 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_lCalcuRight [private]

Definition at line 272 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 115 of file HoughValidUpdate.h.

Referenced by HoughValidUpdate().

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]

Definition at line 95 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 134 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

MdcPrintSvc* HoughValidUpdate::m_mdcPrintSvc [private]

Definition at line 135 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_method [private]

Definition at line 76 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

int HoughValidUpdate::m_minMdcDigi [private]

Definition at line 119 of file HoughValidUpdate.h.

Referenced by HoughValidUpdate().

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]

Definition at line 83 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

int HoughValidUpdate::m_ndev2 [private]

Definition at line 84 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<int> HoughValidUpdate::m_nFitFailure [private]

Definition at line 239 of file HoughValidUpdate.h.

NTuple::Item<int> HoughValidUpdate::m_nHit [private]

Definition at line 172 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nHit_Mc [private]

Definition at line 153 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nHitAxial [private]

Definition at line 174 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nHitAxial_Mc [private]

Definition at line 154 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_nHitAxialSelect [private]

Definition at line 216 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 173 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nhitdis0 [private]

Definition at line 285 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_nhitdis1 [private]

Definition at line 286 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_nhitdis2 [private]

Definition at line 287 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<int> HoughValidUpdate::m_nHitSelect [private]

Definition at line 215 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nHitStereo [private]

Definition at line 175 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nHitStereo_use [private]

Definition at line 270 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_nHitStereoSelect [private]

Definition at line 217 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_noise_intrk1 [private]

Definition at line 306 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_npeak_1 [private]

Definition at line 211 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_npeak_2 [private]

Definition at line 212 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_npeak_after_tracking [private]

Definition at line 214 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_nTrkMC [private]

Definition at line 159 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_omega [private]

Definition at line 231 of file HoughValidUpdate.h.

NTuple::Array<double> HoughValidUpdate::m_omega_final [private]

Definition at line 248 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_omega_tds [private]

Definition at line 258 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 302 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_p [private]

Definition at line 189 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_p_rec_final [private]

Definition at line 245 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_p_rec_tds [private]

Definition at line 255 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_par_use [private]

Definition at line 75 of file HoughValidUpdate.h.

HepPDT::ParticleDataTable* HoughValidUpdate::m_particleTable [private]

Definition at line 128 of file HoughValidUpdate.h.

Referenced by initialize().

std::string HoughValidUpdate::m_pdtFile [private]

Definition at line 87 of file HoughValidUpdate.h.

Referenced by HoughValidUpdate(), and initialize().

int HoughValidUpdate::m_peakCellNumX [private]

Definition at line 81 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

int HoughValidUpdate::m_peakCellNumY [private]

Definition at line 82 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<double> HoughValidUpdate::m_phi0 [private]

Definition at line 230 of file HoughValidUpdate.h.

NTuple::Array<double> HoughValidUpdate::m_phi0_final [private]

Definition at line 247 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_phi0_tds [private]

Definition at line 257 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 88 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 160 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

IMagneticFieldSvc* HoughValidUpdate::m_pIMF [private]

Definition at line 137 of file HoughValidUpdate.h.

Referenced by initialize().

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]

Definition at line 234 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_pt2 [private]

Definition at line 235 of file HoughValidUpdate.h.

NTuple::Array<double> HoughValidUpdate::m_pt2_rec_final [private]

Definition at line 244 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_pt2_rec_tds [private]

Definition at line 254 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_pTruth [private]

Definition at line 164 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_ptTruth [private]

Definition at line 162 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

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]

Definition at line 163 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_Q [private]

Definition at line 252 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_q [private]

Definition at line 120 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<int> HoughValidUpdate::m_Q_tds [private]

Definition at line 262 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_qTruth [private]

Definition at line 165 of file HoughValidUpdate.h.

Referenced by GetMcInfo(), and initialize().

double HoughValidUpdate::m_qualityFactor [private]

Definition at line 142 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<double> HoughValidUpdate::m_r [private]

Definition at line 227 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

RawDataProviderSvc* HoughValidUpdate::m_rawDataProviderSvc [private]

Definition at line 133 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

bool HoughValidUpdate::m_removeBadTrack [private]

Definition at line 143 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

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]

Definition at line 93 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Item<int> HoughValidUpdate::m_runNum [private]

Definition at line 152 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

int HoughValidUpdate::m_setpeak_fit [private]

Definition at line 94 of file HoughValidUpdate.h.

Referenced by execute(), and HoughValidUpdate().

NTuple::Array<int> HoughValidUpdate::m_slant [private]

Definition at line 190 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<int> HoughValidUpdate::m_stereolose [private]

Definition at line 219 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_tanl [private]

Definition at line 233 of file HoughValidUpdate.h.

NTuple::Item<double> HoughValidUpdate::m_tanl_Cal [private]

Definition at line 281 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_tanl_final [private]

Definition at line 250 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 260 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 305 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_trackFailure [private]

Definition at line 213 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_trackNum [private]

Definition at line 222 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_trackNum_fit [private]

Definition at line 223 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<int> HoughValidUpdate::m_trackNum_Mc [private]

Definition at line 221 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_trackNum_tds [private]

Definition at line 224 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

double HoughValidUpdate::m_truthPos[43][288][3] [private]

Definition at line 265 of file HoughValidUpdate.h.

NTuple::Array<double> HoughValidUpdate::m_u [private]

Definition at line 187 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 188 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_x [private]

Definition at line 184 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_x_circle [private]

Definition at line 225 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_x_east [private]

Definition at line 178 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_x_west [private]

Definition at line 180 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<double> HoughValidUpdate::m_xsigma [private]

Definition at line 204 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<int> HoughValidUpdate::m_xybin [private]

Definition at line 203 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Item<int> HoughValidUpdate::m_xybinNum [private]

Definition at line 202 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_y [private]

Definition at line 185 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_y_circle [private]

Definition at line 226 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_y_east [private]

Definition at line 179 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_y_west [private]

Definition at line 181 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Item<double> HoughValidUpdate::m_ysigma [private]

Definition at line 205 of file HoughValidUpdate.h.

Referenced by initialize().

NTuple::Array<double> HoughValidUpdate::m_z [private]

Definition at line 186 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_z0 [private]

Definition at line 232 of file HoughValidUpdate.h.

NTuple::Item<double> HoughValidUpdate::m_z0_Cal [private]

Definition at line 282 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_z0_final [private]

Definition at line 249 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_z0_tds [private]

Definition at line 259 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_z_east [private]

Definition at line 182 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_z_west [private]

Definition at line 183 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_zCalcu [private]

Definition at line 275 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_zCalcuLeft [private]

Definition at line 273 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

NTuple::Array<double> HoughValidUpdate::m_zCalcuRight [private]

Definition at line 274 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

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]

Definition at line 149 of file HoughValidUpdate.h.

Referenced by execute(), and initialize().

double HoughValidUpdate::omega_mc [private]

Definition at line 109 of file HoughValidUpdate.h.

Referenced by execute().

double HoughValidUpdate::phi0_mc [private]

Definition at line 108 of file HoughValidUpdate.h.

Referenced by execute().

int HoughValidUpdate::t_eventNo [private]

Definition at line 136 of file HoughValidUpdate.h.

int HoughValidUpdate::t_eventNum [private]

Definition at line 100 of file HoughValidUpdate.h.

Referenced by execute(), and GetMcInfo().

int HoughValidUpdate::t_runNum [private]

Definition at line 101 of file HoughValidUpdate.h.

Referenced by execute(), and GetMcInfo().

double HoughValidUpdate::tanl_mc [private]

Definition at line 111 of file HoughValidUpdate.h.

Referenced by execute().

int HoughValidUpdate::track_fit [private]

Definition at line 112 of file HoughValidUpdate.h.

Referenced by execute().

HoughValidUpdate Class Reference

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Classes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


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