Loading...
Searching...
No Matches
HoughFinder Class Reference
#include <Hough.h>
Inheritance diagram for HoughFinder:
Public Member Functions | |
| HoughFinder (const std::string &name, ISvcLocator *pSvcLocator) | |
| StatusCode | initialize () |
| StatusCode | execute () |
| StatusCode | finalize () |
| StatusCode | bookTuple () |
| StatusCode | registerTrack (RecMdcTrackCol *&trackList_tds, RecMdcHitCol *&hitList_tds) |
| int | makeHoughHitList () |
| int | separateHoughHitByPhi () |
| int | getMcHitCol () |
| int | getMcParticleCol () |
| int | fillHistogram (HoughHit *hit, TH2D *hitMap, int charge, int vote) |
| int | fillHistogram (HoughHit *hit, TH2D *hitMap, HoughTrack *trkCandi, int vote) |
| int | fillHistogram (vector< HoughHit * > &hitList, TH2D *hitMap, int charge, int vote, HoughTrack *trkCandi=NULL) |
| int | checkHot (vector< HoughTrack > &trackVector) |
| int | checkTrack (vector< HoughTrack > &trackVector) |
| void | clearTrack (vector< HoughTrack > &trackVector) |
| int | storeTrack (vector< HoughTrack > &trackVector, RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol) |
| int | storeTracks (RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector) |
| int | storeRecTracks (RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector) |
| void | printHitList () |
| int | activeUnusedCgemHitsOnly (vector< HoughHit * > &hitPntList) |
| void | solveSharedHits (vector< HoughHit * > &hitList) |
| vector< HoughTrack >::iterator | getHoughTrkIt (vector< HoughTrack > &houghTrkList, int trkId) |
| void | clearMemory () |
| int | printTdsTrack () |
Public Attributes | |
| int | Layer |
| int | m_useCgemInGlobalFit |
| int | m_clearTrack |
| int | storeFlag |
| int | printFlag |
Detailed Description
Constructor & Destructor Documentation
◆ HoughFinder()
| HoughFinder::HoughFinder | ( | const std::string & | name, |
| ISvcLocator * | pSvcLocator ) |
Definition at line 31 of file Hough.cxx.
31 :Algorithm(name, pSvcLocator)
32{
33 declareProperty("debug", m_debug = 0);
34 declareProperty("cgem", m_cgem = 1);
35 declareProperty("mcTruth", m_mcTruth = 0);
36 declareProperty("fillNTuple", m_fillNTuple= 0);
37 declareProperty("trackCharge", m_trackCharge = 0);
38 declareProperty("findPeakMethod",m_findPeakMethod=0);
39
40 //declareProperty("nStep_xy", m_nStep_xy = 1);
41 //declareProperty("xBin_xy", m_xBin_xy);
42 //declareProperty("yBin_xy", m_yBin_xy);
43 //declareProperty("xExtend_xy", m_xExtend_xy);
44 //declareProperty("yExtend_xy", m_yExtend_xy);
45 //declareProperty("nVote_xy", m_nVote_xy);
46 //declareProperty("nRMS_xy", m_nRMS_xy);
47 //declareProperty("xInclude_xy", m_xInclude_xy = 2);
48 //declareProperty("yInclude_xy", m_yInclude_xy = 2);
49
50 //declareProperty("nStep_sz", m_nStep_sz = 1);
51 //declareProperty("xBin_sz", m_xBin_sz);
52 //declareProperty("yBin_sz", m_yBin_sz);
53 //declareProperty("xExtend_sz", m_xExtend_sz);
54 //declareProperty("yExtend_sz", m_yExtend_sz);
55 //declareProperty("nVote_sz", m_nVote_sz);
56 //declareProperty("nRMS_sz", m_nRMS_sz);
57 //declareProperty("xInclude_sz", m_xInclude_sz = 2);
58 //declareProperty("yInclude_sz", m_yInclude_sz = 2);
59
60 //read raw data setup
61 declareProperty("keepBadTdc", m_keepBadTdc = 0);
62 declareProperty("dropHot", m_dropHot = 0);
63 declareProperty("keepUnmatch", m_keepUnmatch = 0);
64 declareProperty("driftTimeUpLimit", m_driftTimeUpLimit = 1500);
65 //declareProperty("",m_ = 0);
66 declareProperty("pdtFile", m_pdtFile = "pdt.table");
67
68 // tests by llwang
69 declareProperty("nBinTheta", m_nBinTheta=100 );
70 declareProperty("nBinRho", m_nBinRho=50 );
71 declareProperty("rhoRange", m_rhoRange = 0.1);
72 declareProperty("ExtPeak_theta",m_ExtPeak_theta= 4);
73 declareProperty("ExtPeak_rho", m_ExtPeak_rho= 4);
74 declareProperty("shareHitRate", m_shareHitRate = 0.7);
75
76 declareProperty("nBinTanl", m_nBinTanl = 100);
77 declareProperty("nBinDz", m_nBinDz = 100);
78 declareProperty("ExtPeak_tanl", m_ExtPeak_tanl = 1);
79 declareProperty("ExtPeak_dz", m_ExtPeak_dz = 1);
80 declareProperty("filter", m_filter= 0);
81 declareProperty("eventFile", m_evtFile= "EventList.txt");
82 declareProperty("fitFlag", m_fitFlag= 3);
84 declareProperty("checkHits", m_checkHits=1);
89 declareProperty("removeNOuterHits", m_removeNOuterHits = 0);
90
91 // --- cuts to combine tracks
92 declareProperty("drCut_combine", m_drCut=2.0);
93 declareProperty("phi0Cut_combine", m_phi0Cut=0.2);
94 declareProperty("kappaCut_combine",m_kappaCut=1.5);
95 declareProperty("dzCut_combine", m_dzCut=5.0);
96 declareProperty("tanlCut_combine", m_tanlCut=0.5);
97
98 // --- cuts for hits saved in RecMdcTrack
99 declareProperty("chi2CutHits", m_chi2CutHits=10);
100
101 m_maxFireLayer = -1;
102 m_totEvtProcessed=0;
103}
Member Function Documentation
◆ activeUnusedCgemHitsOnly()
| int HoughFinder::activeUnusedCgemHitsOnly | ( | vector< HoughHit * > & | hitPntList | ) |
Definition at line 3133 of file Hough.cxx.
3134{
3135 int nCgemLeft(0);
3136 vector<HoughHit*>::iterator iter = hitPntList.begin();
3138 {
3139 //if((*iter)->getPairHit()==NULL) continue;// FIXME
3140 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
3141 int useOld = (*iter)->getUse();
3142 int iLayer = (*iter)->getLayer();
3143 if(iLayer<3) {
3144 if(useOld==-1||useOld==0) {
3145 (*iter)->setUse(0);
3146 nCgemLeft++;
3147 }
3148 }
3149 else
3150 if(useOld==0) (*iter)->setUse(-1);
3151 }
3152
3153 return nCgemLeft;
3154}
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
◆ bookTuple()
| StatusCode HoughFinder::bookTuple | ( | ) |
Definition at line 3474 of file Hough.cxx.
3474 {
3475 MsgStream log(msgSvc(), name());
3476
3478 if ( nt1 ){
3479 ntuple_hit= nt1;
3480 } else {
3482 if(ntuple_hit){
3483 ntuple_hit->addItem("hit_run", m_hit_run);
3484 ntuple_hit->addItem("hit_event", m_hit_event);
3485
3486 ntuple_hit->addItem("hit_nhit", m_hit_nhit, 0, 10000);
3487 ntuple_hit->addItem("hit_hitID", m_hit_nhit, m_hit_hitID);
3488 ntuple_hit->addItem("hit_hitType", m_hit_nhit, m_hit_hitType);
3489 ntuple_hit->addItem("hit_layer", m_hit_nhit, m_hit_layer);
3490 ntuple_hit->addItem("hit_wire", m_hit_nhit, m_hit_wire);
3491 ntuple_hit->addItem("hit_flag", m_hit_nhit, m_hit_flag);
3492 ntuple_hit->addItem("hit_halfCircle", m_hit_nhit, m_hit_halfCircle);
3493 ntuple_hit->addItem("hit_x", m_hit_nhit, m_hit_x);
3494 ntuple_hit->addItem("hit_y", m_hit_nhit, m_hit_y);
3495 ntuple_hit->addItem("hit_z", m_hit_nhit, m_hit_z);
3496 ntuple_hit->addItem("hit_drift", m_hit_nhit, m_hit_drift);
3497
3498 ntuple_hit->addItem("mcHit_hitID", m_hit_nhit, m_mcHit_hitID);
3499 ntuple_hit->addItem("mcHit_hitType", m_hit_nhit, m_mcHit_hitType);
3500 ntuple_hit->addItem("mcHit_layer", m_hit_nhit, m_mcHit_layer);
3501 ntuple_hit->addItem("mcHit_wire", m_hit_nhit, m_mcHit_wire);
3502 ntuple_hit->addItem("mcHit_flag", m_hit_nhit, m_mcHit_flag);
3503 ntuple_hit->addItem("mcHit_halfCircle", m_hit_nhit, m_mcHit_halfCircle);
3504 ntuple_hit->addItem("mcHit_x", m_hit_nhit, m_mcHit_x);
3505 ntuple_hit->addItem("mcHit_y", m_hit_nhit, m_mcHit_y);
3506 ntuple_hit->addItem("mcHit_z", m_hit_nhit, m_mcHit_z);
3507 ntuple_hit->addItem("mcHit_drift", m_hit_nhit, m_mcHit_drift);
3508 } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
3509 return StatusCode::FAILURE;
3510 }
3511 }
3512
3514 if ( nt2 ){
3515 ntuple_track = nt2;
3516 } else {
3518 if(ntuple_track){
3519 ntuple_track->addItem("trk_run", m_trk_run);
3520 ntuple_track->addItem("trk_event", m_trk_event);
3521 ntuple_track->addItem("trk_nTrack", m_trk_nTrack);
3522 ntuple_track->addItem("trk_trackID", m_trk_trackID);
3523 ntuple_track->addItem("trk_charge", m_trk_charge);
3524 ntuple_track->addItem("trk_flag", m_trk_flag);
3525 ntuple_track->addItem("trk_angle", m_trk_angle);
3526 ntuple_track->addItem("trk_rho", m_trk_rho);
3527 ntuple_track->addItem("trk_dAngle", m_trk_dAngle);
3528 ntuple_track->addItem("trk_dRho", m_trk_dRho);
3529 ntuple_track->addItem("trk_dTanl", m_trk_dTanl);
3530 ntuple_track->addItem("trk_dDz", m_trk_dDz);
3531 ntuple_track->addItem("trk_Xc", m_trk_Xc);
3532 ntuple_track->addItem("trk_Yc", m_trk_Yc);
3533 ntuple_track->addItem("trk_R", m_trk_R);
3534 ntuple_track->addItem("trk_dr", m_trk_dr);
3535 ntuple_track->addItem("trk_phi0", m_trk_phi0);
3536 ntuple_track->addItem("trk_kappa", m_trk_kappa);
3537 ntuple_track->addItem("trk_dz", m_trk_dz);
3538 ntuple_track->addItem("trk_tanl", m_trk_tanl);
3539 ntuple_track->addItem("trk_pxy", m_trk_pxy);
3540 ntuple_track->addItem("trk_px", m_trk_px);
3541 ntuple_track->addItem("trk_py", m_trk_py);
3542 ntuple_track->addItem("trk_pz", m_trk_pz);
3543 ntuple_track->addItem("trk_p", m_trk_p);
3544 ntuple_track->addItem("trk_phi", m_trk_phi);
3545 ntuple_track->addItem("trk_theta", m_trk_theta);
3546 ntuple_track->addItem("trk_cosTheta", m_trk_cosTheta);
3547 ntuple_track->addItem("trk_vx", m_trk_vx);
3548 ntuple_track->addItem("trk_vy", m_trk_vy);
3549 ntuple_track->addItem("trk_vz", m_trk_vz);
3550 ntuple_track->addItem("trk_vr", m_trk_vr);
3551 ntuple_track->addItem("trk_chi2", m_trk_chi2);
3552 ntuple_track->addItem("trk_fiTerm", m_trk_fiTerm);
3553 ntuple_track->addItem("trk_nhit", m_trk_nhit);
3554 ntuple_track->addItem("trk_ncluster", m_trk_ncluster);
3555 ntuple_track->addItem("trk_stat", m_trk_stat);
3556 ntuple_track->addItem("trk_ndof", m_trk_ndof);
3557 ntuple_track->addItem("trk_nster", m_trk_nster);
3558 ntuple_track->addItem("trk_nlayer", m_trk_nlayer);
3559 ntuple_track->addItem("trk_firstLayer", m_trk_firstLayer);
3560 ntuple_track->addItem("trk_lastLayer", m_trk_lastLayer);
3561 ntuple_track->addItem("trk_nCgemXClusters", m_trk_nCgemXClusters);
3562 ntuple_track->addItem("trk_nCgemVClusters", m_trk_nCgemVClusters);
3563
3564 ntuple_track->addItem("trk_nHot", m_trk_nHot, 0, 10000);
3565 ntuple_track->addItem("hot_hitID", m_trk_nHot, m_hot_hitID);
3566 ntuple_track->addItem("hot_hitType", m_trk_nHot, m_hot_hitType);
3567 ntuple_track->addItem("hot_layer", m_trk_nHot, m_hot_layer);
3568 ntuple_track->addItem("hot_wire", m_trk_nHot, m_hot_wire);
3569 ntuple_track->addItem("hot_flag", m_trk_nHot, m_hot_flag);
3570 ntuple_track->addItem("hot_halfCircle", m_trk_nHot, m_hot_halfCircle);
3571 ntuple_track->addItem("hot_x", m_trk_nHot, m_hot_x);
3572 ntuple_track->addItem("hot_y", m_trk_nHot, m_hot_y);
3573 ntuple_track->addItem("hot_z", m_trk_nHot, m_hot_z);
3574 ntuple_track->addItem("hot_drift", m_trk_nHot, m_hot_drift);
3575 ntuple_track->addItem("hot_residual", m_trk_nHot, m_hot_residual);
3576
3577 ntuple_track->addItem("mcHot_hitID", m_trk_nHot, m_mcHot_hitID);
3578 ntuple_track->addItem("mcHot_hitType", m_trk_nHot, m_mcHot_hitType);
3579 ntuple_track->addItem("mcHot_layer", m_trk_nHot, m_mcHot_layer);
3580 ntuple_track->addItem("mcHot_wire", m_trk_nHot, m_mcHot_wire);
3581 ntuple_track->addItem("mcHot_flag", m_trk_nHot, m_mcHot_flag);
3582 ntuple_track->addItem("mcHot_halfCircle", m_trk_nHot, m_mcHot_halfCircle);
3583 ntuple_track->addItem("mcHot_x", m_trk_nHot, m_mcHot_x);
3584 ntuple_track->addItem("mcHot_y", m_trk_nHot, m_mcHot_y);
3585 ntuple_track->addItem("mcHot_z", m_trk_nHot, m_mcHot_z);
3586 ntuple_track->addItem("mcHot_drift", m_trk_nHot, m_mcHot_drift);
3587
3588 ntuple_track->addItem("mcTrk_trackID", m_mcTrk_trackID);
3589 ntuple_track->addItem("mcTrk_charge", m_mcTrk_charge);
3590 ntuple_track->addItem("mcTrk_flag", m_mcTrk_flag);
3591 ntuple_track->addItem("mcTrk_angle", m_mcTrk_angle);
3592 ntuple_track->addItem("mcTrk_rho", m_mcTrk_rho);
3593 ntuple_track->addItem("mcTrk_dAngle", m_mcTrk_dAngle);
3594 ntuple_track->addItem("mcTrk_dRho", m_mcTrk_dRho);
3595 ntuple_track->addItem("mcTrk_dTanl", m_mcTrk_dTanl);
3596 ntuple_track->addItem("mcTrk_dDz", m_mcTrk_dDz);
3597 ntuple_track->addItem("mcTrk_Xc", m_mcTrk_Xc);
3598 ntuple_track->addItem("mcTrk_Yc", m_mcTrk_Yc);
3599 ntuple_track->addItem("mcTrk_R", m_mcTrk_R);
3600 ntuple_track->addItem("mcTrk_dr", m_mcTrk_dr);
3601 ntuple_track->addItem("mcTrk_phi0", m_mcTrk_phi0);
3602 ntuple_track->addItem("mcTrk_kappa", m_mcTrk_kappa);
3603 ntuple_track->addItem("mcTrk_dz", m_mcTrk_dz);
3604 ntuple_track->addItem("mcTrk_tanl", m_mcTrk_tanl);
3605 ntuple_track->addItem("mcTrk_pxy", m_mcTrk_pxy);
3606 ntuple_track->addItem("mcTrk_px", m_mcTrk_px);
3607 ntuple_track->addItem("mcTrk_py", m_mcTrk_py);
3608 ntuple_track->addItem("mcTrk_pz", m_mcTrk_pz);
3609 ntuple_track->addItem("mcTrk_p", m_mcTrk_p);
3610 ntuple_track->addItem("mcTrk_phi", m_mcTrk_phi);
3611 ntuple_track->addItem("mcTrk_theta", m_mcTrk_theta);
3612 ntuple_track->addItem("mcTrk_cosTheta", m_mcTrk_cosTheta);
3613 ntuple_track->addItem("mcTrk_vx", m_mcTrk_vx);
3614 ntuple_track->addItem("mcTrk_vy", m_mcTrk_vy);
3615 ntuple_track->addItem("mcTrk_vz", m_mcTrk_vz);
3616 ntuple_track->addItem("mcTrk_vr", m_mcTrk_vr);
3617 ntuple_track->addItem("mcTrk_chi2", m_mcTrk_chi2);
3618 ntuple_track->addItem("mcTrk_fiTerm", m_mcTrk_fiTerm);
3619 ntuple_track->addItem("mcTrk_nhit", m_mcTrk_nhit);
3620 ntuple_track->addItem("mcTrk_ncluster", m_mcTrk_ncluster);
3621 ntuple_track->addItem("mcTrk_stat", m_mcTrk_stat);
3622 ntuple_track->addItem("mcTrk_ndof", m_mcTrk_ndof);
3623 ntuple_track->addItem("mcTrk_nster", m_mcTrk_nster);
3624 ntuple_track->addItem("mcTrk_nlayer", m_mcTrk_nlayer);
3625 ntuple_track->addItem("mcTrk_firstLayer", m_mcTrk_firstLayer);
3626 ntuple_track->addItem("mcTrk_lastLayer", m_mcTrk_lastLayer);
3627 ntuple_track->addItem("mcTrk_nCgemXClusters", m_mcTrk_nCgemXClusters);
3628 ntuple_track->addItem("mcTrk_nCgemVClusters", m_mcTrk_nCgemVClusters);
3629
3630 ntuple_track->addItem("mcTrk_nHot", m_mcTrk_nHot, 0, 10000);
3631 ntuple_track->addItem("mcTrkHot_hitID", m_mcTrk_nHot, m_mcTrkHot_hitID);
3632 ntuple_track->addItem("mcTrkHot_hitType", m_mcTrk_nHot, m_mcTrkHot_hitType);
3633 ntuple_track->addItem("mcTrkHot_layer", m_mcTrk_nHot, m_mcTrkHot_layer);
3634 ntuple_track->addItem("mcTrkHot_wire", m_mcTrk_nHot, m_mcTrkHot_wire);
3635 ntuple_track->addItem("mcTrkHot_flag", m_mcTrk_nHot, m_mcTrkHot_flag);
3636 ntuple_track->addItem("mcTrkHot_halfCircle", m_mcTrk_nHot, m_mcTrkHot_halfCircle);
3637 ntuple_track->addItem("mcTrkHot_x", m_mcTrk_nHot, m_mcTrkHot_x);
3638 ntuple_track->addItem("mcTrkHot_y", m_mcTrk_nHot, m_mcTrkHot_y);
3639 ntuple_track->addItem("mcTrkHot_z", m_mcTrk_nHot, m_mcTrkHot_z);
3640 ntuple_track->addItem("mcTrkHot_drift", m_mcTrk_nHot, m_mcTrkHot_drift);
3641 } else {
3642 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/track" <<endmsg;
3643 return StatusCode::FAILURE;
3644 }
3645 }
3646
3648 if ( nt3 ){
3649 ntuple_event = nt3;
3650 } else {
3652 if(ntuple_event){
3653 ntuple_event->addItem("evt_run", m_evt_run);
3654 ntuple_event->addItem("evt_event", m_evt_event);
3655 ntuple_event->addItem("evt_nXCluster", m_evt_nXCluster);
3656 ntuple_event->addItem("evt_nVCluster", m_evt_nVCluster);
3657 ntuple_event->addItem("evt_nXVCluster", m_evt_nXVCluster);
3658 ntuple_event->addItem("evt_nCgemCluster", m_evt_nCGEMCluster);
3659 ntuple_event->addItem("evt_nAxialHit", m_evt_nAxialHit);
3660 ntuple_event->addItem("evt_nStereoHit", m_evt_nStereoHit);
3661 ntuple_event->addItem("evt_nODCHit", m_evt_nODCHit);
3662 ntuple_event->addItem("evt_nHit", m_evt_nHit);
3663
3664 ntuple_event->addItem("evt_nTrack", m_evt_nTrack, 0, 100);
3665 ntuple_event->addItem("evtTrk_trackID", m_evt_nTrack, m_evtTrk_trackID);
3666 ntuple_event->addItem("evtTrk_charge", m_evt_nTrack, m_evtTrk_charge);
3667 ntuple_event->addItem("evtTrk_flag", m_evt_nTrack, m_evtTrk_flag);
3668 ntuple_event->addItem("evtTrk_angle", m_evt_nTrack, m_evtTrk_angle);
3669 ntuple_event->addItem("evtTrk_rho", m_evt_nTrack, m_evtTrk_rho);
3670 ntuple_event->addItem("evtTrk_dAngle", m_evt_nTrack, m_evtTrk_dAngle);
3671 ntuple_event->addItem("evtTrk_dRho", m_evt_nTrack, m_evtTrk_dRho);
3672 ntuple_event->addItem("evtTrk_dTanl", m_evt_nTrack, m_evtTrk_dTanl);
3673 ntuple_event->addItem("evtTrk_dDz", m_evt_nTrack, m_evtTrk_dDz);
3674 ntuple_event->addItem("evtTrk_Xc", m_evt_nTrack, m_evtTrk_Xc);
3675 ntuple_event->addItem("evtTrk_Yc", m_evt_nTrack, m_evtTrk_Yc);
3676 ntuple_event->addItem("evtTrk_R", m_evt_nTrack, m_evtTrk_R);
3677 ntuple_event->addItem("evtTrk_dr", m_evt_nTrack, m_evtTrk_dr);
3678 ntuple_event->addItem("evtTrk_phi0", m_evt_nTrack, m_evtTrk_phi0);
3679 ntuple_event->addItem("evtTrk_kappa", m_evt_nTrack, m_evtTrk_kappa);
3680 ntuple_event->addItem("evtTrk_dz", m_evt_nTrack, m_evtTrk_dz);
3681 ntuple_event->addItem("evtTrk_tanl", m_evt_nTrack, m_evtTrk_tanl);
3682 ntuple_event->addItem("evtTrk_pxy", m_evt_nTrack, m_evtTrk_pxy);
3683 ntuple_event->addItem("evtTrk_px", m_evt_nTrack, m_evtTrk_px);
3684 ntuple_event->addItem("evtTrk_py", m_evt_nTrack, m_evtTrk_py);
3685 ntuple_event->addItem("evtTrk_pz", m_evt_nTrack, m_evtTrk_pz);
3686 ntuple_event->addItem("evtTrk_p", m_evt_nTrack, m_evtTrk_p);
3687 ntuple_event->addItem("evtTrk_phi", m_evt_nTrack, m_evtTrk_phi);
3688 ntuple_event->addItem("evtTrk_theta", m_evt_nTrack, m_evtTrk_theta);
3689 ntuple_event->addItem("evtTrk_cosTheta", m_evt_nTrack, m_evtTrk_cosTheta);
3690 ntuple_event->addItem("evtTrk_vx", m_evt_nTrack, m_evtTrk_vx);
3691 ntuple_event->addItem("evtTrk_vy", m_evt_nTrack, m_evtTrk_vy);
3692 ntuple_event->addItem("evtTrk_vz", m_evt_nTrack, m_evtTrk_vz);
3693 ntuple_event->addItem("evtTrk_vr", m_evt_nTrack, m_evtTrk_vr);
3694 ntuple_event->addItem("evtTrk_chi2", m_evt_nTrack, m_evtTrk_chi2);
3695 ntuple_event->addItem("evtTrk_fiTerm", m_evt_nTrack, m_evtTrk_fiTerm);
3696 ntuple_event->addItem("evtTrk_nhit", m_evt_nTrack, m_evtTrk_nhit);
3697 ntuple_event->addItem("evtTrk_ncluster", m_evt_nTrack, m_evtTrk_ncluster);
3698 ntuple_event->addItem("evtTrk_stat", m_evt_nTrack, m_evtTrk_stat);
3699 ntuple_event->addItem("evtTrk_ndof", m_evt_nTrack, m_evtTrk_ndof);
3700 ntuple_event->addItem("evtTrk_nster", m_evt_nTrack, m_evtTrk_nster);
3701 ntuple_event->addItem("evtTrk_nlayer", m_evt_nTrack, m_evtTrk_nlayer);
3702 ntuple_event->addItem("evtTrk_firstLayer", m_evt_nTrack, m_evtTrk_firstLayer);
3703 ntuple_event->addItem("evtTrk_lastLayer", m_evt_nTrack, m_evtTrk_lastLayer);
3704 ntuple_event->addItem("evtTrk_nCgemXClusters", m_evt_nTrack, m_evtTrk_nCgemXClusters);
3705 ntuple_event->addItem("evtTrk_nCgemVClusters", m_evt_nTrack, m_evtTrk_nCgemVClusters);
3706
3707 ntuple_event->addItem("mcEvtTrk_trackID", m_evt_nTrack, m_mcEvtTrk_trackID);
3708 ntuple_event->addItem("mcEvtTrk_charge", m_evt_nTrack, m_mcEvtTrk_charge);
3709 ntuple_event->addItem("mcEvtTrk_flag", m_evt_nTrack, m_mcEvtTrk_flag);
3710 ntuple_event->addItem("mcEvtTrk_angle", m_evt_nTrack, m_mcEvtTrk_angle);
3711 ntuple_event->addItem("mcEvtTrk_rho", m_evt_nTrack, m_mcEvtTrk_rho);
3712 ntuple_event->addItem("mcEvtTrk_dAngle", m_evt_nTrack, m_mcEvtTrk_dAngle);
3713 ntuple_event->addItem("mcEvtTrk_dRho", m_evt_nTrack, m_mcEvtTrk_dRho);
3714 ntuple_event->addItem("mcEvtTrk_dTanl", m_evt_nTrack, m_mcEvtTrk_dTanl);
3715 ntuple_event->addItem("mcEvtTrk_dDz", m_evt_nTrack, m_mcEvtTrk_dDz);
3716 ntuple_event->addItem("mcEvtTrk_Xc", m_evt_nTrack, m_mcEvtTrk_Xc);
3717 ntuple_event->addItem("mcEvtTrk_Yc", m_evt_nTrack, m_mcEvtTrk_Yc);
3718 ntuple_event->addItem("mcEvtTrk_R", m_evt_nTrack, m_mcEvtTrk_R);
3719 ntuple_event->addItem("mcEvtTrk_dr", m_evt_nTrack, m_mcEvtTrk_dr);
3720 ntuple_event->addItem("mcEvtTrk_phi0", m_evt_nTrack, m_mcEvtTrk_phi0);
3721 ntuple_event->addItem("mcEvtTrk_kappa", m_evt_nTrack, m_mcEvtTrk_kappa);
3722 ntuple_event->addItem("mcEvtTrk_dz", m_evt_nTrack, m_mcEvtTrk_dz);
3723 ntuple_event->addItem("mcEvtTrk_tanl", m_evt_nTrack, m_mcEvtTrk_tanl);
3724 ntuple_event->addItem("mcEvtTrk_pxy", m_evt_nTrack, m_mcEvtTrk_pxy);
3725 ntuple_event->addItem("mcEvtTrk_px", m_evt_nTrack, m_mcEvtTrk_px);
3726 ntuple_event->addItem("mcEvtTrk_py", m_evt_nTrack, m_mcEvtTrk_py);
3727 ntuple_event->addItem("mcEvtTrk_pz", m_evt_nTrack, m_mcEvtTrk_pz);
3728 ntuple_event->addItem("mcEvtTrk_p", m_evt_nTrack, m_mcEvtTrk_p);
3729 ntuple_event->addItem("mcEvtTrk_phi", m_evt_nTrack, m_mcEvtTrk_phi);
3730 ntuple_event->addItem("mcEvtTrk_theta", m_evt_nTrack, m_mcEvtTrk_theta);
3731 ntuple_event->addItem("mcEvtTrk_cosTheta", m_evt_nTrack, m_mcEvtTrk_cosTheta);
3732 ntuple_event->addItem("mcEvtTrk_vx", m_evt_nTrack, m_mcEvtTrk_vx);
3733 ntuple_event->addItem("mcEvtTrk_vy", m_evt_nTrack, m_mcEvtTrk_vy);
3734 ntuple_event->addItem("mcEvtTrk_vz", m_evt_nTrack, m_mcEvtTrk_vz);
3735 ntuple_event->addItem("mcEvtTrk_vr", m_evt_nTrack, m_mcEvtTrk_vr);
3736 ntuple_event->addItem("mcEvtTrk_chi2", m_evt_nTrack, m_mcEvtTrk_chi2);
3737 ntuple_event->addItem("mcEvtTrk_fiTerm", m_evt_nTrack, m_mcEvtTrk_fiTerm);
3738 ntuple_event->addItem("mcEvtTrk_nhit", m_evt_nTrack, m_mcEvtTrk_nhit);
3739 ntuple_event->addItem("mcEvtTrk_ncluster", m_evt_nTrack, m_mcEvtTrk_ncluster);
3740 ntuple_event->addItem("mcEvtTrk_stat", m_evt_nTrack, m_mcEvtTrk_stat);
3741 ntuple_event->addItem("mcEvtTrk_ndof", m_evt_nTrack, m_mcEvtTrk_ndof);
3742 ntuple_event->addItem("mcEvtTrk_nster", m_evt_nTrack, m_mcEvtTrk_nster);
3743 ntuple_event->addItem("mcEvtTrk_nlayer", m_evt_nTrack, m_mcEvtTrk_nlayer);
3744 ntuple_event->addItem("mcEvtTrk_firstLayer", m_evt_nTrack, m_mcEvtTrk_firstLayer);
3745 ntuple_event->addItem("mcEvtTrk_lastLayer", m_evt_nTrack, m_mcEvtTrk_lastLayer);
3746 ntuple_event->addItem("mcEvtTrk_nCgemXClusters", m_evt_nTrack, m_mcEvtTrk_nCgemXClusters);
3747 ntuple_event->addItem("mcEvtTrk_nCgemVClusters", m_evt_nTrack, m_mcEvtTrk_nCgemVClusters);
3748 } else {
3749 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/event" <<endmsg;
3750 return StatusCode::FAILURE;
3751 }
3752 }
3753 return StatusCode::SUCCESS;
3754}
Referenced by initialize().
◆ checkHot()
| int HoughFinder::checkHot | ( | vector< HoughTrack > & | trackVector | ) |
Definition at line 1629 of file Hough.cxx.
1630{
1631 //m_debug=true;
1632 int nDelete(0);
1633 if(trackVector.size()==0)return 0;
1634 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1635
1636 vector<HoughTrack>::iterator trkIT1 = trackVector.end();
1637 //vector<HoughTrack>::iterator trkIT1 = trackVector.rbegin();
1638 bool loop = true;
1639
1640 while(loop)
1641 {
1642 trkIT1--;
1643 if(trkIT1==trackVector.begin()) loop=false;
1644 if((trkIT1)->getFlag() == 1) continue;
1645 trkIT1->dropRedundentCgemXHits();
1646 //cout<<"hit size "<<hot1->size()<<endl;
1647 int nHits = trkIT1->getVecHitPnt().size();
1648 if(nHits<3){
1649 (trkIT1)->setFlag(1);
1650 trkIT1->releaseSelHits();
1651 if(m_debug) {
1652 cout<<"too less hits ("<<nHits<<") in track "<<(trkIT1)->getTrkID()<<endl;
1653 }
1654 continue;
1655 }
1656 int nShared = trkIT1->getNHitsShared();
1657
1658 //cout<<"shareHit = "<<shareHit<<endl;
1659 if(nShared>m_shareHitRate*nHits){
1660 trkIT1->setFlag(1);
1661 trkIT1->releaseSelHits();
1662 nDelete++;
1663 if(m_debug) {
1664 cout<<"too many shared hits ("<<nShared<<"/"<<nHits
1665 <<") in track "<<(trkIT1)->getTrkID()
1666 <<endl;
1667 }
1668 }
1669 }
1670 //m_debug=false;
1671 return nDelete;
1672}
Referenced by execute().
◆ checkTrack()
| int HoughFinder::checkTrack | ( | vector< HoughTrack > & | trackVector | ) |
Definition at line 1674 of file Hough.cxx.
1675{
1676 int nDelete(0);
1677 double dDr(999);
1678 double dPhi0(999);
1679 double dKappa(999);
1680 double dDz(999);
1681 double dTanl(999);
1682
1683 //m_debug=true;
1684
1685 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1686 for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
1687 if((trkIT1)->getFlag() == 1)continue;
1688 int charge = (trkIT1)->getCharge();
1689 double dr = (trkIT1)->dr();
1690 double phi0 = (trkIT1)->phi0();
1691 double kappa = (trkIT1)->kappa();
1692 double dz = (trkIT1)->dz();
1693 double tanl = (trkIT1)->tanl();
1694 for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
1695 if((trkIT2)->getFlag() == 1)continue;
1696 bool sameTrack(false);
1697 dDr = fabs(dr - (trkIT2)->dr());
1698 dPhi0 = fabs(phi0 - (trkIT2)->phi0());
1699 if((trkIT2)->getCharge() == charge){
1700 dKappa = fabs(kappa - (trkIT2)->kappa());
1701 dTanl = fabs(tanl - (trkIT2)->tanl());
1702 }else{
1703 dKappa = fabs(fabs(kappa) - fabs((trkIT2)->kappa()));
1704 dTanl = fabs(fabs(tanl) - fabs((trkIT2)->tanl()));
1705 }
1706 double r = fabs(dz)<fabs((trkIT2)->dz())?(trkIT1)->radius():(trkIT2)->radius();
1707 double k = dz<fabs((trkIT2)->dz())?tanl:(trkIT2)->tanl();
1708 dDz = fabs(dz - (trkIT2)->dz());
1710 dDz = fabs(dDz - nCircle*(2*M_PI*r*k));
1712 if(dPhi0<m_phi0Cut&&dKappa<m_kappaCut&&dTanl<m_tanlCut) sameTrack=true;
1713 //sameTrack = sameTrack&&(dDr<m_drCut);//FIXME
1714 //sameTrack = sameTrack&&(dPhi0<m_phi0Cut);//FIXME
1715 //sameTrack = sameTrack&&(dKappa<m_kappaCut);//FIXME
1716 //sameTrack = sameTrack&&(dDz<m_dzCut);//FIXME
1717 //sameTrack = sameTrack&&(dTanl<m_tanlCut);//FIXME
1718 if(sameTrack){
1719 //delete *trkIT2;
1720 //trackVector.erase(trkIT2);
1721 //trkIT2--;
1722 (trkIT2)->setFlag(1);
1723 if(m_debug) {
1724 cout<<"dDr, dPhi0, dKappa, dDz, dTanl "
1725 <<setw(10)<<dDr
1726 <<setw(10)<<dPhi0
1727 <<setw(10)<<dKappa
1728 <<setw(10)<<dDz
1729 <<setw(10)<<dTanl
1730 <<endl;
1731 cout<<"similar track "<<trkIT2->getTrkID()<<" is removed from hough track list"<<endl;
1732 }
1733
1734 // --- associate the hits of track 2 to track 1 if within pi/2 // to-do
1735 nDelete++;
1736 }
1737 }
1738 }
1739 if(m_debug) cout<<"checkTrack(): "<<nDelete<<" similar tracks are removed from hough track list"<<endl;
1740 //m_debug=false;
1741 return nDelete;
1742}
◆ clearMemory()
| void HoughFinder::clearMemory | ( | ) |
Definition at line 3117 of file Hough.cxx.
3118{
3119 //cout<<"m_mdcHitCol.size()"<<m_mdcHitCol.size()<<endl;
3120 vector<MdcHit*>::iterator imdcHit = m_mdcHitCol.begin();
3121 for(;imdcHit != m_mdcHitCol.end();imdcHit++){
3122 delete (*imdcHit);
3123 }
3124
3125 for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin();trkIT!=m_houghTrackList.end();trkIT++){
3126 trkIT->clearMemory();
3127 //TrkRecoTrk* trkRecoTrk = trkIT->getTrkRecoTrk();
3128 //delete trkRecoTrk;
3129 }
3130}
Referenced by execute().
◆ clearTrack()
| void HoughFinder::clearTrack | ( | vector< HoughTrack > & | trackVector | ) |
◆ execute()
| StatusCode HoughFinder::execute | ( | ) |
Definition at line 295 of file Hough.cxx.
296{
297 MsgStream log(msgSvc(), name());
298 log << MSG::INFO << "HoughFinder::execute()" << endreq;
299 cout.precision(6);
300
301
302 //cout<<"line "<<__LINE__<<endl;
303 // --- Get event header
304 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
305 if(!eventHeader){
306 log << MSG::ERROR << "Could not find Event Header" << endreq;
307 return StatusCode::FAILURE;
308 }
309 m_run = eventHeader->runNumber();
310 m_event = eventHeader->eventNumber();
311 m_totEvtProcessed++; //if(m_totEvtProcessed%100==0) cout<<"HoughFinder: "<<m_totEvtProcessed<<" events processed."<<endl;
312
313 //cout<<"line "<<__LINE__<<endl;
314 // --- event start time
315 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
316 if(!evTimeCol){
317 log << MSG::ERROR << "Could not retrieve RecEsTimeCol" << endreq;
318 return StatusCode::FAILURE;
319 }
320 m_bunchT0 = 0;
321 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
322 if (iter_evt != evTimeCol->end()){
323 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
324 //cout<<"this event T0 = "<<m_bunchT0<<endl;
325 }
326
327 // --- clear, register and store reconstructed track
328 registerTrack(m_trackList_tds,m_hitList_tds);
329
330
331 //cout<<"line "<<__LINE__<<endl;
332 if(m_filter){
333 ifstream lowPt_Evt;
334 lowPt_Evt.open(m_evtFile.c_str());
335 vector<int> evtlist;
336 int evtNum;
337 while( lowPt_Evt >> evtNum) {
338 evtlist.push_back(evtNum);
339 }
340 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),m_event);
341 if(iter_evt == evtlist.end()){
342 setFilterPassed(false);
343 return StatusCode::SUCCESS;
344 }
345 else{
346 setFilterPassed(true);
347 }
348 }
349 //if(m_event<2000||m_event>3000){
350 //setFilterPassed(false);
351 //return StatusCode::SUCCESS;
352 //}
353
354 if(m_debug)
355 {
356 //cout<<endl;
357 cout<<"===================================================================================================="<<endl;
358 cout<<"run:"<<m_run<<" , event: "<<m_event<<endl;
359 //cout<<m_event<<" , ";
360 //cout<<endl;
361 }
362
363 //cout<<"line "<<__LINE__<<endl;
364 // --- prepare hits for reconstruction
366
367 //cout<<"line "<<__LINE__<<endl;
368 int nMcHit;
369 if(m_run<0&&m_mcTruth){
371 //cout<<"nMcHit "<<nMcHit<<endl;
372 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
373 getMcParticleCol();
374 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
378 //halfCircle(m_mcTrackCol,m_mcHitCol);
379 }
380
381 //cout<<"line "<<__LINE__<<endl;
382 searchCircle();
383 //cout<<"ntracks:"<<m_houghTrackList.size()<<endl;
384 //cout<<"line "<<__LINE__<<endl;
386 //cout<<"line "<<__LINE__<<endl;
388 //checkTrack(m_houghTrackList);
389 //cout<<"line "<<__LINE__<<endl;
390 associateVHits();
391 //cout<<"line "<<__LINE__<<endl;
392 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
394 else
396 else
398 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
402 //cout<<"line "<<__LINE__<<endl;
403
404 if(m_fillNTuple)dumpHit();
405 if(m_fillNTuple==1)dumpHoughTrack();
406 if(m_fillNTuple==1)dumpHoughEvent();
407 if(m_fillNTuple==2)dumpTdsTrack();
408 if(m_fillNTuple==2)dumpTdsEvent();
409 //cout<<"line "<<__LINE__<<endl;
410
411 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
413 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
417
418 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
420 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
424
425 clearTrack(m_houghTrackList);
426 //cout<<"line "<<__LINE__<<endl;
427
428 clearMemory();
429 //cout<<"line "<<__LINE__<<endl;
430 if(m_debug)cout<<endl;
431
432 return StatusCode::SUCCESS;
433}// <-- execute()
int storeRecTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:1896
void solveSharedHits(vector< HoughHit * > &hitList)
Definition Hough.cxx:3157
int storeTrack(vector< HoughTrack > &trackVector, RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition Hough.cxx:1855
int storeTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:2181
StatusCode registerTrack(RecMdcTrackCol *&trackList_tds, RecMdcHitCol *&hitList_tds)
Definition Hough.cxx:1802
◆ fillHistogram() [1/3]
| int HoughFinder::fillHistogram | ( | HoughHit * | hit, |
| TH2D * | hitMap, | ||
| HoughTrack * | trkCandi, | ||
| int | vote ) |
Definition at line 1023 of file Hough.cxx.
1024{
1026 int nBin = hitMap->GetXaxis()->GetNbins();
1027 double xMin = hitMap->GetXaxis()->GetXmin();
1028 double xMax = hitMap->GetXaxis()->GetXmax();
1029 double yMin = hitMap->GetYaxis()->GetXmin();
1030 double yMax = hitMap->GetYaxis()->GetXmax();
1031 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;// step size
1033
1035 int nPoint(0);
1036 if(flag==0) // X-view hits
1037 {
1039 if(!firstHalf) return 0;
1043 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
1044
1045 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
1047 double Y1 = 2*yHit/denominator;
1049 double Y2 = 2*yHit/denominator;
1051
1052 while(x<xMax){
1057 double hitOnCircle1 = charge*slope1*y1;
1058 double hitOnCircle2 = charge*slope2*y2;
1059
1061 cut = yMin<y1 && y1<yMax;
1062 //cut = cut && hitOnCircle1 <=0;
1064 hitMap->Fill(x,y1);
1065 nPoint++;
1066 }
1067
1068 cut = yMin<y2 && y2<yMax;
1069 //cut = cut && hitOnCircle2 <= 0;
1071 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
1072 //int bin = hitMap->FindBin(x,y2);
1073 //if(hitMap->GetBinContent(bin)>0)continue;
1074 hitMap->Fill(x,y2);
1075 nPoint++;
1076 }
1077 x += dx;
1078 }
1079 }else{ // V-view hits
1080 //if(hit.getPairHit()==NULL)return nPoint;// FIXME
1081 //if(hit.getPairHit()->getHalfCircle()!=1)return nPoint;// FIXME
1083 //cout<<nTangency<<endl;
1084 if(nTangency==0){
1085 //hit.setUse(-1);// why? FIXME
1086 }
1087 else
1088 {
1089 vector<HoughHit::S_Z> sz = hit->getSZ();
1090 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1092 {
1095 while(x<xMax)
1096 {
1098 //y = S*cos(x) + Z*sin(x);
1100 cut = yMin<y && y<yMax;
1102 {
1103 hitMap->Fill(x,y);
1104 nPoint++;
1105 }
1106 x += dx;
1107 }
1108 }
1109 }
1110 }
1111 return nPoint;
1112}
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition KarFin.h:27
#define X1
#define X2
Definition DotsConnection.h:26
◆ fillHistogram() [2/3]
| int HoughFinder::fillHistogram | ( | HoughHit * | hit, |
| TH2D * | hitMap, | ||
| int | charge, | ||
| int | vote ) |
Definition at line 947 of file Hough.cxx.
948{
949 int nBin = hitMap->GetXaxis()->GetNbins();
950 double xMin = hitMap->GetXaxis()->GetXmin();
951 double xMax = hitMap->GetXaxis()->GetXmax();
952 double yMin = hitMap->GetYaxis()->GetXmin();
953 double yMax = hitMap->GetYaxis()->GetXmax();
954 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
956
958 int nPoint(0);
959 if(flag==0){
963 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
964
965 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
967 double Y1 = 2*yHit/denominator;
969 double Y2 = 2*yHit/denominator;
971
972 while(x<xMax){
977 double hitOnCircle1 = charge*slope1*y1;
978 double hitOnCircle2 = charge*slope2*y2;
979
981 cut = yMin<y1 && y1<yMax;
984 hitMap->Fill(x,y1);
985 nPoint++;
986 }
987
988 cut = yMin<y2 && y2<yMax;
991 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
992 //int bin = hitMap->FindBin(x,y2);
993 //if(hitMap->GetBinContent(bin)>0)continue;
994 hitMap->Fill(x,y2);
995 nPoint++;
996 }
997 x += dx;
998 }
999 }else{
1000 vector<HoughHit::S_Z> sz = hit->getSZ();
1001 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1005 double dy = -S*dx;
1006 double y = -dy + Z;
1007 while(x<xMax){
1008 x += dx;
1009 y += dy;
1010 //y = S*cos(x) + Z*sin(x);
1012 cut = yMin<y && y<yMax;
1014 hitMap->Fill(x,y);
1015 nPoint++;
1016 }
1017 }
1018 }
1019 }
1020 return nPoint;
1021}
Referenced by fillHistogram().
◆ fillHistogram() [3/3]
| int HoughFinder::fillHistogram | ( | vector< HoughHit * > & | hitList, |
| TH2D * | hitMap, | ||
| int | charge, | ||
| int | vote, | ||
| HoughTrack * | trkCandi = NULL ) |
Definition at line 1114 of file Hough.cxx.
1115{
1116 int nHitsFilled = 0;
1117 m_VHoughHitListOnSZmap.clear();
1118 std::vector<HoughHit*>::iterator iter = hitList.begin();
1120 {
1121 int used = (*iter)->getUse();
1122 //(*iter)->print();
1123 if(used==0) {
1124 int nPoint = 0;
1126 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1127 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1129 }
1130 else{
1131 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1132 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1134 }
1135 if(nPoint>0) {
1136 nHitsFilled++;
1137 m_VHoughHitListOnSZmap.push_back(*iter);
1138 }
1139 }
1140 }
1141 return nHitsFilled;
1142}
int fillHistogram(HoughHit *hit, TH2D *hitMap, int charge, int vote)
Definition Hough.cxx:947
◆ finalize()
| StatusCode HoughFinder::finalize | ( | ) |
◆ getHoughTrkIt()
| vector< HoughTrack >::iterator HoughFinder::getHoughTrkIt | ( | vector< HoughTrack > & | houghTrkList, |
| int | trkId ) |
Definition at line 3202 of file Hough.cxx.
3203{
3204 vector<HoughTrack>::iterator it = houghTrkList.begin();
3205 for(; it!=houghTrkList.end(); it++)
3206 {
3207 if(it->getTrkID()==trkId) break;
3208 }
3209 return it;
3210}
Referenced by solveSharedHits().
◆ getMcHitCol()
| int HoughFinder::getMcHitCol | ( | ) |
— put match RecCgemCluster into CemgMcHit
Definition at line 514 of file Hough.cxx.
515{
516 MsgStream log(msgSvc(), name());
517 int nMcHit(0);
518 if(!(m_mcHitCol.empty()))m_mcHitCol.clear();
519
520 // --- get cgemMcHit
521 if(m_cgem){
522 SmartDataPtr<CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
523 if (!cgemMcHitCol){
524 log << MSG::ERROR << "Could not retrieve CgemMcHitCol" << endreq;
525 }else{
526 CgemMcHitCol::iterator cgemMcHitIt = cgemMcHitCol->begin();
527 for(;cgemMcHitIt != cgemMcHitCol->end();cgemMcHitIt++,nMcHit++){
529 HoughHit hit(cgemMcHit,m_bunchT0,nMcHit);
530 m_mcHitCol.push_back(hit);
531 HoughHit* phit = &(*m_mcHitCol.rbegin());
532
533 // --- match fired strip ID between CemgMcHit and RecCgemCluster
534 vector<int> xFireStripID;
535 vector<int> vFireStripID;
536 int sheet;
538 for(int ii = 0; ii < identifier.GetSize(); ii++){
539 //const Identifier ident = Identifier(identifier.GetAt(ii));
543 sheet = CgemID::sheet(ident);
544 //cout
545 //<< " hitID:" << phit->getHitID()
546 //<< " layerID=" << CgemID::layer(ident)
547 //<< " sheetID=" << CgemID::sheet(ident)
548 //<< " isXstrip=" << CgemID::is_xstrip(ident)
549 //<< " stripID=" << CgemID::strip(ident)
550 //<<endl;
551 }
552 sort(xFireStripID.begin(),xFireStripID.end());
553 sort(vFireStripID.begin(),vFireStripID.end());
554 //vector<int>::iterator uniqueX = unique(xFireStripID.begin(),xFireStrip ID.end());
555 //vector<int>::iterator uniqueV = unique(vFireStripID.begin(),vFireStrip ID.end());
556 //xFireStripID.erase(uniqueX,xFireStripID.end());
557 //vFireStripID.erase(uniqueV,vFireStripID.end());
558 //for(vector<int>::iterator i = vFireStripID.begin();i!=vFireStripID.end();i++)cout<<*i<<endl;
559 //cout<<endl;
560
561 bool findX(false), findV(false);
569 //if(recCgemCluster->getsheetid() != sheet)continue;
570 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
571 if(!findX){
572 findX = recCgemCluster->getflag() == 0
573 && xFireStripID.size() > 0
574 && recCgemCluster->getclusterflagb() == *(xFireStripID.begin())
575 && recCgemCluster->getclusterflage() == *(xFireStripID.rbegin());
576 if(findX){
577 xCluster = recCgemCluster;
578 hitIt->setPairHit(phit);
579 }
580 }
581 if(!findV){
582 findV = recCgemCluster->getflag() == 1
583 && vFireStripID.size() > 0
584 && recCgemCluster->getclusterflagb() == *(vFireStripID.begin())
585 && recCgemCluster->getclusterflage() == *(vFireStripID.rbegin());
586 if(findV){
587 vCluster = recCgemCluster;
588 hitIt->setPairHit(phit);
589 }
590 }
591
592 }
593
594 /// --- put match RecCgemCluster into CemgMcHit
595 if(findX&&findV){
601 //if(recCgemCluster->getsheetid() != sheet)continue;
602 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
605 phit->setCgemCluster(recCgemCluster);
608 phit->setPairHit(&(*hitIt));
609 hitIt->setPairHit(phit);
610 }
611 }
612 }
613 else if(findX&&!findV){
614 phit->setCgemCluster(xCluster);
617 }
618 else if(!findX&&findV){
619 phit->setCgemCluster(vCluster);
622 }
623 }
624 }
625 }
626
627 // --- get mdcMcHit
628 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
629 if(!mdcMcHitCol){
630 log << MSG::ERROR << "Could not retrieve MdcMcHitCol" << endreq;
631 }else{
632 MdcMcHitCol::iterator mdcMcHitIt = mdcMcHitCol->begin();
633 for(;mdcMcHitIt != mdcMcHitCol->end();mdcMcHitIt++,nMcHit++){
635 HoughHit hit(mdcMcHit,m_bunchT0,nMcHit);
636 m_mcHitCol.push_back(hit);
637 HoughHit* phit = &(*m_mcHitCol.rbegin());
642 phit->setDigi(mdcDigi);
643 phit->setPairHit(&(*hitIt));
644 hitIt->setPairHit(phit);
645 }
646 }
647 }
648 }
649
652 if(hitIt->getLayer()!=mcHitIt->getLayer())continue;
654 if(mcHitIt->getCgemCluster()==NULL)continue;
655 int recClusterID(-1);
656 int mcClusterID(-2);
657 int recFlag = hitIt->getCgemCluster()->getflag();
658 int mcFlag = mcHitIt->getCgemCluster()->getflag();;
659 if(recFlag==mcFlag){
660 recClusterID = hitIt->getCgemCluster()->getclusterid();
661 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
662 }
663
664 if(recFlag==1&&mcFlag==2){
665 recClusterID = hitIt->getCgemCluster()->getclusterid();
666 mcClusterID = mcHitIt->getCgemCluster()->getclusterflage();
667 }
668
669 if(recFlag==0&&mcFlag==2){
670 recClusterID = hitIt->getCgemCluster()->getclusterid();
671 mcClusterID = mcHitIt->getCgemCluster()->getclusterflagb();
672 }
673
674 if(recFlag==2&&mcFlag==1){
675 recClusterID = hitIt->getCgemCluster()->getclusterflage();
676 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
677 }
678
679 if(recFlag==2&&mcFlag==0){
680 recClusterID = hitIt->getCgemCluster()->getclusterflagb();
681 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
682 }
683
684 if(recClusterID==mcClusterID){
685 hitIt->setPairHit(&(*mcHitIt));
686 hitIt->setCgemMcHit(mcHitIt->getCgemMcHit());
687 //mcHitIt->addPairHit(&(*hitIt));
688
689 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" "<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
690 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
691 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" ";
692 //cout<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
693 //cout<<recFlag<<" "<<mcFlag<<" ";
694 //cout<<recClusterID<<" "<<mcClusterID<<endl;
695
696 //cout<<hitIt->getHitPosition();
697 //cout<<endl;
698 //cout<<mcHitIt->getHitPosition();
699 //cout<<endl;
700 //cout<<endl;
701 }
702 }
703 else if(hitIt->getHitType() == HoughHit::MDCHIT && mcHitIt->getHitType() == HoughHit::MDCMCHIT){
704 if(mcHitIt->getDigi()==NULL)continue;
705 //if(hitIt->getLayer() == mcHitIt->getLayer() && hitIt->getWire() == mcHitIt->getWire())
706 if(hitIt->getDigi()->identify() == mcHitIt->getDigi()->identify())
707 {
708 hitIt->setPairHit(&(*mcHitIt));
709 hitIt->setMdcMcHit(mcHitIt->getMdcMcHit());
710 //mcHitIt->addPairHit(&(*hitIt));
711 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<endl;
712 //cout<<mcHitIt->getLayer()<<" "<<mcHitIt->getWire()<<endl;
713 //cout<<endl;
714 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
715 }
716 }
717 }
718 }
719 //for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
720 //if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
721 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<" ";
722 //if(hitIt->getPairHit()!=NULL){
723 //cout<<hitIt->getPairHit()->getHitID()<<" ";
724 //cout<<hitIt->getPairHit()->getLayer()<<" "<<hitIt->getPairHit()->getWire()<<" ";
725 //}
726 //cout<<endl;
727 //}
728
729 return nMcHit;
730}
Definition CgemMcHit.h:21
Definition MdcMcHit.h:16
Definition HoughHit.h:25
void setCgemCluster(const RecCgemCluster *cgemCluster)
Definition HoughHit.h:75
Definition Identifier.h:20
Definition MdcDigi.h:15
Definition RecCgemCluster.h:16
Referenced by execute().
◆ getMcParticleCol()
| int HoughFinder::getMcParticleCol | ( | ) |
Definition at line 732 of file Hough.cxx.
732 {
733 const int maxCell[43] = {40,44,48,56, 64,72,80,80,
734 76,76,88,88, 100,100,112,112, 128,128,140,140,
735 160,160,160,160, 176,176,176,176, 208,208,208,208, 240,240,240,240,
736 256,256,256,256, 288,288,288 };
737 m_mcTrackCol.clear();
738 MsgStream log(msgSvc(), name());
739 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
740 if (!mcParticleCol) {
741 log << MSG::ERROR << "Could not find McParticle" << endreq;
742 }else{
743 bool psipDecay(false);
744 McParticleCol::iterator iter_mc = mcParticleCol->begin();
745 for (;iter_mc != mcParticleCol->end(); iter_mc++){
746 int trackIndex = (*iter_mc)->trackIndex();
747 int pid = (*iter_mc)->particleProperty();
748 bool primaryParticle = (*iter_mc)->primaryParticle();
749 bool leafParticle = (*iter_mc)->leafParticle();
750 bool decayFromGenerator = (*iter_mc)->decayFromGenerator();
751 bool decayInFlight = (*iter_mc)->decayInFlight();
752 HepLorentzVector initialPosition = (*iter_mc)->initialPosition();
753 HepLorentzVector initialFourMomentum = (*iter_mc)->initialFourMomentum();
754 //HepLorentzVector finalPosition = (*iter_mc)->finalPosition();
755 //HepLorentzVector finalFourMomentum = (*iter_mc)->finalFourMomentum();
756 unsigned int statusFlags = (*iter_mc)->statusFlags();
757 //string process = (*iter_mc)->getProcess();
758 //if(primaryParticle)continue;
759 //if(!decayFromGenerator)continue;
760 //if(pid==100443)psipDecay = true;
761 //if(!psipDecay) continue;
762 //if(!(fabs(pid)==11||fabs(pid)==211))continue;
763
764 string name;
765 int charge(0);
766 if( pid == 0 ) {
767 //cout << "Wrong particle id" <<endl;
768 continue;
769 }else{
770 IPartPropSvc* p_PartPropSvc;
771 static const bool CREATEIFNOTTHERE(true);
772 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
773 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
774 log << MSG::ERROR << "Could not initialize Particle Properties Service" << endreq;
775 }
776 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
777 if( p_particleTable->particle(pid) ){
778 name = p_particleTable->particle(pid)->name();
779 charge = p_particleTable->particle(pid)->charge();
780 }else if( p_particleTable->particle(-pid) ){
781 name = "anti " + p_particleTable->particle(-pid)->name();
782 charge = (-1)*p_particleTable->particle(-pid)->charge();
783 }
784 }
785
786 //cout
787 //<<setw(8)<<trackIndex
788 //<<setw(12)<<name
789 //<<setw(8)<<pid
790 //<<setw(8)<<primaryParticle
791 //<<setw(8)<<leafParticle
792 //<<setw(8)<<decayFromGenerator
793 //<<setw(8)<<decayInFlight
794 //<<setw(8)<<statusFlags
795 //<<endl;
796
797 HoughTrack track(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
798 //HoughTrack* track = new HoughTrack(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
799 HepPoint3D pivot(0,0,0);
800 track.pivot(pivot);
801
802 int nHot(0);
803 vector<HoughHit*> hotList;
804 vector<HoughHit>& hitList = m_mcHitCol;
806 vector<int> trkID = hitIt->getTrkID();
807 if(trkID[0]==track.getTrkID()){
809 track.addHitPnt(&(*hitIt));
810 track.addVecStereoHitPnt(&(*hitIt));
811 hotList.push_back(&(*hitIt));
812 }
814 if(hitIt->getFlag()==0)track.addHitPnt(&(*hitIt));
815 else track.addVecStereoHitPnt(&(*hitIt));
816 hotList.push_back(&(*hitIt));
817 }
818 nHot++;
819 }
820 }
821 if(nHot>0)m_mcTrackCol.push_back(track);
822
823 ///*
824 bool classifyHotByHotLayer(false);
825 bool classifyHotByTrackCenter(true);
826 if(classifyHotByHotLayer){
827 int lastHitLayer(0);
828 int deltaLayer(0);
829 int halfCircle(0);
830 HoughHit* lastLayerHit;
831 vector<HoughHit*> hitOnLayer;
832 int nMdcHot(0);
833 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
835 if(nMdcHot==0){//the first hit on track
836 hitOnLayer.clear();
837 hitOnLayer.push_back(*hitIt);
838 lastLayerHit = *hitIt;
839 lastHitLayer = (*hitIt)->getLayer();
840 deltaLayer = 0;
841 halfCircle = 1;
842 }else{
843 if((*hitIt)->getLayer()==lastHitLayer){//find hits on the same layer
844 hitOnLayer.push_back(*hitIt);
845 lastHitLayer = (*hitIt)->getLayer();
846 }else{//another layer hit
847 if(deltaLayer*((*hitIt)->getLayer()-lastHitLayer)>=0){//not the turn layer
849 (*iter)->setHalfCircle(halfCircle);
850 lastLayerHit = *iter;
851 }
852 }else{//the turn layer
853 halfCircle++;
854 int maxWireGap(0);
855 int turnPoint(0);
856 int deltaWire(0);
858 deltaWire = fabs((*iter)->getWire()-lastLayerHit->getWire());
859 if(deltaWire>(maxCell[(*iter)->getLayer()])/2) deltaWire = maxCell[(*iter)->getLayer()] - deltaWire;
860 if(deltaWire>maxWireGap){
861 maxWireGap = deltaWire;
862 turnPoint = iter - hitOnLayer.begin();
863 }
864 lastLayerHit = *iter;
865 }
866 deltaWire = fabs((*hitIt)->getWire()-lastLayerHit->getWire());
867 if(deltaWire>(maxCell[(*hitIt)->getLayer()])/2) deltaWire = maxCell[(*hitIt)->getLayer()] - deltaWire;
868 if(deltaWire>maxWireGap){
869 maxWireGap = deltaWire;
870 turnPoint = hitOnLayer.size();
871 }
872
873 if(maxWireGap>2){
876 if(shift<turnPoint)(*iter)->setHalfCircle(halfCircle-1);
877 else (*iter)->setHalfCircle(halfCircle);
878 }
879 }else{
882 if(shift<hitOnLayer.size()/2)(*iter)->setHalfCircle(halfCircle-1);
883 else (*iter)->setHalfCircle(halfCircle);
884 }
885 }
886 }
887 hitOnLayer.clear();
888 hitOnLayer.push_back(*hitIt);
889 lastHitLayer = (*hitIt)->getLayer();
890 deltaLayer = (*hitIt)->getLayer()-lastHitLayer;
891 }
892
893 }
894 nMdcHot++;
895 }
896 }
897 //*/
898 ///*
899 if(classifyHotByTrackCenter){
900 int cgemHalf(0), mdcHalf(0);
901 int cgemSign(0), mdcSign(0);
902 int half(0), sign(0);
903 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
904 //(*hitIt)->print();
906 if(track.judgeHalf(*hitIt) != cgemSign){
907 cgemHalf++;
908 cgemSign = track.judgeHalf(*hitIt);
909 half = cgemHalf;
910 sign = cgemSign;
911 }
912 }
914 if(track.judgeHalf(*hitIt) != mdcSign){
915 mdcHalf++;
916 mdcSign = track.judgeHalf(*hitIt);
917 half = mdcHalf;
918 sign = mdcSign;
919 }
920 }
921 (*hitIt)->setHalfCircle(half);
922 }
923 }
924 if(hotList.size()>0){
926 //track.printHot();
927 //cout<<endl;
928 }
929 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++){
930 HoughHit* hitIt = *hotIt;
932 //cout<<hitIt->getHitID();
933 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
934 //cout<<endl;
935 }
937 }
938 sort(m_mcTrackCol.begin(),m_mcTrackCol.end(),largerPt);
939 //cout<<"N MCtrack: "<<m_mcTrackCol.size()<<endl;
940 for(vector<HoughTrack>::iterator trkIter=m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
941 //trkIter->print();
942 }
943 }
944 return m_mcTrackCol.size();
945}
Definition HoughTrack.h:23
Referenced by execute().
◆ initialize()
| StatusCode HoughFinder::initialize | ( | ) |
Definition at line 105 of file Hough.cxx.
106{
107 MsgStream log(msgSvc(), name());
108 log << MSG::INFO << "HoughFinder::initialize()" << endreq;
109 StatusCode sc;
110
111 // --- RawData
112 IRawDataProviderSvc* irawDataProviderSvc;
113 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
115 if ( sc.isFailure() ){
116 log << MSG::ERROR << name()<<" Could not load RawDataProviderSvc!" << endreq;
117 return StatusCode::FAILURE;
118 }
119
120
121 // --- MDC Geometry
122 IMdcGeomSvc* imdcGeomSvc;
123 sc = service ("MdcGeomSvc", imdcGeomSvc);
125 if ( sc.isFailure() ){
126 log << MSG::ERROR << "Could not load MdcGeoSvc!" << endreq;
127 return StatusCode::FAILURE;
128 }
129
130 // --- MdcCalibFunSvc
131 IMdcCalibFunSvc* imdcCalibSvc;
132 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
134 if ( sc.isFailure() ){
135 log << MSG::ERROR << "Could not load MdcCalibFunSvc!" << endreq;
136 return StatusCode::FAILURE;
137 }
138
139 if(m_cgem){
140 // --- CGEM geometry
141 ICgemGeomSvc* icgemGeomSvc;
142 sc = service ("CgemGeomSvc", icgemGeomSvc);
144 if ( sc.isFailure() ){
145 log << MSG::ERROR << "Could not load CgemGeomSvc!" << endreq;
146 return StatusCode::FAILURE;
147 }
148
149 // --- CgemCalibFunSvc
150 ICgemCalibFunSvc* icgemCalibSvc;
151 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
153 if ( sc.isFailure() ){
154 log << MSG::ERROR << "Could not load CgemCalibFunSvc!" << endreq;
155 return StatusCode::FAILURE;
156 }
157
158 }
159
160 // --- MagneticFieldSvc
161 sc = service ("MagneticFieldSvc",m_pIMF);
162 if( sc.isFailure() ) {
163 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
164 return StatusCode::FAILURE;
165 }
167 //log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
169
171 //read pdt
172 Pdt::readMCppTable(m_pdtFile);
173
174 // --- Mdc Print Service
175 //IMdcPrintSvc* imdcPrintSvc;
176 //sc = service ("MdcPrintSvc", imdcPrintSvc);
177 //m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
178 //if ( sc.isFailure() ){
179 //log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
180 //return StatusCode::FAILURE;
181 //}
182
183 // --- Bes Timer Service
184 //sc = service( "BesTimerSvc", m_timersvc);
185 //if( sc.isFailure() ) {
186 //log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
187 //return StatusCode::FAILURE;
188 //}
189 //m_timer_all = m_timersvc->addItem("Execution");
190 //m_timer_all->propName("nExecution");
191
192 HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
193 HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
194 HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
195 HoughHit::setCgemCalibFunSvc(m_cgemCalibFunSvc);
196 HoughHit::setMdcDetector(m_mdcDetector);
199
200 m_roughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",m_nBinTheta,0.,M_PI, m_nBinRho, -1.*m_rhoRange, m_rhoRange);
201 //m_tanlRange = 0.93/sqrt(1-0.93*0.93);//FIXME
202 m_tanlRange = 3.0;//FIXME
203 m_dzRange = 50;
204 m_roughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",m_nBinTanl,-1.*m_tanlRange,m_tanlRange,m_nBinDz,-1.*m_dzRange,m_dzRange);
205
206 const int nPt=12;
207 double rho[nPt] = {0.07, 0.056,0.045,0.038,0.0335,0.03,0.0198,0.0148,0.0074,0.00376,0.00303,0.00157};
208 double cut1_Cgem_99[12]={-2.14,-1.36, -0.6 , -0.46, -0.35, -0.59, -0.32, -0.26, -0.22, -0.21, -0.21, -0.211};
209 double cut2_Cgem_99[12]={ 0.5 , 1.77, 1.99, 1.94, 1.99, 1.9 , 0.31, 0.27, 0.24, 0.23, 0.23, 0.222};
210 double cut1_ODC1_99[12]={-1.28,-0.71, -0.58, -0.62, -0.64, -0.68, -0.18, -0.14, -0.11, -0.1 , -0.1 , -0.11 };
211 double cut2_ODC1_99[12]={ 0.9 , 0.74, 0.42, 0.37, 0.32, 0.37, 0.28, 0.25, 0.24, 0.24, 0.24, 0.23};
212 double cut1_ODC2_99[12]={ -2.14, -1.25, -1.28, -0.86, -0.47, -0.78, -0.36, -0.44, -0.61, -0.67, -0.64, -0.82 };
213 double cut2_ODC2_99[12]={ -1.35, 0.55 , 0.53 , 0.83 , 0.85 , 0.83 , 0.38 , 0.55 , 0.49 , 0.46 , 0.49 , 0.4};
214 m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
215 m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
216 m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
217 m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
218 m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
219 m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
220
221
223 double momenta[100] = {0};
224 double cuts[2][43][100] = {0};
225 string cutFilePath = getenv("HOUGHTRANSALGROOT");
226 cutFilePath += "/share/cut.txt";
227 ifstream fin(cutFilePath.c_str(), ios::in);
228 //cout<<"cutfiledir:"<<cutFilePath.c_str()<<endl;
229 //ifstream fin("/workfs/bes/huangzhen/workarea/CgemBoss/6.6.5.c/Reconstruction/HoughFinder/HoughFinder-00-00-13/share/cut.txt",ios::in);
230 if (!fin.good())
231 cout << "Error in " << __FILE__<<" when opening "<<cutFilePath.c_str()<<endl;
232 bool firstline(true);
233 string strcom;
234 while(std::getline(fin, strcom)){
235 if(firstline){
236 fin >> npt
237 >> momenta[0] >> momenta[1] >> momenta[2] >> momenta[3] >> momenta[4] >> momenta[5] >> momenta[6] >> momenta[7] >> momenta[8] >> momenta[9]
238 >> momenta[10] >> momenta[11] >> momenta[12] >> momenta[13] >> momenta[14] >> momenta[15] >> momenta[16] >> momenta[17] >> momenta[18] >> momenta[19]
239 >> momenta[20] >> momenta[21] >> momenta[22] >> momenta[23] >> momenta[24] >> momenta[25] >> momenta[26] >> momenta[27] >> momenta[28] >> momenta[29]
240 >> momenta[30] >> momenta[31] >> momenta[32] >> momenta[33] >> momenta[34] >> momenta[35] >> momenta[36] >> momenta[37] >> momenta[38] >> momenta[39]
241 >> momenta[40] >> momenta[41] >> momenta[42] >> momenta[43] >> momenta[44] >> momenta[45] >> momenta[46] >> momenta[47] >> momenta[48] >> momenta[49]
242 >> momenta[50] >> momenta[51] >> momenta[52] >> momenta[53] >> momenta[54] >> momenta[55] >> momenta[56] >> momenta[57] >> momenta[58] >> momenta[59]
243 >> momenta[60] >> momenta[61] >> momenta[62] >> momenta[63] >> momenta[64] >> momenta[65] >> momenta[66] >> momenta[67] >> momenta[68] >> momenta[69]
244 >> momenta[70] >> momenta[71] >> momenta[72] >> momenta[73] >> momenta[74] >> momenta[75] >> momenta[76] >> momenta[77] >> momenta[78] >> momenta[79]
245 >> momenta[80] >> momenta[81] >> momenta[82] >> momenta[83] >> momenta[84] >> momenta[85] >> momenta[86] >> momenta[87] >> momenta[88] >> momenta[89]
246 >> momenta[90] >> momenta[91] >> momenta[92] >> momenta[93] >> momenta[94] >> momenta[95] >> momenta[96] >> momenta[97] >> momenta[98] >> momenta[99]
247 ;
248 firstline = false;
249 continue;
250 }
251 fin >> s >> l;
252 fin >> cuts[s][l][0] >> cuts[s][l][1] >> cuts[s][l][2] >> cuts[s][l][3] >> cuts[s][l][4] >> cuts[s][l][5] >> cuts[s][l][6] >> cuts[s][l][7] >> cuts[s][l][8] >> cuts[s][l][9]
253 >> cuts[s][l][10] >> cuts[s][l][11] >> cuts[s][l][12] >> cuts[s][l][13] >> cuts[s][l][14] >> cuts[s][l][15] >> cuts[s][l][16] >> cuts[s][l][17] >> cuts[s][l][18] >> cuts[s][l][19]
254 >> cuts[s][l][20] >> cuts[s][l][21] >> cuts[s][l][22] >> cuts[s][l][23] >> cuts[s][l][24] >> cuts[s][l][25] >> cuts[s][l][26] >> cuts[s][l][27] >> cuts[s][l][28] >> cuts[s][l][29]
255 >> cuts[s][l][30] >> cuts[s][l][31] >> cuts[s][l][32] >> cuts[s][l][33] >> cuts[s][l][34] >> cuts[s][l][35] >> cuts[s][l][36] >> cuts[s][l][37] >> cuts[s][l][38] >> cuts[s][l][39]
256 >> cuts[s][l][40] >> cuts[s][l][41] >> cuts[s][l][42] >> cuts[s][l][43] >> cuts[s][l][44] >> cuts[s][l][45] >> cuts[s][l][46] >> cuts[s][l][47] >> cuts[s][l][48] >> cuts[s][l][49]
257 >> cuts[s][l][50] >> cuts[s][l][51] >> cuts[s][l][52] >> cuts[s][l][53] >> cuts[s][l][54] >> cuts[s][l][55] >> cuts[s][l][56] >> cuts[s][l][57] >> cuts[s][l][58] >> cuts[s][l][59]
258 >> cuts[s][l][60] >> cuts[s][l][61] >> cuts[s][l][62] >> cuts[s][l][63] >> cuts[s][l][64] >> cuts[s][l][65] >> cuts[s][l][66] >> cuts[s][l][67] >> cuts[s][l][68] >> cuts[s][l][69]
259 >> cuts[s][l][70] >> cuts[s][l][71] >> cuts[s][l][72] >> cuts[s][l][73] >> cuts[s][l][74] >> cuts[s][l][75] >> cuts[s][l][76] >> cuts[s][l][77] >> cuts[s][l][78] >> cuts[s][l][79]
260 >> cuts[s][l][80] >> cuts[s][l][81] >> cuts[s][l][82] >> cuts[s][l][83] >> cuts[s][l][84] >> cuts[s][l][85] >> cuts[s][l][86] >> cuts[s][l][87] >> cuts[s][l][88] >> cuts[s][l][89]
261 >> cuts[s][l][90] >> cuts[s][l][91] >> cuts[s][l][92] >> cuts[s][l][93] >> cuts[s][l][94] >> cuts[s][l][95] >> cuts[s][l][96] >> cuts[s][l][97] >> cuts[s][l][98] >> cuts[s][l][99]
262 ;
263 if(fin.peek()<0)break;
264 }
265 for(int is=0;is<=1;is++){
266 for(int il=0;il<43;il++){
268 int point(0);
269 //cout<<il<<" ";
270 for(int ipt=0;ipt<npt;ipt++){
271 //cout<<cuts[is][il][ipt]<<" ";
272 if(fabs(cuts[is][il][ipt])<1e-6)continue;
273 //cut[is][il]->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
274 gr->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
275 point++;
276 }
278 //cout<<endl;
279 }
280 }
281
282
284 if ( sc.isFailure() ){
285 log << MSG::FATAL << "Could not book Tuple !" << endreq;
286 return StatusCode::FAILURE;
287 }
288
289
290 myDotsHelixFitter.initialize();
291
292 return StatusCode::SUCCESS;
293}
Definition CgemCalibFunSvc.h:23
Definition CgemGeomSvc.h:25
static void setMdcCalibFunSvc(MdcCalibFunSvc *mdcCalibFunSvc)
Definition HoughHit.h:100
static void setCgemCalibFunSvc(CgemCalibFunSvc *cgemCalibFunSvc)
Definition HoughHit.h:102
static void setCgemGeomSvc(CgemGeomSvc *cgemGeomSvc)
Definition HoughHit.h:101
static void setMdcDetector(const MdcDetector *mdcDetector)
Definition HoughHit.h:103
Definition ICgemCalibFunSvc.h:10
Definition ICgemGeomSvc.h:22
Definition IMdcCalibFunSvc.h:17
Definition IRawDataProviderSvc.h:16
Definition MdcCalibFunSvc.h:22
Definition MdcDetector.h:20
Definition RawDataProviderSvc.h:17
Definition TrkContextEv.h:23
◆ makeHoughHitList()
| int HoughFinder::makeHoughHitList | ( | ) |
Definition at line 447 of file Hough.cxx.
448{
449 MsgStream log(msgSvc(), name());
450 int nHit(0);
451 if(!(m_mdcHitCol.empty())) m_mdcHitCol.clear();
452 if(!(m_houghHitList.empty())) m_houghHitList.clear();
453 if(!(m_XHoughHitList.empty())) m_XHoughHitList.clear();
454 if(!(m_VHoughHitList.empty())) m_VHoughHitList.clear();
455
456 if(m_cgem){
457 // --- RecCgemCluster
458 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
459 if(!recCgemClusterCol){
460 log << MSG::ERROR << "Could not retrieve RecCgemClusterCol" << endreq;
461 }else{
462 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterCol->begin();
463 m_recCgemClusterColBegin=cgemClusterIter;
464 for (;cgemClusterIter != recCgemClusterCol->end(); cgemClusterIter++,nHit++){
466 //m_recCgemCluster = (*cgemClusterIter);
467 HoughHit hit(recCgemCluster, m_bunchT0, nHit);
468 m_houghHitList.push_back(hit);
469 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
470 }
471 }
472 }
473
474 uint32_t getDigiFlag = 0;
479
480 MdcDigiVec::iterator mdcDigiIter = mdcDigiVec.begin();
481 for (;mdcDigiIter != mdcDigiVec.end(); mdcDigiIter++,nHit++){
483 HoughHit hit(mdcDigi, m_bunchT0, nHit);
484 if(fabs(hit.driftTime())>m_driftTimeUpLimit) continue;
486 hit.setMdcHit(mdcHit);
487
488 m_mdcHitCol.push_back(mdcHit);
489 m_houghHitList.push_back(hit);
490 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
491 }
492
494 int flag=hitIt->getFlag();
496 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
497 if(flag==2)m_VHoughHitList.push_back(&(*hitIt));
498 }
500 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
501 else m_VHoughHitList.push_back(&(*hitIt));
502 }
503 }
504 //for(vector<HoughHit*>::iterator hitIt = m_XHoughHitList.begin(); hitIt != m_XHoughHitList.end(); hitIt++){
505 //(*hitIt)->print();
506 //}
507 //for(vector<HoughHit*>::iterator hitIt = m_VHoughHitList.begin(); hitIt != m_VHoughHitList.end(); hitIt++){
508 //(*hitIt)->print();
509 //}
510
511 return nHit;
512}
Definition MdcHit.h:44
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition RawDataProviderSvc.h:29
Referenced by execute().
◆ printHitList()
| void HoughFinder::printHitList | ( | ) |
Definition at line 2691 of file Hough.cxx.
2692{
2693 //cout<<"nHit ==========> "<<m_houghHitList.size()<<endl;
2694 cout<<"========== truth =========="<<endl;
2696 hitIt->print();
2697 }
2698 cout<<endl;
2699 cout<<"========== Digi =========="<<endl;
2701 {
2702 hitIt->print();
2703 }
2704 cout<<endl;
2705 /*
2706 {
2707 vector<HoughHit>& hitList = m_mcHitCol;
2708 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2709 cout<<setw(4)<<hitIt->getHitID();
2710 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2711 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2712 vector<int> trkID = hitIt->getTrkID();
2713 cout<<"trkID:"<<setw(4)<<trkID[0];
2714 cout<<"half:"<<setw(4)<<hitIt->getHalfCircle();
2715 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
2716 if(hitIt->getCgemCluster()!=NULL){
2717 //cout<<hitIt->getFlag()<<" ";
2718 //if(hitIt->getFlag()==2)
2719 cout<<setw(5)<<hitIt->getCgemCluster()->getclusterid();
2720 cout<<"["<<setw(3)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(3)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2721 }
2722 //else cout<<"NULL ";
2723 }
2724 else{
2725 if(hitIt->getDigi()!=NULL){
2726 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2727 cout<<hitIt->getDigi()->identify();
2728 }
2729 //else cout<<"NULL ";
2730
2731 }
2732 cout<<endl;
2733 }
2734 }
2735 cout<<endl;
2736 //*/
2737 /*
2738 {
2739 vector<HoughHit>& hitList = m_houghHitList;
2740 cout<<hitList.size()<<endl;
2741 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2742 cout<<setw(4)<<hitIt->getHitID();
2743 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2744 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2745 cout<<setw(4)<<hitIt->getFlag();
2746 if(hitIt->getHitType()==HoughHit::CGEMHIT){
2747 cout<<setw(4)<<hitIt->getCgemCluster()->getclusterid();
2748 cout<<"["<<setw(4)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2749 }
2750 else{
2751 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2752 cout<<hitIt->getDigi()->identify()<<" ";
2753 //cout<<"("<<hitIt->getPairHit()->getLayer()<<", "<<hitIt->getPairHit()->getWire()<<") ";
2754 //if(hitIt->getPairHit()!=NULL)cout<<hitIt->getDriftDist()<<" "<<hitIt->getPairHit()->getDriftDist()<<" "<<hitIt->getDriftDist()-hitIt->getPairHit()->getDriftDist()<<endl;
2755 }
2756 if(hitIt->getPairHit()!=NULL)cout<<hitIt->getPairHit()->getHitID();
2757 //else cout<<"NULL";
2758 cout<<endl;
2759 }
2760 }
2761 cout<<endl;
2762 //*/
2763}
Referenced by execute().
◆ printTdsTrack()
| int HoughFinder::printTdsTrack | ( | ) |
Definition at line 4390 of file Hough.cxx.
4391{
4392 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
4393 if(!recMdcTrackCol)return 0;
4394 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
4395 double dr = (*iter)->helix(0) ;
4396 double phi0 = (*iter)->helix(1) ;
4397 double kappa = (*iter)->helix(2) ;
4398 double dz = (*iter)->helix(3) ;
4399 double tanl = (*iter)->helix(4) ;
4400 cout<<setw(12)<<"dr:" <<setw(15)<<dr
4401 <<setw(12)<<"phi0:" <<setw(15)<<phi0
4402 <<setw(12)<<"kappa:" <<setw(15)<<kappa
4403 <<setw(12)<<"dz:" <<setw(15)<<dz
4404 <<setw(12)<<"tanl:" <<setw(15)<<tanl
4405 <<endl;
4406
4407 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
4408 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
4409 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
4410 //cout<<(*clusterIt)->;
4411 }
4412
4413 HitRefVec hitRefVec = (*iter)->getVecHits();
4414 HitRefVec::iterator hotIt = hitRefVec.begin();
4415 //cout<<"("<<"layer"<<","<<"wire"<<") "
4416 cout<<"("<<"l "<<","<<" w "<<") "
4417 //<<"Stat "
4418 //<<"FlagLR "
4419 <<"S "
4420 <<"F "
4421 <<"DriftLeft "
4422 <<" DriftRight "
4423 <<"ErrDrift_L "
4424 <<"ErrDrift_R "
4425 <<"ChisqAdd "
4426 //<<" Tdc "
4427 //<<" Adc "
4428 <<" DriftT "
4429 <<" Doca "
4430 <<" Entra "
4431 <<" Zhit "
4432 <<" FltLen "
4433 <<endl;
4434 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
4437 cout<<"("<<setw( 2)<<layer<<","<<setw( 3)<<wire<<") "
4438 <<setw( 3)<<(*hotIt)->getStat()
4439 <<setw( 3)<<(*hotIt)->getFlagLR()
4440 <<setw(12)<<(*hotIt)->getDriftDistLeft()
4441 <<setw(12)<<(*hotIt)->getDriftDistRight()
4442 <<setw(12)<<(*hotIt)->getErrDriftDistLeft()
4443 <<setw(12)<<(*hotIt)->getErrDriftDistRight()
4444 <<setw(12)<<(*hotIt)->getChisqAdd()
4445 //<<setw( 6)<<(*hotIt)->getTdc()
4446 //<<setw( 6)<<(*hotIt)->getAdc()
4447 <<setw(10)<<(*hotIt)->getDriftT()
4448 <<setw(12)<<(*hotIt)->getDoca()
4449 <<setw(12)<<(*hotIt)->getEntra()
4450 <<setw(10)<<(*hotIt)->getZhit()
4451 <<setw(8 )<<(*hotIt)->getFltLen()
4452 <<endl;
4453 }
4454 cout<<endl;
4455 }
4456 return recMdcTrackCol->size();
4457}
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
Referenced by execute().
◆ registerTrack()
| StatusCode HoughFinder::registerTrack | ( | RecMdcTrackCol *& | trackList_tds, |
| RecMdcHitCol *& | hitList_tds ) |
Definition at line 1802 of file Hough.cxx.
1803{
1804 MsgStream log(msgSvc(), name());
1805 StatusCode sc;
1806 IDataProviderSvc* eventSvc = NULL;
1807 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1808 if (eventSvc) {
1809 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1810 } else {
1811 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1812 return StatusCode::FAILURE;
1813 //return StatusCode::SUCCESS;
1814 }
1815 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1816
1817
1818 // --- register RecMdcTrack
1820 DataObject *aRecMdcTrackCol;
1821 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1822 if(aRecMdcTrackCol != NULL) {
1823 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1824 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1825 }
1826
1827 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1828 if(sc.isFailure()) {
1829 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1830 return StatusCode::FAILURE;
1831 }
1832 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1833
1834
1835
1837 DataObject *aRecMdcHitCol;
1838 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1839 if(aRecMdcHitCol != NULL) {
1840 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1841 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1842 }
1843
1844 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1845 if(sc.isFailure()) {
1846 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1847 return StatusCode::FAILURE;
1848 }
1849 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1850
1851
1852 return StatusCode::SUCCESS;
1853}//end of stroeTracks
Referenced by execute().
◆ separateHoughHitByPhi()
| int HoughFinder::separateHoughHitByPhi | ( | ) |
◆ solveSharedHits()
| void HoughFinder::solveSharedHits | ( | vector< HoughHit * > & | hitList | ) |
Definition at line 3157 of file Hough.cxx.
3158{
3159 vector<HoughHit*>::iterator itHit = hitList.begin();
3160 for(; itHit!=hitList.end(); itHit++)
3161 {
3162 //vector<HoughTrack*> trkPntList = it->getTrkPntVec();
3163 //int nTrkSharing = trkPntList.size();
3164 vector<int> trkIdList = (*itHit)->getTrkID();
3165 int nTrkSharing = trkIdList.size();
3166 //cout<<"nTrkSharing = "<<nTrkSharing<<endl;
3167 if(nTrkSharing>1) {
3168 //cout<<"found one hit shared by "<<nTrkSharing<<" trks: layer "<<(*itHit)->getLayer()
3169 //<<", pos("<<(*itHit)->getHitPosition().x()<<","<<(*itHit)->getHitPosition().y()<<") with ";
3170 int trkId_minDelD=-1;
3171 double minResid = 99.0;
3172 //vector<int> trkIdList_toDel;
3173 vector<double> vecRes = (*itHit)->getVecResid();
3174 int nTrk = vecRes.size();
3175 //cout<<nTrk<<" residuals: ";
3176 for(int i=0; i<nTrk; i++)// find out the minimum residual
3177 {
3178 //cout<<" trk "<<trkIdList[i]<<":"<<vecRes[i];
3179 if(minResid>fabs(vecRes[i])) {
3180 minResid=fabs(vecRes[i]);
3181 trkId_minDelD=trkIdList[i];
3182 }
3183 }
3184 //cout<<endl;
3185 //cout<<" trkId_minDelD, minResid = "<<trkId_minDelD<<", "<<minResid<<endl;
3186 for(int i=0; i<nTrk; i++)// remove other associations
3187 {
3188 if(trkIdList[i]!=trkId_minDelD) {
3189 (*itHit)->dropTrkID(trkIdList[i]);// remove the trk from the hit
3190 vector<HoughTrack>::iterator itTrk = getHoughTrkIt(m_houghTrackList,trkIdList[i]);
3191 if(itTrk!=m_houghTrackList.end()) itTrk->dropHitPnt(&(*(*itHit)));// remove the hit from the trk
3192 }
3193 }
3194 trkIdList = (*itHit)->getTrkID();
3195 nTrkSharing = trkIdList.size();
3196 //cout<<"new nTrkSharing = "<<nTrkSharing<<endl;
3197 }// if(nTrkSharing>1)
3198 }// loop hits
3199}
vector< HoughTrack >::iterator getHoughTrkIt(vector< HoughTrack > &houghTrkList, int trkId)
Definition Hough.cxx:3202
Referenced by execute().
◆ storeRecTracks()
| int HoughFinder::storeRecTracks | ( | RecMdcTrackCol * | trackList, |
| RecMdcHitCol * | hitList, | ||
| vector< HoughTrack > & | trackVector ) |
Definition at line 1896 of file Hough.cxx.
1897{
1898 bool debug=false; //debug=true;
1899 int trackId(0);
1900 int tkStat =4;
1901 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
1902 sort(trackVector.begin(),trackVector.end(),largerPt);
1903
1904 for(vector<HoughTrack>::iterator trkIT = trackVector.begin(); trkIT!=trackVector.end(); trkIT++)
1905 {
1906
1907 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
1908
1909 //if((trkIT)->getFlag() == 1)continue;
1910 //if((trkIT)->getFlag() == -1)continue;
1911 //if((trkIT)->getFlag() == -2)continue;
1912 if((trkIT)->getFlag()!=0) continue;
1913
1914 if(debug)
1915 {
1916 cout<<"the following track will be saved: "<<endl;
1917 trkIT->print();
1918 }
1919
1920 //new a Rec. Track MdcTrack
1922
1923 recMdcTrack->setTrackId(trackId);
1924
1925 double helixPar[5];
1926 helixPar[0]=trkIT->getDr();
1927 helixPar[1]=trkIT->getPhi0();
1928 helixPar[2]=trkIT->getKappa();
1929 helixPar[3]=trkIT->getDz();
1930 helixPar[4]=trkIT->getTanl();//helixPar[4]+=0.1;
1931 recMdcTrack->setHelix(helixPar);
1932
1934 double pxy = trkIT->pt();
1935 double px = trkIT->momentum(0).x();
1936 double py = trkIT->momentum(0).y();
1937 double pz = trkIT->momentum(0).z();
1938 double p = trkIT->momentum(0).mag();
1939 double theta = trkIT->direction(0).theta();
1940 double phi = trkIT->direction(0).phi();
1941 HepPoint3D poca = trkIT->x(0);
1942 HepPoint3D pivot = trkIT->pivot();
1943 double r = poca.perp();
1944 HepSymMatrix Ea = trkIT->Ea();
1945 //cout<<"Ea="<<Ea<<endl;
1946 double errorMat[15];
1947 int k = 0;
1948 for (int ie = 0 ; ie < 5 ; ie ++){
1949 for (int je = ie ; je < 5 ; je ++){
1950 errorMat[k] = Ea[ie][je];
1951 k++;
1952 }
1953 }
1954 // ---> a Helix for DotsHelixFitter
1955 HepVector a(5,0);
1956 a(1) =trkIT->getDr();
1957 a(2) =trkIT->getPhi0();
1958 a(3) =trkIT->getKappa();
1959 a(4) =trkIT->getDz();
1960 a(5) =trkIT->getTanl();
1962 //myDotsHelixFitter.setInitialHelix(aHelix);// set a Helix for DotsHelixFitter
1963 double chisq = trkIT->getChi2();
1965 recMdcTrack->setPxy(pxy);
1966 recMdcTrack->setPx(px);
1967 recMdcTrack->setPy(py);
1968 recMdcTrack->setPz(pz);
1969 recMdcTrack->setP(p);
1970 recMdcTrack->setTheta(theta);
1971 recMdcTrack->setPhi(phi);
1972 recMdcTrack->setPoca(poca);
1974 recMdcTrack->setY(poca.y());
1975 recMdcTrack->setZ(poca.z());
1976 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
1977 HepPoint3D apivot(0.,0.,0.);
1978 recMdcTrack->setPivot(apivot);
1980 recMdcTrack->setVY0(0.);
1981 recMdcTrack->setVZ0(0.);
1982 recMdcTrack->setError(errorMat);
1983 recMdcTrack->setError(Ea);
1984 recMdcTrack->setChi2(chisq);
1985 recMdcTrack->setStat(tkStat);
1986
1987 //-----------hit list-------------
1988 int hitId = 0;
1989 HitRefVec hitRefVec;
1990 ClusterRefVec clusterRefVec;
1991 //vector<int> vecHits;
1992 map<int,int> clusterFitStat;
1993 //int maxLayer = 0;
1994 int maxLayerId = -1;
1995 int minLayerId = 43;
1996 double fiTerm = 999.;
1997 HoughHit* fiTermHot = NULL;
1998
1999 //vector<HoughHit*> hotList = trkIT->getVecHitPnt();
2000 vector<HoughHit*> hotList = trkIT->getHotList(2);
2001 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2002 {
2003 //cout<<"handle hit "<<hitId<<endl;
2004 HoughHit* hot(*hotIt);
2005 //int layer = (**hotIt).getLayer();
2006 //int wire = (**hotIt).getWire();
2007 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2008
2009 /*
2010 if(hot->getHitType()==HoughHit::MDCHIT)
2011 {
2012 //if(hot->getUse()!=1) continue;
2013 //cout<<"create a new RecMdcHit"<<endl;
2014
2015 // --- make a RecMdcHit from digi, aHelix using DotsHelixFitter
2016 RecMdcHit* recMdcHit = new RecMdcHit(myDotsHelixFitter.makeRecMdcHit(hot->getDigi(), aHelix));
2017 recMdcHit->setId(hitId);
2018 recMdcHit->setTrkId(trackId);
2019
2020 //
2021 // --> make a RecMdcHit from a HoughHit (parameters from Hough)
2022 RecMdcHit* recMdcHit = new RecMdcHit;
2023 recMdcHit->setId(hitId);
2024 recMdcHit->setTrkId(trackId);
2025 HepPoint3D hitPoint = hot->getHitPosition();
2026 double dPhi = trkIT->dPhi(hitPoint);
2027 HepPoint3D position = trkIT->x(dPhi);
2028 Hep3Vector direction = trkIT->direction(dPhi);
2029 //cout<<"get 1"<<endl;
2030
2031 double ddl = -1*hot->getDriftDist();
2032 double ddr = hot->getDriftDist();
2033 //double erddl = hot->hitSigma();
2034 //double erddr = hot->hitSigma();
2035 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2036 int lr = 0;
2037 int stat = hot->getUse();
2038 stat = 1; // llwang
2039 //cout<<"hit stat = "<<stat<<endl;
2040 //cout<<"get 2"<<endl;
2041 Identifier mdcId = hot->getDigi()->identify();
2042 //cout<<"get 2.1"<<endl;
2043 double tdc = hot->getDigi()->getTimeChannel();
2044 //cout<<"get 2.2"<<endl;
2045 double adc = hot->getDigi()->getChargeChannel();
2046 //cout<<"get 2.3"<<endl;
2047 double driftT = hot->driftTime();
2048 //cout<<"get 2.4"<<endl;
2049 double doca = hot->getDriftDist() +hot->getResidual();
2050 //cout<<"get 2.5"<<endl;
2051 double entra = direction.phi()-position.phi();//FIXME
2052 //cout<<"get 2.6"<<endl;
2053 while(entra<-0.5*M_PI)entra+= M_PI;
2054 while(entra> 0.5*M_PI)entra -= M_PI;
2055 //cout<<"get 3"<<endl;
2056 double zhit = hitPoint.z();
2057 double fltLen = trkIT->flightLength(hitPoint);
2058 //cout<<"after get "<<hitId<<endl;
2059
2060 recMdcHit->setDriftDistLeft(ddl);
2061 recMdcHit->setDriftDistRight(ddr);
2062 //recMdcHit->setErrDriftDistLeft(erddl);
2063 //recMdcHit->setErrDriftDistRight(erddr);
2064 recMdcHit->setChisqAdd(pChisq);
2065 recMdcHit->setFlagLR(lr);
2066 recMdcHit->setStat(stat);
2067 recMdcHit->setMdcId(mdcId);
2068 recMdcHit->setTdc(tdc);
2069 recMdcHit->setAdc(adc);
2070 recMdcHit->setDriftT(driftT);
2071 recMdcHit->setDoca(doca);
2072 recMdcHit->setEntra(entra);
2073 recMdcHit->setZhit(zhit);
2074 recMdcHit->setFltLen(fltLen);
2075 //cout<<"after set "<<hitId<<endl;
2076 // <--- make a RecMdcHit from a HoughHit
2077
2078 hitList->push_back(recMdcHit);
2079 SmartRef<RecMdcHit> refHit(recMdcHit);
2080 hitRefVec.push_back(refHit);
2081
2082 if(hot->getLayer()>maxLayerId){
2083 maxLayerId = hot->getLayer();
2084 fiTermHot = hot;
2085 }
2086 if(hot->getLayer()<minLayerId){
2087 minLayerId = hot->getLayer();
2088 }
2089 //if(hot->getUse()==1){
2090 recMdcHit->setStat(1);
2091 //}else{
2092 // recMdcHit->setStat(0);
2093 //}
2094 //cout<<"finish a DC hit "<<hitId<<endl;
2095 hitId++;
2096 }*/
2098 {
2099 //hot->print();
2100 //cout<<endl;
2103 //int stat = hot->getUse();
2104 int stat = 1;
2105 clusterFitStat[clusterid] = stat;
2106 clusterRefVec.push_back(recCgemCluster);
2107 if(debug) cout<<" cluster "<<clusterid<<" kept! "<<endl;
2108 if(hot->getLayer()>maxLayerId)
2109 {
2110 maxLayerId = hot->getLayer();
2111 fiTermHot = hot;
2112 }
2113 if(hot->getLayer()<minLayerId)
2114 {
2115 minLayerId = hot->getLayer();
2116 }
2117 }
2118 }// loop hotList (vector<HoughHit*>) of a hough track
2119
2120 // --- loop RecMdcHitVec in Hough track
2121 double fltLen = -1;
2122 vector<RecMdcHit>* RecMdcHitVecPnt = trkIT->getRecMdcHitVec();
2123 int nMdcHits=RecMdcHitVecPnt->size();
2124 int nMdcHitsKept=0;
2125 vector<RecMdcHit>::iterator iter_recMdcHit = RecMdcHitVecPnt->begin();
2126 for(; iter_recMdcHit!=RecMdcHitVecPnt->end(); iter_recMdcHit++)
2127 {
2130 {
2131 //cout<<"skip hit with chi2="<<recMdcHit->getChisqAdd()<<endl;
2132 delete recMdcHit;
2133 continue;
2134 }
2135 recMdcHit->setId(hitId);
2136 recMdcHit->setTrkId(trackId);
2137 recMdcHit->setStat(1);
2138 hitList->push_back(recMdcHit);
2139 SmartRef<RecMdcHit> refHit(recMdcHit);
2140 hitRefVec.push_back(refHit);
2141 nMdcHitsKept++;
2142
2146 if(debug) cout<<" hit at layer "<<layer<<" wire "<<wire<<" kept! "<<endl;
2147 if(layer>maxLayerId)
2148 {
2149 maxLayerId = layer;
2150 fltLen=recMdcHit->getFltLen();
2151 //fiTermHot = hot;
2152 }
2153 if(layer<minLayerId)
2154 {
2155 minLayerId = layer;
2156 }
2157 hitId++;
2158 }
2159 if(debug) cout<<"track "<<trackId<<", "<<nMdcHitsKept<<"/"<<nMdcHits<<" hits kept"<<endl;
2160
2161 //fi terminal angle
2162 if (fiTermHot!=NULL){
2163 //fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
2165 fiTerm = trkIT->flightArc(point)/trkIT->radius();
2166 }
2168 recMdcTrack->setFiTerm(fiTerm);
2169
2170 // --- setN* functions called in setVec* functions
2171 recMdcTrack->setVecHits(hitRefVec);
2172 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2173 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2174 trackList->push_back(recMdcTrack);
2175 trackId++;
2176 }// loop HoughTrack list
2177 return trackList->size();
2178
2179}
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition CgemMdcFitAlg.cxx:768
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
Helix parameter class.
Definition RecMdcHit.h:22
Definition RecMdcTrack.h:30
void setVecClusters(ClusterRefVec vecclusters)
Definition RecMdcTrack.cxx:86
Referenced by execute().
◆ storeTrack()
| int HoughFinder::storeTrack | ( | vector< HoughTrack > & | trackVector, |
| RecMdcTrackCol *& | recMdcTrackCol, | ||
| RecMdcHitCol *& | recMdcHitCol ) |
Definition at line 1855 of file Hough.cxx.
1856{
1857 MdcTrackParams mdcTrackParams;
1858 //mdcTrackParams.lPrint=8;
1859 mdcTrackParams.lRemoveInActive=1;
1860 mdcTrackParams.maxGapLength=99;
1861 mdcTrackParams.nHitDeleted=99;
1862 MdcTrackList mdcTrackList(mdcTrackParams);
1863 mdcTrackList.setNoInner(true);
1864 vector<MdcTrack*> mdcTrackVector;
1865 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
1866 if((trkIT)->getFlag() == 1)continue;
1867 if((trkIT)->getFlag() == -1)continue;
1868 if((trkIT)->getFlag() == -2)continue;
1869 //(trkIT)->sortHot();
1871 mdcTrackList.append(mdcTrack);
1872 mdcTrackVector.push_back(mdcTrack);
1873 }
1874 int nDeleted(0);
1875 mdcTrackList.setNoInner(true);
1876 if(m_checkHits)nDeleted = mdcTrackList.arbitrateHits();
1877 if(nDeleted>0)cout<<"nDeleted "<<nDeleted<<endl;
1878
1879 int trkID=0;
1880 for(int l=0;l<mdcTrackList.length();l++){
1881 mdcTrackList[l]->storeTrack(trkID,recMdcTrackCol,recMdcHitCol,4);
1882 trkID++;
1883 }
1884
1885 vector<MdcTrack*>::iterator iter = mdcTrackVector.begin();
1888 }
1889 return trkID;
1890}
Definition MdcTrackList.h:31
Definition MdcTrackParams.h:13
Definition MdcTrack.h:19
Referenced by execute().
◆ storeTracks()
| int HoughFinder::storeTracks | ( | RecMdcTrackCol * | trackList, |
| RecMdcHitCol * | hitList, | ||
| vector< HoughTrack > & | trackVector ) |
Definition at line 2181 of file Hough.cxx.
2181 {
2182 bool debug=false; //debug=true;
2183 int trackId(0);
2184 int tkStat =4;
2185 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
2186 sort(trackVector.begin(),trackVector.end(),largerPt);
2187
2188 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2189 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
2190 if((trkIT)->getFlag() == 1)continue;
2191 if((trkIT)->getFlag() == -1)continue;
2192 if((trkIT)->getFlag() == -2)continue;
2193 if(trkIT->getTrkRecoTrk()==NULL)
2194 {
2195 if(debug) cout<<"trk getTrkRecoTrk NULL!"<<endl;
2196 continue;
2197 }
2198 //get the results of the fit to this track
2200 //check the fit worked
2201 if (fitresult == 0)
2202 {
2203 if(debug) cout<<"no fit result!"<<endl;
2204 continue;
2205 }
2206 if(debug)
2207 {
2208 cout<<"the following track will be saved: "<<endl;
2209 trkIT->print();
2210 }
2213 //cout<<endl;
2214
2215 //new a Rec. Track MdcTrack
2218 //some track info
2221 //Fit info
2222 int hitId = 0;
2223 int nHits=0;
2224 int nSt=0;
2225 //int nAct=0; int nSeg=0;
2226 //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
2227 //****************************
2228 //* active flag open this
2229 //****************************
2230 //* if (allHit>0){
2231 //* nHits = aList->nHit();//add inActive hit also
2232 //* }else{
2233 //* nHits = fitresult->nMdc();// = nActive()
2234 //* }
2235 //****************************
2236 //for 5.1.0 use all hits
2237 nHits = aList->nHit();
2238 //****************************
2239 // use active only
2240 // nHits = fitresult->nMdc();// = nActive()
2241 //****************************
2242
2246
2247 //if(nHits-nDof!=5) cout<<"HoughFinder::storeTracks: nHits-nDof!=5, nHits="<<nHits<<" , nDof="<<nDof<<", chi2="<<chisq<<" , evt="<<m_event<<endl;
2248
2249 //track parameters at closest approach to beamline
2250 double fltLenPoca = 0.0;
2252 //std::cout<< __FILE__ << " phi0 " << helix.phi0() << " omega " <<helix.omega()<<std::endl;
2255 //double z0 = helix.z0();
2257 //momenta and position
2258 //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
2260
2261 double helixPar[5];
2262 //dro =-d0
2263 helixPar[0] = -d0;
2264 //phi0 = phi0 - pi/2 [0,2pi)
2267 helixPar[1] = tphi0;
2268 //kappa = q/pxy
2270 if (pxy == 0.) helixPar[2] = 9999.;
2272 if(pxy>9999.) helixPar[2] = 0.00001;
2273 //dz = z0
2274 helixPar[3] = helix.z0();
2275 //tanl
2276 helixPar[4] = tanDip;
2277 //error
2278 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
2279 HepSymMatrix mS(helix.params().num_row(),0);
2280 mS[0][0]=-1.;//dr0=-d0
2281 mS[1][1]=1.;
2282 mS[2][2]=-333.567/Bz;//pxy -> cpa
2283 mS[3][3]=1.;
2284 mS[4][4]=1.;
2285 HepSymMatrix mVy = helix.covariance().similarity(mS);
2286 double errorMat[15];
2287 int k = 0;
2288 for (int ie = 0 ; ie < 5 ; ie ++){
2289 for (int je = ie ; je < 5 ; je ++) {
2290 errorMat[k] = mVy[ie][je];
2291 k++;
2292 }
2293 }
2294 double p,px,py,pz;
2295 px = pxy * (-sin(helixPar[1]));
2296 py = pxy * cos(helixPar[1]);
2297 pz = pxy * helixPar[4];
2298 p = sqrt(pxy*pxy + pz*pz);
2299 //theta, phi
2300 double theta = acos(pz/p);
2301 double phi = atan2(py,px);
2302 recMdcTrack->setTrackId(trackId);
2303 recMdcTrack->setHelix(helixPar);
2305 recMdcTrack->setPxy(pxy);
2306 recMdcTrack->setPx(px);
2307 recMdcTrack->setPy(py);
2308 recMdcTrack->setPz(pz);
2309 recMdcTrack->setP(p);
2310 recMdcTrack->setTheta(theta);
2311 recMdcTrack->setPhi(phi);
2312 recMdcTrack->setPoca(poca);
2314 recMdcTrack->setY(poca.y());
2315 recMdcTrack->setZ(poca.z());
2316 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
2317 HepPoint3D pivot(0.,0.,0.);
2318 recMdcTrack->setPivot(pivot);
2320 recMdcTrack->setVY0(0.);
2321 recMdcTrack->setVZ0(0.);
2322 recMdcTrack->setError(errorMat);
2323 recMdcTrack->setError(mVy);
2324 recMdcTrack->setChi2(chisq);
2325 recMdcTrack->setStat(tkStat);
2326 //recMdcTrack->setNhits(nHits);
2327 //-----------hit list-------------
2328 HitRefVec hitRefVec;
2329 ClusterRefVec clusterRefVec;
2330 vector<int> vecHits;
2331 map<int,int> clusterFitStat;
2332 //for fiterm
2333 int maxLayer = 0;
2334 //cout<<__FILE__ <<" "<<__LINE__ <<endl;
2336 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2341 if(maxLayer < layerId)maxLayer = layerId;
2342 }
2343 }
2344 //cout<<"maxLayer:"<<maxLayer<<endl;
2345 int maxLayerId = -1;
2346 int minLayerId = 43;
2347 double fiTerm = 999.;
2351 //cout<<__FILE__ <<" "<<__LINE__ <<" "<< typeid(*hot).name()<<endl;
2352 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2357 if(layerId>(maxLayer - m_removeNOuterHits))continue;
2358 //if (!recoHot->mdcHit()) return;
2360 //id
2361 recMdcHit->setId(hitId);
2362 //---------BESIII----------
2363 //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
2364 //phiWire - phiHit >0 flagLR=1 right driftright>0;
2365 //flag = 2 ambig;
2366 //-----Babar wireAmbig----
2367 //-1 -> right, 0 -> don't know, +1 -> left
2368 //+1 phiWire-phiHit<0 Left,
2369 //-1 phiWire-phiHit>0 Right,
2370 //0 don't know
2371 //ambig() w.r.t track direction
2372 //wireAmbig() w.r.t. wire location
2377 int flagLR=2;
2378 if ( hotWireAmbig == 1){
2379 flagLR = 0; //left driftDist <0
2380 doca *= -1.;//2012-07-19
2381 }else if( hotWireAmbig == -1){
2382 flagLR = 1; //right driftDist >0
2383 }else if( hotWireAmbig == 0){
2384 flagLR = 2; //don't know
2385 }
2386 recMdcHit->setFlagLR(flagLR);
2387 recMdcHit->setDriftDistLeft((-1)*driftDist);
2388 recMdcHit->setErrDriftDistLeft(sigma);
2389 recMdcHit->setDriftDistRight(driftDist);
2390 recMdcHit->setErrDriftDistRight(sigma);
2391 //Mdc Id
2393 recMdcHit->setMdcId(mdcId);
2394 //corrected ADC
2395
2396 //contribution to chisquare
2397 //fill chisq
2398 double res=999.,rese=999.;
2400 }else{}
2401 double deltaChi=0;
2403 recMdcHit->setChisqAdd( deltaChi * deltaChi );
2404 //set tracks containing this hit
2405 recMdcHit->setTrkId(trackId);
2406 //doca
2408 //entrance angle
2409
2411 //z of hit
2412 HepPoint3D pos; Hep3Vector dir;
2414 recMdcHit->setFltLen(fltLen);
2416 recMdcHit->setZhit(pos.z());
2420 //drift time
2422 //for fiterm
2424 //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig "<<hotWireAmbig<<" y "<<pos.y()<<std::endl;
2425 // <<recoHot->entranceAngle()<<std::endl;
2426 if (layerId >= maxLayerId){
2427 maxLayerId = layerId;
2428 fiTermHot = recoHot;
2429 }
2430 if (layerId < minLayerId){
2431 minLayerId = layerId;
2432 }
2433 // status flag
2435 recMdcHit->setStat(1);
2436 }else{
2437 recMdcHit->setStat(0);
2438 }
2439 // for 5.1.0 use all hits
2441 ++nSt;
2442 }
2443 hitList->push_back(recMdcHit);
2444 SmartRef<RecMdcHit> refHit(recMdcHit);
2445 hitRefVec.push_back(refHit);
2446 vecHits.push_back(mdcId.get_value());
2447 ++hitId;
2453 if(stat==0) continue;// added 2020-05-26 wangll
2454 //cout<<"stat:"<<stat<<endl;
2455 clusterFitStat[clusterid] = stat;
2456 //cout<<clusterid<<" "<<stat<<endl;
2457 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2458 clusterRefVec.push_back(recCgemCluster);
2459 }
2460 }
2461 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2462 //fi terminal angle
2463 if (fiTermHot!=NULL){
2465 }
2466 recMdcTrack->setFiTerm(fiTerm);
2467 // number of stereo hits contained
2468 //recMdcTrack->setNhits(hitId);
2469 //recMdcTrack->setNster(nSt);
2470 //recMdcTrack->setNdof(nDof);
2471 //recMdcTrack->setFirstLayer(maxLayerId);
2472 //recMdcTrack->setLastLayer(minLayerId);
2473 //recMdcTrack->setFirstLayer(minLayerId);
2474 //recMdcTrack->setLastLayer(maxLayerId);
2475 //recMdcTrack->setClusterFitStat(clusterFitStat);
2476 //recMdcTrack->setVecClusters(clusterRefVec);
2477 //recMdcTrack->setNcluster(clusterRefVec.size());
2478
2479 // --- setN* functions called in setVec* functions
2480 recMdcTrack->setVecHits(hitRefVec);
2481 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2482 trackList->push_back(recMdcTrack);
2483 trackId++;
2484 }
2485 //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
2486 /*
2487 int trackId(0);
2488 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2489//TrkErrCode trkErrCode = trkIT->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
2490HoughTrack* houghTrack = &(*trkIT);
2491if(houghTrack->getFlag() == 1)continue;
2492if(houghTrack->getFlag() < 0)continue;
2493//houghTrack->print();
2494//houghTrack->printHot();
2495trackId=trackList->size();
2496//cout<<"start store trk "<<trackId;
2497RecMdcTrack* recMdcTrack = new RecMdcTrack();
2498//vector<HoughHit>* hotList = houghTrack->getHot();
2499//vector<HoughHit*> hotList = houghTrack->getVecHitPnt();
2500vector<HoughHit*> hotList = houghTrack->getHotList();
2501//cout<<" with "<<hotList.size()<<" hits";//<<endl;
2502double helixPar[5];
2503helixPar[0] = houghTrack->dr();
2504helixPar[1] = houghTrack->phi0();
2505helixPar[2] = houghTrack->kappa();
2506helixPar[3] = houghTrack->dz();
2507helixPar[4] = houghTrack->tanl();
2508int q = houghTrack->getCharge();
2509double pxy = houghTrack->pt();
2510double px = houghTrack->momentum(0).x();
2511double py = houghTrack->momentum(0).y();
2512double pz = houghTrack->momentum(0).z();
2513double p = houghTrack->momentum(0).mag();
2514double theta = houghTrack->direction(0).theta();
2515double phi = houghTrack->direction(0).phi();
2516HepPoint3D poca = houghTrack->x(0);
2517HepPoint3D pivot = houghTrack->pivot();
2518double r = poca.perp();
2519HepSymMatrix Ea = houghTrack->Ea();
2520//cout<<" phi, pxy, pz = "<<phi<<", "<<pxy<<", "<<pz<<endl;
2521double errorMat[15];
2522int k = 0;
2523for (int ie = 0 ; ie < 5 ; ie ++){
2524for (int je = ie ; je < 5 ; je ++){
2525errorMat[k] = Ea[ie][je];
2526k++;
2527}
2528}
2529double chisq = houghTrack->getChi2();
2530int tkStat = 4+houghTrack->getCircleFitStat();
2531//if(tkStat==4) cout<<" track "<<trackId<<" has no good circle fit"<<endl;
2532int nHits = hotList.size();
2533
2534recMdcTrack->setTrackId(trackId);
2535recMdcTrack->setHelix(helixPar);
2536recMdcTrack->setCharge(q);
2537recMdcTrack->setPxy(pxy);
2538recMdcTrack->setPx(px);
2539recMdcTrack->setPy(py);
2540recMdcTrack->setPz(pz);
2541recMdcTrack->setP(p);
2542recMdcTrack->setTheta(theta);
2543recMdcTrack->setPhi(phi);
2544recMdcTrack->setPoca(poca);
2545recMdcTrack->setX(poca.x());//poca
2546recMdcTrack->setY(poca.y());
2547recMdcTrack->setZ(poca.z());
2548recMdcTrack->setR(r);
2549recMdcTrack->setPivot(pivot);
2550recMdcTrack->setVX0(0.);//pivot
2551recMdcTrack->setVY0(0.);
2552recMdcTrack->setVZ0(0.);
2553recMdcTrack->setError(errorMat);
2554recMdcTrack->setError(Ea);
2555recMdcTrack->setChi2(chisq);
2556recMdcTrack->setStat(tkStat);
2557recMdcTrack->setNhits(nHits);
2558//-----------hit list-------------
2559int hitId = 0;
2560int nSt(0);
2561int maxLayerId = -1;
2562int minLayerId = 43;
2563double fiTerm = 999.;
2564HoughHit fiTermHot;
2565
2566HitRefVec hitRefVec;
2567ClusterRefVec clusterRefVec;
2568map<int,int> clusterFitStat;
2569for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2570{
2571 //cout<<"handle hit "<<hitId<<endl;
2572 HoughHit* hot(*hotIt);
2573 //int layer = (**hotIt).getLayer();
2574 //int wire = (**hotIt).getWire();
2575 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2576 if(hot->getHitType()==HoughHit::MDCHIT){
2577 //cout<<"create a new RecMdcHit"<<endl;
2578 RecMdcHit* recMdcHit = new RecMdcHit;
2579
2580 recMdcHit->setId(hitId);
2581 recMdcHit->setTrkId(trackId);
2582 HepPoint3D hitPoint = hot->getHitPosition();
2583 double dPhi = houghTrack->dPhi(hitPoint);
2584 HepPoint3D position = houghTrack->x(dPhi);
2585 Hep3Vector direction = houghTrack->direction(dPhi);
2586 //cout<<"get 1"<<endl;
2587
2588 double ddl = -1*hot->getDriftDist();
2589 double ddr = hot->getDriftDist();
2590 //double erddl = hot->hitSigma();
2591 //double erddr = hot->hitSigma();
2592 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2593 int lr = 0;
2594 int stat = hot->getUse();
2595 stat = 1; // llwang
2596 //cout<<"hit stat = "<<stat<<endl;
2597 //cout<<"get 2"<<endl;
2598 Identifier mdcId = hot->getDigi()->identify();
2599 //cout<<"get 2.1"<<endl;
2600 double tdc = hot->getDigi()->getTimeChannel();
2601 //cout<<"get 2.2"<<endl;
2602 double adc = hot->getDigi()->getChargeChannel();
2603 //cout<<"get 2.3"<<endl;
2604 double driftT = hot->driftTime();
2605 //cout<<"get 2.4"<<endl;
2606 double doca = hot->getDriftDist() +hot->getResidual();
2607 //cout<<"get 2.5"<<endl;
2608 double entra = direction.phi()-position.phi();//FIXME
2609 //cout<<"get 2.6"<<endl;
2610 while(entra<-0.5*M_PI)entra+= M_PI;
2611 while(entra> 0.5*M_PI)entra -= M_PI;
2612 //cout<<"get 3"<<endl;
2613 double zhit = hitPoint.z();
2614 double fltLen = houghTrack->flightLength(hitPoint);
2615 //cout<<"after get "<<hitId<<endl;
2616
2617 recMdcHit->setDriftDistLeft(ddl);
2618 recMdcHit->setDriftDistRight(ddr);
2619 //recMdcHit->setErrDriftDistLeft(erddl);
2620 //recMdcHit->setErrDriftDistRight(erddr);
2621 recMdcHit->setChisqAdd(pChisq);
2622 recMdcHit->setFlagLR(lr);
2623 recMdcHit->setStat(stat);
2624 recMdcHit->setMdcId(mdcId);
2625 recMdcHit->setTdc(tdc);
2626 recMdcHit->setAdc(adc);
2627 recMdcHit->setDriftT(driftT);
2628 recMdcHit->setDoca(doca);
2629 recMdcHit->setEntra(entra);
2630 recMdcHit->setZhit(zhit);
2631 recMdcHit->setFltLen(fltLen);
2632 //cout<<"after set "<<hitId<<endl;
2633
2634 hitList->push_back(recMdcHit);
2635 SmartRef<RecMdcHit> refHit(recMdcHit);
2636 hitRefVec.push_back(refHit);
2637
2638 if(hot->getLayer()>maxLayerId){
2639 maxLayerId = hot->getLayer();
2640 fiTermHot = *hot;
2641 }
2642 if(hot->getLayer()<minLayerId){
2643 minLayerId = hot->getLayer();
2644 }
2645 if(hot->getUse()==1){
2646 recMdcHit->setStat(1);
2647 }else{
2648 recMdcHit->setStat(0);
2649 }
2650 if(hot->getFlag()!=0)++nSt;
2651 //cout<<"finish a DC hit "<<hitId<<endl;
2652 }else if(hot->getHitType()==HoughHit::CGEMHIT){
2653 //hot->print();
2654 //cout<<endl;
2655 const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
2656 int clusterid = hot->getCgemCluster()->getclusterid();
2657 //int stat = hot->getUse();
2658 int stat = 1;
2659 clusterFitStat[clusterid] = stat;
2660 clusterRefVec.push_back(recCgemCluster);
2661 }
2662 hitId++;
2663}
2664
2665//HepPoint3D point = fiTermHot.getHitPosition();
2666//fiTerm = houghTrack->flightArc(point)/houghTrack->radius();
2667//recMdcTrack->setFiTerm(fiTerm);
2668recMdcTrack->setNdof(nHits-5);
2669recMdcTrack->setNster(nSt);
2670recMdcTrack->setFirstLayer(minLayerId);
2671recMdcTrack->setLastLayer(maxLayerId);
2672//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2673recMdcTrack->setVecHits(hitRefVec);
2674//cout<<"before setVecClusters"<<endl;
2675//recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2676//recMdcTrack->setVecClusters(clusterRefVec);
2677recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2678//cout<<"clusterRefVec.size = "<<clusterRefVec.size()<<endl;
2679//cout<<"end setVecClusters"<<endl;
2680recMdcTrack->setNcluster(clusterRefVec.size());
2681//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2682trackList->push_back(recMdcTrack);
2683//trackId++;
2684//cout<<helixPar[2]<<endl;
2685}
2686return trackList->size();
2687//*/
2688return trackList->size();
2689}
double bFieldZ() const
Definition CgemHitOnTrack.h:44
Definition MdcHitOnTrack.h:40
Definition MdcRecoHitOnTrack.h:27
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &direction) const =0
virtual void print(std::ostream &ostr) const =0
virtual double pt(double fltL=0.) const =0
virtual double chisq() const =0
virtual int charge() const =0
virtual int nDof() const =0
virtual const TrkDifTraj & traj() const =0
virtual HepPoint3D position(double fltL) const =0
Definition TrkExchangePar.h:30
Definition TrkFit.h:23
virtual TrkExchangePar helix(double fltL) const =0
Definition TrkHitList.h:35
Definition TrkHitOnTrkIter.h:34
TrkErrCode getFitStuff(HepVector &derivs, double &deltaChi) const
Definition TrkHitOnTrk.cxx:248
Referenced by execute().
Member Data Documentation
◆ Layer
| int HoughFinder::Layer |
Definition at line 63 of file Hough.h.
Referenced by HoughFinder().
◆ m_clearTrack
| int HoughFinder::m_clearTrack |
Definition at line 65 of file Hough.h.
Referenced by HoughFinder(), initialize(), and storeTrack().
◆ m_useCgemInGlobalFit
| int HoughFinder::m_useCgemInGlobalFit |
Definition at line 64 of file Hough.h.
Referenced by HoughFinder(), and initialize().
◆ printFlag
| int HoughFinder::printFlag |
Definition at line 68 of file Hough.h.
Referenced by execute(), getMcParticleCol(), HoughFinder(), and storeTracks().
◆ storeFlag
| int HoughFinder::storeFlag |
Definition at line 66 of file Hough.h.
Referenced by execute(), and HoughFinder().
The documentation for this class was generated from the following files:
