#include <MdcTrkRecon.h>
Definition at line 31 of file MdcTrkRecon.h.
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 }
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;}
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] |
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] |
int MdcTrkRecon::isPrimaryOfMcTk[100] [private] |
int MdcTrkRecon::m_allHit [private] |
bool MdcTrkRecon::m_arbitrateHits [private] |
BField* MdcTrkRecon::m_bfield [private] |
bool MdcTrkRecon::m_combineTracking [private] |
std::string MdcTrkRecon::m_configFile [private] |
double MdcTrkRecon::m_d0Cut [private] |
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] |
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] |
bool MdcTrkRecon::m_fieldCosmic [private] |
int MdcTrkRecon::m_fittingMethod [private] |
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] |
std::vector<float> MdcTrkRecon::m_helixHitsSigma [private] |
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] |
int MdcTrkRecon::m_minMdcDigi [private] |
Definition at line 81 of file MdcTrkRecon.h.
Referenced by fillEvent(), fillMcTruth(), and MdcTrkRecon().
bool MdcTrkRecon::m_onlyUnusedHits [private] |
std::string MdcTrkRecon::m_paramFile [private] |
std::string MdcTrkRecon::m_pdtFile [private] |
IMagneticFieldSvc* MdcTrkRecon::m_pIMF [private] |
bool MdcTrkRecon::m_poisonHits [private] |
Definition at line 63 of file MdcTrkRecon.h.
Referenced by dumpDigi(), fillEvent(), fillMcTruth(), and initialize().
bool MdcTrkRecon::m_recForEsTime [private] |
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] |
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] |
double MdcTrkRecon::m_z0Cut [private] |
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] |
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] |
double MdcTrkRecon::t_mcTkNum [private] |
Definition at line 109 of file MdcTrkRecon.h.
Referenced by fillEvent(), fillMcTruth(), and fillTrackList().
int MdcTrkRecon::t_nDigi [private] |
int MdcTrkRecon::t_nHitInTk[100] [private] |
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().