#include <MdcTrkRecon.h>
|
||||||||||||
|
00068 : 00069 Algorithm(name, pSvcLocator), 00070 m_hitData(0), m_segs(0), m_tracks(0), m_segFinder(0) 00071 { 00072 //Declare the properties-------------------------------- 00073 declareProperty("FittingMethod", m_fittingMethod = 2); 00074 declareProperty("ConfigFile", m_configFile = "MDCConfig.xml"); 00075 declareProperty("OnlyUnusedHits", m_onlyUnusedHits = 0); 00076 declareProperty("PoisonHits", m_poisonHits = 0); 00077 declareProperty("doLineFit", m_doLineFit = false); 00078 declareProperty("paramFile", m_paramFile = "params"); 00079 declareProperty("pdtFile", m_pdtFile = "pdt.table"); 00080 declareProperty("allHit", m_allHit = 1); 00081 declareProperty("hist", m_hist= false); 00082 declareProperty("recForEsTime", m_recForEsTime= false); 00083 declareProperty("tryBunch", m_tryBunch = false); 00084 declareProperty("d0Cut", m_d0Cut= -999.); 00085 declareProperty("z0Cut", m_z0Cut= -999.); 00086 declareProperty("dropTrkPt", m_dropTrkPt = -999.); 00087 declareProperty("debug", m_debug= 0); 00088 declareProperty("mcHist", m_mcHist= false); 00089 declareProperty("fieldCosmic", m_fieldCosmic= false); 00090 declareProperty("doSag", m_doSagFlag= false); 00091 declareProperty("arbitrateHits", m_arbitrateHits = true); 00092 00093 declareProperty("helixHitsSigma", m_helixHitsSigma); 00094 //for fill hist 00095 declareProperty("getDigiFlag", m_getDigiFlag = 0); 00096 declareProperty("maxMdcDigi", m_maxMdcDigi= 0); 00097 declareProperty("keepBadTdc", m_keepBadTdc= 0); 00098 declareProperty("dropHot", m_dropHot= 0); 00099 declareProperty("keepUnmatch", m_keepUnmatch= 0); 00100 declareProperty("minMdcDigi", m_minMdcDigi = 0); 00101 declareProperty("selEvtNo", m_selEvtNo); 00102 declareProperty("combineTracking",m_combineTracking = false);//yzhang 2010-05-28 00103 00104 #ifdef MDCPATREC_RESLAYER 00105 declareProperty("resLayer",m_resLayer= -1); 00106 #endif 00107 00108 }
|
|
|
00111 {
00112 m_segs.reset(0);
00113 m_segFinder.reset(0);
00114 m_hitData.reset(0);
00115 m_tracks.reset(0);
00116 }
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
00119 {
00120 MsgStream log(msgSvc(), name());
00121 log << MSG::INFO << "in beginRun()" << endreq;
00122
00123 //Detector geometry
00124 _gm = MdcDetector::instance(m_doSagFlag);
00125
00126 if(NULL == _gm) return StatusCode::FAILURE;
00127
00128 //Initialize segList
00129 m_segs.reset( new MdcSegList(_gm->nSuper(), m_flags.segPar) );
00130
00131 return StatusCode::SUCCESS;
00132 }
|
|
|
|
|
|
00467 {
00468 MsgStream log(msgSvc(), name());
00469 StatusCode sc = StatusCode::SUCCESS;
00470
00471 //Initialize
00472 for (int ii=0;ii<43;ii++){
00473 for (int jj=0;jj<288;jj++){
00474 haveDigi[ii][jj] = -2;
00475 recHitMap[ii][jj] = -2;
00476 mcDrift[ii][jj] = -99.;
00477 mcX[ii][jj] = -99.;
00478 mcY[ii][jj] = -99.;
00479 mcZ[ii][jj] = -99.;
00480 mcLR[ii][jj] = -99.;
00481 }
00482 }
00483
00484 //--------------book histogram------------------
00485 h_sfHit = histoSvc()->book( "sfHit", "Hits after segment finding",46,-1.,44. );
00486 //g_actHit = histoSvc()->book( "actHit", "Active hits",46,-1.,44. );
00487 //g_hitEff = histoSvc()->book( "hitEff", "Hit efficiency in the event",100,-0.5,1.2 );
00488 g_residVsLayer = histoSvc()->book( "residVsLayer", "Residual vs. layer",60,-5,50.,100,-0.5,1.2 );
00489 g_residVsDoca = histoSvc()->book( "residVsDoca", "Residial vs. doca",50,-2.,2.,100,-0.5,1.2 );
00490 g_segChi2 = histoSvc()->book( "segChi2", "chi2 per dof in segment finding",710,-10.,700. );
00491 g_cirTkChi2 = histoSvc()->book( "cirTkChi2", "chi2 per dof in circle finding",210,-10. ,200. );
00492 g_3dTkChi2 = histoSvc()->book( "chi2Helix", "maxChisq dof helix fit" ,710,-10,700 );
00493 g_nSigAdd = histoSvc()->book( "nSigAdd", "max allowed n sigma to add hit to existing segment",50,0,50 );
00494 g_z0Cut = histoSvc()->book( "z0Cut", "z0 for including stereo segs in cluster",1200,0,1200 );
00495 g_ctCut = histoSvc()->book( "ctCut", "ct for including stereo segs in cluster",100,-50,50 );
00496 g_delCt = histoSvc()->book( "delCt", "del ct cut forming seg group",100,-50,50 );
00497 g_delZ0 = histoSvc()->book( "delZ0", "del z0 cut forming seg group",200,-100,100 );
00498 g_phiDiff = histoSvc()->book( "phiDiff", "phidiff mult for drop dup seg",100,-0.5,6.5 );
00499 g_slopeDiff = histoSvc()->book( "slopeDiff", "slopediff for drop dup seg",100,-0.4,0.4 );
00500 g_maxSegChisqO = histoSvc()->book( "maxSegChisqO", "max chisqDof comb segs from origin" ,200,-10.,200.);
00501
00502 for(int i=0;i<11;i++){
00503 char title1[80],title2[80];
00504 sprintf(title1,"segChi2Pat3%d",i);
00505 sprintf(title2,"segChi2Pat4%d",i);
00506 g_segChi2SlayPat3[i] = histoSvc()->book( title1, "chi2/dof per layer of pat3",710,-10.,700. );
00507 g_segChi2SlayPat4[i] = histoSvc()->book( title2, "chi2/dof per layer of pat4",710,-10.,700. );
00508 }
00509
00510 //book tuple of test
00511 NTuplePtr ntt(ntupleSvc(), "FILE101/test");
00512 if ( ntt ) { m_tuplet = ntt;}
00513 else {
00514 m_tuplet = ntupleSvc()->book ("FILE101/test", CLID_ColumnWiseTuple, "MdcTrkRecon t particle");
00515 if ( m_tuplet ) {
00516 sc = m_tuplet->addItem ("layer", m_tl);
00517 sc = m_tuplet->addItem ("wire", m_tw);
00518 sc = m_tuplet->addItem ("phi", m_tphi);
00519 sc = m_tuplet->addItem ("dphi", m_tdphi);
00520 sc = m_tuplet->addItem ("dncell", m_tdncell);
00521 } else {
00522 log << MSG::ERROR << "Cannot book tuple FILE101/test" << endmsg;
00523 return StatusCode::FAILURE;
00524 }
00525 }
00526
00527 //book tuple of MC truth
00528 NTuplePtr ntMc(ntupleSvc(), "FILE101/mcTk");
00529 if ( ntMc ) { m_tupleMc = ntMc;}
00530 else {
00531 m_tupleMc = ntupleSvc()->book ("FILE101/mcTk", CLID_ColumnWiseTuple, "MdcTrkRecon Mc particle");
00532 if ( m_tupleMc ) {
00533 sc = m_tupleMc->addItem ("evtNo", m_tMcEvtNo);
00534 sc = m_tupleMc->addItem ("nDigi", m_tMcNTk);
00535 } else {
00536 log << MSG::ERROR << "Cannot book tuple FILE101/mc" << endmsg;
00537 return StatusCode::FAILURE;
00538 }
00539 }
00540
00541 //book tuple of MC truth
00542 NTuplePtr ntMcHit(ntupleSvc(), "FILE101/mcHit");
00543 if ( ntMcHit ) { m_tupleMcHit = ntMcHit;}
00544 else {
00545 m_tupleMcHit = ntupleSvc()->book ("FILE101/mcHit", CLID_ColumnWiseTuple, "MdcTrkRecon Mc hit");
00546 if ( m_tupleMcHit ) {
00547 sc = m_tupleMcHit->addItem ("tkId", m_tMcTkId);
00548 sc = m_tupleMcHit->addItem ("pid", m_tMcPid);
00549 sc = m_tupleMcHit->addItem ("px", m_tMcPx);
00550 sc = m_tupleMcHit->addItem ("py", m_tMcPy);
00551 sc = m_tupleMcHit->addItem ("pz", m_tMcPz);
00552 sc = m_tupleMcHit->addItem ("d0", m_tMcD0);
00553 sc = m_tupleMcHit->addItem ("z0", m_tMcZ0);
00554 sc = m_tupleMcHit->addItem ("theta",m_tMcTheta);
00555 sc = m_tupleMcHit->addItem ("fiterm",m_tMcFiTerm);
00556 sc = m_tupleMcHit->addItem ("nMcHit", m_tMcNHit, 0, 6796);
00557 sc = m_tupleMcHit->addIndexedItem ("layer",m_tMcNHit, m_tMcLayer);
00558 sc = m_tupleMcHit->addIndexedItem ("wire", m_tMcNHit, m_tMcWire);
00559 sc = m_tupleMcHit->addIndexedItem ("drift",m_tMcNHit, m_tMcDrift);
00560 sc = m_tupleMcHit->addIndexedItem ("x", m_tMcNHit, m_tMcX);
00561 sc = m_tupleMcHit->addIndexedItem ("y", m_tMcNHit, m_tMcY);
00562 sc = m_tupleMcHit->addIndexedItem ("z", m_tMcNHit, m_tMcZ);
00563 sc = m_tupleMcHit->addIndexedItem ("lr", m_tMcNHit, m_tMcLR);
00564 } else {
00565 log << MSG::ERROR << "Cannot book tuple FILE101/mcHit" << endmsg;
00566 return StatusCode::FAILURE;
00567 }
00568 }
00569 //book tuple of recnstruction results
00570 NTuplePtr nt1(ntupleSvc(), "FILE101/rec");
00571 if ( nt1 ) { m_tuple1 = nt1;}
00572 else {
00573 m_tuple1 = ntupleSvc()->book ("FILE101/rec", CLID_ColumnWiseTuple, "MdcTrkRecon results");
00574 if ( m_tuple1 ) {
00575 sc = m_tuple1->addItem ("t0", m_t0);
00576 sc = m_tuple1->addItem ("t0Stat", m_t0Stat);
00577 sc = m_tuple1->addItem ("t0Truth", m_t0Truth);
00578
00579 sc = m_tuple1->addItem ("nTdsTk", m_nTdsTk);
00580 sc = m_tuple1->addItem ("nMcTk", m_nMcTk);
00581
00582 sc = m_tuple1->addItem ("p", m_p);
00583 sc = m_tuple1->addItem ("pt", m_pt);
00584 sc = m_tuple1->addItem ("nSlay", m_nSlay);
00585 sc = m_tuple1->addItem ("pz", m_pz);
00586 sc = m_tuple1->addItem ("d0", m_d0);
00587 sc = m_tuple1->addItem ("phi0", m_phi0);
00588 sc = m_tuple1->addItem ("cpa", m_cpa);
00589 sc = m_tuple1->addItem ("z0", m_z0);
00590 sc = m_tuple1->addItem ("tanl", m_tanl);
00591 sc = m_tuple1->addItem ("q", m_q);
00592 sc = m_tuple1->addItem ("pocax", m_pocax);
00593 sc = m_tuple1->addItem ("pocay", m_pocay);
00594 sc = m_tuple1->addItem ("pocaz", m_pocaz);
00595
00596 sc = m_tuple1->addItem ("evtNo", m_evtNo);
00597 sc = m_tuple1->addItem ("nSt", m_nSt);
00598 sc = m_tuple1->addItem ("nAct", m_nAct);
00599 sc = m_tuple1->addItem ("nDof", m_nDof);
00600 sc = m_tuple1->addItem ("chi2", m_chi2);
00601 sc = m_tuple1->addItem ("tkId", m_tkId);
00602 sc = m_tuple1->addItem ("nHit", m_nHit, 0, 6796);
00603 sc = m_tuple1->addIndexedItem ("resid", m_nHit, m_resid);
00604 sc = m_tuple1->addIndexedItem ("sigma", m_nHit, m_sigma);
00605 sc = m_tuple1->addIndexedItem ("driftD", m_nHit, m_driftD);
00606 sc = m_tuple1->addIndexedItem ("driftT", m_nHit, m_driftT);
00607 sc = m_tuple1->addIndexedItem ("doca", m_nHit, m_doca);
00608 sc = m_tuple1->addIndexedItem ("entra", m_nHit, m_entra);
00609 sc = m_tuple1->addIndexedItem ("ambig", m_nHit, m_ambig);
00610 sc = m_tuple1->addIndexedItem ("fltLen", m_nHit, m_fltLen);
00611 sc = m_tuple1->addIndexedItem ("tof", m_nHit, m_tof);
00612 sc = m_tuple1->addIndexedItem ("act", m_nHit, m_act);
00613 sc = m_tuple1->addIndexedItem ("tdc", m_nHit, m_tdc);
00614 sc = m_tuple1->addIndexedItem ("adc", m_nHit, m_adc);
00615 sc = m_tuple1->addIndexedItem ("layer", m_nHit, m_layer);
00616 sc = m_tuple1->addIndexedItem ("wire", m_nHit, m_wire);
00617 sc = m_tuple1->addIndexedItem ("x", m_nHit, m_x);
00618 sc = m_tuple1->addIndexedItem ("y", m_nHit, m_y);
00619 sc = m_tuple1->addIndexedItem ("z", m_nHit, m_z);
00620 sc = m_tuple1->addIndexedItem ("dx", m_nHit, m_dx);
00621 sc = m_tuple1->addIndexedItem ("dy", m_nHit, m_dy);
00622 sc = m_tuple1->addIndexedItem ("dz", m_nHit, m_dz);
00623 sc = m_tuple1->addIndexedItem ("dDriftD", m_nHit, m_dDriftD);
00624 sc = m_tuple1->addIndexedItem ("dlr", m_nHit, m_dlr);
00625 } else {
00626 log << MSG::ERROR << "Cannot book tuple FILE101/rec" << endmsg;
00627 return StatusCode::FAILURE;
00628 }
00629 }
00630 //book tuple of segment
00631 NTuplePtr ntSeg(ntupleSvc(), "FILE101/seg");
00632 if ( ntSeg ) { m_tupleSeg = ntSeg;}
00633 else {
00634 m_tupleSeg = ntupleSvc()->book ("FILE101/seg", CLID_ColumnWiseTuple, "MdcTrkRecon segment data");
00635 if ( m_tupleSeg ) {
00636 sc = m_tupleSeg->addItem ("evtNo", m_tsEvtNo);
00637 sc = m_tupleSeg->addItem ("nDigi", m_tsNDigi, 0, 6796);
00638 sc = m_tupleSeg->addIndexedItem ("layer", m_tsNDigi, m_tsLayer);
00639 sc = m_tupleSeg->addIndexedItem ("wire", m_tsNDigi, m_tsWire);
00640 sc = m_tupleSeg->addIndexedItem ("inSeg", m_tsNDigi, m_tsInSeg);
00641 sc = m_tupleSeg->addIndexedItem ("mcTkId", m_tsNDigi, m_tsMcTkId);
00642 } else {
00643 log << MSG::ERROR << "Cannot book tuple FILE101/seg" << endmsg;
00644 return StatusCode::FAILURE;
00645 }
00646 }
00647
00648 //book tuple of event
00649 NTuplePtr nt4(ntupleSvc(), "FILE101/evt");
00650 if ( nt4 ) { m_tupleEvt = nt4;}
00651 else {
00652 m_tupleEvt = ntupleSvc()->book ("FILE101/evt", CLID_ColumnWiseTuple, "MdcTrkRecon event data");
00653 if ( m_tupleEvt ) {
00654 sc = m_tupleEvt->addItem ("evtNo", m_t4EvtNo);
00655 sc = m_tupleEvt->addItem ("nMcTk", m_t4nMcTk );
00656 sc = m_tupleEvt->addItem ("nTdsTk", m_t4nRecTk );
00657 sc = m_tupleEvt->addItem ("t0", m_t4t0);
00658 sc = m_tupleEvt->addItem ("t0Stat", m_t4t0Stat);
00659 sc = m_tupleEvt->addItem ("t0Truth", m_t4t0Truth);
00660 sc = m_tupleEvt->addItem ("nDigiUn", m_t4nDigiUnmatch);
00661 sc = m_tupleEvt->addItem ("nDigi", m_t4nDigi, 0, 6796);
00662 sc = m_tupleEvt->addIndexedItem ("layer", m_t4nDigi, m_t4Layer);
00663 sc = m_tupleEvt->addIndexedItem ("wire", m_t4nDigi, m_t4Wire);
00664 sc = m_tupleEvt->addIndexedItem ("rt", m_t4nDigi, m_t4Time);
00665 sc = m_tupleEvt->addIndexedItem ("rc", m_t4nDigi, m_t4Charge);
00666 sc = m_tupleEvt->addIndexedItem ("phiMid", m_t4nDigi, m_t4PhiMid);
00667 sc = m_tupleEvt->addIndexedItem ("hot", m_t4nDigi, m_t4Hot);
00668 //sc = m_tupleEvt->addIndexedItem ("recHit", m_t4nDigi, m_t4recHit);
00669 //sc = m_tupleEvt->addIndexedItem ("rawHit", m_t4nDigi, m_t4rawHit);
00670 } else {
00671 log << MSG::ERROR << "Cannot book N-tuple: FILE101/mc" << endmsg;
00672 return StatusCode::FAILURE;
00673 }
00674 }
00675
00676 //book tuple of residula for every layer
00677 NTuplePtr ntCombAx(ntupleSvc(), "FILE101/combAx");
00678 if ( ntCombAx ) { g_tupleCombAx= ntCombAx;}
00679 else {
00680 g_tupleCombAx = ntupleSvc()->book ("FILE101/combAx", CLID_RowWiseTuple, "MdcTrkRecon cuts in combine ax seg");
00681 if ( g_tupleCombAx) {
00682 sc = g_tupleCombAx->addItem ("dPhi0", g_combAxdPhi0);
00683 sc = g_tupleCombAx->addItem ("dCurv", g_combAxdCurv);
00684 sc = g_tupleCombAx->addItem ("sigPhi0",g_combAxSigPhi0);
00685 sc = g_tupleCombAx->addItem ("sigCurv",g_combAxSigCurv);
00686 sc = g_tupleCombAx->addItem ("slSeed", g_combAxSlSeed);
00687 sc = g_tupleCombAx->addItem ("slTest", g_combAxSlTest);
00688 sc = g_tupleCombAx->addItem ("qSeed", g_combAxQualitySeed);
00689 sc = g_tupleCombAx->addItem ("qTest", g_combAxQualityTest);
00690 sc = g_tupleCombAx->addItem ("mc", g_combAxMc);
00691 } else {
00692 log << MSG::ERROR << "Cannot book N-tuple: FILE101/combAx" << endmsg;
00693 return StatusCode::FAILURE;
00694 }
00695 }
00696 #ifdef MDCPATREC_TIMETEST
00697 //book tuple of time test
00698 NTuplePtr nt5(ntupleSvc(), "FILE101/ti");
00699 if ( nt5 ) { m_tupleTime = nt5;}
00700 else {
00701 m_tupleTime = ntupleSvc()->book ("FILE101/ti", CLID_RowWiseTuple, "MdcTrkRecon time");
00702 if ( m_tupleTime ) {
00703 sc = m_tupleTime->addItem ("eventtime", m_eventTime);
00704 sc = m_tupleTime->addItem ("recTkNum", m_t5recTkNum);
00705 sc = m_tupleTime->addItem ("mcTkNum", m_mcTkNum);
00706 sc = m_tupleTime->addItem ("evtNo", m_t5EvtNo);
00707 sc = m_tupleTime->addItem ("nHit", m_t5nHit);
00708 sc = m_tupleTime->addItem ("nDigi", m_t5nDigi);
00709 } else {
00710 log << MSG::ERROR << "Cannot book N-tuple: FILE101/time" << endmsg;
00711 return StatusCode::FAILURE;
00712 }
00713 }
00714 sc = service( "BesTimerSvc", m_timersvc);
00715 if( sc.isFailure() ) {
00716 log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
00717 return StatusCode::FAILURE;
00718 }
00719 m_timer[1] = m_timersvc->addItem("Execution");
00720 m_timer[1]->propName("nExecution");
00721 #endif
00722
00723 #ifdef MDCPATREC_RESLAYER
00724 //book tuple of residula for every layer
00725 NTuplePtr nt6(ntupleSvc(), "FILE101/res");
00726 if ( nt6 ) { g_tupleres = nt6;}
00727 else {
00728 g_tupleres= ntupleSvc()->book ("FILE101/res", CLID_RowWiseTuple, "MdcTrkRecon res");
00729 if ( g_tupleres ) {
00730 sc = g_tupleres->addItem ("resid", g_resLayer);
00731 sc = g_tupleres->addItem ("layer", g_t6Layer);
00732 } else {
00733 log << MSG::ERROR << "Cannot book N-tuple: FILE104/res" << endmsg;
00734 return StatusCode::FAILURE;
00735 }
00736 }
00737 #endif
00738 return sc;
00739 }// end of bookNTuple()
|
|
|
|
|
|
01107 {
01108 uint32_t getDigiFlag = 0;
01109 getDigiFlag += m_maxMdcDigi;
01110 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
01111 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
01112 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
01113 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
01114 MdcDigiCol::iterator iter = mdcDigiVec.begin();
01115 std::cout<<name()<<" dump MdcDigiVec, nDigi="<<mdcDigiVec.size()<<std::endl;
01116 for (int iDigi=0;iter!= mdcDigiVec.end(); iter++,iDigi++ ) {
01117 long l = MdcID::layer((*iter)->identify());
01118 long w = MdcID::wire((*iter)->identify());
01119 int tkTruth = (*iter)->getTrackIndex();
01120 std::cout<<"("<<l<<","<<w<<";"<<tkTruth<<")";
01121 if(iDigi%4==0) std::cout<<std::endl;
01122 }
01123 std::cout<<std::endl;
01124 }
|
|
|
|
|
|
00892 {
00893 std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
00894 RecMdcTrackCol::iterator it = trackList->begin();
00895 for (;it!= trackList->end();it++){
00896 RecMdcTrack *tk = *it;
00897 std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
00898 cout <<" d0 "<<tk->helix(0)
00899 <<" phi0 "<<tk->helix(1)
00900 <<" cpa "<<tk->helix(2)
00901 <<" z0 "<<tk->helix(3)
00902 <<" tanl "<<tk->helix(4)
00903 <<endl;
00904 std::cout<<" q "<<tk->charge()
00905 <<" theta "<<tk->theta()
00906 <<" phi "<<tk->phi()
00907 <<" x0 "<<tk->x()
00908 <<" y0 "<<tk->y()
00909 <<" z0 "<<tk->z()
00910 <<" r0 "<<tk->r()
00911 <<endl;
00912 std::cout <<" p "<<tk->p()
00913 <<" pt "<<tk->pxy()
00914 <<" px "<<tk->px()
00915 <<" py "<<tk->py()
00916 <<" pz "<<tk->pz()
00917 <<endl;
00918 std::cout<<" tkStat "<<tk->stat()
00919 <<" chi2 "<<tk->chi2()
00920 <<" ndof "<<tk->ndof()
00921 <<" nhit "<<tk->getNhits()
00922 <<" nst "<<tk->nster()
00923 <<endl;
00924 std::cout<< "errmat mat " << std::endl;
00925 std::cout<< tk->err()<<std::endl;
00926 std::cout<< "errmat array " << std::endl;
00927 for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
00928 std::cout<< " " << std::endl;
00929
00930 int nhits = tk->getVecHits().size();
00931 std::cout<<nhits <<" Hits: " << std::endl;
00932 for(int ii=0; ii <nhits ; ii++){
00933 Identifier id(tk->getVecHits()[ii]->getMdcId());
00934 int layer = MdcID::layer(id);
00935 int wire = MdcID::wire(id);
00936 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
00937 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
00938 }//end of hit list
00939 std::cout << " "<< std::endl;
00940 }//end of tk list
00941 std::cout << " "<< std::endl;
00942 }//end of dumpTdsTrack
|
|
|
|
|
|
00202 {
00203 setFilterPassed(false);
00204
00205 MsgStream log(msgSvc(), name());
00206 log << MSG::INFO << "in execute()" << endreq;
00207
00208 TrkContextEv context(m_bfield);
00209 #ifdef MDCPATREC_TIMETEST
00210 m_timer[1]->start();
00211 ti_nHit=0;
00212 #endif
00213 //------------------read event header--------------
00214 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00215 if (!evtHead) {
00216 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00217 return( StatusCode::FAILURE);
00218 }
00219 t_eventNo = evtHead->eventNumber();
00220
00221 if (m_selEvtNo.size() >0){
00222 std::vector<int>::iterator it = m_selEvtNo.begin();
00223 for (; it!=m_selEvtNo.end(); it++){
00224 if((*it)== t_eventNo) setFilterPassed(true);
00225 }
00226 }
00227 //------------------get event start time-----------
00228 double t0 = 0.;
00229 t_t0 = -1.;
00230 t_t0Stat = -1;
00231 bool iterateT0 = false;
00232 const int nBunch = 3;//3 bunch/event
00233 double iBunch = 0 ; //form 0 to 2
00234 double BeamDeltaTime = 24. / nBunch;// 8ns
00235 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00236 if (!evTimeCol || evTimeCol->size()==0) {
00237 log << MSG::WARNING <<t_eventNo<< " Could not find RecEsTimeCol" << endreq;
00238 if(m_fieldCosmic){//yzhang temp for bes3 csmc test
00239 return StatusCode::SUCCESS;
00240 }
00241 }
00242 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00243 // skip RecEsTimeCol no rec data
00244 if (iter_evt != evTimeCol->end()){
00245 t_t0Stat = (*iter_evt)->getStat();
00246 t_t0 = (*iter_evt)->getTest();
00247
00248 if (m_doLineFit){
00249 //t0 -= m_t0Const;
00250 if (abs(t_t0)<0.00001 || (t_t0 < 0) || (t_t0 > 9999.0) ){
00251 log << MSG::WARNING << "Skip with t0 = "<<t_t0 << endreq;
00252 return StatusCode::SUCCESS;//for bes3 cosmic test
00253 }
00254 }
00255 t0 = t_t0*1.e-9;
00256 }
00257
00258
00259 restart:
00260 if ( !m_recForEsTime && m_tryBunch) {
00261 if ( t_t0Stat < 10 ) return StatusCode::SUCCESS;
00262 }
00263 if ( m_tryBunch ){
00264 if ( iBunch > (nBunch - 1) ) return StatusCode::SUCCESS;
00265 if ( t_t0Stat < 0 ){ iterateT0 = true; }
00266 if ( iterateT0 ){
00267 //double EventDeltaTime = 24.;//24ns/event
00268 if ( (t_t0Stat > -1) && (fabs( iBunch * BeamDeltaTime - t_t0) < 2.) ){
00269 ++iBunch;
00270 goto restart;
00271 }else{ t0 = iBunch * BeamDeltaTime *1.e-9;}
00272 ++iBunch;
00273 }
00274 }
00275
00276 SmartDataPtr<MdcHitMap> hitMap(eventSvc(),"/Event/Hit/MdcHitMap");
00277 if (!hitMap) {
00278 log << MSG::FATAL << "Could not retrieve MDC hit map" << endreq;
00279 return( StatusCode::FAILURE);
00280 }
00281 //---------------------Read an event--------------
00282 SmartDataPtr<MdcHitCol> hitCol(eventSvc(),"/Event/Hit/MdcHitCol");
00283 if (!hitCol) {
00284 log << MSG::FATAL << "Could not retrieve MDC hit list" << endreq;
00285 return( StatusCode::FAILURE);
00286 }
00287 StatusCode sc;
00288
00289 //---------- register RecMdcTrackCol and RecMdcHitCol to tds---------
00290 DataObject *aReconEvent;
00291 eventSvc()->findObject("/Event/Recon",aReconEvent);
00292 if(aReconEvent==NULL) {
00293 aReconEvent = new ReconEvent();
00294 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00295 if(sc!=StatusCode::SUCCESS) {
00296 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00297 return( StatusCode::FAILURE);
00298 }
00299 }
00300
00301 DataObject *aTrackCol;
00302 DataObject *aRecHitCol;
00303 //yzhang 2010-05-28
00304 if(!m_combineTracking){
00305 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00306 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00307 if(aTrackCol != NULL) {
00308 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00309 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00310 }
00311 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00312 if(aRecHitCol != NULL) {
00313 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00314 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00315 }
00316 }
00317
00318 RecMdcTrackCol* trackList;
00319 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00320 if (aTrackCol) {
00321 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00322 }else{
00323 trackList = new RecMdcTrackCol;
00324 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00325 if(!sc.isSuccess()) {
00326 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00327 return StatusCode::FAILURE;
00328 }
00329 }
00330 RecMdcHitCol* hitList;
00331 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00332 if (aRecHitCol) {
00333 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00334 }else{
00335 hitList = new RecMdcHitCol;
00336 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00337 if(!sc.isSuccess()) {
00338 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00339 return StatusCode::FAILURE;
00340 }
00341 }
00342
00343 m_hitData->loadevent(hitCol, hitMap, t0);
00344 t_nDigi = hitCol->size();
00345
00346 if (poisoningHits()) {
00347 m_hitData->poisonHits(_gm);
00348 }
00349 #ifdef MDCPATREC_TIMETEST
00350 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00351 if(!mcParticleCol){
00352 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
00353 }
00354 m_mcTkNum = mcParticleCol->size();
00355 #endif
00356
00357 if(m_debug>1) dumpDigi();
00358 if(m_flags.lHist ) fillMcTruth();
00359 //--------------------------------------------------------------------------
00360 // Segment finding
00361 //--------------------------------------------------------------------------
00362 int nSegs = 0;
00363 if (m_flags.findSegs) {
00364 // Form track segments in superlayers
00365 nSegs = m_segFinder->createSegs(_gm, *m_segs, m_hitData->segUsage(),
00366 m_hitData->hitMap(), t0);
00367 }
00368 log << MSG::INFO << " Found " << nSegs << " segments" << endreq;
00369
00370 if (m_debug>1){
00371 std::cout<<"------------------------------------------------"<<std::endl;
00372 std::cout<<"==event "<<t_eventNo<< " Found " << nSegs << " segments" << std::endl;
00373 m_segs->plot();
00374 }
00375
00376 if (m_flags.lHist||m_flags.segPar.lPrint) {
00377 fillSegList();
00378 }
00379
00380 //--------------------------------------------------------------------------
00381 // Combine segments to form track
00382 //--------------------------------------------------------------------------
00383 int nTracks = 0;
00384 int nDeleted = 0;
00385 if (m_flags.findTracks && m_flags.findSegs) {
00386 if (nSegs > 2) {
00387 nTracks = m_tracks->createFromSegs(m_segs.get(), m_hitData->hitMap(),
00388 _gm, context, t0);
00389 }
00390
00391 if(m_arbitrateHits) nDeleted = m_tracks->arbitrateHits();
00392 int nKeep = nTracks - nDeleted;
00393
00394 if (m_debug>0 && (nTracks - nDeleted)>0){
00395 std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
00396 <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
00397 << " keep "<< nKeep <<" track(s)";
00398 if (nDeleted!=0){ std::cout <<",deleted " <<nDeleted; }
00399 std::cout<< std::endl;//yzhang debug
00400 }else{
00401 if (m_debug>0){
00402 std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
00403 <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
00404 <<" found 0 track "<< endl;
00405 }
00406 }
00407
00408 if (m_flags.lHist) t_nRecTk = nKeep;
00409
00410 #ifdef MDCPATREC_RESLAYER
00411 m_tracks->setResLayer(m_resLayer);
00412 #endif
00413 if (m_flags.lHist){ fillTrackList(); }
00414 // Stick the found tracks into the list in TDS
00415 m_tracks->store(trackList, hitList);
00416 if (m_debug>1) { dumpTdsTrack(trackList); }
00417
00418 //if(nKeep != (int)trackList->size()) std::cout<<__FILE__<<" "<<__LINE__<<"!!!!!!!!!!!!!!!!! nKeep != nTdsTk"<< std::endl;
00419 #ifdef MDCPATREC_TIMETEST
00420 RecMdcTrackCol::iterator iter = trackList->begin();
00421 for (;iter != trackList->end(); iter++) {
00422 MdcTrack* tk = *iter;
00423 ti_nHit += tk->getNhits();
00424 }
00425 #endif
00426
00427 }
00428 //-------------End track-finding------------------
00429 m_segs->destroySegs();
00430
00431 //if ( t_eventNo% 1000 == 0) {
00432 //std::cout<<"evt "<<t_eventNo<< " Found " << trackList->size() << " track(s)" <<endl;//yzhang debug
00433 //}
00434 #ifdef MDCPATREC_TIMETEST
00435 m_timer[1]->stop();
00436 //cout << "m_timer[1]->elapsed()::" << m_timer[1]->elapsed() << endl;
00437 //cout << "m_timer[1]->mean()::" << m_timer[1]->mean() << endl;
00438 m_eventTime = m_timer[1]->elapsed();
00439 m_t5recTkNum = m_tracks->length();
00440 m_t5EvtNo = t_eventNo;
00441 m_t5nHit = ti_nHit;
00442 m_t5nDigi = t_nDigi;
00443 sc = m_tupleTime->write();
00444 #endif
00445 // for event start time, if track not found try another t0
00446 if ( m_tryBunch ){
00447 if ( nTracks <1 ){ iterateT0 = true;
00448 goto restart;
00449 }
00450 }
00451 if (m_flags.lHist ) fillEvent();
00452 return StatusCode::SUCCESS;
00453 }
|
|
|
|
|
|
01064 {
01065 if (m_tupleEvt!=NULL){
01066
01067 uint32_t getDigiFlag = 0;
01068 getDigiFlag += m_maxMdcDigi;
01069 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
01070 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
01071 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
01072 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
01073 if ((int)mdcDigiVec.size()<m_minMdcDigi){
01074 std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
01075 }
01076
01077 m_t4nDigi = mdcDigiVec.size();
01078 int iDigi=0;
01079 MdcDigiCol::iterator iter = mdcDigiVec.begin();
01080 for (;iDigi<m_t4nDigi; iter++ ) {
01081 double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
01082 double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
01083 m_t4Time[iDigi] = tdc;
01084 m_t4Charge[iDigi] = adc;
01085 long l = MdcID::layer((*iter)->identify());
01086 long w = MdcID::wire((*iter)->identify());
01087 m_t4Layer[iDigi] = l;
01088 m_t4Wire[iDigi] = w;
01089 m_t4PhiMid[iDigi] = _gm->Layer(l)->phiWire(w);
01090 m_t4Hot[iDigi] = recHitMap[l][w];
01091 iDigi++;
01092 }
01093 m_t4t0 = t_t0;
01094 m_t4t0Stat = t_t0Stat;
01095 m_t4t0Truth = t_t0Truth;
01096 m_t4EvtNo = t_eventNo;
01097 m_t4nRecTk = t_nRecTk;
01098 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
01099 if (mdcDigiCol){
01100 m_t4nDigiUnmatch = mdcDigiCol->size();
01101 }
01102 m_t4nMcTk= t_mcTkNum;
01103 m_tupleEvt->write();
01104 }
01105 }//end of fillEvent
|
|
|
|
|
|
00741 {
00742 MsgStream log(msgSvc(), name());
00743 StatusCode sc;
00744 t_mcTkNum = 0;
00745 if (m_mcHist){
00746 //------------------Retrieve MC truth MdcMcHit------------
00747 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
00748 if (!mcMdcMcHitCol) {
00749 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
00750 }else{
00751 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
00752 for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
00753 const Identifier id= (*iter_mchit)->identify();
00754 int layer = MdcID::layer(id);
00755 int wire = MdcID::wire(id);
00756 int iMcTk = (*iter_mchit)->getTrackIndex();
00757 //int iHit = t_nHitInTk[iMcTk];
00758 mcDrift[layer][wire]= (*iter_mchit)->getDriftDistance()/10.; //drift in MC.
00759 mcX[layer][wire] = (*iter_mchit)->getPositionX()/10.;
00760 mcY[layer][wire] = (*iter_mchit)->getPositionY()/10.;
00761 mcZ[layer][wire] = (*iter_mchit)->getPositionZ()/10.;
00762 mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
00763 //std::cout<< "lr== "<< (*iter_mchit)->getPositionFlag() << std::endl;//yzhang debug
00764 if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
00765 t_nHitInTk[iMcTk]++;
00766 }
00767 }
00768
00769
00770 //------------------get event start time truth-----------
00771 t_t0Truth = -10.;
00772 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00773 if(!mcParticleCol){
00774 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
00775 }else {
00776 int nMcTk = 0;
00777 Event::McParticleCol::iterator it = mcParticleCol->begin();
00778 for (;it != mcParticleCol->end(); it++){
00779 int pdg_code = (*it)->particleProperty();
00780 if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
00781 t_mcTkNum++;
00782 nMcTk++;
00783 }
00784 m_tMcNTk= nMcTk;
00785 it = mcParticleCol->begin();
00786 for (;it != mcParticleCol->end(); it++){
00787 if ((*it)->primaryParticle()){
00788 t_t0Truth = (*it)->initialPosition().t();
00789 }
00790 //fill charged particle
00791 int pdg_code = (*it)->particleProperty();
00792 //FIXME skip charged track and track outside MDC
00793
00794 if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
00795 if (m_tupleMcHit!=NULL){
00796 const CLHEP::Hep3Vector& true_mom = (*it)->initialFourMomentum().vect();
00797 const CLHEP::HepLorentzVector& posIni = (*it)->initialPosition();
00798 const CLHEP::HepLorentzVector& posFin = (*it)->finalPosition();
00799 m_tMcPx = true_mom.x();
00800 m_tMcPy = true_mom.y();
00801 m_tMcPz = true_mom.z();
00802 m_tMcD0 = posIni.perp()/10.;
00803 m_tMcZ0 = posIni.z()/10.;
00804 m_tMcTheta = (*it)->initialFourMomentum().theta();
00805 m_tMcFiTerm= posFin.phi();
00806 m_tMcPid = pdg_code;
00807 int tkId = (*it)->trackIndex();
00808 m_tMcTkId = tkId;
00809 m_tMcNHit = t_nHitInTk[tkId];
00810 //std::cout<< "t_nHitInTk "<<t_nHitInTk[tkId]<<" "<<tkId << std::endl;//yzhang debug
00811 //fill mcHit tuples
00812 /*
00813 for (int i =0; i<t_nHitInTk[tkId]; i++){
00814 m_tMcLayer[i] = mcLayer[tkId][i];
00815 m_tMcWire[i] = mcWire[tkId][i];
00816 m_tMcDrift[i] = mcDrift[tkId][i];
00817 m_tMcX[i] = mcX[tkId][i];
00818 m_tMcY[i] = mcY[tkId][i];
00819 m_tMcZ[i] = mcZ[tkId][i];
00820 m_tMcLR[i] = mcLR[tkId][i];
00821 }
00822 */
00823 }
00824 m_tupleMcHit->write();
00825 }//end of loop mcParticleCol
00826 m_tMcEvtNo = t_eventNo;
00827 m_tupleMc->write();
00828 }
00829 }
00830
00831 uint32_t getDigiFlag = 0;
00832 getDigiFlag += m_maxMdcDigi;
00833 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00834 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00835 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00836 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00837 if ((int)mdcDigiVec.size()<m_minMdcDigi){
00838 std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
00839 }
00840
00841 if (mdcDigiVec.size()<=0) return;
00842 //fill raw digi and t0 status
00843 //int nDigi = mdcDigiVec.size();
00844 MdcDigiCol::iterator iter = mdcDigiVec.begin();
00845 for (;iter!= mdcDigiVec.end(); iter++ ) {
00846 long l = MdcID::layer((*iter)->identify());
00847 long w = MdcID::wire((*iter)->identify());
00848 haveDigi[l][w]=(*iter)->getTrackIndex();
00849 }
00850 }//end of fillMcTruth()
|
|
|
|
|
|
00852 {
00853 if (2 == m_flags.segPar.lPrint) {
00854 std::cout << "print after Segment finding " << std::endl;
00855 //m_segs->plot(10);
00856 }
00857 if (!m_flags.lHist) return;
00858 // Fill hits of every layer after segment finding
00859 int hitInSegList[43][288];
00860 for (int ii=0;ii<43;ii++){
00861 for (int jj=0;jj<288;jj++){ hitInSegList[ii][jj] = 0; }
00862 }
00863 for (int i = 0; i < _gm->nSuper(); i++) {
00864 MdcSeg* aSeg = (MdcSeg *) m_segs->oneList(i)->first();
00865 while (aSeg != NULL) {
00866 for (int ihit=0 ; ihit< aSeg->nHit() ; ihit++){
00867 const MdcHit* hit = aSeg->hit(ihit)->mdcHit();
00868 int layer = hit->layernumber();
00869 int wire = hit->wirenumber();
00870 hitInSegList[layer][wire]++;
00871 }
00872 aSeg = (MdcSeg *) aSeg->next();
00873 }
00874 }//end for slayer
00875 m_tsEvtNo = t_eventNo;
00876 m_tsNDigi = t_nDigi;
00877 int iDigi=0;
00878 for (int ii=0;ii<43;ii++){
00879 for (int jj=0;jj<288;jj++){
00880 if (haveDigi[ii][jj] > -2){
00881 m_tsLayer[iDigi] = ii;
00882 m_tsWire[iDigi] = jj;
00883 m_tsInSeg[iDigi] = hitInSegList[ii][jj];
00884 m_tsMcTkId[iDigi] = haveDigi[ii][jj];
00885 iDigi++;
00886 }
00887 }
00888 }
00889 m_tupleSeg->write();
00890 }//end of fillSegList
|
|
|
|
|
|
00944 {
00945 if (m_flags.debugFlag()>1) {
00946 std::cout << "=======Print after Track finding=======" << std::endl;
00947 m_tracks->plot();
00948 }
00949
00950 if (!m_flags.lHist) return;
00951 //----------- fill hit -----------
00952 int nTk = (*m_tracks).nTrack();
00953 for (int itrack = 0; itrack < nTk; itrack++) {
00954 MdcTrack *atrack = (*m_tracks)[itrack];
00955 if (atrack==NULL) continue;
00956 const TrkFit* fitResult = atrack->track().fitResult();
00957 if (fitResult == 0) continue;//check the fit worked
00958
00959 //for fill ntuples
00960 int nSt=0; //int nSeg=0;
00961 int seg[11] = {0}; int segme;
00962 //-----------hit list-------------
00963 const TrkHitList* hitlist = atrack->track().hits();
00964 TrkHitList::hot_iterator hot = hitlist->begin();
00965 int layerCount[43]={0};
00966 for(int iii=0;iii<43;iii++){layerCount[iii]=0;}
00967 int i=0;
00968 for (;hot!= hitlist->end();hot++){
00969 const MdcRecoHitOnTrack* recoHot;
00970 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
00971 int layer = recoHot->mdcHit()->layernumber();
00972 int wire = recoHot->mdcHit()->wirenumber();
00973 m_layer[i] = layer;
00974 m_wire[i] = wire;
00975 layerCount[layer]++;
00976 recHitMap[layer][wire]++;
00977 m_ambig[i] = recoHot->wireAmbig();// wire ambig
00978 //fltLen
00979 double fltLen = recoHot->fltLen();
00980 m_fltLen[i] = fltLen;
00981 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
00982 //position
00983 HepPoint3D pos; Hep3Vector dir;
00984 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
00985 m_x[i] = pos.x();
00986 m_y[i] = pos.y();
00987 m_z[i] = pos.z();
00988 m_driftT[i] = recoHot->mdcHit()->driftTime(tof,pos.z());
00989 m_tof[i] = tof;
00990 m_driftD[i] = recoHot->drift();
00991 m_sigma[i] = recoHot->hitRms();
00992 m_tdc[i] = recoHot->rawTime();
00993 m_adc[i] = recoHot->mdcHit()->charge();
00994 m_doca[i] = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME
00995 m_entra[i] = recoHot->entranceAngle();
00996 m_act[i] = recoHot->isActive();
00997 //diff with truth
00998 m_dx[i] = m_x[i] - mcX[layer][wire];
00999 m_dy[i] = m_y[i] - mcY[layer][wire];
01000 m_dz[i] = m_z[i] - mcZ[layer][wire];
01001 m_dDriftD[i] = m_driftD[i] - mcDrift[layer][wire];
01002 m_dlr[i] = m_ambig[i] - mcLR[layer][wire];
01003 //resid
01004 double res=999.,rese=999.;
01005 if (recoHot->resid(res,rese,false)){
01006 }else{}
01007 m_resid[i] = res;
01008 //for n seg
01009 segme=0;
01010 if ( layer >0 ) {segme=(layer-1)/4;}
01011 seg[segme]++;
01012 if (recoHot->layer()->view()) { ++nSt; }
01013 i++;
01014 }// end fill hit
01015 int nSlay=0;
01016 for(int i=0;i<11;i++){
01017 if (seg[i]>0) nSlay++;
01018 }
01019 //------------fill track result-------------
01020 m_tkId = itrack;
01021 //track parameters at closest approach to beamline
01022 double fltLenPoca = 0.0;
01023 TrkExchangePar helix = fitResult->helix(fltLenPoca);
01024 m_q = fitResult->charge();
01025 m_phi0 = helix.phi0();
01026 m_tanl = helix.tanDip();
01027 m_z0 = helix.z0();
01028 m_d0 = helix.d0();
01029 m_pt = fitResult->pt();
01030 m_nSlay = nSlay;
01031 if (m_pt > 0.00001){
01032 m_cpa = fitResult->charge()/fitResult->pt();
01033 }
01034 //momenta and position
01035 CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
01036 double px,py,pz,pxy;
01037 pxy = fitResult->pt();
01038 px = p1.x();
01039 py = p1.y();
01040 pz = p1.z();
01041 m_p = p1.mag();
01042 m_pz = pz;
01043 HepPoint3D poca = fitResult->position(fltLenPoca);
01044 m_pocax = poca.x();
01045 m_pocay = poca.y();
01046 m_pocaz = poca.z();
01047 m_nAct = fitResult->nActive();
01048 m_nHit = hitlist->nHit();
01049 m_nDof = fitResult->nDof();
01050 m_nSt = nSt;
01051 m_chi2 = fitResult->chisq();
01052 //for (int l=0;l<43;l++) m_layerCount[l] = layerCount[l];
01053 m_t0 = t_t0;
01054 m_t0Stat = t_t0Stat;
01055 m_t0Truth = t_t0Truth;
01056 m_nTdsTk = t_nRecTk;
01057 m_nMcTk = t_mcTkNum;
01058 m_evtNo = t_eventNo;
01059 m_tuple1->write();
01060 }//end of loop rec trk list
01061 }//end of fillTrackList
|
|
|
|
|
|
00456 {
00457 MsgStream log(msgSvc(), name());
00458 log << MSG::INFO << "in finalize()" << endreq;
00459
00460 m_tracks.reset(0);
00461 #ifdef MDCPATREC_TIMETEST
00462 m_timersvc->print();
00463 #endif
00464 return StatusCode::SUCCESS;
00465 }
|
|
|
00040 {return m_onlyUnusedHits;}
|
|
|
00040 {return m_onlyUnusedHits;}
|
|
|
|
|
|
00135 {
00136
00137 MsgStream log(msgSvc(), name());
00138 log << MSG::INFO << "in initialize()" << endreq;
00139 StatusCode sc;
00140
00141 //Pdt
00142 Pdt::readMCppTable(m_pdtFile);
00143
00144 //Flag and Pars
00145 m_flags.readPar(m_paramFile);
00146 m_flags.lHist = m_hist;
00147 m_flags.setDebug(m_debug);
00148 for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
00149 if (!m_doLineFit) MdcTrackListBase::m_d0Cut = m_d0Cut;
00150 if (!m_doLineFit) MdcTrackListBase::m_z0Cut = m_z0Cut;
00151 MdcTrackListBase::m_ptCut = m_dropTrkPt;
00152 if (m_debug>0) {
00153 std::cout<< "MdcTrkRecon d0Cut "<<m_d0Cut<< " z0Cut "<<
00154 m_z0Cut<<" ptCut "<<m_dropTrkPt << std::endl;
00155 m_flags.printPar();
00156 }
00157
00158 //Initailize magnetic filed
00159 sc = service ("MagneticFieldSvc",m_pIMF);
00160 if(sc!=StatusCode::SUCCESS) {
00161 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00162 }
00163 m_bfield = new BField(m_pIMF);
00164 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00165
00166 //Initialize RawDataProviderSvc
00167 IRawDataProviderSvc* irawDataProviderSvc;
00168 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00169 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00170 if ( sc.isFailure() ){
00171 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00172 return StatusCode::FAILURE;
00173 }
00174
00175 //Initialize hitdata
00176 m_hitData.reset( new MdcSegData( ignoringUsedHits() ));
00177
00178 //Initialize segFinder
00179 if (m_flags.findSegs) {
00180 m_segFinder.reset( new MdcSegFinder(m_flags.segPar.useAllAmbig) );
00181 }
00182
00183 //Initialize trackList
00184 if ( m_doLineFit ){
00185 m_tracks.reset( new MdcTrackListCsmc(m_flags.tkParTight) );
00186 } else {
00187 m_tracks.reset( new MdcTrackList(m_flags.tkParTight) );
00188 }
00189
00190 //Book NTuple and Histograms
00191 if (m_flags.lHist){
00192 StatusCode sc = bookNTuple();
00193 if (!sc.isSuccess()) { m_flags.lHist = 0;}
00194 }
00195
00196 //yzhang for HelixFitter debug
00197 if(7 == m_flags.debugFlag())TrkHelixFitter::m_debug = true;
00198 return StatusCode::SUCCESS;
00199 }
|
|
|
00041 {return m_poisonHits;}
|
|
|
00041 {return m_poisonHits;}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1