MdcTrkRecon Class Reference

#include <MdcTrkRecon.h>

List of all members.

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

Private Attributes

int t_iExecute
int m_fittingMethod
int m_allHit
std::string m_configFile
std::string m_paramFile
std::string m_pdtFile
const MdcDetector_gm
MdcFlagHold m_flags
std::auto_ptr< MdcSegDatam_hitData
std::auto_ptr< MdcSegListm_segs
std::auto_ptr< MdcTrackListBasem_tracks
std::auto_ptr< MdcSegFinderm_segFinder
RawDataProviderSvcm_rawDataProviderSvc
MdcPrintSvcm_mdcPrintSvc
BFieldm_bfield
std::vector< int > m_selEvtNo
bool m_onlyUnusedHits
bool m_poisonHits
bool m_doLineFit
bool m_tryBunch
bool m_recForEsTime
std::vector< float > m_helixHitsSigma
uint32_t m_getDigiFlag
int m_maxMdcDigi
bool m_keepBadTdc
bool m_dropHot
bool m_keepUnmatch
int m_minMdcDigi
bool m_fieldCosmic
double m_d0Cut
double m_z0Cut
double m_dropTrkPt
bool m_mcHist
int m_hist
int m_debug
bool m_doSagFlag
bool m_arbitrateHits
bool m_combineTracking
double mcDrift [43][288]
double mcLayer [43][288]
double mcWire [43][288]
double mcLR [43][288]
double mcX [43][288]
double mcY [43][288]
double mcZ [43][288]
double hitOnSeg [43]
double hitOnMcTk [43]
int hitPoisoned [43][288]
int hitInSegList [43][288]
int t_nHitInTk [100]
int t_t0Stat
double t_t0
double t_t0Truth
double t_mcTkNum
int t_nDigi
int t_eventNo
int t_nRecTk
int t_iExexute
int isPrimaryOfMcTk [100]
int pdgOfMcTk [100]
IMagneticFieldSvcm_pIMF


Detailed Description

Definition at line 31 of file MdcTrkRecon.h.


Constructor & Destructor Documentation

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

Definition at line 71 of file MdcTrkRecon.cxx.

References m_allHit, m_arbitrateHits, m_combineTracking, m_configFile, m_d0Cut, m_debug, m_doLineFit, m_doSagFlag, m_dropHot, m_dropTrkPt, m_fieldCosmic, m_fittingMethod, m_getDigiFlag, m_helixHitsSigma, m_hist, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_mcHist, m_minMdcDigi, m_onlyUnusedHits, m_paramFile, m_pdtFile, m_poisonHits, m_recForEsTime, m_selEvtNo, m_tryBunch, and m_z0Cut.

00071                                                                         :
00072   Algorithm(name, pSvcLocator),
00073   m_hitData(0), m_segs(0), m_tracks(0), m_segFinder(0)  
00074 { 
00075   //Declare the properties--------------------------------
00076   declareProperty("FittingMethod",  m_fittingMethod = 2);
00077   declareProperty("ConfigFile",     m_configFile = "MDCConfig.xml");
00078   declareProperty("OnlyUnusedHits", m_onlyUnusedHits = 0);
00079   declareProperty("PoisonHits",     m_poisonHits = 0);
00080   declareProperty("doLineFit",      m_doLineFit = false);       
00081   declareProperty("paramFile",      m_paramFile = "params"); 
00082   declareProperty("pdtFile",        m_pdtFile = "pdt.table");
00083   declareProperty("allHit",         m_allHit = 1);
00084   declareProperty("hist",           m_hist= false);
00085   declareProperty("recForEsTime",   m_recForEsTime= false);
00086   declareProperty("tryBunch",       m_tryBunch = false);
00087   declareProperty("d0Cut",          m_d0Cut= -999.);
00088   declareProperty("z0Cut",          m_z0Cut= -999.);
00089   declareProperty("dropTrkPt",      m_dropTrkPt = -999.);
00090   declareProperty("debug",          m_debug= 0); 
00091   declareProperty("mcHist",         m_mcHist= false); 
00092   declareProperty("fieldCosmic",    m_fieldCosmic= false); 
00093   declareProperty("doSag",          m_doSagFlag= false);
00094   declareProperty("arbitrateHits",  m_arbitrateHits = true);
00095 
00096   declareProperty("helixHitsSigma", m_helixHitsSigma);
00097   //for fill hist 
00098   declareProperty("getDigiFlag",    m_getDigiFlag = 0);
00099   declareProperty("maxMdcDigi",     m_maxMdcDigi = 0);
00100   declareProperty("keepBadTdc",     m_keepBadTdc = 0);
00101   declareProperty("dropHot",        m_dropHot= 0);
00102   declareProperty("keepUnmatch",    m_keepUnmatch = 0);
00103   declareProperty("minMdcDigi",     m_minMdcDigi = 0);
00104   declareProperty("selEvtNo",       m_selEvtNo);
00105   declareProperty("combineTracking",m_combineTracking = false);//yzhang 2010-05-28 
00106 
00107 #ifdef MDCPATREC_RESLAYER
00108   declareProperty("resLayer",m_resLayer= -1);
00109 #endif
00110 
00111 }

MdcTrkRecon::~MdcTrkRecon (  ) 

Definition at line 113 of file MdcTrkRecon.cxx.

References m_hitData, m_segFinder, m_segs, and m_tracks.

00114 { 
00115   m_segs.reset(0);
00116   m_segFinder.reset(0);
00117   m_hitData.reset(0);
00118   m_tracks.reset(0); 
00119 }


Member Function Documentation

StatusCode MdcTrkRecon::beginRun (  ) 

Definition at line 122 of file MdcTrkRecon.cxx.

References _gm, Bes_Common::INFO, MdcDetector::instance(), m_doSagFlag, m_flags, m_segs, msgSvc(), MdcDetector::nSuper(), and MdcFlagHold::segPar.

00122                                 {
00123   MsgStream log(msgSvc(), name());  
00124   log << MSG::INFO << "in beginRun()" << endreq;
00125 
00126   //Detector geometry 
00127   _gm = MdcDetector::instance(m_doSagFlag);
00128 
00129   if(NULL == _gm) return StatusCode::FAILURE;
00130 
00131   //Initialize segList
00132   m_segs.reset( new MdcSegList(_gm->nSuper(), m_flags.segPar) );
00133 
00134   return StatusCode::SUCCESS;
00135 }

StatusCode MdcTrkRecon::bookNTuple (  )  [protected]

Definition at line 523 of file MdcTrkRecon.cxx.

References calibUtil::ERROR, g_3dTkChi2, g_cirTkChi2, g_combAxdCurv, g_combAxdPhi0, g_combAxMc, g_combAxMcAmbigSeed, g_combAxMcAmbigTest, g_combAxMcPhi, g_combAxMcPt, g_combAxMcTheta, g_combAxNhitSeed, g_combAxNhitTest, g_combAxPatSeed, g_combAxPatTest, g_combAxQualitySeed, g_combAxQualityTest, g_combAxSigCurv, g_combAxSigPhi0, g_combAxSlSeed, g_combAxSlTest, g_ctCut, g_delCt, g_delZ0, g_findSegAmbig, g_findSegChi2, g_findSegChi2Add, g_findSegChi2Refit, g_findSegIntercept, g_findSegMc, g_findSegNhit, g_findSegPat, g_findSegPat34, g_findSegSl, g_findSegSlope, g_maxSegChisqAx, g_maxSegChisqSt, g_nSigAdd, g_phiDiff, g_pickHitWire, g_residVsDoca, g_residVsLayer, g_slopeDiff, g_tupleCombAx, g_tupleFindSeg, g_z0Cut, h_sfHit, histoSvc(), m_act, m_adc, m_ambig, m_cellWidth, m_chi2, m_cpa, m_d0, m_dDriftD, m_dlr, m_doca, m_driftD, m_driftT, m_dx, m_dy, m_dz, m_entra, m_evtNo, m_fltLen, m_layer, m_nAct, m_nDof, m_nHit, m_nMcTk, m_nSlay, m_nSt, m_nTdsTk, m_p, m_phi0, m_pickCurv, m_pickDrift, m_pickDriftCut, m_pickDriftTruth, m_pickIs2D, m_pickLayer, m_pickMcTk, m_pickNCell, m_pickNCellWithDigi, m_pickPhiHighCut, m_pickPhiLowCut, m_pickPhizOk, m_pickPredDrift, m_pickPt, m_pickResid, m_pickSigma, m_pickWire, m_pickZ, m_pocax, m_pocay, m_pocaz, m_pt, m_pz, m_q, m_resid, m_sigma, m_t0, m_t0Stat, m_t0Truth, m_t4Charge, m_t4EvtNo, m_t4Hot, m_t4Layer, m_t4nDigi, m_t4nDigiUnmatch, m_t4nMcTk, m_t4nRecTk, m_t4PhiMid, m_t4t0, m_t4t0Stat, m_t4t0Truth, m_t4Time, m_t4Wire, m_tanl, m_tdc, m_tdncell, m_tdphi, m_tkId, m_tl, m_tMcD0, m_tMcDrift, m_tMcEvtNo, m_tMcFiTerm, m_tMcLayer, m_tMcLR, m_tMcNHit, m_tMcNTk, m_tMcPid, m_tMcPx, m_tMcPy, m_tMcPz, m_tMcTheta, m_tMcTkId, m_tMcWire, m_tMcX, m_tMcY, m_tMcZ, m_tMcZ0, m_tof, m_tphi, m_tsEvtNo, m_tsInSeg, m_tsLayer, m_tsMcTkId, m_tsNDigi, m_tsNSeg, m_tsWire, m_tuple1, m_tupleEvt, m_tupleMc, m_tupleMcHit, m_tuplePick, m_tupleSeg, m_tuplet, m_tupleWireEff, m_tw, m_we_gwire, m_we_layer, m_we_pdg, m_we_phi, m_we_poisoned, m_we_primary, m_we_pt, m_we_seg, m_we_theta, m_we_tkId, m_we_track, m_we_wire, m_wire, m_x, m_y, m_z, m_z0, msgSvc(), ntupleSvc(), and Bes_Common::WARNING.

Referenced by initialize().

00523                                   {
00524   MsgStream log(msgSvc(), name());
00525   StatusCode sc = StatusCode::SUCCESS;
00526 
00527 
00528   //--------------book histogram------------------  
00529   h_sfHit  = histoSvc()->book( "sfHit", "Hits after segment finding",46,-1.,44. );
00530   //g_actHit = histoSvc()->book( "actHit", "Active hits",46,-1.,44. );
00531   //g_hitEff = histoSvc()->book( "hitEff", "Hit efficiency in the event",100,-0.5,1.2 );
00532   g_residVsLayer = histoSvc()->book( "residVsLayer", "Residual vs. layer",60,-5,50.,100,-0.5,1.2 );
00533   g_residVsDoca  = histoSvc()->book( "residVsDoca", "Residial vs. doca",50,-2.,2.,100,-0.5,1.2 );
00534   //segment
00535   //g_segChi2  = histoSvc()->book( "segChi2", "chi2 per dof in segment finding",101,-1.,100. );
00536   g_cirTkChi2  = histoSvc()->book( "cirTkChi2", "chi2 per dof in circle finding",21,-1. ,20. );
00537   g_3dTkChi2 = histoSvc()->book( "chi2Helix", "maxChisq dof helix fit" ,51.,-1.,50);
00538   g_nSigAdd = histoSvc()->book( "nSigAdd", "max allowed n sigma to add hit to existing segment",50,0,50 );
00539   g_z0Cut = histoSvc()->book( "z0Cut", "z0 for including stereo segs in cluster",100,0,600 );
00540   g_ctCut = histoSvc()->book( "ctCut", "ct for including stereo segs in cluster",100,-50,50 );
00541   g_delCt = histoSvc()->book( "delCt", "del ct cut forming seg group",100,-20,20 );
00542   g_delZ0 = histoSvc()->book( "delZ0", "del z0 cut forming seg group",100,-60,60 );
00543   g_phiDiff = histoSvc()->book( "phiDiff", "phidiff mult for drop dup seg",100,-0.5,6.5 );
00544   g_slopeDiff = histoSvc()->book( "slopeDiff", "slopediff for drop dup seg",100,-0.3,0.3 );
00545   g_maxSegChisqAx = histoSvc()->book( "maxSegChisqAx", "max chisqDof ax seg combine" ,1000,0.,1000.);
00546   g_maxSegChisqSt = histoSvc()->book( "maxSegChisqSt", "max chisqDof st seg combine" ,200,-10.,200.);
00547   g_pickHitWire= histoSvc()->book( "pickHitWire", "nWire of pickHit" ,600,-288,288 );
00548 
00549   //for(int i=0;i<11;i++){
00550   //char title1[80],title2[80];
00551   //sprintf(title1,"segChi2Pat3%d",i);
00552   //sprintf(title2,"segChi2Pat4%d",i);
00553   //g_segChi2SlayPat3[i]  = histoSvc()->book( title1, "chi2/dof per layer of pat3",710,-10.,700. );
00554   //g_segChi2SlayPat4[i]  = histoSvc()->book( title2, "chi2/dof per layer of pat4",710,-10.,700. );
00555   //}
00556 
00557   //book tuple of wire efficiency
00558   NTuplePtr ntWireEff(ntupleSvc(), "MdcWire/mc");
00559   if ( ntWireEff ) { m_tupleWireEff = ntWireEff;}
00560   else {
00561     m_tupleWireEff = ntupleSvc()->book ("MdcWire/mc", CLID_ColumnWiseTuple, "MdcWire");
00562     if ( m_tupleWireEff ) {
00563       sc = m_tupleWireEff->addItem ("tkId",      m_we_tkId);
00564       sc = m_tupleWireEff->addItem ("pdg",       m_we_pdg);
00565       sc = m_tupleWireEff->addItem ("primary",   m_we_primary);
00566       sc = m_tupleWireEff->addItem ("layer",     m_we_layer);
00567       sc = m_tupleWireEff->addItem ("wire",      m_we_wire);
00568       sc = m_tupleWireEff->addItem ("gwire",     m_we_gwire);
00569       sc = m_tupleWireEff->addItem ("poisoned",  m_we_poisoned);
00570       sc = m_tupleWireEff->addItem ("seg",       m_we_seg);
00571       sc = m_tupleWireEff->addItem ("track",     m_we_track);
00572       sc = m_tupleWireEff->addItem ("pt",        m_we_pt);
00573       sc = m_tupleWireEff->addItem ("theta",     m_we_theta);
00574       sc = m_tupleWireEff->addItem ("phi",       m_we_phi);
00575       //sc = m_tupleWireEff->addItem ("tdc",       m_we_tdc);
00576     } else {   
00577       log << MSG::ERROR << "Cannot book tuple MdcWire/mc" << endmsg;
00578     }
00579   }
00580 
00581   //book tuple of test
00582   NTuplePtr ntt(ntupleSvc(), "MdcRec/test");
00583   if ( ntt ) { m_tuplet = ntt;}
00584   else {
00585     m_tuplet = ntupleSvc()->book ("MdcRec/test", CLID_ColumnWiseTuple, "MdcTrkRecon t particle");
00586     if ( m_tuplet ) {
00587       sc = m_tuplet->addItem ("layer",   m_tl);
00588       sc = m_tuplet->addItem ("wire",    m_tw);
00589       sc = m_tuplet->addItem ("phi",     m_tphi);
00590       sc = m_tuplet->addItem ("dphi",    m_tdphi);
00591       sc = m_tuplet->addItem ("dncell",  m_tdncell);
00592     } else {   
00593       log << MSG::ERROR << "Cannot book tuple MdcRec/test" << endmsg;
00594       //return StatusCode::FAILURE;
00595     }
00596   }
00597 
00598   //book tuple of MC truth 
00599   NTuplePtr ntMc(ntupleSvc(), "MdcMc/mcTk");
00600   if ( ntMc ) { m_tupleMc = ntMc;}
00601   else {
00602     m_tupleMc = ntupleSvc()->book ("MdcMc/mcTk", CLID_ColumnWiseTuple, "MdcTrkRecon Mc particle");
00603     if ( m_tupleMc ) {
00604       sc = m_tupleMc->addItem ("evtNo",  m_tMcEvtNo);
00605       sc = m_tupleMc->addItem ("nDigi",  m_tMcNTk);
00606     } else {   
00607       log << MSG::ERROR << "Cannot book tuple MdcMc/mcTk" << endmsg;
00608       //return StatusCode::FAILURE;
00609     }
00610   }
00611 
00612   //book tuple of MC truth 
00613   NTuplePtr ntMcHit(ntupleSvc(), "MdcMc/mcHit");
00614   if ( ntMcHit ) { m_tupleMcHit = ntMcHit;}
00615   else {
00616     m_tupleMcHit = ntupleSvc()->book ("MdcMc/mcHit", CLID_ColumnWiseTuple, "MdcTrkRecon Mc hit");
00617     if ( m_tupleMcHit ) {
00618       sc = m_tupleMcHit->addItem ("tkId", m_tMcTkId);
00619       sc = m_tupleMcHit->addItem ("pid",  m_tMcPid);
00620       sc = m_tupleMcHit->addItem ("px",   m_tMcPx);
00621       sc = m_tupleMcHit->addItem ("py",   m_tMcPy);
00622       sc = m_tupleMcHit->addItem ("pz",   m_tMcPz);
00623       sc = m_tupleMcHit->addItem ("d0",   m_tMcD0);
00624       sc = m_tupleMcHit->addItem ("z0",   m_tMcZ0);
00625       sc = m_tupleMcHit->addItem ("theta",m_tMcTheta);
00626       sc = m_tupleMcHit->addItem ("fiterm",m_tMcFiTerm);
00627       sc = m_tupleMcHit->addItem ("nMcHit",      m_tMcNHit, 0, 6796);
00628       sc = m_tupleMcHit->addIndexedItem ("layer",m_tMcNHit, m_tMcLayer);
00629       sc = m_tupleMcHit->addIndexedItem ("wire", m_tMcNHit, m_tMcWire);
00630       sc = m_tupleMcHit->addIndexedItem ("drift",m_tMcNHit, m_tMcDrift);
00631       sc = m_tupleMcHit->addIndexedItem ("x",    m_tMcNHit, m_tMcX);
00632       sc = m_tupleMcHit->addIndexedItem ("y",    m_tMcNHit, m_tMcY);
00633       sc = m_tupleMcHit->addIndexedItem ("z",    m_tMcNHit, m_tMcZ);
00634       sc = m_tupleMcHit->addIndexedItem ("lr",    m_tMcNHit, m_tMcLR);
00635     } else {   
00636       log << MSG::ERROR << "Cannot book tuple MdcMc/mcHit" << endmsg;
00637       //return StatusCode::FAILURE;
00638     }
00639   }
00640   //book tuple of recnstruction results
00641   NTuplePtr nt1(ntupleSvc(), "MdcRec/rec");
00642   if ( nt1 ) { m_tuple1 = nt1;}
00643   else {
00644     m_tuple1 = ntupleSvc()->book ("MdcRec/rec", CLID_ColumnWiseTuple, "MdcTrkRecon results");
00645     if ( m_tuple1 ) {
00646       sc = m_tuple1->addItem ("t0",      m_t0);
00647       sc = m_tuple1->addItem ("t0Stat",  m_t0Stat);
00648       sc = m_tuple1->addItem ("t0Truth", m_t0Truth);
00649 
00650       sc = m_tuple1->addItem ("nTdsTk",  m_nTdsTk);
00651       sc = m_tuple1->addItem ("nMcTk",   m_nMcTk);
00652 
00653       sc = m_tuple1->addItem ("p",       m_p);
00654       sc = m_tuple1->addItem ("pt",      m_pt);
00655       sc = m_tuple1->addItem ("nSlay",   m_nSlay);
00656       sc = m_tuple1->addItem ("pz",      m_pz);
00657       sc = m_tuple1->addItem ("d0",      m_d0);
00658       sc = m_tuple1->addItem ("phi0",    m_phi0);
00659       sc = m_tuple1->addItem ("cpa",     m_cpa);
00660       sc = m_tuple1->addItem ("z0",      m_z0);
00661       sc = m_tuple1->addItem ("tanl",    m_tanl);
00662       sc = m_tuple1->addItem ("q",       m_q);
00663       sc = m_tuple1->addItem ("pocax",   m_pocax);
00664       sc = m_tuple1->addItem ("pocay",   m_pocay);
00665       sc = m_tuple1->addItem ("pocaz",   m_pocaz);
00666 
00667       sc = m_tuple1->addItem ("evtNo",   m_evtNo);
00668       sc = m_tuple1->addItem ("nSt",     m_nSt);
00669       sc = m_tuple1->addItem ("nAct",    m_nAct);
00670       sc = m_tuple1->addItem ("nDof",    m_nDof);
00671       sc = m_tuple1->addItem ("chi2",    m_chi2);
00672       sc = m_tuple1->addItem ("tkId",    m_tkId);
00673       sc = m_tuple1->addItem ("nHit",    m_nHit, 0, 6796);
00674       sc = m_tuple1->addIndexedItem ("resid",    m_nHit, m_resid);
00675       sc = m_tuple1->addIndexedItem ("sigma",    m_nHit, m_sigma);
00676       sc = m_tuple1->addIndexedItem ("driftD",   m_nHit, m_driftD);
00677       sc = m_tuple1->addIndexedItem ("driftT",   m_nHit, m_driftT);
00678       sc = m_tuple1->addIndexedItem ("doca",     m_nHit, m_doca);
00679       sc = m_tuple1->addIndexedItem ("entra",    m_nHit, m_entra);
00680       sc = m_tuple1->addIndexedItem ("ambig",    m_nHit, m_ambig);
00681       sc = m_tuple1->addIndexedItem ("fltLen",   m_nHit, m_fltLen);
00682       sc = m_tuple1->addIndexedItem ("tof",      m_nHit, m_tof);
00683       sc = m_tuple1->addIndexedItem ("act",      m_nHit, m_act);
00684       sc = m_tuple1->addIndexedItem ("tdc",      m_nHit, m_tdc);
00685       sc = m_tuple1->addIndexedItem ("adc",      m_nHit, m_adc);
00686       sc = m_tuple1->addIndexedItem ("layer",    m_nHit, m_layer);
00687       sc = m_tuple1->addIndexedItem ("wire",     m_nHit, m_wire);
00688       sc = m_tuple1->addIndexedItem ("x",        m_nHit, m_x);
00689       sc = m_tuple1->addIndexedItem ("y",        m_nHit, m_y);
00690       sc = m_tuple1->addIndexedItem ("z",        m_nHit, m_z);
00691       sc = m_tuple1->addIndexedItem ("dx",       m_nHit, m_dx);
00692       sc = m_tuple1->addIndexedItem ("dy",       m_nHit, m_dy);
00693       sc = m_tuple1->addIndexedItem ("dz",       m_nHit, m_dz);
00694       sc = m_tuple1->addIndexedItem ("dDriftD",  m_nHit, m_dDriftD);
00695       sc = m_tuple1->addIndexedItem ("dlr",      m_nHit, m_dlr);
00696       sc = m_tuple1->addIndexedItem ("cellWidth",m_nHit, m_cellWidth);//for pickHits tuning
00697     } else {
00698       log << MSG::ERROR << "Cannot book tuple MdcRec/rec" << endmsg;
00699       //return StatusCode::FAILURE;
00700     }
00701   }
00702   //book tuple of segment
00703   NTuplePtr ntSeg(ntupleSvc(), "MdcSeg/seg");
00704   if ( ntSeg ) { m_tupleSeg = ntSeg;}
00705   else {
00706     m_tupleSeg = ntupleSvc()->book ("MdcSeg/seg", CLID_ColumnWiseTuple, "MdcTrkRecon segment data");
00707     if ( m_tupleSeg ) {
00708       sc = m_tupleSeg->addItem ("evtNo",         m_tsEvtNo);
00709       sc = m_tupleSeg->addItem ("nSeg",          m_tsNSeg);
00710       sc = m_tupleSeg->addItem ("nDigi",         m_tsNDigi, 0, 6796);
00711       sc = m_tupleSeg->addIndexedItem ("layer",  m_tsNDigi, m_tsLayer);
00712       sc = m_tupleSeg->addIndexedItem ("wire",   m_tsNDigi, m_tsWire);
00713       sc = m_tupleSeg->addIndexedItem ("inSeg",  m_tsNDigi, m_tsInSeg);
00714       sc = m_tupleSeg->addIndexedItem ("mcTkId", m_tsNDigi, m_tsMcTkId);
00715     } else {   
00716       log << MSG::ERROR << "Cannot book tuple MdcSeg/seg" << endmsg;
00717       //return StatusCode::FAILURE;
00718     }
00719   }
00720 
00721   //book tuple of event
00722   NTuplePtr nt4(ntupleSvc(), "MdcRec/evt");
00723   if ( nt4 ) { m_tupleEvt = nt4;}
00724   else {
00725     m_tupleEvt = ntupleSvc()->book ("MdcRec/evt", CLID_ColumnWiseTuple, "MdcTrkRecon event data");
00726     if ( m_tupleEvt ) {
00727       sc = m_tupleEvt->addItem ("evtNo",         m_t4EvtNo);
00728       sc = m_tupleEvt->addItem ("nMcTk",         m_t4nMcTk );
00729       sc = m_tupleEvt->addItem ("nTdsTk",        m_t4nRecTk );
00730       sc = m_tupleEvt->addItem ("t0",            m_t4t0);
00731       sc = m_tupleEvt->addItem ("t0Stat",        m_t4t0Stat);
00732       sc = m_tupleEvt->addItem ("t0Truth",       m_t4t0Truth);
00733       sc = m_tupleEvt->addItem ("nDigiUn",       m_t4nDigiUnmatch);
00734       sc = m_tupleEvt->addItem ("nDigi",         m_t4nDigi, 0, 6796);
00735       sc = m_tupleEvt->addIndexedItem ("layer",  m_t4nDigi, m_t4Layer);
00736       sc = m_tupleEvt->addIndexedItem ("wire",   m_t4nDigi, m_t4Wire);
00737       sc = m_tupleEvt->addIndexedItem ("rt",     m_t4nDigi, m_t4Time);
00738       sc = m_tupleEvt->addIndexedItem ("rc",     m_t4nDigi, m_t4Charge);
00739       sc = m_tupleEvt->addIndexedItem ("phiMid", m_t4nDigi, m_t4PhiMid);
00740       sc = m_tupleEvt->addIndexedItem ("hot",    m_t4nDigi, m_t4Hot);
00741       //sc = m_tupleEvt->addIndexedItem ("recHit", m_t4nDigi, m_t4recHit);
00742       //sc = m_tupleEvt->addIndexedItem ("rawHit", m_t4nDigi, m_t4rawHit);
00743     } else {
00744       log << MSG::ERROR << "Cannot book N-tuple: MdcRec/evt"   << endmsg;
00745       //return StatusCode::FAILURE;
00746     }
00747   }
00748 
00749   //book tuple of residula for every layer
00750   NTuplePtr ntCombAx(ntupleSvc(), "MdcCombAx/combAx");
00751   if ( ntCombAx ) { g_tupleCombAx= ntCombAx;}
00752   else { 
00753     g_tupleCombAx = ntupleSvc()->book ("MdcCombAx/combAx", CLID_RowWiseTuple, "MdcTrkRecon cuts in combine ax seg");
00754     if ( g_tupleCombAx) {
00755       sc = g_tupleCombAx->addItem ("dPhi0",     g_combAxdPhi0);
00756       sc = g_tupleCombAx->addItem ("dCurv",     g_combAxdCurv);
00757       sc = g_tupleCombAx->addItem ("sigPhi0",   g_combAxSigPhi0);
00758       sc = g_tupleCombAx->addItem ("sigCurv",   g_combAxSigCurv);
00759       sc = g_tupleCombAx->addItem ("slSeed",    g_combAxSlSeed);
00760       sc = g_tupleCombAx->addItem ("slTest",    g_combAxSlTest);
00761       sc = g_tupleCombAx->addItem ("qSeed",     g_combAxQualitySeed);
00762       sc = g_tupleCombAx->addItem ("qTest",     g_combAxQualityTest);
00763       sc = g_tupleCombAx->addItem ("pSeed",     g_combAxPatSeed);
00764       sc = g_tupleCombAx->addItem ("pTest",     g_combAxPatTest);
00765       sc = g_tupleCombAx->addItem ("nHitSeed",  g_combAxNhitSeed);
00766       sc = g_tupleCombAx->addItem ("nHitTest",  g_combAxNhitTest);
00767       sc = g_tupleCombAx->addItem ("mc",        g_combAxMc);
00768       sc = g_tupleCombAx->addItem ("ptTruth",   g_combAxMcPt);
00769       sc = g_tupleCombAx->addItem ("thetaTruth",g_combAxMcTheta);
00770       sc = g_tupleCombAx->addItem ("phiTruth",  g_combAxMcPhi);
00771       sc = g_tupleCombAx->addItem ("ambigSeed", g_combAxMcAmbigSeed);
00772       sc = g_tupleCombAx->addItem ("ambigTest", g_combAxMcAmbigTest);
00773     } else {   
00774       log << MSG::ERROR << "Cannot book N-tuple: FILE101/combAx"   << endmsg;
00775       //return StatusCode::FAILURE;
00776     }
00777   }
00778 
00779   //book tuple of residula for every layer
00780   NTuplePtr ntFindSeg(ntupleSvc(), "MdcSeg/findSeg");
00781   if ( ntFindSeg ) { g_tupleFindSeg = ntFindSeg;}
00782   else { 
00783     g_tupleFindSeg = ntupleSvc()->book ("MdcSeg/findSeg", CLID_RowWiseTuple, "MdcTrkRecon cuts in segment finder");
00784     if ( g_tupleFindSeg) {
00785       sc = g_tupleFindSeg->addItem ("chi2",     g_findSegChi2);
00786       sc = g_tupleFindSeg->addItem ("slope",    g_findSegSlope);
00787       sc = g_tupleFindSeg->addItem ("intercept",g_findSegIntercept);
00788       sc = g_tupleFindSeg->addItem ("chi2Refit",g_findSegChi2Refit);
00789       sc = g_tupleFindSeg->addItem ("chi2Add",  g_findSegChi2Add);
00790       sc = g_tupleFindSeg->addItem ("pat",      g_findSegPat);
00791       sc = g_tupleFindSeg->addItem ("pat34",    g_findSegPat34);
00792       sc = g_tupleFindSeg->addItem ("nhit",     g_findSegNhit);
00793       sc = g_tupleFindSeg->addItem ("slayer",   g_findSegSl);
00794       sc = g_tupleFindSeg->addItem ("mc",       g_findSegMc);
00795       sc = g_tupleFindSeg->addItem ("ambig",    g_findSegAmbig);
00796     } else {   
00797       log << MSG::ERROR << "Cannot book N-tuple: FILE101/findSeg"   << endmsg;
00798       //return StatusCode::FAILURE;
00799     }
00800   }
00801 
00802   //book tuple of event
00803   NTuplePtr ntPick(ntupleSvc(), "MdcTrk/pick");
00804   if ( ntPick ) { m_tuplePick = ntPick;}
00805   else {
00806     m_tuplePick = ntupleSvc()->book ("MdcTrk/pick", CLID_ColumnWiseTuple, "MdcTrkRecon pickHits");
00807     if ( m_tuplePick ) {
00808       sc = m_tuplePick->addItem ("layer",               m_pickLayer);
00809       sc = m_tuplePick->addItem ("nCell",               m_pickNCell, 0, 288);
00810       sc = m_tuplePick->addItem ("nCellWithDigi",       m_pickNCellWithDigi, 0, 288);
00811       sc = m_tuplePick->addIndexedItem ("wire",    m_pickNCellWithDigi, m_pickWire);
00812       sc = m_tuplePick->addIndexedItem ("predDrift",    m_pickNCellWithDigi, m_pickPredDrift);
00813       sc = m_tuplePick->addIndexedItem ("drift",        m_pickNCellWithDigi, m_pickDrift);
00814       sc = m_tuplePick->addIndexedItem ("driftTruth",   m_pickNCellWithDigi, m_pickDriftTruth);
00815       sc = m_tuplePick->addIndexedItem ("phiZOk",       m_pickNCellWithDigi, m_pickPhizOk);
00816       sc = m_tuplePick->addIndexedItem ("z",            m_pickNCellWithDigi, m_pickZ);
00817       sc = m_tuplePick->addIndexedItem ("resid",        m_pickNCellWithDigi, m_pickResid);
00818       sc = m_tuplePick->addIndexedItem ("sigma",        m_pickNCellWithDigi, m_pickSigma);
00819       sc = m_tuplePick->addIndexedItem ("phiLowCut",    m_pickNCellWithDigi, m_pickPhiLowCut);
00820       sc = m_tuplePick->addIndexedItem ("phiHighCut",   m_pickNCellWithDigi, m_pickPhiHighCut);
00821       sc = m_tuplePick->addIndexedItem ("goodDriftCut", m_pickNCellWithDigi, m_pickDriftCut);
00822       sc = m_tuplePick->addIndexedItem ("mcTk",         m_pickNCellWithDigi, m_pickMcTk);
00823       sc = m_tuplePick->addIndexedItem ("is2d",         m_pickNCellWithDigi, m_pickIs2D);
00824       sc = m_tuplePick->addIndexedItem ("pt",           m_pickNCellWithDigi, m_pickPt);
00825       sc = m_tuplePick->addIndexedItem ("curv",           m_pickNCellWithDigi, m_pickCurv);
00826     } else {   
00827       log << MSG::ERROR << "Cannot book N-tuple: MdcTrk/pick"   << endmsg;
00828       //return StatusCode::FAILURE;
00829     }
00830   }
00831 
00832 #ifdef MDCPATREC_TIMETEST
00833   //book tuple of time test
00834   NTuplePtr nt5(ntupleSvc(), "MdcTime/ti");
00835   if ( nt5 ) { m_tupleTime = nt5;}
00836   else {
00837     m_tupleTime = ntupleSvc()->book ("MdcTime/ti", CLID_RowWiseTuple, "MdcTrkRecon time");
00838     if ( m_tupleTime ) {
00839       sc = m_tupleTime->addItem ("eventtime",  m_eventTime);
00840       sc = m_tupleTime->addItem ("recTkNum",   m_t5recTkNum);
00841       sc = m_tupleTime->addItem ("mcTkNum",    m_mcTkNum);
00842       sc = m_tupleTime->addItem ("evtNo",      m_t5EvtNo);
00843       sc = m_tupleTime->addItem ("nHit",       m_t5nHit);
00844       sc = m_tupleTime->addItem ("nDigi",      m_t5nDigi);
00845     } else {
00846       log << MSG::ERROR << "Cannot book N-tuple: FILE101/time"   << endmsg;
00847       //return StatusCode::FAILURE;
00848     }
00849   }
00850   sc = service( "BesTimerSvc", m_timersvc);
00851   if( sc.isFailure() ) {
00852     log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
00853     //return StatusCode::FAILURE;
00854   } 
00855   m_timer[1] = m_timersvc->addItem("Execution");
00856   m_timer[1]->propName("nExecution");
00857 #endif
00858 
00859 #ifdef MDCPATREC_RESLAYER
00860   //book tuple of residula for every layer
00861   NTuplePtr nt6(ntupleSvc(), "MdcRes/res");
00862   if ( nt6 ) { g_tupleres = nt6;}
00863   else {
00864     g_tupleres= ntupleSvc()->book ("MdcRes/res", CLID_RowWiseTuple, "MdcTrkRecon res");
00865     if ( g_tupleres ) {
00866       sc = g_tupleres->addItem ("resid",  g_resLayer);
00867       sc = g_tupleres->addItem ("layer",  g_t6Layer); 
00868     } else {   
00869       log << MSG::ERROR << "Cannot book N-tuple: FILE101/res"   << endmsg;
00870       return StatusCode::FAILURE;
00871     }
00872   }
00873 #endif
00874 
00875   return sc;
00876 }// end of bookNTuple()

void MdcTrkRecon::dumpDigi (  )  [protected]

Definition at line 1294 of file MdcTrkRecon.cxx.

References MdcRawDataProvider::b_dropHot, MdcRawDataProvider::b_keepBadTdc, MdcRawDataProvider::b_keepUnmatch, RawDataProviderSvc::getMdcDigiVec(), iter(), MdcID::layer(), m_dropHot, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_rawDataProviderSvc, RawDataUtil::MdcCharge(), RawDataUtil::MdcTime(), w, and MdcID::wire().

Referenced by execute().

01294                           {
01295   uint32_t getDigiFlag = 0;
01296   getDigiFlag += m_maxMdcDigi;
01297   if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
01298   if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
01299   if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
01300   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
01301   MdcDigiCol::iterator iter = mdcDigiVec.begin();
01302   std::cout<<name()<<" dump MdcDigiVec, nDigi="<<mdcDigiVec.size()<<std::endl;
01303   for (int iDigi=0;iter!= mdcDigiVec.end(); iter++,iDigi++ ) {
01304     long l = MdcID::layer((*iter)->identify());
01305     long w = MdcID::wire((*iter)->identify());
01306     int tkTruth = (*iter)->getTrackIndex();
01307     double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
01308     double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
01309     std::cout<<"("<<l<<","<<w<<";"<<tkTruth<<",t "<<tdc<<",c "<<adc<<")";
01310     if(iDigi%4==0) std::cout<<std::endl;
01311   }
01312   std::cout<<std::endl;
01313 }

void MdcTrkRecon::dumpTdsTrack ( RecMdcTrackCol trackList  )  [protected]

Definition at line 1072 of file MdcTrkRecon.cxx.

References MdcID::layer(), and MdcID::wire().

01072                                                        {
01073   std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
01074   RecMdcTrackCol::iterator it = trackList->begin();
01075   for (;it!= trackList->end();it++){
01076     RecMdcTrack *tk = *it;
01077     std::cout<< endl<<"//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
01078     cout <<" d0 "<<tk->helix(0)
01079       <<" phi0 "<<tk->helix(1)
01080       <<" cpa "<<tk->helix(2)
01081       <<" z0 "<<tk->helix(3)
01082       <<" tanl "<<tk->helix(4)
01083       <<endl;
01084     std::cout<<" q "<<tk->charge() 
01085       <<" theta "<<tk->theta()
01086       <<" phi "<<tk->phi()
01087       <<" x0 "<<tk->x()
01088       <<" y0 "<<tk->y()
01089       <<" z0 "<<tk->z()
01090       <<" r0 "<<tk->r()
01091       <<endl;
01092     std::cout <<" p "<<tk->p()
01093       <<" pt "<<tk->pxy()
01094       <<" px "<<tk->px()
01095       <<" py "<<tk->py()
01096       <<" pz "<<tk->pz()
01097       <<endl;
01098     std::cout<<" tkStat "<<tk->stat()
01099       <<" chi2 "<<tk->chi2()
01100       <<" ndof "<<tk->ndof()
01101       <<" nhit "<<tk->getNhits()
01102       <<" nst "<<tk->nster()
01103       <<endl;
01104     std::cout<< "errmat mat  " << std::endl;
01105     std::cout<< tk->err()<<std::endl;
01106     //std::cout<< "errmat array  " << std::endl;
01107     //for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
01108     //std::cout<< "   " << std::endl;
01109 
01110     int nhits = tk->getVecHits().size();
01111     std::cout<<nhits <<" Hits: " << std::endl;
01112     for(int ii=0; ii <nhits ; ii++){
01113       Identifier id(tk->getVecHits()[ii]->getMdcId());
01114       int layer = MdcID::layer(id);
01115       int wire = MdcID::wire(id);
01116       cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
01117         <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
01118     }//end of hit list
01119     std::cout << "  "<< std::endl;
01120   }//end of tk list
01121   std::cout << "  "<< std::endl;
01122 }//end of dumpTdsTrack

StatusCode MdcTrkRecon::execute (  ) 

Definition at line 215 of file MdcTrkRecon.cxx.

References _gm, abs, dumpDigi(), Bes_Common::FATAL, fillEvent(), fillMcTruth(), fillSegList(), fillTrackList(), MdcFlagHold::findSegs, MdcFlagHold::findTracks, haveDigiAmbig, haveDigiDrift, haveDigiPhi, haveDigiPt, haveDigiTheta, haveDigiTk, hitPoisoned, genRecEmupikp::i, Bes_Common::INFO, iter(), MdcFlagHold::lHist, MdcSegParams::lPrint, m_arbitrateHits, m_bfield, m_combineTracking, m_debug, m_doLineFit, m_fieldCosmic, m_flags, m_hist, m_hitData, m_mcHist, m_mdcPrintSvc, m_recForEsTime, m_segFinder, m_segs, m_selEvtNo, m_tracks, m_tryBunch, m_tupleWireEff, mcDrift, mcLR, mcX, mcY, mcZ, msgSvc(), poisoningHits(), MdcPrintSvc::printRecMdcTrackCol(), recHitMap, EventModel::Recon::RecMdcHitCol, EventModel::Recon::RecMdcTrackCol, MdcFlagHold::segPar, t_eventNo, t_iExexute, t_nDigi, t_nRecTk, t_t0, t_t0Stat, w, and Bes_Common::WARNING.

00215                                 {
00216   setFilterPassed(false); 
00217 
00218   MsgStream log(msgSvc(), name());
00219   log << MSG::INFO << "in execute()" << endreq;
00220 
00221   //Initialize
00222   if(m_hist>0){
00223     for (int ii=0;ii<43;ii++){ 
00224       for (int jj=0;jj<288;jj++){ 
00225         haveDigiTk[ii][jj] = -2; 
00226         haveDigiPt[ii][jj] = -2; 
00227         haveDigiTheta[ii][jj] = -999.; 
00228         haveDigiPhi[ii][jj] = -999.; 
00229         haveDigiDrift[ii][jj] = -999.; 
00230         haveDigiAmbig[ii][jj] = -999.; 
00231         recHitMap[ii][jj] = 0; 
00232         mcDrift[ii][jj] = -99.; 
00233         mcX[ii][jj] = -99.; 
00234         mcY[ii][jj] = -99.; 
00235         mcZ[ii][jj] = -99.; 
00236         mcLR[ii][jj] = -99.; 
00237         hitPoisoned[ii][jj]=-99;
00238       } 
00239     }
00240   }
00241 
00242 
00243   TrkContextEv context(m_bfield);
00244 #ifdef MDCPATREC_TIMETEST
00245   m_timer[1]->start();
00246   ti_nHit=0;
00247 #endif
00248   //------------------read event header--------------
00249   SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00250   if (!evtHead) {
00251     log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00252     return( StatusCode::FAILURE);
00253   }
00254   t_eventNo = evtHead->eventNumber(); 
00255   if(m_debug!=0)std::cout<<t_iExexute<<" run "<<evtHead->runNumber()<<" evt "<<t_eventNo<<std::endl;
00256   t_iExexute++;
00257 
00258   if (m_selEvtNo.size() >0){
00259     std::vector<int>::iterator it = m_selEvtNo.begin();
00260     for (; it!=m_selEvtNo.end(); it++){
00261       if((*it)== t_eventNo) setFilterPassed(true);
00262     }
00263   }
00264   //------------------get event start time-----------
00265   double t0 = 0.;
00266   t_t0 = -1.;
00267   t_t0Stat = -1;
00268   bool iterateT0 = false;
00269   const int nBunch = 3;//3 bunch/event
00270   double iBunch = 0 ; //form 0 to 2
00271   double BeamDeltaTime = 24. / nBunch;// 8ns
00272   SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00273   if (!evTimeCol || evTimeCol->size()==0) { 
00274     log << MSG::WARNING <<" run "<<evtHead->runNumber() <<" evt "<<t_eventNo<< " Could not find RecEsTimeCol" << endreq;
00275     if(m_fieldCosmic){//yzhang temp for bes3 csmc test 
00276       return StatusCode::SUCCESS;
00277     } 
00278   }
00279   RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00280   // skip RecEsTimeCol no rec data
00281   if (iter_evt != evTimeCol->end()){
00282     t_t0Stat = (*iter_evt)->getStat();
00283     t_t0 = (*iter_evt)->getTest();
00284 
00285     if (m_doLineFit){
00286       //t0 -= m_t0Const;
00287       if (abs(t_t0)<0.00001 || (t_t0 < 0) || (t_t0 > 9999.0) ){
00288         log << MSG::WARNING << "Skip with t0 = "<<t_t0 << endreq;
00289         return StatusCode::SUCCESS;//for bes3 cosmic test
00290       }
00291     }
00292     t0 =  t_t0*1.e-9;
00293   }
00294 
00295 
00296 restart:
00297   if ( !m_recForEsTime && m_tryBunch) {
00298     if ( t_t0Stat < 10 ) return StatusCode::SUCCESS;
00299   }
00300   if ( m_tryBunch ){
00301     if ( iBunch > (nBunch - 1) ) return StatusCode::SUCCESS;
00302     if ( t_t0Stat < 0 ){ iterateT0 = true; }
00303     if ( iterateT0 ){
00304       //double EventDeltaTime = 24.;//24ns/event
00305       if ( (t_t0Stat > -1) && (fabs( iBunch * BeamDeltaTime - t_t0) < 2.) ){
00306         ++iBunch;
00307         goto restart;
00308       }else{ t0 = iBunch * BeamDeltaTime *1.e-9;}
00309       ++iBunch;
00310     }
00311   } 
00312 
00313   SmartDataPtr<MdcHitMap> hitMap(eventSvc(),"/Event/Hit/MdcHitMap");
00314   if (!hitMap) {
00315     log << MSG::FATAL << "Could not retrieve MDC hit map" << endreq;
00316     return( StatusCode::FAILURE);
00317   }
00318   //---------------------Read an event--------------
00319   SmartDataPtr<MdcHitCol> hitCol(eventSvc(),"/Event/Hit/MdcHitCol");
00320   if (!hitCol) {
00321     log << MSG::FATAL << "Could not retrieve MDC hit list" << endreq;
00322     return( StatusCode::FAILURE);
00323   }
00324   StatusCode sc;
00325 
00326   //---------- register RecMdcTrackCol and RecMdcHitCol to tds---------
00327   DataObject *aReconEvent;
00328   eventSvc()->findObject("/Event/Recon",aReconEvent);
00329   if(aReconEvent==NULL) {
00330     aReconEvent = new ReconEvent();
00331     StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00332     if(sc!=StatusCode::SUCCESS) {
00333       log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00334       return( StatusCode::FAILURE);
00335     } 
00336   }
00337 
00338   DataObject *aTrackCol;
00339   DataObject *aRecHitCol;
00340   //yzhang 2010-05-28 
00341   SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
00342   if(!m_combineTracking){
00343     eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00344     if(aTrackCol != NULL) {
00345       dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00346       eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00347     }
00348     eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00349     if(aRecHitCol != NULL) {
00350       dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00351       eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00352     }
00353   }
00354 
00355   RecMdcTrackCol* trackList;
00356   eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00357   if (aTrackCol) {
00358     trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00359   }else{
00360     trackList = new RecMdcTrackCol;
00361     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00362     if(!sc.isSuccess()) {
00363       log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00364       return StatusCode::FAILURE;
00365     }
00366   }
00367   RecMdcHitCol* hitList;
00368   eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00369   if (aRecHitCol) {
00370     hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00371   }else{
00372     hitList = new RecMdcHitCol;
00373     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00374     if(!sc.isSuccess()) {
00375       log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00376       return StatusCode::FAILURE;
00377     }
00378   }
00379 
00380   m_hitData->loadevent(hitCol, hitMap, t0);
00381   t_nDigi = hitCol->size();
00382 
00383   if (poisoningHits()) {
00384     m_hitData->poisonHits(_gm,m_debug);
00385   }
00386 
00387   if((m_hist>0) && m_tupleWireEff){
00388     for(int i=0;i<m_hitData->nhits();i++){
00389       const MdcHit* thisHit = m_hitData->hit(i);
00390       int l = thisHit->layernumber();
00391       int w = thisHit->wirenumber();
00392       if(m_hitData->segUsage().get(thisHit)->deadHit()){
00393         hitPoisoned[l][w] = 1;
00394       }else{
00395         hitPoisoned[l][w] = 0;
00396       }
00397     }
00398   }
00399 
00400 #ifdef MDCPATREC_TIMETEST
00401   SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00402   if(!mcParticleCol){
00403     log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
00404   }
00405   m_mcTkNum = mcParticleCol->size();
00406 #endif
00407 
00408   
00409   if(m_debug>1) dumpDigi();
00410   //--------------------------------------------------------------------------
00411   // Segment finding
00412   //--------------------------------------------------------------------------
00413   int nSegs = 0;
00414   if (m_flags.findSegs) {
00415     // Form track segments in superlayers
00416     nSegs = m_segFinder->createSegs(_gm, *m_segs, m_hitData->segUsage(), 
00417         m_hitData->hitMap(), t0);
00418   }
00419   log << MSG::INFO << " Found " << nSegs << " segments" << endreq;
00420   if (m_debug>1){
00421     std::cout<<"------------------------------------------------"<<std::endl;
00422     std::cout<<"==event "<<t_eventNo<< " Found " << nSegs << " segments" << std::endl;
00423     m_segs->plot();
00424   }
00425 
00426   if (m_flags.lHist||m_flags.segPar.lPrint) {
00427     fillSegList();
00428   } 
00429 
00430   //--------------------------------------------------------------------------
00431   // Combine segments to form track
00432   //-------------------------------------------------------------------------- 
00433   int nTracks  = 0;
00434   int nDeleted = 0;
00435   if (m_flags.findTracks && m_flags.findSegs) {
00436     if (nSegs > 2) {
00437       nTracks = m_tracks->createFromSegs(m_segs.get(), m_hitData->hitMap(),
00438           _gm, context, t0); 
00439     }
00440 
00441     if(m_arbitrateHits) nDeleted = m_tracks->arbitrateHits();
00442     int nKeep = nTracks - nDeleted;
00443 
00444     if (m_debug>0 && (nTracks - nDeleted)>0){
00445       std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
00446         <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
00447         << " keep  "<< nKeep <<" track(s)"; 
00448       if (nDeleted!=0){ std::cout <<",deleted " <<nDeleted; }
00449       std::cout<< std::endl;//yzhang debug
00450     }else{
00451       if (m_debug>0){
00452         std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
00453           <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
00454           <<" found 0 track "<< endl;
00455       }
00456     }
00457 
00458     if (m_flags.lHist) t_nRecTk = nKeep;
00459 
00460 #ifdef MDCPATREC_RESLAYER
00461     m_tracks->setResLayer(m_resLayer);
00462 #endif 
00463     if (m_flags.lHist){ fillTrackList(); }
00464     // Stick the found tracks into the list in TDS
00465     m_tracks->store(trackList, hitList);
00466     //if (m_debug>1) { dumpTdsTrack(trackList); }
00467     if(m_debug>1){
00468       std::cout<<name()<<" print track list  "<<std::endl;
00469       m_mdcPrintSvc->printRecMdcTrackCol();
00470     }
00471 
00472     //if(nKeep != (int)trackList->size()) std::cout<<__FILE__<<" "<<__LINE__<<"!!!!!!!!!!!!!!!!! nKeep != nTdsTk"<<  std::endl;
00473 #ifdef MDCPATREC_TIMETEST
00474     RecMdcTrackCol::iterator iter = trackList->begin();
00475     for (;iter != trackList->end(); iter++) {
00476       MdcTrack* tk = *iter;
00477       ti_nHit += tk->getNhits(); 
00478     }
00479 #endif
00480 
00481   } 
00482   //-------------End track-finding------------------ 
00483   m_segs->destroySegs();
00484 
00485   //if ( t_eventNo% 1000 == 0) { 
00486   //std::cout<<"evt "<<t_eventNo<< " Found " << trackList->size() << " track(s)" <<endl;//yzhang debug
00487   //}
00488   if(m_flags.lHist && m_mcHist) fillMcTruth();
00489 #ifdef MDCPATREC_TIMETEST
00490   m_timer[1]->stop(); 
00491   //cout << "m_timer[1]->elapsed()::" << m_timer[1]->elapsed() << endl;
00492   //cout << "m_timer[1]->mean()::" << m_timer[1]->mean() << endl;
00493   m_eventTime = m_timer[1]->elapsed();
00494   m_t5recTkNum = m_tracks->length();
00495   m_t5EvtNo = t_eventNo;
00496   m_t5nHit = ti_nHit;
00497   m_t5nDigi = t_nDigi;
00498   sc = m_tupleTime->write();
00499 #endif
00500   // for event start time, if track not found try another t0
00501   if ( m_tryBunch ){
00502     if ( nTracks <1 ){ iterateT0 = true; 
00503       goto restart;
00504     }
00505   }
00506   if (m_flags.lHist ) fillEvent();
00507 
00508   return StatusCode::SUCCESS;
00509 }

void MdcTrkRecon::fillEvent (  )  [protected]

Definition at line 1251 of file MdcTrkRecon.cxx.

References _gm, MdcRawDataProvider::b_dropHot, MdcRawDataProvider::b_keepBadTdc, MdcRawDataProvider::b_keepUnmatch, RawDataProviderSvc::getMdcDigiVec(), iter(), MdcDetector::Layer(), MdcID::layer(), m_dropHot, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_minMdcDigi, m_rawDataProviderSvc, m_t4Charge, m_t4EvtNo, m_t4Hot, m_t4Layer, m_t4nDigi, m_t4nDigiUnmatch, m_t4nMcTk, m_t4nRecTk, m_t4PhiMid, m_t4t0, m_t4t0Stat, m_t4t0Truth, m_t4Time, m_t4Wire, m_tupleEvt, RawDataUtil::MdcCharge(), RawDataUtil::MdcTime(), MdcLayer::phiWire(), recHitMap, t_eventNo, t_mcTkNum, t_nRecTk, t_t0, t_t0Stat, t_t0Truth, w, and MdcID::wire().

Referenced by execute().

01251                            {
01252   if (m_tupleEvt!=NULL){
01253 
01254     uint32_t getDigiFlag = 0;
01255     getDigiFlag += m_maxMdcDigi;
01256     if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
01257     if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
01258     if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
01259     MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
01260     if ((int)mdcDigiVec.size()<m_minMdcDigi){
01261       std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
01262     }
01263 
01264     m_t4nDigi = mdcDigiVec.size();
01265     int iDigi=0;
01266     MdcDigiCol::iterator iter = mdcDigiVec.begin();
01267     for (;iDigi<m_t4nDigi; iter++ ) {
01268       double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
01269       double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
01270       m_t4Time[iDigi] = tdc;
01271       m_t4Charge[iDigi] = adc;
01272       long l = MdcID::layer((*iter)->identify());
01273       long w = MdcID::wire((*iter)->identify());
01274       m_t4Layer[iDigi] = l;
01275       m_t4Wire[iDigi] = w;
01276       m_t4PhiMid[iDigi] = _gm->Layer(l)->phiWire(w);
01277       m_t4Hot[iDigi] = recHitMap[l][w];
01278       iDigi++;
01279     }
01280     m_t4t0 = t_t0;
01281     m_t4t0Stat = t_t0Stat;
01282     m_t4t0Truth = t_t0Truth;
01283     m_t4EvtNo =  t_eventNo;
01284     m_t4nRecTk = t_nRecTk;
01285     SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
01286     if (mdcDigiCol){
01287       m_t4nDigiUnmatch = mdcDigiCol->size();
01288     }
01289     m_t4nMcTk= t_mcTkNum;
01290     m_tupleEvt->write();
01291   }
01292 }//end of fillEvent 

void MdcTrkRecon::fillMcTruth (  )  [protected]

Definition at line 878 of file MdcTrkRecon.cxx.

References abs, MdcRawDataProvider::b_dropHot, MdcRawDataProvider::b_keepBadTdc, MdcRawDataProvider::b_keepUnmatch, RawDataProviderSvc::getMdcDigiVec(), haveDigiAmbig, haveDigiDrift, haveDigiPhi, haveDigiPt, haveDigiTheta, haveDigiTk, hitInSegList, hitPoisoned, genRecEmupikp::i, Bes_Common::INFO, isPrimaryOfMcTk, iter(), MdcID::layer(), m_dropHot, m_keepBadTdc, m_keepUnmatch, m_maxMdcDigi, m_mcHist, m_minMdcDigi, m_rawDataProviderSvc, m_tMcD0, m_tMcEvtNo, m_tMcFiTerm, m_tMcNHit, m_tMcNTk, m_tMcPid, m_tMcPx, m_tMcPy, m_tMcPz, m_tMcTheta, m_tMcTkId, m_tMcZ0, m_tupleMc, m_tupleMcHit, m_tupleWireEff, m_we_gwire, m_we_layer, m_we_pdg, m_we_phi, m_we_poisoned, m_we_primary, m_we_pt, m_we_seg, m_we_theta, m_we_tkId, m_we_track, m_we_wire, mcDrift, mcLR, mcX, mcY, mcZ, msgSvc(), Constants::nWireBeforeLayer, pdgOfMcTk, recHitMap, t_eventNo, t_mcTkNum, t_nHitInTk, t_t0Truth, w, and MdcID::wire().

Referenced by execute().

00878                              {
00879   MsgStream log(msgSvc(), name());
00880   StatusCode sc;
00881   t_mcTkNum = 0;
00882   for(int i=0;i<100;i++){
00883     isPrimaryOfMcTk[i]=-9999;
00884     pdgOfMcTk[i]=-9999;
00885   }
00886   double ptOfMcTk[100]={0.};
00887   double thetaOfMcTk[100]={0.};
00888   double phiOfMcTk[100]={0.};
00889   double nDigiOfMcTk[100]={0.};
00890   if (m_mcHist){
00891     //------------------Retrieve MC truth MdcMcHit------------
00892     SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol"); 
00893     if (!mcMdcMcHitCol) {
00894       log << MSG::INFO << "Could not find MdcMcHit" << endreq; 
00895     }else{
00896       Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
00897       for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
00898         const Identifier id= (*iter_mchit)->identify();
00899         int layer = MdcID::layer(id);
00900         int wire = MdcID::wire(id);
00901         int iMcTk =  (*iter_mchit)->getTrackIndex();
00902         if(iMcTk<100&&iMcTk>0) nDigiOfMcTk[iMcTk]++;
00903         mcDrift[layer][wire]= (*iter_mchit)->getDriftDistance()/10.;  //drift in MC.
00904         //std::cout<<"   "<<__FILE__<<" mcDrift  "<<mcDrift[layer][wire]<<"   "<<std::endl;
00905         mcX[layer][wire] = (*iter_mchit)->getPositionX()/10.;
00906         mcY[layer][wire] = (*iter_mchit)->getPositionY()/10.;
00907         mcZ[layer][wire] = (*iter_mchit)->getPositionZ()/10.;
00908         mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
00909         if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
00910         t_nHitInTk[iMcTk]++; 
00911         haveDigiAmbig[layer][wire] =  mcLR[layer][wire];
00912         haveDigiDrift[layer][wire] = mcDrift[layer][wire];
00913       }
00914     }
00915 
00916     //------------------get event start time truth-----------
00917     t_t0Truth = -10.;
00918     SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00919     if(!mcParticleCol){
00920       log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
00921     }else {
00922       int nMcTk = 0;
00923       Event::McParticleCol::iterator it = mcParticleCol->begin();
00924       for (;it != mcParticleCol->end(); it++){ 
00925         //int pdg_code = (*it)->particleProperty();
00926         //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
00927         t_mcTkNum++;
00928         nMcTk++;
00929       }
00930       it = mcParticleCol->begin();
00931       for (;it != mcParticleCol->end(); it++){
00932         int tkId   = (*it)->trackIndex();
00933         if ((*it)->primaryParticle()){
00934           t_t0Truth = (*it)->initialPosition().t();
00935           isPrimaryOfMcTk[tkId] = 1;
00936         }else{
00937           isPrimaryOfMcTk[tkId] = 0;
00938           //continue;
00939         }
00940         //fill charged particle
00941         int pdg_code = (*it)->particleProperty();
00942         pdgOfMcTk[tkId] = pdg_code;
00943         //std::cout<<" tkId  "<<tkId<<" pdg_code  "<<pdg_code<<"   "<<std::endl;
00944         //FIXME skip charged track and track outside MDC
00945         //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
00946         const CLHEP::Hep3Vector& true_mom = (*it)->initialFourMomentum().vect();
00947         const CLHEP::HepLorentzVector& posIni = (*it)->initialPosition();
00948         const CLHEP::HepLorentzVector& posFin = (*it)->finalPosition();
00949         if(tkId<100&&tkId>=0) {
00950           ptOfMcTk[tkId] = true_mom.perp();
00951           thetaOfMcTk[tkId] = (*it)->initialFourMomentum().theta();
00952           phiOfMcTk[tkId] = (*it)->initialFourMomentum().phi();
00953         }
00954 
00955         if ( m_tupleMcHit){
00956           m_tMcPx    = true_mom.x();
00957           m_tMcPy    = true_mom.y();
00958           m_tMcPz    = true_mom.z();
00959           m_tMcD0    = posIni.perp()/10.;
00960           m_tMcZ0    = posIni.z()/10.;
00961           m_tMcTheta = (*it)->initialFourMomentum().theta();
00962           m_tMcFiTerm= posFin.phi();
00963           m_tMcPid   = pdg_code;
00964           m_tMcTkId  = tkId;
00965           m_tMcNHit  = t_nHitInTk[tkId];
00966           m_tupleMcHit->write();
00967         }
00968       }//end of loop mcParticleCol
00969       if(m_tupleMc) {
00970         m_tMcNTk= nMcTk;
00971         m_tMcEvtNo = t_eventNo;
00972         m_tupleMc->write();
00973       }
00974     }
00975   }//end of m_mcHist
00976 
00977   uint32_t getDigiFlag = 0;
00978   getDigiFlag += m_maxMdcDigi;
00979   if(m_dropHot)     getDigiFlag |= MdcRawDataProvider::b_dropHot;
00980   if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00981   if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00982   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00983   if ((int)mdcDigiVec.size()<m_minMdcDigi){
00984     std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
00985   }
00986 
00987   if (mdcDigiVec.size()<=0) return;
00988   //fill raw digi and t0 status
00989   MdcDigiCol::iterator iter = mdcDigiVec.begin();
00990   for (;iter!= mdcDigiVec.end(); iter++ ) {
00991     long l = MdcID::layer((*iter)->identify());
00992     long w = MdcID::wire((*iter)->identify());
00993     int tkId = (*iter)->getTrackIndex();
00994     haveDigiTk[l][w]= tkId;
00995     //std::cout<<" l  "<<l<<" w  "<<w<<"  tk "<<tkId<<std::endl;
00996     haveDigiPt[l][w] = ptOfMcTk[(*iter)->getTrackIndex()];
00997     haveDigiTheta[l][w] = thetaOfMcTk[(*iter)->getTrackIndex()];
00998     haveDigiPhi[l][w] = phiOfMcTk[(*iter)->getTrackIndex()];
00999     if(m_tupleWireEff){
01000       m_we_tkId = tkId;
01001       if(tkId>=0){
01002         if(tkId>=1000) tkId-=1000;
01003         m_we_pdg = pdgOfMcTk[tkId];
01004         m_we_primary = isPrimaryOfMcTk[tkId];
01005         //if(pdgOfMcTk[tkId]==-22) cout<<" primaryParticle ? "<< isPrimaryOfMcTk[tkId]<< " tkId "<<tkId <<" layer "<<l<<" wire "<<w<<" pdg="<<pdgOfMcTk[tkId]<<endl;
01006       }else{
01007         m_we_pdg = -999;
01008         m_we_primary = -999;
01009       }
01010       m_we_layer = l;
01011       m_we_wire = w;
01012       int gwire = w+Constants::nWireBeforeLayer[l];
01013       m_we_gwire = gwire;
01014       m_we_poisoned = hitPoisoned[l][w];
01015       if(hitInSegList[l][w]>0) m_we_seg = 1;
01016       else m_we_seg = 0;
01017       if(recHitMap[l][w]>0) m_we_track = 1;
01018       else m_we_track = 0;
01019       if(l==35&&tkId>=0&&abs(pdgOfMcTk[tkId])==11&&hitInSegList[l][w]<=0) {
01020         cout<<"layer "<<l<<" cell "<<w<<" gwire "<<gwire<<" inseg "<<hitInSegList[l][w]<<endl;
01021       }
01022       m_we_pt = ptOfMcTk[tkId];
01023       m_we_theta = thetaOfMcTk[tkId];
01024       m_we_phi = phiOfMcTk[tkId];
01025       m_tupleWireEff->write();
01026     }
01027   }
01028 }//end of fillMcTruth()

void MdcTrkRecon::fillSegList (  )  [protected]

Definition at line 1030 of file MdcTrkRecon.cxx.

References _gm, haveDigiTk, hitInSegList, genRecEmupikp::i, MdcFlagHold::lHist, MdcSegParams::lPrint, m_flags, m_segs, m_tsEvtNo, m_tsInSeg, m_tsLayer, m_tsMcTkId, m_tsNDigi, m_tsNSeg, m_tsWire, m_tupleSeg, MdcDetector::nSuper(), MdcFlagHold::segPar, t_eventNo, and t_nDigi.

Referenced by execute().

01030                               {
01031   if (2 == m_flags.segPar.lPrint) {
01032     std::cout << "print after Segment finding " << std::endl;
01033   }
01034   if (!m_flags.lHist) return;
01035   // Fill hits of every layer after segment finding
01036   for (int ii=0;ii<43;ii++){
01037     for (int jj=0;jj<288;jj++){ hitInSegList[ii][jj] = 0; }
01038   }
01039   int nSeg=0;
01040   for (int i = 0; i < _gm->nSuper(); i++) {
01041     MdcSeg* aSeg = (MdcSeg *) m_segs->oneList(i)->first();
01042     while (aSeg != NULL) {
01043       nSeg++;
01044       for (int ihit=0 ; ihit< aSeg->nHit() ; ihit++){
01045         const MdcHit* hit = aSeg->hit(ihit)->mdcHit();
01046         int layer = hit->layernumber();
01047         int wire = hit->wirenumber();
01048         hitInSegList[layer][wire]++;
01049       }
01050       aSeg = (MdcSeg *) aSeg->next();
01051     }
01052   }//end for slayer
01053   if (!m_tupleSeg) return;
01054   m_tsEvtNo =  t_eventNo;
01055   m_tsNDigi = t_nDigi;
01056   int iDigi=0;
01057   for (int ii=0;ii<43;ii++){
01058     for (int jj=0;jj<288;jj++){ 
01059       if (haveDigiTk[ii][jj] > -2){
01060         m_tsLayer[iDigi] = ii;
01061         m_tsWire[iDigi] = jj;
01062         m_tsInSeg[iDigi] = hitInSegList[ii][jj];
01063         m_tsMcTkId[iDigi] = haveDigiTk[ii][jj];
01064         iDigi++;
01065       }
01066     }
01067   }
01068   m_tsNSeg=nSeg;
01069   m_tupleSeg->write();
01070 }//end of fillSegList

void MdcTrkRecon::fillTrackList (  )  [protected]

Definition at line 1124 of file MdcTrkRecon.cxx.

References _gm, TrkHitList::begin(), TrkAbsFit::charge(), TrkAbsFit::chisq(), TrkExchangePar::d0(), MdcFlagHold::debugFlag(), TrkHitList::end(), TrkFit::helix(), genRecEmupikp::i, MdcDetector::Layer(), MdcFlagHold::lHist, m_act, m_adc, m_ambig, m_cellWidth, m_chi2, m_cpa, m_d0, m_dDriftD, m_dlr, m_doca, m_driftD, m_driftT, m_dx, m_dy, m_dz, m_entra, m_evtNo, m_flags, m_fltLen, m_layer, m_nAct, m_nDof, m_nHit, m_nMcTk, m_nSlay, m_nSt, m_nTdsTk, m_p, m_phi0, m_pocax, m_pocay, m_pocaz, m_pt, m_pz, m_q, m_resid, m_sigma, m_t0, m_t0Stat, m_t0Truth, m_tanl, m_tdc, m_tkId, m_tof, m_tracks, m_tuple1, m_wire, m_x, m_y, m_z, m_z0, mcDrift, mcLR, mcX, mcY, mcZ, TrkAbsFit::momentum(), TrkFit::nActive(), TrkAbsFit::nDof(), TrkHitList::nHit(), MdcLayer::nWires(), TrkExchangePar::phi0(), boss::pos, TrkAbsFit::position(), TrkAbsFit::pt(), recHitMap, MdcLayer::rMid(), t_eventNo, t_mcTkNum, t_nRecTk, t_t0, t_t0Stat, t_t0Truth, TrkExchangePar::tanDip(), Constants::twoPi, and TrkExchangePar::z0().

Referenced by execute().

01124                                { 
01125   if (m_flags.debugFlag()>1) {
01126     std::cout << "=======Print after Track finding=======" << std::endl;
01127     m_tracks->plot();
01128   }
01129 
01130   if (!m_flags.lHist) return; 
01131   //----------- fill hit -----------
01132   int nTk = (*m_tracks).nTrack();
01133   for (int itrack = 0; itrack < nTk; itrack++) {
01134     MdcTrack *atrack = (*m_tracks)[itrack];
01135     if (atrack==NULL) continue;
01136     const TrkFit* fitResult = atrack->track().fitResult();
01137     if (fitResult == 0) continue;//check the fit worked
01138 
01139     //for fill ntuples
01140     int nSt=0; //int nSeg=0; 
01141     int seg[11] = {0}; int segme;
01142     //-----------hit list-------------
01143     const TrkHitList* hitlist = atrack->track().hits();
01144     TrkHitList::hot_iterator hot = hitlist->begin();
01145     int layerCount[43]={0};
01146     for(int iii=0;iii<43;iii++){layerCount[iii]=0;}
01147     int i=0;
01148     for (;hot!= hitlist->end();hot++){
01149       const MdcRecoHitOnTrack* recoHot;                   
01150       recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
01151       int layer = recoHot->mdcHit()->layernumber();
01152       int wire  = recoHot->mdcHit()->wirenumber();
01153       layerCount[layer]++;
01154       recHitMap[layer][wire]++;
01155       //for n seg
01156       segme=0;
01157       if ( layer >0 ) {segme=(layer-1)/4;}
01158       seg[segme]++;
01159       if (recoHot->layer()->view()) { ++nSt; }
01160 
01161       if(!m_tuple1) continue;
01162       m_layer[i]    = layer;
01163       m_wire[i]     = wire;
01164       m_ambig[i]    = recoHot->wireAmbig();// wire ambig
01165       //fltLen
01166       double fltLen = recoHot->fltLen(); 
01167       m_fltLen[i]   = fltLen;
01168       double tof    = recoHot->getParentRep()->arrivalTime(fltLen);
01169       //position
01170       HepPoint3D pos; Hep3Vector dir;    
01171       recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
01172       m_x[i] = pos.x();      
01173       m_y[i] = pos.y();      
01174       m_z[i] = pos.z();      
01175       m_driftT[i]   = recoHot->mdcHit()->driftTime(tof,pos.z());
01176       m_tof[i]      = tof;
01177       m_driftD[i]   = recoHot->drift();
01178       m_sigma[i]    = recoHot->hitRms();            
01179       m_tdc[i]      = recoHot->rawTime();
01180       m_adc[i]      = recoHot->mdcHit()->charge(); 
01181       m_doca[i]     = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME 
01182       m_entra[i]    = recoHot->entranceAngle();
01183       m_act[i]      = recoHot->isActive();
01184       //diff with truth
01185       m_dx[i] = m_x[i] - mcX[layer][wire];
01186       m_dy[i] = m_y[i] - mcY[layer][wire];
01187       m_dz[i] = m_z[i] - mcZ[layer][wire];
01188       m_dDriftD[i] = m_driftD[i] - mcDrift[layer][wire];
01189       m_dlr[i] = m_ambig[i] - mcLR[layer][wire];
01190       //yzhang for pickHits debug
01191       m_cellWidth[i] = Constants::twoPi*_gm->Layer(layer)->rMid()/_gm->Layer(layer)->nWires();
01192       //zhangy
01193       //resid
01194       double res=999.,rese=999.;
01195       if (recoHot->resid(res,rese,false)){
01196       }else{}
01197       m_resid[i] = res;
01198       i++;
01199     }// end fill hit
01200     int nSlay=0;
01201     for(int i=0;i<11;i++){
01202       if (seg[i]>0) nSlay++;
01203     }
01204     if(m_tuple1){
01205       //------------fill track result-------------
01206       m_tkId = itrack;
01207       //track parameters at closest approach to beamline
01208       double fltLenPoca = 0.0;
01209       TrkExchangePar helix = fitResult->helix(fltLenPoca);
01210       m_q    = fitResult->charge();
01211       m_phi0 = helix.phi0();
01212       m_tanl = helix.tanDip();
01213       m_z0 = helix.z0();
01214       m_d0 = helix.d0();
01215       m_pt = fitResult->pt();
01216       m_nSlay = nSlay;
01217       if (m_pt > 0.00001){
01218         m_cpa = fitResult->charge()/fitResult->pt();
01219       }
01220       //momenta and position
01221       CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
01222       double px,py,pz,pxy;
01223       pxy = fitResult->pt();
01224       px = p1.x();
01225       py = p1.y();
01226       pz = p1.z();
01227       m_p = p1.mag();
01228       m_pz = pz;
01229       HepPoint3D poca = fitResult->position(fltLenPoca);
01230       m_pocax = poca.x();
01231       m_pocay = poca.y();
01232       m_pocaz = poca.z(); 
01233       m_nAct = fitResult->nActive();
01234       m_nHit = hitlist->nHit();
01235       m_nDof = fitResult->nDof();
01236       m_nSt  = nSt;
01237       m_chi2 = fitResult->chisq();
01238       //for (int l=0;l<43;l++) m_layerCount[l] = layerCount[l];
01239       m_t0 = t_t0;
01240       m_t0Stat = t_t0Stat;
01241       m_t0Truth = t_t0Truth;
01242       m_nTdsTk = t_nRecTk;
01243       m_nMcTk = t_mcTkNum;
01244       m_evtNo = t_eventNo;
01245       m_tuple1->write();
01246     }
01247   }//end of loop rec trk list
01248 }//end of fillTrackList

StatusCode MdcTrkRecon::finalize (  ) 

Definition at line 512 of file MdcTrkRecon.cxx.

References Bes_Common::INFO, m_tracks, and msgSvc().

00512                                  {
00513   MsgStream log(msgSvc(), name());
00514   log << MSG::INFO << "in finalize()" << endreq;
00515 
00516   m_tracks.reset(0);
00517 #ifdef MDCPATREC_TIMETEST
00518   m_timersvc->print();
00519 #endif
00520   return StatusCode::SUCCESS;
00521 }

bool MdcTrkRecon::ignoringUsedHits (  )  const [inline, protected]

Definition at line 41 of file MdcTrkRecon.h.

References m_onlyUnusedHits.

Referenced by initialize().

00041 {return m_onlyUnusedHits;}

StatusCode MdcTrkRecon::initialize (  ) 

Definition at line 138 of file MdcTrkRecon.cxx.

References BField::bFieldNominal(), bookNTuple(), MdcFlagHold::debugFlag(), calibUtil::ERROR, Bes_Common::FATAL, MdcFlagHold::findSegs, genRecEmupikp::i, ignoringUsedHits(), Bes_Common::INFO, MdcFlagHold::lHist, m_bfield, m_d0Cut, MdcTrackListBase::m_d0Cut, TrkHelixFitter::m_debug, m_debug, m_doLineFit, m_dropTrkPt, m_flags, m_helixHitsSigma, m_hist, m_hitData, m_mdcPrintSvc, m_paramFile, m_pdtFile, m_pIMF, MdcTrackListBase::m_ptCut, m_rawDataProviderSvc, m_segFinder, m_tracks, m_z0Cut, MdcTrackListBase::m_z0Cut, msgSvc(), TrkHelixFitter::nSigmaCut, Pdt::readMCppTable(), MdcFlagHold::readPar(), MdcFlagHold::segPar, MdcFlagHold::setDebug(), t_iExexute, MdcFlagHold::tkParTight, and MdcSegParams::useAllAmbig.

00138                                   {
00139 
00140   MsgStream log(msgSvc(), name());
00141   log << MSG::INFO << "in initialize()" << endreq;
00142   StatusCode sc;
00143 
00144   //Pdt
00145   Pdt::readMCppTable(m_pdtFile);
00146 
00147   //Flag and Pars
00148   m_flags.readPar(m_paramFile);
00149   m_flags.lHist = m_hist;
00150   m_flags.setDebug(m_debug);
00151   for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
00152   if (!m_doLineFit) MdcTrackListBase::m_d0Cut = m_d0Cut;
00153   if (!m_doLineFit) MdcTrackListBase::m_z0Cut = m_z0Cut;
00154   MdcTrackListBase::m_ptCut = m_dropTrkPt;
00155   if (m_debug>0) {
00156     std::cout<< "MdcTrkRecon d0Cut "<<m_d0Cut<< " z0Cut "<< 
00157       m_z0Cut<<" ptCut "<<m_dropTrkPt << std::endl;
00158   }
00159 
00160   //Initailize magnetic filed 
00161   sc = service ("MagneticFieldSvc",m_pIMF);
00162   if(sc!=StatusCode::SUCCESS) {
00163     log << MSG::ERROR << name() << " Unable to open magnetic field service"<<endreq;
00164   }
00165   m_bfield = new BField(m_pIMF);
00166   log << MSG::INFO <<name() << " field z = "<<m_bfield->bFieldNominal()<< endreq;
00167 
00168   //Initialize RawDataProviderSvc
00169   IRawDataProviderSvc* irawDataProviderSvc;
00170   sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00171   m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00172   if ( sc.isFailure() ){
00173     log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
00174     return StatusCode::FAILURE;
00175   } 
00176 
00177   //Initailize MdcPrintSvc
00178   IMdcPrintSvc* imdcPrintSvc;
00179   sc = service ("MdcPrintSvc", imdcPrintSvc);
00180   m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
00181   if ( sc.isFailure() ){
00182     log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
00183     return StatusCode::FAILURE;
00184   }
00185 
00186   //Initialize hitdata
00187   m_hitData.reset( new MdcSegData( ignoringUsedHits() ));
00188 
00189   //Initialize segFinder
00190   if (m_flags.findSegs) {
00191     m_segFinder.reset( new MdcSegFinder(m_flags.segPar.useAllAmbig) );
00192   }
00193 
00194   //Initialize trackList
00195   if ( m_doLineFit ){
00196     m_tracks.reset( new MdcTrackListCsmc(m_flags.tkParTight) );
00197   } else {
00198     m_tracks.reset( new MdcTrackList(m_flags.tkParTight) );
00199   }
00200 
00201   //Book NTuple and Histograms
00202   if (m_flags.lHist){ 
00203     StatusCode sc = bookNTuple();
00204     if (!sc.isSuccess()) { m_flags.lHist = 0;}
00205   }
00206 
00207   //yzhang for HelixFitter debug
00208   if(7 == m_flags.debugFlag())TrkHelixFitter::m_debug = true;
00209 
00210   t_iExexute = 0;
00211   return StatusCode::SUCCESS;
00212 }

bool MdcTrkRecon::poisoningHits (  )  const [inline, protected]

Definition at line 42 of file MdcTrkRecon.h.

References m_poisonHits.

Referenced by execute().

00042 {return m_poisonHits;}


Member Data Documentation

const MdcDetector* MdcTrkRecon::_gm [private]

Definition at line 57 of file MdcTrkRecon.h.

Referenced by beginRun(), execute(), fillEvent(), fillSegList(), and fillTrackList().

int MdcTrkRecon::hitInSegList[43][288] [private]

Definition at line 104 of file MdcTrkRecon.h.

Referenced by fillMcTruth(), and fillSegList().

double MdcTrkRecon::hitOnMcTk[43] [private]

Definition at line 102 of file MdcTrkRecon.h.

double MdcTrkRecon::hitOnSeg[43] [private]

Definition at line 101 of file MdcTrkRecon.h.

int MdcTrkRecon::hitPoisoned[43][288] [private]

Definition at line 103 of file MdcTrkRecon.h.

Referenced by execute(), and fillMcTruth().

int MdcTrkRecon::isPrimaryOfMcTk[100] [private]

Definition at line 114 of file MdcTrkRecon.h.

Referenced by fillMcTruth().

int MdcTrkRecon::m_allHit [private]

Definition at line 53 of file MdcTrkRecon.h.

Referenced by MdcTrkRecon().

bool MdcTrkRecon::m_arbitrateHits [private]

Definition at line 91 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

BField* MdcTrkRecon::m_bfield [private]

Definition at line 66 of file MdcTrkRecon.h.

Referenced by execute(), and initialize().

bool MdcTrkRecon::m_combineTracking [private]

Definition at line 92 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

std::string MdcTrkRecon::m_configFile [private]

Definition at line 54 of file MdcTrkRecon.h.

Referenced by MdcTrkRecon().

double MdcTrkRecon::m_d0Cut [private]

Definition at line 84 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

int MdcTrkRecon::m_debug [private]

Definition at line 89 of file MdcTrkRecon.h.

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

bool MdcTrkRecon::m_doLineFit [private]

Definition at line 71 of file MdcTrkRecon.h.

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

bool MdcTrkRecon::m_doSagFlag [private]

Definition at line 90 of file MdcTrkRecon.h.

Referenced by beginRun(), and MdcTrkRecon().

bool MdcTrkRecon::m_dropHot [private]

Definition at line 79 of file MdcTrkRecon.h.

Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and MdcTrkRecon().

double MdcTrkRecon::m_dropTrkPt [private]

Definition at line 86 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

bool MdcTrkRecon::m_fieldCosmic [private]

Definition at line 83 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

int MdcTrkRecon::m_fittingMethod [private]

Definition at line 52 of file MdcTrkRecon.h.

Referenced by MdcTrkRecon().

MdcFlagHold MdcTrkRecon::m_flags [private]

Definition at line 58 of file MdcTrkRecon.h.

Referenced by beginRun(), execute(), fillSegList(), fillTrackList(), and initialize().

uint32_t MdcTrkRecon::m_getDigiFlag [private]

Definition at line 76 of file MdcTrkRecon.h.

Referenced by MdcTrkRecon().

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

Definition at line 74 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

int MdcTrkRecon::m_hist [private]

Definition at line 88 of file MdcTrkRecon.h.

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

std::auto_ptr<MdcSegData> MdcTrkRecon::m_hitData [private]

Definition at line 59 of file MdcTrkRecon.h.

Referenced by execute(), initialize(), and ~MdcTrkRecon().

bool MdcTrkRecon::m_keepBadTdc [private]

Definition at line 78 of file MdcTrkRecon.h.

Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and MdcTrkRecon().

bool MdcTrkRecon::m_keepUnmatch [private]

Definition at line 80 of file MdcTrkRecon.h.

Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and MdcTrkRecon().

int MdcTrkRecon::m_maxMdcDigi [private]

Definition at line 77 of file MdcTrkRecon.h.

Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and MdcTrkRecon().

bool MdcTrkRecon::m_mcHist [private]

Definition at line 87 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and MdcTrkRecon().

MdcPrintSvc* MdcTrkRecon::m_mdcPrintSvc [private]

Definition at line 64 of file MdcTrkRecon.h.

Referenced by execute(), and initialize().

int MdcTrkRecon::m_minMdcDigi [private]

Definition at line 81 of file MdcTrkRecon.h.

Referenced by fillEvent(), fillMcTruth(), and MdcTrkRecon().

bool MdcTrkRecon::m_onlyUnusedHits [private]

Definition at line 69 of file MdcTrkRecon.h.

Referenced by ignoringUsedHits(), and MdcTrkRecon().

std::string MdcTrkRecon::m_paramFile [private]

Definition at line 55 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

std::string MdcTrkRecon::m_pdtFile [private]

Definition at line 56 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

IMagneticFieldSvc* MdcTrkRecon::m_pIMF [private]

Definition at line 123 of file MdcTrkRecon.h.

Referenced by initialize().

bool MdcTrkRecon::m_poisonHits [private]

Definition at line 70 of file MdcTrkRecon.h.

Referenced by MdcTrkRecon(), and poisoningHits().

RawDataProviderSvc* MdcTrkRecon::m_rawDataProviderSvc [private]

Definition at line 63 of file MdcTrkRecon.h.

Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and initialize().

bool MdcTrkRecon::m_recForEsTime [private]

Definition at line 73 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

std::auto_ptr<MdcSegFinder> MdcTrkRecon::m_segFinder [private]

Definition at line 62 of file MdcTrkRecon.h.

Referenced by execute(), initialize(), and ~MdcTrkRecon().

std::auto_ptr<MdcSegList> MdcTrkRecon::m_segs [private]

Definition at line 60 of file MdcTrkRecon.h.

Referenced by beginRun(), execute(), fillSegList(), and ~MdcTrkRecon().

std::vector<int> MdcTrkRecon::m_selEvtNo [private]

Definition at line 68 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

std::auto_ptr<MdcTrackListBase> MdcTrkRecon::m_tracks [private]

Definition at line 61 of file MdcTrkRecon.h.

Referenced by execute(), fillTrackList(), finalize(), initialize(), and ~MdcTrkRecon().

bool MdcTrkRecon::m_tryBunch [private]

Definition at line 72 of file MdcTrkRecon.h.

Referenced by execute(), and MdcTrkRecon().

double MdcTrkRecon::m_z0Cut [private]

Definition at line 85 of file MdcTrkRecon.h.

Referenced by initialize(), and MdcTrkRecon().

double MdcTrkRecon::mcDrift[43][288] [private]

Definition at line 94 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and fillTrackList().

double MdcTrkRecon::mcLayer[43][288] [private]

Definition at line 95 of file MdcTrkRecon.h.

double MdcTrkRecon::mcLR[43][288] [private]

Definition at line 97 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and fillTrackList().

double MdcTrkRecon::mcWire[43][288] [private]

Definition at line 96 of file MdcTrkRecon.h.

double MdcTrkRecon::mcX[43][288] [private]

Definition at line 98 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and fillTrackList().

double MdcTrkRecon::mcY[43][288] [private]

Definition at line 99 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and fillTrackList().

double MdcTrkRecon::mcZ[43][288] [private]

Definition at line 100 of file MdcTrkRecon.h.

Referenced by execute(), fillMcTruth(), and fillTrackList().

int MdcTrkRecon::pdgOfMcTk[100] [private]

Definition at line 115 of file MdcTrkRecon.h.

Referenced by fillMcTruth().

int MdcTrkRecon::t_eventNo [private]

Definition at line 111 of file MdcTrkRecon.h.

Referenced by execute(), fillEvent(), fillMcTruth(), fillSegList(), and fillTrackList().

int MdcTrkRecon::t_iExecute [private]

Definition at line 51 of file MdcTrkRecon.h.

int MdcTrkRecon::t_iExexute [private]

Definition at line 113 of file MdcTrkRecon.h.

Referenced by execute(), and initialize().

double MdcTrkRecon::t_mcTkNum [private]

Definition at line 109 of file MdcTrkRecon.h.

Referenced by fillEvent(), fillMcTruth(), and fillTrackList().

int MdcTrkRecon::t_nDigi [private]

Definition at line 110 of file MdcTrkRecon.h.

Referenced by execute(), and fillSegList().

int MdcTrkRecon::t_nHitInTk[100] [private]

Definition at line 105 of file MdcTrkRecon.h.

Referenced by fillMcTruth().

int MdcTrkRecon::t_nRecTk [private]

Definition at line 112 of file MdcTrkRecon.h.

Referenced by execute(), fillEvent(), and fillTrackList().

double MdcTrkRecon::t_t0 [private]

Definition at line 107 of file MdcTrkRecon.h.

Referenced by execute(), fillEvent(), and fillTrackList().

int MdcTrkRecon::t_t0Stat [private]

Definition at line 106 of file MdcTrkRecon.h.

Referenced by execute(), fillEvent(), and fillTrackList().

double MdcTrkRecon::t_t0Truth [private]

Definition at line 108 of file MdcTrkRecon.h.

Referenced by fillEvent(), fillMcTruth(), and fillTrackList().


Generated on Tue Nov 29 23:20:19 2016 for BOSS_7.0.2 by  doxygen 1.4.7