MdcTrkRecon Class Reference

#include <MdcTrkRecon.h>

Inheritance diagram for MdcTrkRecon:

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

Protected Member Functions
bool	ignoringUsedHits () const
bool	poisoningHits () const
void	fillSegList ()
void	fillTrackList ()
void	dumpDigi ()
void	dumpTdsTrack (RecMdcTrackCol *trackList)
void	fillMcTruth ()
void	fillEvent ()
StatusCode	bookNTuple ()

Detailed Description

Definition at line 31 of file MdcTrkRecon.h.

Constructor & Destructor Documentation

MdcTrkRecon()

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

Definition at line 71 of file MdcTrkRecon.cxx.

                                                                        :
  Algorithm(name, pSvcLocator),
  m_hitData(0), m_segs(0), m_tracks(0), m_segFinder(0)  
{ 
  //Declare the properties--------------------------------
  declareProperty("FittingMethod",  m_fittingMethod = 2);
  declareProperty("ConfigFile",     m_configFile = "MDCConfig.xml");
  declareProperty("OnlyUnusedHits", m_onlyUnusedHits = 0);
  declareProperty("PoisonHits",     m_poisonHits = 0);
  declareProperty("doLineFit",      m_doLineFit = false);   
  declareProperty("paramFile",      m_paramFile = "params"); 
  declareProperty("pdtFile",        m_pdtFile = "pdt.table");
  declareProperty("allHit",         m_allHit = 1);
  declareProperty("hist",           m_hist= false);
  declareProperty("recForEsTime",   m_recForEsTime= false);
  declareProperty("tryBunch",       m_tryBunch = false);
  declareProperty("d0Cut",          m_d0Cut= -999.);
  declareProperty("z0Cut",          m_z0Cut= -999.);
  declareProperty("dropTrkPt",      m_dropTrkPt = -999.);
  declareProperty("debug",          m_debug= 0); 
  declareProperty("mcHist",         m_mcHist= false); 
  declareProperty("fieldCosmic",    m_fieldCosmic= false); 
  declareProperty("doSag",      m_doSagFlag= false);
  declareProperty("arbitrateHits",  m_arbitrateHits = true);
 
  declareProperty("helixHitsSigma", m_helixHitsSigma);
  //for fill hist 
  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);
  declareProperty("selEvtNo",       m_selEvtNo);
  declareProperty("combineTracking",m_combineTracking = false);//yzhang 2010-05-28 
 
#ifdef MDCPATREC_RESLAYER
  declareProperty("resLayer",m_resLayer= -1);
#endif
 
}

~MdcTrkRecon()

MdcTrkRecon::~MdcTrkRecon

(

)

Definition at line 113 of file MdcTrkRecon.cxx.

{ 
  m_segs.reset(0);
  m_segFinder.reset(0);
  m_hitData.reset(0);
  m_tracks.reset(0); 
}

Member Function Documentation

beginRun()

StatusCode MdcTrkRecon::beginRun

(

)

Definition at line 122 of file MdcTrkRecon.cxx.

                                {
  MsgStream log(msgSvc(), name());  
  log << MSG::INFO << "in beginRun()" << endreq;
 
  //Detector geometry 
  _gm = MdcDetector::instance(m_doSagFlag);
 
  if(NULL == _gm) return StatusCode::FAILURE;
 
  //Initialize segList
  m_segs.reset( new MdcSegList(_gm->nSuper(), m_flags.segPar) );
 
  return StatusCode::SUCCESS;
}

bookNTuple()

StatusCode MdcTrkRecon::bookNTuple ( )

protected

Definition at line 523 of file MdcTrkRecon.cxx.

                                  {
  MsgStream log(msgSvc(), name());
  StatusCode sc = StatusCode::SUCCESS;
 
 
  //--------------book histogram------------------  
  h_sfHit  = histoSvc()->book( "sfHit", "Hits after segment finding",46,-1.,44. );
  //g_actHit = histoSvc()->book( "actHit", "Active hits",46,-1.,44. );
  //g_hitEff = histoSvc()->book( "hitEff", "Hit efficiency in the event",100,-0.5,1.2 );
  g_residVsLayer = histoSvc()->book( "residVsLayer", "Residual vs. layer",60,-5,50.,100,-0.5,1.2 );
  g_residVsDoca  = histoSvc()->book( "residVsDoca", "Residial vs. doca",50,-2.,2.,100,-0.5,1.2 );
  //segment
  //g_segChi2  = histoSvc()->book( "segChi2", "chi2 per dof in segment finding",101,-1.,100. );
  g_cirTkChi2  = histoSvc()->book( "cirTkChi2", "chi2 per dof in circle finding",21,-1. ,20. );
  g_3dTkChi2 = histoSvc()->book( "chi2Helix", "maxChisq dof helix fit" ,51.,-1.,50);
  g_nSigAdd = histoSvc()->book( "nSigAdd", "max allowed n sigma to add hit to existing segment",50,0,50 );
  g_z0Cut = histoSvc()->book( "z0Cut", "z0 for including stereo segs in cluster",100,0,600 );
  g_ctCut = histoSvc()->book( "ctCut", "ct for including stereo segs in cluster",100,-50,50 );
  g_delCt = histoSvc()->book( "delCt", "del ct cut forming seg group",100,-20,20 );
  g_delZ0 = histoSvc()->book( "delZ0", "del z0 cut forming seg group",100,-60,60 );
  g_phiDiff = histoSvc()->book( "phiDiff", "phidiff mult for drop dup seg",100,-0.5,6.5 );
  g_slopeDiff = histoSvc()->book( "slopeDiff", "slopediff for drop dup seg",100,-0.3,0.3 );
  g_maxSegChisqAx = histoSvc()->book( "maxSegChisqAx", "max chisqDof ax seg combine" ,1000,0.,1000.);
  g_maxSegChisqSt = histoSvc()->book( "maxSegChisqSt", "max chisqDof st seg combine" ,200,-10.,200.);
  g_pickHitWire= histoSvc()->book( "pickHitWire", "nWire of pickHit" ,600,-288,288 );
 
  //for(int i=0;i<11;i++){
  //char title1[80],title2[80];
  //sprintf(title1,"segChi2Pat3%d",i);
  //sprintf(title2,"segChi2Pat4%d",i);
  //g_segChi2SlayPat3[i]  = histoSvc()->book( title1, "chi2/dof per layer of pat3",710,-10.,700. );
  //g_segChi2SlayPat4[i]  = histoSvc()->book( title2, "chi2/dof per layer of pat4",710,-10.,700. );
  //}
 
  //book tuple of wire efficiency
  NTuplePtr ntWireEff(ntupleSvc(), "MdcWire/mc");
  if ( ntWireEff ) { m_tupleWireEff = ntWireEff;}
  else {
    m_tupleWireEff = ntupleSvc()->book ("MdcWire/mc", CLID_ColumnWiseTuple, "MdcWire");
    if ( m_tupleWireEff ) {
      sc = m_tupleWireEff->addItem ("tkId",  m_we_tkId);
      sc = m_tupleWireEff->addItem ("pdg",   m_we_pdg);
      sc = m_tupleWireEff->addItem ("primary",   m_we_primary);
      sc = m_tupleWireEff->addItem ("layer",     m_we_layer);
      sc = m_tupleWireEff->addItem ("wire",      m_we_wire);
      sc = m_tupleWireEff->addItem ("gwire",     m_we_gwire);
      sc = m_tupleWireEff->addItem ("poisoned",  m_we_poisoned);
      sc = m_tupleWireEff->addItem ("seg",       m_we_seg);
      sc = m_tupleWireEff->addItem ("track",     m_we_track);
      sc = m_tupleWireEff->addItem ("pt",    m_we_pt);
      sc = m_tupleWireEff->addItem ("theta",     m_we_theta);
      sc = m_tupleWireEff->addItem ("phi",   m_we_phi);
      //sc = m_tupleWireEff->addItem ("tdc",       m_we_tdc);
    } else {   
      log << MSG::ERROR << "Cannot book tuple MdcWire/mc" << endmsg;
    }
  }
 
  //book tuple of test
  NTuplePtr ntt(ntupleSvc(), "MdcRec/test");
  if ( ntt ) { m_tuplet = ntt;}
  else {
    m_tuplet = ntupleSvc()->book ("MdcRec/test", CLID_ColumnWiseTuple, "MdcTrkRecon t particle");
    if ( m_tuplet ) {
      sc = m_tuplet->addItem ("layer",   m_tl);
      sc = m_tuplet->addItem ("wire",    m_tw);
      sc = m_tuplet->addItem ("phi",     m_tphi);
      sc = m_tuplet->addItem ("dphi",    m_tdphi);
      sc = m_tuplet->addItem ("dncell",  m_tdncell);
    } else {   
      log << MSG::ERROR << "Cannot book tuple MdcRec/test" << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of MC truth 
  NTuplePtr ntMc(ntupleSvc(), "MdcMc/mcTk");
  if ( ntMc ) { m_tupleMc = ntMc;}
  else {
    m_tupleMc = ntupleSvc()->book ("MdcMc/mcTk", CLID_ColumnWiseTuple, "MdcTrkRecon Mc particle");
    if ( m_tupleMc ) {
      sc = m_tupleMc->addItem ("evtNo",  m_tMcEvtNo);
      sc = m_tupleMc->addItem ("nDigi",  m_tMcNTk);
    } else {   
      log << MSG::ERROR << "Cannot book tuple MdcMc/mcTk" << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of MC truth 
  NTuplePtr ntMcHit(ntupleSvc(), "MdcMc/mcHit");
  if ( ntMcHit ) { m_tupleMcHit = ntMcHit;}
  else {
    m_tupleMcHit = ntupleSvc()->book ("MdcMc/mcHit", CLID_ColumnWiseTuple, "MdcTrkRecon Mc hit");
    if ( m_tupleMcHit ) {
      sc = m_tupleMcHit->addItem ("tkId", m_tMcTkId);
      sc = m_tupleMcHit->addItem ("pid",  m_tMcPid);
      sc = m_tupleMcHit->addItem ("px",   m_tMcPx);
      sc = m_tupleMcHit->addItem ("py",   m_tMcPy);
      sc = m_tupleMcHit->addItem ("pz",   m_tMcPz);
      sc = m_tupleMcHit->addItem ("d0",   m_tMcD0);
      sc = m_tupleMcHit->addItem ("z0",   m_tMcZ0);
      sc = m_tupleMcHit->addItem ("theta",m_tMcTheta);
      sc = m_tupleMcHit->addItem ("fiterm",m_tMcFiTerm);
      sc = m_tupleMcHit->addItem ("nMcHit",  m_tMcNHit, 0, 6796);
      sc = m_tupleMcHit->addIndexedItem ("layer",m_tMcNHit, m_tMcLayer);
      sc = m_tupleMcHit->addIndexedItem ("wire", m_tMcNHit, m_tMcWire);
      sc = m_tupleMcHit->addIndexedItem ("drift",m_tMcNHit, m_tMcDrift);
      sc = m_tupleMcHit->addIndexedItem ("x",    m_tMcNHit, m_tMcX);
      sc = m_tupleMcHit->addIndexedItem ("y",    m_tMcNHit, m_tMcY);
      sc = m_tupleMcHit->addIndexedItem ("z",    m_tMcNHit, m_tMcZ);
      sc = m_tupleMcHit->addIndexedItem ("lr",    m_tMcNHit, m_tMcLR);
    } else {   
      log << MSG::ERROR << "Cannot book tuple MdcMc/mcHit" << endmsg;
      //return StatusCode::FAILURE;
    }
  }
  //book tuple of recnstruction results
  NTuplePtr nt1(ntupleSvc(), "MdcRec/rec");
  if ( nt1 ) { m_tuple1 = nt1;}
  else {
    m_tuple1 = ntupleSvc()->book ("MdcRec/rec", CLID_ColumnWiseTuple, "MdcTrkRecon results");
    if ( m_tuple1 ) {
      sc = m_tuple1->addItem ("t0",      m_t0);
      sc = m_tuple1->addItem ("t0Stat",  m_t0Stat);
      sc = m_tuple1->addItem ("t0Truth", m_t0Truth);
 
      sc = m_tuple1->addItem ("nTdsTk",  m_nTdsTk);
      sc = m_tuple1->addItem ("nMcTk",   m_nMcTk);
 
      sc = m_tuple1->addItem ("p",       m_p);
      sc = m_tuple1->addItem ("pt",      m_pt);
      sc = m_tuple1->addItem ("nSlay",   m_nSlay);
      sc = m_tuple1->addItem ("pz",      m_pz);
      sc = m_tuple1->addItem ("d0",      m_d0);
      sc = m_tuple1->addItem ("phi0",    m_phi0);
      sc = m_tuple1->addItem ("cpa",     m_cpa);
      sc = m_tuple1->addItem ("z0",      m_z0);
      sc = m_tuple1->addItem ("tanl",    m_tanl);
      sc = m_tuple1->addItem ("q",       m_q);
      sc = m_tuple1->addItem ("pocax",   m_pocax);
      sc = m_tuple1->addItem ("pocay",   m_pocay);
      sc = m_tuple1->addItem ("pocaz",   m_pocaz);
 
      sc = m_tuple1->addItem ("evtNo",   m_evtNo);
      sc = m_tuple1->addItem ("nSt",     m_nSt);
      sc = m_tuple1->addItem ("nAct",    m_nAct);
      sc = m_tuple1->addItem ("nDof",    m_nDof);
      sc = m_tuple1->addItem ("chi2",    m_chi2);
      sc = m_tuple1->addItem ("tkId",    m_tkId);
      sc = m_tuple1->addItem ("nHit",    m_nHit, 0, 6796);
      sc = m_tuple1->addIndexedItem ("resid",    m_nHit, m_resid);
      sc = m_tuple1->addIndexedItem ("sigma",    m_nHit, m_sigma);
      sc = m_tuple1->addIndexedItem ("driftD",   m_nHit, m_driftD);
      sc = m_tuple1->addIndexedItem ("driftT",   m_nHit, m_driftT);
      sc = m_tuple1->addIndexedItem ("doca",     m_nHit, m_doca);
      sc = m_tuple1->addIndexedItem ("entra",    m_nHit, m_entra);
      sc = m_tuple1->addIndexedItem ("ambig",    m_nHit, m_ambig);
      sc = m_tuple1->addIndexedItem ("fltLen",   m_nHit, m_fltLen);
      sc = m_tuple1->addIndexedItem ("tof",      m_nHit, m_tof);
      sc = m_tuple1->addIndexedItem ("act",      m_nHit, m_act);
      sc = m_tuple1->addIndexedItem ("tdc",      m_nHit, m_tdc);
      sc = m_tuple1->addIndexedItem ("adc",      m_nHit, m_adc);
      sc = m_tuple1->addIndexedItem ("layer",    m_nHit, m_layer);
      sc = m_tuple1->addIndexedItem ("wire",     m_nHit, m_wire);
      sc = m_tuple1->addIndexedItem ("x",        m_nHit, m_x);
      sc = m_tuple1->addIndexedItem ("y",        m_nHit, m_y);
      sc = m_tuple1->addIndexedItem ("z",        m_nHit, m_z);
      sc = m_tuple1->addIndexedItem ("dx",       m_nHit, m_dx);
      sc = m_tuple1->addIndexedItem ("dy",       m_nHit, m_dy);
      sc = m_tuple1->addIndexedItem ("dz",       m_nHit, m_dz);
      sc = m_tuple1->addIndexedItem ("dDriftD",  m_nHit, m_dDriftD);
      sc = m_tuple1->addIndexedItem ("dlr",  m_nHit, m_dlr);
      sc = m_tuple1->addIndexedItem ("cellWidth",m_nHit, m_cellWidth);//for pickHits tuning
    } else {
      log << MSG::ERROR << "Cannot book tuple MdcRec/rec" << endmsg;
      //return StatusCode::FAILURE;
    }
  }
  //book tuple of segment
  NTuplePtr ntSeg(ntupleSvc(), "MdcSeg/seg");
  if ( ntSeg ) { m_tupleSeg = ntSeg;}
  else {
    m_tupleSeg = ntupleSvc()->book ("MdcSeg/seg", CLID_ColumnWiseTuple, "MdcTrkRecon segment data");
    if ( m_tupleSeg ) {
      sc = m_tupleSeg->addItem ("evtNo",     m_tsEvtNo);
      sc = m_tupleSeg->addItem ("nSeg",      m_tsNSeg);
      sc = m_tupleSeg->addItem ("nDigi",     m_tsNDigi, 0, 6796);
      sc = m_tupleSeg->addIndexedItem ("layer",  m_tsNDigi, m_tsLayer);
      sc = m_tupleSeg->addIndexedItem ("wire",   m_tsNDigi, m_tsWire);
      sc = m_tupleSeg->addIndexedItem ("inSeg",  m_tsNDigi, m_tsInSeg);
      sc = m_tupleSeg->addIndexedItem ("mcTkId", m_tsNDigi, m_tsMcTkId);
    } else {   
      log << MSG::ERROR << "Cannot book tuple MdcSeg/seg" << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of event
  NTuplePtr nt4(ntupleSvc(), "MdcRec/evt");
  if ( nt4 ) { m_tupleEvt = nt4;}
  else {
    m_tupleEvt = ntupleSvc()->book ("MdcRec/evt", CLID_ColumnWiseTuple, "MdcTrkRecon event data");
    if ( m_tupleEvt ) {
      sc = m_tupleEvt->addItem ("evtNo",     m_t4EvtNo);
      sc = m_tupleEvt->addItem ("nMcTk",     m_t4nMcTk );
      sc = m_tupleEvt->addItem ("nTdsTk",    m_t4nRecTk );
      sc = m_tupleEvt->addItem ("t0",        m_t4t0);
      sc = m_tupleEvt->addItem ("t0Stat",        m_t4t0Stat);
      sc = m_tupleEvt->addItem ("t0Truth",       m_t4t0Truth);
      sc = m_tupleEvt->addItem ("nDigiUn",       m_t4nDigiUnmatch);
      sc = m_tupleEvt->addItem ("nDigi",     m_t4nDigi, 0, 6796);
      sc = m_tupleEvt->addIndexedItem ("layer",  m_t4nDigi, m_t4Layer);
      sc = m_tupleEvt->addIndexedItem ("wire",   m_t4nDigi, m_t4Wire);
      sc = m_tupleEvt->addIndexedItem ("rt",     m_t4nDigi, m_t4Time);
      sc = m_tupleEvt->addIndexedItem ("rc",     m_t4nDigi, m_t4Charge);
      sc = m_tupleEvt->addIndexedItem ("phiMid", m_t4nDigi, m_t4PhiMid);
      sc = m_tupleEvt->addIndexedItem ("hot",    m_t4nDigi, m_t4Hot);
      //sc = m_tupleEvt->addIndexedItem ("recHit", m_t4nDigi, m_t4recHit);
      //sc = m_tupleEvt->addIndexedItem ("rawHit", m_t4nDigi, m_t4rawHit);
    } else {
      log << MSG::ERROR << "Cannot book N-tuple: MdcRec/evt"   << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of residula for every layer
  NTuplePtr ntCombAx(ntupleSvc(), "MdcCombAx/combAx");
  if ( ntCombAx ) { g_tupleCombAx= ntCombAx;}
  else { 
    g_tupleCombAx = ntupleSvc()->book ("MdcCombAx/combAx", CLID_RowWiseTuple, "MdcTrkRecon cuts in combine ax seg");
    if ( g_tupleCombAx) {
      sc = g_tupleCombAx->addItem ("dPhi0",     g_combAxdPhi0);
      sc = g_tupleCombAx->addItem ("dCurv",     g_combAxdCurv);
      sc = g_tupleCombAx->addItem ("sigPhi0",   g_combAxSigPhi0);
      sc = g_tupleCombAx->addItem ("sigCurv",   g_combAxSigCurv);
      sc = g_tupleCombAx->addItem ("slSeed",    g_combAxSlSeed);
      sc = g_tupleCombAx->addItem ("slTest",    g_combAxSlTest);
      sc = g_tupleCombAx->addItem ("qSeed",     g_combAxQualitySeed);
      sc = g_tupleCombAx->addItem ("qTest",     g_combAxQualityTest);
      sc = g_tupleCombAx->addItem ("pSeed",     g_combAxPatSeed);
      sc = g_tupleCombAx->addItem ("pTest",     g_combAxPatTest);
      sc = g_tupleCombAx->addItem ("nHitSeed",  g_combAxNhitSeed);
      sc = g_tupleCombAx->addItem ("nHitTest",  g_combAxNhitTest);
      sc = g_tupleCombAx->addItem ("mc",        g_combAxMc);
      sc = g_tupleCombAx->addItem ("ptTruth",   g_combAxMcPt);
      sc = g_tupleCombAx->addItem ("thetaTruth",g_combAxMcTheta);
      sc = g_tupleCombAx->addItem ("phiTruth",  g_combAxMcPhi);
      sc = g_tupleCombAx->addItem ("ambigSeed", g_combAxMcAmbigSeed);
      sc = g_tupleCombAx->addItem ("ambigTest", g_combAxMcAmbigTest);
    } else {   
      log << MSG::ERROR << "Cannot book N-tuple: FILE101/combAx"   << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of residula for every layer
  NTuplePtr ntFindSeg(ntupleSvc(), "MdcSeg/findSeg");
  if ( ntFindSeg ) { g_tupleFindSeg = ntFindSeg;}
  else { 
    g_tupleFindSeg = ntupleSvc()->book ("MdcSeg/findSeg", CLID_RowWiseTuple, "MdcTrkRecon cuts in segment finder");
    if ( g_tupleFindSeg) {
      sc = g_tupleFindSeg->addItem ("chi2",     g_findSegChi2);
      sc = g_tupleFindSeg->addItem ("slope",    g_findSegSlope);
      sc = g_tupleFindSeg->addItem ("intercept",g_findSegIntercept);
      sc = g_tupleFindSeg->addItem ("chi2Refit",g_findSegChi2Refit);
      sc = g_tupleFindSeg->addItem ("chi2Add",  g_findSegChi2Add);
      sc = g_tupleFindSeg->addItem ("pat",      g_findSegPat);
      sc = g_tupleFindSeg->addItem ("pat34",    g_findSegPat34);
      sc = g_tupleFindSeg->addItem ("nhit",     g_findSegNhit);
      sc = g_tupleFindSeg->addItem ("slayer",   g_findSegSl);
      sc = g_tupleFindSeg->addItem ("mc",   g_findSegMc);
      sc = g_tupleFindSeg->addItem ("ambig",    g_findSegAmbig);
    } else {   
      log << MSG::ERROR << "Cannot book N-tuple: FILE101/findSeg"   << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
  //book tuple of event
  NTuplePtr ntPick(ntupleSvc(), "MdcTrk/pick");
  if ( ntPick ) { m_tuplePick = ntPick;}
  else {
    m_tuplePick = ntupleSvc()->book ("MdcTrk/pick", CLID_ColumnWiseTuple, "MdcTrkRecon pickHits");
    if ( m_tuplePick ) {
      sc = m_tuplePick->addItem ("layer",           m_pickLayer);
      sc = m_tuplePick->addItem ("nCell",       m_pickNCell, 0, 288);
      sc = m_tuplePick->addItem ("nCellWithDigi",   m_pickNCellWithDigi, 0, 288);
      sc = m_tuplePick->addIndexedItem ("wire",    m_pickNCellWithDigi, m_pickWire);
      sc = m_tuplePick->addIndexedItem ("predDrift",    m_pickNCellWithDigi, m_pickPredDrift);
      sc = m_tuplePick->addIndexedItem ("drift",        m_pickNCellWithDigi, m_pickDrift);
      sc = m_tuplePick->addIndexedItem ("driftTruth",   m_pickNCellWithDigi, m_pickDriftTruth);
      sc = m_tuplePick->addIndexedItem ("phiZOk",       m_pickNCellWithDigi, m_pickPhizOk);
      sc = m_tuplePick->addIndexedItem ("z",            m_pickNCellWithDigi, m_pickZ);
      sc = m_tuplePick->addIndexedItem ("resid",        m_pickNCellWithDigi, m_pickResid);
      sc = m_tuplePick->addIndexedItem ("sigma",        m_pickNCellWithDigi, m_pickSigma);
      sc = m_tuplePick->addIndexedItem ("phiLowCut",    m_pickNCellWithDigi, m_pickPhiLowCut);
      sc = m_tuplePick->addIndexedItem ("phiHighCut",   m_pickNCellWithDigi, m_pickPhiHighCut);
      sc = m_tuplePick->addIndexedItem ("goodDriftCut", m_pickNCellWithDigi, m_pickDriftCut);
      sc = m_tuplePick->addIndexedItem ("mcTk",         m_pickNCellWithDigi, m_pickMcTk);
      sc = m_tuplePick->addIndexedItem ("is2d",         m_pickNCellWithDigi, m_pickIs2D);
      sc = m_tuplePick->addIndexedItem ("pt",           m_pickNCellWithDigi, m_pickPt);
      sc = m_tuplePick->addIndexedItem ("curv",           m_pickNCellWithDigi, m_pickCurv);
    } else {   
      log << MSG::ERROR << "Cannot book N-tuple: MdcTrk/pick"   << endmsg;
      //return StatusCode::FAILURE;
    }
  }
 
#ifdef MDCPATREC_TIMETEST
  //book tuple of time test
  NTuplePtr nt5(ntupleSvc(), "MdcTime/ti");
  if ( nt5 ) { m_tupleTime = nt5;}
  else {
    m_tupleTime = ntupleSvc()->book ("MdcTime/ti", CLID_RowWiseTuple, "MdcTrkRecon time");
    if ( m_tupleTime ) {
      sc = m_tupleTime->addItem ("eventtime",  m_eventTime);
      sc = m_tupleTime->addItem ("recTkNum",   m_t5recTkNum);
      sc = m_tupleTime->addItem ("mcTkNum",    m_mcTkNum);
      sc = m_tupleTime->addItem ("evtNo",      m_t5EvtNo);
      sc = m_tupleTime->addItem ("nHit",       m_t5nHit);
      sc = m_tupleTime->addItem ("nDigi",      m_t5nDigi);
    } else {
      log << MSG::ERROR << "Cannot book N-tuple: FILE101/time"   << endmsg;
      //return StatusCode::FAILURE;
    }
  }
  sc = service( "BesTimerSvc", m_timersvc);
  if( sc.isFailure() ) {
    log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
    //return StatusCode::FAILURE;
  } 
  m_timer[1] = m_timersvc->addItem("Execution");
  m_timer[1]->propName("nExecution");
#endif
 
#ifdef MDCPATREC_RESLAYER
  //book tuple of residula for every layer
  NTuplePtr nt6(ntupleSvc(), "MdcRes/res");
  if ( nt6 ) { g_tupleres = nt6;}
  else {
    g_tupleres= ntupleSvc()->book ("MdcRes/res", CLID_RowWiseTuple, "MdcTrkRecon res");
    if ( g_tupleres ) {
      sc = g_tupleres->addItem ("resid",  g_resLayer);
      sc = g_tupleres->addItem ("layer",  g_t6Layer); 
    } else {   
      log << MSG::ERROR << "Cannot book N-tuple: FILE101/res"   << endmsg;
      return StatusCode::FAILURE;
    }
  }
#endif
 
  return sc;
}// end of bookNTuple()

Referenced by initialize().

dumpDigi()

void MdcTrkRecon::dumpDigi ( )

protected

Definition at line 1294 of file MdcTrkRecon.cxx.

                          {
  uint32_t getDigiFlag = 0;
  getDigiFlag += m_maxMdcDigi;
  if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
  if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
  if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
  MdcDigiCol::iterator iter = mdcDigiVec.begin();
  std::cout<<name()<<" dump MdcDigiVec, nDigi="<<mdcDigiVec.size()<<std::endl;
  for (int iDigi=0;iter!= mdcDigiVec.end(); iter++,iDigi++ ) {
    long l = MdcID::layer((*iter)->identify());
    long w = MdcID::wire((*iter)->identify());
    int tkTruth = (*iter)->getTrackIndex();
    double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
    double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
    std::cout<<"("<<l<<","<<w<<";"<<tkTruth<<",t "<<tdc<<",c "<<adc<<")";
    if(iDigi%4==0) std::cout<<std::endl;
  }
  std::cout<<std::endl;
}

Referenced by execute().

dumpTdsTrack()

void MdcTrkRecon::dumpTdsTrack ( RecMdcTrackCol *  trackList )

protected

Definition at line 1072 of file MdcTrkRecon.cxx.

                                                       {
  std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
  RecMdcTrackCol::iterator it = trackList->begin();
  for (;it!= trackList->end();it++){
    RecMdcTrack *tk = *it;
    std::cout<< endl<<"//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
    cout <<" d0 "<<tk->helix(0)
      <<" phi0 "<<tk->helix(1)
      <<" cpa "<<tk->helix(2)
      <<" z0 "<<tk->helix(3)
      <<" tanl "<<tk->helix(4)
      <<endl;
    std::cout<<" q "<<tk->charge() 
      <<" theta "<<tk->theta()
      <<" phi "<<tk->phi()
      <<" x0 "<<tk->x()
      <<" y0 "<<tk->y()
      <<" z0 "<<tk->z()
      <<" r0 "<<tk->r()
      <<endl;
    std::cout <<" p "<<tk->p()
      <<" pt "<<tk->pxy()
      <<" px "<<tk->px()
      <<" py "<<tk->py()
      <<" pz "<<tk->pz()
      <<endl;
    std::cout<<" tkStat "<<tk->stat()
      <<" chi2 "<<tk->chi2()
      <<" ndof "<<tk->ndof()
      <<" nhit "<<tk->getNhits()
      <<" nst "<<tk->nster()
      <<endl;
    std::cout<< "errmat mat  " << std::endl;
    std::cout<< tk->err()<<std::endl;
    //std::cout<< "errmat array  " << std::endl;
    //for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
    //std::cout<< "   " << std::endl;
 
    int nhits = tk->getVecHits().size();
    std::cout<<nhits <<" Hits: " << std::endl;
    for(int ii=0; ii <nhits ; ii++){
      Identifier id(tk->getVecHits()[ii]->getMdcId());
      int layer = MdcID::layer(id);
      int wire = MdcID::wire(id);
      cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
    <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
    }//end of hit list
    std::cout << "  "<< std::endl;
  }//end of tk list
  std::cout << "  "<< std::endl;
}//end of dumpTdsTrack

execute()

StatusCode MdcTrkRecon::execute

(

)

Definition at line 215 of file MdcTrkRecon.cxx.

                                {
  setFilterPassed(false); 
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
 
  //Initialize
  if(m_hist>0){
    for (int ii=0;ii<43;ii++){ 
      for (int jj=0;jj<288;jj++){ 
    haveDigiTk[ii][jj] = -2; 
    haveDigiPt[ii][jj] = -2; 
    haveDigiTheta[ii][jj] = -999.; 
    haveDigiPhi[ii][jj] = -999.; 
    haveDigiDrift[ii][jj] = -999.; 
    haveDigiAmbig[ii][jj] = -999.; 
    recHitMap[ii][jj] = 0; 
    mcDrift[ii][jj] = -99.; 
    mcX[ii][jj] = -99.; 
    mcY[ii][jj] = -99.; 
    mcZ[ii][jj] = -99.; 
    mcLR[ii][jj] = -99.; 
    hitPoisoned[ii][jj]=-99;
      } 
    }
  }
 
 
  TrkContextEv context(m_bfield);
#ifdef MDCPATREC_TIMETEST
  m_timer[1]->start();
  ti_nHit=0;
#endif
  //------------------read event header--------------
  SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
  if (!evtHead) {
    log << MSG::FATAL<< "Could not retrieve event header" << endreq;
    return( StatusCode::FAILURE);
  }
  t_eventNo = evtHead->eventNumber(); 
  if(m_debug!=0)std::cout<<t_iExexute<<" run "<<evtHead->runNumber()<<" evt "<<t_eventNo<<std::endl;
  t_iExexute++;
 
  if (m_selEvtNo.size() >0){
    std::vector<int>::iterator it = m_selEvtNo.begin();
    for (; it!=m_selEvtNo.end(); it++){
      if((*it)== t_eventNo) setFilterPassed(true);
    }
  }
  //------------------get event start time-----------
  double t0 = 0.;
  t_t0 = -1.;
  t_t0Stat = -1;
  bool iterateT0 = false;
  const int nBunch = 3;//3 bunch/event
  double iBunch = 0 ; //form 0 to 2
  double BeamDeltaTime = 24. / nBunch;// 8ns
  SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
  if (!evTimeCol || evTimeCol->size()==0) { 
    log << MSG::WARNING <<" run "<<evtHead->runNumber() <<" evt "<<t_eventNo<< " Could not find RecEsTimeCol" << endreq;
    if(m_fieldCosmic){//yzhang temp for bes3 csmc test 
      return StatusCode::SUCCESS;
    } 
  }
  RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
  // skip RecEsTimeCol no rec data
  if (iter_evt != evTimeCol->end()){
    t_t0Stat = (*iter_evt)->getStat();
    t_t0 = (*iter_evt)->getTest();
 
    if (m_doLineFit){
      //t0 -= m_t0Const;
      if (abs(t_t0)<0.00001 || (t_t0 < 0) || (t_t0 > 9999.0) ){
    log << MSG::WARNING << "Skip with t0 = "<<t_t0 << endreq;
    return StatusCode::SUCCESS;//for bes3 cosmic test
      }
    }
    t0 =  t_t0*1.e-9;
  }
 
 
restart:
  if ( !m_recForEsTime && m_tryBunch) {
    if ( t_t0Stat < 10 ) return StatusCode::SUCCESS;
  }
  if ( m_tryBunch ){
    if ( iBunch > (nBunch - 1) ) return StatusCode::SUCCESS;
    if ( t_t0Stat < 0 ){ iterateT0 = true; }
    if ( iterateT0 ){
      //double EventDeltaTime = 24.;//24ns/event
      if ( (t_t0Stat > -1) && (fabs( iBunch * BeamDeltaTime - t_t0) < 2.) ){
    ++iBunch;
    goto restart;
      }else{ t0 = iBunch * BeamDeltaTime *1.e-9;}
      ++iBunch;
    }
  } 
 
  SmartDataPtr<MdcHitMap> hitMap(eventSvc(),"/Event/Hit/MdcHitMap");
  if (!hitMap) {
    log << MSG::FATAL << "Could not retrieve MDC hit map" << endreq;
    return( StatusCode::FAILURE);
  }
  //---------------------Read an event--------------
  SmartDataPtr<MdcHitCol> hitCol(eventSvc(),"/Event/Hit/MdcHitCol");
  if (!hitCol) {
    log << MSG::FATAL << "Could not retrieve MDC hit list" << endreq;
    return( StatusCode::FAILURE);
  }
  StatusCode sc;
 
  //---------- register RecMdcTrackCol and RecMdcHitCol to tds---------
  DataObject *aReconEvent;
  eventSvc()->findObject("/Event/Recon",aReconEvent);
  if(aReconEvent==NULL) {
    aReconEvent = new ReconEvent();
    StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
    if(sc!=StatusCode::SUCCESS) {
      log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
      return( StatusCode::FAILURE);
    } 
  }
 
  DataObject *aTrackCol;
  DataObject *aRecHitCol;
  //yzhang 2010-05-28 
  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;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
  RecMdcHitCol* hitList;
  eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
  if (aRecHitCol) {
    hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
  }else{
    hitList = new RecMdcHitCol;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
 
  m_hitData->loadevent(hitCol, hitMap, t0);
  t_nDigi = hitCol->size();
 
  if (poisoningHits()) {
    m_hitData->poisonHits(_gm,m_debug);
  }
 
  if((m_hist>0) && m_tupleWireEff){
    for(int i=0;i<m_hitData->nhits();i++){
      const MdcHit* thisHit = m_hitData->hit(i);
      int l = thisHit->layernumber();
      int w = thisHit->wirenumber();
      if(m_hitData->segUsage().get(thisHit)->deadHit()){
    hitPoisoned[l][w] = 1;
      }else{
    hitPoisoned[l][w] = 0;
      }
    }
  }
 
#ifdef MDCPATREC_TIMETEST
  SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
  if(!mcParticleCol){
    log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
  }
  m_mcTkNum = mcParticleCol->size();
#endif
 
  
  if(m_debug>1) dumpDigi();
  //--------------------------------------------------------------------------
  // Segment finding
  //--------------------------------------------------------------------------
  int nSegs = 0;
  if (m_flags.findSegs) {
    // Form track segments in superlayers
    nSegs = m_segFinder->createSegs(_gm, *m_segs, m_hitData->segUsage(), 
    m_hitData->hitMap(), t0);
  }
  log << MSG::INFO << " Found " << nSegs << " segments" << endreq;
  if (m_debug>1){
    std::cout<<"------------------------------------------------"<<std::endl;
    std::cout<<"==event "<<t_eventNo<< " Found " << nSegs << " segments" << std::endl;
    m_segs->plot();
  }
 
  if (m_flags.lHist||m_flags.segPar.lPrint) {
    fillSegList();
  } 
 
  //--------------------------------------------------------------------------
  // Combine segments to form track
  //-------------------------------------------------------------------------- 
  int nTracks  = 0;
  int nDeleted = 0;
  if (m_flags.findTracks && m_flags.findSegs) {
    if (nSegs > 2) {
      nTracks = m_tracks->createFromSegs(m_segs.get(), m_hitData->hitMap(),
      _gm, context, t0); 
    }
 
    if(m_arbitrateHits) nDeleted = m_tracks->arbitrateHits();
    int nKeep = nTracks - nDeleted;
 
    if (m_debug>0 && (nTracks - nDeleted)>0){
      std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
    <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
    << " keep  "<< nKeep <<" track(s)"; 
      if (nDeleted!=0){ std::cout <<",deleted " <<nDeleted; }
      std::cout<< std::endl;//yzhang debug
    }else{
      if (m_debug>0){
    std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
      <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
      <<" found 0 track "<< endl;
      }
    }
 
    if (m_flags.lHist) t_nRecTk = nKeep;
 
#ifdef MDCPATREC_RESLAYER
    m_tracks->setResLayer(m_resLayer);
#endif 
    if (m_flags.lHist){ fillTrackList(); }
    // Stick the found tracks into the list in TDS
    m_tracks->store(trackList, hitList);
    //if (m_debug>1) { dumpTdsTrack(trackList); }
    if(m_debug>1){
      std::cout<<name()<<" print track list  "<<std::endl;
      m_mdcPrintSvc->printRecMdcTrackCol();
    }
 
    //if(nKeep != (int)trackList->size()) std::cout<<__FILE__<<" "<<__LINE__<<"!!!!!!!!!!!!!!!!! nKeep != nTdsTk"<<  std::endl;
#ifdef MDCPATREC_TIMETEST
    RecMdcTrackCol::iterator iter = trackList->begin();
    for (;iter != trackList->end(); iter++) {
      MdcTrack* tk = *iter;
      ti_nHit += tk->getNhits(); 
    }
#endif
 
  } 
  //-------------End track-finding------------------ 
  m_segs->destroySegs();
 
  //if ( t_eventNo% 1000 == 0) { 
  //std::cout<<"evt "<<t_eventNo<< " Found " << trackList->size() << " track(s)" <<endl;//yzhang debug
  //}
  if(m_flags.lHist && m_mcHist) fillMcTruth();
#ifdef MDCPATREC_TIMETEST
  m_timer[1]->stop(); 
  //cout << "m_timer[1]->elapsed()::" << m_timer[1]->elapsed() << endl;
  //cout << "m_timer[1]->mean()::" << m_timer[1]->mean() << endl;
  m_eventTime = m_timer[1]->elapsed();
  m_t5recTkNum = m_tracks->length();
  m_t5EvtNo = t_eventNo;
  m_t5nHit = ti_nHit;
  m_t5nDigi = t_nDigi;
  sc = m_tupleTime->write();
#endif
  // for event start time, if track not found try another t0
  if ( m_tryBunch ){
    if ( nTracks <1 ){ iterateT0 = true; 
      goto restart;
    }
  }
  if (m_flags.lHist ) fillEvent();
 
  return StatusCode::SUCCESS;
}

fillEvent()

void MdcTrkRecon::fillEvent ( )

protected

Definition at line 1251 of file MdcTrkRecon.cxx.

                           {
  if (m_tupleEvt!=NULL){
 
    uint32_t getDigiFlag = 0;
    getDigiFlag += m_maxMdcDigi;
    if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
    if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
    if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
    MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
    if ((int)mdcDigiVec.size()<m_minMdcDigi){
      std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
    }
 
    m_t4nDigi = mdcDigiVec.size();
    int iDigi=0;
    MdcDigiCol::iterator iter = mdcDigiVec.begin();
    for (;iDigi<m_t4nDigi; iter++ ) {
      double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
      double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
      m_t4Time[iDigi] = tdc;
      m_t4Charge[iDigi] = adc;
      long l = MdcID::layer((*iter)->identify());
      long w = MdcID::wire((*iter)->identify());
      m_t4Layer[iDigi] = l;
      m_t4Wire[iDigi] = w;
      m_t4PhiMid[iDigi] = _gm->Layer(l)->phiWire(w);
      m_t4Hot[iDigi] = recHitMap[l][w];
      iDigi++;
    }
    m_t4t0 = t_t0;
    m_t4t0Stat = t_t0Stat;
    m_t4t0Truth = t_t0Truth;
    m_t4EvtNo =  t_eventNo;
    m_t4nRecTk = t_nRecTk;
    SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
    if (mdcDigiCol){
      m_t4nDigiUnmatch = mdcDigiCol->size();
    }
    m_t4nMcTk= t_mcTkNum;
    m_tupleEvt->write();
  }
}//end of fillEvent 

Referenced by execute().

fillMcTruth()

void MdcTrkRecon::fillMcTruth ( )

protected

Definition at line 878 of file MdcTrkRecon.cxx.

                             {
  MsgStream log(msgSvc(), name());
  StatusCode sc;
  t_mcTkNum = 0;
  for(int i=0;i<100;i++){
    isPrimaryOfMcTk[i]=-9999;
    pdgOfMcTk[i]=-9999;
  }
  double ptOfMcTk[100]={0.};
  double thetaOfMcTk[100]={0.};
  double phiOfMcTk[100]={0.};
  double nDigiOfMcTk[100]={0.};
  if (m_mcHist){
    //------------------Retrieve MC truth MdcMcHit------------
    SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol"); 
    if (!mcMdcMcHitCol) {
      log << MSG::INFO << "Could not find MdcMcHit" << endreq; 
    }else{
      Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
      for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
    const Identifier id= (*iter_mchit)->identify();
    int layer = MdcID::layer(id);
    int wire = MdcID::wire(id);
    int iMcTk =  (*iter_mchit)->getTrackIndex();
    if(iMcTk<100&&iMcTk>0) nDigiOfMcTk[iMcTk]++;
    mcDrift[layer][wire]= (*iter_mchit)->getDriftDistance()/10.;  //drift in MC.
    //std::cout<<"   "<<__FILE__<<" mcDrift  "<<mcDrift[layer][wire]<<"   "<<std::endl;
    mcX[layer][wire] = (*iter_mchit)->getPositionX()/10.;
    mcY[layer][wire] = (*iter_mchit)->getPositionY()/10.;
    mcZ[layer][wire] = (*iter_mchit)->getPositionZ()/10.;
    mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
    if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
    t_nHitInTk[iMcTk]++; 
    haveDigiAmbig[layer][wire] =  mcLR[layer][wire];
    haveDigiDrift[layer][wire] = mcDrift[layer][wire];
      }
    }
 
    //------------------get event start time truth-----------
    t_t0Truth = -10.;
    SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
    if(!mcParticleCol){
      log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
    }else {
      int nMcTk = 0;
      Event::McParticleCol::iterator it = mcParticleCol->begin();
      for (;it != mcParticleCol->end(); it++){ 
    //int pdg_code = (*it)->particleProperty();
    //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
    t_mcTkNum++;
    nMcTk++;
      }
      it = mcParticleCol->begin();
      for (;it != mcParticleCol->end(); it++){
    int tkId   = (*it)->trackIndex();
    if ((*it)->primaryParticle()){
      t_t0Truth = (*it)->initialPosition().t();
      isPrimaryOfMcTk[tkId] = 1;
    }else{
      isPrimaryOfMcTk[tkId] = 0;
      //continue;
    }
    //fill charged particle
    int pdg_code = (*it)->particleProperty();
    pdgOfMcTk[tkId] = pdg_code;
    //std::cout<<" tkId  "<<tkId<<" pdg_code  "<<pdg_code<<"   "<<std::endl;
    //FIXME skip charged track and track outside MDC
    //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
    const CLHEP::Hep3Vector& true_mom = (*it)->initialFourMomentum().vect();
    const CLHEP::HepLorentzVector& posIni = (*it)->initialPosition();
    const CLHEP::HepLorentzVector& posFin = (*it)->finalPosition();
    if(tkId<100&&tkId>=0) {
      ptOfMcTk[tkId] = true_mom.perp();
      thetaOfMcTk[tkId] = (*it)->initialFourMomentum().theta();
      phiOfMcTk[tkId] = (*it)->initialFourMomentum().phi();
    }
 
    if ( m_tupleMcHit){
      m_tMcPx    = true_mom.x();
      m_tMcPy    = true_mom.y();
      m_tMcPz    = true_mom.z();
      m_tMcD0    = posIni.perp()/10.;
      m_tMcZ0    = posIni.z()/10.;
      m_tMcTheta = (*it)->initialFourMomentum().theta();
      m_tMcFiTerm= posFin.phi();
      m_tMcPid   = pdg_code;
      m_tMcTkId  = tkId;
      m_tMcNHit  = t_nHitInTk[tkId];
      m_tupleMcHit->write();
    }
      }//end of loop mcParticleCol
      if(m_tupleMc) {
    m_tMcNTk= nMcTk;
    m_tMcEvtNo = t_eventNo;
    m_tupleMc->write();
      }
    }
  }//end of m_mcHist
 
  uint32_t getDigiFlag = 0;
  getDigiFlag += m_maxMdcDigi;
  if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
  if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
  if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
  if ((int)mdcDigiVec.size()<m_minMdcDigi){
    std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
  }
 
  if (mdcDigiVec.size()<=0) return;
  //fill raw digi and t0 status
  MdcDigiCol::iterator iter = mdcDigiVec.begin();
  for (;iter!= mdcDigiVec.end(); iter++ ) {
    long l = MdcID::layer((*iter)->identify());
    long w = MdcID::wire((*iter)->identify());
    int tkId = (*iter)->getTrackIndex();
    haveDigiTk[l][w]= tkId;
    //std::cout<<" l  "<<l<<" w  "<<w<<"  tk "<<tkId<<std::endl;
    haveDigiPt[l][w] = ptOfMcTk[(*iter)->getTrackIndex()];
    haveDigiTheta[l][w] = thetaOfMcTk[(*iter)->getTrackIndex()];
    haveDigiPhi[l][w] = phiOfMcTk[(*iter)->getTrackIndex()];
    if(m_tupleWireEff){
      m_we_tkId = tkId;
      if(tkId>=0){
    if(tkId>=1000) tkId-=1000;
    m_we_pdg = pdgOfMcTk[tkId];
    m_we_primary = isPrimaryOfMcTk[tkId];
    //if(pdgOfMcTk[tkId]==-22) cout<<" primaryParticle ? "<< isPrimaryOfMcTk[tkId]<< " tkId "<<tkId <<" layer "<<l<<" wire "<<w<<" pdg="<<pdgOfMcTk[tkId]<<endl;
      }else{
    m_we_pdg = -999;
    m_we_primary = -999;
      }
      m_we_layer = l;
      m_we_wire = w;
      int gwire = w+Constants::nWireBeforeLayer[l];
      m_we_gwire = gwire;
      m_we_poisoned = hitPoisoned[l][w];
      if(hitInSegList[l][w]>0) m_we_seg = 1;
      else m_we_seg = 0;
      if(recHitMap[l][w]>0) m_we_track = 1;
      else m_we_track = 0;
      if(l==35&&tkId>=0&&abs(pdgOfMcTk[tkId])==11&&hitInSegList[l][w]<=0) {
    cout<<"layer "<<l<<" cell "<<w<<" gwire "<<gwire<<" inseg "<<hitInSegList[l][w]<<endl;
      }
      m_we_pt = ptOfMcTk[tkId];
      m_we_theta = thetaOfMcTk[tkId];
      m_we_phi = phiOfMcTk[tkId];
      m_tupleWireEff->write();
    }
  }
}//end of fillMcTruth()

Referenced by execute().

fillSegList()

void MdcTrkRecon::fillSegList ( )

protected

Definition at line 1030 of file MdcTrkRecon.cxx.

                              {
  if (2 == m_flags.segPar.lPrint) {
    std::cout << "print after Segment finding " << std::endl;
  }
  if (!m_flags.lHist) return;
  // Fill hits of every layer after segment finding
  for (int ii=0;ii<43;ii++){
    for (int jj=0;jj<288;jj++){ hitInSegList[ii][jj] = 0; }
  }
  int nSeg=0;
  for (int i = 0; i < _gm->nSuper(); i++) {
    MdcSeg* aSeg = (MdcSeg *) m_segs->oneList(i)->first();
    while (aSeg != NULL) {
      nSeg++;
      for (int ihit=0 ; ihit< aSeg->nHit() ; ihit++){
    const MdcHit* hit = aSeg->hit(ihit)->mdcHit();
    int layer = hit->layernumber();
    int wire = hit->wirenumber();
    hitInSegList[layer][wire]++;
      }
      aSeg = (MdcSeg *) aSeg->next();
    }
  }//end for slayer
  if (!m_tupleSeg) return;
  m_tsEvtNo =  t_eventNo;
  m_tsNDigi = t_nDigi;
  int iDigi=0;
  for (int ii=0;ii<43;ii++){
    for (int jj=0;jj<288;jj++){ 
      if (haveDigiTk[ii][jj] > -2){
    m_tsLayer[iDigi] = ii;
    m_tsWire[iDigi] = jj;
    m_tsInSeg[iDigi] = hitInSegList[ii][jj];
    m_tsMcTkId[iDigi] = haveDigiTk[ii][jj];
    iDigi++;
      }
    }
  }
  m_tsNSeg=nSeg;
  m_tupleSeg->write();
}//end of fillSegList

Referenced by execute().

fillTrackList()

void MdcTrkRecon::fillTrackList ( )

protected

Definition at line 1124 of file MdcTrkRecon.cxx.

                               { 
  if (m_flags.debugFlag()>1) {
    std::cout << "=======Print after Track finding=======" << std::endl;
    m_tracks->plot();
  }
 
  if (!m_flags.lHist) return; 
  //----------- fill hit -----------
  int nTk = (*m_tracks).nTrack();
  for (int itrack = 0; itrack < nTk; itrack++) {
    MdcTrack *atrack = (*m_tracks)[itrack];
    if (atrack==NULL) continue;
    const TrkFit* fitResult = atrack->track().fitResult();
    if (fitResult == 0) continue;//check the fit worked
 
    //for fill ntuples
    int nSt=0; //int nSeg=0; 
    int seg[11] = {0}; int segme;
    //-----------hit list-------------
    const TrkHitList* hitlist = atrack->track().hits();
    TrkHitList::hot_iterator hot = hitlist->begin();
    int layerCount[43]={0};
    for(int iii=0;iii<43;iii++){layerCount[iii]=0;}
    int i=0;
    for (;hot!= hitlist->end();hot++){
      const MdcRecoHitOnTrack* recoHot;           
      recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
      int layer = recoHot->mdcHit()->layernumber();
      int wire  = recoHot->mdcHit()->wirenumber();
      layerCount[layer]++;
      recHitMap[layer][wire]++;
      //for n seg
      segme=0;
      if ( layer >0 ) {segme=(layer-1)/4;}
      seg[segme]++;
      if (recoHot->layer()->view()) { ++nSt; }
 
      if(!m_tuple1) continue;
      m_layer[i]    = layer;
      m_wire[i]     = wire;
      m_ambig[i]    = recoHot->wireAmbig();// wire ambig
      //fltLen
      double fltLen = recoHot->fltLen(); 
      m_fltLen[i]   = fltLen;
      double tof    = recoHot->getParentRep()->arrivalTime(fltLen);
      //position
      HepPoint3D pos; Hep3Vector dir;    
      recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
      m_x[i] = pos.x();      
      m_y[i] = pos.y();      
      m_z[i] = pos.z();      
      m_driftT[i]   = recoHot->mdcHit()->driftTime(tof,pos.z());
      m_tof[i]      = tof;
      m_driftD[i]   = recoHot->drift();
      m_sigma[i]    = recoHot->hitRms();        
      m_tdc[i]      = recoHot->rawTime();
      m_adc[i]      = recoHot->mdcHit()->charge(); 
      m_doca[i]     = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME 
      m_entra[i]    = recoHot->entranceAngle();
      m_act[i]      = recoHot->isActive();
      //diff with truth
      m_dx[i] = m_x[i] - mcX[layer][wire];
      m_dy[i] = m_y[i] - mcY[layer][wire];
      m_dz[i] = m_z[i] - mcZ[layer][wire];
      m_dDriftD[i] = m_driftD[i] - mcDrift[layer][wire];
      m_dlr[i] = m_ambig[i] - mcLR[layer][wire];
      //yzhang for pickHits debug
      m_cellWidth[i] = Constants::twoPi*_gm->Layer(layer)->rMid()/_gm->Layer(layer)->nWires();
      //zhangy
      //resid
      double res=999.,rese=999.;
      if (recoHot->resid(res,rese,false)){
      }else{}
      m_resid[i] = res;
      i++;
    }// end fill hit
    int nSlay=0;
    for(int i=0;i<11;i++){
      if (seg[i]>0) nSlay++;
    }
    if(m_tuple1){
      //------------fill track result-------------
      m_tkId = itrack;
      //track parameters at closest approach to beamline
      double fltLenPoca = 0.0;
      TrkExchangePar helix = fitResult->helix(fltLenPoca);
      m_q    = fitResult->charge();
      m_phi0 = helix.phi0();
      m_tanl = helix.tanDip();
      m_z0 = helix.z0();
      m_d0 = helix.d0();
      m_pt = fitResult->pt();
      m_nSlay = nSlay;
      if (m_pt > 0.00001){
    m_cpa = fitResult->charge()/fitResult->pt();
      }
      //momenta and position
      CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
      double px,py,pz,pxy;
      pxy = fitResult->pt();
      px = p1.x();
      py = p1.y();
      pz = p1.z();
      m_p = p1.mag();
      m_pz = pz;
      HepPoint3D poca = fitResult->position(fltLenPoca);
      m_pocax = poca.x();
      m_pocay = poca.y();
      m_pocaz = poca.z(); 
      m_nAct = fitResult->nActive();
      m_nHit = hitlist->nHit();
      m_nDof = fitResult->nDof();
      m_nSt  = nSt;
      m_chi2 = fitResult->chisq();
      //for (int l=0;l<43;l++) m_layerCount[l] = layerCount[l];
      m_t0 = t_t0;
      m_t0Stat = t_t0Stat;
      m_t0Truth = t_t0Truth;
      m_nTdsTk = t_nRecTk;
      m_nMcTk = t_mcTkNum;
      m_evtNo = t_eventNo;
      m_tuple1->write();
    }
  }//end of loop rec trk list
}//end of fillTrackList

Referenced by execute().

finalize()

StatusCode MdcTrkRecon::finalize

(

)

Definition at line 512 of file MdcTrkRecon.cxx.

                                 {
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in finalize()" << endreq;
 
  m_tracks.reset(0);
#ifdef MDCPATREC_TIMETEST
  m_timersvc->print();
#endif
  return StatusCode::SUCCESS;
}

ignoringUsedHits()

bool MdcTrkRecon::ignoringUsedHits ( ) const

inlineprotected

Definition at line 41 of file MdcTrkRecon.h.

{return m_onlyUnusedHits;}

Referenced by initialize().

initialize()

StatusCode MdcTrkRecon::initialize

(

)

Definition at line 138 of file MdcTrkRecon.cxx.

                                  {
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in initialize()" << endreq;
  StatusCode sc;
 
  //Pdt
  Pdt::readMCppTable(m_pdtFile);
 
  //Flag and Pars
  m_flags.readPar(m_paramFile);
  m_flags.lHist = m_hist;
  m_flags.setDebug(m_debug);
  for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
  if (!m_doLineFit) MdcTrackListBase::m_d0Cut = m_d0Cut;
  if (!m_doLineFit) MdcTrackListBase::m_z0Cut = m_z0Cut;
  MdcTrackListBase::m_ptCut = m_dropTrkPt;
  if (m_debug>0) {
    std::cout<< "MdcTrkRecon d0Cut "<<m_d0Cut<< " z0Cut "<< 
      m_z0Cut<<" ptCut "<<m_dropTrkPt << std::endl;
  }
 
  //Initailize magnetic filed 
  sc = service ("MagneticFieldSvc",m_pIMF);
  if(sc!=StatusCode::SUCCESS) {
    log << MSG::ERROR << name() << " Unable to open magnetic field service"<<endreq;
  }
  m_bfield = new BField(m_pIMF);
  log << MSG::INFO <<name() << " field z = "<<m_bfield->bFieldNominal()<< endreq;
 
  //Initialize RawDataProviderSvc
  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;
  } 
 
  //Initailize 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;
  }
 
  //Initialize hitdata
  m_hitData.reset( new MdcSegData( ignoringUsedHits() ));
 
  //Initialize segFinder
  if (m_flags.findSegs) {
    m_segFinder.reset( new MdcSegFinder(m_flags.segPar.useAllAmbig) );
  }
 
  //Initialize trackList
  if ( m_doLineFit ){
    m_tracks.reset( new MdcTrackListCsmc(m_flags.tkParTight) );
  } else {
    m_tracks.reset( new MdcTrackList(m_flags.tkParTight) );
  }
 
  //Book NTuple and Histograms
  if (m_flags.lHist){ 
    StatusCode sc = bookNTuple();
    if (!sc.isSuccess()) { m_flags.lHist = 0;}
  }
 
  //yzhang for HelixFitter debug
  if(7 == m_flags.debugFlag())TrkHelixFitter::m_debug = true;
 
  t_iExexute = 0;
  return StatusCode::SUCCESS;
}

poisoningHits()

bool MdcTrkRecon::poisoningHits ( ) const

inlineprotected

Definition at line 42 of file MdcTrkRecon.h.

{return m_poisonHits;}

Referenced by execute().

---

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

• MdcTrkRecon.h
• MdcTrkRecon.cxx