BOSS_7.0.2: MdcHoughFinder Class Reference

MdcHoughFinder::MdcHoughFinder	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 56 of file MdcHoughFinder.cxx.

References m_binx, m_biny, m_combineTracking, m_data, m_debug, m_dropHot, m_fpro, m_getDigiFlag, m_hit_pro, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_method, m_minMdcDigi, m_ndev, m_pdtFile, m_peakCellNum, m_pickHits, m_pid, and m_trackNum_Mc_set.

00056                                                                               :
00057   Algorithm(name, pSvcLocator)
00058 {
00059   // Declare the properties  
00060   declareProperty("trackNumMc", m_trackNum_Mc_set=1);
00061   declareProperty("method", m_method=0);
00062   declareProperty("debug",          m_debug = 0);
00063   declareProperty("data",          m_data= 0);
00064   declareProperty("binx",          m_binx= 100);
00065   declareProperty("biny",          m_biny= 100);
00066   declareProperty("peakCellNum",          m_peakCellNum);
00067   declareProperty("ndev",          m_ndev= 3);
00068   declareProperty("f_pro",          m_fpro= 0.8);
00069   declareProperty("f_hit_pro",          m_hit_pro= 0.8);
00070   //declareProperty("debug",          m_debug = 0);
00071   declareProperty("pdtFile",        m_pdtFile = "pdt.table");
00072   declareProperty("pickHits",       m_pickHits = true);
00073   declareProperty("pid",        m_pid = 0);
00074   //declareProperty("helixHitsSigma", m_helixHitsSigma);
00075   declareProperty("combineTracking",m_combineTracking = true);
00076   declareProperty("getDigiFlag",    m_getDigiFlag = 0);
00077   declareProperty("maxMdcDigi",     m_maxMdcDigi = 0);
00078   declareProperty("keepBadTdc",     m_keepBadTdc = 0);
00079   declareProperty("dropHot",        m_dropHot= 0);
00080   declareProperty("keepUnmatch",    m_keepUnmatch = 0);
00081   declareProperty("minMdcDigi",     m_minMdcDigi = 0);
00082 
00083 }

StatusCode MdcHoughFinder::beginRun ( )

Definition at line 85 of file MdcHoughFinder.cxx.

References MdcDetector::instance(), and m_gm.

00085                                    {
00086 
00087   //Initailize MdcDetector
00088   m_gm = MdcDetector::instance(0);
00089   if(NULL == m_gm) return StatusCode::FAILURE;
00090 
00091 
00092   return StatusCode::SUCCESS;
00093 }

double MdcHoughFinder::CFtrans	(	double	x,
		double	y
	)			`[private]`

Definition at line 1667 of file MdcHoughFinder.cxx.

Referenced by execute().

01667                                                {
01668   return x/(x*x+y*y);
01669 }

int MdcHoughFinder::digiToHots	(	TrkRecoTrk *	newTrack,
		vector< bool >	vec_truthHit
	)			`[private]`

Definition at line 1606 of file MdcHoughFinder.cxx.

References TrkHitList::appendHot(), MdcHit::driftDist(), MdcHit::driftTime(), RawData::getChargeChannel(), RawDataProviderSvc::getMdcDigiVec(), MdcCalibFunSvc::getT0(), MdcCalibFunSvc::getTimeWalk(), TrkRecoTrk::hits(), RawData::identify(), MdcDetector::Layer(), MdcID::layer(), m_bunchT0, m_debug, m_gm, m_mdcCalibFunSvc, m_rawDataProviderSvc, TrkHitList::nHit(), MdcHit::rawTime(), MdcLayer::rMid(), MdcHit::setCalibSvc(), MdcHit::setCountPropTime(), TrkHitOnTrk::setFltLen(), subSeperate::temp, and MdcID::wire().

Referenced by execute().

01606                                                                             {
01607   TrkHitList* trkHitList = newTrack->hits();
01608   //HepPoint3D poca; Hep3Vector tempDir;
01609   //getInfo(helix,0,poca,tempDir);
01610 
01611   //if(m_debug>3)std::cout<<"track poca  "<<poca<<std::endl;
01612   // new MdcRecoHitOnTracks
01613   //SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
01614   //if (!mdcDigiCol) {
01615   //  cout << "Could not find MdcDigiCol" << endl;
01616   //  return 0;
01617   //}
01618 
01619   //if(mdcDigiCol->size()<5) return mdcDigiCol->size();
01620 
01621   //if(m_debug>0) std::cout<<" nDigi = "<<mdcDigiCol->size()<<std::endl;
01622   int digiId=0;
01623   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
01624   MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
01625   for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
01626         //      MdcDigiVec::iterator iter = mdcDigiCol->begin();
01627         //for (int i=0;iter != mdcDigiCol->end(); iter++,i++ ) {
01628         if(vec_truthHit[digiId]==false) continue;
01629         const MdcDigi* aDigi = *iter1;
01630         MdcHit *hit = new MdcHit(aDigi, m_gm);
01631         Identifier id = aDigi->identify();
01632         int layer = MdcID::layer(id);
01633         int wire  = MdcID::wire(id);
01634 
01635         if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}   
01636         //      if ((layer>=0&&layer<=7)||(layer>=20)) {continue;}      
01637         //if (layer<=35) {continue;}    
01638 
01639         if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
01640         hit->setCalibSvc(m_mdcCalibFunSvc);
01641         hit->setCountPropTime(true); //hit->setCountPropTime(m_countPropTime);
01642         //      if(m_debug>0)   cout<<"hitx:            "<<hit->x()<<"   "<<"hity:           "<<hit->y()<<"    "<<endl;
01643         //hit->setZTruth(t_mcZ[layer][wire]);
01644 
01645         // new fltLen and ambig
01646         int newAmbig = 0;
01647         // calc. position of this point
01648         //m_hitCol->push_back(hit);
01649         MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here
01650         MdcHitOnTrack* newhot = &temp;
01651         //      cout<<"ambig:    "<<newhot->ambig()<<"   xxxxxxx"<<endl;
01652         double fltLen = m_gm->Layer(layer)->rMid();
01653         newhot->setFltLen(fltLen);
01654         //newhot->setFltLen(0);
01655         //newhot->setHitRms(0.02);
01656         //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl;
01657         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())<<endl;
01658         //if(hit->driftTime(m_bunchT0,0)>800) continue; 
01659         if(hit->driftDist(m_bunchT0,0)>1.2) continue;
01660         trkHitList->appendHot(newhot);
01661         //}//end of loop hits
01662   }
01663   if(m_debug>1) std::cout<<std::endl;
01664   return trkHitList->nHit();
01665 }

int MdcHoughFinder::digiToHots2	(	TrkRecoTrk *	newTrack,
		vector< bool >	vec_truthHit
	)			`[private]`

Definition at line 2010 of file MdcHoughFinder.cxx.

References TrkHitList::appendHot(), MdcHit::driftDist(), MdcHit::driftTime(), RawData::getChargeChannel(), RawDataProviderSvc::getMdcDigiVec(), MdcCalibFunSvc::getT0(), MdcCalibFunSvc::getTimeWalk(), TrkRecoTrk::hits(), RawData::identify(), MdcDetector::Layer(), MdcID::layer(), m_bunchT0, m_debug, m_gm, m_mdcCalibFunSvc, m_rawDataProviderSvc, TrkHitList::nHit(), MdcHit::rawTime(), MdcLayer::rMid(), MdcHit::setCalibSvc(), MdcHit::setCountPropTime(), TrkHitOnTrk::setFltLen(), subSeperate::temp, and MdcID::wire().

Referenced by execute().

02010                                                                              {
02011   TrkHitList* trkHitList = newTrack->hits();
02012   //HepPoint3D poca; Hep3Vector tempDir;
02013   //getInfo(helix,0,poca,tempDir);
02014 
02015   //if(m_debug>3)std::cout<<"track poca  "<<poca<<std::endl;
02016   // new MdcRecoHitOnTracks
02017   //SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
02018   //if (!mdcDigiCol) {
02019   //  cout << "Could not find MdcDigiCol" << endl;
02020   //  return 0;
02021   //}
02022 
02023   //if(mdcDigiCol->size()<5) return mdcDigiCol->size();
02024 
02025   //if(m_debug>0) std::cout<<" nDigi = "<<mdcDigiCol->size()<<std::endl;
02026   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
02027   MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
02028   int digiId=0;
02029   for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
02030         //MdcDigiVec::iterator iter = mdcDigiCol->begin();
02031         //for (int i=0;iter != mdcDigiCol->end(); iter++,i++ ) {
02032         if(vec_truthHit[digiId]==false) continue;
02033         const MdcDigi* aDigi = *iter1;
02034         MdcHit *hit = new MdcHit(aDigi, m_gm);
02035         Identifier id = aDigi->identify();
02036         int layer = MdcID::layer(id);
02037         int wire  = MdcID::wire(id);
02038 
02039         //      if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}   
02040         //      if ((layer>=0&&layer<=7)||(layer>=20)) {continue;}      
02041         //if (layer<=35) {continue;}    
02042 
02043         if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
02044         hit->setCalibSvc(m_mdcCalibFunSvc);
02045         hit->setCountPropTime(true); //hit->setCountPropTime(m_countPropTime);
02046         //      if(m_debug>0)   cout<<"hitx:            "<<hit->x()<<"   "<<"hity:           "<<hit->y()<<"    "<<endl;
02047         //hit->setZTruth(t_mcZ[layer][wire]);
02048 
02049         // new fltLen and ambig
02050         int newAmbig = 0;
02051         // calc. position of this point
02052         //m_hitCol->push_back(hit);
02053         MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here
02054         MdcHitOnTrack* newhot = &temp;
02055         //      cout<<"ambig:    "<<newhot->ambig()<<"   xxxxxxx"<<endl;
02056         double fltLen = m_gm->Layer(layer)->rMid();
02057         newhot->setFltLen(fltLen);
02058         //newhot->setFltLen(0);
02059         //newhot->setHitRms(0.02);
02060         //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl;
02061         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())<<endl;
02062         //if(hit->driftTime(m_bunchT0,0)>800) continue; 
02063         if(hit->driftDist(m_bunchT0,0)>1.2) continue;
02064         trkHitList->appendHot(newhot);
02065         //}//end of loop hits
02066   }
02067   if(m_debug>1) std::cout<<std::endl;
02068   return trkHitList->nHit();
02069 }

StatusCode MdcHoughFinder::execute ( )

Definition at line 275 of file MdcHoughFinder.cxx.

References abs, MdcRawDataProvider::b_dropHot, MdcRawDataProvider::b_keepBadTdc, MdcRawDataProvider::b_keepUnmatch, MdcGeoWire::Backward(), TrkHotList::begin(), binX, binY, MdcGeoWire::Cell(), CFtrans(), cos(), count, TrkExchangePar::d0(), digiToHots(), digiToHots2(), dz_mc, TrkHotList::end(), showlog::err, TrkErrCode::failure(), Bes_Common::FATAL, FillCells(), FillHist(), FillRhoTheta(), TrkHitList::fit(), TrkRecoTrk::fitResult(), MdcGeoWire::Forward(), GetMcInfo(), RawDataProviderSvc::getMdcDigiVec(), TrkFit::helix(), TrkRecoTrk::hits(), TrkHitList::hotList(), genRecEmupikp::i, Bes_Common::INFO, ganga-rec::j, MdcGeoWire::Layer(), MdcID::layer(), LeastSquare(), Linefit(), m_2d_nFit, m_3d_nFit, m_binx, m_biny, m_bunchT0, m_cell, m_cell_Mc, m_combineTracking, m_context, m_cosCut, m_costaTruth, m_d0, m_data, m_debug, m_dropHot, m_eventNum, m_gm, m_hit_pro, m_hitSignal, m_keepBadTdc, m_keepUnmatch, m_l, m_layer, m_layer_Mc, m_layerNhit, m_mdcGeomSvc, m_method, m_ndev, m_nEvtSuccess, m_nFitFailure, m_nFitSucess, m_nHit, m_nHit_Mc, m_nHitAxial, m_nHitAxialSignal, m_nHitAxialSignal_select, m_nHitSignal, m_nHitSignal_select, m_omega, m_p, m_peakCellNum, m_phi0, M_PI, m_pickHits, m_pt, m_pt2, m_pxyz, m_pz, m_rawDataProviderSvc, m_runNum, m_slant, m_tanl, m_trackNum, m_trackNum_Mc, m_trackNum_Mc_set, m_u, m_v, m_x, m_x_east, m_x_west, m_y, m_y_east, m_y_west, m_z, m_z0, m_z_east, m_z_west, m_zStereo, m_zStereoNum, TrkSimpleMaker< T >::makeTrack(), msgSvc(), TrkFit::nActive(), MdcDetector::nLayer(), ntuplehit, TrkExchangePar::omega(), phi0, TrkRecoTrk::print(), push_back(), EventModel::Recon::RecMdcHitCol, EventModel::Recon::RecMdcTrackCol, RhoTheta(), TrkSimpleMaker< T >::setFlipAndDrop(), sin(), SlantId(), MdcTrack::storeTrack(), swap, t_eventNum, t_maxP, t_minP, t_runNum, tanl_mc, track_fit, v, Bes_Common::WARNING, MdcGeomSvc::Wire(), MdcID::wire(), x, and Zposition().

00275                                    {
00276 
00277   MsgStream log(msgSvc(), name());
00278   log << MSG::INFO << "in execute()" << endreq;
00279 
00280   //event start time
00281   SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00282   if (!evTimeCol) {
00283         log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
00284         m_bunchT0=0.;
00285   }else{
00286         RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00287         if (iter_evt != evTimeCol->end()){
00288           m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
00289         }
00290   }
00291 
00292   //if(m_debug)TrkHelixFitter::m_debug = true;
00293   // Get the event header, print out event and run number
00294 
00295   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00296   if (!eventHeader) {
00297         log << MSG::FATAL << "Could not find Event Header" << endreq;
00298         return( StatusCode::FAILURE);
00299   }
00300 
00301 
00302   DataObject *aTrackCol;
00303   DataObject *aRecHitCol;
00304   SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
00305   if(!m_combineTracking){
00306         eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00307         if(aTrackCol != NULL) {
00308           dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00309           eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00310         }
00311         eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00312         if(aRecHitCol != NULL) {
00313           dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00314           eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00315         }
00316   }
00317 
00318   RecMdcTrackCol* trackList;
00319   eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00320   if (aTrackCol) {
00321         trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00322   }else{
00323         trackList = new RecMdcTrackCol;
00324         StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00325         if(!sc.isSuccess()) {
00326           log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00327           return StatusCode::FAILURE;
00328         }
00329   }
00330   int nTdsTk = (int) trackList->size();
00331   int nFitSucess = 0;
00332 
00333   RecMdcHitCol* hitList;
00334   eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00335   if (aRecHitCol) {
00336         hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00337   }else{
00338         hitList = new RecMdcHitCol;
00339         StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00340         if(!sc.isSuccess()) {
00341           log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00342           return StatusCode::FAILURE;
00343         }
00344   }
00345 
00346   //pick the first layer
00347 
00348   // m_nLayerNum=43;
00349   // for(int i=0;i<43;i++){
00350   //     const MdcGeoWire *geowir_first_layer =m_mdcGeomSvc->Wire(i,1);
00351   //     HepPoint3D eastP = geowir_first_layer->Backward()/10.0;//mm 2 cm
00352   //     HepPoint3D westP = geowir_first_layer->Forward() /10.0;//mm 2 cm
00353   //     m_layer1X[i]=(eastP.x()+westP.x())/2;
00354   //     m_layer1Y[i]=(eastP.y()+westP.y())/2;
00355   //     m_layer1U[i]=CFtrans((eastP.x()+westP.x())/2.,(eastP.y()+westP.y())/2.);
00356   //     m_layer1V[i]=CFtrans((eastP.y()+westP.y())/2.,(eastP.x()+westP.x())/2.);
00357   //     
00358   //     cout<<"("<<i<<","<<"1"<<")"<<":"<<"("<<m_layer1X[i]<<","<<m_layer1Y[i]<<")"<<endl;
00359   // }
00360 
00361   t_eventNum=eventHeader->eventNumber();
00362   t_runNum=eventHeader->runNumber();
00363   m_eventNum=t_eventNum;
00364   m_runNum=t_runNum;
00365   log << MSG::INFO << "MdcHoughFinder: retrieved event: " << eventHeader->eventNumber()  << "  run: " <<  eventHeader->runNumber() << endreq;
00366 
00367   // get mc particle digi
00368   int mc_particle=GetMcInfo();
00369 
00370 
00371   if(abs(m_costaTruth)<=0.83){m_cosCut=1;}
00372   else{m_cosCut=0;}
00373 
00374   //int n_curve=0;
00375   //double vec_ptTruth=m_ptTruth;
00376   //double n_costaTruth=m_costaTruth;
00377   //double z_flight=2*M_PI*(771.4/2)*(vec_ptTruth/tan(n_costaTruth));
00378   //if(vec_ptTruth<0.12&&(z_flight<1193.)){
00381   //}
00382   vector<double> vec_u;
00383   vector<double> vec_v;
00384   vector<double> vec_p;
00385   vector<double> vec_x_east;
00386   vector<double> vec_y_east;
00387   vector<double> vec_z_east;
00388   vector<double> vec_x_west;
00389   vector<double> vec_y_west;
00390   vector<double> vec_z_west;
00391   vector<double> vec_z_Mc;
00392   vector<double> vec_x;
00393   vector<double> vec_y;
00394   vector<double> vec_z;
00395   vector<int> vec_layer;
00396   vector<int> vec_wire;
00397   vector<int> vec_slant;
00398   vector<double> vec_zStereo;
00399   vector<double> l;
00400 
00401   vector<int> vec_layer_Mc;
00402   vector<int> vec_wire_Mc;
00403   t_maxP=-999.;
00404   t_minP=999.;
00405 
00406   // retrieve Mdc truth hits
00407   int t_nHit_Mc;
00408   int digiId_Mc=0;
00409   int nHitAxial_Mc=0;
00410 
00411   SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol"); 
00412   if (!mcMdcMcHitCol) {
00413         log << MSG::INFO << "Could not find MdcMcHit" << endreq; 
00414   }else{
00415         Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin();
00416 
00417         for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) {
00418           Identifier mdcid = (*iterMcHit)->identify();
00419           int layerId_Mc = MdcID::layer(mdcid);
00420           int wireId_Mc = MdcID::wire(mdcid);
00421 
00422           vec_layer_Mc.push_back(layerId_Mc);
00423           vec_wire_Mc.push_back(wireId_Mc);
00424 
00425           //      cout<<"mcHit:  "<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<endl;
00426           //vec_layer.push_back(layerId_Mc);
00427           //vec_wire.push_back(wireId_Mc);
00428           m_layer_Mc[digiId_Mc]=layerId_Mc;
00429           m_cell_Mc[digiId_Mc]=wireId_Mc;
00430 
00431           if(m_data==0){
00432                 if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) {        
00433                   nHitAxial_Mc++;}
00434                 if(m_debug>0) {cout<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<"nHitAxial_Mc:  "<<nHitAxial_Mc<<endl;}
00435                 m_layer[digiId_Mc]=layerId_Mc;
00436                 m_cell[digiId_Mc]=wireId_Mc;
00437                 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc);
00438                 HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
00439                 HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
00440 
00441                 m_x_east[digiId_Mc]=eastP.x();
00442                 m_y_east[digiId_Mc]=eastP.y();
00443                 m_z_east[digiId_Mc]=eastP.z();
00444                 m_x_west[digiId_Mc]=westP.x();
00445                 m_y_west[digiId_Mc]=westP.y();
00446                 m_z_west[digiId_Mc]=westP.z();
00447 
00448                 vec_x_east.push_back(eastP.x());
00449                 vec_y_east.push_back(eastP.y());
00450                 vec_z_east.push_back(eastP.z());
00451                 vec_x_west.push_back(westP.x());
00452                 vec_y_west.push_back(westP.y());
00453                 vec_z_west.push_back(westP.z());
00454 
00455                 double x = (*iterMcHit)->getPositionX()/10.;//mm 2 cm
00456                 double y = (*iterMcHit)->getPositionY()/10.;//mm 2 cm
00457                 double z = (*iterMcHit)->getPositionZ()/10.;//mm 2 cm
00458 
00459                 double u=CFtrans(x,y);
00460                 double v=CFtrans(y,x);
00461                 double p=sqrt(u*u+v*v);
00462 
00463                 vec_x.push_back(x);
00464                 vec_y.push_back(y);
00465                 vec_z.push_back(z);
00466                 vec_u.push_back(u);
00467                 vec_v.push_back(v);
00468                 vec_p.push_back(p);
00469 
00470                 m_x[digiId_Mc]=x;
00471                 m_y[digiId_Mc]=y;
00472                 m_z[digiId_Mc]=z;
00473                 m_u[digiId_Mc]=u;
00474                 m_v[digiId_Mc]=v;
00475                 m_p[digiId_Mc]=p;
00476 
00477                 int layer= layerId_Mc;
00478                 vec_slant.push_back(SlantId(layer));
00479                 m_slant[digiId_Mc]=SlantId(layer);
00480                 if (m_slant!=0)
00481                 {cout<<"layer:  "<<layerId_Mc<<"wire:  "<<wireId_Mc<<"x_east:  "<<m_x_east[digiId_Mc]<<"y_east:  "<<m_y_east[digiId_Mc]<<endl;}
00482           }
00483         }
00484   }
00485 
00486   t_nHit_Mc=digiId_Mc;
00487   m_nHit_Mc=digiId_Mc;
00488   if(m_data==0) {m_nHit=digiId_Mc;}
00489 
00490   //m_nHitAxial=nHitAxial_Mc;
00491   //    m_nFitFailure=0;
00492 
00493   //    for(int i =0;i<t_nHit_Mc;i++){
00494   //      if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];} 
00495   //    }
00496   //cout<<"t_maxP:"<<t_maxP<<"  "<<endl;
00497 
00498   //    for(int i =0;i<t_nHit_Mc;i++){
00499   //      if(t_minP>vec_p[i]) {t_minP=vec_p[i];} 
00500   //    }
00501   //cout<<"t_minP:"<<t_minP<<"  "<<endl;
00502 
00503 
00504   //  retrieve Mdc digi vector form RawDataProviderSvc 
00505   vector<int> vec_hitSignal;
00506   vector<int> vec_track_index;
00507   int track_index_max=0;
00508   uint32_t getDigiFlag = 0;
00509   if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
00510   if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00511   if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00512   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00513   //  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
00514   cout<<"MdcDigiVec:  "<<mdcDigiVec.size()<<endl;
00515   MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00516   int t_nHit;
00517   int digiId = 0;
00518   int nHitAxial = 0;
00519   int nHitLayer[43];
00520   int nHitSignal=0;
00521   int nHitAxialSignal=0;
00522 
00523   for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00524         int track_index=(*iter1)->getTrackIndex();
00525 //      cout<<"mc track_index: "<<track_index<<endl;
00526         if(track_index>=1000) track_index-=1000;
00527         vec_track_index.push_back(track_index);
00528         if(track_index>=0){
00529           nHitSignal++;
00530           m_hitSignal[digiId]=1;
00531         }
00532         //m_hitCol[digiId]=digiId; 
00533         Identifier mdcId= (*iter1)->identify();
00534         int     layerId = MdcID::layer(mdcId);
00535         int     wireId = MdcID::wire(mdcId);
00536 
00537 
00538         //cout<<"layer:  "<<layerId<<"  "<<"wire:  "<<wireId<<"hitSignal:  "<<vec_hitSignal[digiId]<<endl;
00539 
00540         if ((layerId>=8&&layerId<=19)||(layerId>=36)) { 
00541           nHitAxial++;}
00542         if (((layerId>=8&&layerId<=19)||(layerId>=36))&&track_index>=0) {       
00543           nHitAxialSignal++;}
00544 
00545         vec_layer.push_back(layerId);
00546         vec_wire.push_back(wireId);
00547 
00548         nHitLayer[layerId]++;
00549         m_layerNhit[layerId]=nHitLayer[layerId];        
00550 
00551         const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId);
00552         HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
00553         HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
00554 
00555         vec_x_east.push_back(eastP.x());
00556         vec_y_east.push_back(eastP.y());
00557         vec_z_east.push_back(eastP.z());
00558         vec_x_west.push_back(westP.x());
00559         vec_y_west.push_back(westP.y());
00560         vec_z_west.push_back(westP.z());
00561 
00562         m_x_east[digiId]=eastP.x();
00563         m_y_east[digiId]=eastP.y();
00564         m_z_east[digiId]=eastP.z();
00565         m_x_west[digiId]=westP.x();
00566         m_y_west[digiId]=westP.y();
00567         m_z_west[digiId]=westP.z();
00568         if(m_debug>0) {cout<<"event:  "<<t_eventNum<<"  "<<"layer:  "<<layerId<<" "<<"wireId:  "<<wireId<<" "<<"zeast:  "<<vec_z_east.at(digiId)<<"  "<<"zwest:  "<<vec_z_west.at(digiId)<<"  "<<endl;}
00569         //conformal trans:
00570         double x =(eastP.x()+westP.x())/2.;
00571         double y =(eastP.y()+westP.y())/2.;
00572         vec_x.push_back((vec_x_east[digiId]+vec_x_west[digiId])/2); 
00573         vec_y.push_back((vec_y_east[digiId]+vec_y_west[digiId])/2); 
00574         m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2; 
00575         m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2; 
00576         double u=CFtrans(x,y);
00577         double v=CFtrans(y,x);
00578         //      cout<< "u:   "<<u<<"v:   "<<v<<"digiId:  "<<digiId<<endl;
00579         //cout<<"layer:  "<<layerId<<" "<<"wireId:  "<<wireId<<" "<<"x:      "<<x<<"   "<<"y:     "<<y<<"     "<<endl;
00580         vec_p.push_back(sqrt(u*u+v*v));
00581         vec_u.push_back(u);
00582         vec_v.push_back(v);
00583         m_u[digiId]=u;
00584         m_v[digiId]=v;
00585         m_p[digiId]=sqrt(u*u+v*v);
00586 
00587         m_layer[digiId]=geowir->Layer();
00588         m_cell[digiId]=geowir->Cell();
00589         int layer= layerId;
00590         vec_slant.push_back(SlantId(layer));
00591         m_slant[digiId]=SlantId(layer);
00592         //cout<<"layer:  "<<layer<<"  "<<"slant:  "<<m_slant[digiId]<<endl;
00593         //   if (m_slant!=0)
00594         //      {cout<<"layer:  "<<layerId<<"wire:  "<<wireId<<"x_east:  "<<m_x_east[digiId]<<"y_east:  "<<m_y_east[digiId]<<endl;}
00595   }
00596   for(int i=0;i<digiId;i++){
00597         if(track_index_max<=vec_track_index[i]) track_index_max=vec_track_index[i];
00598   }
00599   track_index_max++;
00600   m_trackNum_Mc=track_index_max;
00601 //  cout<<"mc truth has :"<<track_index_max<<"  mc track."<<endl;
00602   vector< vector <int> > mc_track_hit(track_index_max,vector<int>() );
00603   for(int i=0;i<track_index_max;i++){ 
00604         for(int j=0;j<digiId;j++){
00605           if(vec_track_index[j]==i) mc_track_hit[i].push_back(j);
00606         }
00607   }
00608   //composory set mc_track number
00609   track_index_max=m_trackNum_Mc_set;
00610   //debug mc track hit
00611   //for(int i=0;i<track_index_max;i++){ 
00612   //  for(int j=0;j<mc_track_hit[i].size();j++){
00613   //    int hitId=mc_track_hit[i][j];
00614   //    cout<<"track: "<<i<<"  has hit:  "<<mc_track_hit[i][j]<<"("<<vec_layer[hitId]<<","<<vec_wire[hitId]<<")"<<endl;
00615   //  }
00616   //}
00617 
00618   //t_nHit=digiId;
00619   m_nHit=digiId;
00620   m_nHitAxial=nHitAxial;
00621   //m_nFitFailure=0;
00622 
00623   m_nHitSignal=nHitSignal;
00624   m_nHitAxialSignal=nHitAxialSignal;
00625 
00626   cout<<"hit number:  "<<digiId<<endl;
00627 
00628 
00629   // get the max&min amplitude of the curve 
00630   for(int i =0;i<m_nHit;i++){
00631         if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];} 
00632   }
00633 
00634   for(int i =0;i<m_nHit;i++){
00635         if(t_minP>vec_p[i]) {t_minP=vec_p[i];} 
00636   }
00637   t_maxP=t_maxP+0.01;
00638 
00639 
00640   //double peakArea[100];
00641   //m_npeak=100;
00642   //for( int bin_i=0;bin_i<10;bin_i++){
00643   //    for(int bin_j=0;bin_j<10;bin_j++){
00644   //      int binx=m_binx+100*bin_i;
00645   //      int biny=m_biny+100*bin_j;
00646   int nhit;
00647   if(m_data==0) {
00648         nhit=m_nHit_Mc;}
00649   else{
00650         nhit=m_nHit;}
00651 
00652   binX=m_binx;
00653   binY=m_biny;
00654   double binwidth=M_PI/binX;
00655   double binhigh=2*t_maxP/binY;
00656   TH2D *h1=new TH2D("h1","h1",binX,0,M_PI,binY,-t_maxP,t_maxP);
00657 
00658   //method 0
00659   //2-point combine
00660   vector<double> vec_rho;
00661   vector<double> vec_theta;
00662   vector< vector<int> > vec_hitNum(2,vector<int>());      //store 2 hits in a cross point
00663   int numCross=(int)(nhit*(nhit-1)/2);
00664   if(m_method==0){
00665         RhoTheta(numCross,nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum);    //calculate cross point
00666         FillRhoTheta(h1,vec_theta,vec_rho,numCross);         //fill histgram method by rhotheta
00667   }
00668   // method  1
00669   vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) ); 
00670   if(m_method==1){
00671         FillHist(h1,vec_u,vec_v,nhit,vec_selectNum);
00672   }
00673 
00674 
00675   vector<bool> vec_truthHit(nhit,false);
00676   vector< vector <int> > countij(binX,vector <int> (binY,0) ); 
00677   //  vector< vector < vector<int> > > vec_selectNum(binX,vector < vector<int> >(binY,vector<int>() )); 
00678   vector< vector< vector<int> > > vec_selectHit(binX,vector< vector<int> >(binY,vector<int>() ) ); 
00679   if(m_method==2){
00680         FillCells(h1,nhit,binX,binY,vec_u,vec_v,vec_layer,vec_wire,countij,vec_selectNum,vec_selectHit);
00681   }
00682 
00683   //cout h1 
00684   //for(int i=0;i<binX;i++){
00685   //  for(int j=0;j<binY;j++){
00686   //    int count =h1->GetBinContent(i+1,j+1);
00688   //    cout<<"("<<i<<","<<j<<")"<<":  "<<count<<endl;
00690   //  }
00691   //}
00692 
00693 
00694   //int nHitSelect=0;
00695   //int nHitSignal_select=0;
00696   //int nHitAxialSignal_select=0;
00697   //find peak
00698   int max_count=0;
00699   int max_binx=1;
00700   int max_biny=1;
00701   for(int i=0;i<binX;i++){
00702         for(int j=0;j<binY;j++){
00703           int count=h1->GetBinContent(i+1,j+1);
00704           if(max_count<count) {
00705                 max_count=count;
00706                 max_binx=i+1;
00707                 max_biny=j+1;
00708           }
00709         }
00710   }
00711   //            cout<<"("<<max_binx<<","<<max_biny<<","<<max_count<<")"<<"numCross:  "<<numCross<<endl;
00712   //  h1->Draw("lego2");
00713 
00714   //int peak_cellNum=1;         // to determine how many cell should be embraced in one orientation
00715   //int count_sum[1000];
00716   //for(int i=0;i<1000;i++){
00717   //  count_sum[i]=0;
00718   //  for(int alphax=max_binx-i;alphax<=max_binx+i;alphax++){
00719   //    for(int rhox=max_biny-i;rhox<=max_biny+i;rhox++){
00720   //    int ax;
00721   //    if (alphax<0) { ax=alphax+binX; }
00722   //    else if (alphax>=binX) { ax=alphax-binX; }
00723   //    else { ax=alphax; }
00724   //    count_sum[i]+=h1->GetBinContent(ax,rhox);
00725   //    }
00726   //  }
00727   //  if(count_sum[i]>m_fpro*numCross){
00728   //    peak_cellNum=i;
00729   //                      cout<<"i: "<<i<<endl;
00730   //                    cout<<"count_sum:  "<<count_sum[i]<<endl;
00731   //                    cout<<"pro: "<<(double)count_sum[i]/(double)numCross<<endl;
00732   //    break;
00733   //  }
00734   //}
00735 
00736 
00737   /*
00738   //return hit number
00739   vector<bool> vec_hitSelect(m_nHit,false);
00740   int bin_left=max_binx-peak_cellNum;
00741   int bin_right=max_binx+peak_cellNum;
00742   int bin_up=max_biny+peak_cellNum;
00743   int bin_down=max_biny-peak_cellNum;
00744   double area_left=(bin_left-1)*binwidth;
00745   double area_right=(bin_right)*binwidth;
00746   double area_down=(bin_down-1)*binhigh-t_maxP;
00747   double area_up=(bin_up)*binhigh-t_maxP;
00748   // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
00749   double peak_width=(2*peak_cellNum+1)*binwidth;
00750   double peak_high=(2*peak_cellNum+1)*binhigh;
00751   cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
00752   //test if the hit is in this area
00753   for(int i=0;i<numCross;i++){
00754   if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
00755   int hitNum_1=vec_hitNum[0][i];
00756   int hitNum_2=vec_hitNum[1][i];
00757   vec_hitSelect[hitNum_1]=true;
00758   vec_hitSelect[hitNum_2]=true;
00759   }
00760   }
00761   int selectHitNum=0;
00762   for(int i=0;i<m_nHit;i++){
00763   if(vec_hitSelect[i]==1) {selectHitNum++;}
00764   //    cout<<"if hit is select: "<<i<<":  "<<vec_hitSelect[i]<<endl;
00765   }
00766   cout<<"how many Hit have been selected: "<<selectHitNum<<"  "<<(double)selectHitNum/(double)m_nHit;
00767   m_areaSelectHit=(double)selectHitNum/(double)m_nHit;
00768 
00769 
00770   //    // the smallest cell has been found
00771   //      //cout<<"width: "<<peak_width<<"  "<<"high: "<<peak_high<<endl;
00772   //      int column=(int)(bin_i*10+bin_j);
00773   //      peakArea[column]=peak_width*peak_high;
00774   //      m_peakWidth[column]=peak_width;
00775   //      m_peakHigh[column]=peak_high;
00777   //      int columnx_split=(int)(column/10)*100+100;
00778   //      int columny_split=(int)(column%10)*100+100;
00779   //      //cout<<"column: "<<column<<"  "<<"("<<columnx_split<<","<<columny_split<<")"<<peakArea[column]<<endl;
00780   delete h1;
00781   //    }
00782   //  }
00783   //  double area_least=peakArea[0];
00784   //  for(int i=0;i<100;i++){
00785   //    if(area_least>peakArea[i]){
00786   //      area_least=peakArea[i];
00787   //    }
00788   //  }
00789   //  int num_area_least=0;
00790   //  for(int i=0;i<100;i++){
00791   //    if(area_least==peakArea[i]){
00792   //      num_area_least=i;
00793   //    }
00794   //  }
00795   //  int binx_split=(int)(num_area_least/10)*100+100;
00796   //  int biny_split=(int)(num_area_least%10)*100+100;
00797   //  m_areaLeast=area_least;
00798   //  m_areaLeastNum=num_area_least;
00799   //cout<<"the max peakArea: "<<num_area_least<<"("<<binx_split<<","<<biny_split<<")"<<"peakarea is : "<<area_least<<endl;
00800    */
00801 
00802 
00803 
00804 
00805   // select hit
00806   int count_use=0;
00807   double sum=0;
00808   double sum2=0;        
00809   for(int i=0;i<binX;i++){
00810         for(int j=0;j<binY;j++){
00811           int count=h1->GetBinContent(i+1,j+1);
00812           //            cout<<"("<<i<<","<<j<<")"<<"  "<<"count:  "<<count<<"          ";
00813           if(count!=0){
00814                 count_use++;
00815                 sum+=count;
00816                 sum2+=count*count;
00817           }
00818         }
00819   }
00820   double f_n=m_ndev;
00821   cout<<"count_use:  "<<count_use<<endl;
00822   double aver=sum/count_use;
00823   double dev=sqrt(sum2/count_use-aver*aver);
00824   int min_counts=(int)(aver+f_n*dev);
00825   cout<<"aver:  "<<aver<<"  "<<"dev:  "<<dev<<"  min: "<<min_counts<<endl;
00826 
00827 
00828   vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) ); 
00829   vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) ); 
00830   //  int hough_trans_CS[binX][binY];
00831   // int hough_trans_CS_peak[binX][binY];
00832   for(int i=0;i<binX;i++){
00833         for(int j=0;j<binY;j++){
00834           hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1);
00835         }
00836   }
00837 
00838   int peak_num=0;
00839   for (int r=0; r<binY; r++) {
00840         for (int a=0; a<binX; a++) {
00841           long max_around=0;            
00842           for (int ar=a-1; ar<=a+1; ar++) {
00843                 for (int rx=r-1; rx<=r+1; rx++) {
00844                   int ax;
00845                   if (ar<0) { ax=ar+binX; }
00846                   else if (ar>=binX) { ax=ar-binX; }
00847                   else { ax=ar; }
00848                   if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00849                         if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; }
00850                   }
00851                 }
00852           }
00853           if (hough_trans_CS[a][r]<=min_counts) { hough_trans_CS_peak[a][r]=0; }
00854           else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; }
00855           else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; }
00856           else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; }
00857           if(hough_trans_CS_peak[a][r]!=0) { 
00858                 //  cout<<" peak:  "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<"  "<<hough_trans_CS[a][r]<<endl;
00859                 peak_num++;
00860           }
00861         }
00862   }
00863   //cout<<"peak_num:   "<<peak_num<<endl;
00864 
00865   //  //find double-point-peaks in parameter space
00866   int peak_num_2=0;
00867   for (int r=0; r<binY; r++) {
00868         for (int a=0; a<binX; a++) {
00869           if (hough_trans_CS_peak[a][r]==1) {
00870                 hough_trans_CS_peak[a][r]=2;
00871                 for (int ar=a-1; ar<=a+1; ar++) {
00872                   for (int rx=r-1; rx<=r+1; rx++) {
00873                         int ax;
00874                         if (ar<0) { ax=ar+binX; }
00875                         else if (ar>=binX) { ax=ar-binX; }
00876                         else { ax=ar; }
00877                         if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00878                           if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; }
00879                         }
00880                   }
00881                 }
00882           }
00883           if(hough_trans_CS_peak[a][r]!=0) {
00884                 peak_num_2++;
00885                 //cout<<" peak:  "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<"  "<<hough_trans_CS[a][r]<<endl;
00886                 //              for(int i=0;i<hough_trans_CS[a][r];i++){
00887                 //                int hitNum=selectNum[a][r][i];
00888                 //                cout<<"hit:  "<<hitNum<<"("<<layer[hitNum]<<","<<cell[hitNum]<<")"<<endl;
00889                 //              }
00890           }
00891         }
00892   }
00893   //cout<<"peak_mum_2:   "<<peak_num_2<<endl;
00894 
00895   //fill the histogram again
00896   TH2D *h2 = new TH2D("h2","test2",binX,0,3.14,binY,-t_maxP,t_maxP);
00897   for(int i=0;i<binX;i++){
00898         for(int j=0;j<binY;j++){
00899           if(hough_trans_CS_peak[i][j]==2){
00900                 h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]);
00901           }
00902         }
00903   }
00904 
00905   vector<int> peakList1[3];
00906   for(int n=0;n<peak_num_2;n++){
00907         for (int r=0; r<binY; r++) {
00908           for (int a=0; a<binX; a++) {
00909                 if (hough_trans_CS_peak[a][r]==2) {
00910                   peakList1[0].push_back(a+1);
00911                   peakList1[1].push_back(r+1);
00912                   peakList1[2].push_back(hough_trans_CS[a][r]);
00913                 }
00914           }
00915         }
00916   }
00917   cout<<"finish peak?"<<endl;
00918   if(m_method==0) cout<<"numCross:  "<<numCross<<endl;
00919   //  //sort peaks by height
00920   int npeak1=peak_num_2;
00921   int n_max;
00922   for (int na=0; na<npeak1-1; na++) {
00923         n_max=na;
00924         for (int nb=na+1; nb<npeak1; nb++) {
00925           if (peakList1[2][n_max]<peakList1[2][nb]) { n_max=nb; }
00926         }   
00927         if (n_max!=na) {        // swap
00928           float swap[3];
00929           for (int i=0; i<3; i++) {
00930                 swap[i]=peakList1[i][na];
00931                 peakList1[i][na]=peakList1[i][n_max];
00932                 peakList1[i][n_max]=swap[i];
00933           }
00934         }
00935   }
00936   cout<<"npeak1:  "<<npeak1<<endl;
00937   for(int n=0;n<npeak1;n++){
00938         int bintheta=peakList1[0][n];
00939         int binrho=peakList1[1][n];
00940         int count=peakList1[2][n];
00941         for(int i=0;i<count;i++){
00942           //cout<<"("<<bintheta<<","<<binrho<<","<<count<<")"<<endl;
00943           //cout<<"n"<<n<<":  "<<"("<<bintheta<<","<<binrho<<","<<count<<"):"<<vec_selectNum[bintheta-1][binrho-1][i]<<endl;
00944         }
00945         //cout<<"peakList1:  "<<n<<"  "<<"bina:  "<<peakList1[0][n]<<"  "<<"binr:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl; 
00946   }
00947 
00948   typedef std::map<int, int> M2;
00949   typedef std::map<int, M2> M1;
00950   M2 peak_combine_num;
00951   M1 hit_combine;
00952 
00953   int peak_track=0;
00954 
00955   if(m_method==2){
00956         //combine aroud
00957         int m=0;
00958         int n=0;
00959         int addnum=999;
00960         vector<bool> peaktrack(npeak1,true);
00961         while(addnum!=0)
00962         {
00963           addnum=0;
00964           double peak_cellNum[43];
00965           int peakX[43];
00966           int peakY[43];
00967           double  bin_left[43];
00968           double  bin_right[43];
00969           double  bin_up[43];
00970           double  bin_down[43];
00971           double area_left[43];
00972           double area_right[43];
00973           //double area_mid[43];
00974           double area_down[43];
00975           double area_up[43];
00976           // for(int i=0;i<npeak1;i++)  { peaktrack[i]=true; }
00977           for(int iLayer=0; iLayer<m_gm->nLayer(); iLayer++){
00978                 peak_cellNum[iLayer]=m_peakCellNum[iLayer];
00979                 peakX[iLayer]=peakList1[0][n];        //start from the highest peak
00980                 peakY[iLayer]=peakList1[1][n];
00981                 bin_left[iLayer]=peakX[iLayer]-peak_cellNum[iLayer];
00982                 bin_right[iLayer]=peakX[iLayer]+peak_cellNum[iLayer];
00983                 bin_up[iLayer]=peakY[iLayer]+peak_cellNum[iLayer];
00984                 bin_down[iLayer]=peakY[iLayer]-peak_cellNum[iLayer];
00985                 area_left[iLayer]=(bin_left[iLayer]-1)*binwidth;
00986                 area_right[iLayer]=(bin_right[iLayer])*binwidth;
00987                 // area_mid=(area_left+area_right)/2.;
00988                 area_down[iLayer]=(bin_down[iLayer]-1)*binhigh-t_maxP;
00989                 area_up[iLayer]=(bin_up[iLayer])*binhigh-t_maxP;
00990           }
00991           int count_peak=0;
00992           for(int k=0;k<nhit;k++){
00993                 int layer=vec_layer[k];
00994                 double lineLeft=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
00995                 double lineRight=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
00996                 // double lineMid=vec_u[k]*cos(area_mid)+vec_v[k]*sin(area_mid);
00997                 if((lineLeft<area_up[layer])&&(lineLeft>area_down[layer])||(lineRight<area_up[layer])&&(lineRight>area_down[layer])||((lineLeft-area_up[layer])*(area_down[layer]-lineRight)>0)){ 
00998                   //      hit_combine[m].push_back(k);
00999                   //            if(m==0){
01000                   //                    nHitSelect++;
01001                   //              if(vec_track_index[k]>=0){nHitSignal_select++;}
01002                   //              if(vec_track_index[k]>=0&&vec_slant[k]==0){nHitAxialSignal_select++;}
01003                   //            }
01004                   hit_combine[m][count_peak]=k;
01005                   //vec_truthHit[m][k]=1;
01006                   count_peak++;
01007                 }
01008                 //else {vec_truthHit[m][k]=0;}
01009                 peak_combine_num[m]=count_peak;
01010                 //int sizepeak=hit_combine[0].size();
01011                 //cout<<"peak1size:  "<<sizepeak<<endl;
01012           }
01013           for(int i=n+1;i<npeak1;i++){
01014                 if(peaktrack[i]==false) continue;
01015                 int peaktheta=peakList1[0][i];
01016                 int peakrho=peakList1[1][i];
01017                 int peakhitNum=peakList1[2][i];
01018                 int count_same_hit=0;
01019                 for(int j=0;j<peakhitNum;j++){
01020                   //for(int k=0;k<hit_combine[m].size();k++)
01021                   for(int k=0;k<peak_combine_num[m];k++){
01022                         if(vec_selectNum[peaktheta-1][peakrho-1][j]==hit_combine[m][k]){
01023                           count_same_hit++;
01024                         }
01025                   }
01026                 }
01027                 double f_hit=m_hit_pro;
01028                 if(count_same_hit>f_hit*peakhitNum){
01029                   peaktrack[i]=false;
01030                 }
01031           }
01032           for(int i=n+1;i<npeak1;i++){
01033                 if(peaktrack[i]==true){
01034                   addnum=i;
01035                   break;
01036                 }
01037           }
01038           if(addnum!=0) m++; 
01039           cout<<"peak_m:  "<<m+1<<endl;
01040           n=n+addnum;
01041         }
01042 
01043 
01044         for(int i=0;i<npeak1;i++){
01045           cout<<"track"<<i<<":  "<<"("<<peakList1[0][i]<<","<<peakList1[1][i]<<","<<peakList1[2][i]<<")"<<"  "<<"truth:  "<<peaktrack[i]<<endl;
01046           if(peaktrack[i]==true){
01047                 peak_track++;
01048           }
01049         }
01050         for( int i=0;i<peak_track;i++){
01051           for(int j =0;j<peak_combine_num[i];j++){
01052                 int hit_number=hit_combine[i][j];
01053                 cout<<" peak "<<i<<"  has select hits: "<<vec_layer[hit_number]<<"   "<<vec_wire[hit_number]<<endl;
01054           }
01055           cout<<"has collect :"<<peak_combine_num[i]<<" hits "<<endl;
01056         }
01057         cout<<"event"<<t_eventNum<<":  "<<"has found: "<<peak_track<<" track."<<endl;
01058         m_trackNum=peak_track;
01059         //m_nHitSelect=nHitSelect;
01060         //m_nHitSignal_select=nHitSignal_select;
01061         //m_nHitAxialSignal_select=nHitAxialSignal_select;
01062 
01063         //according mc and track
01064         vector< vector<int> > rec_mc_num(peak_track,vector<int>(track_index_max,0) );
01065         vector< vector<int> > rec_mc_num_axial(peak_track,vector<int>(track_index_max,0) );
01066         if(track_index_max!=peak_track)  cout<<"not match track number !"<<endl;
01067         for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
01068           for(int i=0;i<peak_track;i++){
01069                 for(int j=0;j<peak_combine_num[i];j++){
01070                   for(int k=0;k<mc_track_hit[mc_track_num].size();k++){
01071                         if(hit_combine[i][j]==mc_track_hit[mc_track_num][k]){
01072                           rec_mc_num[i][mc_track_num]++; 
01073                           int hit_num=mc_track_hit[mc_track_num][k];
01074                           if(vec_slant[hit_num]==0) rec_mc_num_axial[i][mc_track_num]++; 
01075                         }
01076                   }
01077                 }
01078           }
01079         }
01080         vector<int> rec_mc(peak_track,999);
01081         for(int i=0;i<peak_track;i++){
01082           for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
01083                 if(rec_mc_num[i][mc_track_num]>0.5*peak_combine_num[i])  rec_mc[i]=mc_track_num; 
01084                 cout<<"trackrec: "<<i<<" trackmc:  "<<mc_track_num<<"   "<<rec_mc_num[i][mc_track_num]<<"  "<<peak_combine_num[i]<<endl;
01085           }
01086         }
01087         for(int i=0;i<peak_track;i++){
01088           cout<<"rec :"<<i<<"belong to mc track: "<<rec_mc[i]<<endl;
01089         }
01090         for(int i=0;i<peak_track;i++){
01091           int mc_track_num=rec_mc[i];
01092           if(mc_track_num!=999) {
01093                 cout<<"debug mc_track_num: "<<mc_track_num<<endl;
01094                 m_nHitSignal_select[i]=rec_mc_num[i][mc_track_num];          //????????? mc_track_num=999?
01095                 m_nHitAxialSignal_select[i]=rec_mc_num_axial[i][mc_track_num];          //????????? mc_track_num=999?
01096                 cout<<"m_nHitSignal_select:  "<<m_nHitSignal_select[i]<<endl;
01097           }
01098         }
01099 
01100         //cout<<"peak1: "<<endl;
01101         //for(int i=0;i<peak_track;i++){
01102         //  for(int j=0;j<peak_combine_num[i];j++){
01103         //      int hitNumtemp=hit_combine[i][j];
01104         //      vec_truthHit[hitNumtemp]=i;
01105         //      cout<<"         "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<":  "<<vec_track_index[hitNumtemp]<<endl;
01106         //  }
01107         //}
01108 
01109         //      if(peak_track==2){
01110         //        cout<<"peak2: "<<endl;
01111         //        for(int i=0;i<hit_combine[1].size();i++){
01112         //              int hitNumtemp=hit_combine[1][i];
01113         //              vec_truthHit[hitNumtemp]=1;
01114         //              cout<<"         "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<":  "<<vec_track_index[hitNumtemp]<<endl;
01115         //        }
01116         //      }
01117 
01118 
01119   }
01120 
01121 
01122   //yzhang 2015-03-03 delete
01123 //  //return hit number
01124 //  if(m_method==0){
01125 //      vector<int> numCross_select;
01126 //      vector< vector<bool> > vec_hitSelect(npeak1,vector<bool> (nhit,false) );
01127 //      vector<int> vec_selectHitNum;
01128 //      vector<int> vec_selectHitNum_signal;
01129 //      for(int peakNum=0;peakNum<npeak1;peakNum++){
01130 //        int peak_cellNum=m_peakCellNum;
01131 //        //    vector<bool> vec_hitSelect(nhit,false);
01132 //        int peakX=peakList1[0][peakNum];
01133 //        int peakY=peakList1[1][peakNum];
01134 //        int bin_left=peakX-peak_cellNum;
01135 //        int bin_right=peakX+peak_cellNum;
01136 //        int bin_up=peakY+peak_cellNum;
01137 //        int bin_down=peakY-peak_cellNum;
01138 //        double area_left=(bin_left-1)*binwidth;
01139 //        double area_right=(bin_right)*binwidth;
01140 //        double area_down=(bin_down-1)*binhigh-t_maxP;
01141 //        double area_up=(bin_up)*binhigh-t_maxP;
01142 //        // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
01143 //        double peak_width=(2*peak_cellNum+1)*binwidth;
01144 //        double peak_high=(2*peak_cellNum+1)*binhigh;
01145 //        //cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
01146 //        //test if the hit is in this area
01147 //        int numCross_select_temp=0;
01148 //        for(int i=0;i<numCross;i++){
01149 //              if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
01150 //                numCross_select_temp++;
01151 //                int hitNum_1=vec_hitNum[0][i];
01152 //                int hitNum_2=vec_hitNum[1][i];
01153 //                vec_hitSelect[peakNum][hitNum_1]=true;
01154 //                vec_hitSelect[peakNum][hitNum_2]=true;
01155 //              }
01156 //        }
01157 //        numCross_select.push_back(numCross_select_temp);
01158 //        int selectHitNum=0;
01159 //        int selectHitNum_signal=0;
01160 //        for(int i=0;i<nhit;i++){
01161 //              if(vec_hitSelect[peakNum][i]==true) {
01162 //                selectHitNum++;
01163 //                if(vec_hitSignal[i]==1){
01164 //                      selectHitNum_signal++;
01165 //                }
01166 //                //    cout<<"if hit is select: "<<i<<":  "<<vec_hitSelect[i]<<endl;
01167 //              }
01168 //        }
01169 //        vec_selectHitNum.push_back(selectHitNum);
01170 //        vec_selectHitNum_signal.push_back(selectHitNum_signal);
01171 //        //  cout<<"hit: "<<selectHitNum<<"  "<<(double)selectHitNum/(double)nhit<<endl;
01172 //        //cout<<"true hit : "<<selectHitNum_signal<<"  "<<(double)selectHitNum_signal/(double)nhit<<endl;
01173 //        m_areaSelectHit=(double)selectHitNum/(double)nhit;
01174 //        m_areaSelectHit_signal=(double)selectHitNum_signal/(double)nhit;
01175 //      }
01176 //
01177 //      for(int i=0;i<npeak1;i++){
01178 //        int selectHitNum=0;
01179 //        for(int j=0;j<nhit;j++){
01180 //              if(vec_hitSelect[i][j]==true){
01181 //                selectHitNum++;
01182 //                //  cout<<"event"<<t_eventNum<<": "<<"peak"<<i<<": "<<"hit: "<<"("<<vec_layer[j]<<","<<vec_wire[j]<<")"<<endl;    //cout picked hit
01183 //              }
01184 //        }
01185 //        cout<<"peak"<<i<<": "<<selectHitNum<<"  hits"<<endl;
01186 //      }
01187 //      //  for(int i=0;i<3;i++){
01188 //      //      int bina=peakList1[0][i];
01189 //      //      int binr=peakList1[1][i];
01190 //      //      int hit_peak=peakList1[2][i];
01191 //      //      for(int j=0;j<hit_peak;j++){
01192 //      //        int hitList=vec_selectHit[bina][binr][j];
01193 //      //        cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
01194 //      //        //      cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
01195 //      //      }
01196 //      //  }
01197 //      //cout<<npeak1<<endl;
01198 //
01199 //      for(int n=0;n<npeak1;n++){
01200 //        cout<<"peakList1:  "<<n<<"  "<<"bina:  "<<peakList1[0][n]<<"  "<<"binr:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<"  "<<"numCross_select:"<<numCross_select[n]<<"  "<<" cross:  "<<(double)numCross_select[n]/(double)numCross<<"  "<<"hit_selectNum: "<<vec_selectHitNum_signal[n]<<"    "<<100*(double)vec_selectHitNum_signal[n]/(double)m_nHitSignal<<"%    "<<"noise:  "<<(double)(vec_selectHitNum[n]-vec_selectHitNum_signal[n])<<endl;
01201 //      }
01202 //  }
01203 
01204   /*
01205   //combine the neighbor cell
01206   //vector<int> peak_hitList[npeak1];
01207   vector< vector <int> > peak_hitList(npeak1,vector <int> () ); 
01208   for(int i=0;i<npeak1;i++){
01209   int bina=peakList1[0][i];
01210   int binr=peakList1[1][i];
01211   int hit_peak=peakList1[2][i];
01212   for(int j=0;j<hit_peak;j++){
01213   peak_hitList[i].push_back(vec_selectHit[bina][binr][j]);
01214 //  cout<<"j:  "<<j<<"  "<<vec_selectHit[bina][binr][j]<<endl;
01215 }
01216 
01217 for(int a=bina-1;a<=bina+1;a++){
01218 for(int rx=binr-1;rx<=binr+1;rx++){
01219 int ax;
01220 if (a<0) { ax=a+binX; }
01221 else if (a>=binX) { ax=a-binX; }
01222 else { ax=a; }
01223 if ( (ax!=bina || rx!=binr) && rx>=0 && rx<binY) {
01224 int hit_around=count[ax][rx]; 
01225 for(int k=0;k<hit_around;k++){
01226 bool combine=1;
01227 for(int m=0;m<hit_peak;m++){
01228 if(vec_selectHit[ax][rx][k]==peak_hitList[i][m]){
01229 //                              cout<<i<<"  "<<"  "<<ax<<"  "<<rx<<"  "<<k<<"  "<<vec_selectHit[ax][rx][k]<<endl;
01230 combine=0;
01231 continue;
01232 }
01233 }
01234 if(combine==1)  {
01235 peak_hitList[i].push_back(vec_selectHit[ax][rx][k]);
01236 peakList1[2][i]++;
01237 hit_peak++;
01238 //                cout<<" there is a combine"<<"  "<<vec_selectHit[ax][rx][k]<<endl;
01239 }
01240 }
01241 }
01242 }
01243 }
01244 }
01245 
01246 
01247 //for(int n=0;n<npeak1;n++){
01248 //  cout<<"peakList2:  "<<n<<"  "<<"a:  "<<peakList1[0][n]<<"  "<<"r:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl;
01249 //}
01250 
01251 
01252 for(int i=0;i<1;i++){
01253 int bina=peakList1[0][i];
01254 int binr=peakList1[1][i];
01255 int hit_peak=peakList1[2][i];
01256 for(int j=0;j<hit_peak;j++){
01257 int hitList=peak_hitList[i][j];
01258 //        cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
01259 //        cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
01260 }
01261 }
01262 
01263 //sort another time 
01264 int n_max_2;
01265 for (int na=0; na<npeak1-1; na++) {
01266 n_max_2=na;
01267 for (int nb=na+1; nb<npeak1; nb++) {
01268 if (peakList1[2][n_max_2]<peakList1[2][nb]) { n_max_2=nb; }
01269 }   
01270 if (n_max_2!=na) {        // swap
01271 float swap[3];
01272 for (int i=0; i<3; i++) {
01273 swap[i]=peakList1[i][na];
01274 peakList1[i][na]=peakList1[i][n_max_2];
01275 peakList1[i][n_max_2]=swap[i];
01276 }
01277 }
01278 }
01279 */
01280 
01281 //  for(int n=0;n<npeak1;n++){
01282 //      cout<<"peakList3:  "<<n<<"  "<<"a:  "<<peakList1[0][n]<<"  "<<"r:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl;
01283 //  }
01284 
01285 //  int max_count=-999;
01286 //  for(int i=0;i<binX;i++){
01287 //      for(int j=0;j<binY;j++){
01288 //        if(max_count<count[i][j]) {max_count=count[i][j];} 
01289 //      }
01290 //  }
01291 
01292 //  m_maxCount=max_count;
01293 //if(m_debug>0) cout<<"maxcount:   "<<max_count<<endl;
01294 
01295 //  for(int i=0;i<binX;i++){
01296 //      for(int j=0;j<binY;j++){
01297 //        //            cout<<"("<<i<<","<<j<<")"<<"  "<<"count:  "<<count[i][j]<<"  "<<endl;
01298 //        for(int l=0,m=vec_selectNum[i][j].size();l<m;++l){
01299 //              if(count[i][j]==max_count){
01300 //                //int layerId=vec_layer.at(vec_selectNum[i][j].at(l));
01301 //                //if ((layerId>=8&&layerId<=19)||(layerId>=36)) {     
01302 //                //    nHitAxialSelect++;
01303 //                //}
01304 //                vec_truthHit[vec_selectNum[i][j][l]]=true;  
01305 //
01306 //              }
01307 //        }
01308 //      }
01309 //}
01310 //  for(int i=0;i<n;i++){
01311 //      cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<endl;
01312 //  }
01313 //  cout<<"selected true hits is:    "<<endl;
01314 
01315 //for(int i=0;i<n;i++){
01316 //  if(vec_truthHit[i]==true){
01317 //    cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<"  "<<"hitSignal:"<<vec_hitSignal[i]<<endl;
01318 //  }
01319 //}
01320 
01321 
01322 //int nHitAxialSelect=0;
01323 //int nHitSelect=0;
01324 //int nHitSignal_select=0;
01325 //int nHitAxialSignal_select=0;
01326 //for(int i=0;i<n;i++){
01327 //  if(vec_truthHit[i]==false) continue;
01328 //  nHitSelect++;
01329 //  if(vec_hitSignal[i]==1) {nHitSignal_select++;}
01330 //  int layerId=vec_layer[i];
01331 //  if ((layerId>=8&&layerId<=19)||(layerId>=36)) {     
01332 //      nHitAxialSelect++;
01333 //      if(vec_hitSignal[i]==1) {nHitAxialSignal_select++;}
01334 //  }   
01335 //}
01336 //
01337 //m_nHitAxialSelect=nHitAxialSelect;
01338 //m_nHitSelect=nHitSelect;
01339 //
01340 //m_nHitSignal_select=nHitSignal_select;
01341 //m_nHitAxialSignal_select=nHitAxialSignal_select;
01342 
01343 //  backtransform all track-
01344 
01345 double d0=-999.;
01346 double phi0=-999.;
01347 double omega=-999.;
01348 double z0=0;
01349 double tanl=0;
01350 
01351 for(track_fit=0;track_fit<peak_track;track_fit++){
01352   for(int i=0;i<nhit;i++){
01353         vec_truthHit[i]=false;
01354   }
01355   cout<<"track: "<<track_fit<<"  has select: "<<peak_combine_num[track_fit]<<"  hit ."<<endl;
01356   for(int i=0;i<peak_combine_num[track_fit];i++){
01357         int hitNum=hit_combine[track_fit][i];
01358         vec_truthHit[hitNum]=true; 
01359   }
01360 
01361   int nHitAxialSelect_temp=10;
01362   int leastSquare=LeastSquare(nhit,nHitAxialSelect_temp,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0,phi0,omega);
01363 
01364   if(leastSquare==0){
01365 
01366         TrkExchangePar tt(d0,phi0,omega,z0,tanl);
01367         TrkCircleMaker circlefactory;
01368         float chisum =0.;
01369         TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
01370         bool permitFlips = true;
01371         bool lPickHits = m_pickHits;
01372         circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
01373         // combine hit list
01374         int nDigi = digiToHots(newTrack,vec_truthHit);
01375         int nRecTk = 0;
01376         //if(m_debug>0) newTrack->printAll(std::cout);
01377         //fit 
01378         TrkHitList* qhits = newTrack->hits();
01379         TrkErrCode err=qhits->fit();
01380         //cout<<"++++++++++++++++++"<<err.failure()<<endl;
01381         const TrkFit* newFit = newTrack->fitResult();
01382         bool fitSuccess = false;
01383         //test fit result
01384         if (!newFit || (newFit->nActive()<3)) {
01385           if(m_debug>0){
01386                 cout << "evt "<<t_eventNum<<" global fit failed ";
01387                 if(newFit) cout <<" nAct "<<newFit->nActive();
01388                 cout<<endl;
01389           }
01390           //return StatusCode::SUCCESS;
01391         }else{
01392           nRecTk = 1;
01393           fitSuccess = true;
01394           m_nEvtSuccess++;
01395           //    if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
01396           if(m_debug>0) newTrack->print(std::cout);
01397 
01398         }
01399 
01400         vector<int>  nfit2d(peak_track,0);
01401         if(m_debug>0) cout<<" hot list:"<<endl;
01402         TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
01403         while (hotIter!=qhits->hotList().end()) {
01404           nfit2d[track_fit]++;
01405           cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
01406                 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 
01407                 <<":"<<hotIter->isActive()<<")  ";
01408           hotIter++;
01409         }
01410 
01411         m_2d_nFit[track_fit]=nfit2d[track_fit];
01412 
01413         //-------------try to out put parameters---
01414         double x_cirtemp;
01415         double y_cirtemp;
01416         double r_temp;
01417         if(err.failure()==0){
01418           //cout<<"++++++++++++++++++"<<err.failure()<<endl;
01419           TrkExchangePar par = newFit->helix(0.);
01420           d0=par.d0();
01421           phi0=par.phi0();
01422           omega=par.omega();
01423           //    cout<<"d0:                  "<<par.d0()<<"           "<<"d0temp:                     "<<-d0_temp<<endl;
01424           //    cout<<"phi0:                 "<<par.phi0()<<"         "<<"phi0temp:                  "<<phi0_temp<<endl;
01425           //    cout<<"w:               "<<par.omega()<<"         "<<"omegatemp:                  "<<omega_temp<<endl;
01426           // vector<double> b; // vector<double> x_stereo;
01427           r_temp=1./-par.omega();
01428           x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
01429           y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
01430 
01431           m_pt[track_fit]=r_temp/333.567;
01432           cout<<"pt_rec:   "<<m_pt[track_fit]<<endl;
01433         }
01434         else{
01435           m_nFitFailure[track_fit]=2;
01436         }
01437         //caculate z position
01438 
01439         int z_stereoNum= Zposition(nhit,vec_slant,x_cirtemp,y_cirtemp,r_temp,vec_x_east,vec_x_west,vec_y_east, vec_y_west,vec_z_east,vec_z_west, vec_layer, vec_wire,vec_z,vec_zStereo,l,vec_truthHit);
01440         //cout<<"z_stereoNum: "<<z_stereoNum<<endl;
01441         //cout<<"if zposition is finish: "<<endl;
01442 
01443         m_zStereoNum=z_stereoNum;
01444         for(int i=0;i<z_stereoNum;i++){
01445 
01446           m_zStereoNum=z_stereoNum;
01447           m_zStereo[i]=vec_zStereo[i];
01448           //      cout<<"eventNum:   "<<t_eventNum<<"   "<<"1111111111111111"<<endl;
01449           m_l[i]=l[i];
01450           if(m_debug>0) cout<<" l: "<<m_l[i]<<"  "<<"z:  "<<vec_zStereo[i]<<endl;
01451         }
01452         //      cout<<"333333333333333333"<<endl;
01453 
01454         Linefit(z_stereoNum,vec_zStereo,l,tanl,z0);
01455         
01456         cout<<"tanl:   "<<tanl<<endl;
01457         cout<<"z0:   "<<z0<<endl;
01458         z0=dz_mc;
01459         tanl=tanl_mc;
01460         cout<<"tanltruth:   "<<tanl<<endl;
01461         cout<<"z0truth:   "<<z0<<endl;
01462         //-------------------------------------------5 parameter fit-----------------------
01463 
01464         TrkExchangePar tt2(d0,phi0,omega,z0,tanl);
01465         TrkHelixMaker helixfactory;
01466         chisum =0.;
01467         TrkRecoTrk* newTrack2 = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0);
01468         permitFlips = true;
01469         lPickHits = m_pickHits;
01470         helixfactory.setFlipAndDrop(*newTrack2, permitFlips, lPickHits);
01471         // combine hit list
01472         int nDigi2 = digiToHots2(newTrack2,vec_truthHit);
01473         //  int nRecTk = 0;
01474         //if(m_debug>0) newTrack->printAll(std::cout);
01475         //fit 
01476         TrkHitList* qhits2 = newTrack2->hits();
01477         TrkErrCode err2=qhits2->fit();
01478         //cout<<"++++++++++++++++++"<<err.failure()<<endl;
01479         const TrkFit* newFit2 = newTrack2->fitResult();
01480         //  bool fitSuccess = false;
01481 
01482         //test fit result
01483         if (!newFit2 || (newFit2->nActive()<5)) {
01484           //fit failed
01485           if(m_debug>0){
01486                 cout << "evt "<<t_eventNum<<" global fit failed ";
01487                 if(newFit2) cout <<" nAct "<<newFit2->nActive();
01488                 cout<<endl;
01489           }
01490           //return StatusCode::SUCCESS;
01491         }else{
01492           //fit success
01493           nFitSucess++;
01494           MdcTrack* mdcTrack = new MdcTrack(newTrack2);
01495           int tkStat = 4;//track find by Hough set stat=4
01496           int tkId = nTdsTk + track_fit;
01497           mdcTrack->storeTrack(tkId, trackList, hitList, tkStat);
01498           nRecTk = 1;
01499           fitSuccess = true;
01500           m_nEvtSuccess++;
01501           //    if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
01502           if(m_debug>0) newTrack2->print(std::cout);
01503 
01504         }
01505 
01506         vector<int>  nfit3d(peak_track,0);
01507         if(m_debug>0) cout<<" hot list:"<<endl;
01508         TrkHotList::hot_iterator hotIter2= qhits2->hotList().begin();
01509         while (hotIter2!=qhits2->hotList().end()) {
01510           nfit3d[track_fit]++;
01511           cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber()
01512                 <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber() 
01513                 <<":"<<hotIter2->isActive()<<")  ";
01514           hotIter2++;
01515         }
01516 
01517         m_3d_nFit[track_fit]=nfit3d[track_fit];
01518 
01519         // -------------try to out put parameters---
01520         if(err2.failure()==0){
01521           //cout<<"++++++++++++++++++"<<err2.failure()<<endl;
01522           TrkExchangePar par2 = newFit2->helix(0.);
01523 
01524           //cout<<"d0:                  "<<par2.d0()<<endl;
01525           //cout<<"phi0:                 "<<par2.phi0()<<endl;
01526           //cout<<"w:               "<<par2.omega()<<endl;
01527           //cout<<"z:               "<<par2.z0()<<endl;
01528           //cout<<"tanl:               "<<par2.tanDip()<<endl;
01529 
01530           r_temp=1./-par2.omega();
01531           x_cirtemp=(r_temp-par2.d0())*cos(par2.phi0()-M_PI/2.);
01532           y_cirtemp=(r_temp-par2.d0())*sin(par2.phi0()-M_PI/2.);
01533 
01534           m_d0[track_fit]=par2.d0();
01535           m_phi0[track_fit]=par2.phi0();
01536           m_omega[track_fit]=par2.omega();
01537           m_z0[track_fit]=par2.z0();
01538           m_tanl[track_fit]=par2.tanDip();
01539 
01540           double pxy=r_temp/333.567;
01541           double pz=pxy*par2.tanDip();
01542           double pxyz=sqrt(pxy*pxy+pz*pz);
01543           m_pt2[track_fit]=pxy;
01544           m_pz[track_fit]=pxy*par2.tanDip();
01545           m_pxyz[track_fit]=pxyz;
01546           cout<<"track"<<track_fit<<":  "<<"pt_rec2:   "<<m_pt2[track_fit]<<endl;
01547           cout<<"eventNum:  "<<t_eventNum<<"   "<<"track"<<track_fit<<":  "<<"pz_rec:   "<<m_pz[track_fit]<<endl;
01548         }
01549         else{
01550           m_nFitFailure[track_fit]=3;
01551         }
01552         // vector<double> b; // vector<double> x_stereo;
01553         //      r_temp=1./-par.omega();
01554         //      x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
01555         //      y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
01556         //      if(m_data==0) {
01557         //        m_mc_pt=r_temp/333.567;
01558         //        cout<<"pt_rec:   "<<m_mc_pt<<endl;
01559         //      }
01560         //      if(m_data==1) {m_pt=r_temp/333.567;}
01561         //  
01562         //  else{
01563         //      if(m_data==0) {m_mc_nfitfailure=2;}
01564         //      if(m_data==1){m_nFitFailure=2;}
01565         //  }
01566 
01567 
01568         //------------------------------------clean-------------------------
01569   }
01570   //vector<double>().swap(vec_u);
01571   //vector<double>().swap(vec_v);
01572   //vector<double>().swap(vec_p);
01573   //vector<double>().swap(vec_x);
01574   //vector<double>().swap(vec_y);
01575   //vector<double>().swap(vec_z);
01576   //vector<double>().swap(vec_x_east);
01577   //vector<double>().swap(vec_x_west);
01578   //vector<double>().swap(vec_y_east);
01579   //vector<double>().swap(vec_y_west);
01580   //vector<double>().swap(vec_z_east);
01581   //vector<double>().swap(vec_z_west);
01582   //vector<double>().swap(vec_zStereo);
01583   //vector<double>().swap(l);
01584   //vector<int>().swap(vec_layer);
01585   //vector<int>().swap(vec_wire);
01586 }
01587 m_nFitSucess=nFitSucess;
01588 cout<<"  Hough finder find "<<nFitSucess<<" tot "<<nTdsTk<< endl;
01589 
01590 delete h1;
01591 delete h2;
01592 cout<<"break:  ????"<<endl;
01593 ntuplehit->write();
01594 cout<<"finish event  "<<t_eventNum<<endl;
01595 cout<<endl;
01596 cout<<endl;
01597 return StatusCode::SUCCESS;
01598 }

void MdcHoughFinder::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]`

Definition at line 2085 of file MdcHoughFinder.cxx.

References cos(), genRecEmupikp::i, ganga-rec::j, m_debug, M_PI, sin(), and t_maxP.

Referenced by execute().

02085                                                                                                                                                                                                                                                                                    {
02086   double binWidth=M_PI/binX; 
02087   double binHigh=2*t_maxP/binY;
02088 
02089   //  vector<int>** selectNum=new vector<int>*[binX];
02090   //  for(int i=0;i<binX;i++){
02091   //    selectNum[i]=new vector<int>[binY];
02092   //    for(int j=0;j<binY;j++){
02093   //      selectNum[i][j]=new vector<int>;
02094   //    }
02095   //  }
02096   //  TH2D *h = new TH2D("hough","hough",binX,0,3.14,binY,-max_p,max_p);
02097 
02098   for(int i=0;i<binX;i++){
02099         for(int j=0;j<binY;j++){
02100           int count=0;
02101           double binLeft=i*binWidth;
02102           double binRight=(i+1)*binWidth;
02103           double binDown=-t_maxP+j*binHigh;
02104           double binUp=-t_maxP+(j+1)*binHigh;
02105           double binMid=(binLeft+binRight)/2;
02106           for(int k=0;k<n;k++){
02107                 double lineLeft=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
02108                 double lineRight=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
02109                 double lineMid=vec_u[k]*cos(binMid)+vec_v[k]*sin(binMid);
02110                 if(((lineLeft<binUp)&&(lineLeft>binDown)) || ((lineRight<binUp)&&(lineRight>binDown)) || (((lineLeft-binUp)*(binDown-lineRight)>0))){ 
02111                   vec_selectNum[i][j].push_back(k);
02112                   int layerId=vec_layer[k];
02113                   int cellId=vec_wire[k];
02114                   int hit_temp=(int)layerId*1000+cellId;
02115                   vec_selectHit[i][j].push_back(hit_temp);
02116                   count++;
02117                 }
02118           }
02119           countij[i][j]=count;
02120 
02121           //for(int l=0,m=vec_selectNum[i][j].size();l<m;++l){
02122           //  {//         cout << "section " << "(" << i <<  ", " << j << ") statistics: " << count[i][j] << endl;
02123           //       cout << "("<<i<<","<<j<<")"<<": "<<vec_selectNum[i][j][l] << ", "<<endl;
02124           //    // trueTrk[vec_selectNum[i][j].at(l)]=true;  
02125           //  }
02126           //}
02127           if(m_debug>0)  cout<<"countij:  "<<countij[i][j]<<endl;
02128         }
02129   }
02130   // for(int i=0;i<binX;i++){
02131   //   for(int j=0;j<binY;j++){
02132   //     for(int k=0;k<countij[i][j];k++){
02133   //    cout<<"("<<i<<","<<j<<")"<<vec_selectNum[i][j][k]<<endl;
02134   //     }
02135   //   }
02136   // }
02137 
02138   for(int i=0;i<binX;i++){
02139         for(int j=0;j<binY;j++){
02140           h1->SetBinContent(i+1,j+1,countij[i][j]);
02141         }
02142   }
02143 }

void MdcHoughFinder::FillHist	(	TH2D *	h1,
		vector< double >	vec_u,
		vector< double >	vec_v,
		int	m_nHit,
		vector< vector< vector< int > > >	vec_selectNum
	)			`[private]`

Definition at line 2175 of file MdcHoughFinder.cxx.

References cos(), genRecEmupikp::i, ganga-rec::j, M_PI, sin(), and t_maxP.

Referenced by execute().

02175                                                                                                                                          {
02176   for(int n=0;n<nhit;n++){
02177         for(int i=0;i<binX;i++){
02178           double binwidth=M_PI/binX;
02179           double binhigh=2*t_maxP/binY;
02180           double binMid=(i+0.5)*binwidth;
02181           double rho=vec_u[n]*cos(binMid)+vec_v[n]*sin(binMid);
02182           int j=(int)((rho+t_maxP)/binhigh);
02183           int count=0;
02184           count=h1->GetBinContent(i+1,j+1);
02185           count++;
02186           h1->SetBinContent(i+1,j+1,count);
02187         }
02188   }
02189   //for(int i=0;i<binX;i++){
02190   //  double binwidth=M_PI/binX;
02191   //  double binhigh=2*t_maxP/binY;
02192   //  double binMid=(i+0.5)*binwidth;
02193   //  for(int n=0;n<nhit;n++){
02194   //    double rho=vec_u[n]*cos(binMid)+vec_v[n]*sin(binMid);
02195   //    int j=(int)((rho+t_maxP)/binhigh);
02196   //    int count=h1->GetBinContent(i+1,j+1);
02197   //    count++;
02198   //    h1->SetBinContent(i+1,j+1,count);
02199   //    //        cout<<"i: "<<i<<"  "<<"j: "<<j<<"  "<<"n: "<<n<<endl;
02200   //    vec_selectNum[i][j].push_back(n);
02201   //  }
02202   //}
02203 }

void MdcHoughFinder::FillRhoTheta	(	TH2D *	h1,
		vector< double >	vec_theta,
		vector< double >	vec_rho,
		int	numCross
	)			`[private]`

Definition at line 2204 of file MdcHoughFinder.cxx.

References genRecEmupikp::i, M_PI, and t_maxP.

Referenced by execute().

02204                                                                                                        {
02205   for(int i=0;i<numCross;i++){
02206         double binwidth=M_PI/binX;
02207         double binhigh=2*t_maxP/binY;
02208         int binxNum=vec_theta[i]/binwidth+1;
02209         int binyNum=(vec_rho[i]+t_maxP)/binhigh+1; 
02210         int count =h1->GetBinContent(binxNum,binyNum);
02211         count++;
02212         h1->SetBinContent(binxNum,binyNum,count);
02213   }
02214 }

StatusCode MdcHoughFinder::finalize ( )

Definition at line 1601 of file MdcHoughFinder.cxx.

References Bes_Common::INFO, and msgSvc().

01601                                     {
01602   MsgStream log(msgSvc(), name());
01603   log << MSG::INFO << "in finalize()" << endreq;
01604   return StatusCode::SUCCESS;
01605 }

int MdcHoughFinder::GetMcInfo ( ) [private]

Definition at line 1692 of file MdcHoughFinder.cxx.

References Helix::a(), Helix::cosTheta(), Helix::dr(), Helix::dz(), dz_mc, calibUtil::ERROR, Helix, m_costaTruth, m_debug, m_drTruth, m_dzTruth, m_phi0Truth, m_pid, m_pidTruth, m_pTruth, m_ptTruth, m_pzTruth, m_qTruth, Helix::momentum(), msgSvc(), Helix::phi0(), pid, Helix::pivot(), Helix::pt(), deljobs::string, t_eventNum, t_nTrkMC, t_runNum, Helix::tanl(), tanl_mc, and Bes_Common::WARNING.

Referenced by execute().

01692                              {
01693   //t_t0Truth=-1;
01694   MsgStream log(msgSvc(), name());
01695   StatusCode sc;
01696   SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01697   if (!mcParticleCol) {
01698         log << MSG::ERROR << "Could not find McParticle" << endreq;
01699         return -999;
01700   }
01701   int itk = 0;
01702   int t_qTruth = 0;
01703   int t_pidTruth = -999;
01704   int t_McTkId = -999;
01705   t_nTrkMC=0;
01706   Helix* mchelix;
01707   McParticleCol::iterator iter_mc = mcParticleCol->begin();
01708   for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
01709         if(!(*iter_mc)->primaryParticle()) continue;
01710         t_pidTruth = (*iter_mc)->particleProperty();
01711 
01712         if(m_debug>2) cout<< "tk "<<itk<< " pid="<< t_pidTruth<< endl;
01713         if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
01714 
01715         //t_t0Truth=(*iter_mc)->initialPosition().t();
01716 
01717         if(itk>=10) break;
01718 
01719         int pid = t_pidTruth; 
01720         if( pid == 0 ) { 
01721           log << MSG::WARNING << "Wrong particle id" <<endreq;
01722           continue;
01723         }else{
01724           IPartPropSvc* p_PartPropSvc;
01725           static const bool CREATEIFNOTTHERE(true);
01726           StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
01727           if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
01728                 std::cout<< " Could not initialize Particle Properties Service" << std::endl;
01729           }
01730           HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT(); 
01731           std::string name;
01732           if( p_particleTable->particle(pid) ){
01733                 name = p_particleTable->particle(pid)->name();
01734                 t_qTruth = p_particleTable->particle(pid)->charge();
01735           }else if( p_particleTable->particle(-pid) ){
01736                 name = "anti " + p_particleTable->particle(-pid)->name();
01737                 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
01738           }
01739           if(m_debug>2) std::cout << " particle "<<name<<" charge= "<<t_qTruth<<std::endl;
01740         }
01741 
01742         t_McTkId = (*iter_mc)->trackIndex();
01743 
01744         double px = (*iter_mc)->initialFourMomentum().px();//GeV
01745         double py = (*iter_mc)->initialFourMomentum().py();//GeV
01746         double pz = (*iter_mc)->initialFourMomentum().pz();//GeV
01747 
01748         m_pzTruth=pz;
01749 
01750         mchelix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);//cm
01751         mchelix->pivot( HepPoint3D(0.,0.,0.) );
01752 
01753         if(m_debug>0){
01754           std::cout<<"Truth tk "<<itk<<" pid "<<pid<<" charge "<<t_qTruth
01755                 << " helix "<< mchelix->a()<<" p "<<mchelix->momentum(0.)<<std::endl;
01756         }
01757         t_nTrkMC++;
01758         itk++;
01759   }
01760   if(t_nTrkMC!=1) {
01761         std::cout<<"WARNING run "<<t_runNum<<" evt "<<t_eventNum<<" not single event. nTrkMc="<<t_nTrkMC<<std::endl;
01762         return -999;
01763   }else{
01764         if(m_debug>2) std::cout<<"nTrkMc=1"<<std::endl;
01765   }
01766   //--fill truth in evt 
01767 
01768   m_pidTruth   = t_pidTruth;
01769   m_costaTruth = mchelix->cosTheta();
01770   m_phi0Truth  = mchelix->phi0();
01771   m_drTruth    = mchelix->dr();
01772   m_dzTruth    = mchelix->dz();
01773   m_ptTruth    = mchelix->pt();
01774   m_pTruth     = mchelix->momentum(0.).mag();
01775   m_qTruth     = t_qTruth;
01776   dz_mc        =mchelix->dz();
01777   tanl_mc      =mchelix->tanl();
01778   return 1;
01779 }

StatusCode MdcHoughFinder::initialize ( )

Definition at line 96 of file MdcHoughFinder.cxx.

References BField::bFieldNominal(), calibUtil::ERROR, Bes_Common::FATAL, Bes_Common::INFO, m_2d_nFit, m_3d_nFit, m_areaLeast, m_areaLeastNum, m_areaSelectHit, m_areaSelectHit_signal, m_bfield, m_cell, m_cell_Mc, m_context, m_cosCut, m_costaTruth, m_d0, TrkHelixFitter::m_debug, m_debug, m_drTruth, m_dzTruth, m_eventNum, m_hitSignal, m_l, m_layer, m_layer_Mc, m_layerNhit, m_maxCount, m_mdcCalibFunSvc, m_mdcGeomSvc, m_nCross, m_nFitFailure, m_nFitSucess, m_nHit, m_nHit_Mc, m_nHitAxial, m_nHitAxialSignal, m_nHitAxialSignal_select, m_nHitSignal, m_nHitSignal_select, m_npeak, m_omega, m_p, m_particleTable, m_pdtFile, m_peakArea, m_peakHigh, m_peakWidth, m_phi0, m_phi0Truth, m_pidTruth, m_pIMF, m_pt, m_pt2, m_pTruth, m_ptTruth, m_pxyz, m_pz, m_pzTruth, m_qTruth, m_r, m_rawDataProviderSvc, m_rho, m_runNum, m_slant, m_tanl, m_theta, m_trackNum, m_trackNum_Mc, m_u, m_v, m_x, m_x_circle, m_x_east, m_x_west, m_y, m_y_circle, m_y_east, m_y_west, m_z, m_z0, m_z_east, m_z_west, m_zStereo, m_zStereoNum, msgSvc(), ntuplehit, ntupleSvc(), and Pdt::readMCppTable().

00096                                      {
00097 
00098   MsgStream log(msgSvc(), name());
00099   log << MSG::INFO << "in initialize()" << endreq;
00100 
00101   StatusCode sc;
00102 
00103   //for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
00104 
00105   //nfailure=0;
00106   IPartPropSvc* p_PartPropSvc;
00107   static const bool CREATEIFNOTTHERE(true);
00108   sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00109   if (!sc.isSuccess() || 0 == p_PartPropSvc) {
00110         log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
00111         return sc;
00112   }
00113   m_particleTable = p_PartPropSvc->PDT();
00114 
00115   // RawData
00116   IRawDataProviderSvc* irawDataProviderSvc;
00117   sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00118   m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00119   if ( sc.isFailure() ){
00120         log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
00121         return StatusCode::FAILURE;
00122 
00123   }
00124 
00125   //  Geometry
00126   IMdcGeomSvc*   imdcGeomSvc;
00127   sc = service ("MdcGeomSvc", imdcGeomSvc);
00128   m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
00129   if ( sc.isFailure() ){
00130         log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
00131         return StatusCode::FAILURE;
00132   }
00133 
00134   sc = service ("MagneticFieldSvc",m_pIMF);
00135   if(sc!=StatusCode::SUCCESS) {
00136         log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00137   }
00138   m_bfield = new BField(m_pIMF);
00139   log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00140 
00141   m_context = new TrkContextEv(m_bfield);
00142 
00143   Pdt::readMCppTable(m_pdtFile);
00144 
00145   //Get MdcCalibFunSvc
00146   IMdcCalibFunSvc* imdcCalibSvc; 
00147   sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00148   m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00149   if ( sc.isFailure() ){
00150         log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00151         return StatusCode::FAILURE;
00152   }
00153 
00154   //initialize ntuplehit
00155   NTuplePtr nt(ntupleSvc(), "MdcHoughFinder/hit");
00156   if ( nt ){
00157         ntuplehit = nt;
00158   } else {
00159         ntuplehit = ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
00160         if(ntuplehit){
00161           ntuplehit->addItem("eventNum",           m_eventNum);
00162           ntuplehit->addItem("runNum",             m_runNum);
00163           ntuplehit->addItem("costaCut",        m_cosCut);
00164           //    test the first layer in hough space
00165           //      ntuplehit->addItem("Layer1Num",               m_nLayerNum,0, 43);
00166           //      ntuplehit->addItem("layer1X",    m_nLayerNum,   m_layer1X);
00167           //      ntuplehit->addItem("layer1Y",    m_nLayerNum,   m_layer1Y);
00168           //      ntuplehit->addItem("layer1U",    m_nLayerNum,   m_layer1U);
00169           //      ntuplehit->addItem("layer1V",    m_nLayerNum,   m_layer1V);
00170 
00171 
00172           ntuplehit->addItem("nHitMc",              m_nHit_Mc,0, 6796);
00173           ntuplehit->addItem("layerMc",    m_nHit_Mc,   m_layer_Mc);
00174           ntuplehit->addItem("cellMc",     m_nHit_Mc,   m_cell_Mc);
00175 
00176 
00177           ntuplehit->addItem("nHit",               m_nHit,0, 6796);
00178           ntuplehit->addItem("layerNhit",   43,    m_layerNhit);
00179           //      ntuplehit->addItem("hitCol",   m_nHit,   m_hitCol);
00180           ntuplehit->addItem("nCross",           m_nCross,0, 125000);
00181           ntuplehit->addItem("rho",    m_nCross,   m_rho);
00182           ntuplehit->addItem("theta",    m_nCross,   m_theta);
00183 
00184           ntuplehit->addItem("hitSignal",    m_nHit,   m_hitSignal);
00185           ntuplehit->addItem("layer",    m_nHit,   m_layer);
00186           ntuplehit->addItem("cell",     m_nHit,   m_cell);
00187           ntuplehit->addItem("x_east",   m_nHit,   m_x_east);
00188           ntuplehit->addItem("y_east",   m_nHit,   m_y_east);
00189           ntuplehit->addItem("z_east",   m_nHit,   m_z_east);
00190           ntuplehit->addItem("x_west",   m_nHit,   m_x_west);
00191           ntuplehit->addItem("y_west",   m_nHit,   m_y_west);
00192           ntuplehit->addItem("z_west",   m_nHit,   m_z_west);
00193           ntuplehit->addItem("x",        m_nHit,   m_x);
00194           ntuplehit->addItem("y",        m_nHit,   m_y);
00195           ntuplehit->addItem("z",        m_nHit,   m_z);
00196           ntuplehit->addItem("u",        m_nHit,   m_u);
00197           ntuplehit->addItem("v",        m_nHit,   m_v);
00198           ntuplehit->addItem("p",        m_nHit,   m_p);
00199           ntuplehit->addItem("slant",    m_nHit,   m_slant);
00200           //            ntuplehit->addItem("z_stereo",    m_stereohit,   m_zStereo);
00201           //            ntuplehit->addItem("z_num",    m_stereohit,   m_z_num);
00202           ntuplehit->addItem("maxCount",          m_maxCount);
00203 
00204           ntuplehit->addItem("peakColumn",           m_npeak,0,100);
00205           ntuplehit->addItem("peakWidth",  m_npeak,       m_peakWidth);
00206           ntuplehit->addItem("peakHigh",   m_npeak,       m_peakHigh);
00207           ntuplehit->addItem("peakArea",   m_npeak,       m_peakArea);
00208           ntuplehit->addItem("peakAreaLeast",       m_areaLeast);
00209           ntuplehit->addItem("peakAreaLeastNum",       m_areaLeastNum);
00210           ntuplehit->addItem("peakHitSelect",       m_areaSelectHit);
00211           ntuplehit->addItem("peakHitSelectSignal",       m_areaSelectHit_signal);
00212           //            ntuplehit->addItem("mc_x",    m_Mcnhit,  mc_x);
00213           //            ntuplehit->addItem("mc_y",    m_Mcnhit,  mc_y);
00214           //            ntuplehit->addItem("mc_z",    m_Mcnhit,  mc_z);
00215           ntuplehit->addItem("circleCenterX",      m_x_circle);
00216           ntuplehit->addItem("circleCenterY",      m_y_circle);
00217           ntuplehit->addItem("circleR",            m_r);
00218 
00219 
00220           ntuplehit->addItem("zStereoNum",           m_zStereoNum,0,1000);
00221           ntuplehit->addItem("zStereo",     m_zStereoNum,      m_zStereo );
00222           ntuplehit->addItem("l",            m_zStereoNum,      m_l);
00223 
00224           ntuplehit->addItem("trackNum_Mc",     m_trackNum_Mc,  0,100);
00225           ntuplehit->addItem("trackNum",     m_trackNum,  0,100);
00226           ntuplehit->addItem("d0",           m_trackNum,      m_d0);
00227           ntuplehit->addItem("phi0",         m_trackNum,      m_phi0);
00228           ntuplehit->addItem("omega",        m_trackNum,      m_omega);
00229           ntuplehit->addItem("z0",           m_trackNum,      m_z0);
00230           ntuplehit->addItem("tanl",         m_trackNum,      m_tanl);
00231 
00232           ntuplehit->addItem("pt_rec",     m_trackNum,     m_pt);
00233           ntuplehit->addItem("pt2_rec",    m_trackNum,     m_pt2);
00234           ntuplehit->addItem("pz_rec",     m_trackNum,     m_pz);
00235           ntuplehit->addItem("p_rec",      m_trackNum,     m_pxyz);
00236 
00237           ntuplehit->addItem("nHitSignal",                  m_nHitSignal);
00238           ntuplehit->addItem("nHitAxialSignal",             m_nHitAxialSignal);
00239           ntuplehit->addItem("nHitSignal_select",      m_trackNum,     m_nHitSignal_select);
00240           ntuplehit->addItem("nHitAxialSignal_selcet", m_trackNum,     m_nHitAxialSignal_select);
00241           ntuplehit->addItem("fitFailure",    m_trackNum,      m_nFitFailure);
00242           ntuplehit->addItem("nFit",             m_trackNum,   m_2d_nFit);
00243           ntuplehit->addItem("3dnFit",           m_trackNum,   m_3d_nFit);
00244           ntuplehit->addItem("nHitAxial",          m_nHitAxial);
00245           ntuplehit->addItem("nFitSucess",          m_nFitSucess);
00246           //      ntuplehit->addItem("nHitAxialSelect",   m_nHitAxialSelect);
00247           //      ntuplehit->addItem("nHitSelect",        m_trackNum,  m_nHitSelect);
00248 
00249           ntuplehit->addItem("pidTruth",   m_pidTruth);
00250           ntuplehit->addItem("costaTruth", m_costaTruth);
00251           ntuplehit->addItem("phiTruth",   m_phi0Truth);
00252           ntuplehit->addItem("drTruth",    m_drTruth);
00253           ntuplehit->addItem("dzTruth",    m_dzTruth);
00254           ntuplehit->addItem("ptTruth",    m_ptTruth);
00255           ntuplehit->addItem("pzTruth",    m_pzTruth);
00256           ntuplehit->addItem("pTruth",     m_pTruth);
00257           ntuplehit->addItem("qTruth",     m_qTruth);
00258 
00259 
00260         } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hit" << endmsg;
00261           return StatusCode::FAILURE;
00262         }
00263   }
00264 
00265 
00266   if(m_debug>2)TrkHelixFitter::m_debug = true;
00267 
00268   return StatusCode::SUCCESS;
00269 
00270 
00271 }

int MdcHoughFinder::LeastSquare	(	int	n,
		int	nselecthit_axial,
		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
	)			`[private]`

Definition at line 1782 of file MdcHoughFinder.cxx.

References genRecEmupikp::i, m_debug, M_PI, m_r, m_x_circle, and m_y_circle.

Referenced by execute().

01782                                                                                                                                                                                                                 {
01783   //if(nHitAxialSelect<3) {
01784   //m_nFitFailure=1;
01785   //    return 1;
01786   // }
01787   //else{
01788   double x_sum=0;
01789   double y_sum=0;
01790   double x2_sum=0;
01791   double y2_sum=0;
01792   double x3_sum=0;
01793   double y3_sum=0;
01794   double x2y_sum=0;
01795   double xy2_sum=0;
01796   double xy_sum=0;
01797   double a1=0;
01798   double a2=0;
01799   double b1=0;
01800   double b2=0;
01801   double c1=0;
01802   double c2=0;
01803   double x_aver,y_aver,r2;
01804   for(int i=0;i<n;i++){
01805         if(vec_truthHit[i]==false) continue;
01806         if(vec_slant[i]!=0) continue;
01807         //cout<<"hitNum:  "<<i<<"  "<<"layer:"<<"  "<<vec_layer[i]<<"  "<<endl;
01808         x_sum=x_sum+vec_x[i];
01809         y_sum=y_sum+vec_y[i];
01810         x2_sum=x2_sum+vec_x[i]*vec_x[i];
01811         y2_sum=y2_sum+vec_y[i]*vec_y[i];
01812         x3_sum=x3_sum+vec_x[i]*vec_x[i]*vec_x[i];
01813         y3_sum=y3_sum+vec_y[i]*vec_y[i]*vec_y[i];
01814         x2y_sum=x2y_sum+vec_x[i]*vec_x[i]*vec_y[i];
01815         xy2_sum=xy2_sum+vec_x[i]*vec_y[i]*vec_y[i];
01816         xy_sum=xy_sum+vec_x[i]*vec_y[i];
01817   }
01818   a1=x2_sum-x_sum*x_sum/n;
01819   a2=xy_sum-x_sum*y_sum/n;
01820   b1=a2;
01821   b2=y2_sum-y_sum*y_sum/n;
01822   c1=(x3_sum+xy2_sum-x_sum*(x2_sum+y2_sum)/n)/2.;
01823   c2=(y3_sum+x2y_sum-y_sum*(x2_sum+y2_sum)/n)/2.;
01824 
01825   x_aver=x_sum/n;
01826   y_aver=y_sum/n;
01827 
01828   m_x_circle =(b1*c2-b2*c1)/(b1*b1-a1*b2);
01829   m_y_circle =(b1*c1-a1*c2)/(b1*b1-a1*b2);
01830 
01831   double  x_cirtemp =(b1*c2-b2*c1)/(b1*b1-a1*b2);
01832   double  y_cirtemp =(b1*c1-a1*c2)/(b1*b1-a1*b2);
01833 
01834   r2=(x2_sum+y2_sum-2*x_cirtemp*x_sum-2*y_cirtemp*y_sum)/n+x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp; 
01835 
01836   m_r=sqrt(r2);
01837 
01838   double r_temp=sqrt(r2);
01839   double d0_temp=sqrt(x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp)-r_temp;
01840   //cout<<"rtemp:   "<<r_temp<<"  "<<endl;
01841   //cout<<"xtemp:   "<<x_cirtemp<<"  "<<endl;
01842   //cout<<"ytemp:   "<<y_cirtemp<<"  "<<endl;
01843   //cout<<"d0temp:  "<<-d0_temp<<"    "<<endl;
01844   double  pt_temp=r_temp/333.567;
01845   cout<<"pt_temp: "<<pt_temp<<endl;
01846   //--------------------------------------------------add drift and fit---------------------------
01847 
01848   d0=-d0_temp;
01849   // double phi0=atan2(y_cirtemp,x_cirtemp);
01850   // double omega=-1./r_temp;
01851   phi0=atan2(y_cirtemp,x_cirtemp)+M_PI/2.;
01852   omega=-1./r_temp;
01853   double z0=0.;
01854   double tanl=0.;
01855 
01856   //    double phi0_temp=phi0;
01857   //    double  omega_temp=omega;
01858   //d0=0;
01859   //phi0=-3.08301;
01860   //omega=-0.0117116;
01861   //double z0=0.0;
01862   //double tanl=0.0;
01863   if(m_debug<0)std::cout<<" before global fit  Helix PatPar :"<<d0<<","<<phi0<<","<<omega<<","<<z0<<","<<tanl<<  std::endl;
01864   return 0;
01865   //}
01866 }

void MdcHoughFinder::Linefit	(	int	z_stereoNum,
		vector< double >	z_stereo_aver,
		vector< double >	l,
		double &	tanl,
		double &	z0
	)			`[private]`

Definition at line 1985 of file MdcHoughFinder.cxx.

References genRecEmupikp::i.

Referenced by execute().

01985                                                                                                                {
01986 
01987   double l_sum=0;
01988   double z_sum=0;
01989   double l2_sum=0;
01990   double z2_sum=0;
01991   double lz_sum=0;
01992   for(int i=0;i<z_stereoNum;i++){
01993         //      cout<<"event:  "<<t_eventNum<<"  "<<"z:  "<<vec_zStereo[i]<<"  "<<"l:  "<<l[i]<<endl;
01994         l_sum=l_sum+l[i];
01995         z_sum=z_sum+vec_zStereo[i];
01996         l2_sum=l2_sum+l[i]*l[i];
01997         //      cout<<l2_sum<<endl;
01998         z2_sum=z2_sum+vec_zStereo[i]*vec_zStereo[i];
01999         lz_sum=lz_sum+l[i]*vec_zStereo[i];
02000   }
02001   double b=(l2_sum*z_sum-lz_sum*l_sum)/(z_stereoNum*l2_sum-l_sum);
02002   double k=(z_stereoNum*lz_sum-l_sum*z_sum)/(z_stereoNum*l2_sum-l_sum*l_sum);
02003   z0=b;
02004   tanl=k;
02005   cout<<k<<"  "<<b<<endl;
02006 }

void MdcHoughFinder::Multi_array	(	int	binX,
		int	binY
	)			`[private]`

Definition at line 2070 of file MdcHoughFinder.cxx.

References genRecEmupikp::i, and ganga-rec::j.

02070                                                   {
02071   int **count=new int*[binX];
02072   for(int i=0;i<binX;i++){
02073         count[i]=new int [binY];
02074   }
02075   //vector<int> vec_selectNum[binX][binY];
02076   int ***selectNum=new int**[binX];
02077   for(int i=0;i<binX;i++){
02078         selectNum[i]=new int* [binY];
02079         for(int j=0;j<binY;j++){
02080           //selectNum[i][j]=new int[];
02081           selectNum[i][j]=new int [count[i][j]];
02082         }
02083   }
02084 }

void MdcHoughFinder::RhoTheta	(	int	numCross,
		int	m_nHit,
		vector< double >	vec_u,
		vector< double >	vec_v,
		vector< double > &	vec_rho,
		vector< double > &	vec_theta,
		vector< vector< int > > &	vec_hitNum
	)			`[private]`

Definition at line 2144 of file MdcHoughFinder.cxx.

References genRecEmupikp::i, ganga-rec::j, m_nCross, M_PI, m_rho, and m_theta.

Referenced by execute().

02144                                                                                                                                                                                 {
02145   for(int i=0;i<nhit;i++){
02146         for(int j=i+1;j<nhit;j++){
02147           double k,b,x_cross,y_cross;
02148           double rho_temp,theta_temp; 
02149           k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]);
02150           b=vec_v[i]-k*vec_u[i];
02151           x_cross=-b/(k+1/k);
02152           y_cross=b/(k*k+1);
02153           rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross);
02154           theta_temp=atan2(y_cross,x_cross);
02155           if(theta_temp<0)  {
02156                 theta_temp=theta_temp+M_PI;
02157                 rho_temp=-rho_temp;
02158           }
02159           vec_rho.push_back(rho_temp);
02160           vec_theta.push_back(theta_temp);
02161           vec_hitNum[0].push_back(i);
02162           vec_hitNum[1].push_back(j);
02163           //cout<<"u:"<<vec_u[i]<<" "<<"v:"<<vec_v[i]<<"  "<<"("<<vec_theta[i]<<","<<vec_rho[i]<<")"<<endl;
02164         }
02165   }
02166   //cout<<"numCross: "<<numCross<<endl;
02167   m_nCross=numCross;
02168 
02169   for(int i=0;i<numCross;i++){
02170         m_rho[i]=vec_rho[i];
02171         m_theta[i]=vec_theta[i];
02172         //cout<<" "<<vec_theta[i]<<" "<<vec_rho[i]<<endl;
02173   }     
02174 }

int MdcHoughFinder::SlantId ( int layer ) [private]

Definition at line 1678 of file MdcHoughFinder.cxx.

Referenced by execute().

01678                                     {
01679   if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) { 
01680         if ((layer>=0&&layer<=3)||(layer>=20&&layer<=23)||(layer>=28&&layer<=31)){
01681           //      m_slant[digiId]=-1;
01682           return -1;
01683         }else{ 
01684           return 1;
01685         }
01686   } else{
01687         return 0;
01688   }
01689 }

int MdcHoughFinder::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 1868 of file MdcHoughFinder.cxx.

References abs, genRecEmupikp::i, and M_PI.

Referenced by execute().

01868                                                                                                                                                                                                                                                                                                                                                                                                             {
01869   double x_stereo;
01870   double y_stereo;
01871   //  double theta;
01872   int stereoId=0;
01873   //vector<double> z_stereo;
01874   // vector<double> vec_zStereo;
01875   int nDropDelta=0;
01876   for(int i=0;i<n;i++){
01877         if(vec_truthHit[i]==false) continue;
01878         double k;
01879         double b;
01880         if(vec_slant[i]!=0){
01881           if(vec_x_east[i]!=vec_x_west[i]){
01882                 k=(vec_y_east[i]-vec_y_west[i])/(vec_x_east[i]-vec_x_west[i]);
01883                 b=-k*vec_x_east[i]+vec_y_east[i];
01884                 //        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;
01885                 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);
01886                 //cout<<"layer:  "<<vec_layer[i]<<"wire:  "<<vec_wire[i]<<endl;
01887                 //cout<<"x_east:  "<<vec_x_east[i]<<"  "<<"x_west:   "<<vec_x_west[i]<<endl;
01888                 //cout<<"k:   "<<k<<"  "<<"b:  "<<b<<endl;
01889                 //cout<<"x_cirtemp:  "<<x_cirtemp<<"  "<<"y_cirtemp:  "<<y_cirtemp<<endl;
01890                 if(delta<0) {
01891                   nDropDelta++;
01892                   continue;
01893                 }
01894                 //cout<<"delta:     "<<delta<<endl;
01895                 double x1=(2*x_cirtemp+2*k*y_cirtemp-2*k*b+sqrt(delta))/(2*(k*k+1));
01896                 double x2=(2*x_cirtemp+2*k*y_cirtemp-2*k*b-sqrt(delta))/(2*(k*k+1));
01897                 double x_middle=(vec_x_east[i]+vec_x_west[i])/2.;
01898                 //        if(abs(x1-vec_x[i])>abs(x2-vec_x[i])) {x_stereo=x2;} 
01899                 if(abs(x1-x_middle)>abs(x2-x_middle)) {x_stereo=x2;} 
01900                 else {x_stereo=x1;}
01901                 y_stereo=k*x_stereo+b;
01902           }
01903           else{
01904                 //cout<<"layer:  "<<vec_layer[i]<<"wire:  "<<vec_wire[i]<<endl;
01905                 //cout<<"x_east:  "<<vec_x_east[i]<<"  "<<"x_west:   "<<vec_x_west[i]<<endl;
01906                 x_stereo=vec_x_east[i];
01907                 double y1=sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo))+y_cirtemp;
01908                 double y2=y_cirtemp-sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo));
01909                 double y_middle=(vec_y_east[i]+vec_y_west[i])/2.;
01910                 if(abs(y1-y_middle)>abs(y2-y_middle)) {y_stereo=y2;}
01911                 else{y_stereo=y1;}
01912           }
01913           //      cout<<"layer:  "<<vec_layer[i]<<"wire:  "<<vec_wire[i]<<"x_stereo:  "<<x_stereo<<"  "<<"y_stereo:  "<<y_stereo<<endl;
01914 
01915           double theta_temp;
01916           if(x_cirtemp==0||x_stereo-x_cirtemp==0){
01917                 theta_temp=0;   
01918           }
01919           else{
01920                 theta_temp=M_PI-atan2(y_stereo-y_cirtemp,x_stereo-x_cirtemp)+atan2(y_cirtemp,x_cirtemp);
01921                 if(theta_temp>2*M_PI){
01922                   theta_temp=theta_temp-2*M_PI;
01923                 }
01924                 if(theta_temp<0){
01925                   theta_temp=theta_temp+2*M_PI;
01926                 }
01927           }
01928           //cout<<"thetatemp:   "<<theta_temp<<endl;
01929           //if(theta_temp>2*M_PI) {
01930           //  theta_temp=theta_temp-2*M_PI;
01931           //}
01932           //if(theta_temp<0.) {
01933           //  theta_temp=theta_temp+2*M_PI;
01934           //}
01935           //cout<<"debug: if thetatemp is cout over"<<endl;
01936           l.push_back(r_temp*theta_temp);
01937           //cout<<"debug: if l is pushback"<<endl;
01938           //l.push_back(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp));
01939           //l.push_back(r_temp*(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp)));
01940 
01941           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]));
01942           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]));
01943           vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
01944           //cout<<"debug: if vec_zStereo is pushback"<<endl;
01945           //if(stereoId==0) {vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);}
01946           //if(stereoId>0){
01947           //  if(vec_layer[i]==vec_layer.at(i-1)) {
01948           //    vec_zStereo[(stereoId-1)]=((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.);
01949           //    vec_zStereo.push_back((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.);
01950           //  }
01951           //  else{
01952           //    vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
01953           //  }
01954           //}
01955 
01956           //      double delphi=atan2(y_stereo,x_stereo)-par.phi0();
01957           //cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<" "<<m_cell[i]<<" "<<"z_stereo: "<<z_stereo<<"  "<<d1<<"   "<<d2<<"    "<<endl;
01958           //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;
01959           //      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;
01960 
01961           //      cout<<"layer:  "<<vec_layer[i]<<"wire:  "<<vec_wire[i]<<"z:  "<<vec_zStereo[stereoId]<<"  "<<"l:  "<<l[stereoId]<<"  "<<vec_truthHit[i]<<endl;
01962           stereoId=stereoId+1;
01963           //cout<<"theta:  "<<theta<<endl;
01964         }       
01965         else {k=b=0;}
01966   }
01967   cout<<"nDropDelta:  "<<nDropDelta<<endl;
01968   //  cout<<"if for is finished :  "<<endl;
01969   //int digiId=0;
01970   //for(int i=0;i<n;i++){
01971   //  if(vec_slant[i]!=0){
01972   //    cout<<"layerid: "<<"   "<<vec_layer[i]<<"   "<<"cellid: "<<vec_wire[i]<<"    "<<"zstereo:  "<<z_stereo.at(digiId)<<"   "<<vec_zStereo.at(digiId)<<"    "<<"ztruth:  "<<vec_z[i]<<endl;
01973   //    digiId++;
01974   //  } 
01975   //}
01976   //double sum_zstereo=0;
01977   //for(int i=0;i<stereoId;i++){
01978   //  sum_zstereo=sum_zstereo+z_stereo[i];
01979   //}
01980   //cout<<"sum_zstereo:  "<<sum_zstereo<<endl;
01981   return stereoId;
01982 }

int MdcHoughFinder::binX [private]

Definition at line 73 of file MdcHoughFinder.h.

Referenced by execute().

int MdcHoughFinder::binY [private]

Definition at line 74 of file MdcHoughFinder.h.

Referenced by execute().

double MdcHoughFinder::dz_mc [private]

Definition at line 76 of file MdcHoughFinder.h.

Referenced by execute(), and GetMcInfo().

NTuple::Array<int> MdcHoughFinder::m_2d_nFit [private]

Definition at line 187 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<int> MdcHoughFinder::m_3d_nFit [private]

Definition at line 188 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_areaLeast [private]

Definition at line 156 of file MdcHoughFinder.h.

Referenced by initialize().

NTuple::Item<double> MdcHoughFinder::m_areaLeastNum [private]

Definition at line 157 of file MdcHoughFinder.h.

Referenced by initialize().

NTuple::Item<double> MdcHoughFinder::m_areaSelectHit [private]

Definition at line 158 of file MdcHoughFinder.h.

Referenced by initialize().

NTuple::Item<double> MdcHoughFinder::m_areaSelectHit_signal [private]

Definition at line 159 of file MdcHoughFinder.h.

Referenced by initialize().

BField* MdcHoughFinder::m_bfield [private]

Definition at line 91 of file MdcHoughFinder.h.

Referenced by initialize().

int MdcHoughFinder::m_binx [private]

Definition at line 59 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

int MdcHoughFinder::m_biny [private]

Definition at line 60 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

double MdcHoughFinder::m_bunchT0 [private]

Definition at line 68 of file MdcHoughFinder.h.

Referenced by digiToHots(), digiToHots2(), and execute().

NTuple::Array<int> MdcHoughFinder::m_cell [private]

Definition at line 111 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 113 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

bool MdcHoughFinder::m_combineTracking [private]

Definition at line 103 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

std::string MdcHoughFinder::m_configFile [private]

Definition at line 87 of file MdcHoughFinder.h.

TrkContextEv* MdcHoughFinder::m_context [private]

Definition at line 90 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_cosCut [private]

Definition at line 148 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_costaTruth [private]

Definition at line 201 of file MdcHoughFinder.h.

Referenced by execute(), GetMcInfo(), and initialize().

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

Definition at line 168 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

int MdcHoughFinder::m_data [private]

Definition at line 58 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

int MdcHoughFinder::m_debug [private]

Definition at line 57 of file MdcHoughFinder.h.

Referenced by digiToHots(), digiToHots2(), execute(), FillCells(), GetMcInfo(), initialize(), LeastSquare(), and MdcHoughFinder().

bool MdcHoughFinder::m_dropHot [private]

Definition at line 83 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

NTuple::Item<double> MdcHoughFinder::m_drTruth [private]

Definition at line 203 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_dzTruth [private]

Definition at line 204 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

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

Definition at line 140 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

double MdcHoughFinder::m_fpro [private]

Definition at line 63 of file MdcHoughFinder.h.

Referenced by MdcHoughFinder().

uint32_t MdcHoughFinder::m_getDigiFlag [private]

Definition at line 80 of file MdcHoughFinder.h.

Referenced by MdcHoughFinder().

const MdcDetector* MdcHoughFinder::m_gm [private]

Definition at line 88 of file MdcHoughFinder.h.

Referenced by beginRun(), digiToHots(), digiToHots2(), and execute().

std::vector<float> MdcHoughFinder::m_helixHitsSigma [private]

Definition at line 104 of file MdcHoughFinder.h.

double MdcHoughFinder::m_hit_pro [private]

Definition at line 64 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

NTuple::Array<int> MdcHoughFinder::m_hitCol [private]

Definition at line 107 of file MdcHoughFinder.h.

NTuple::Array<int> MdcHoughFinder::m_hitSignal [private]

Definition at line 109 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

bool MdcHoughFinder::m_keepBadTdc [private]

Definition at line 82 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

bool MdcHoughFinder::m_keepUnmatch [private]

Definition at line 84 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

NTuple::Array<double> MdcHoughFinder::m_l [private]

Definition at line 194 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 110 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 112 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 108 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_maxCount [private]

Definition at line 150 of file MdcHoughFinder.h.

Referenced by initialize().

int MdcHoughFinder::m_maxMdcDigi [private]

Definition at line 81 of file MdcHoughFinder.h.

Referenced by MdcHoughFinder().

const MdcCalibFunSvc* MdcHoughFinder::m_mdcCalibFunSvc [private]

Definition at line 92 of file MdcHoughFinder.h.

Referenced by digiToHots(), digiToHots2(), and initialize().

MdcGeomSvc* MdcHoughFinder::m_mdcGeomSvc [private]

Definition at line 95 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

int MdcHoughFinder::m_method [private]

Definition at line 56 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

int MdcHoughFinder::m_minMdcDigi [private]

Definition at line 85 of file MdcHoughFinder.h.

Referenced by MdcHoughFinder().

NTuple::Item<int> MdcHoughFinder::m_nCross [private]

Definition at line 142 of file MdcHoughFinder.h.

Referenced by initialize(), and RhoTheta().

int MdcHoughFinder::m_ndev [private]

Definition at line 62 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

int MdcHoughFinder::m_nEvtSuccess [private]

Definition at line 97 of file MdcHoughFinder.h.

Referenced by execute().

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

Definition at line 179 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_nFitSucess [private]

Definition at line 190 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 143 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 144 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 181 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_nHitAxialSignal [private]

Definition at line 182 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<int> MdcHoughFinder::m_nHitAxialSignal_select [private]

Definition at line 186 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 183 of file MdcHoughFinder.h.

NTuple::Item<int> MdcHoughFinder::m_nHitSignal [private]

Definition at line 180 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<int> MdcHoughFinder::m_nHitSignal_select [private]

Definition at line 184 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_npeak [private]

Definition at line 152 of file MdcHoughFinder.h.

Referenced by initialize().

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

Definition at line 170 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 127 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 89 of file MdcHoughFinder.h.

Referenced by initialize().

std::string MdcHoughFinder::m_pdtFile [private]

Definition at line 65 of file MdcHoughFinder.h.

Referenced by initialize(), and MdcHoughFinder().

NTuple::Array<double> MdcHoughFinder::m_peakArea [private]

Definition at line 155 of file MdcHoughFinder.h.

Referenced by initialize().

std::vector<float> MdcHoughFinder::m_peakCellNum [private]

Definition at line 61 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

NTuple::Array<double> MdcHoughFinder::m_peakHigh [private]

Definition at line 154 of file MdcHoughFinder.h.

Referenced by initialize().

NTuple::Array<double> MdcHoughFinder::m_peakWidth [private]

Definition at line 153 of file MdcHoughFinder.h.

Referenced by initialize().

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

Definition at line 169 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_phi0Truth [private]

Definition at line 202 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

bool MdcHoughFinder::m_pickHits [private]

Definition at line 66 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

int MdcHoughFinder::m_pid [private]

Definition at line 101 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and MdcHoughFinder().

NTuple::Item<double> MdcHoughFinder::m_pidTruth [private]

Definition at line 200 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

IMagneticFieldSvc* MdcHoughFinder::m_pIMF [private]

Definition at line 98 of file MdcHoughFinder.h.

Referenced by initialize().

NTuple::Array<double> MdcHoughFinder::m_pt [private]

Definition at line 174 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 175 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_pTruth [private]

Definition at line 207 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_ptTruth [private]

Definition at line 205 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_pxyz [private]

Definition at line 177 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_pz [private]

Definition at line 176 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_pzTruth [private]

Definition at line 206 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_qTruth [private]

Definition at line 208 of file MdcHoughFinder.h.

Referenced by GetMcInfo(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_r [private]

Definition at line 164 of file MdcHoughFinder.h.

Referenced by initialize(), and LeastSquare().

RawDataProviderSvc* MdcHoughFinder::m_rawDataProviderSvc [private]

Definition at line 94 of file MdcHoughFinder.h.

Referenced by digiToHots(), digiToHots2(), execute(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_rho [private]

Definition at line 125 of file MdcHoughFinder.h.

Referenced by initialize(), and RhoTheta().

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

Definition at line 141 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 128 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 172 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_theta [private]

Definition at line 126 of file MdcHoughFinder.h.

Referenced by initialize(), and RhoTheta().

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

Definition at line 167 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 166 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

int MdcHoughFinder::m_trackNum_Mc_set [private]

Definition at line 55 of file MdcHoughFinder.h.

Referenced by execute(), and MdcHoughFinder().

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

Definition at line 197 of file MdcHoughFinder.h.

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

Definition at line 123 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_v [private]

Definition at line 124 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 120 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_x_circle [private]

Definition at line 162 of file MdcHoughFinder.h.

Referenced by initialize(), and LeastSquare().

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

Definition at line 114 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 116 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 121 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MdcHoughFinder::m_y_circle [private]

Definition at line 163 of file MdcHoughFinder.h.

Referenced by initialize(), and LeastSquare().

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

Definition at line 115 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 117 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 122 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 171 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 118 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

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

Definition at line 119 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Array<double> MdcHoughFinder::m_zStereo [private]

Definition at line 193 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<int> MdcHoughFinder::m_zStereoNum [private]

Definition at line 192 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

int MdcHoughFinder::nfailure [private]

Definition at line 102 of file MdcHoughFinder.h.

NTuple::Tuple* MdcHoughFinder::ntuplehit [private]

Definition at line 106 of file MdcHoughFinder.h.

Referenced by execute(), and initialize().

int MdcHoughFinder::t_eventNo [private]

Definition at line 96 of file MdcHoughFinder.h.

int MdcHoughFinder::t_eventNum [private]

Definition at line 69 of file MdcHoughFinder.h.

Referenced by execute(), and GetMcInfo().

double MdcHoughFinder::t_maxP [private]

Definition at line 71 of file MdcHoughFinder.h.

Referenced by execute(), FillCells(), FillHist(), and FillRhoTheta().

double MdcHoughFinder::t_minP [private]

Definition at line 72 of file MdcHoughFinder.h.

Referenced by execute().

double MdcHoughFinder::t_nTrkMC [private]

Definition at line 100 of file MdcHoughFinder.h.

Referenced by GetMcInfo().

int MdcHoughFinder::t_runNum [private]

Definition at line 70 of file MdcHoughFinder.h.

Referenced by execute(), and GetMcInfo().

double MdcHoughFinder::t_t0Truth [private]

Definition at line 99 of file MdcHoughFinder.h.

double MdcHoughFinder::tanl_mc [private]

Definition at line 77 of file MdcHoughFinder.h.

Referenced by execute(), and GetMcInfo().

int MdcHoughFinder::track_fit [private]

Definition at line 78 of file MdcHoughFinder.h.

Referenced by execute().

MdcHoughFinder Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	MdcHoughFinder (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	beginRun ()
Private Member Functions
int	GetMcInfo ()
int	digiToHots (TrkRecoTrk *newTrack, vector< bool > vec_truthHit)
int	digiToHots2 (TrkRecoTrk *newTrack, vector< bool > vec_truthHit)
double	CFtrans (double x, double y)
int	SlantId (int layer)
int	LeastSquare (int n, int nselecthit_axial, 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)
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 numCross, int m_nHit, vector< double > vec_u, vector< double > vec_v, vector< double > &vec_rho, vector< double > &vec_theta, vector< vector< int > > &vec_hitNum)
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 numCross)
Private Attributes
int	m_trackNum_Mc_set
int	m_method
int	m_debug
int	m_data
int	m_binx
int	m_biny
std::vector< float >	m_peakCellNum
int	m_ndev
double	m_fpro
double	m_hit_pro
std::string	m_pdtFile
bool	m_pickHits
double	m_bunchT0
int	t_eventNum
int	t_runNum
double	t_maxP
double	t_minP
int	binX
int	binY
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
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
int	t_eventNo
int	m_nEvtSuccess
IMagneticFieldSvc *	m_pIMF
double	t_t0Truth
double	t_nTrkMC
int	m_pid
int	nfailure
bool	m_combineTracking
std::vector< float >	m_helixHitsSigma
NTuple::Tuple *	ntuplehit
NTuple::Array< int >	m_hitCol
NTuple::Array< int >	m_layerNhit
NTuple::Array< int >	m_hitSignal
NTuple::Array< int >	m_layer
NTuple::Array< int >	m_cell
NTuple::Array< int >	m_layer_Mc
NTuple::Array< int >	m_cell_Mc
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_rho
NTuple::Array< double >	m_theta
NTuple::Array< double >	m_p
NTuple::Array< int >	m_slant
NTuple::Item< int >	m_eventNum
NTuple::Item< int >	m_runNum
NTuple::Item< int >	m_nCross
NTuple::Item< int >	m_nHit
NTuple::Item< int >	m_nHit_Mc
NTuple::Item< int >	m_cosCut
NTuple::Item< int >	m_maxCount
NTuple::Item< int >	m_npeak
NTuple::Array< double >	m_peakWidth
NTuple::Array< double >	m_peakHigh
NTuple::Array< double >	m_peakArea
NTuple::Item< double >	m_areaLeast
NTuple::Item< double >	m_areaLeastNum
NTuple::Item< double >	m_areaSelectHit
NTuple::Item< double >	m_areaSelectHit_signal
NTuple::Item< double >	m_x_circle
NTuple::Item< double >	m_y_circle
NTuple::Item< double >	m_r
NTuple::Item< int >	m_trackNum_Mc
NTuple::Item< int >	m_trackNum
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::Item< int >	m_nHitSignal
NTuple::Item< int >	m_nHitAxial
NTuple::Item< int >	m_nHitAxialSignal
NTuple::Array< int >	m_nHitSelect
NTuple::Array< int >	m_nHitSignal_select
NTuple::Array< int >	m_nHitAxialSignal_select
NTuple::Array< int >	m_2d_nFit
NTuple::Array< int >	m_3d_nFit
NTuple::Item< int >	m_nFitSucess
NTuple::Item< int >	m_zStereoNum
NTuple::Array< double >	m_zStereo
NTuple::Array< double >	m_l
double	m_truthPos [43][288][3]
NTuple::Item< double >	m_pidTruth
NTuple::Item< double >	m_costaTruth
NTuple::Item< double >	m_phi0Truth
NTuple::Item< double >	m_drTruth
NTuple::Item< double >	m_dzTruth
NTuple::Item< double >	m_ptTruth
NTuple::Item< double >	m_pzTruth
NTuple::Item< double >	m_pTruth
NTuple::Item< double >	m_qTruth