MdcHoughFinder Class Reference

#include <Hough.h>

Inheritance diagram for MdcHoughFinder:

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

Detailed Description

Definition at line 45 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough.h.

Constructor & Destructor Documentation

MdcHoughFinder() [1/2]

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

Definition at line 76 of file Hough.cxx.

                                                                              :
  Algorithm(name, pSvcLocator)
{
  // Declare the properties  
  declareProperty("debug",          m_debug = 0);
  declareProperty("debugMap",       m_debugMap = 0);
  declareProperty("debug2D",          m_debug2D = 0);
  declareProperty("debugTrack",          m_debugTrack = 0);
  declareProperty("debugStereo",          m_debugStereo= 0);
  declareProperty("debugZs",          m_debugZs= 0);
  declareProperty("debug3D",          m_debug3D= 0);
  declareProperty("debugArbHit",          m_debugArbHit= 0);
  declareProperty("hist",           m_hist  = 1);
  declareProperty("filter",           m_filter= 0);
  //read raw data setup
  declareProperty("keepBadTdc",     m_keepBadTdc = 0);
  declareProperty("dropHot",        m_dropHot= 0);
  declareProperty("keepUnmatch",    m_keepUnmatch = 0);
  // combine pattsf
  declareProperty("combineTracking",m_combineTracking = false);
  declareProperty("removeBadTrack", m_removeBadTrack = 0);
  declareProperty("dropTrkDrCut",   m_dropTrkDrCut= 1.);
  declareProperty("dropTrkDzCut",   m_dropTrkDzCut= 10.);
  declareProperty("dropTrkPtCut",   m_dropTrkPtCut= 0.12);
  declareProperty("dropTrkChi2Cut", m_dropTrkChi2Cut = 10000.);
  //input setup
  declareProperty("inputType",      m_inputType = 0);
  //set MdcHoughFinder map 
  declareProperty("mapCharge",      m_mapCharge= -1);    //0 use whole ; 1 only half 
  declareProperty("useHalf",        m_useHalf= 0);    //0 use whole ; 1 only half 
  declareProperty("mapHitStyle",    m_mapHit= 0);    //0 : all ; 1 :axial
  declareProperty("nbinTheta",      m_nBinTheta = 100);
  declareProperty("nbinRho",        m_nBinRho = 100);
  declareProperty("rhoRange",       m_rhoRange = 0.1);
  declareProperty("peakWidth",      m_peakWidth= 3);
  declareProperty("peakHigh",       m_peakHigh= 1);
  declareProperty("hitPro",         m_hitPro= 0.4);
 
  declareProperty("recpos",         m_recpos= 1);
  declareProperty("recneg",         m_recneg= 1);
  declareProperty("combineSelf",    m_combine= 1);
  declareProperty("z0CutCompareHough",    m_z0Cut_compareHough= 10);
 
  //split drift circle
  declareProperty("n1",             m_npart= 100);
  declareProperty("n2",             m_n2= 16);
 
  declareProperty("pdtFile",        m_pdtFile = "pdt.table");
  declareProperty("eventFile",      m_evtFile= "EventList");
 
}

MdcHoughFinder() [2/2]

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

Member Function Documentation

beginRun() [1/2]

StatusCode MdcHoughFinder::beginRun

(

)

Definition at line 128 of file Hough.cxx.

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

beginRun() [2/2]

StatusCode MdcHoughFinder::beginRun

(

)

bookTuple() [1/2]

StatusCode MdcHoughFinder::bookTuple

(

)

Definition at line 1343 of file Hough.cxx.

                                    {
  MsgStream log(msgSvc(), name());
  NTuplePtr nt1(ntupleSvc(), "MdcHoughFinder/evt");
  if ( nt1 ){
    ntuple_evt = nt1;
  } else {
    ntuple_evt = ntupleSvc()->book("MdcHoughFinder/evt", CLID_ColumnWiseTuple, "evt");
    if(ntuple_evt){
      ntuple_evt->addItem("eventNum",           m_eventNum);
      ntuple_evt->addItem("runNum",             m_runNum);
 
      ntuple_evt->addItem("nHit",          m_nHit,0, 6796);
 
      //ntuple_evt->addItem("layer",         m_nHit,   m_layer);
      //ntuple_evt->addItem("cell",          m_nHit,   m_cell);
      //ntuple_evt->addItem("x",             m_nHit,   m_x);
      //ntuple_evt->addItem("y",             m_nHit,   m_y);
      //ntuple_evt->addItem("u",             m_nHit,   m_u);
      //ntuple_evt->addItem("v",             m_nHit,   m_v);
      //ntuple_evt->addItem("drift",         m_nHit,   m_drift);
      //ntuple_evt->addItem("r",             m_nHit,   m_r);
 
      //ntuple_evt->addItem("uTruth",        m_nHit,   m_uTruth);
      //ntuple_evt->addItem("vTruth",        m_nHit,   m_vTruth);
      //ntuple_evt->addItem("xTruth",        m_nHit,   m_xTruth);
      //ntuple_evt->addItem("yTruth",        m_nHit,   m_yTruth);
      //ntuple_evt->addItem("driftTruth",    m_nHit,   m_driftTruth);
      //ntuple_evt->addItem("rTruth",        m_nHit,   m_rTruth);
      //ntuple_evt->addItem("type",          m_nHit,   m_type);
 
      //ntuple_evt->addItem("deltaD",        m_nHit,   m_deltaD);
      //ntuple_evt->addItem("flt",           m_nHit,   m_flt);
      //ntuple_evt->addItem("slant",         m_nHit,   m_slant);
 
      //ntuple_evt->addItem("nSigAxial",     m_nSig_Axial, 0, 6796);
      //ntuple_evt->addItem("xybinNum",      m_xybinNum,   0,10000000); 
      //ntuple_evt->addItem("xybinMax",         m_xybinMax); 
      //ntuple_evt->addItem("xybinNL",       m_xybinNum,   m_xybinNL); 
      //ntuple_evt->addItem("xybinRL",       m_xybinNum,   m_xybinRL); 
      //ntuple_evt->addItem("xybinS",        m_xybinNum,   m_xybinS); 
      //ntuple_evt->addItem("theta_left",    m_theta_left); 
      //ntuple_evt->addItem("theta_right",   m_theta_right); 
      //ntuple_evt->addItem("rho_down",      m_rho_down); 
      //ntuple_evt->addItem("rho_up",        m_rho_up); 
      ////by houghmap maxbin value  0~M_PI
      //ntuple_evt->addItem("theta",         m_theta); 
      //ntuple_evt->addItem("rho",           m_rho); 
      //ntuple_evt->addItem("height",        m_height); 
      //ntuple_evt->addItem("aver",          m_aver); 
      //ntuple_evt->addItem("sigma",         m_sigma); 
 
      ////multiturn
      //ntuple_evt->addItem("multiturn",         m_exit_multiturn); 
 
      ////diff two charge
      //ntuple_evt->addItem("nMap1pk",          m_nMap1Pk,  0,100); 
      //ntuple_evt->addItem("nMap1tk",          m_nMap1Tk,  0,100); 
      //ntuple_evt->addItem("nMap2pk",          m_nMap2Pk,  0,100); 
      //ntuple_evt->addItem("nMap2tk",          m_nMap2Tk,  0,100); 
      //ntuple_evt->addItem("map1_pkrho",       m_nMap1Pk,  m_PkRho1); 
      //ntuple_evt->addItem("map1_pktheta",     m_nMap1Pk,  m_PkTheta1); 
      //ntuple_evt->addItem("map1_tkrho",       m_nMap1Tk,  m_TkRho1); 
      //ntuple_evt->addItem("map1_tktheta",     m_nMap1Tk,  m_TkTheta1); 
      //ntuple_evt->addItem("map2_pkrho",       m_nMap2Pk,  m_PkRho2); 
      //ntuple_evt->addItem("map2_pktheta",     m_nMap2Pk,  m_PkTheta2); 
      //ntuple_evt->addItem("map2_tkrho",       m_nMap2Tk,  m_TkRho2); 
      //ntuple_evt->addItem("map2_tktheta",     m_nMap2Tk,  m_TkTheta2); 
 
      ////MC truth
      //ntuple_evt->addItem("thetaline",     m_theta_line); 
      //ntuple_evt->addItem("rholine",       m_rho_line); 
 
      //      ntuple_evt->addItem("nmax",          m_nMaxLayerSlant,0,200); 
      //      ntuple_evt->addItem("maxslant",      m_nMaxLayerSlant,   m_MaxLayerSlant); 
      //      ntuple_evt->addItem("maxslantlayer", m_nMaxLayerSlant,   m_MaxLayer); 
      //      ntuple_evt->addItem("nnomax",        m_nNoMaxLayerSlant,0,1000); 
      //      ntuple_evt->addItem("nomaxslant",    m_nNoMaxLayerSlant, m_NoMaxLayerSlant); 
      //      ntuple_evt->addItem("nomaxlayerid",  m_nNoMaxLayerSlant, m_NoMaxLayerid); 
 
      //      ntuple_evt->addItem("MapMax",        m_MapMax ); 
      ntuple_evt->addItem("nMapPk",        m_nMapPk,  0,100); 
      //      ntuple_evt->addItem("nPeakNum",      m_nMapPk,  m_PeakOrder); 
      //      ntuple_evt->addItem("nPeakRho",      m_nMapPk,  m_PeakRho); 
      //      ntuple_evt->addItem("nPeakTheta",    m_nMapPk,  m_PeakTheta); 
      //      ntuple_evt->addItem("nPeakHeight",      m_nMapPk,  m_PeakHeight); 
      //      ntuple_evt->addItem("nPeakHit",      m_nMapPk,  m_PeakHit); 
      //      ntuple_evt->addItem("nPeakHitA",     m_nMapPk,  m_PeakHitA); 
      //      ntuple_evt->addItem("nMapTrk",       m_nMapTrk, 0,1000); 
      //      ntuple_evt->addItem("nTrkRho",       m_nMapTrk,  m_TrackRho); 
      //      ntuple_evt->addItem("nTrkTheta",     m_nMapTrk,  m_TrackTheta); 
      //      ntuple_evt->addItem("nTrkHitA",      m_nMapTrk,  m_TrackHitA); 
 
      //rec - charge
      //global 2d fit 
      //      ntuple_evt->addItem("nTrk2D_neg",        m_nTrk2D_neg,  0,100); 
      //      ntuple_evt->addItem("pt2D_neg",          m_nTrk2D_neg,  m_pt2D_neg); 
      //      ntuple_evt->addItem("nHit2D_neg",        m_nTrk2D_neg,  m_nHit2D_neg); 
      //      ntuple_evt->addItem("chi2_2D_neg",       m_nTrk2D_neg,  m_chi2_2D_neg); 
      //      //global 3d fit 
      //      ntuple_evt->addItem("nTrk3D_neg",        m_nTrk3D_neg,  0,100); 
      //      ntuple_evt->addItem("tanl_neg",          m_nTrk3D_neg,  m_tanl_neg); 
      //      ntuple_evt->addItem("tanl3D_neg",        m_nTrk3D_neg,  m_tanl3D_neg); 
      //      ntuple_evt->addItem("z0_neg",            m_nTrk3D_neg,  m_z0_neg); 
      //      ntuple_evt->addItem("z0_3D_neg",         m_nTrk3D_neg,  m_z0_3D_neg); 
      //      ntuple_evt->addItem("pro_neg",              m_nTrk3D_neg,  m_pro_neg); 
      //      ntuple_evt->addItem("pt3D_neg",          m_nTrk3D_neg,  m_pt3D_neg); 
      //      ntuple_evt->addItem("nHit3D_neg",        m_nTrk3D_neg,  m_nHit3D_neg); 
      //      ntuple_evt->addItem("chi2_3D_neg",       m_nTrk3D_neg,  m_chi2_3D_neg); 
      //
      //      //truth inf
      //      //      ntuple_evt->addItem("nTrkMC",     m_nTrkMC,0,10);
      //      ntuple_evt->addItem("pidTruth",   m_pidTruth);
      //      ntuple_evt->addItem("costaTruth", m_costaTruth);
      //      ntuple_evt->addItem("ptTruth",    m_ptTruth);
      //      ntuple_evt->addItem("pzTruth",    m_pzTruth);
      //      ntuple_evt->addItem("pTruth",     m_pTruth);
      //      ntuple_evt->addItem("qTruth",     m_qTruth);
      //      ntuple_evt->addItem("drTruth",    m_drTruth);
      //      ntuple_evt->addItem("phiTruth",   m_phi0Truth);
      //      ntuple_evt->addItem("omegaTruth", m_omegaTruth);
      //      ntuple_evt->addItem("dzTruth",    m_dzTruth);
      //      ntuple_evt->addItem("tanlTruth",  m_tanl_mc);
      //      ntuple_evt->addItem("rhoTruth",   m_rho_mc);
      //      ntuple_evt->addItem("thetaTruth", m_theta_mc);
 
 
      ntuple_evt->addItem("time", m_time);
 
    } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hough" <<endmsg;
      return StatusCode::FAILURE;
    }
  }
  /*
     NTuplePtr nt2(ntupleSvc(), "MdcHoughFinder/hit");
     if ( nt2 ){
     ntuplehit= nt2;
     } else {
     ntuplehit= ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
     if(ntuplehit){
     ntuplehit->addItem("evt_stereo",           m_evt_stereo);
     ntuplehit->addItem("run_stereo",           m_run_stereo);
     ntuplehit->addItem("cos_stereo",           m_cos_stereo);
     ntuplehit->addItem("tanlTruth_stereo",     m_tanlTruth_stereo);
     ntuplehit->addItem("ncir_stereo",          m_ncir_stereo);
     ntuplehit->addItem("npair_stereo",         m_npair_stereo);
     ntuplehit->addItem("tanl_stereo",          m_tanl_stereo);
     ntuplehit->addItem("tanl3D_stereo",        m_tanl3D_stereo);
     ntuplehit->addItem("z0_stereo",            m_z0_stereo);
     ntuplehit->addItem("z03D_stereo",          m_z03D_stereo);
     ntuplehit->addItem("nHit_axial",           m_nHit_axial, 0, 6796);
     ntuplehit->addItem("axial_layer",          m_nHit_axial, m_axial_layer);
     ntuplehit->addItem("axial_wire",           m_nHit_axial, m_axial_wire);
     ntuplehit->addItem("axial_deltaD",         m_nHit_axial, m_axial_deltaD);
     ntuplehit->addItem("axial_flt",            m_nHit_axial, m_axial_flt);
     ntuplehit->addItem("nHit_stereo",          m_nHit_stereo, 0, 6796);
     ntuplehit->addItem("stereo_layer",         m_nHit_stereo, m_stereo_layer);
     ntuplehit->addItem("stereo_wire",          m_nHit_stereo, m_stereo_wire);
     ntuplehit->addItem("stereo_style",         m_nHit_stereo, m_stereo_style);
     ntuplehit->addItem("stereo_z0",            m_nHit_stereo, m_stereo_z0);
     ntuplehit->addItem("stereo_z1",            m_nHit_stereo, m_stereo_z1);
     ntuplehit->addItem("stereo_s0",            m_nHit_stereo, m_stereo_s0);
     ntuplehit->addItem("stereo_s1",            m_nHit_stereo, m_stereo_s1);
     ntuplehit->addItem("stereo_z",             m_nHit_stereo, m_stereo_z);
     ntuplehit->addItem("stereo_s",             m_nHit_stereo, m_stereo_s);
     ntuplehit->addItem("stereo_zTruth",             m_nHit_stereo, m_stereo_zTruth);
     ntuplehit->addItem("stereo_sTruth",             m_nHit_stereo, m_stereo_sTruth);
     ntuplehit->addItem("stereo_deltaZ",             m_nHit_stereo, m_stereo_deltaZ);
     ntuplehit->addItem("stereo_nsol",          m_nHit_stereo, m_stereo_nsol);
     ntuplehit->addItem("stereo_ambig",         m_nHit_stereo, m_stereo_ambig);
     ntuplehit->addItem("stereo_ambig_truth",   m_nHit_stereo, m_stereo_ambig_truth);
     ntuplehit->addItem("stereo_disToCir",      m_nHit_stereo, m_stereo_disToCir);
     ntuplehit->addItem("stereo_cirlist",       m_nHit_stereo, m_stereo_cirlist);
     }
     }
     */
 
  return StatusCode::SUCCESS;
}

Referenced by initialize().

bookTuple() [2/2]

StatusCode MdcHoughFinder::bookTuple

(

)

execute() [1/2]

StatusCode MdcHoughFinder::execute

(

)

Definition at line 243 of file Hough.cxx.

                                   {
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
  cout.precision(6);
 
 
  m_timer_all->start();
 
  //event start time
  SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
  if (!evTimeCol) {
    log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
    return StatusCode::SUCCESS;
  }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
    }
  }
  HoughHit::setBunchTime(m_bunchT0);
 
  //-- 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;
  }
 
  t_eventNum=eventHeader->eventNumber();
  t_runNum=eventHeader->runNumber();
  if(m_hist){
    m_eventNum=eventHeader->eventNumber();
    m_runNum=eventHeader->runNumber();
  }
  if(m_debug>0) cout<<"begin evt "<<eventHeader->eventNumber()<<endl;
 
  //prepare tds 
  RecMdcTrackCol* trackList_tds;
  RecMdcHitCol* hitList_tds;
  vector<RecMdcTrack*> vec_trackPatTds;
  int nTdsTk = storeTracks(trackList_tds,hitList_tds,vec_trackPatTds);
  //print track in pattsf with bad vertex
  if(m_debug>0){
    RecMdcTrackCol::iterator iteritrk_pattsf = vec_trackPatTds.begin();
    for(;iteritrk_pattsf!=vec_trackPatTds.end();iteritrk_pattsf++){
      cout<<"in PATTSF LOST: "<<(*iteritrk_pattsf)->helix(0)<<" "<<(*iteritrk_pattsf)->helix(1)<<" "<<(*iteritrk_pattsf)->helix(2)<<" "<<(*iteritrk_pattsf)->helix(3)<<" "<<(*iteritrk_pattsf)->helix(4)<<" chi2 "<<(*iteritrk_pattsf)->chi2()<<endl;
    }
  }
 
  //for arbi hits
  MdcTrackParams  m_par;
  if(m_debugArbHit>0 ) m_par.lPrint=8;
  m_par.lRemoveInActive=1;
  //m_par.lUseQualCuts=0;
  //m_par.maxGapLength=99;
  //m_par.maxChisq=99;
  //m_par.minHits=99;
 
  // if filter read eventNum in file
  if(m_filter){
    ifstream lowPt_Evt;
    lowPt_Evt.open(m_evtFile.c_str());
    vector<int> evtlist;
    int evtNum;
    while( lowPt_Evt >> evtNum) {
      evtlist.push_back(evtNum);
    }
    vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),t_eventNum);
    if( iter_evt == evtlist.end() )  { setFilterPassed(false); return StatusCode::SUCCESS; }
    else setFilterPassed(true);
  }
 
  if(m_inputType == -1) GetMcInfo();
 
  //-- Get MDC digi vector
  if(m_debug>0) cout<<"step1 :  prepare digi "<<endl;
  MdcDigiVec mdcDigiVec = prepareDigi();
 
  //-- Create MdcHoughFinder hit list
  bool debugTruth = false;
  if(m_inputType == -1) debugTruth= true;
  if(m_debug>0) cout<<"step2 : hits-> hough hit list "<<endl;
  HoughHitList houghHitList(mdcDigiVec);
  //  if( mdcDigiVec.size() < 10 ) return StatusCode::SUCCESS;
  if( houghHitList.nHit() < 10 || houghHitList.nHit()>500 ) return StatusCode::SUCCESS;
  if(m_debug>0) houghHitList.printAll();
  if(debugTruth) houghHitList.addTruthInfo(g_tkParTruth);
 
  vector<MdcHit*>  mdcHitCol_neg;
  vector<MdcHit*>  mdcHitCol_pos;
  MdcDigiVec::iterator iter = mdcDigiVec.begin();
  for (;iter != mdcDigiVec.end(); iter++) {
    const MdcDigi* digi = (*iter);
    if( HoughHit(digi).driftTime()>1000 || HoughHit(digi).driftTime()<=0 ) continue;
    MdcHit *mdcHit= new MdcHit(digi, Global::m_gm);
    MdcHit *mdcHit_pos= new MdcHit(digi, Global::m_gm);
    mdcHitCol_neg.push_back(mdcHit);
    mdcHitCol_pos.push_back(mdcHit_pos);
  }
 
  HoughMap *m_map = new HoughMap(-1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
 
  //  if(m_hist) dumpHoughHitList(houghHitList);
  if(m_hist) m_nHit = houghHitList.nHit();
  if(m_hist) dumpHoughMap(*m_map);
 
  //track
  if(m_debug>0) cout<<"step3 : neg track list "<<endl;
  vector<HoughTrack*> vec_track_neg; 
  vector<HoughTrack*> vec_track2D_neg; 
  MdcTrackList  mdcTrackList_neg(m_par) ;
  if(m_recneg){
    HoughTrackList trackList_neg(*m_map);
    int trackList_size = trackList_neg.getTrackNum();
    vector< vector<int> > stat_2d(2,vector<int>(trackList_size,0) );
    vector< vector<int> > stat_3d(2,vector<int>(trackList_size,0) );
    int ifit=0;
    int ifit3D=0;
    for(int itrack=0;itrack<trackList_size;itrack++){
      if(m_debug>0) cout<<"begin track: "<<itrack<<endl;
      //suppose charge -
      if(m_debug>0) cout<<" suppose charge -1 "<<endl;
      HoughTrack &track_neg = trackList_neg.getTrack(itrack);
      track_neg.setMdcHit( &mdcHitCol_neg);
      track_neg.setHoughHitList(houghHitList.getList());
      track_neg.setCharge(-1);
      stat_2d[0][itrack] = track_neg.fit2D(m_bunchT0);
      int track_charge_2d = track_neg.trackCharge2D();
      if(m_debug>0) cout<<" charge -1 stat2d "<<stat_2d[0][itrack]<<" "<<track_charge_2d<<endl; 
      if( stat_2d[0][itrack] == 0 || track_charge_2d ==0  )  continue;
      //vec_track2D_neg.push_back( &track_neg );
      //fill 2D inf
      ifit++;
      //3D fit
      int nHit3d = track_neg.find_stereo_hit();
      int npair= track_neg.find_pair_hit();
      //cout<<"npair "<<npair<<endl;
 
      int track_charge_3d = track_neg.trackCharge3D();
      if(m_debug>0) cout<<" nhitstereo -1 "<<nHit3d<<" "<<track_charge_3d<<endl; 
      if( nHit3d <3 || track_charge_3d==0 )  continue;
 
      //choose fit method
      if( npair==0 ) stat_3d[0][itrack] = track_neg.fit3D();  
      else stat_3d[0][itrack] = track_neg.fit3D_inner();  
      //else stat_3d[0][itrack] = track_neg.fit3D();  
 
      if(m_debug>0) cout<<" charge -1 stat3d "<<stat_3d[0][itrack]<<" "<<endl;
      if( stat_3d[0][itrack]==0 ) continue;
      vec_track_neg.push_back( &track_neg );
      //fill 3D inf
      /*
         if(m_hist){
         m_tanl_neg[ifit3D] = track_neg.getTanl_zs(); 
         m_tanl3D_neg[ifit3D] = track_neg.getTanl(); 
         m_z0_neg[ifit3D] = track_neg.getZ0_zs(); 
         m_z0_3D_neg[ifit3D] = track_neg.getZ0(); 
         m_pro_neg[ifit3D] = track_neg.getPro(); 
         m_pt3D_neg[ifit3D] = track_neg.getPt3D(); 
         m_nHit3D_neg[ifit3D] = track_neg.getNfit3D(); 
         m_chi2_3D_neg[ifit3D] = track_neg.getChi2_3D(); 
         }
         ifit3D++;
         int axial_use=0;
         int stereo_use=0;
         int cir0=0;
         int cirmore=0;
         int ncir_track;
         for(int ister =0;ister<track_neg.getHoughHitList().size();ister++){
         HoughRecHit *rechit = &(track_neg.getHoughHitList().at(ister));
         if( rechit->getSlayerType()==0 ) {
         double deltaD  = rechit->getDeltaD();
         double flt= rechit->getFltLen();
         if(m_hist){
         m_axial_layer[axial_use]=rechit->getLayerId();
         m_axial_wire[axial_use]=rechit->getWireId();
         m_axial_deltaD[axial_use] = deltaD;
         m_axial_flt[axial_use] = flt;
         }
      //judge which circle
      int cirlist = rechit->getCirList();
      if(cirlist == 0)  cir0++;
      else cirmore++;
      axial_use++;
      }
      if( rechit->getSlayerType()!=0 ) {
      //stereo hit info
      double z0= rechit->getzAmb(0);
      double z1= rechit->getzAmb(1);
      double s0= rechit->getsAmb(0);
      double s1= rechit->getsAmb(1);
      int ambig = rechit->getAmbig();
      int ambigTruth = rechit->getLrTruth();
      int style = rechit->getStyle();
      //cout<<"ambig ambigTruth "<<ambig<<" "<<ambigTruth<<endl;
      double z = rechit->getzPos();
      double s = rechit->getsPos();
      double zTruth;
      double sTruth;
      if(ambigTruth == 1 ) {zTruth = z0;sTruth=s0;}
      else if(ambigTruth == -1 ) {zTruth = z1;sTruth=s1;}
      //calcu dist of hitz to zsline
      double deltaZ = zTruth - ( sTruth * track_neg.getTanl_zs()+track_neg.getZ0_zs());
      //else cout<<" not get ambig truth info "<<endl;
      int nsol = rechit->getnsol();
      double disToCir= rechit->getDisToCir();
      int multicir = rechit->getCirList();
      if(m_hist){
      m_stereo_layer[stereo_use]=rechit->getLayerId();
      m_stereo_wire[stereo_use]=rechit->getWireId();
      m_stereo_style[stereo_use]=style;
      m_stereo_z0[stereo_use]=z0;
      m_stereo_z1[stereo_use]=z1;
      m_stereo_s0[stereo_use]=s0;
      m_stereo_s1[stereo_use]=s1;
      m_stereo_z[stereo_use]=z;
      m_stereo_s[stereo_use]=s;
      m_stereo_zTruth[stereo_use]=zTruth;
      m_stereo_sTruth[stereo_use]=sTruth;
      m_stereo_deltaZ[stereo_use]=deltaZ;
      m_stereo_nsol[stereo_use]=nsol;
      m_stereo_ambig[stereo_use]=ambig;
      m_stereo_ambig_truth[stereo_use]=ambigTruth;
      m_stereo_disToCir[stereo_use]=disToCir;
      m_stereo_cirlist[stereo_use]=multicir;
    }
    stereo_use++;
    }
    }
    if( cir0>cirmore ) ncir_track = 0;
    else ncir_track = 1;
    if(m_hist){
      m_evt_stereo= t_eventNum;
      m_run_stereo= t_runNum;
      m_cos_stereo= t_cos;
      m_tanlTruth_stereo= t_tanl;
      m_ncir_stereo= ncir_track;
      m_npair_stereo= npair;
      m_tanl_stereo= track_neg.getTanl_zs();
      m_z0_stereo= track_neg.getZ0_zs();
      m_tanl3D_stereo= track_neg.getTanl();
      m_z03D_stereo= track_neg.getZ0();
      m_nHit_axial= axial_use;
      m_nHit_stereo = stereo_use;
      ntuplehit->write();
    }
    */
    }
 
    std::sort ( vec_track_neg.begin(),vec_track_neg.end(),more_pt); 
    //track for clear
    vector<MdcTrack*> vec_MdcTrack_neg;
    for( unsigned int i =0;i<vec_track_neg.size();i++){
      MdcTrack *mdcTrack = new MdcTrack(vec_track_neg[i]->getTrk());
      vec_MdcTrack_neg.push_back(mdcTrack);
      if(m_debug>0) cout<<"trackneg: "<<i<<" pt "<<vec_track_neg[i]->getPt3D()<<endl;
      if(m_debug>0) vec_track_neg[i]->print();
    }
    if(m_debug>0) cout<<"step4 : judge neg track list "<<endl;
    judgeHit(mdcTrackList_neg,vec_MdcTrack_neg);
    if(m_debug>0) cout<<"finish - charge "<<endl;
  }
 
  //do rec pos charge
 
  HoughMap *m_map2 = new HoughMap(1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
  if(m_debug>0) cout<<"step5 : pos track list "<<endl;
  vector<HoughTrack*> vec_track_pos; 
  MdcTrackList  mdcTrackList_pos(m_par) ;
  if(m_recpos){
    HoughTrackList tracklist_pos(*m_map2);
    int tracklist_size2 = tracklist_pos.getTrackNum();
    vector< vector<int> > stat_2d2(2,vector<int>(tracklist_size2,0) );
    vector< vector<int> > stat_3d2(2,vector<int>(tracklist_size2,0) );
    int ifit=0;
    int ifit3D=0;
    for(int itrack=0;itrack<tracklist_size2;itrack++){
      //suppose charge +
      if(m_debug>0) cout<<" suppose charge +1 "<<endl;
      HoughTrack &track_pos = tracklist_pos.getTrack(itrack);
      track_pos.setMdcHit( &mdcHitCol_pos);
      track_pos.setHoughHitList(houghHitList.getList());
      track_pos.setCharge(1);
      stat_2d2[0][itrack] = track_pos.fit2D(m_bunchT0);
      int track_charge_2d = track_pos.trackCharge2D();
      if(m_debug>0) cout<<" charge +1 stat2d "<<stat_2d2[1][itrack]<<" "<<track_charge_2d<<endl; 
      if( stat_2d2[0][itrack] == 0 || track_charge_2d==0 )  continue;
      //fill 2d inf
      ifit++;
      //3d fit
      int nHit3d = track_pos.find_stereo_hit();
      int npair= track_pos.find_pair_hit();
      //cout<<"npair "<<npair<<endl;
 
      int track_charge_3d = track_pos.trackCharge3D();
      if(m_debug>0) cout<<" nhitstereo +1 "<<nHit3d<<" "<<track_charge_3d<<endl; 
      if( nHit3d <3 || track_pos.trackCharge3D()==0 )  continue;
      //choose fit method
      if( npair==0 ) stat_3d2[0][itrack] = track_pos.fit3D();  
      else stat_3d2[0][itrack] = track_pos.fit3D_inner();  
      //else stat_3d2[0][itrack] = track_neg.fit3D();  
 
      if(m_debug>0) cout<<" charge +1 stat3d "<<stat_3d2[1][itrack]<<" "<<endl;
      if( stat_3d2[0][itrack]==0 ) continue;
      vec_track_pos.push_back( &track_pos );
      //fill 3d inf
      ifit3D++;
    }
 
    //sort pos tracklist
    std::sort ( vec_track_pos.begin(),vec_track_pos.end(),more_pt); 
 
    // clear pos track
    vector<MdcTrack*> vec_MdcTrack_pos;
    for( unsigned int i =0;i<vec_track_pos.size();i++){
      MdcTrack *mdcTrack = new MdcTrack(vec_track_pos[i]->getTrk());
      vec_MdcTrack_pos.push_back(mdcTrack);
      if(m_debug>0) cout<<"trackpos : "<<i<<" pt "<<vec_track_pos[i]->getPt3D()<<endl;
      if(m_debug>0) vec_track_pos[i]->print();
    }
    if(m_debug>0) cout<<"step6 : judge pos track list "<<endl;
    judgeHit(mdcTrackList_pos,vec_MdcTrack_pos);
  }
 
  //combine hough itself
  //  if(m_debug>0) cout<<"step7 : combine neg&pos track list "<<endl;
  if(m_combine){
    compareHough(mdcTrackList_neg);
    compareHough(mdcTrackList_pos);
  }
  if( mdcTrackList_neg.length()!=0 && mdcTrackList_pos.length()!=0 ) judgeChargeTrack(mdcTrackList_neg,mdcTrackList_pos);
  //add tracklist
  mdcTrackList_neg+=mdcTrackList_pos;   //neg -> all charge
 
  //compare pattsf&hough self
  if(m_combineTracking) nTdsTk = comparePatTsf(mdcTrackList_neg,trackList_tds);
 
  // store tds
  if(m_debug>0) cout<<"step8 : store tds  "<<endl;
  if(m_debug>0) cout<<"store tds "<<endl;
  int tkId = nTdsTk ;   //trkid
  for( unsigned int i =0;i<mdcTrackList_neg.length();i++){
    if(m_debug>0) cout<<"- charge size i : "<<i<<" "<<mdcTrackList_neg.length()<<endl;
    int tkStat = 4;//track find by houghspace set stat=4
    mdcTrackList_neg[i]->storeTrack(tkId, trackList_tds, hitList_tds, tkStat);
    tkId++;
    delete  mdcTrackList_neg[i];
  }
  if(m_debug>0)  m_mdcPrintSvc->printRecMdcTrackCol();
 
  m_timer_all->stop();
  double t_teventTime = m_timer_all->elapsed();
  if(m_hist) m_time = t_teventTime;
 
  if( m_hist )  ntuple_evt->write();
 
  delete m_map;
  delete m_map2;
  if(m_debug>0) cout<<"after delete map "<<endl;
  for(int ihit=0;ihit<mdcHitCol_neg.size();ihit++){
    delete mdcHitCol_neg[ihit];
    delete mdcHitCol_pos[ihit];
  }
 
  if(m_debug>0) cout<<"after delete hit"<<endl;
  //clearMem(mdcTrackList_neg,mdcTrackList_pos);
  if(m_debug>0) cout<<"end event "<<endl;
 
  return StatusCode::SUCCESS;
}

execute() [2/2]

StatusCode MdcHoughFinder::execute

(

)

finalize() [1/2]

StatusCode MdcHoughFinder::finalize

(

)

Definition at line 616 of file Hough.cxx.

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

finalize() [2/2]

StatusCode MdcHoughFinder::finalize

(

)

initialize() [1/2]

StatusCode MdcHoughFinder::initialize

(

)

Definition at line 137 of file Hough.cxx.

                                     {
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in initialize()" << endreq;
 
  StatusCode sc;
 
  //particle
  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;
  }
 
  //  magneticfield
  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);
 
  //read pdt
  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 MdcPrintSvc
  IMdcPrintSvc* imdcPrintSvc;
  sc = service ("MdcPrintSvc", imdcPrintSvc);
  m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
  if ( sc.isFailure() ){
    log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
    return StatusCode::FAILURE;
  }
 
  //time 
  sc = service( "BesTimerSvc", m_timersvc);
  if( sc.isFailure() ) {
    log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
    return StatusCode::FAILURE;
  }
  m_timer_all = m_timersvc->addItem("Execution");
  m_timer_all->propName("nExecution");
 
  //initialize static 
  HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
  HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
  HoughHit::setBunchTime(m_bunchT0);
  Hough2D::setContext(m_context);
  Hough2D::setCalib(m_mdcCalibFunSvc);
  Hough3D::setContext(m_context);
  Hough3D::setCalib(m_mdcCalibFunSvc);
 
  HoughHit::_npart=m_npart;
  HoughMap::m_useHalfCir=m_useHalf;
  HoughMap::m_N1=m_npart;
  HoughMap::m_N2=m_n2;
 
 
  if(m_debugMap> 0)   HoughMap::m_debug = 1;
  if(m_debug2D > 0)   Hough2D::m_debug = 1;
  if(m_debug3D > 0)   Hough3D::m_debug = 1;
  if(m_debugTrack> 0) HoughTrack::m_debug = 1;
  if(m_debugStereo> 0) HoughStereo::m_debug = 1;
  if(m_debugZs> 0)    HoughZsFit::m_debug = 1;
  if(m_debug3D > 4)   TrkHelixFitter::m_debug = 1;
 
  if(m_hist) sc = bookTuple();
  if ( sc.isFailure() ){
    log << MSG::FATAL << "Could not book Tuple !" << endreq;
    return StatusCode::FAILURE;
  }
  return StatusCode::SUCCESS;
}

initialize() [2/2]

StatusCode MdcHoughFinder::initialize

(

)

---

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

• InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough.h
• Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-14/MdcHoughFinder/Hough.h
• Hough.cxx