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 for circle
92 declareProperty("circle_chiCut", m_circle_chiCut=25.0);
93
94 // --- cuts to combine tracks
95 declareProperty("drCut_combine", m_drCut=2.0);
96 declareProperty("phi0Cut_combine", m_phi0Cut=0.2);
97 declareProperty("kappaCut_combine",m_kappaCut=1.5);
98 declareProperty("dzCut_combine", m_dzCut=5.0);
99 declareProperty("tanlCut_combine", m_tanlCut=0.5);
100
101 // --- cuts for hits saved in RecMdcTrack
102 declareProperty("chi2CutHits", m_chi2CutHits=25);
103
104 m_maxFireLayer = -1;
105 m_totEvtProcessed=0;
106}
Member Function Documentation
◆ activeUnusedCgemHitsOnly()
| int HoughFinder::activeUnusedCgemHitsOnly | ( | vector< HoughHit * > & | hitPntList | ) |
Definition at line 3228 of file Hough.cxx.
3229{
3230 int nCgemLeft(0);
3231 vector<HoughHit*>::iterator iter = hitPntList.begin();
3233 {
3234 //if((*iter)->getPairHit()==NULL) continue;// FIXME
3235 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
3236 int useOld = (*iter)->getUse();
3237 int iLayer = (*iter)->getLayer();
3238 if(iLayer<3) {
3239 if(useOld==-1||useOld==0) {
3240 (*iter)->setUse(0);
3241 nCgemLeft++;
3242 }
3243 }
3244 else
3245 if(useOld==0) (*iter)->setUse(-1);
3246 }
3247
3248 return nCgemLeft;
3249}
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
◆ bookTuple()
| StatusCode HoughFinder::bookTuple | ( | ) |
Definition at line 3603 of file Hough.cxx.
3603 {
3604 MsgStream log(msgSvc(), name());
3605
3607 if ( nt1 ){
3608 ntuple_hit= nt1;
3609 } else {
3611 if(ntuple_hit){
3612 ntuple_hit->addItem("hit_run", m_hit_run);
3613 ntuple_hit->addItem("hit_event", m_hit_event);
3614
3615 ntuple_hit->addItem("hit_nhit", m_hit_nhit, 0, 10000);
3616 ntuple_hit->addItem("hit_hitID", m_hit_nhit, m_hit_hitID);
3617 ntuple_hit->addItem("hit_hitType", m_hit_nhit, m_hit_hitType);
3618 ntuple_hit->addItem("hit_layer", m_hit_nhit, m_hit_layer);
3619 ntuple_hit->addItem("hit_wire", m_hit_nhit, m_hit_wire);
3620 ntuple_hit->addItem("hit_flag", m_hit_nhit, m_hit_flag);
3621 ntuple_hit->addItem("hit_halfCircle", m_hit_nhit, m_hit_halfCircle);
3622 ntuple_hit->addItem("hit_x", m_hit_nhit, m_hit_x);
3623 ntuple_hit->addItem("hit_y", m_hit_nhit, m_hit_y);
3624 ntuple_hit->addItem("hit_z", m_hit_nhit, m_hit_z);
3625 ntuple_hit->addItem("hit_drift", m_hit_nhit, m_hit_drift);
3626
3627 ntuple_hit->addItem("mcHit_hitID", m_hit_nhit, m_mcHit_hitID);
3628 ntuple_hit->addItem("mcHit_hitType", m_hit_nhit, m_mcHit_hitType);
3629 ntuple_hit->addItem("mcHit_layer", m_hit_nhit, m_mcHit_layer);
3630 ntuple_hit->addItem("mcHit_wire", m_hit_nhit, m_mcHit_wire);
3631 ntuple_hit->addItem("mcHit_flag", m_hit_nhit, m_mcHit_flag);
3632 ntuple_hit->addItem("mcHit_halfCircle", m_hit_nhit, m_mcHit_halfCircle);
3633 ntuple_hit->addItem("mcHit_x", m_hit_nhit, m_mcHit_x);
3634 ntuple_hit->addItem("mcHit_y", m_hit_nhit, m_mcHit_y);
3635 ntuple_hit->addItem("mcHit_z", m_hit_nhit, m_mcHit_z);
3636 ntuple_hit->addItem("mcHit_drift", m_hit_nhit, m_mcHit_drift);
3637 } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
3638 return StatusCode::FAILURE;
3639 }
3640 }
3641
3643 if ( nt2 ){
3644 ntuple_track = nt2;
3645 } else {
3647 if(ntuple_track){
3648 ntuple_track->addItem("trk_run", m_trk_run);
3649 ntuple_track->addItem("trk_event", m_trk_event);
3650 ntuple_track->addItem("trk_nTrack", m_trk_nTrack);
3651 ntuple_track->addItem("trk_trackID", m_trk_trackID);
3652 ntuple_track->addItem("trk_charge", m_trk_charge);
3653 ntuple_track->addItem("trk_flag", m_trk_flag);
3654 ntuple_track->addItem("trk_angle", m_trk_angle);
3655 ntuple_track->addItem("trk_rho", m_trk_rho);
3656 ntuple_track->addItem("trk_dAngle", m_trk_dAngle);
3657 ntuple_track->addItem("trk_dRho", m_trk_dRho);
3658 ntuple_track->addItem("trk_dTanl", m_trk_dTanl);
3659 ntuple_track->addItem("trk_dDz", m_trk_dDz);
3660 ntuple_track->addItem("trk_Xc", m_trk_Xc);
3661 ntuple_track->addItem("trk_Yc", m_trk_Yc);
3662 ntuple_track->addItem("trk_R", m_trk_R);
3663 ntuple_track->addItem("trk_dr", m_trk_dr);
3664 ntuple_track->addItem("trk_phi0", m_trk_phi0);
3665 ntuple_track->addItem("trk_kappa", m_trk_kappa);
3666 ntuple_track->addItem("trk_dz", m_trk_dz);
3667 ntuple_track->addItem("trk_tanl", m_trk_tanl);
3668 ntuple_track->addItem("trk_pxy", m_trk_pxy);
3669 ntuple_track->addItem("trk_px", m_trk_px);
3670 ntuple_track->addItem("trk_py", m_trk_py);
3671 ntuple_track->addItem("trk_pz", m_trk_pz);
3672 ntuple_track->addItem("trk_p", m_trk_p);
3673 ntuple_track->addItem("trk_phi", m_trk_phi);
3674 ntuple_track->addItem("trk_theta", m_trk_theta);
3675 ntuple_track->addItem("trk_cosTheta", m_trk_cosTheta);
3676 ntuple_track->addItem("trk_vx", m_trk_vx);
3677 ntuple_track->addItem("trk_vy", m_trk_vy);
3678 ntuple_track->addItem("trk_vz", m_trk_vz);
3679 ntuple_track->addItem("trk_vr", m_trk_vr);
3680 ntuple_track->addItem("trk_chi2", m_trk_chi2);
3681 ntuple_track->addItem("trk_fiTerm", m_trk_fiTerm);
3682 ntuple_track->addItem("trk_nhit", m_trk_nhit);
3683 ntuple_track->addItem("trk_ncluster", m_trk_ncluster);
3684 ntuple_track->addItem("trk_stat", m_trk_stat);
3685 ntuple_track->addItem("trk_ndof", m_trk_ndof);
3686 ntuple_track->addItem("trk_nster", m_trk_nster);
3687 ntuple_track->addItem("trk_nlayer", m_trk_nlayer);
3688 ntuple_track->addItem("trk_firstLayer", m_trk_firstLayer);
3689 ntuple_track->addItem("trk_lastLayer", m_trk_lastLayer);
3690 ntuple_track->addItem("trk_nCgemXClusters", m_trk_nCgemXClusters);
3691 ntuple_track->addItem("trk_nCgemVClusters", m_trk_nCgemVClusters);
3692
3693 ntuple_track->addItem("trk_nHot", m_trk_nHot, 0, 10000);
3694 ntuple_track->addItem("hot_hitID", m_trk_nHot, m_hot_hitID);
3695 ntuple_track->addItem("hot_hitType", m_trk_nHot, m_hot_hitType);
3696 ntuple_track->addItem("hot_layer", m_trk_nHot, m_hot_layer);
3697 ntuple_track->addItem("hot_wire", m_trk_nHot, m_hot_wire);
3698 ntuple_track->addItem("hot_flag", m_trk_nHot, m_hot_flag);
3699 ntuple_track->addItem("hot_halfCircle", m_trk_nHot, m_hot_halfCircle);
3700 ntuple_track->addItem("hot_x", m_trk_nHot, m_hot_x);
3701 ntuple_track->addItem("hot_y", m_trk_nHot, m_hot_y);
3702 ntuple_track->addItem("hot_z", m_trk_nHot, m_hot_z);
3703 ntuple_track->addItem("hot_drift", m_trk_nHot, m_hot_drift);
3704 ntuple_track->addItem("hot_residual", m_trk_nHot, m_hot_residual);
3705
3706 ntuple_track->addItem("mcHot_hitID", m_trk_nHot, m_mcHot_hitID);
3707 ntuple_track->addItem("mcHot_hitType", m_trk_nHot, m_mcHot_hitType);
3708 ntuple_track->addItem("mcHot_layer", m_trk_nHot, m_mcHot_layer);
3709 ntuple_track->addItem("mcHot_wire", m_trk_nHot, m_mcHot_wire);
3710 ntuple_track->addItem("mcHot_flag", m_trk_nHot, m_mcHot_flag);
3711 ntuple_track->addItem("mcHot_halfCircle", m_trk_nHot, m_mcHot_halfCircle);
3712 ntuple_track->addItem("mcHot_x", m_trk_nHot, m_mcHot_x);
3713 ntuple_track->addItem("mcHot_y", m_trk_nHot, m_mcHot_y);
3714 ntuple_track->addItem("mcHot_z", m_trk_nHot, m_mcHot_z);
3715 ntuple_track->addItem("mcHot_drift", m_trk_nHot, m_mcHot_drift);
3716
3717 ntuple_track->addItem("mcTrk_trackID", m_mcTrk_trackID);
3718 ntuple_track->addItem("mcTrk_charge", m_mcTrk_charge);
3719 ntuple_track->addItem("mcTrk_flag", m_mcTrk_flag);
3720 ntuple_track->addItem("mcTrk_angle", m_mcTrk_angle);
3721 ntuple_track->addItem("mcTrk_rho", m_mcTrk_rho);
3722 ntuple_track->addItem("mcTrk_dAngle", m_mcTrk_dAngle);
3723 ntuple_track->addItem("mcTrk_dRho", m_mcTrk_dRho);
3724 ntuple_track->addItem("mcTrk_dTanl", m_mcTrk_dTanl);
3725 ntuple_track->addItem("mcTrk_dDz", m_mcTrk_dDz);
3726 ntuple_track->addItem("mcTrk_Xc", m_mcTrk_Xc);
3727 ntuple_track->addItem("mcTrk_Yc", m_mcTrk_Yc);
3728 ntuple_track->addItem("mcTrk_R", m_mcTrk_R);
3729 ntuple_track->addItem("mcTrk_dr", m_mcTrk_dr);
3730 ntuple_track->addItem("mcTrk_phi0", m_mcTrk_phi0);
3731 ntuple_track->addItem("mcTrk_kappa", m_mcTrk_kappa);
3732 ntuple_track->addItem("mcTrk_dz", m_mcTrk_dz);
3733 ntuple_track->addItem("mcTrk_tanl", m_mcTrk_tanl);
3734 ntuple_track->addItem("mcTrk_pxy", m_mcTrk_pxy);
3735 ntuple_track->addItem("mcTrk_px", m_mcTrk_px);
3736 ntuple_track->addItem("mcTrk_py", m_mcTrk_py);
3737 ntuple_track->addItem("mcTrk_pz", m_mcTrk_pz);
3738 ntuple_track->addItem("mcTrk_p", m_mcTrk_p);
3739 ntuple_track->addItem("mcTrk_phi", m_mcTrk_phi);
3740 ntuple_track->addItem("mcTrk_theta", m_mcTrk_theta);
3741 ntuple_track->addItem("mcTrk_cosTheta", m_mcTrk_cosTheta);
3742 ntuple_track->addItem("mcTrk_vx", m_mcTrk_vx);
3743 ntuple_track->addItem("mcTrk_vy", m_mcTrk_vy);
3744 ntuple_track->addItem("mcTrk_vz", m_mcTrk_vz);
3745 ntuple_track->addItem("mcTrk_vr", m_mcTrk_vr);
3746 ntuple_track->addItem("mcTrk_chi2", m_mcTrk_chi2);
3747 ntuple_track->addItem("mcTrk_fiTerm", m_mcTrk_fiTerm);
3748 ntuple_track->addItem("mcTrk_nhit", m_mcTrk_nhit);
3749 ntuple_track->addItem("mcTrk_ncluster", m_mcTrk_ncluster);
3750 ntuple_track->addItem("mcTrk_stat", m_mcTrk_stat);
3751 ntuple_track->addItem("mcTrk_ndof", m_mcTrk_ndof);
3752 ntuple_track->addItem("mcTrk_nster", m_mcTrk_nster);
3753 ntuple_track->addItem("mcTrk_nlayer", m_mcTrk_nlayer);
3754 ntuple_track->addItem("mcTrk_firstLayer", m_mcTrk_firstLayer);
3755 ntuple_track->addItem("mcTrk_lastLayer", m_mcTrk_lastLayer);
3756 ntuple_track->addItem("mcTrk_nCgemXClusters", m_mcTrk_nCgemXClusters);
3757 ntuple_track->addItem("mcTrk_nCgemVClusters", m_mcTrk_nCgemVClusters);
3758
3759 ntuple_track->addItem("mcTrk_nHot", m_mcTrk_nHot, 0, 10000);
3760 ntuple_track->addItem("mcTrkHot_hitID", m_mcTrk_nHot, m_mcTrkHot_hitID);
3761 ntuple_track->addItem("mcTrkHot_hitType", m_mcTrk_nHot, m_mcTrkHot_hitType);
3762 ntuple_track->addItem("mcTrkHot_layer", m_mcTrk_nHot, m_mcTrkHot_layer);
3763 ntuple_track->addItem("mcTrkHot_wire", m_mcTrk_nHot, m_mcTrkHot_wire);
3764 ntuple_track->addItem("mcTrkHot_flag", m_mcTrk_nHot, m_mcTrkHot_flag);
3765 ntuple_track->addItem("mcTrkHot_halfCircle", m_mcTrk_nHot, m_mcTrkHot_halfCircle);
3766 ntuple_track->addItem("mcTrkHot_x", m_mcTrk_nHot, m_mcTrkHot_x);
3767 ntuple_track->addItem("mcTrkHot_y", m_mcTrk_nHot, m_mcTrkHot_y);
3768 ntuple_track->addItem("mcTrkHot_z", m_mcTrk_nHot, m_mcTrkHot_z);
3769 ntuple_track->addItem("mcTrkHot_drift", m_mcTrk_nHot, m_mcTrkHot_drift);
3770 } else {
3771 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/track" <<endmsg;
3772 return StatusCode::FAILURE;
3773 }
3774 }
3775
3777 if ( nt3 ){
3778 ntuple_event = nt3;
3779 } else {
3781 if(ntuple_event){
3782 ntuple_event->addItem("evt_run", m_evt_run);
3783 ntuple_event->addItem("evt_event", m_evt_event);
3784 ntuple_event->addItem("evt_nXCluster", m_evt_nXCluster);
3785 ntuple_event->addItem("evt_nVCluster", m_evt_nVCluster);
3786 ntuple_event->addItem("evt_nXVCluster", m_evt_nXVCluster);
3787 ntuple_event->addItem("evt_nCgemCluster", m_evt_nCGEMCluster);
3788 ntuple_event->addItem("evt_nAxialHit", m_evt_nAxialHit);
3789 ntuple_event->addItem("evt_nStereoHit", m_evt_nStereoHit);
3790 ntuple_event->addItem("evt_nODCHit", m_evt_nODCHit);
3791 ntuple_event->addItem("evt_nHit", m_evt_nHit);
3792
3793 ntuple_event->addItem("evt_nTrack", m_evt_nTrack, 0, 100);
3794 ntuple_event->addItem("evtTrk_trackID", m_evt_nTrack, m_evtTrk_trackID);
3795 ntuple_event->addItem("evtTrk_charge", m_evt_nTrack, m_evtTrk_charge);
3796 ntuple_event->addItem("evtTrk_flag", m_evt_nTrack, m_evtTrk_flag);
3797 ntuple_event->addItem("evtTrk_angle", m_evt_nTrack, m_evtTrk_angle);
3798 ntuple_event->addItem("evtTrk_rho", m_evt_nTrack, m_evtTrk_rho);
3799 ntuple_event->addItem("evtTrk_dAngle", m_evt_nTrack, m_evtTrk_dAngle);
3800 ntuple_event->addItem("evtTrk_dRho", m_evt_nTrack, m_evtTrk_dRho);
3801 ntuple_event->addItem("evtTrk_dTanl", m_evt_nTrack, m_evtTrk_dTanl);
3802 ntuple_event->addItem("evtTrk_dDz", m_evt_nTrack, m_evtTrk_dDz);
3803 ntuple_event->addItem("evtTrk_Xc", m_evt_nTrack, m_evtTrk_Xc);
3804 ntuple_event->addItem("evtTrk_Yc", m_evt_nTrack, m_evtTrk_Yc);
3805 ntuple_event->addItem("evtTrk_R", m_evt_nTrack, m_evtTrk_R);
3806 ntuple_event->addItem("evtTrk_dr", m_evt_nTrack, m_evtTrk_dr);
3807 ntuple_event->addItem("evtTrk_phi0", m_evt_nTrack, m_evtTrk_phi0);
3808 ntuple_event->addItem("evtTrk_kappa", m_evt_nTrack, m_evtTrk_kappa);
3809 ntuple_event->addItem("evtTrk_dz", m_evt_nTrack, m_evtTrk_dz);
3810 ntuple_event->addItem("evtTrk_tanl", m_evt_nTrack, m_evtTrk_tanl);
3811 ntuple_event->addItem("evtTrk_pxy", m_evt_nTrack, m_evtTrk_pxy);
3812 ntuple_event->addItem("evtTrk_px", m_evt_nTrack, m_evtTrk_px);
3813 ntuple_event->addItem("evtTrk_py", m_evt_nTrack, m_evtTrk_py);
3814 ntuple_event->addItem("evtTrk_pz", m_evt_nTrack, m_evtTrk_pz);
3815 ntuple_event->addItem("evtTrk_p", m_evt_nTrack, m_evtTrk_p);
3816 ntuple_event->addItem("evtTrk_phi", m_evt_nTrack, m_evtTrk_phi);
3817 ntuple_event->addItem("evtTrk_theta", m_evt_nTrack, m_evtTrk_theta);
3818 ntuple_event->addItem("evtTrk_cosTheta", m_evt_nTrack, m_evtTrk_cosTheta);
3819 ntuple_event->addItem("evtTrk_vx", m_evt_nTrack, m_evtTrk_vx);
3820 ntuple_event->addItem("evtTrk_vy", m_evt_nTrack, m_evtTrk_vy);
3821 ntuple_event->addItem("evtTrk_vz", m_evt_nTrack, m_evtTrk_vz);
3822 ntuple_event->addItem("evtTrk_vr", m_evt_nTrack, m_evtTrk_vr);
3823 ntuple_event->addItem("evtTrk_chi2", m_evt_nTrack, m_evtTrk_chi2);
3824 ntuple_event->addItem("evtTrk_fiTerm", m_evt_nTrack, m_evtTrk_fiTerm);
3825 ntuple_event->addItem("evtTrk_nhit", m_evt_nTrack, m_evtTrk_nhit);
3826 ntuple_event->addItem("evtTrk_ncluster", m_evt_nTrack, m_evtTrk_ncluster);
3827 ntuple_event->addItem("evtTrk_stat", m_evt_nTrack, m_evtTrk_stat);
3828 ntuple_event->addItem("evtTrk_ndof", m_evt_nTrack, m_evtTrk_ndof);
3829 ntuple_event->addItem("evtTrk_nster", m_evt_nTrack, m_evtTrk_nster);
3830 ntuple_event->addItem("evtTrk_nlayer", m_evt_nTrack, m_evtTrk_nlayer);
3831 ntuple_event->addItem("evtTrk_firstLayer", m_evt_nTrack, m_evtTrk_firstLayer);
3832 ntuple_event->addItem("evtTrk_lastLayer", m_evt_nTrack, m_evtTrk_lastLayer);
3833 ntuple_event->addItem("evtTrk_nCgemXClusters", m_evt_nTrack, m_evtTrk_nCgemXClusters);
3834 ntuple_event->addItem("evtTrk_nCgemVClusters", m_evt_nTrack, m_evtTrk_nCgemVClusters);
3835
3836 ntuple_event->addItem("mcEvtTrk_trackID", m_evt_nTrack, m_mcEvtTrk_trackID);
3837 ntuple_event->addItem("mcEvtTrk_charge", m_evt_nTrack, m_mcEvtTrk_charge);
3838 ntuple_event->addItem("mcEvtTrk_flag", m_evt_nTrack, m_mcEvtTrk_flag);
3839 ntuple_event->addItem("mcEvtTrk_angle", m_evt_nTrack, m_mcEvtTrk_angle);
3840 ntuple_event->addItem("mcEvtTrk_rho", m_evt_nTrack, m_mcEvtTrk_rho);
3841 ntuple_event->addItem("mcEvtTrk_dAngle", m_evt_nTrack, m_mcEvtTrk_dAngle);
3842 ntuple_event->addItem("mcEvtTrk_dRho", m_evt_nTrack, m_mcEvtTrk_dRho);
3843 ntuple_event->addItem("mcEvtTrk_dTanl", m_evt_nTrack, m_mcEvtTrk_dTanl);
3844 ntuple_event->addItem("mcEvtTrk_dDz", m_evt_nTrack, m_mcEvtTrk_dDz);
3845 ntuple_event->addItem("mcEvtTrk_Xc", m_evt_nTrack, m_mcEvtTrk_Xc);
3846 ntuple_event->addItem("mcEvtTrk_Yc", m_evt_nTrack, m_mcEvtTrk_Yc);
3847 ntuple_event->addItem("mcEvtTrk_R", m_evt_nTrack, m_mcEvtTrk_R);
3848 ntuple_event->addItem("mcEvtTrk_dr", m_evt_nTrack, m_mcEvtTrk_dr);
3849 ntuple_event->addItem("mcEvtTrk_phi0", m_evt_nTrack, m_mcEvtTrk_phi0);
3850 ntuple_event->addItem("mcEvtTrk_kappa", m_evt_nTrack, m_mcEvtTrk_kappa);
3851 ntuple_event->addItem("mcEvtTrk_dz", m_evt_nTrack, m_mcEvtTrk_dz);
3852 ntuple_event->addItem("mcEvtTrk_tanl", m_evt_nTrack, m_mcEvtTrk_tanl);
3853 ntuple_event->addItem("mcEvtTrk_pxy", m_evt_nTrack, m_mcEvtTrk_pxy);
3854 ntuple_event->addItem("mcEvtTrk_px", m_evt_nTrack, m_mcEvtTrk_px);
3855 ntuple_event->addItem("mcEvtTrk_py", m_evt_nTrack, m_mcEvtTrk_py);
3856 ntuple_event->addItem("mcEvtTrk_pz", m_evt_nTrack, m_mcEvtTrk_pz);
3857 ntuple_event->addItem("mcEvtTrk_p", m_evt_nTrack, m_mcEvtTrk_p);
3858 ntuple_event->addItem("mcEvtTrk_phi", m_evt_nTrack, m_mcEvtTrk_phi);
3859 ntuple_event->addItem("mcEvtTrk_theta", m_evt_nTrack, m_mcEvtTrk_theta);
3860 ntuple_event->addItem("mcEvtTrk_cosTheta", m_evt_nTrack, m_mcEvtTrk_cosTheta);
3861 ntuple_event->addItem("mcEvtTrk_vx", m_evt_nTrack, m_mcEvtTrk_vx);
3862 ntuple_event->addItem("mcEvtTrk_vy", m_evt_nTrack, m_mcEvtTrk_vy);
3863 ntuple_event->addItem("mcEvtTrk_vz", m_evt_nTrack, m_mcEvtTrk_vz);
3864 ntuple_event->addItem("mcEvtTrk_vr", m_evt_nTrack, m_mcEvtTrk_vr);
3865 ntuple_event->addItem("mcEvtTrk_chi2", m_evt_nTrack, m_mcEvtTrk_chi2);
3866 ntuple_event->addItem("mcEvtTrk_fiTerm", m_evt_nTrack, m_mcEvtTrk_fiTerm);
3867 ntuple_event->addItem("mcEvtTrk_nhit", m_evt_nTrack, m_mcEvtTrk_nhit);
3868 ntuple_event->addItem("mcEvtTrk_ncluster", m_evt_nTrack, m_mcEvtTrk_ncluster);
3869 ntuple_event->addItem("mcEvtTrk_stat", m_evt_nTrack, m_mcEvtTrk_stat);
3870 ntuple_event->addItem("mcEvtTrk_ndof", m_evt_nTrack, m_mcEvtTrk_ndof);
3871 ntuple_event->addItem("mcEvtTrk_nster", m_evt_nTrack, m_mcEvtTrk_nster);
3872 ntuple_event->addItem("mcEvtTrk_nlayer", m_evt_nTrack, m_mcEvtTrk_nlayer);
3873 ntuple_event->addItem("mcEvtTrk_firstLayer", m_evt_nTrack, m_mcEvtTrk_firstLayer);
3874 ntuple_event->addItem("mcEvtTrk_lastLayer", m_evt_nTrack, m_mcEvtTrk_lastLayer);
3875 ntuple_event->addItem("mcEvtTrk_nCgemXClusters", m_evt_nTrack, m_mcEvtTrk_nCgemXClusters);
3876 ntuple_event->addItem("mcEvtTrk_nCgemVClusters", m_evt_nTrack, m_mcEvtTrk_nCgemVClusters);
3877 } else {
3878 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/event" <<endmsg;
3879 return StatusCode::FAILURE;
3880 }
3881 }
3882 return StatusCode::SUCCESS;
3883}
Referenced by initialize().
◆ checkHot()
| int HoughFinder::checkHot | ( | vector< HoughTrack > & | trackVector | ) |
Definition at line 1655 of file Hough.cxx.
1656{
1657 //m_debug=true;
1658 int nDelete(0);
1659 if(trackVector.size()==0)return 0;
1660 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1661
1662 vector<HoughTrack>::iterator trkIT1 = trackVector.end();
1663 //vector<HoughTrack>::iterator trkIT1 = trackVector.rbegin();
1664 bool loop = true;
1665
1666 while(loop)
1667 {
1668 trkIT1--;
1669 if(trkIT1==trackVector.begin()) loop=false;
1670 if((trkIT1)->getFlag() == 1) continue;
1671 //trkIT1->dropRedundentCgemXHits();
1672 //cout<<"hit size "<<hot1->size()<<endl;
1673 int nHits = trkIT1->getVecHitPnt().size();
1674 if(nHits<3){
1675 (trkIT1)->setFlag(1);
1676 trkIT1->releaseSelHits();
1677 if(m_debug==1) {
1678 cout<<"too less hits ("<<nHits<<") in track "<<(trkIT1)->getTrkID()<<endl;
1679 }
1680 continue;
1681 }
1682 int nShared = trkIT1->getNHitsShared();
1683
1684 //cout<<"shareHit = "<<shareHit<<endl;
1685 if(nShared>m_shareHitRate*nHits){
1686 trkIT1->setFlag(1);
1687 trkIT1->releaseSelHits();
1688 nDelete++;
1689 if(m_debug==1) {
1690 cout<<"too many shared hits ("<<nShared<<"/"<<nHits
1691 <<") in track "<<(trkIT1)->getTrkID()
1692 <<endl;
1693 }
1694 }
1695 }
1696 //m_debug=false;
1697 return nDelete;
1698}
Referenced by execute().
◆ checkTrack()
| int HoughFinder::checkTrack | ( | vector< HoughTrack > & | trackVector | ) |
Definition at line 1700 of file Hough.cxx.
1701{
1702 int nDelete(0);
1703 double dDr(999);
1704 double dPhi0(999);
1705 double dKappa(999);
1706 double dDz(999);
1707 double dTanl(999);
1708
1709 //m_debug=true;
1710
1711 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1712 for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
1713 if((trkIT1)->getFlag() == 1)continue;
1714 int charge = (trkIT1)->getCharge();
1715 double dr = (trkIT1)->dr();
1716 double phi0 = (trkIT1)->phi0();
1717 double kappa = (trkIT1)->kappa();
1718 double dz = (trkIT1)->dz();
1719 double tanl = (trkIT1)->tanl();
1720 for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
1721 if((trkIT2)->getFlag() == 1)continue;
1722 bool sameTrack(false);
1723 dDr = fabs(dr - (trkIT2)->dr());
1724 dPhi0 = fabs(phi0 - (trkIT2)->phi0());
1725 if((trkIT2)->getCharge() == charge){
1726 dKappa = fabs(kappa - (trkIT2)->kappa());
1727 dTanl = fabs(tanl - (trkIT2)->tanl());
1728 }else{
1729 dKappa = fabs(fabs(kappa) - fabs((trkIT2)->kappa()));
1730 dTanl = fabs(fabs(tanl) - fabs((trkIT2)->tanl()));
1731 }
1732 double r = fabs(dz)<fabs((trkIT2)->dz())?(trkIT1)->radius():(trkIT2)->radius();
1733 double k = dz<fabs((trkIT2)->dz())?tanl:(trkIT2)->tanl();
1734 dDz = fabs(dz - (trkIT2)->dz());
1736 dDz = fabs(dDz - nCircle*(2*M_PI*r*k));
1738 if(dPhi0<m_phi0Cut&&dKappa<m_kappaCut&&dTanl<m_tanlCut) sameTrack=true;
1739 //sameTrack = sameTrack&&(dDr<m_drCut);//FIXME
1740 //sameTrack = sameTrack&&(dPhi0<m_phi0Cut);//FIXME
1741 //sameTrack = sameTrack&&(dKappa<m_kappaCut);//FIXME
1742 //sameTrack = sameTrack&&(dDz<m_dzCut);//FIXME
1743 //sameTrack = sameTrack&&(dTanl<m_tanlCut);//FIXME
1744 if(sameTrack){
1745 //delete *trkIT2;
1746 //trackVector.erase(trkIT2);
1747 //trkIT2--;
1748 (trkIT2)->setFlag(1);
1749 if(m_debug==5) {
1750 cout<<"dDr, dPhi0, dKappa, dDz, dTanl "
1751 <<setw(10)<<dDr
1752 <<setw(10)<<dPhi0
1753 <<setw(10)<<dKappa
1754 <<setw(10)<<dDz
1755 <<setw(10)<<dTanl
1756 <<endl;
1757 cout<<"similar track "<<trkIT2->getTrkID()<<" is removed from hough track list"<<endl;
1758 }
1759
1760 // --- associate the hits of track 2 to track 1 if within pi/2 // to-do
1761 nDelete++;
1762 }
1763 }
1764 }
1765 if(m_debug==5) cout<<"checkTrack(): "<<nDelete<<" similar tracks are removed from hough track list"<<endl;
1766 //m_debug=false;
1767 return nDelete;
1768}
◆ clearMemory()
| void HoughFinder::clearMemory | ( | ) |
Definition at line 3212 of file Hough.cxx.
3213{
3214 //cout<<"m_mdcHitCol.size()"<<m_mdcHitCol.size()<<endl;
3215 vector<MdcHit*>::iterator imdcHit = m_mdcHitCol.begin();
3216 for(;imdcHit != m_mdcHitCol.end();imdcHit++){
3217 delete (*imdcHit);
3218 }
3219
3220 for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin();trkIT!=m_houghTrackList.end();trkIT++){
3221 trkIT->clearMemory();
3222 //TrkRecoTrk* trkRecoTrk = trkIT->getTrkRecoTrk();
3223 //delete trkRecoTrk;
3224 }
3225}
Referenced by execute().
◆ clearTrack()
| void HoughFinder::clearTrack | ( | vector< HoughTrack > & | trackVector | ) |
◆ execute()
| StatusCode HoughFinder::execute | ( | ) |
Definition at line 299 of file Hough.cxx.
300{
301 MsgStream log(msgSvc(), name());
302 log << MSG::INFO << "HoughFinder::execute()" << endreq;
303 cout.precision(6);
304
305
306 //cout<<"line "<<__LINE__<<endl;
307 // --- Get event header
308 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
309 if(!eventHeader){
310 log << MSG::ERROR << "Could not find Event Header" << endreq;
311 return StatusCode::FAILURE;
312 }
313 m_run = eventHeader->runNumber();
314 m_event = eventHeader->eventNumber();
315 m_totEvtProcessed++; //if(m_totEvtProcessed%100==0) cout<<"HoughFinder: "<<m_totEvtProcessed<<" events processed."<<endl;
316
317 myDotsHelixFitter.updateBField();
318
319 //cout<<"line "<<__LINE__<<endl;
320
321 // --- event start time
322 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
323 if(!evTimeCol){
324 log << MSG::ERROR << "Could not retrieve RecEsTimeCol" << endreq;
325 return StatusCode::FAILURE;
326 }
327 m_bunchT0 = 0;
328 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
329 if (iter_evt != evTimeCol->end()){
330 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
331 //cout<<"this event T0 = "<<m_bunchT0<<endl;
332 }
333
334 // --- clear, register and store reconstructed track
335 registerTrack(m_trackList_tds,m_hitList_tds);
336
337
338 //cout<<"line "<<__LINE__<<endl;
339 if(m_filter){
340 ifstream lowPt_Evt;
341 lowPt_Evt.open(m_evtFile.c_str());
342 vector<int> evtlist;
343 int evtNum;
344 while( lowPt_Evt >> evtNum) {
345 evtlist.push_back(evtNum);
346 }
347 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),m_event);
348 if(iter_evt == evtlist.end()){
349 setFilterPassed(false);
350 return StatusCode::SUCCESS;
351 }
352 else{
353 setFilterPassed(true);
354 }
355 }
356 //if(m_event<2000||m_event>3000){
357 //setFilterPassed(false);
358 //return StatusCode::SUCCESS;
359 //}
360
361 if(m_debug==5)
362 {
363 //cout<<endl;
364 cout<<"===================================================================================================="<<endl;
365 cout<<"run:"<<m_run<<" , event: "<<m_event<<endl;
367 //cout<<m_event<<" , ";
368 //cout<<endl;
369 }
370
371 //cout<<"line "<<__LINE__<<endl;
372 // --- prepare hits for reconstruction
374
375 //cout<<"line "<<__LINE__<<endl;
376 int nMcHit;
377 if(m_run<0&&m_mcTruth){
379 //cout<<"nMcHit "<<nMcHit<<endl;
380 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
381 getMcParticleCol();
382 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
386 //halfCircle(m_mcTrackCol,m_mcHitCol);
387 }
388
389 //cout<<"line "<<__LINE__<<endl;
390 searchCircle();
391 //cout<<"ntracks:"<<m_houghTrackList.size()<<endl;
392 //cout<<"line "<<__LINE__<<endl;
394 //cout<<"line "<<__LINE__<<endl;
396 //updateXHitsStatistics(m_houghTrackList);
397 //checkTrack(m_houghTrackList);
398 //cout<<"line "<<__LINE__<<endl;
399 associateVHits();
400 //cout<<"line "<<__LINE__<<endl;
401 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
403 else
405 else
407 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
411 //cout<<"line "<<__LINE__<<endl;
412
413 if(m_fillNTuple)dumpHit();
414 if(m_fillNTuple==1)dumpHoughTrack();
415 if(m_fillNTuple==1)dumpHoughEvent();
416 if(m_fillNTuple==2)dumpTdsTrack();
417 if(m_fillNTuple==2)dumpTdsEvent();
418 //cout<<"line "<<__LINE__<<endl;
419
420 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
422 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
426
427 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
429 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
433
434 clearTrack(m_houghTrackList);
435 //cout<<"line "<<__LINE__<<endl;
436
437 clearMemory();
438 //cout<<"line "<<__LINE__<<endl;
439 if(m_debug==5) cout<<endl;
440
441 return StatusCode::SUCCESS;
442}// <-- execute()
double bFieldNominal() const
int storeRecTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:1922
void solveSharedHits(vector< HoughHit * > &hitList)
Definition Hough.cxx:3252
int storeTrack(vector< HoughTrack > &trackVector, RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition Hough.cxx:1881
int storeTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:2207
StatusCode registerTrack(RecMdcTrackCol *&trackList_tds, RecMdcHitCol *&hitList_tds)
Definition Hough.cxx:1828
◆ fillHistogram() [1/3]
| int HoughFinder::fillHistogram | ( | HoughHit * | hit, |
| TH2D * | hitMap, | ||
| HoughTrack * | trkCandi, | ||
| int | vote ) |
Definition at line 1037 of file Hough.cxx.
1038{
1040 int nBin = hitMap->GetXaxis()->GetNbins();
1041 double xMin = hitMap->GetXaxis()->GetXmin();
1042 double xMax = hitMap->GetXaxis()->GetXmax();
1043 double yMin = hitMap->GetYaxis()->GetXmin();
1044 double yMax = hitMap->GetYaxis()->GetXmax();
1045 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;// step size
1047
1049 int nPoint(0);
1050 if(flag==0) // X-view hits
1051 {
1053 if(!firstHalf) return 0;
1057 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
1058
1059 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
1061 double Y1 = 2*yHit/denominator;
1063 double Y2 = 2*yHit/denominator;
1065
1066 while(x<xMax){
1071 double hitOnCircle1 = charge*slope1*y1;
1072 double hitOnCircle2 = charge*slope2*y2;
1073
1075 cut = yMin<y1 && y1<yMax;
1076 //cut = cut && hitOnCircle1 <=0;
1078 hitMap->Fill(x,y1);
1079 nPoint++;
1080 }
1081
1082 cut = yMin<y2 && y2<yMax;
1083 //cut = cut && hitOnCircle2 <= 0;
1085 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
1086 //int bin = hitMap->FindBin(x,y2);
1087 //if(hitMap->GetBinContent(bin)>0)continue;
1088 hitMap->Fill(x,y2);
1089 nPoint++;
1090 }
1091 x += dx;
1092 }
1093 }else{ // V-view hits
1094 //if(hit.getPairHit()==NULL)return nPoint;// FIXME
1095 //if(hit.getPairHit()->getHalfCircle()!=1)return nPoint;// FIXME
1097 if(nTangency==0){
1098 //hit.setUse(-1);// why? FIXME
1099 }
1100 else
1101 {
1102 vector<HoughHit::S_Z> sz = hit->getSZ();
1103 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1105 {
1108 while(x<xMax)
1109 {
1111 //y = S*cos(x) + Z*sin(x);
1113 cut = yMin<y && y<yMax;
1115 {
1116 hitMap->Fill(x,y);
1117 nPoint++;
1118 }
1119 x += dx;
1120 }
1121 }
1122 }
1123 //cout<<setw(8)<<""<<"nTangency="<<nTangency<<", nPointFilled="<<nPoint<<endl;
1124 }
1125 return nPoint;
1126}
*********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 956 of file Hough.cxx.
957{
958 //cout<<"fillHistogram "<<__LINE__<<endl;
959 int nBin = hitMap->GetXaxis()->GetNbins();
960 double xMin = hitMap->GetXaxis()->GetXmin();
961 double xMax = hitMap->GetXaxis()->GetXmax();
962 double yMin = hitMap->GetYaxis()->GetXmin();
963 double yMax = hitMap->GetYaxis()->GetXmax();
964 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
966
968 int nPoint(0);
969 if(flag==0){
973 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
974
975 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
977 double Y1 = 2*yHit/denominator;
979 double Y2 = 2*yHit/denominator;
981
982 while(x<xMax){
987 double hitOnCircle1 = charge*slope1*y1;
988 double hitOnCircle2 = charge*slope2*y2;
989
991 cut = yMin<y1 && y1<yMax;
994 hitMap->Fill(x,y1);
995 nPoint++;
996 }
997
999 cut = yMin<y2 && y2<yMax;
1002 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
1003 //int bin = hitMap->FindBin(x,y2);
1004 //if(hitMap->GetBinContent(bin)>0)continue;
1005 hitMap->Fill(x,y2);
1006 nPoint++;
1007 }
1008 }
1009
1010 x += dx;
1011 }
1012 }else{
1013 vector<HoughHit::S_Z> sz = hit->getSZ();
1014 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1018 double dy = -S*dx;
1019 double y = -dy + Z;
1020 while(x<xMax){
1021 x += dx;
1022 y += dy;
1023 //y = S*cos(x) + Z*sin(x);
1025 cut = yMin<y && y<yMax;
1027 hitMap->Fill(x,y);
1028 nPoint++;
1029 }
1030 }
1031 }
1032 }
1033 //cout<<"fillHistogram "<<__LINE__<<endl;
1034 return nPoint;
1035}
Referenced by fillHistogram().
◆ fillHistogram() [3/3]
| int HoughFinder::fillHistogram | ( | vector< HoughHit * > & | hitList, |
| TH2D * | hitMap, | ||
| int | charge, | ||
| int | vote, | ||
| HoughTrack * | trkCandi = NULL ) |
Definition at line 1128 of file Hough.cxx.
1129{
1130 bool nearStereoHits= false;
1132
1133 int nHitsFilled = 0;
1134 m_nVClusterOnSZmap=0;
1135 m_VHoughHitListOnSZmap.clear();
1136 std::vector<HoughHit*>::iterator iter = hitList.begin();
1138 {
1139 int used = (*iter)->getUse();
1140 if(used==0) // too tight? FIXME
1141 {
1142 int nPoint = 0;
1143 //(*iter)->print();
1145 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1146 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1147 if(!nearStereoHits&&(*iter)->getHitType()==HoughHit::MDCHIT&&(*iter)->getFlag()!=0) continue;// if no axial hits around, the stereo hits should not be filled
1148
1150 }
1151 else{
1152 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1153 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1155 }
1156 if(nPoint>0) {
1157 nHitsFilled++;
1159 {
1161 m_VHoughHitListOnSZmap.push_back(*iter);
1162 }
1163 }
1164 //cout<<nPoint<<" nPoint filled"<<(*iter)->print();
1165 }
1166 }
1167 return nHitsFilled;
1168}
int fillHistogram(HoughHit *hit, TH2D *hitMap, int charge, int vote)
Definition Hough.cxx:956
◆ finalize()
| StatusCode HoughFinder::finalize | ( | ) |
◆ getHoughTrkIt()
| vector< HoughTrack >::iterator HoughFinder::getHoughTrkIt | ( | vector< HoughTrack > & | houghTrkList, |
| int | trkId ) |
Definition at line 3297 of file Hough.cxx.
3298{
3299 vector<HoughTrack>::iterator it = houghTrkList.begin();
3300 for(; it!=houghTrkList.end(); it++)
3301 {
3302 if(it->getTrkID()==trkId) break;
3303 }
3304 return it;
3305}
Referenced by solveSharedHits().
◆ getMcHitCol()
| int HoughFinder::getMcHitCol | ( | ) |
— put match RecCgemCluster into CemgMcHit
Definition at line 523 of file Hough.cxx.
524{
525 MsgStream log(msgSvc(), name());
526 int nMcHit(0);
527 if(!(m_mcHitCol.empty()))m_mcHitCol.clear();
528
529 // --- get cgemMcHit
530 if(m_cgem){
531 SmartDataPtr<CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
532 if (!cgemMcHitCol){
533 log << MSG::ERROR << "Could not retrieve CgemMcHitCol" << endreq;
534 }else{
535 CgemMcHitCol::iterator cgemMcHitIt = cgemMcHitCol->begin();
536 for(;cgemMcHitIt != cgemMcHitCol->end();cgemMcHitIt++,nMcHit++){
538 HoughHit hit(cgemMcHit,m_bunchT0,nMcHit);
539 m_mcHitCol.push_back(hit);
540 HoughHit* phit = &(*m_mcHitCol.rbegin());
541
542 // --- match fired strip ID between CemgMcHit and RecCgemCluster
543 vector<int> xFireStripID;
544 vector<int> vFireStripID;
545 int sheet;
547 for(int ii = 0; ii < identifier.GetSize(); ii++){
548 //const Identifier ident = Identifier(identifier.GetAt(ii));
552 sheet = CgemID::sheet(ident);
553 //cout
554 //<< " hitID:" << phit->getHitID()
555 //<< " layerID=" << CgemID::layer(ident)
556 //<< " sheetID=" << CgemID::sheet(ident)
557 //<< " isXstrip=" << CgemID::is_xstrip(ident)
558 //<< " stripID=" << CgemID::strip(ident)
559 //<<endl;
560 }
561 sort(xFireStripID.begin(),xFireStripID.end());
562 sort(vFireStripID.begin(),vFireStripID.end());
563 //vector<int>::iterator uniqueX = unique(xFireStripID.begin(),xFireStrip ID.end());
564 //vector<int>::iterator uniqueV = unique(vFireStripID.begin(),vFireStrip ID.end());
565 //xFireStripID.erase(uniqueX,xFireStripID.end());
566 //vFireStripID.erase(uniqueV,vFireStripID.end());
567 //for(vector<int>::iterator i = vFireStripID.begin();i!=vFireStripID.end();i++)cout<<*i<<endl;
568 //cout<<endl;
569
570 bool findX(false), findV(false);
578 //if(recCgemCluster->getsheetid() != sheet)continue;
579 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
580 if(!findX){
581 findX = recCgemCluster->getflag() == 0
582 && xFireStripID.size() > 0
583 && recCgemCluster->getclusterflagb() == *(xFireStripID.begin())
584 && recCgemCluster->getclusterflage() == *(xFireStripID.rbegin());
585 if(findX){
586 xCluster = recCgemCluster;
587 hitIt->setPairHit(phit);
588 }
589 }
590 if(!findV){
591 findV = recCgemCluster->getflag() == 1
592 && vFireStripID.size() > 0
593 && recCgemCluster->getclusterflagb() == *(vFireStripID.begin())
594 && recCgemCluster->getclusterflage() == *(vFireStripID.rbegin());
595 if(findV){
596 vCluster = recCgemCluster;
597 hitIt->setPairHit(phit);
598 }
599 }
600
601 }
602
603 /// --- put match RecCgemCluster into CemgMcHit
604 if(findX&&findV){
610 //if(recCgemCluster->getsheetid() != sheet)continue;
611 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
614 phit->setCgemCluster(recCgemCluster);
617 phit->setPairHit(&(*hitIt));
618 hitIt->setPairHit(phit);
619 }
620 }
621 }
622 else if(findX&&!findV){
623 phit->setCgemCluster(xCluster);
626 }
627 else if(!findX&&findV){
628 phit->setCgemCluster(vCluster);
631 }
632 }
633 }
634 }
635
636 // --- get mdcMcHit
637 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
638 if(!mdcMcHitCol){
639 log << MSG::ERROR << "Could not retrieve MdcMcHitCol" << endreq;
640 }else{
641 MdcMcHitCol::iterator mdcMcHitIt = mdcMcHitCol->begin();
642 for(;mdcMcHitIt != mdcMcHitCol->end();mdcMcHitIt++,nMcHit++){
644 HoughHit hit(mdcMcHit,m_bunchT0,nMcHit);
645 m_mcHitCol.push_back(hit);
646 HoughHit* phit = &(*m_mcHitCol.rbegin());
651 phit->setDigi(mdcDigi);
652 phit->setPairHit(&(*hitIt));
653 hitIt->setPairHit(phit);
654 }
655 }
656 }
657 }
658
661 if(hitIt->getLayer()!=mcHitIt->getLayer())continue;
663 if(mcHitIt->getCgemCluster()==NULL)continue;
664 int recClusterID(-1);
665 int mcClusterID(-2);
666 int recFlag = hitIt->getCgemCluster()->getflag();
667 int mcFlag = mcHitIt->getCgemCluster()->getflag();;
668 if(recFlag==mcFlag){
669 recClusterID = hitIt->getCgemCluster()->getclusterid();
670 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
671 }
672
673 if(recFlag==1&&mcFlag==2){
674 recClusterID = hitIt->getCgemCluster()->getclusterid();
675 mcClusterID = mcHitIt->getCgemCluster()->getclusterflage();
676 }
677
678 if(recFlag==0&&mcFlag==2){
679 recClusterID = hitIt->getCgemCluster()->getclusterid();
680 mcClusterID = mcHitIt->getCgemCluster()->getclusterflagb();
681 }
682
683 if(recFlag==2&&mcFlag==1){
684 recClusterID = hitIt->getCgemCluster()->getclusterflage();
685 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
686 }
687
688 if(recFlag==2&&mcFlag==0){
689 recClusterID = hitIt->getCgemCluster()->getclusterflagb();
690 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
691 }
692
693 if(recClusterID==mcClusterID){
694 hitIt->setPairHit(&(*mcHitIt));
695 hitIt->setCgemMcHit(mcHitIt->getCgemMcHit());
696 //mcHitIt->addPairHit(&(*hitIt));
697
698 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" "<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
699 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
700 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" ";
701 //cout<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
702 //cout<<recFlag<<" "<<mcFlag<<" ";
703 //cout<<recClusterID<<" "<<mcClusterID<<endl;
704
705 //cout<<hitIt->getHitPosition();
706 //cout<<endl;
707 //cout<<mcHitIt->getHitPosition();
708 //cout<<endl;
709 //cout<<endl;
710 }
711 }
712 else if(hitIt->getHitType() == HoughHit::MDCHIT && mcHitIt->getHitType() == HoughHit::MDCMCHIT){
713 if(mcHitIt->getDigi()==NULL)continue;
714 //if(hitIt->getLayer() == mcHitIt->getLayer() && hitIt->getWire() == mcHitIt->getWire())
715 if(hitIt->getDigi()->identify() == mcHitIt->getDigi()->identify())
716 {
717 hitIt->setPairHit(&(*mcHitIt));
718 hitIt->setMdcMcHit(mcHitIt->getMdcMcHit());
719 //mcHitIt->addPairHit(&(*hitIt));
720 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<endl;
721 //cout<<mcHitIt->getLayer()<<" "<<mcHitIt->getWire()<<endl;
722 //cout<<endl;
723 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
724 }
725 }
726 }
727 }
728 //for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
729 //if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
730 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<" ";
731 //if(hitIt->getPairHit()!=NULL){
732 //cout<<hitIt->getPairHit()->getHitID()<<" ";
733 //cout<<hitIt->getPairHit()->getLayer()<<" "<<hitIt->getPairHit()->getWire()<<" ";
734 //}
735 //cout<<endl;
736 //}
737
738 return nMcHit;
739}
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 741 of file Hough.cxx.
741 {
742 const int maxCell[43] = {40,44,48,56, 64,72,80,80,
743 76,76,88,88, 100,100,112,112, 128,128,140,140,
744 160,160,160,160, 176,176,176,176, 208,208,208,208, 240,240,240,240,
745 256,256,256,256, 288,288,288 };
746 m_mcTrackCol.clear();
747 MsgStream log(msgSvc(), name());
748 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
749 if (!mcParticleCol) {
750 log << MSG::ERROR << "Could not find McParticle" << endreq;
751 }else{
752 bool psipDecay(false);
753 McParticleCol::iterator iter_mc = mcParticleCol->begin();
754 for (;iter_mc != mcParticleCol->end(); iter_mc++){
755 int trackIndex = (*iter_mc)->trackIndex();
756 int pid = (*iter_mc)->particleProperty();
757 bool primaryParticle = (*iter_mc)->primaryParticle();
758 bool leafParticle = (*iter_mc)->leafParticle();
759 bool decayFromGenerator = (*iter_mc)->decayFromGenerator();
760 bool decayInFlight = (*iter_mc)->decayInFlight();
761 HepLorentzVector initialPosition = (*iter_mc)->initialPosition();
762 HepLorentzVector initialFourMomentum = (*iter_mc)->initialFourMomentum();
763 //HepLorentzVector finalPosition = (*iter_mc)->finalPosition();
764 //HepLorentzVector finalFourMomentum = (*iter_mc)->finalFourMomentum();
765 unsigned int statusFlags = (*iter_mc)->statusFlags();
766 //string process = (*iter_mc)->getProcess();
767 //if(primaryParticle)continue;
768 //if(!decayFromGenerator)continue;
769 //if(pid==100443)psipDecay = true;
770 //if(!psipDecay) continue;
771 //if(!(fabs(pid)==11||fabs(pid)==211))continue;
772
773 string name;
774 int charge(0);
775 if( pid == 0 ) {
776 //cout << "Wrong particle id" <<endl;
777 continue;
778 }else{
779 IPartPropSvc* p_PartPropSvc;
780 static const bool CREATEIFNOTTHERE(true);
781 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
782 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
783 log << MSG::ERROR << "Could not initialize Particle Properties Service" << endreq;
784 }
785 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
786 if( p_particleTable->particle(pid) ){
787 name = p_particleTable->particle(pid)->name();
788 charge = p_particleTable->particle(pid)->charge();
789 }else if( p_particleTable->particle(-pid) ){
790 name = "anti " + p_particleTable->particle(-pid)->name();
791 charge = (-1)*p_particleTable->particle(-pid)->charge();
792 }
793 }
794
795 //cout
796 //<<setw(8)<<trackIndex
797 //<<setw(12)<<name
798 //<<setw(8)<<pid
799 //<<setw(8)<<primaryParticle
800 //<<setw(8)<<leafParticle
801 //<<setw(8)<<decayFromGenerator
802 //<<setw(8)<<decayInFlight
803 //<<setw(8)<<statusFlags
804 //<<endl;
805
806 HoughTrack track(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
807 //HoughTrack* track = new HoughTrack(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
808 HepPoint3D pivot(0,0,0);
809 track.pivot(pivot);
810
811 int nHot(0);
812 vector<HoughHit*> hotList;
813 vector<HoughHit>& hitList = m_mcHitCol;
815 vector<int> trkID = hitIt->getTrkID();
816 if(trkID[0]==track.getTrkID()){
818 track.addHitPnt(&(*hitIt));
819 track.addVecStereoHitPnt(&(*hitIt));
820 hotList.push_back(&(*hitIt));
821 }
823 if(hitIt->getFlag()==0)track.addHitPnt(&(*hitIt));
824 else track.addVecStereoHitPnt(&(*hitIt));
825 hotList.push_back(&(*hitIt));
826 }
827 nHot++;
828 }
829 }
830 if(nHot>0)m_mcTrackCol.push_back(track);
831
832 ///*
833 bool classifyHotByHotLayer(false);
834 bool classifyHotByTrackCenter(true);
835 if(classifyHotByHotLayer){
836 int lastHitLayer(0);
837 int deltaLayer(0);
838 int halfCircle(0);
839 HoughHit* lastLayerHit;
840 vector<HoughHit*> hitOnLayer;
841 int nMdcHot(0);
842 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
844 if(nMdcHot==0){//the first hit on track
845 hitOnLayer.clear();
846 hitOnLayer.push_back(*hitIt);
847 lastLayerHit = *hitIt;
848 lastHitLayer = (*hitIt)->getLayer();
849 deltaLayer = 0;
850 halfCircle = 1;
851 }else{
852 if((*hitIt)->getLayer()==lastHitLayer){//find hits on the same layer
853 hitOnLayer.push_back(*hitIt);
854 lastHitLayer = (*hitIt)->getLayer();
855 }else{//another layer hit
856 if(deltaLayer*((*hitIt)->getLayer()-lastHitLayer)>=0){//not the turn layer
858 (*iter)->setHalfCircle(halfCircle);
859 lastLayerHit = *iter;
860 }
861 }else{//the turn layer
862 halfCircle++;
863 int maxWireGap(0);
864 int turnPoint(0);
865 int deltaWire(0);
867 deltaWire = fabs((*iter)->getWire()-lastLayerHit->getWire());
868 if(deltaWire>(maxCell[(*iter)->getLayer()])/2) deltaWire = maxCell[(*iter)->getLayer()] - deltaWire;
869 if(deltaWire>maxWireGap){
870 maxWireGap = deltaWire;
871 turnPoint = iter - hitOnLayer.begin();
872 }
873 lastLayerHit = *iter;
874 }
875 deltaWire = fabs((*hitIt)->getWire()-lastLayerHit->getWire());
876 if(deltaWire>(maxCell[(*hitIt)->getLayer()])/2) deltaWire = maxCell[(*hitIt)->getLayer()] - deltaWire;
877 if(deltaWire>maxWireGap){
878 maxWireGap = deltaWire;
879 turnPoint = hitOnLayer.size();
880 }
881
882 if(maxWireGap>2){
885 if(shift<turnPoint)(*iter)->setHalfCircle(halfCircle-1);
886 else (*iter)->setHalfCircle(halfCircle);
887 }
888 }else{
891 if(shift<hitOnLayer.size()/2)(*iter)->setHalfCircle(halfCircle-1);
892 else (*iter)->setHalfCircle(halfCircle);
893 }
894 }
895 }
896 hitOnLayer.clear();
897 hitOnLayer.push_back(*hitIt);
898 lastHitLayer = (*hitIt)->getLayer();
899 deltaLayer = (*hitIt)->getLayer()-lastHitLayer;
900 }
901
902 }
903 nMdcHot++;
904 }
905 }
906 //*/
907 ///*
908 if(classifyHotByTrackCenter){
909 int cgemHalf(0), mdcHalf(0);
910 int cgemSign(0), mdcSign(0);
911 int half(0), sign(0);
912 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
913 //(*hitIt)->print();
915 if(track.judgeHalf(*hitIt) != cgemSign){
916 cgemHalf++;
917 cgemSign = track.judgeHalf(*hitIt);
918 half = cgemHalf;
919 sign = cgemSign;
920 }
921 }
923 if(track.judgeHalf(*hitIt) != mdcSign){
924 mdcHalf++;
925 mdcSign = track.judgeHalf(*hitIt);
926 half = mdcHalf;
927 sign = mdcSign;
928 }
929 }
930 (*hitIt)->setHalfCircle(half);
931 }
932 }
933 if(hotList.size()>0){
935 //track.printHot();
936 //cout<<endl;
937 }
938 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++){
939 HoughHit* hitIt = *hotIt;
941 //cout<<hitIt->getHitID();
942 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
943 //cout<<endl;
944 }
946 }
947 sort(m_mcTrackCol.begin(),m_mcTrackCol.end(),largerPt);
948 //cout<<"N MCtrack: "<<m_mcTrackCol.size()<<endl;
949 for(vector<HoughTrack>::iterator trkIter=m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
950 //trkIter->print();
951 }
952 }
953 return m_mcTrackCol.size();
954}
Definition HoughTrack.h:23
Referenced by execute().
◆ initialize()
| StatusCode HoughFinder::initialize | ( | ) |
Definition at line 108 of file Hough.cxx.
109{
110 MsgStream log(msgSvc(), name());
111 log << MSG::INFO << "HoughFinder::initialize()" << endreq;
112 StatusCode sc;
113
114 // --- RawData
115 IRawDataProviderSvc* irawDataProviderSvc;
116 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
118 if ( sc.isFailure() ){
119 log << MSG::ERROR << name()<<" Could not load RawDataProviderSvc!" << endreq;
120 return StatusCode::FAILURE;
121 }
122
123
124 // --- MDC Geometry
125 IMdcGeomSvc* imdcGeomSvc;
126 sc = service ("MdcGeomSvc", imdcGeomSvc);
128 if ( sc.isFailure() ){
129 log << MSG::ERROR << "Could not load MdcGeoSvc!" << endreq;
130 return StatusCode::FAILURE;
131 }
132
133 // --- MdcCalibFunSvc
134 IMdcCalibFunSvc* imdcCalibSvc;
135 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
137 if ( sc.isFailure() ){
138 log << MSG::ERROR << "Could not load MdcCalibFunSvc!" << endreq;
139 return StatusCode::FAILURE;
140 }
141
142 if(m_cgem){
143 // --- CGEM geometry
144 ICgemGeomSvc* icgemGeomSvc;
145 sc = service ("CgemGeomSvc", icgemGeomSvc);
147 if ( sc.isFailure() ){
148 log << MSG::ERROR << "Could not load CgemGeomSvc!" << endreq;
149 return StatusCode::FAILURE;
150 }
151
152 // --- CgemCalibFunSvc
153 ICgemCalibFunSvc* icgemCalibSvc;
154 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
156 if ( sc.isFailure() ){
157 log << MSG::ERROR << "Could not load CgemCalibFunSvc!" << endreq;
158 return StatusCode::FAILURE;
159 }
160
161 }
162
163 // --- MagneticFieldSvc
164 sc = service ("MagneticFieldSvc",m_pIMF);
165 if( sc.isFailure() ) {
166 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
167 return StatusCode::FAILURE;
168 }
170 //log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
173
175 //read pdt
176 Pdt::readMCppTable(m_pdtFile);
177
178 // --- Mdc Print Service
179 //IMdcPrintSvc* imdcPrintSvc;
180 //sc = service ("MdcPrintSvc", imdcPrintSvc);
181 //m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
182 //if ( sc.isFailure() ){
183 //log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
184 //return StatusCode::FAILURE;
185 //}
186
187 // --- Bes Timer Service
188 //sc = service( "BesTimerSvc", m_timersvc);
189 //if( sc.isFailure() ) {
190 //log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
191 //return StatusCode::FAILURE;
192 //}
193 //m_timer_all = m_timersvc->addItem("Execution");
194 //m_timer_all->propName("nExecution");
195
196 HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
197 HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
198 HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
199 HoughHit::setCgemCalibFunSvc(m_cgemCalibFunSvc);
200 HoughHit::setMdcDetector(m_mdcDetector);
203
204 m_roughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",m_nBinTheta,0.,M_PI, m_nBinRho, -1.*m_rhoRange, m_rhoRange);
205 //m_tanlRange = 0.93/sqrt(1-0.93*0.93);//FIXME
206 m_tanlRange = 3.0;//FIXME
207 m_dzRange = 50;
208 m_roughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",m_nBinTanl,-1.*m_tanlRange,m_tanlRange,m_nBinDz,-1.*m_dzRange,m_dzRange);
209
210 const int nPt=12;
211 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};
212 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};
213 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};
214 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 };
215 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};
216 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 };
217 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};
218 m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
219 m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
220 m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
221 m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
222 m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
223 m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
224
225
227 double momenta[100] = {0};
228 double cuts[2][43][100] = {0};
229 string cutFilePath = getenv("HOUGHTRANSALGROOT");
230 cutFilePath += "/share/cut.txt";
231 ifstream fin(cutFilePath.c_str(), ios::in);
232 //cout<<"cutfiledir:"<<cutFilePath.c_str()<<endl;
233 //ifstream fin("/workfs/bes/huangzhen/workarea/CgemBoss/6.6.5.c/Reconstruction/HoughFinder/HoughFinder-00-00-13/share/cut.txt",ios::in);
234 if (!fin.good())
235 cout << "Error in " << __FILE__<<" when opening "<<cutFilePath.c_str()<<endl;
236 bool firstline(true);
237 string strcom;
238 while(std::getline(fin, strcom)){
239 if(firstline){
240 fin >> npt
241 >> momenta[0] >> momenta[1] >> momenta[2] >> momenta[3] >> momenta[4] >> momenta[5] >> momenta[6] >> momenta[7] >> momenta[8] >> momenta[9]
242 >> momenta[10] >> momenta[11] >> momenta[12] >> momenta[13] >> momenta[14] >> momenta[15] >> momenta[16] >> momenta[17] >> momenta[18] >> momenta[19]
243 >> momenta[20] >> momenta[21] >> momenta[22] >> momenta[23] >> momenta[24] >> momenta[25] >> momenta[26] >> momenta[27] >> momenta[28] >> momenta[29]
244 >> momenta[30] >> momenta[31] >> momenta[32] >> momenta[33] >> momenta[34] >> momenta[35] >> momenta[36] >> momenta[37] >> momenta[38] >> momenta[39]
245 >> momenta[40] >> momenta[41] >> momenta[42] >> momenta[43] >> momenta[44] >> momenta[45] >> momenta[46] >> momenta[47] >> momenta[48] >> momenta[49]
246 >> momenta[50] >> momenta[51] >> momenta[52] >> momenta[53] >> momenta[54] >> momenta[55] >> momenta[56] >> momenta[57] >> momenta[58] >> momenta[59]
247 >> momenta[60] >> momenta[61] >> momenta[62] >> momenta[63] >> momenta[64] >> momenta[65] >> momenta[66] >> momenta[67] >> momenta[68] >> momenta[69]
248 >> momenta[70] >> momenta[71] >> momenta[72] >> momenta[73] >> momenta[74] >> momenta[75] >> momenta[76] >> momenta[77] >> momenta[78] >> momenta[79]
249 >> momenta[80] >> momenta[81] >> momenta[82] >> momenta[83] >> momenta[84] >> momenta[85] >> momenta[86] >> momenta[87] >> momenta[88] >> momenta[89]
250 >> momenta[90] >> momenta[91] >> momenta[92] >> momenta[93] >> momenta[94] >> momenta[95] >> momenta[96] >> momenta[97] >> momenta[98] >> momenta[99]
251 ;
252 firstline = false;
253 continue;
254 }
255 fin >> s >> l;
256 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]
257 >> 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]
258 >> 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]
259 >> 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]
260 >> 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]
261 >> 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]
262 >> 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]
263 >> 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]
264 >> 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]
265 >> 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]
266 ;
267 if(fin.peek()<0)break;
268 }
269 for(int is=0;is<=1;is++){
270 for(int il=0;il<43;il++){
272 int point(0);
273 //cout<<il<<" ";
274 for(int ipt=0;ipt<npt;ipt++){
275 //cout<<cuts[is][il][ipt]<<" ";
276 if(fabs(cuts[is][il][ipt])<1e-6)continue;
277 //cut[is][il]->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
278 gr->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
279 point++;
280 }
282 //cout<<endl;
283 }
284 }
285
286
288 if ( sc.isFailure() ){
289 log << MSG::FATAL << "Could not book Tuple !" << endreq;
290 return StatusCode::FAILURE;
291 }
292
293
294 myDotsHelixFitter.initialize();
295
296 return StatusCode::SUCCESS;
297}
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 456 of file Hough.cxx.
457{
458 MsgStream log(msgSvc(), name());
459 int nHit(0);
460 if(!(m_mdcHitCol.empty())) m_mdcHitCol.clear();
461 if(!(m_houghHitList.empty())) m_houghHitList.clear();
462 if(!(m_XHoughHitList.empty())) m_XHoughHitList.clear();
463 if(!(m_VHoughHitList.empty())) m_VHoughHitList.clear();
464
465 if(m_cgem){
466 // --- RecCgemCluster
467 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
468 if(!recCgemClusterCol){
469 log << MSG::ERROR << "Could not retrieve RecCgemClusterCol" << endreq;
470 }else{
471 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterCol->begin();
472 m_recCgemClusterColBegin=cgemClusterIter;
473 for (;cgemClusterIter != recCgemClusterCol->end(); cgemClusterIter++,nHit++){
475 //m_recCgemCluster = (*cgemClusterIter);
476 HoughHit hit(recCgemCluster, m_bunchT0, nHit);
477 m_houghHitList.push_back(hit);
478 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
479 }
480 }
481 }
482
483 uint32_t getDigiFlag = 0;
488
489 MdcDigiVec::iterator mdcDigiIter = mdcDigiVec.begin();
490 for (;mdcDigiIter != mdcDigiVec.end(); mdcDigiIter++,nHit++){
492 HoughHit hit(mdcDigi, m_bunchT0, nHit);
493 if(fabs(hit.driftTime())>m_driftTimeUpLimit) continue;
495 hit.setMdcHit(mdcHit);
496
497 m_mdcHitCol.push_back(mdcHit);
498 m_houghHitList.push_back(hit);
499 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
500 }
501
503 int flag=hitIt->getFlag();
505 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
506 if(flag==2)m_VHoughHitList.push_back(&(*hitIt));
507 }
509 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
510 else m_VHoughHitList.push_back(&(*hitIt));
511 }
512 }
513 //for(vector<HoughHit*>::iterator hitIt = m_XHoughHitList.begin(); hitIt != m_XHoughHitList.end(); hitIt++){
514 //(*hitIt)->print();
515 //}
516 //for(vector<HoughHit*>::iterator hitIt = m_VHoughHitList.begin(); hitIt != m_VHoughHitList.end(); hitIt++){
517 //(*hitIt)->print();
518 //}
519
520 return nHit;
521}
Definition MdcHit.h:44
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition RawDataProviderSvc.h:29
Referenced by execute().
◆ printHitList()
| void HoughFinder::printHitList | ( | ) |
Definition at line 2717 of file Hough.cxx.
2718{
2719 //cout<<"nHit ==========> "<<m_houghHitList.size()<<endl;
2720 cout<<"========== truth =========="<<endl;
2722 hitIt->print();
2723 }
2724 cout<<endl;
2725 cout<<"========== Digi =========="<<endl;
2727 {
2728 hitIt->print();
2729 }
2730 cout<<endl;
2731 /*
2732 {
2733 vector<HoughHit>& hitList = m_mcHitCol;
2734 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2735 cout<<setw(4)<<hitIt->getHitID();
2736 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2737 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2738 vector<int> trkID = hitIt->getTrkID();
2739 cout<<"trkID:"<<setw(4)<<trkID[0];
2740 cout<<"half:"<<setw(4)<<hitIt->getHalfCircle();
2741 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
2742 if(hitIt->getCgemCluster()!=NULL){
2743 //cout<<hitIt->getFlag()<<" ";
2744 //if(hitIt->getFlag()==2)
2745 cout<<setw(5)<<hitIt->getCgemCluster()->getclusterid();
2746 cout<<"["<<setw(3)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(3)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2747 }
2748 //else cout<<"NULL ";
2749 }
2750 else{
2751 if(hitIt->getDigi()!=NULL){
2752 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2753 cout<<hitIt->getDigi()->identify();
2754 }
2755 //else cout<<"NULL ";
2756
2757 }
2758 cout<<endl;
2759 }
2760 }
2761 cout<<endl;
2762 //*/
2763 /*
2764 {
2765 vector<HoughHit>& hitList = m_houghHitList;
2766 cout<<hitList.size()<<endl;
2767 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2768 cout<<setw(4)<<hitIt->getHitID();
2769 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2770 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2771 cout<<setw(4)<<hitIt->getFlag();
2772 if(hitIt->getHitType()==HoughHit::CGEMHIT){
2773 cout<<setw(4)<<hitIt->getCgemCluster()->getclusterid();
2774 cout<<"["<<setw(4)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2775 }
2776 else{
2777 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2778 cout<<hitIt->getDigi()->identify()<<" ";
2779 //cout<<"("<<hitIt->getPairHit()->getLayer()<<", "<<hitIt->getPairHit()->getWire()<<") ";
2780 //if(hitIt->getPairHit()!=NULL)cout<<hitIt->getDriftDist()<<" "<<hitIt->getPairHit()->getDriftDist()<<" "<<hitIt->getDriftDist()-hitIt->getPairHit()->getDriftDist()<<endl;
2781 }
2782 if(hitIt->getPairHit()!=NULL)cout<<hitIt->getPairHit()->getHitID();
2783 //else cout<<"NULL";
2784 cout<<endl;
2785 }
2786 }
2787 cout<<endl;
2788 // */
2789}
Referenced by execute().
◆ printTdsTrack()
| int HoughFinder::printTdsTrack | ( | ) |
Definition at line 4519 of file Hough.cxx.
4520{
4521 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
4522 if(!recMdcTrackCol)return 0;
4523 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
4524 double dr = (*iter)->helix(0) ;
4525 double phi0 = (*iter)->helix(1) ;
4526 double kappa = (*iter)->helix(2) ;
4527 double dz = (*iter)->helix(3) ;
4528 double tanl = (*iter)->helix(4) ;
4529 cout<<setw(12)<<"dr:" <<setw(15)<<dr
4530 <<setw(12)<<"phi0:" <<setw(15)<<phi0
4531 <<setw(12)<<"kappa:" <<setw(15)<<kappa
4532 <<setw(12)<<"dz:" <<setw(15)<<dz
4533 <<setw(12)<<"tanl:" <<setw(15)<<tanl
4534 <<endl;
4535
4536 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
4537 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
4538 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
4539 //cout<<(*clusterIt)->;
4540 }
4541
4542 HitRefVec hitRefVec = (*iter)->getVecHits();
4543 HitRefVec::iterator hotIt = hitRefVec.begin();
4544 //cout<<"("<<"layer"<<","<<"wire"<<") "
4545 cout<<"("<<"l "<<","<<" w "<<") "
4546 //<<"Stat "
4547 //<<"FlagLR "
4548 <<"S "
4549 <<"F "
4550 <<"DriftLeft "
4551 <<" DriftRight "
4552 <<"ErrDrift_L "
4553 <<"ErrDrift_R "
4554 <<"ChisqAdd "
4555 //<<" Tdc "
4556 //<<" Adc "
4557 <<" DriftT "
4558 <<" Doca "
4559 <<" Entra "
4560 <<" Zhit "
4561 <<" FltLen "
4562 <<endl;
4563 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
4566 cout<<"("<<setw( 2)<<layer<<","<<setw( 3)<<wire<<") "
4567 <<setw( 3)<<(*hotIt)->getStat()
4568 <<setw( 3)<<(*hotIt)->getFlagLR()
4569 <<setw(12)<<(*hotIt)->getDriftDistLeft()
4570 <<setw(12)<<(*hotIt)->getDriftDistRight()
4571 <<setw(12)<<(*hotIt)->getErrDriftDistLeft()
4572 <<setw(12)<<(*hotIt)->getErrDriftDistRight()
4573 <<setw(12)<<(*hotIt)->getChisqAdd()
4574 //<<setw( 6)<<(*hotIt)->getTdc()
4575 //<<setw( 6)<<(*hotIt)->getAdc()
4576 <<setw(10)<<(*hotIt)->getDriftT()
4577 <<setw(12)<<(*hotIt)->getDoca()
4578 <<setw(12)<<(*hotIt)->getEntra()
4579 <<setw(10)<<(*hotIt)->getZhit()
4580 <<setw(8 )<<(*hotIt)->getFltLen()
4581 <<endl;
4582 }
4583 cout<<endl;
4584 }
4585 return recMdcTrackCol->size();
4586}
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 1828 of file Hough.cxx.
1829{
1830 MsgStream log(msgSvc(), name());
1831 StatusCode sc;
1832 IDataProviderSvc* eventSvc = NULL;
1833 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1834 if (eventSvc) {
1835 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1836 } else {
1837 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1838 return StatusCode::FAILURE;
1839 //return StatusCode::SUCCESS;
1840 }
1841 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1842
1843
1844 // --- register RecMdcTrack
1846 DataObject *aRecMdcTrackCol;
1847 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1848 if(aRecMdcTrackCol != NULL) {
1849 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1850 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1851 }
1852
1853 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1854 if(sc.isFailure()) {
1855 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1856 return StatusCode::FAILURE;
1857 }
1858 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1859
1860
1861
1863 DataObject *aRecMdcHitCol;
1864 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1865 if(aRecMdcHitCol != NULL) {
1866 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1867 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1868 }
1869
1870 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1871 if(sc.isFailure()) {
1872 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1873 return StatusCode::FAILURE;
1874 }
1875 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1876
1877
1878 return StatusCode::SUCCESS;
1879}//end of stroeTracks
Referenced by execute().
◆ separateHoughHitByPhi()
| int HoughFinder::separateHoughHitByPhi | ( | ) |
◆ solveSharedHits()
| void HoughFinder::solveSharedHits | ( | vector< HoughHit * > & | hitList | ) |
Definition at line 3252 of file Hough.cxx.
3253{
3254 vector<HoughHit*>::iterator itHit = hitList.begin();
3255 for(; itHit!=hitList.end(); itHit++)
3256 {
3257 //vector<HoughTrack*> trkPntList = it->getTrkPntVec();
3258 //int nTrkSharing = trkPntList.size();
3259 vector<int> trkIdList = (*itHit)->getTrkID();
3260 int nTrkSharing = trkIdList.size();
3261 //cout<<"nTrkSharing = "<<nTrkSharing<<endl;
3262 if(nTrkSharing>1) {
3263 //cout<<"found one hit shared by "<<nTrkSharing<<" trks: layer "<<(*itHit)->getLayer()
3264 //<<", pos("<<(*itHit)->getHitPosition().x()<<","<<(*itHit)->getHitPosition().y()<<") with ";
3265 int trkId_minDelD=-1;
3266 double minResid = 99.0;
3267 //vector<int> trkIdList_toDel;
3268 vector<double> vecRes = (*itHit)->getVecResid();
3269 int nTrk = vecRes.size();
3270 //cout<<nTrk<<" residuals: ";
3271 for(int i=0; i<nTrk; i++)// find out the minimum residual
3272 {
3273 //cout<<" trk "<<trkIdList[i]<<":"<<vecRes[i];
3274 if(minResid>fabs(vecRes[i])) {
3275 minResid=fabs(vecRes[i]);
3276 trkId_minDelD=trkIdList[i];
3277 }
3278 }
3279 //cout<<endl;
3280 //cout<<" trkId_minDelD, minResid = "<<trkId_minDelD<<", "<<minResid<<endl;
3281 for(int i=0; i<nTrk; i++)// remove other associations
3282 {
3283 if(trkIdList[i]!=trkId_minDelD) {
3284 (*itHit)->dropTrkID(trkIdList[i]);// remove the trk from the hit
3285 vector<HoughTrack>::iterator itTrk = getHoughTrkIt(m_houghTrackList,trkIdList[i]);
3286 if(itTrk!=m_houghTrackList.end()) itTrk->dropHitPnt(&(*(*itHit)));// remove the hit from the trk
3287 }
3288 }
3289 trkIdList = (*itHit)->getTrkID();
3290 nTrkSharing = trkIdList.size();
3291 //cout<<"new nTrkSharing = "<<nTrkSharing<<endl;
3292 }// if(nTrkSharing>1)
3293 }// loop hits
3294}
vector< HoughTrack >::iterator getHoughTrkIt(vector< HoughTrack > &houghTrkList, int trkId)
Definition Hough.cxx:3297
Referenced by execute().
◆ storeRecTracks()
| int HoughFinder::storeRecTracks | ( | RecMdcTrackCol * | trackList, |
| RecMdcHitCol * | hitList, | ||
| vector< HoughTrack > & | trackVector ) |
Definition at line 1922 of file Hough.cxx.
1923{
1924 bool debug=false; //debug=true;
1925 int trackId(0);
1926 int tkStat =4;
1927 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
1928 sort(trackVector.begin(),trackVector.end(),largerPt);
1929
1930 for(vector<HoughTrack>::iterator trkIT = trackVector.begin(); trkIT!=trackVector.end(); trkIT++)
1931 {
1932
1933 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
1934
1935 //if((trkIT)->getFlag() == 1)continue;
1936 //if((trkIT)->getFlag() == -1)continue;
1937 //if((trkIT)->getFlag() == -2)continue;
1938 if((trkIT)->getFlag()!=0) continue;
1939
1940 if(debug)
1941 {
1942 cout<<"the following track will be saved: "<<endl;
1943 trkIT->print();
1944 }
1945
1946 //new a Rec. Track MdcTrack
1948
1949 recMdcTrack->setTrackId(trackId);
1950
1951 double helixPar[5];
1952 helixPar[0]=trkIT->getDr();
1953 helixPar[1]=trkIT->getPhi0();
1954 helixPar[2]=trkIT->getKappa();
1955 helixPar[3]=trkIT->getDz();
1956 helixPar[4]=trkIT->getTanl();//helixPar[4]+=0.1;
1957 recMdcTrack->setHelix(helixPar);
1958
1960 double pxy = trkIT->pt();
1961 double px = trkIT->momentum(0).x();
1962 double py = trkIT->momentum(0).y();
1963 double pz = trkIT->momentum(0).z();
1964 double p = trkIT->momentum(0).mag();
1965 double theta = trkIT->direction(0).theta();
1966 double phi = trkIT->direction(0).phi();
1967 HepPoint3D poca = trkIT->x(0);
1968 HepPoint3D pivot = trkIT->pivot();
1969 double r = poca.perp();
1970 HepSymMatrix Ea = trkIT->Ea();
1971 //cout<<"Ea="<<Ea<<endl;
1972 double errorMat[15];
1973 int k = 0;
1974 for (int ie = 0 ; ie < 5 ; ie ++){
1975 for (int je = ie ; je < 5 ; je ++){
1976 errorMat[k] = Ea[ie][je];
1977 k++;
1978 }
1979 }
1980 // ---> a Helix for DotsHelixFitter
1981 HepVector a(5,0);
1982 a(1) =trkIT->getDr();
1983 a(2) =trkIT->getPhi0();
1984 a(3) =trkIT->getKappa();
1985 a(4) =trkIT->getDz();
1986 a(5) =trkIT->getTanl();
1988 //myDotsHelixFitter.setInitialHelix(aHelix);// set a Helix for DotsHelixFitter
1989 double chisq = trkIT->getChi2();
1991 recMdcTrack->setPxy(pxy);
1992 recMdcTrack->setPx(px);
1993 recMdcTrack->setPy(py);
1994 recMdcTrack->setPz(pz);
1995 recMdcTrack->setP(p);
1996 recMdcTrack->setTheta(theta);
1997 recMdcTrack->setPhi(phi);
1998 recMdcTrack->setPoca(poca);
2000 recMdcTrack->setY(poca.y());
2001 recMdcTrack->setZ(poca.z());
2002 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
2003 HepPoint3D apivot(0.,0.,0.);
2004 recMdcTrack->setPivot(apivot);
2006 recMdcTrack->setVY0(0.);
2007 recMdcTrack->setVZ0(0.);
2008 recMdcTrack->setError(errorMat);
2009 recMdcTrack->setError(Ea);
2010 recMdcTrack->setChi2(chisq);
2011 recMdcTrack->setStat(tkStat);
2012
2013 //-----------hit list-------------
2014 int hitId = 0;
2015 HitRefVec hitRefVec;
2016 ClusterRefVec clusterRefVec;
2017 //vector<int> vecHits;
2018 map<int,int> clusterFitStat;
2019 //int maxLayer = 0;
2020 int maxLayerId = -1;
2021 int minLayerId = 43;
2022 double fiTerm = 999.;
2023 HoughHit* fiTermHot = NULL;
2024
2025 //vector<HoughHit*> hotList = trkIT->getVecHitPnt();
2026 vector<HoughHit*> hotList = trkIT->getHotList(2);
2027 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2028 {
2029 //cout<<"handle hit "<<hitId<<endl;
2030 HoughHit* hot(*hotIt);
2031 //int layer = (**hotIt).getLayer();
2032 //int wire = (**hotIt).getWire();
2033 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2034
2035 /*
2036 if(hot->getHitType()==HoughHit::MDCHIT)
2037 {
2038 //if(hot->getUse()!=1) continue;
2039 //cout<<"create a new RecMdcHit"<<endl;
2040
2041 // --- make a RecMdcHit from digi, aHelix using DotsHelixFitter
2042 RecMdcHit* recMdcHit = new RecMdcHit(myDotsHelixFitter.makeRecMdcHit(hot->getDigi(), aHelix));
2043 recMdcHit->setId(hitId);
2044 recMdcHit->setTrkId(trackId);
2045
2046 //
2047 // --> make a RecMdcHit from a HoughHit (parameters from Hough)
2048 RecMdcHit* recMdcHit = new RecMdcHit;
2049 recMdcHit->setId(hitId);
2050 recMdcHit->setTrkId(trackId);
2051 HepPoint3D hitPoint = hot->getHitPosition();
2052 double dPhi = trkIT->dPhi(hitPoint);
2053 HepPoint3D position = trkIT->x(dPhi);
2054 Hep3Vector direction = trkIT->direction(dPhi);
2055 //cout<<"get 1"<<endl;
2056
2057 double ddl = -1*hot->getDriftDist();
2058 double ddr = hot->getDriftDist();
2059 //double erddl = hot->hitSigma();
2060 //double erddr = hot->hitSigma();
2061 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2062 int lr = 0;
2063 int stat = hot->getUse();
2064 stat = 1; // llwang
2065 //cout<<"hit stat = "<<stat<<endl;
2066 //cout<<"get 2"<<endl;
2067 Identifier mdcId = hot->getDigi()->identify();
2068 //cout<<"get 2.1"<<endl;
2069 double tdc = hot->getDigi()->getTimeChannel();
2070 //cout<<"get 2.2"<<endl;
2071 double adc = hot->getDigi()->getChargeChannel();
2072 //cout<<"get 2.3"<<endl;
2073 double driftT = hot->driftTime();
2074 //cout<<"get 2.4"<<endl;
2075 double doca = hot->getDriftDist() +hot->getResidual();
2076 //cout<<"get 2.5"<<endl;
2077 double entra = direction.phi()-position.phi();//FIXME
2078 //cout<<"get 2.6"<<endl;
2079 while(entra<-0.5*M_PI)entra+= M_PI;
2080 while(entra> 0.5*M_PI)entra -= M_PI;
2081 //cout<<"get 3"<<endl;
2082 double zhit = hitPoint.z();
2083 double fltLen = trkIT->flightLength(hitPoint);
2084 //cout<<"after get "<<hitId<<endl;
2085
2086 recMdcHit->setDriftDistLeft(ddl);
2087 recMdcHit->setDriftDistRight(ddr);
2088 //recMdcHit->setErrDriftDistLeft(erddl);
2089 //recMdcHit->setErrDriftDistRight(erddr);
2090 recMdcHit->setChisqAdd(pChisq);
2091 recMdcHit->setFlagLR(lr);
2092 recMdcHit->setStat(stat);
2093 recMdcHit->setMdcId(mdcId);
2094 recMdcHit->setTdc(tdc);
2095 recMdcHit->setAdc(adc);
2096 recMdcHit->setDriftT(driftT);
2097 recMdcHit->setDoca(doca);
2098 recMdcHit->setEntra(entra);
2099 recMdcHit->setZhit(zhit);
2100 recMdcHit->setFltLen(fltLen);
2101 //cout<<"after set "<<hitId<<endl;
2102 // <--- make a RecMdcHit from a HoughHit
2103
2104 hitList->push_back(recMdcHit);
2105 SmartRef<RecMdcHit> refHit(recMdcHit);
2106 hitRefVec.push_back(refHit);
2107
2108 if(hot->getLayer()>maxLayerId){
2109 maxLayerId = hot->getLayer();
2110 fiTermHot = hot;
2111 }
2112 if(hot->getLayer()<minLayerId){
2113 minLayerId = hot->getLayer();
2114 }
2115 //if(hot->getUse()==1){
2116 recMdcHit->setStat(1);
2117 //}else{
2118 // recMdcHit->setStat(0);
2119 //}
2120 //cout<<"finish a DC hit "<<hitId<<endl;
2121 hitId++;
2122 }*/
2124 {
2125 //hot->print();
2126 //cout<<endl;
2129 //int stat = hot->getUse();
2130 int stat = 1;
2131 clusterFitStat[clusterid] = stat;
2132 clusterRefVec.push_back(recCgemCluster);
2133 if(debug) cout<<" cluster "<<clusterid<<" kept! "<<endl;
2134 if(hot->getLayer()>maxLayerId)
2135 {
2136 maxLayerId = hot->getLayer();
2137 fiTermHot = hot;
2138 }
2139 if(hot->getLayer()<minLayerId)
2140 {
2141 minLayerId = hot->getLayer();
2142 }
2143 }
2144 }// loop hotList (vector<HoughHit*>) of a hough track
2145
2146 // --- loop RecMdcHitVec in Hough track
2147 double fltLen = -1;
2148 vector<RecMdcHit>* RecMdcHitVecPnt = trkIT->getRecMdcHitVec();
2149 int nMdcHits=RecMdcHitVecPnt->size();
2150 int nMdcHitsKept=0;
2151 vector<RecMdcHit>::iterator iter_recMdcHit = RecMdcHitVecPnt->begin();
2152 for(; iter_recMdcHit!=RecMdcHitVecPnt->end(); iter_recMdcHit++)
2153 {
2156 {
2157 //cout<<"skip hit with chi2="<<recMdcHit->getChisqAdd()<<endl;
2158 delete recMdcHit;
2159 continue;
2160 }
2161 recMdcHit->setId(hitId);
2162 recMdcHit->setTrkId(trackId);
2163 recMdcHit->setStat(1);
2164 hitList->push_back(recMdcHit);
2165 SmartRef<RecMdcHit> refHit(recMdcHit);
2166 hitRefVec.push_back(refHit);
2167 nMdcHitsKept++;
2168
2172 if(debug) cout<<" hit at layer "<<layer<<" wire "<<wire<<" kept! "<<endl;
2173 if(layer>maxLayerId)
2174 {
2175 maxLayerId = layer;
2176 fltLen=recMdcHit->getFltLen();
2177 //fiTermHot = hot;
2178 }
2179 if(layer<minLayerId)
2180 {
2181 minLayerId = layer;
2182 }
2183 hitId++;
2184 }
2185 if(debug) cout<<"track "<<trackId<<", "<<nMdcHitsKept<<"/"<<nMdcHits<<" hits kept"<<endl;
2186
2187 //fi terminal angle
2188 if (fiTermHot!=NULL){
2189 //fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
2191 fiTerm = trkIT->flightArc(point)/trkIT->radius();
2192 }
2194 recMdcTrack->setFiTerm(fiTerm);
2195
2196 // --- setN* functions called in setVec* functions
2197 recMdcTrack->setVecHits(hitRefVec);
2198 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2199 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2200 trackList->push_back(recMdcTrack);
2201 trackId++;
2202 }// loop HoughTrack list
2203 return trackList->size();
2204
2205}
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 1881 of file Hough.cxx.
1882{
1883 MdcTrackParams mdcTrackParams;
1884 //mdcTrackParams.lPrint=8;
1885 mdcTrackParams.lRemoveInActive=1;
1886 mdcTrackParams.maxGapLength=99;
1887 mdcTrackParams.nHitDeleted=99;
1888 MdcTrackList mdcTrackList(mdcTrackParams);
1889 mdcTrackList.setNoInner(true);
1890 vector<MdcTrack*> mdcTrackVector;
1891 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
1892 if((trkIT)->getFlag() == 1)continue;
1893 if((trkIT)->getFlag() == -1)continue;
1894 if((trkIT)->getFlag() == -2)continue;
1895 //(trkIT)->sortHot();
1897 mdcTrackList.append(mdcTrack);
1898 mdcTrackVector.push_back(mdcTrack);
1899 }
1900 int nDeleted(0);
1901 mdcTrackList.setNoInner(true);
1902 if(m_checkHits)nDeleted = mdcTrackList.arbitrateHits();
1903 if(nDeleted>0)cout<<"nDeleted "<<nDeleted<<endl;
1904
1905 int trkID=0;
1906 for(int l=0;l<mdcTrackList.length();l++){
1907 mdcTrackList[l]->storeTrack(trkID,recMdcTrackCol,recMdcHitCol,4);
1908 trkID++;
1909 }
1910
1911 vector<MdcTrack*>::iterator iter = mdcTrackVector.begin();
1914 }
1915 return trkID;
1916}
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 2207 of file Hough.cxx.
2207 {
2208 bool debug=false; //debug=true;
2209 int trackId(0);
2210 int tkStat =4;
2211 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
2212 sort(trackVector.begin(),trackVector.end(),largerPt);
2213
2214 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2215 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
2216 if((trkIT)->getFlag() == 1)continue;
2217 if((trkIT)->getFlag() == -1)continue;
2218 if((trkIT)->getFlag() == -2)continue;
2219 if(trkIT->getTrkRecoTrk()==NULL)
2220 {
2221 if(debug) cout<<"trk getTrkRecoTrk NULL!"<<endl;
2222 continue;
2223 }
2224 //get the results of the fit to this track
2226 //check the fit worked
2227 if (fitresult == 0)
2228 {
2229 if(debug) cout<<"no fit result!"<<endl;
2230 continue;
2231 }
2232 if(debug)
2233 {
2234 cout<<"the following track will be saved: "<<endl;
2235 trkIT->print();
2236 }
2239 //cout<<endl;
2240
2241 //new a Rec. Track MdcTrack
2244 //some track info
2247 //Fit info
2248 int hitId = 0;
2249 int nHits=0;
2250 int nSt=0;
2251 //int nAct=0; int nSeg=0;
2252 //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
2253 //****************************
2254 //* active flag open this
2255 //****************************
2256 //* if (allHit>0){
2257 //* nHits = aList->nHit();//add inActive hit also
2258 //* }else{
2259 //* nHits = fitresult->nMdc();// = nActive()
2260 //* }
2261 //****************************
2262 //for 5.1.0 use all hits
2263 nHits = aList->nHit();
2264 //****************************
2265 // use active only
2266 // nHits = fitresult->nMdc();// = nActive()
2267 //****************************
2268
2272
2273 //if(nHits-nDof!=5) cout<<"HoughFinder::storeTracks: nHits-nDof!=5, nHits="<<nHits<<" , nDof="<<nDof<<", chi2="<<chisq<<" , evt="<<m_event<<endl;
2274
2275 //track parameters at closest approach to beamline
2276 double fltLenPoca = 0.0;
2278 //std::cout<< __FILE__ << " phi0 " << helix.phi0() << " omega " <<helix.omega()<<std::endl;
2281 //double z0 = helix.z0();
2283 //momenta and position
2284 //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
2286
2287 double helixPar[5];
2288 //dro =-d0
2289 helixPar[0] = -d0;
2290 //phi0 = phi0 - pi/2 [0,2pi)
2293 helixPar[1] = tphi0;
2294 //kappa = q/pxy
2296 if (pxy == 0.) helixPar[2] = 9999.;
2298 if(pxy>9999.) helixPar[2] = 0.00001;
2299 //dz = z0
2300 helixPar[3] = helix.z0();
2301 //tanl
2302 helixPar[4] = tanDip;
2303 //error
2304 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
2305 HepSymMatrix mS(helix.params().num_row(),0);
2306 mS[0][0]=-1.;//dr0=-d0
2307 mS[1][1]=1.;
2308 mS[2][2]=-333.567/Bz;//pxy -> cpa
2309 mS[3][3]=1.;
2310 mS[4][4]=1.;
2311 HepSymMatrix mVy = helix.covariance().similarity(mS);
2312 double errorMat[15];
2313 int k = 0;
2314 for (int ie = 0 ; ie < 5 ; ie ++){
2315 for (int je = ie ; je < 5 ; je ++) {
2316 errorMat[k] = mVy[ie][je];
2317 k++;
2318 }
2319 }
2320 double p,px,py,pz;
2321 px = pxy * (-sin(helixPar[1]));
2322 py = pxy * cos(helixPar[1]);
2323 pz = pxy * helixPar[4];
2324 p = sqrt(pxy*pxy + pz*pz);
2325 //theta, phi
2326 double theta = acos(pz/p);
2327 double phi = atan2(py,px);
2328 recMdcTrack->setTrackId(trackId);
2329 recMdcTrack->setHelix(helixPar);
2331 recMdcTrack->setPxy(pxy);
2332 recMdcTrack->setPx(px);
2333 recMdcTrack->setPy(py);
2334 recMdcTrack->setPz(pz);
2335 recMdcTrack->setP(p);
2336 recMdcTrack->setTheta(theta);
2337 recMdcTrack->setPhi(phi);
2338 recMdcTrack->setPoca(poca);
2340 recMdcTrack->setY(poca.y());
2341 recMdcTrack->setZ(poca.z());
2342 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
2343 HepPoint3D pivot(0.,0.,0.);
2344 recMdcTrack->setPivot(pivot);
2346 recMdcTrack->setVY0(0.);
2347 recMdcTrack->setVZ0(0.);
2348 recMdcTrack->setError(errorMat);
2349 recMdcTrack->setError(mVy);
2350 recMdcTrack->setChi2(chisq);
2351 recMdcTrack->setStat(tkStat);
2352 //recMdcTrack->setNhits(nHits);
2353 //-----------hit list-------------
2354 HitRefVec hitRefVec;
2355 ClusterRefVec clusterRefVec;
2356 vector<int> vecHits;
2357 map<int,int> clusterFitStat;
2358 //for fiterm
2359 int maxLayer = 0;
2360 //cout<<__FILE__ <<" "<<__LINE__ <<endl;
2362 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2367 if(maxLayer < layerId)maxLayer = layerId;
2368 }
2369 }
2370 //cout<<"maxLayer:"<<maxLayer<<endl;
2371 int maxLayerId = -1;
2372 int minLayerId = 43;
2373 double fiTerm = 999.;
2377 //cout<<__FILE__ <<" "<<__LINE__ <<" "<< typeid(*hot).name()<<endl;
2378 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2383 if(layerId>(maxLayer - m_removeNOuterHits))continue;
2384 //if (!recoHot->mdcHit()) return;
2386 //id
2387 recMdcHit->setId(hitId);
2388 //---------BESIII----------
2389 //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
2390 //phiWire - phiHit >0 flagLR=1 right driftright>0;
2391 //flag = 2 ambig;
2392 //-----Babar wireAmbig----
2393 //-1 -> right, 0 -> don't know, +1 -> left
2394 //+1 phiWire-phiHit<0 Left,
2395 //-1 phiWire-phiHit>0 Right,
2396 //0 don't know
2397 //ambig() w.r.t track direction
2398 //wireAmbig() w.r.t. wire location
2403 int flagLR=2;
2404 if ( hotWireAmbig == 1){
2405 flagLR = 0; //left driftDist <0
2406 doca *= -1.;//2012-07-19
2407 }else if( hotWireAmbig == -1){
2408 flagLR = 1; //right driftDist >0
2409 }else if( hotWireAmbig == 0){
2410 flagLR = 2; //don't know
2411 }
2412 recMdcHit->setFlagLR(flagLR);
2413 recMdcHit->setDriftDistLeft((-1)*driftDist);
2414 recMdcHit->setErrDriftDistLeft(sigma);
2415 recMdcHit->setDriftDistRight(driftDist);
2416 recMdcHit->setErrDriftDistRight(sigma);
2417 //Mdc Id
2419 recMdcHit->setMdcId(mdcId);
2420 //corrected ADC
2421
2422 //contribution to chisquare
2423 //fill chisq
2424 double res=999.,rese=999.;
2426 }else{}
2427 double deltaChi=0;
2429 recMdcHit->setChisqAdd( deltaChi * deltaChi );
2430 //set tracks containing this hit
2431 recMdcHit->setTrkId(trackId);
2432 //doca
2434 //entrance angle
2435
2437 //z of hit
2438 HepPoint3D pos; Hep3Vector dir;
2440 recMdcHit->setFltLen(fltLen);
2442 recMdcHit->setZhit(pos.z());
2446 //drift time
2448 //for fiterm
2450 //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig "<<hotWireAmbig<<" y "<<pos.y()<<std::endl;
2451 // <<recoHot->entranceAngle()<<std::endl;
2452 if (layerId >= maxLayerId){
2453 maxLayerId = layerId;
2454 fiTermHot = recoHot;
2455 }
2456 if (layerId < minLayerId){
2457 minLayerId = layerId;
2458 }
2459 // status flag
2461 recMdcHit->setStat(1);
2462 }else{
2463 recMdcHit->setStat(0);
2464 }
2465 // for 5.1.0 use all hits
2467 ++nSt;
2468 }
2469 hitList->push_back(recMdcHit);
2470 SmartRef<RecMdcHit> refHit(recMdcHit);
2471 hitRefVec.push_back(refHit);
2472 vecHits.push_back(mdcId.get_value());
2473 ++hitId;
2479 if(stat==0) continue;// added 2020-05-26 wangll
2480 //cout<<"stat:"<<stat<<endl;
2481 clusterFitStat[clusterid] = stat;
2482 //cout<<clusterid<<" "<<stat<<endl;
2483 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2484 clusterRefVec.push_back(recCgemCluster);
2485 }
2486 }
2487 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2488 //fi terminal angle
2489 if (fiTermHot!=NULL){
2491 }
2492 recMdcTrack->setFiTerm(fiTerm);
2493 // number of stereo hits contained
2494 //recMdcTrack->setNhits(hitId);
2495 //recMdcTrack->setNster(nSt);
2496 //recMdcTrack->setNdof(nDof);
2497 //recMdcTrack->setFirstLayer(maxLayerId);
2498 //recMdcTrack->setLastLayer(minLayerId);
2499 //recMdcTrack->setFirstLayer(minLayerId);
2500 //recMdcTrack->setLastLayer(maxLayerId);
2501 //recMdcTrack->setClusterFitStat(clusterFitStat);
2502 //recMdcTrack->setVecClusters(clusterRefVec);
2503 //recMdcTrack->setNcluster(clusterRefVec.size());
2504
2505 // --- setN* functions called in setVec* functions
2506 recMdcTrack->setVecHits(hitRefVec);
2507 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2508 trackList->push_back(recMdcTrack);
2509 trackId++;
2510 }
2511 //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
2512 /*
2513 int trackId(0);
2514 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2515//TrkErrCode trkErrCode = trkIT->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
2516HoughTrack* houghTrack = &(*trkIT);
2517if(houghTrack->getFlag() == 1)continue;
2518if(houghTrack->getFlag() < 0)continue;
2519//houghTrack->print();
2520//houghTrack->printHot();
2521trackId=trackList->size();
2522//cout<<"start store trk "<<trackId;
2523RecMdcTrack* recMdcTrack = new RecMdcTrack();
2524//vector<HoughHit>* hotList = houghTrack->getHot();
2525//vector<HoughHit*> hotList = houghTrack->getVecHitPnt();
2526vector<HoughHit*> hotList = houghTrack->getHotList();
2527//cout<<" with "<<hotList.size()<<" hits";//<<endl;
2528double helixPar[5];
2529helixPar[0] = houghTrack->dr();
2530helixPar[1] = houghTrack->phi0();
2531helixPar[2] = houghTrack->kappa();
2532helixPar[3] = houghTrack->dz();
2533helixPar[4] = houghTrack->tanl();
2534int q = houghTrack->getCharge();
2535double pxy = houghTrack->pt();
2536double px = houghTrack->momentum(0).x();
2537double py = houghTrack->momentum(0).y();
2538double pz = houghTrack->momentum(0).z();
2539double p = houghTrack->momentum(0).mag();
2540double theta = houghTrack->direction(0).theta();
2541double phi = houghTrack->direction(0).phi();
2542HepPoint3D poca = houghTrack->x(0);
2543HepPoint3D pivot = houghTrack->pivot();
2544double r = poca.perp();
2545HepSymMatrix Ea = houghTrack->Ea();
2546//cout<<" phi, pxy, pz = "<<phi<<", "<<pxy<<", "<<pz<<endl;
2547double errorMat[15];
2548int k = 0;
2549for (int ie = 0 ; ie < 5 ; ie ++){
2550for (int je = ie ; je < 5 ; je ++){
2551errorMat[k] = Ea[ie][je];
2552k++;
2553}
2554}
2555double chisq = houghTrack->getChi2();
2556int tkStat = 4+houghTrack->getCircleFitStat();
2557//if(tkStat==4) cout<<" track "<<trackId<<" has no good circle fit"<<endl;
2558int nHits = hotList.size();
2559
2560recMdcTrack->setTrackId(trackId);
2561recMdcTrack->setHelix(helixPar);
2562recMdcTrack->setCharge(q);
2563recMdcTrack->setPxy(pxy);
2564recMdcTrack->setPx(px);
2565recMdcTrack->setPy(py);
2566recMdcTrack->setPz(pz);
2567recMdcTrack->setP(p);
2568recMdcTrack->setTheta(theta);
2569recMdcTrack->setPhi(phi);
2570recMdcTrack->setPoca(poca);
2571recMdcTrack->setX(poca.x());//poca
2572recMdcTrack->setY(poca.y());
2573recMdcTrack->setZ(poca.z());
2574recMdcTrack->setR(r);
2575recMdcTrack->setPivot(pivot);
2576recMdcTrack->setVX0(0.);//pivot
2577recMdcTrack->setVY0(0.);
2578recMdcTrack->setVZ0(0.);
2579recMdcTrack->setError(errorMat);
2580recMdcTrack->setError(Ea);
2581recMdcTrack->setChi2(chisq);
2582recMdcTrack->setStat(tkStat);
2583recMdcTrack->setNhits(nHits);
2584//-----------hit list-------------
2585int hitId = 0;
2586int nSt(0);
2587int maxLayerId = -1;
2588int minLayerId = 43;
2589double fiTerm = 999.;
2590HoughHit fiTermHot;
2591
2592HitRefVec hitRefVec;
2593ClusterRefVec clusterRefVec;
2594map<int,int> clusterFitStat;
2595for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2596{
2597 //cout<<"handle hit "<<hitId<<endl;
2598 HoughHit* hot(*hotIt);
2599 //int layer = (**hotIt).getLayer();
2600 //int wire = (**hotIt).getWire();
2601 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2602 if(hot->getHitType()==HoughHit::MDCHIT){
2603 //cout<<"create a new RecMdcHit"<<endl;
2604 RecMdcHit* recMdcHit = new RecMdcHit;
2605
2606 recMdcHit->setId(hitId);
2607 recMdcHit->setTrkId(trackId);
2608 HepPoint3D hitPoint = hot->getHitPosition();
2609 double dPhi = houghTrack->dPhi(hitPoint);
2610 HepPoint3D position = houghTrack->x(dPhi);
2611 Hep3Vector direction = houghTrack->direction(dPhi);
2612 //cout<<"get 1"<<endl;
2613
2614 double ddl = -1*hot->getDriftDist();
2615 double ddr = hot->getDriftDist();
2616 //double erddl = hot->hitSigma();
2617 //double erddr = hot->hitSigma();
2618 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2619 int lr = 0;
2620 int stat = hot->getUse();
2621 stat = 1; // llwang
2622 //cout<<"hit stat = "<<stat<<endl;
2623 //cout<<"get 2"<<endl;
2624 Identifier mdcId = hot->getDigi()->identify();
2625 //cout<<"get 2.1"<<endl;
2626 double tdc = hot->getDigi()->getTimeChannel();
2627 //cout<<"get 2.2"<<endl;
2628 double adc = hot->getDigi()->getChargeChannel();
2629 //cout<<"get 2.3"<<endl;
2630 double driftT = hot->driftTime();
2631 //cout<<"get 2.4"<<endl;
2632 double doca = hot->getDriftDist() +hot->getResidual();
2633 //cout<<"get 2.5"<<endl;
2634 double entra = direction.phi()-position.phi();//FIXME
2635 //cout<<"get 2.6"<<endl;
2636 while(entra<-0.5*M_PI)entra+= M_PI;
2637 while(entra> 0.5*M_PI)entra -= M_PI;
2638 //cout<<"get 3"<<endl;
2639 double zhit = hitPoint.z();
2640 double fltLen = houghTrack->flightLength(hitPoint);
2641 //cout<<"after get "<<hitId<<endl;
2642
2643 recMdcHit->setDriftDistLeft(ddl);
2644 recMdcHit->setDriftDistRight(ddr);
2645 //recMdcHit->setErrDriftDistLeft(erddl);
2646 //recMdcHit->setErrDriftDistRight(erddr);
2647 recMdcHit->setChisqAdd(pChisq);
2648 recMdcHit->setFlagLR(lr);
2649 recMdcHit->setStat(stat);
2650 recMdcHit->setMdcId(mdcId);
2651 recMdcHit->setTdc(tdc);
2652 recMdcHit->setAdc(adc);
2653 recMdcHit->setDriftT(driftT);
2654 recMdcHit->setDoca(doca);
2655 recMdcHit->setEntra(entra);
2656 recMdcHit->setZhit(zhit);
2657 recMdcHit->setFltLen(fltLen);
2658 //cout<<"after set "<<hitId<<endl;
2659
2660 hitList->push_back(recMdcHit);
2661 SmartRef<RecMdcHit> refHit(recMdcHit);
2662 hitRefVec.push_back(refHit);
2663
2664 if(hot->getLayer()>maxLayerId){
2665 maxLayerId = hot->getLayer();
2666 fiTermHot = *hot;
2667 }
2668 if(hot->getLayer()<minLayerId){
2669 minLayerId = hot->getLayer();
2670 }
2671 if(hot->getUse()==1){
2672 recMdcHit->setStat(1);
2673 }else{
2674 recMdcHit->setStat(0);
2675 }
2676 if(hot->getFlag()!=0)++nSt;
2677 //cout<<"finish a DC hit "<<hitId<<endl;
2678 }else if(hot->getHitType()==HoughHit::CGEMHIT){
2679 //hot->print();
2680 //cout<<endl;
2681 const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
2682 int clusterid = hot->getCgemCluster()->getclusterid();
2683 //int stat = hot->getUse();
2684 int stat = 1;
2685 clusterFitStat[clusterid] = stat;
2686 clusterRefVec.push_back(recCgemCluster);
2687 }
2688 hitId++;
2689}
2690
2691//HepPoint3D point = fiTermHot.getHitPosition();
2692//fiTerm = houghTrack->flightArc(point)/houghTrack->radius();
2693//recMdcTrack->setFiTerm(fiTerm);
2694recMdcTrack->setNdof(nHits-5);
2695recMdcTrack->setNster(nSt);
2696recMdcTrack->setFirstLayer(minLayerId);
2697recMdcTrack->setLastLayer(maxLayerId);
2698//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2699recMdcTrack->setVecHits(hitRefVec);
2700//cout<<"before setVecClusters"<<endl;
2701//recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2702//recMdcTrack->setVecClusters(clusterRefVec);
2703recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2704//cout<<"clusterRefVec.size = "<<clusterRefVec.size()<<endl;
2705//cout<<"end setVecClusters"<<endl;
2706recMdcTrack->setNcluster(clusterRefVec.size());
2707//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2708trackList->push_back(recMdcTrack);
2709//trackId++;
2710//cout<<helixPar[2]<<endl;
2711}
2712return trackList->size();
2713// */
2714return trackList->size();
2715}
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 66 of file Hough.h.
Referenced by HoughFinder().
◆ m_clearTrack
| int HoughFinder::m_clearTrack |
Definition at line 68 of file Hough.h.
Referenced by HoughFinder(), initialize(), and storeTrack().
◆ m_useCgemInGlobalFit
| int HoughFinder::m_useCgemInGlobalFit |
Definition at line 67 of file Hough.h.
Referenced by HoughFinder(), and initialize().
◆ printFlag
| int HoughFinder::printFlag |
Definition at line 71 of file Hough.h.
Referenced by execute(), getMcParticleCol(), HoughFinder(), and storeTracks().
◆ storeFlag
| int HoughFinder::storeFlag |
Definition at line 69 of file Hough.h.
Referenced by execute(), and HoughFinder().
The documentation for this class was generated from the following files:
