MdcHoughFinder Class Reference

#include <MdcHoughFinder.h>

Inheritance diagram for MdcHoughFinder:

Public Member Functions
	MdcHoughFinder (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	beginRun ()
	MdcHoughFinder (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	beginRun ()

Detailed Description

Definition at line 31 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/MdcHoughFinder.h.

Constructor & Destructor Documentation

MdcHoughFinder() [1/2]

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

Definition at line 56 of file MdcHoughFinder.cxx.

                                                                              :
  Algorithm(name, pSvcLocator)
{
  // Declare the properties  
  declareProperty("trackNumMc", m_trackNum_Mc_set=1);
  declareProperty("method", m_method=0);
  declareProperty("debug",          m_debug = 0);
  declareProperty("data",          m_data= 0);
  declareProperty("binx",          m_binx= 100);
  declareProperty("biny",          m_biny= 100);
  declareProperty("peakCellNum",          m_peakCellNum);
  declareProperty("ndev",          m_ndev= 3);
  declareProperty("f_pro",          m_fpro= 0.8);
  declareProperty("f_hit_pro",          m_hit_pro= 0.8);
  //declareProperty("debug",          m_debug = 0);
  declareProperty("pdtFile",        m_pdtFile = "pdt.table");
  declareProperty("pickHits",       m_pickHits = true);
  declareProperty("pid",        m_pid = 0);
  //declareProperty("helixHitsSigma", m_helixHitsSigma);
  declareProperty("combineTracking",m_combineTracking = true);
  declareProperty("getDigiFlag",    m_getDigiFlag = 0);
  declareProperty("maxMdcDigi",     m_maxMdcDigi = 0);
  declareProperty("keepBadTdc",     m_keepBadTdc = 0);
  declareProperty("dropHot",        m_dropHot= 0);
  declareProperty("keepUnmatch",    m_keepUnmatch = 0);
  declareProperty("minMdcDigi",     m_minMdcDigi = 0);
 
}

MdcHoughFinder() [2/2]

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

Member Function Documentation

beginRun() [1/2]

StatusCode MdcHoughFinder::beginRun

(

)

Definition at line 85 of file MdcHoughFinder.cxx.

                                   {
 
  //Initailize MdcDetector
  m_gm = MdcDetector::instance(0);
  if(NULL == m_gm) return StatusCode::FAILURE;
 
 
  return StatusCode::SUCCESS;
}

beginRun() [2/2]

StatusCode MdcHoughFinder::beginRun

(

)

execute() [1/2]

StatusCode MdcHoughFinder::execute

(

)

Definition at line 275 of file MdcHoughFinder.cxx.

                                   {
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
 
  //event start time
  SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
  if (!evTimeCol) {
    log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
    m_bunchT0=0.;
  }else{
    RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
    if (iter_evt != evTimeCol->end()){
      m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
    }
  }
 
  //if(m_debug)TrkHelixFitter::m_debug = true;
  // Get the event header, print out event and run number
 
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  if (!eventHeader) {
    log << MSG::FATAL << "Could not find Event Header" << endreq;
    return( StatusCode::FAILURE);
  }
 
 
  DataObject *aTrackCol;
  DataObject *aRecHitCol;
  SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
  if(!m_combineTracking){
    eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
    if(aTrackCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
    }
    eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
    if(aRecHitCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
    }
  }
 
  RecMdcTrackCol* trackList;
  eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
  if (aTrackCol) {
    trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
  }else{
    trackList = new RecMdcTrackCol;
    StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
  int nTdsTk = (int) trackList->size();
  int nFitSucess = 0;
 
  RecMdcHitCol* hitList;
  eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
  if (aRecHitCol) {
    hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
  }else{
    hitList = new RecMdcHitCol;
    StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
 
  //pick the first layer
 
  // m_nLayerNum=43;
  // for(int i=0;i<43;i++){
  //     const MdcGeoWire *geowir_first_layer =m_mdcGeomSvc->Wire(i,1);
  //     HepPoint3D eastP = geowir_first_layer->Backward()/10.0;//mm 2 cm
  //     HepPoint3D westP = geowir_first_layer->Forward() /10.0;//mm 2 cm
  //     m_layer1X[i]=(eastP.x()+westP.x())/2;
  //     m_layer1Y[i]=(eastP.y()+westP.y())/2;
  //     m_layer1U[i]=CFtrans((eastP.x()+westP.x())/2.,(eastP.y()+westP.y())/2.);
  //     m_layer1V[i]=CFtrans((eastP.y()+westP.y())/2.,(eastP.x()+westP.x())/2.);
  //     
  //     cout<<"("<<i<<","<<"1"<<")"<<":"<<"("<<m_layer1X[i]<<","<<m_layer1Y[i]<<")"<<endl;
  // }
 
  t_eventNum=eventHeader->eventNumber();
  t_runNum=eventHeader->runNumber();
  m_eventNum=t_eventNum;
  m_runNum=t_runNum;
  log << MSG::INFO << "MdcHoughFinder: retrieved event: " << eventHeader->eventNumber()  << "  run: " <<  eventHeader->runNumber() << endreq;
 
  // get mc particle digi
  int mc_particle=GetMcInfo();
 
 
  if(abs(m_costaTruth)<=0.83){m_cosCut=1;}
  else{m_cosCut=0;}
 
  //int n_curve=0;
  //double vec_ptTruth=m_ptTruth;
  //double n_costaTruth=m_costaTruth;
  //double z_flight=2*M_PI*(771.4/2)*(vec_ptTruth/tan(n_costaTruth));
  //if(vec_ptTruth<0.12&&(z_flight<1193.)){
  ///// n_curve++;
  ////  m_curve=n_curve;
  //}
  vector<double> vec_u;
  vector<double> vec_v;
  vector<double> vec_p;
  vector<double> vec_x_east;
  vector<double> vec_y_east;
  vector<double> vec_z_east;
  vector<double> vec_x_west;
  vector<double> vec_y_west;
  vector<double> vec_z_west;
  vector<double> vec_z_Mc;
  vector<double> vec_x;
  vector<double> vec_y;
  vector<double> vec_z;
  vector<int> vec_layer;
  vector<int> vec_wire;
  vector<int> vec_slant;
  vector<double> vec_zStereo;
  vector<double> l;
 
  vector<int> vec_layer_Mc;
  vector<int> vec_wire_Mc;
  t_maxP=-999.;
  t_minP=999.;
 
  // retrieve Mdc truth hits
  int t_nHit_Mc;
  int digiId_Mc=0;
  int nHitAxial_Mc=0;
 
  SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol"); 
  if (!mcMdcMcHitCol) {
    log << MSG::INFO << "Could not find MdcMcHit" << endreq; 
  }else{
    Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin();
 
    for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) {
      Identifier mdcid = (*iterMcHit)->identify();
      int layerId_Mc = MdcID::layer(mdcid);
      int wireId_Mc = MdcID::wire(mdcid);
 
      vec_layer_Mc.push_back(layerId_Mc);
      vec_wire_Mc.push_back(wireId_Mc);
 
      //      cout<<"mcHit:  "<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<endl;
      //vec_layer.push_back(layerId_Mc);
      //vec_wire.push_back(wireId_Mc);
      m_layer_Mc[digiId_Mc]=layerId_Mc;
      m_cell_Mc[digiId_Mc]=wireId_Mc;
 
      if(m_data==0){
        if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) {    
          nHitAxial_Mc++;}
        if(m_debug>0) {cout<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<"nHitAxial_Mc:  "<<nHitAxial_Mc<<endl;}
        m_layer[digiId_Mc]=layerId_Mc;
        m_cell[digiId_Mc]=wireId_Mc;
        const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc);
        HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
        HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
 
        m_x_east[digiId_Mc]=eastP.x();
        m_y_east[digiId_Mc]=eastP.y();
        m_z_east[digiId_Mc]=eastP.z();
        m_x_west[digiId_Mc]=westP.x();
        m_y_west[digiId_Mc]=westP.y();
        m_z_west[digiId_Mc]=westP.z();
 
        vec_x_east.push_back(eastP.x());
        vec_y_east.push_back(eastP.y());
        vec_z_east.push_back(eastP.z());
        vec_x_west.push_back(westP.x());
        vec_y_west.push_back(westP.y());
        vec_z_west.push_back(westP.z());
 
        double x = (*iterMcHit)->getPositionX()/10.;//mm 2 cm
        double y = (*iterMcHit)->getPositionY()/10.;//mm 2 cm
        double z = (*iterMcHit)->getPositionZ()/10.;//mm 2 cm
 
        double u=CFtrans(x,y);
        double v=CFtrans(y,x);
        double p=sqrt(u*u+v*v);
 
        vec_x.push_back(x);
        vec_y.push_back(y);
        vec_z.push_back(z);
        vec_u.push_back(u);
        vec_v.push_back(v);
        vec_p.push_back(p);
 
        m_x[digiId_Mc]=x;
        m_y[digiId_Mc]=y;
        m_z[digiId_Mc]=z;
        m_u[digiId_Mc]=u;
        m_v[digiId_Mc]=v;
        m_p[digiId_Mc]=p;
 
        int layer= layerId_Mc;
        vec_slant.push_back(SlantId(layer));
        m_slant[digiId_Mc]=SlantId(layer);
        if (m_slant!=0)
        {cout<<"layer:  "<<layerId_Mc<<"wire:  "<<wireId_Mc<<"x_east:  "<<m_x_east[digiId_Mc]<<"y_east:  "<<m_y_east[digiId_Mc]<<endl;}
      }
    }
  }
 
  t_nHit_Mc=digiId_Mc;
  m_nHit_Mc=digiId_Mc;
  if(m_data==0) {m_nHit=digiId_Mc;}
 
  //m_nHitAxial=nHitAxial_Mc;
  //    m_nFitFailure=0;
 
  //    for(int i =0;i<t_nHit_Mc;i++){
  //      if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];} 
  //    }
  //cout<<"t_maxP:"<<t_maxP<<"  "<<endl;
 
  //    for(int i =0;i<t_nHit_Mc;i++){
  //      if(t_minP>vec_p[i]) {t_minP=vec_p[i];} 
  //    }
  //cout<<"t_minP:"<<t_minP<<"  "<<endl;
 
 
  //  retrieve Mdc digi vector form RawDataProviderSvc 
  vector<int> vec_hitSignal;
  vector<int> vec_track_index;
  int track_index_max=0;
  uint32_t getDigiFlag = 0;
  if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
  if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
  if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
  //  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
  cout<<"MdcDigiVec:  "<<mdcDigiVec.size()<<endl;
  MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
  int t_nHit;
  int digiId = 0;
  int nHitAxial = 0;
  int nHitLayer[43];
  int nHitSignal=0;
  int nHitAxialSignal=0;
 
  for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
    int track_index=(*iter1)->getTrackIndex();
//  cout<<"mc track_index: "<<track_index<<endl;
    if(track_index>=1000) track_index-=1000;
    vec_track_index.push_back(track_index);
    if(track_index>=0){
      nHitSignal++;
      m_hitSignal[digiId]=1;
    }
    //m_hitCol[digiId]=digiId; 
    Identifier mdcId= (*iter1)->identify();
    int layerId = MdcID::layer(mdcId);
    int wireId = MdcID::wire(mdcId);
 
 
    //cout<<"layer:  "<<layerId<<"  "<<"wire:  "<<wireId<<"hitSignal:  "<<vec_hitSignal[digiId]<<endl;
 
    if ((layerId>=8&&layerId<=19)||(layerId>=36)) { 
      nHitAxial++;}
    if (((layerId>=8&&layerId<=19)||(layerId>=36))&&track_index>=0) {   
      nHitAxialSignal++;}
 
    vec_layer.push_back(layerId);
    vec_wire.push_back(wireId);
 
    nHitLayer[layerId]++;
    m_layerNhit[layerId]=nHitLayer[layerId];    
 
    const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId);
    HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
    HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
 
    vec_x_east.push_back(eastP.x());
    vec_y_east.push_back(eastP.y());
    vec_z_east.push_back(eastP.z());
    vec_x_west.push_back(westP.x());
    vec_y_west.push_back(westP.y());
    vec_z_west.push_back(westP.z());
 
    m_x_east[digiId]=eastP.x();
    m_y_east[digiId]=eastP.y();
    m_z_east[digiId]=eastP.z();
    m_x_west[digiId]=westP.x();
    m_y_west[digiId]=westP.y();
    m_z_west[digiId]=westP.z();
    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;}
    //conformal trans:
    double x =(eastP.x()+westP.x())/2.;
    double y =(eastP.y()+westP.y())/2.;
    vec_x.push_back((vec_x_east[digiId]+vec_x_west[digiId])/2); 
    vec_y.push_back((vec_y_east[digiId]+vec_y_west[digiId])/2); 
    m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2; 
    m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2; 
    double u=CFtrans(x,y);
    double v=CFtrans(y,x);
    //  cout<< "u:   "<<u<<"v:   "<<v<<"digiId:  "<<digiId<<endl;
    //cout<<"layer:  "<<layerId<<" "<<"wireId:  "<<wireId<<" "<<"x:      "<<x<<"   "<<"y:     "<<y<<"     "<<endl;
    vec_p.push_back(sqrt(u*u+v*v));
    vec_u.push_back(u);
    vec_v.push_back(v);
    m_u[digiId]=u;
    m_v[digiId]=v;
    m_p[digiId]=sqrt(u*u+v*v);
 
    m_layer[digiId]=geowir->Layer();
    m_cell[digiId]=geowir->Cell();
    int layer= layerId;
    vec_slant.push_back(SlantId(layer));
    m_slant[digiId]=SlantId(layer);
    //cout<<"layer:  "<<layer<<"  "<<"slant:  "<<m_slant[digiId]<<endl;
    //   if (m_slant!=0)
    //  {cout<<"layer:  "<<layerId<<"wire:  "<<wireId<<"x_east:  "<<m_x_east[digiId]<<"y_east:  "<<m_y_east[digiId]<<endl;}
  }
  for(int i=0;i<digiId;i++){
    if(track_index_max<=vec_track_index[i]) track_index_max=vec_track_index[i];
  }
  track_index_max++;
  m_trackNum_Mc=track_index_max;
//  cout<<"mc truth has :"<<track_index_max<<"  mc track."<<endl;
  vector< vector <int> > mc_track_hit(track_index_max,vector<int>() );
  for(int i=0;i<track_index_max;i++){ 
    for(int j=0;j<digiId;j++){
      if(vec_track_index[j]==i) mc_track_hit[i].push_back(j);
    }
  }
  //composory set mc_track number
  track_index_max=m_trackNum_Mc_set;
  //debug mc track hit
  //for(int i=0;i<track_index_max;i++){ 
  //  for(int j=0;j<mc_track_hit[i].size();j++){
  //    int hitId=mc_track_hit[i][j];
  //    cout<<"track: "<<i<<"  has hit:  "<<mc_track_hit[i][j]<<"("<<vec_layer[hitId]<<","<<vec_wire[hitId]<<")"<<endl;
  //  }
  //}
 
  //t_nHit=digiId;
  m_nHit=digiId;
  m_nHitAxial=nHitAxial;
  //m_nFitFailure=0;
 
  m_nHitSignal=nHitSignal;
  m_nHitAxialSignal=nHitAxialSignal;
 
  cout<<"hit number:  "<<digiId<<endl;
 
 
  // get the max&min amplitude of the curve 
  for(int i =0;i<m_nHit;i++){
    if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];} 
  }
 
  for(int i =0;i<m_nHit;i++){
    if(t_minP>vec_p[i]) {t_minP=vec_p[i];} 
  }
  t_maxP=t_maxP+0.01;
 
 
  //double peakArea[100];
  //m_npeak=100;
  //for( int bin_i=0;bin_i<10;bin_i++){
  //    for(int bin_j=0;bin_j<10;bin_j++){
  //      int binx=m_binx+100*bin_i;
  //      int biny=m_biny+100*bin_j;
  int nhit;
  if(m_data==0) {
    nhit=m_nHit_Mc;}
  else{
    nhit=m_nHit;}
 
  binX=m_binx;
  binY=m_biny;
  double binwidth=M_PI/binX;
  double binhigh=2*t_maxP/binY;
  TH2D *h1=new TH2D("h1","h1",binX,0,M_PI,binY,-t_maxP,t_maxP);
 
  //method 0
  //2-point combine
  vector<double> vec_rho;
  vector<double> vec_theta;
  vector< vector<int> > vec_hitNum(2,vector<int>());      //store 2 hits in a cross point
  int numCross=(int)(nhit*(nhit-1)/2);
  if(m_method==0){
    RhoTheta(numCross,nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum);    //calculate cross point
    FillRhoTheta(h1,vec_theta,vec_rho,numCross);         //fill histgram method by rhotheta
  }
  // method  1
  vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) ); 
  if(m_method==1){
    FillHist(h1,vec_u,vec_v,nhit,vec_selectNum);
  }
 
 
  vector<bool> vec_truthHit(nhit,false);
  vector< vector <int> > countij(binX,vector <int> (binY,0) ); 
  //  vector< vector < vector<int> > > vec_selectNum(binX,vector < vector<int> >(binY,vector<int>() )); 
  vector< vector< vector<int> > > vec_selectHit(binX,vector< vector<int> >(binY,vector<int>() ) ); 
  if(m_method==2){
    FillCells(h1,nhit,binX,binY,vec_u,vec_v,vec_layer,vec_wire,countij,vec_selectNum,vec_selectHit);
  }
 
  //cout h1 
  //for(int i=0;i<binX;i++){
  //  for(int j=0;j<binY;j++){
  //    int count =h1->GetBinContent(i+1,j+1);
  ////    for(int k=0;k<count;k++){
  //    cout<<"("<<i<<","<<j<<")"<<":  "<<count<<endl;
  ////    }
  //  }
  //}
 
 
  //int nHitSelect=0;
  //int nHitSignal_select=0;
  //int nHitAxialSignal_select=0;
  //find peak
  int max_count=0;
  int max_binx=1;
  int max_biny=1;
  for(int i=0;i<binX;i++){
    for(int j=0;j<binY;j++){
      int count=h1->GetBinContent(i+1,j+1);
      if(max_count<count) {
        max_count=count;
        max_binx=i+1;
        max_biny=j+1;
      }
    }
  }
  //        cout<<"("<<max_binx<<","<<max_biny<<","<<max_count<<")"<<"numCross:  "<<numCross<<endl;
  //  h1->Draw("lego2");
 
  //int peak_cellNum=1;         // to determine how many cell should be embraced in one orientation
  //int count_sum[1000];
  //for(int i=0;i<1000;i++){
  //  count_sum[i]=0;
  //  for(int alphax=max_binx-i;alphax<=max_binx+i;alphax++){
  //    for(int rhox=max_biny-i;rhox<=max_biny+i;rhox++){
  //    int ax;
  //    if (alphax<0) { ax=alphax+binX; }
  //    else if (alphax>=binX) { ax=alphax-binX; }
  //    else { ax=alphax; }
  //    count_sum[i]+=h1->GetBinContent(ax,rhox);
  //    }
  //  }
  //  if(count_sum[i]>m_fpro*numCross){
  //    peak_cellNum=i;
  //              cout<<"i: "<<i<<endl;
  //            cout<<"count_sum:  "<<count_sum[i]<<endl;
  //            cout<<"pro: "<<(double)count_sum[i]/(double)numCross<<endl;
  //    break;
  //  }
  //}
 
 
  /*
  //return hit number
  vector<bool> vec_hitSelect(m_nHit,false);
  int bin_left=max_binx-peak_cellNum;
  int bin_right=max_binx+peak_cellNum;
  int bin_up=max_biny+peak_cellNum;
  int bin_down=max_biny-peak_cellNum;
  double area_left=(bin_left-1)*binwidth;
  double area_right=(bin_right)*binwidth;
  double area_down=(bin_down-1)*binhigh-t_maxP;
  double area_up=(bin_up)*binhigh-t_maxP;
  // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
  double peak_width=(2*peak_cellNum+1)*binwidth;
  double peak_high=(2*peak_cellNum+1)*binhigh;
  cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
  //test if the hit is in this area
  for(int i=0;i<numCross;i++){
  if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
  int hitNum_1=vec_hitNum[0][i];
  int hitNum_2=vec_hitNum[1][i];
  vec_hitSelect[hitNum_1]=true;
  vec_hitSelect[hitNum_2]=true;
  }
  }
  int selectHitNum=0;
  for(int i=0;i<m_nHit;i++){
  if(vec_hitSelect[i]==1) {selectHitNum++;}
  //    cout<<"if hit is select: "<<i<<":  "<<vec_hitSelect[i]<<endl;
  }
  cout<<"how many Hit have been selected: "<<selectHitNum<<"  "<<(double)selectHitNum/(double)m_nHit;
  m_areaSelectHit=(double)selectHitNum/(double)m_nHit;
 
 
  //    // the smallest cell has been found
  //      //cout<<"width: "<<peak_width<<"  "<<"high: "<<peak_high<<endl;
  //      int column=(int)(bin_i*10+bin_j);
  //      peakArea[column]=peak_width*peak_high;
  //      m_peakWidth[column]=peak_width;
  //      m_peakHigh[column]=peak_high;
  ///     m_peakArea[column]=peak_width*peak_high;
  //      int columnx_split=(int)(column/10)*100+100;
  //      int columny_split=(int)(column%10)*100+100;
  //      //cout<<"column: "<<column<<"  "<<"("<<columnx_split<<","<<columny_split<<")"<<peakArea[column]<<endl;
  delete h1;
  //    }
  //  }
  //  double area_least=peakArea[0];
  //  for(int i=0;i<100;i++){
  //    if(area_least>peakArea[i]){
  //      area_least=peakArea[i];
  //    }
  //  }
  //  int num_area_least=0;
  //  for(int i=0;i<100;i++){
  //    if(area_least==peakArea[i]){
  //      num_area_least=i;
  //    }
  //  }
  //  int binx_split=(int)(num_area_least/10)*100+100;
  //  int biny_split=(int)(num_area_least%10)*100+100;
  //  m_areaLeast=area_least;
  //  m_areaLeastNum=num_area_least;
  //cout<<"the max peakArea: "<<num_area_least<<"("<<binx_split<<","<<biny_split<<")"<<"peakarea is : "<<area_least<<endl;
   */
 
 
 
 
  // select hit
  int count_use=0;
  double sum=0;
  double sum2=0;    
  for(int i=0;i<binX;i++){
    for(int j=0;j<binY;j++){
      int count=h1->GetBinContent(i+1,j+1);
      //        cout<<"("<<i<<","<<j<<")"<<"  "<<"count:  "<<count<<"          ";
      if(count!=0){
        count_use++;
        sum+=count;
        sum2+=count*count;
      }
    }
  }
  double f_n=m_ndev;
  cout<<"count_use:  "<<count_use<<endl;
  double aver=sum/count_use;
  double dev=sqrt(sum2/count_use-aver*aver);
  int min_counts=(int)(aver+f_n*dev);
  cout<<"aver:  "<<aver<<"  "<<"dev:  "<<dev<<"  min: "<<min_counts<<endl;
 
 
  vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) ); 
  vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) ); 
  //  int hough_trans_CS[binX][binY];
  // int hough_trans_CS_peak[binX][binY];
  for(int i=0;i<binX;i++){
    for(int j=0;j<binY;j++){
      hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1);
    }
  }
 
  int peak_num=0;
  for (int r=0; r<binY; r++) {
    for (int a=0; a<binX; a++) {
      long max_around=0;            
      for (int ar=a-1; ar<=a+1; ar++) {
        for (int rx=r-1; rx<=r+1; rx++) {
          int ax;
          if (ar<0) { ax=ar+binX; }
          else if (ar>=binX) { ax=ar-binX; }
          else { ax=ar; }
          if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
            if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; }
          }
        }
      }
      if (hough_trans_CS[a][r]<=min_counts) { hough_trans_CS_peak[a][r]=0; }
      else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; }
      else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; }
      else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; }
      if(hough_trans_CS_peak[a][r]!=0) { 
        //  cout<<" peak:  "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<"  "<<hough_trans_CS[a][r]<<endl;
        peak_num++;
      }
    }
  }
  //cout<<"peak_num:   "<<peak_num<<endl;
 
  //  //find double-point-peaks in parameter space
  int peak_num_2=0;
  for (int r=0; r<binY; r++) {
    for (int a=0; a<binX; a++) {
      if (hough_trans_CS_peak[a][r]==1) {
        hough_trans_CS_peak[a][r]=2;
        for (int ar=a-1; ar<=a+1; ar++) {
          for (int rx=r-1; rx<=r+1; rx++) {
            int ax;
            if (ar<0) { ax=ar+binX; }
            else if (ar>=binX) { ax=ar-binX; }
            else { ax=ar; }
            if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
              if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; }
            }
          }
        }
      }
      if(hough_trans_CS_peak[a][r]!=0) {
        peak_num_2++;
        //cout<<" peak:  "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<"  "<<hough_trans_CS[a][r]<<endl;
        //      for(int i=0;i<hough_trans_CS[a][r];i++){
        //        int hitNum=selectNum[a][r][i];
        //        cout<<"hit:  "<<hitNum<<"("<<layer[hitNum]<<","<<cell[hitNum]<<")"<<endl;
        //      }
      }
    }
  }
  //cout<<"peak_mum_2:   "<<peak_num_2<<endl;
 
  //fill the histogram again
  TH2D *h2 = new TH2D("h2","test2",binX,0,3.14,binY,-t_maxP,t_maxP);
  for(int i=0;i<binX;i++){
    for(int j=0;j<binY;j++){
      if(hough_trans_CS_peak[i][j]==2){
        h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]);
      }
    }
  }
 
  vector<int> peakList1[3];
  for(int n=0;n<peak_num_2;n++){
    for (int r=0; r<binY; r++) {
      for (int a=0; a<binX; a++) {
        if (hough_trans_CS_peak[a][r]==2) {
          peakList1[0].push_back(a+1);
          peakList1[1].push_back(r+1);
          peakList1[2].push_back(hough_trans_CS[a][r]);
        }
      }
    }
  }
  cout<<"finish peak?"<<endl;
  if(m_method==0) cout<<"numCross:  "<<numCross<<endl;
  //  //sort peaks by height
  int npeak1=peak_num_2;
  int n_max;
  for (int na=0; na<npeak1-1; na++) {
    n_max=na;
    for (int nb=na+1; nb<npeak1; nb++) {
      if (peakList1[2][n_max]<peakList1[2][nb]) { n_max=nb; }
    }   
    if (n_max!=na) {        // swap
      float swap[3];
      for (int i=0; i<3; i++) {
        swap[i]=peakList1[i][na];
        peakList1[i][na]=peakList1[i][n_max];
        peakList1[i][n_max]=swap[i];
      }
    }
  }
  cout<<"npeak1:  "<<npeak1<<endl;
  for(int n=0;n<npeak1;n++){
    int bintheta=peakList1[0][n];
    int binrho=peakList1[1][n];
    int count=peakList1[2][n];
    for(int i=0;i<count;i++){
      //cout<<"("<<bintheta<<","<<binrho<<","<<count<<")"<<endl;
      //cout<<"n"<<n<<":  "<<"("<<bintheta<<","<<binrho<<","<<count<<"):"<<vec_selectNum[bintheta-1][binrho-1][i]<<endl;
    }
    //cout<<"peakList1:  "<<n<<"  "<<"bina:  "<<peakList1[0][n]<<"  "<<"binr:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl; 
  }
 
  typedef std::map<int, int> M2;
  typedef std::map<int, M2> M1;
  M2 peak_combine_num;
  M1 hit_combine;
 
  int peak_track=0;
 
  if(m_method==2){
    //combine aroud
    int m=0;
    int n=0;
    int addnum=999;
    vector<bool> peaktrack(npeak1,true);
    while(addnum!=0)
    {
      addnum=0;
      double peak_cellNum[43];
      int peakX[43];
      int peakY[43];
      double  bin_left[43];
      double  bin_right[43];
      double  bin_up[43];
      double  bin_down[43];
      double area_left[43];
      double area_right[43];
      //double area_mid[43];
      double area_down[43];
      double area_up[43];
      // for(int i=0;i<npeak1;i++)  { peaktrack[i]=true; }
      for(int iLayer=0; iLayer<m_gm->nLayer(); iLayer++){
        peak_cellNum[iLayer]=m_peakCellNum[iLayer];
        peakX[iLayer]=peakList1[0][n];        //start from the highest peak
        peakY[iLayer]=peakList1[1][n];
        bin_left[iLayer]=peakX[iLayer]-peak_cellNum[iLayer];
        bin_right[iLayer]=peakX[iLayer]+peak_cellNum[iLayer];
        bin_up[iLayer]=peakY[iLayer]+peak_cellNum[iLayer];
        bin_down[iLayer]=peakY[iLayer]-peak_cellNum[iLayer];
        area_left[iLayer]=(bin_left[iLayer]-1)*binwidth;
        area_right[iLayer]=(bin_right[iLayer])*binwidth;
        // area_mid=(area_left+area_right)/2.;
        area_down[iLayer]=(bin_down[iLayer]-1)*binhigh-t_maxP;
        area_up[iLayer]=(bin_up[iLayer])*binhigh-t_maxP;
      }
      int count_peak=0;
      for(int k=0;k<nhit;k++){
        int layer=vec_layer[k];
        double lineLeft=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
        double lineRight=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
        // double lineMid=vec_u[k]*cos(area_mid)+vec_v[k]*sin(area_mid);
        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)){ 
          //      hit_combine[m].push_back(k);
          //        if(m==0){
          //            nHitSelect++;
          //          if(vec_track_index[k]>=0){nHitSignal_select++;}
          //          if(vec_track_index[k]>=0&&vec_slant[k]==0){nHitAxialSignal_select++;}
          //        }
          hit_combine[m][count_peak]=k;
          //vec_truthHit[m][k]=1;
          count_peak++;
        }
        //else {vec_truthHit[m][k]=0;}
        peak_combine_num[m]=count_peak;
        //int sizepeak=hit_combine[0].size();
        //cout<<"peak1size:  "<<sizepeak<<endl;
      }
      for(int i=n+1;i<npeak1;i++){
        if(peaktrack[i]==false) continue;
        int peaktheta=peakList1[0][i];
        int peakrho=peakList1[1][i];
        int peakhitNum=peakList1[2][i];
        int count_same_hit=0;
        for(int j=0;j<peakhitNum;j++){
          //for(int k=0;k<hit_combine[m].size();k++)
          for(int k=0;k<peak_combine_num[m];k++){
            if(vec_selectNum[peaktheta-1][peakrho-1][j]==hit_combine[m][k]){
              count_same_hit++;
            }
          }
        }
        double f_hit=m_hit_pro;
        if(count_same_hit>f_hit*peakhitNum){
          peaktrack[i]=false;
        }
      }
      for(int i=n+1;i<npeak1;i++){
        if(peaktrack[i]==true){
          addnum=i;
          break;
        }
      }
      if(addnum!=0) m++; 
      cout<<"peak_m:  "<<m+1<<endl;
      n=n+addnum;
    }
 
 
    for(int i=0;i<npeak1;i++){
      cout<<"track"<<i<<":  "<<"("<<peakList1[0][i]<<","<<peakList1[1][i]<<","<<peakList1[2][i]<<")"<<"  "<<"truth:  "<<peaktrack[i]<<endl;
      if(peaktrack[i]==true){
        peak_track++;
      }
    }
    for( int i=0;i<peak_track;i++){
      for(int j =0;j<peak_combine_num[i];j++){
        int hit_number=hit_combine[i][j];
        cout<<" peak "<<i<<"  has select hits: "<<vec_layer[hit_number]<<"   "<<vec_wire[hit_number]<<endl;
      }
      cout<<"has collect :"<<peak_combine_num[i]<<" hits "<<endl;
    }
    cout<<"event"<<t_eventNum<<":  "<<"has found: "<<peak_track<<" track."<<endl;
    m_trackNum=peak_track;
    //m_nHitSelect=nHitSelect;
    //m_nHitSignal_select=nHitSignal_select;
    //m_nHitAxialSignal_select=nHitAxialSignal_select;
 
    //according mc and track
    vector< vector<int> > rec_mc_num(peak_track,vector<int>(track_index_max,0) );
    vector< vector<int> > rec_mc_num_axial(peak_track,vector<int>(track_index_max,0) );
    if(track_index_max!=peak_track)  cout<<"not match track number !"<<endl;
    for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
      for(int i=0;i<peak_track;i++){
        for(int j=0;j<peak_combine_num[i];j++){
          for(int k=0;k<mc_track_hit[mc_track_num].size();k++){
            if(hit_combine[i][j]==mc_track_hit[mc_track_num][k]){
              rec_mc_num[i][mc_track_num]++; 
              int hit_num=mc_track_hit[mc_track_num][k];
              if(vec_slant[hit_num]==0) rec_mc_num_axial[i][mc_track_num]++; 
            }
          }
        }
      }
    }
    vector<int> rec_mc(peak_track,999);
    for(int i=0;i<peak_track;i++){
      for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
        if(rec_mc_num[i][mc_track_num]>0.5*peak_combine_num[i])  rec_mc[i]=mc_track_num; 
        cout<<"trackrec: "<<i<<" trackmc:  "<<mc_track_num<<"   "<<rec_mc_num[i][mc_track_num]<<"  "<<peak_combine_num[i]<<endl;
      }
    }
    for(int i=0;i<peak_track;i++){
      cout<<"rec :"<<i<<"belong to mc track: "<<rec_mc[i]<<endl;
    }
    for(int i=0;i<peak_track;i++){
      int mc_track_num=rec_mc[i];
      if(mc_track_num!=999) {
        cout<<"debug mc_track_num: "<<mc_track_num<<endl;
        m_nHitSignal_select[i]=rec_mc_num[i][mc_track_num];          //????????? mc_track_num=999?
        m_nHitAxialSignal_select[i]=rec_mc_num_axial[i][mc_track_num];          //????????? mc_track_num=999?
        cout<<"m_nHitSignal_select:  "<<m_nHitSignal_select[i]<<endl;
      }
    }
 
    //cout<<"peak1: "<<endl;
    //for(int i=0;i<peak_track;i++){
    //  for(int j=0;j<peak_combine_num[i];j++){
    //  int hitNumtemp=hit_combine[i][j];
    //  vec_truthHit[hitNumtemp]=i;
    //  cout<<"         "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<":  "<<vec_track_index[hitNumtemp]<<endl;
    //  }
    //}
 
    //  if(peak_track==2){
    //    cout<<"peak2: "<<endl;
    //    for(int i=0;i<hit_combine[1].size();i++){
    //      int hitNumtemp=hit_combine[1][i];
    //      vec_truthHit[hitNumtemp]=1;
    //      cout<<"         "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<":  "<<vec_track_index[hitNumtemp]<<endl;
    //    }
    //  }
 
 
  }
 
 
  //yzhang 2015-03-03 delete
//  //return hit number
//  if(m_method==0){
//  vector<int> numCross_select;
//  vector< vector<bool> > vec_hitSelect(npeak1,vector<bool> (nhit,false) );
//  vector<int> vec_selectHitNum;
//  vector<int> vec_selectHitNum_signal;
//  for(int peakNum=0;peakNum<npeak1;peakNum++){
//    int peak_cellNum=m_peakCellNum;
//    //    vector<bool> vec_hitSelect(nhit,false);
//    int peakX=peakList1[0][peakNum];
//    int peakY=peakList1[1][peakNum];
//    int bin_left=peakX-peak_cellNum;
//    int bin_right=peakX+peak_cellNum;
//    int bin_up=peakY+peak_cellNum;
//    int bin_down=peakY-peak_cellNum;
//    double area_left=(bin_left-1)*binwidth;
//    double area_right=(bin_right)*binwidth;
//    double area_down=(bin_down-1)*binhigh-t_maxP;
//    double area_up=(bin_up)*binhigh-t_maxP;
//    // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
//    double peak_width=(2*peak_cellNum+1)*binwidth;
//    double peak_high=(2*peak_cellNum+1)*binhigh;
//    //cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
//    //test if the hit is in this area
//    int numCross_select_temp=0;
//    for(int i=0;i<numCross;i++){
//      if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
//        numCross_select_temp++;
//        int hitNum_1=vec_hitNum[0][i];
//        int hitNum_2=vec_hitNum[1][i];
//        vec_hitSelect[peakNum][hitNum_1]=true;
//        vec_hitSelect[peakNum][hitNum_2]=true;
//      }
//    }
//    numCross_select.push_back(numCross_select_temp);
//    int selectHitNum=0;
//    int selectHitNum_signal=0;
//    for(int i=0;i<nhit;i++){
//      if(vec_hitSelect[peakNum][i]==true) {
//        selectHitNum++;
//        if(vec_hitSignal[i]==1){
//          selectHitNum_signal++;
//        }
//        //    cout<<"if hit is select: "<<i<<":  "<<vec_hitSelect[i]<<endl;
//      }
//    }
//    vec_selectHitNum.push_back(selectHitNum);
//    vec_selectHitNum_signal.push_back(selectHitNum_signal);
//    //  cout<<"hit: "<<selectHitNum<<"  "<<(double)selectHitNum/(double)nhit<<endl;
//    //cout<<"true hit : "<<selectHitNum_signal<<"  "<<(double)selectHitNum_signal/(double)nhit<<endl;
//    m_areaSelectHit=(double)selectHitNum/(double)nhit;
//    m_areaSelectHit_signal=(double)selectHitNum_signal/(double)nhit;
//  }
//
//  for(int i=0;i<npeak1;i++){
//    int selectHitNum=0;
//    for(int j=0;j<nhit;j++){
//      if(vec_hitSelect[i][j]==true){
//        selectHitNum++;
//        //  cout<<"event"<<t_eventNum<<": "<<"peak"<<i<<": "<<"hit: "<<"("<<vec_layer[j]<<","<<vec_wire[j]<<")"<<endl;    //cout picked hit
//      }
//    }
//    cout<<"peak"<<i<<": "<<selectHitNum<<"  hits"<<endl;
//  }
//  //  for(int i=0;i<3;i++){
//  //  int bina=peakList1[0][i];
//  //  int binr=peakList1[1][i];
//  //  int hit_peak=peakList1[2][i];
//  //  for(int j=0;j<hit_peak;j++){
//  //    int hitList=vec_selectHit[bina][binr][j];
//  //    cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
//  //    //      cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
//  //  }
//  //  }
//  //cout<<npeak1<<endl;
//
//  for(int n=0;n<npeak1;n++){
//    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;
//  }
//  }
 
  /*
  //combine the neighbor cell
  //vector<int> peak_hitList[npeak1];
  vector< vector <int> > peak_hitList(npeak1,vector <int> () ); 
  for(int i=0;i<npeak1;i++){
  int bina=peakList1[0][i];
  int binr=peakList1[1][i];
  int hit_peak=peakList1[2][i];
  for(int j=0;j<hit_peak;j++){
  peak_hitList[i].push_back(vec_selectHit[bina][binr][j]);
//  cout<<"j:  "<<j<<"  "<<vec_selectHit[bina][binr][j]<<endl;
}
 
for(int a=bina-1;a<=bina+1;a++){
for(int rx=binr-1;rx<=binr+1;rx++){
int ax;
if (a<0) { ax=a+binX; }
else if (a>=binX) { ax=a-binX; }
else { ax=a; }
if ( (ax!=bina || rx!=binr) && rx>=0 && rx<binY) {
int hit_around=count[ax][rx]; 
for(int k=0;k<hit_around;k++){
bool combine=1;
for(int m=0;m<hit_peak;m++){
if(vec_selectHit[ax][rx][k]==peak_hitList[i][m]){
//              cout<<i<<"  "<<"  "<<ax<<"  "<<rx<<"  "<<k<<"  "<<vec_selectHit[ax][rx][k]<<endl;
combine=0;
continue;
}
}
if(combine==1)  {
peak_hitList[i].push_back(vec_selectHit[ax][rx][k]);
peakList1[2][i]++;
hit_peak++;
//        cout<<" there is a combine"<<"  "<<vec_selectHit[ax][rx][k]<<endl;
}
}
}
}
}
}
 
 
//for(int n=0;n<npeak1;n++){
//  cout<<"peakList2:  "<<n<<"  "<<"a:  "<<peakList1[0][n]<<"  "<<"r:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl;
//}
 
 
for(int i=0;i<1;i++){
int bina=peakList1[0][i];
int binr=peakList1[1][i];
int hit_peak=peakList1[2][i];
for(int j=0;j<hit_peak;j++){
int hitList=peak_hitList[i][j];
//    cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
//    cout<<"peaknum:  "<<i<<"("<<bina<<","<<binr<<")"<<"  "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
}
}
 
//sort another time 
int n_max_2;
for (int na=0; na<npeak1-1; na++) {
n_max_2=na;
for (int nb=na+1; nb<npeak1; nb++) {
if (peakList1[2][n_max_2]<peakList1[2][nb]) { n_max_2=nb; }
}   
if (n_max_2!=na) {        // swap
float swap[3];
for (int i=0; i<3; i++) {
swap[i]=peakList1[i][na];
peakList1[i][na]=peakList1[i][n_max_2];
peakList1[i][n_max_2]=swap[i];
}
}
}
*/
 
//  for(int n=0;n<npeak1;n++){
//  cout<<"peakList3:  "<<n<<"  "<<"a:  "<<peakList1[0][n]<<"  "<<"r:  "<<peakList1[1][n]<<"  "<<"houghpeak:  "<<peakList1[2][n]<<endl;
//  }
 
//  int max_count=-999;
//  for(int i=0;i<binX;i++){
//  for(int j=0;j<binY;j++){
//    if(max_count<count[i][j]) {max_count=count[i][j];} 
//  }
//  }
 
//  m_maxCount=max_count;
//if(m_debug>0) cout<<"maxcount:   "<<max_count<<endl;
 
//  for(int i=0;i<binX;i++){
//  for(int j=0;j<binY;j++){
//    //        cout<<"("<<i<<","<<j<<")"<<"  "<<"count:  "<<count[i][j]<<"  "<<endl;
//    for(int l=0,m=vec_selectNum[i][j].size();l<m;++l){
//      if(count[i][j]==max_count){
//        //int layerId=vec_layer.at(vec_selectNum[i][j].at(l));
//        //if ((layerId>=8&&layerId<=19)||(layerId>=36)) { 
//        //    nHitAxialSelect++;
//        //}
//        vec_truthHit[vec_selectNum[i][j][l]]=true;  
//
//      }
//    }
//  }
//}
//  for(int i=0;i<n;i++){
//  cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<endl;
//  }
//  cout<<"selected true hits is:    "<<endl;
 
//for(int i=0;i<n;i++){
//  if(vec_truthHit[i]==true){
//    cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<"  "<<"hitSignal:"<<vec_hitSignal[i]<<endl;
//  }
//}
 
 
//int nHitAxialSelect=0;
//int nHitSelect=0;
//int nHitSignal_select=0;
//int nHitAxialSignal_select=0;
//for(int i=0;i<n;i++){
//  if(vec_truthHit[i]==false) continue;
//  nHitSelect++;
//  if(vec_hitSignal[i]==1) {nHitSignal_select++;}
//  int layerId=vec_layer[i];
//  if ((layerId>=8&&layerId<=19)||(layerId>=36)) { 
//  nHitAxialSelect++;
//  if(vec_hitSignal[i]==1) {nHitAxialSignal_select++;}
//  }   
//}
//
//m_nHitAxialSelect=nHitAxialSelect;
//m_nHitSelect=nHitSelect;
//
//m_nHitSignal_select=nHitSignal_select;
//m_nHitAxialSignal_select=nHitAxialSignal_select;
 
//  backtransform all track-
 
double d0=-999.;
double phi0=-999.;
double omega=-999.;
double z0=0;
double tanl=0;
 
for(track_fit=0;track_fit<peak_track;track_fit++){
  for(int i=0;i<nhit;i++){
    vec_truthHit[i]=false;
  }
  cout<<"track: "<<track_fit<<"  has select: "<<peak_combine_num[track_fit]<<"  hit ."<<endl;
  for(int i=0;i<peak_combine_num[track_fit];i++){
    int hitNum=hit_combine[track_fit][i];
    vec_truthHit[hitNum]=true; 
  }
 
  int nHitAxialSelect_temp=10;
  int leastSquare=LeastSquare(nhit,nHitAxialSelect_temp,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0,phi0,omega);
 
  if(leastSquare==0){
 
    TrkExchangePar tt(d0,phi0,omega,z0,tanl);
    TrkCircleMaker circlefactory;
    float chisum =0.;
    TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
    bool permitFlips = true;
    bool lPickHits = m_pickHits;
    circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
    // combine hit list
    int nDigi = digiToHots(newTrack,vec_truthHit);
    int nRecTk = 0;
    //if(m_debug>0) newTrack->printAll(std::cout);
    //fit 
    TrkHitList* qhits = newTrack->hits();
    TrkErrCode err=qhits->fit();
    //cout<<"++++++++++++++++++"<<err.failure()<<endl;
    const TrkFit* newFit = newTrack->fitResult();
    bool fitSuccess = false;
    //test fit result
    if (!newFit || (newFit->nActive()<3)) {
      if(m_debug>0){
        cout << "evt "<<t_eventNum<<" global fit failed ";
        if(newFit) cout <<" nAct "<<newFit->nActive();
        cout<<endl;
      }
      //return StatusCode::SUCCESS;
    }else{
      nRecTk = 1;
      fitSuccess = true;
      m_nEvtSuccess++;
      //    if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
      if(m_debug>0) newTrack->print(std::cout);
 
    }
 
    vector<int>  nfit2d(peak_track,0);
    if(m_debug>0) cout<<" hot list:"<<endl;
    TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
    while (hotIter!=qhits->hotList().end()) {
      nfit2d[track_fit]++;
      cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
        <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 
        <<":"<<hotIter->isActive()<<")  ";
      hotIter++;
    }
 
    m_2d_nFit[track_fit]=nfit2d[track_fit];
 
    //-------------try to out put parameters---
    double x_cirtemp;
    double y_cirtemp;
    double r_temp;
    if(err.failure()==0){
      //cout<<"++++++++++++++++++"<<err.failure()<<endl;
      TrkExchangePar par = newFit->helix(0.);
      d0=par.d0();
      phi0=par.phi0();
      omega=par.omega();
      //    cout<<"d0:                  "<<par.d0()<<"           "<<"d0temp:                     "<<-d0_temp<<endl;
      //    cout<<"phi0:                 "<<par.phi0()<<"         "<<"phi0temp:                  "<<phi0_temp<<endl;
      //    cout<<"w:               "<<par.omega()<<"         "<<"omegatemp:                  "<<omega_temp<<endl;
      // vector<double> b; // vector<double> x_stereo;
      r_temp=1./-par.omega();
      x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
      y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
 
      m_pt[track_fit]=r_temp/333.567;
      cout<<"pt_rec:   "<<m_pt[track_fit]<<endl;
    }
    else{
      m_nFitFailure[track_fit]=2;
    }
    //caculate z position
 
    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);
    //cout<<"z_stereoNum: "<<z_stereoNum<<endl;
    //cout<<"if zposition is finish: "<<endl;
 
    m_zStereoNum=z_stereoNum;
    for(int i=0;i<z_stereoNum;i++){
 
      m_zStereoNum=z_stereoNum;
      m_zStereo[i]=vec_zStereo[i];
      //      cout<<"eventNum:   "<<t_eventNum<<"   "<<"1111111111111111"<<endl;
      m_l[i]=l[i];
      if(m_debug>0) cout<<" l: "<<m_l[i]<<"  "<<"z:  "<<vec_zStereo[i]<<endl;
    }
    //  cout<<"333333333333333333"<<endl;
 
    Linefit(z_stereoNum,vec_zStereo,l,tanl,z0);
    
    cout<<"tanl:   "<<tanl<<endl;
    cout<<"z0:   "<<z0<<endl;
    z0=dz_mc;
    tanl=tanl_mc;
    cout<<"tanltruth:   "<<tanl<<endl;
    cout<<"z0truth:   "<<z0<<endl;
    //-------------------------------------------5 parameter fit-----------------------
 
    TrkExchangePar tt2(d0,phi0,omega,z0,tanl);
    TrkHelixMaker helixfactory;
    chisum =0.;
    TrkRecoTrk* newTrack2 = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0);
    permitFlips = true;
    lPickHits = m_pickHits;
    helixfactory.setFlipAndDrop(*newTrack2, permitFlips, lPickHits);
    // combine hit list
    int nDigi2 = digiToHots2(newTrack2,vec_truthHit);
    //  int nRecTk = 0;
    //if(m_debug>0) newTrack->printAll(std::cout);
    //fit 
    TrkHitList* qhits2 = newTrack2->hits();
    TrkErrCode err2=qhits2->fit();
    //cout<<"++++++++++++++++++"<<err.failure()<<endl;
    const TrkFit* newFit2 = newTrack2->fitResult();
    //  bool fitSuccess = false;
 
    //test fit result
    if (!newFit2 || (newFit2->nActive()<5)) {
      //fit failed
      if(m_debug>0){
        cout << "evt "<<t_eventNum<<" global fit failed ";
        if(newFit2) cout <<" nAct "<<newFit2->nActive();
        cout<<endl;
      }
      //return StatusCode::SUCCESS;
    }else{
      //fit success
      nFitSucess++;
      MdcTrack* mdcTrack = new MdcTrack(newTrack2);
      int tkStat = 4;//track find by Hough set stat=4
      int tkId = nTdsTk + track_fit;
      mdcTrack->storeTrack(tkId, trackList, hitList, tkStat);
      nRecTk = 1;
      fitSuccess = true;
      m_nEvtSuccess++;
      //    if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
      if(m_debug>0) newTrack2->print(std::cout);
 
    }
 
    vector<int>  nfit3d(peak_track,0);
    if(m_debug>0) cout<<" hot list:"<<endl;
    TrkHotList::hot_iterator hotIter2= qhits2->hotList().begin();
    while (hotIter2!=qhits2->hotList().end()) {
      nfit3d[track_fit]++;
      cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber()
        <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber() 
        <<":"<<hotIter2->isActive()<<")  ";
      hotIter2++;
    }
 
    m_3d_nFit[track_fit]=nfit3d[track_fit];
 
    // -------------try to out put parameters---
    if(err2.failure()==0){
      //cout<<"++++++++++++++++++"<<err2.failure()<<endl;
      TrkExchangePar par2 = newFit2->helix(0.);
 
      //cout<<"d0:                  "<<par2.d0()<<endl;
      //cout<<"phi0:                 "<<par2.phi0()<<endl;
      //cout<<"w:               "<<par2.omega()<<endl;
      //cout<<"z:               "<<par2.z0()<<endl;
      //cout<<"tanl:               "<<par2.tanDip()<<endl;
 
      r_temp=1./-par2.omega();
      x_cirtemp=(r_temp-par2.d0())*cos(par2.phi0()-M_PI/2.);
      y_cirtemp=(r_temp-par2.d0())*sin(par2.phi0()-M_PI/2.);
 
      m_d0[track_fit]=par2.d0();
      m_phi0[track_fit]=par2.phi0();
      m_omega[track_fit]=par2.omega();
      m_z0[track_fit]=par2.z0();
      m_tanl[track_fit]=par2.tanDip();
 
      double pxy=r_temp/333.567;
      double pz=pxy*par2.tanDip();
      double pxyz=sqrt(pxy*pxy+pz*pz);
      m_pt2[track_fit]=pxy;
      m_pz[track_fit]=pxy*par2.tanDip();
      m_pxyz[track_fit]=pxyz;
      cout<<"track"<<track_fit<<":  "<<"pt_rec2:   "<<m_pt2[track_fit]<<endl;
      cout<<"eventNum:  "<<t_eventNum<<"   "<<"track"<<track_fit<<":  "<<"pz_rec:   "<<m_pz[track_fit]<<endl;
    }
    else{
      m_nFitFailure[track_fit]=3;
    }
    // vector<double> b; // vector<double> x_stereo;
    //  r_temp=1./-par.omega();
    //  x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
    //  y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
    //  if(m_data==0) {
    //    m_mc_pt=r_temp/333.567;
    //    cout<<"pt_rec:   "<<m_mc_pt<<endl;
    //  }
    //  if(m_data==1) {m_pt=r_temp/333.567;}
    //  
    //  else{
    //  if(m_data==0) {m_mc_nfitfailure=2;}
    //  if(m_data==1){m_nFitFailure=2;}
    //  }
 
 
    //------------------------------------clean-------------------------
  }
  //vector<double>().swap(vec_u);
  //vector<double>().swap(vec_v);
  //vector<double>().swap(vec_p);
  //vector<double>().swap(vec_x);
  //vector<double>().swap(vec_y);
  //vector<double>().swap(vec_z);
  //vector<double>().swap(vec_x_east);
  //vector<double>().swap(vec_x_west);
  //vector<double>().swap(vec_y_east);
  //vector<double>().swap(vec_y_west);
  //vector<double>().swap(vec_z_east);
  //vector<double>().swap(vec_z_west);
  //vector<double>().swap(vec_zStereo);
  //vector<double>().swap(l);
  //vector<int>().swap(vec_layer);
  //vector<int>().swap(vec_wire);
}
m_nFitSucess=nFitSucess;
cout<<"  Hough finder find "<<nFitSucess<<" tot "<<nTdsTk<< endl;
 
delete h1;
delete h2;
cout<<"break:  ????"<<endl;
ntuplehit->write();
cout<<"finish event  "<<t_eventNum<<endl;
cout<<endl;
cout<<endl;
return StatusCode::SUCCESS;
}

execute() [2/2]

StatusCode MdcHoughFinder::execute

(

)

finalize() [1/2]

StatusCode MdcHoughFinder::finalize

(

)

Definition at line 1601 of file MdcHoughFinder.cxx.

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

finalize() [2/2]

StatusCode MdcHoughFinder::finalize

(

)

initialize() [1/2]

StatusCode MdcHoughFinder::initialize

(

)

Definition at line 96 of file MdcHoughFinder.cxx.

                                     {
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in initialize()" << endreq;
 
  StatusCode sc;
 
  //for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
 
  //nfailure=0;
  IPartPropSvc* p_PartPropSvc;
  static const bool CREATEIFNOTTHERE(true);
  sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
  if (!sc.isSuccess() || 0 == p_PartPropSvc) {
    log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
    return sc;
  }
  m_particleTable = p_PartPropSvc->PDT();
 
  // RawData
  IRawDataProviderSvc* irawDataProviderSvc;
  sc = service ("RawDataProviderSvc", irawDataProviderSvc);
  m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
  if ( sc.isFailure() ){
    log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
    return StatusCode::FAILURE;
 
  }
 
  //  Geometry
  IMdcGeomSvc*   imdcGeomSvc;
  sc = service ("MdcGeomSvc", imdcGeomSvc);
  m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
  if ( sc.isFailure() ){
    log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
    return StatusCode::FAILURE;
  }
 
  sc = service ("MagneticFieldSvc",m_pIMF);
  if(sc!=StatusCode::SUCCESS) {
    log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
  }
  m_bfield = new BField(m_pIMF);
  log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
 
  m_context = new TrkContextEv(m_bfield);
 
  Pdt::readMCppTable(m_pdtFile);
 
  //Get MdcCalibFunSvc
  IMdcCalibFunSvc* imdcCalibSvc; 
  sc = service ("MdcCalibFunSvc", imdcCalibSvc);
  m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
  if ( sc.isFailure() ){
    log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
    return StatusCode::FAILURE;
  }
 
  //initialize ntuplehit
  NTuplePtr nt(ntupleSvc(), "MdcHoughFinder/hit");
  if ( nt ){
    ntuplehit = nt;
  } else {
    ntuplehit = ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
    if(ntuplehit){
      ntuplehit->addItem("eventNum",           m_eventNum);
      ntuplehit->addItem("runNum",             m_runNum);
      ntuplehit->addItem("costaCut",        m_cosCut);
      //    test the first layer in hough space
      //      ntuplehit->addItem("Layer1Num",               m_nLayerNum,0, 43);
      //      ntuplehit->addItem("layer1X",    m_nLayerNum,   m_layer1X);
      //      ntuplehit->addItem("layer1Y",    m_nLayerNum,   m_layer1Y);
      //      ntuplehit->addItem("layer1U",    m_nLayerNum,   m_layer1U);
      //      ntuplehit->addItem("layer1V",    m_nLayerNum,   m_layer1V);
 
 
      ntuplehit->addItem("nHitMc",              m_nHit_Mc,0, 6796);
      ntuplehit->addItem("layerMc",    m_nHit_Mc,   m_layer_Mc);
      ntuplehit->addItem("cellMc",     m_nHit_Mc,   m_cell_Mc);
 
 
      ntuplehit->addItem("nHit",               m_nHit,0, 6796);
      ntuplehit->addItem("layerNhit",   43,    m_layerNhit);
      //      ntuplehit->addItem("hitCol",   m_nHit,   m_hitCol);
      ntuplehit->addItem("nCross",           m_nCross,0, 125000);
      ntuplehit->addItem("rho",    m_nCross,   m_rho);
      ntuplehit->addItem("theta",    m_nCross,   m_theta);
 
      ntuplehit->addItem("hitSignal",    m_nHit,   m_hitSignal);
      ntuplehit->addItem("layer",    m_nHit,   m_layer);
      ntuplehit->addItem("cell",     m_nHit,   m_cell);
      ntuplehit->addItem("x_east",   m_nHit,   m_x_east);
      ntuplehit->addItem("y_east",   m_nHit,   m_y_east);
      ntuplehit->addItem("z_east",   m_nHit,   m_z_east);
      ntuplehit->addItem("x_west",   m_nHit,   m_x_west);
      ntuplehit->addItem("y_west",   m_nHit,   m_y_west);
      ntuplehit->addItem("z_west",   m_nHit,   m_z_west);
      ntuplehit->addItem("x",        m_nHit,   m_x);
      ntuplehit->addItem("y",        m_nHit,   m_y);
      ntuplehit->addItem("z",        m_nHit,   m_z);
      ntuplehit->addItem("u",        m_nHit,   m_u);
      ntuplehit->addItem("v",        m_nHit,   m_v);
      ntuplehit->addItem("p",        m_nHit,   m_p);
      ntuplehit->addItem("slant",    m_nHit,   m_slant);
      //        ntuplehit->addItem("z_stereo",    m_stereohit,   m_zStereo);
      //        ntuplehit->addItem("z_num",    m_stereohit,   m_z_num);
      ntuplehit->addItem("maxCount",          m_maxCount);
 
      ntuplehit->addItem("peakColumn",           m_npeak,0,100);
      ntuplehit->addItem("peakWidth",  m_npeak,       m_peakWidth);
      ntuplehit->addItem("peakHigh",   m_npeak,       m_peakHigh);
      ntuplehit->addItem("peakArea",   m_npeak,       m_peakArea);
      ntuplehit->addItem("peakAreaLeast",       m_areaLeast);
      ntuplehit->addItem("peakAreaLeastNum",       m_areaLeastNum);
      ntuplehit->addItem("peakHitSelect",       m_areaSelectHit);
      ntuplehit->addItem("peakHitSelectSignal",       m_areaSelectHit_signal);
      //        ntuplehit->addItem("mc_x",    m_Mcnhit,  mc_x);
      //        ntuplehit->addItem("mc_y",    m_Mcnhit,  mc_y);
      //        ntuplehit->addItem("mc_z",    m_Mcnhit,  mc_z);
      ntuplehit->addItem("circleCenterX",      m_x_circle);
      ntuplehit->addItem("circleCenterY",      m_y_circle);
      ntuplehit->addItem("circleR",            m_r);
 
 
      ntuplehit->addItem("zStereoNum",           m_zStereoNum,0,1000);
      ntuplehit->addItem("zStereo",     m_zStereoNum,      m_zStereo );
      ntuplehit->addItem("l",            m_zStereoNum,      m_l);
 
      ntuplehit->addItem("trackNum_Mc",     m_trackNum_Mc,  0,100);
      ntuplehit->addItem("trackNum",     m_trackNum,  0,100);
      ntuplehit->addItem("d0",           m_trackNum,      m_d0);
      ntuplehit->addItem("phi0",         m_trackNum,      m_phi0);
      ntuplehit->addItem("omega",        m_trackNum,      m_omega);
      ntuplehit->addItem("z0",           m_trackNum,      m_z0);
      ntuplehit->addItem("tanl",         m_trackNum,      m_tanl);
 
      ntuplehit->addItem("pt_rec",     m_trackNum,     m_pt);
      ntuplehit->addItem("pt2_rec",    m_trackNum,     m_pt2);
      ntuplehit->addItem("pz_rec",     m_trackNum,     m_pz);
      ntuplehit->addItem("p_rec",      m_trackNum,     m_pxyz);
 
      ntuplehit->addItem("nHitSignal",                  m_nHitSignal);
      ntuplehit->addItem("nHitAxialSignal",             m_nHitAxialSignal);
      ntuplehit->addItem("nHitSignal_select",      m_trackNum,     m_nHitSignal_select);
      ntuplehit->addItem("nHitAxialSignal_selcet", m_trackNum,     m_nHitAxialSignal_select);
      ntuplehit->addItem("fitFailure",    m_trackNum,      m_nFitFailure);
      ntuplehit->addItem("nFit",             m_trackNum,   m_2d_nFit);
      ntuplehit->addItem("3dnFit",           m_trackNum,   m_3d_nFit);
      ntuplehit->addItem("nHitAxial",          m_nHitAxial);
      ntuplehit->addItem("nFitSucess",          m_nFitSucess);
      //      ntuplehit->addItem("nHitAxialSelect",   m_nHitAxialSelect);
      //      ntuplehit->addItem("nHitSelect",        m_trackNum,  m_nHitSelect);
 
      ntuplehit->addItem("pidTruth",   m_pidTruth);
      ntuplehit->addItem("costaTruth", m_costaTruth);
      ntuplehit->addItem("phiTruth",   m_phi0Truth);
      ntuplehit->addItem("drTruth",    m_drTruth);
      ntuplehit->addItem("dzTruth",    m_dzTruth);
      ntuplehit->addItem("ptTruth",    m_ptTruth);
      ntuplehit->addItem("pzTruth",    m_pzTruth);
      ntuplehit->addItem("pTruth",     m_pTruth);
      ntuplehit->addItem("qTruth",     m_qTruth);
 
 
    } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hit" << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
 
  if(m_debug>2)TrkHelixFitter::m_debug = true;
 
  return StatusCode::SUCCESS;
 
 
}

initialize() [2/2]

StatusCode MdcHoughFinder::initialize

(

)

---

The documentation for this class was generated from the following files:

• InstallArea/include/MdcHoughFinder/MdcHoughFinder/MdcHoughFinder.h
• Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-02/MdcHoughFinder/MdcHoughFinder.h
• MdcHoughFinder.cxx