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 ()

Detailed Description

Definition at line 40 of file Hough.h.

Constructor & Destructor Documentation

MdcHoughFinder()

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

Definition at line 49 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" );
}

Referenced by MdcHoughFinder().

Member Function Documentation

beginRun()

StatusCode MdcHoughFinder::beginRun

(

)

Definition at line 99 of file Hough.cxx.

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

Referenced by execute().

bookTuple()

StatusCode MdcHoughFinder::bookTuple

(

)

Definition at line 1486 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().

execute()

StatusCode MdcHoughFinder::execute

(

)

Definition at line 219 of file Hough.cxx.

                                   {
  if ( !m_beginRun )
  {
    StatusCode sc = beginRun();
    if ( sc.isFailure() )
    {
      error() << "beginRun failed" << endmsg;
      return StatusCode::FAILURE;
    }
    m_beginRun = true;
  }
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
  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" << endmsg;
    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" << endmsg;
    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;
}

finalize()

StatusCode MdcHoughFinder::finalize

(

)

Definition at line 643 of file Hough.cxx.

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

initialize()

StatusCode MdcHoughFinder::initialize

(

)

Definition at line 108 of file Hough.cxx.

                                      {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  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" << endmsg;
    return sc;
  }
  m_particleTable = p_PartPropSvc->PDT();
 
  // RawData
  IRawDataProviderSvc* irawDataProviderSvc;
  sc                   = service( "RawDataProviderSvc", irawDataProviderSvc );
  m_rawDataProviderSvc = irawDataProviderSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << name() << " Could not load RawDataProviderSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  //  Geometry
  IMdcGeomSvc* imdcGeomSvc;
  sc           = service( "MdcGeomSvc", imdcGeomSvc );
  m_mdcGeomSvc = imdcGeomSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load MdcGeoSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  //  magneticfield
  sc = service( "MagneticFieldSvc", m_pIMF );
  if ( sc != StatusCode::SUCCESS )
  { log << MSG::ERROR << "Unable to open Magnetic field service" << endmsg; }
  m_bfield = new BField( m_pIMF );
  log << MSG::INFO << "field z = " << m_bfield->bFieldNominal() << endmsg;
  m_context = new TrkContextEv( m_bfield );
 
  // read pdt
  Pdt::readMCppTable( m_pdtFile );
 
  // Get MdcCalibFunSvc
  IMdcCalibFunSvc* imdcCalibSvc;
  sc               = service( "MdcCalibFunSvc", imdcCalibSvc );
  m_mdcCalibFunSvc = imdcCalibSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  // initialize MdcPrintSvc
  IMdcPrintSvc* imdcPrintSvc;
  sc            = service( "MdcPrintSvc", imdcPrintSvc );
  m_mdcPrintSvc = imdcPrintSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load MdcPrintSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  // time
  sc = service( "BesTimerSvc", m_timersvc );
  if ( sc.isFailure() )
  {
    log << MSG::WARNING << " Unable to locate BesTimerSvc" << endmsg;
    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 !" << endmsg;
    return StatusCode::FAILURE;
  }
  return StatusCode::SUCCESS;
}

---

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

• Hough.h
• Hough.cxx