CGEM BOSS 6.6.5.f
BESIII Offline Software System
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)
 
void printHitList ()
 
int activeUnusedCgemHitsOnly (vector< HoughHit * > &hitPntList)
 
void checkSharedHits (vector< HoughHit * > &hitList)
 
vector< HoughTrack >::iterator getHoughTrkIt (vector< HoughTrack > &houghTrkList, int trkId)
 
void clearMemory ()
 
int printTdsTrack ()
 
 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)
 
void printHitList ()
 
int activeUnusedCgemHitsOnly (vector< HoughHit * > &hitPntList)
 
void checkSharedHits (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() [1/2]

HoughFinder::HoughFinder ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 31 of file HoughTransAlg/HoughTransAlg-00-00-14/src/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("",m_ = 0);
65 declareProperty("pdtFile", m_pdtFile = "pdt.table");
66
67 // tests by llwang
68 declareProperty("nBinTheta", m_nBinTheta=100 );
69 declareProperty("nBinRho", m_nBinRho=50 );
70 declareProperty("rhoRange", m_rhoRange = 0.1);
71 declareProperty("ExtPeak_theta",m_ExtPeak_theta= 4);
72 declareProperty("ExtPeak_rho", m_ExtPeak_rho= 4);
73 declareProperty("shareHitRate", m_shareHitRate = 0.7);
74
75 declareProperty("nBinTanl", m_nBinTanl = 100);
76 declareProperty("nBinDz", m_nBinDz = 100);
77 declareProperty("ExtPeak_tanl", m_ExtPeak_tanl = 1);
78 declareProperty("ExtPeak_dz", m_ExtPeak_dz = 1);
79 declareProperty("filter", m_filter= 0);
80 declareProperty("eventFile", m_evtFile= "EventList.txt");
81 declareProperty("fitFlag", m_fitFlag= 3);
82 declareProperty("storeFlag", storeFlag=1);
83 declareProperty("checkHits", m_checkHits=1);
84 declareProperty("Layer", Layer=20);
85 declareProperty("clearTrack", m_clearTrack=1);
86 declareProperty("useCgemInGlobalFit",m_useCgemInGlobalFit=3);
87 declareProperty("printFlag", printFlag=0);
88 declareProperty("removeNOuterHits", m_removeNOuterHits = 0);
89
90 m_maxFireLayer = -1;
91 m_totEvtProcessed=0;
92}

◆ HoughFinder() [2/2]

HoughFinder::HoughFinder ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Member Function Documentation

◆ activeUnusedCgemHitsOnly() [1/2]

int HoughFinder::activeUnusedCgemHitsOnly ( vector< HoughHit * > &  hitPntList)

Definition at line 2734 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2735{
2736 int nCgemLeft(0);
2737 vector<HoughHit*>::iterator iter = hitPntList.begin();
2738 for(; iter!=hitPntList.end(); iter++)
2739 {
2740 //if((*iter)->getPairHit()==NULL) continue;// FIXME
2741 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
2742 int useOld = (*iter)->getUse();
2743 int iLayer = (*iter)->getLayer();
2744 if(iLayer<3) {
2745 if(useOld==-1||useOld==0) {
2746 (*iter)->setUse(0);
2747 nCgemLeft++;
2748 }
2749 }
2750 else
2751 if(useOld==0) (*iter)->setUse(-1);
2752 }
2753
2754 return nCgemLeft;
2755}
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)

◆ activeUnusedCgemHitsOnly() [2/2]

int HoughFinder::activeUnusedCgemHitsOnly ( vector< HoughHit * > &  hitPntList)

◆ bookTuple() [1/2]

StatusCode HoughFinder::bookTuple ( )

Definition at line 2981 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2981 {
2982 MsgStream log(msgSvc(), name());
2983
2984 NTuplePtr nt1(ntupleSvc(), "mdcHoughFinder/hit");
2985 if ( nt1 ){
2986 ntuple_hit= nt1;
2987 } else {
2988 ntuple_hit= ntupleSvc()->book("mdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
2989 if(ntuple_hit){
2990 ntuple_hit->addItem("hit_run", m_hit_run);
2991 ntuple_hit->addItem("hit_event", m_hit_event);
2992
2993 ntuple_hit->addItem("hit_nhit", m_hit_nhit, 0, 10000);
2994 ntuple_hit->addItem("hit_hitID", m_hit_nhit, m_hit_hitID);
2995 ntuple_hit->addItem("hit_hitType", m_hit_nhit, m_hit_hitType);
2996 ntuple_hit->addItem("hit_layer", m_hit_nhit, m_hit_layer);
2997 ntuple_hit->addItem("hit_wire", m_hit_nhit, m_hit_wire);
2998 ntuple_hit->addItem("hit_flag", m_hit_nhit, m_hit_flag);
2999 ntuple_hit->addItem("hit_halfCircle", m_hit_nhit, m_hit_halfCircle);
3000 ntuple_hit->addItem("hit_x", m_hit_nhit, m_hit_x);
3001 ntuple_hit->addItem("hit_y", m_hit_nhit, m_hit_y);
3002 ntuple_hit->addItem("hit_z", m_hit_nhit, m_hit_z);
3003 ntuple_hit->addItem("hit_drift", m_hit_nhit, m_hit_drift);
3004
3005 ntuple_hit->addItem("mcHit_hitID", m_hit_nhit, m_mcHit_hitID);
3006 ntuple_hit->addItem("mcHit_hitType", m_hit_nhit, m_mcHit_hitType);
3007 ntuple_hit->addItem("mcHit_layer", m_hit_nhit, m_mcHit_layer);
3008 ntuple_hit->addItem("mcHit_wire", m_hit_nhit, m_mcHit_wire);
3009 ntuple_hit->addItem("mcHit_flag", m_hit_nhit, m_mcHit_flag);
3010 ntuple_hit->addItem("mcHit_halfCircle", m_hit_nhit, m_mcHit_halfCircle);
3011 ntuple_hit->addItem("mcHit_x", m_hit_nhit, m_mcHit_x);
3012 ntuple_hit->addItem("mcHit_y", m_hit_nhit, m_mcHit_y);
3013 ntuple_hit->addItem("mcHit_z", m_hit_nhit, m_mcHit_z);
3014 ntuple_hit->addItem("mcHit_drift", m_hit_nhit, m_mcHit_drift);
3015 } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
3016 return StatusCode::FAILURE;
3017 }
3018 }
3019
3020 NTuplePtr nt2(ntupleSvc(), "mdcHoughFinder/track");
3021 if ( nt2 ){
3022 ntuple_track = nt2;
3023 } else {
3024 ntuple_track = ntupleSvc()->book("mdcHoughFinder/track", CLID_ColumnWiseTuple, "track");
3025 if(ntuple_track){
3026 ntuple_track->addItem("trk_run", m_trk_run);
3027 ntuple_track->addItem("trk_event", m_trk_event);
3028 ntuple_track->addItem("trk_nTrack", m_trk_nTrack);
3029 ntuple_track->addItem("trk_trackID", m_trk_trackID);
3030 ntuple_track->addItem("trk_charge", m_trk_charge);
3031 ntuple_track->addItem("trk_flag", m_trk_flag);
3032 ntuple_track->addItem("trk_angle", m_trk_angle);
3033 ntuple_track->addItem("trk_rho", m_trk_rho);
3034 ntuple_track->addItem("trk_dAngle", m_trk_dAngle);
3035 ntuple_track->addItem("trk_dRho", m_trk_dRho);
3036 ntuple_track->addItem("trk_dTanl", m_trk_dTanl);
3037 ntuple_track->addItem("trk_dDz", m_trk_dDz);
3038 ntuple_track->addItem("trk_Xc", m_trk_Xc);
3039 ntuple_track->addItem("trk_Yc", m_trk_Yc);
3040 ntuple_track->addItem("trk_R", m_trk_R);
3041 ntuple_track->addItem("trk_dr", m_trk_dr);
3042 ntuple_track->addItem("trk_phi0", m_trk_phi0);
3043 ntuple_track->addItem("trk_kappa", m_trk_kappa);
3044 ntuple_track->addItem("trk_dz", m_trk_dz);
3045 ntuple_track->addItem("trk_tanl", m_trk_tanl);
3046 ntuple_track->addItem("trk_pxy", m_trk_pxy);
3047 ntuple_track->addItem("trk_px", m_trk_px);
3048 ntuple_track->addItem("trk_py", m_trk_py);
3049 ntuple_track->addItem("trk_pz", m_trk_pz);
3050 ntuple_track->addItem("trk_p", m_trk_p);
3051 ntuple_track->addItem("trk_phi", m_trk_phi);
3052 ntuple_track->addItem("trk_theta", m_trk_theta);
3053 ntuple_track->addItem("trk_cosTheta", m_trk_cosTheta);
3054 ntuple_track->addItem("trk_vx", m_trk_vx);
3055 ntuple_track->addItem("trk_vy", m_trk_vy);
3056 ntuple_track->addItem("trk_vz", m_trk_vz);
3057 ntuple_track->addItem("trk_vr", m_trk_vr);
3058 ntuple_track->addItem("trk_chi2", m_trk_chi2);
3059 ntuple_track->addItem("trk_fiTerm", m_trk_fiTerm);
3060 ntuple_track->addItem("trk_nhit", m_trk_nhit);
3061 ntuple_track->addItem("trk_ncluster", m_trk_ncluster);
3062 ntuple_track->addItem("trk_stat", m_trk_stat);
3063 ntuple_track->addItem("trk_ndof", m_trk_ndof);
3064 ntuple_track->addItem("trk_nster", m_trk_nster);
3065 ntuple_track->addItem("trk_nlayer", m_trk_nlayer);
3066 ntuple_track->addItem("trk_firstLayer", m_trk_firstLayer);
3067 ntuple_track->addItem("trk_lastLayer", m_trk_lastLayer);
3068 ntuple_track->addItem("trk_nCgemXClusters", m_trk_nCgemXClusters);
3069 ntuple_track->addItem("trk_nCgemVClusters", m_trk_nCgemVClusters);
3070
3071 ntuple_track->addItem("trk_nHot", m_trk_nHot, 0, 10000);
3072 ntuple_track->addItem("hot_hitID", m_trk_nHot, m_hot_hitID);
3073 ntuple_track->addItem("hot_hitType", m_trk_nHot, m_hot_hitType);
3074 ntuple_track->addItem("hot_layer", m_trk_nHot, m_hot_layer);
3075 ntuple_track->addItem("hot_wire", m_trk_nHot, m_hot_wire);
3076 ntuple_track->addItem("hot_flag", m_trk_nHot, m_hot_flag);
3077 ntuple_track->addItem("hot_halfCircle", m_trk_nHot, m_hot_halfCircle);
3078 ntuple_track->addItem("hot_x", m_trk_nHot, m_hot_x);
3079 ntuple_track->addItem("hot_y", m_trk_nHot, m_hot_y);
3080 ntuple_track->addItem("hot_z", m_trk_nHot, m_hot_z);
3081 ntuple_track->addItem("hot_drift", m_trk_nHot, m_hot_drift);
3082 ntuple_track->addItem("hot_residual", m_trk_nHot, m_hot_residual);
3083
3084 ntuple_track->addItem("mcHot_hitID", m_trk_nHot, m_mcHot_hitID);
3085 ntuple_track->addItem("mcHot_hitType", m_trk_nHot, m_mcHot_hitType);
3086 ntuple_track->addItem("mcHot_layer", m_trk_nHot, m_mcHot_layer);
3087 ntuple_track->addItem("mcHot_wire", m_trk_nHot, m_mcHot_wire);
3088 ntuple_track->addItem("mcHot_flag", m_trk_nHot, m_mcHot_flag);
3089 ntuple_track->addItem("mcHot_halfCircle", m_trk_nHot, m_mcHot_halfCircle);
3090 ntuple_track->addItem("mcHot_x", m_trk_nHot, m_mcHot_x);
3091 ntuple_track->addItem("mcHot_y", m_trk_nHot, m_mcHot_y);
3092 ntuple_track->addItem("mcHot_z", m_trk_nHot, m_mcHot_z);
3093 ntuple_track->addItem("mcHot_drift", m_trk_nHot, m_mcHot_drift);
3094
3095 ntuple_track->addItem("mcTrk_trackID", m_mcTrk_trackID);
3096 ntuple_track->addItem("mcTrk_charge", m_mcTrk_charge);
3097 ntuple_track->addItem("mcTrk_flag", m_mcTrk_flag);
3098 ntuple_track->addItem("mcTrk_angle", m_mcTrk_angle);
3099 ntuple_track->addItem("mcTrk_rho", m_mcTrk_rho);
3100 ntuple_track->addItem("mcTrk_dAngle", m_mcTrk_dAngle);
3101 ntuple_track->addItem("mcTrk_dRho", m_mcTrk_dRho);
3102 ntuple_track->addItem("mcTrk_dTanl", m_mcTrk_dTanl);
3103 ntuple_track->addItem("mcTrk_dDz", m_mcTrk_dDz);
3104 ntuple_track->addItem("mcTrk_Xc", m_mcTrk_Xc);
3105 ntuple_track->addItem("mcTrk_Yc", m_mcTrk_Yc);
3106 ntuple_track->addItem("mcTrk_R", m_mcTrk_R);
3107 ntuple_track->addItem("mcTrk_dr", m_mcTrk_dr);
3108 ntuple_track->addItem("mcTrk_phi0", m_mcTrk_phi0);
3109 ntuple_track->addItem("mcTrk_kappa", m_mcTrk_kappa);
3110 ntuple_track->addItem("mcTrk_dz", m_mcTrk_dz);
3111 ntuple_track->addItem("mcTrk_tanl", m_mcTrk_tanl);
3112 ntuple_track->addItem("mcTrk_pxy", m_mcTrk_pxy);
3113 ntuple_track->addItem("mcTrk_px", m_mcTrk_px);
3114 ntuple_track->addItem("mcTrk_py", m_mcTrk_py);
3115 ntuple_track->addItem("mcTrk_pz", m_mcTrk_pz);
3116 ntuple_track->addItem("mcTrk_p", m_mcTrk_p);
3117 ntuple_track->addItem("mcTrk_phi", m_mcTrk_phi);
3118 ntuple_track->addItem("mcTrk_theta", m_mcTrk_theta);
3119 ntuple_track->addItem("mcTrk_cosTheta", m_mcTrk_cosTheta);
3120 ntuple_track->addItem("mcTrk_vx", m_mcTrk_vx);
3121 ntuple_track->addItem("mcTrk_vy", m_mcTrk_vy);
3122 ntuple_track->addItem("mcTrk_vz", m_mcTrk_vz);
3123 ntuple_track->addItem("mcTrk_vr", m_mcTrk_vr);
3124 ntuple_track->addItem("mcTrk_chi2", m_mcTrk_chi2);
3125 ntuple_track->addItem("mcTrk_fiTerm", m_mcTrk_fiTerm);
3126 ntuple_track->addItem("mcTrk_nhit", m_mcTrk_nhit);
3127 ntuple_track->addItem("mcTrk_ncluster", m_mcTrk_ncluster);
3128 ntuple_track->addItem("mcTrk_stat", m_mcTrk_stat);
3129 ntuple_track->addItem("mcTrk_ndof", m_mcTrk_ndof);
3130 ntuple_track->addItem("mcTrk_nster", m_mcTrk_nster);
3131 ntuple_track->addItem("mcTrk_nlayer", m_mcTrk_nlayer);
3132 ntuple_track->addItem("mcTrk_firstLayer", m_mcTrk_firstLayer);
3133 ntuple_track->addItem("mcTrk_lastLayer", m_mcTrk_lastLayer);
3134 ntuple_track->addItem("mcTrk_nCgemXClusters", m_mcTrk_nCgemXClusters);
3135 ntuple_track->addItem("mcTrk_nCgemVClusters", m_mcTrk_nCgemVClusters);
3136
3137 ntuple_track->addItem("mcTrk_nHot", m_mcTrk_nHot, 0, 10000);
3138 ntuple_track->addItem("mcTrkHot_hitID", m_mcTrk_nHot, m_mcTrkHot_hitID);
3139 ntuple_track->addItem("mcTrkHot_hitType", m_mcTrk_nHot, m_mcTrkHot_hitType);
3140 ntuple_track->addItem("mcTrkHot_layer", m_mcTrk_nHot, m_mcTrkHot_layer);
3141 ntuple_track->addItem("mcTrkHot_wire", m_mcTrk_nHot, m_mcTrkHot_wire);
3142 ntuple_track->addItem("mcTrkHot_flag", m_mcTrk_nHot, m_mcTrkHot_flag);
3143 ntuple_track->addItem("mcTrkHot_halfCircle", m_mcTrk_nHot, m_mcTrkHot_halfCircle);
3144 ntuple_track->addItem("mcTrkHot_x", m_mcTrk_nHot, m_mcTrkHot_x);
3145 ntuple_track->addItem("mcTrkHot_y", m_mcTrk_nHot, m_mcTrkHot_y);
3146 ntuple_track->addItem("mcTrkHot_z", m_mcTrk_nHot, m_mcTrkHot_z);
3147 ntuple_track->addItem("mcTrkHot_drift", m_mcTrk_nHot, m_mcTrkHot_drift);
3148 } else {
3149 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/track" <<endmsg;
3150 return StatusCode::FAILURE;
3151 }
3152 }
3153
3154 NTuplePtr nt3(ntupleSvc(), "mdcHoughFinder/event");
3155 if ( nt3 ){
3156 ntuple_event = nt3;
3157 } else {
3158 ntuple_event = ntupleSvc()->book("mdcHoughFinder/event", CLID_ColumnWiseTuple, "event");
3159 if(ntuple_event){
3160 ntuple_event->addItem("evt_run", m_evt_run);
3161 ntuple_event->addItem("evt_event", m_evt_event);
3162 ntuple_event->addItem("evt_nXCluster", m_evt_nXCluster);
3163 ntuple_event->addItem("evt_nVCluster", m_evt_nVCluster);
3164 ntuple_event->addItem("evt_nXVCluster", m_evt_nXVCluster);
3165 ntuple_event->addItem("evt_nCgemCluster", m_evt_nCGEMCluster);
3166 ntuple_event->addItem("evt_nAxialHit", m_evt_nAxialHit);
3167 ntuple_event->addItem("evt_nStereoHit", m_evt_nStereoHit);
3168 ntuple_event->addItem("evt_nODCHit", m_evt_nODCHit);
3169 ntuple_event->addItem("evt_nHit", m_evt_nHit);
3170
3171 ntuple_event->addItem("evt_nTrack", m_evt_nTrack, 0, 100);
3172 ntuple_event->addItem("evtTrk_trackID", m_evt_nTrack, m_evtTrk_trackID);
3173 ntuple_event->addItem("evtTrk_charge", m_evt_nTrack, m_evtTrk_charge);
3174 ntuple_event->addItem("evtTrk_flag", m_evt_nTrack, m_evtTrk_flag);
3175 ntuple_event->addItem("evtTrk_angle", m_evt_nTrack, m_evtTrk_angle);
3176 ntuple_event->addItem("evtTrk_rho", m_evt_nTrack, m_evtTrk_rho);
3177 ntuple_event->addItem("evtTrk_dAngle", m_evt_nTrack, m_evtTrk_dAngle);
3178 ntuple_event->addItem("evtTrk_dRho", m_evt_nTrack, m_evtTrk_dRho);
3179 ntuple_event->addItem("evtTrk_dTanl", m_evt_nTrack, m_evtTrk_dTanl);
3180 ntuple_event->addItem("evtTrk_dDz", m_evt_nTrack, m_evtTrk_dDz);
3181 ntuple_event->addItem("evtTrk_Xc", m_evt_nTrack, m_evtTrk_Xc);
3182 ntuple_event->addItem("evtTrk_Yc", m_evt_nTrack, m_evtTrk_Yc);
3183 ntuple_event->addItem("evtTrk_R", m_evt_nTrack, m_evtTrk_R);
3184 ntuple_event->addItem("evtTrk_dr", m_evt_nTrack, m_evtTrk_dr);
3185 ntuple_event->addItem("evtTrk_phi0", m_evt_nTrack, m_evtTrk_phi0);
3186 ntuple_event->addItem("evtTrk_kappa", m_evt_nTrack, m_evtTrk_kappa);
3187 ntuple_event->addItem("evtTrk_dz", m_evt_nTrack, m_evtTrk_dz);
3188 ntuple_event->addItem("evtTrk_tanl", m_evt_nTrack, m_evtTrk_tanl);
3189 ntuple_event->addItem("evtTrk_pxy", m_evt_nTrack, m_evtTrk_pxy);
3190 ntuple_event->addItem("evtTrk_px", m_evt_nTrack, m_evtTrk_px);
3191 ntuple_event->addItem("evtTrk_py", m_evt_nTrack, m_evtTrk_py);
3192 ntuple_event->addItem("evtTrk_pz", m_evt_nTrack, m_evtTrk_pz);
3193 ntuple_event->addItem("evtTrk_p", m_evt_nTrack, m_evtTrk_p);
3194 ntuple_event->addItem("evtTrk_phi", m_evt_nTrack, m_evtTrk_phi);
3195 ntuple_event->addItem("evtTrk_theta", m_evt_nTrack, m_evtTrk_theta);
3196 ntuple_event->addItem("evtTrk_cosTheta", m_evt_nTrack, m_evtTrk_cosTheta);
3197 ntuple_event->addItem("evtTrk_vx", m_evt_nTrack, m_evtTrk_vx);
3198 ntuple_event->addItem("evtTrk_vy", m_evt_nTrack, m_evtTrk_vy);
3199 ntuple_event->addItem("evtTrk_vz", m_evt_nTrack, m_evtTrk_vz);
3200 ntuple_event->addItem("evtTrk_vr", m_evt_nTrack, m_evtTrk_vr);
3201 ntuple_event->addItem("evtTrk_chi2", m_evt_nTrack, m_evtTrk_chi2);
3202 ntuple_event->addItem("evtTrk_fiTerm", m_evt_nTrack, m_evtTrk_fiTerm);
3203 ntuple_event->addItem("evtTrk_nhit", m_evt_nTrack, m_evtTrk_nhit);
3204 ntuple_event->addItem("evtTrk_ncluster", m_evt_nTrack, m_evtTrk_ncluster);
3205 ntuple_event->addItem("evtTrk_stat", m_evt_nTrack, m_evtTrk_stat);
3206 ntuple_event->addItem("evtTrk_ndof", m_evt_nTrack, m_evtTrk_ndof);
3207 ntuple_event->addItem("evtTrk_nster", m_evt_nTrack, m_evtTrk_nster);
3208 ntuple_event->addItem("evtTrk_nlayer", m_evt_nTrack, m_evtTrk_nlayer);
3209 ntuple_event->addItem("evtTrk_firstLayer", m_evt_nTrack, m_evtTrk_firstLayer);
3210 ntuple_event->addItem("evtTrk_lastLayer", m_evt_nTrack, m_evtTrk_lastLayer);
3211 ntuple_event->addItem("evtTrk_nCgemXClusters", m_evt_nTrack, m_evtTrk_nCgemXClusters);
3212 ntuple_event->addItem("evtTrk_nCgemVClusters", m_evt_nTrack, m_evtTrk_nCgemVClusters);
3213
3214 ntuple_event->addItem("mcEvtTrk_trackID", m_evt_nTrack, m_mcEvtTrk_trackID);
3215 ntuple_event->addItem("mcEvtTrk_charge", m_evt_nTrack, m_mcEvtTrk_charge);
3216 ntuple_event->addItem("mcEvtTrk_flag", m_evt_nTrack, m_mcEvtTrk_flag);
3217 ntuple_event->addItem("mcEvtTrk_angle", m_evt_nTrack, m_mcEvtTrk_angle);
3218 ntuple_event->addItem("mcEvtTrk_rho", m_evt_nTrack, m_mcEvtTrk_rho);
3219 ntuple_event->addItem("mcEvtTrk_dAngle", m_evt_nTrack, m_mcEvtTrk_dAngle);
3220 ntuple_event->addItem("mcEvtTrk_dRho", m_evt_nTrack, m_mcEvtTrk_dRho);
3221 ntuple_event->addItem("mcEvtTrk_dTanl", m_evt_nTrack, m_mcEvtTrk_dTanl);
3222 ntuple_event->addItem("mcEvtTrk_dDz", m_evt_nTrack, m_mcEvtTrk_dDz);
3223 ntuple_event->addItem("mcEvtTrk_Xc", m_evt_nTrack, m_mcEvtTrk_Xc);
3224 ntuple_event->addItem("mcEvtTrk_Yc", m_evt_nTrack, m_mcEvtTrk_Yc);
3225 ntuple_event->addItem("mcEvtTrk_R", m_evt_nTrack, m_mcEvtTrk_R);
3226 ntuple_event->addItem("mcEvtTrk_dr", m_evt_nTrack, m_mcEvtTrk_dr);
3227 ntuple_event->addItem("mcEvtTrk_phi0", m_evt_nTrack, m_mcEvtTrk_phi0);
3228 ntuple_event->addItem("mcEvtTrk_kappa", m_evt_nTrack, m_mcEvtTrk_kappa);
3229 ntuple_event->addItem("mcEvtTrk_dz", m_evt_nTrack, m_mcEvtTrk_dz);
3230 ntuple_event->addItem("mcEvtTrk_tanl", m_evt_nTrack, m_mcEvtTrk_tanl);
3231 ntuple_event->addItem("mcEvtTrk_pxy", m_evt_nTrack, m_mcEvtTrk_pxy);
3232 ntuple_event->addItem("mcEvtTrk_px", m_evt_nTrack, m_mcEvtTrk_px);
3233 ntuple_event->addItem("mcEvtTrk_py", m_evt_nTrack, m_mcEvtTrk_py);
3234 ntuple_event->addItem("mcEvtTrk_pz", m_evt_nTrack, m_mcEvtTrk_pz);
3235 ntuple_event->addItem("mcEvtTrk_p", m_evt_nTrack, m_mcEvtTrk_p);
3236 ntuple_event->addItem("mcEvtTrk_phi", m_evt_nTrack, m_mcEvtTrk_phi);
3237 ntuple_event->addItem("mcEvtTrk_theta", m_evt_nTrack, m_mcEvtTrk_theta);
3238 ntuple_event->addItem("mcEvtTrk_cosTheta", m_evt_nTrack, m_mcEvtTrk_cosTheta);
3239 ntuple_event->addItem("mcEvtTrk_vx", m_evt_nTrack, m_mcEvtTrk_vx);
3240 ntuple_event->addItem("mcEvtTrk_vy", m_evt_nTrack, m_mcEvtTrk_vy);
3241 ntuple_event->addItem("mcEvtTrk_vz", m_evt_nTrack, m_mcEvtTrk_vz);
3242 ntuple_event->addItem("mcEvtTrk_vr", m_evt_nTrack, m_mcEvtTrk_vr);
3243 ntuple_event->addItem("mcEvtTrk_chi2", m_evt_nTrack, m_mcEvtTrk_chi2);
3244 ntuple_event->addItem("mcEvtTrk_fiTerm", m_evt_nTrack, m_mcEvtTrk_fiTerm);
3245 ntuple_event->addItem("mcEvtTrk_nhit", m_evt_nTrack, m_mcEvtTrk_nhit);
3246 ntuple_event->addItem("mcEvtTrk_ncluster", m_evt_nTrack, m_mcEvtTrk_ncluster);
3247 ntuple_event->addItem("mcEvtTrk_stat", m_evt_nTrack, m_mcEvtTrk_stat);
3248 ntuple_event->addItem("mcEvtTrk_ndof", m_evt_nTrack, m_mcEvtTrk_ndof);
3249 ntuple_event->addItem("mcEvtTrk_nster", m_evt_nTrack, m_mcEvtTrk_nster);
3250 ntuple_event->addItem("mcEvtTrk_nlayer", m_evt_nTrack, m_mcEvtTrk_nlayer);
3251 ntuple_event->addItem("mcEvtTrk_firstLayer", m_evt_nTrack, m_mcEvtTrk_firstLayer);
3252 ntuple_event->addItem("mcEvtTrk_lastLayer", m_evt_nTrack, m_mcEvtTrk_lastLayer);
3253 ntuple_event->addItem("mcEvtTrk_nCgemXClusters", m_evt_nTrack, m_mcEvtTrk_nCgemXClusters);
3254 ntuple_event->addItem("mcEvtTrk_nCgemVClusters", m_evt_nTrack, m_mcEvtTrk_nCgemVClusters);
3255 } else {
3256 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/event" <<endmsg;
3257 return StatusCode::FAILURE;
3258 }
3259 }
3260 return StatusCode::SUCCESS;
3261}

Referenced by initialize().

◆ bookTuple() [2/2]

StatusCode HoughFinder::bookTuple ( )

◆ checkHot() [1/2]

int HoughFinder::checkHot ( vector< HoughTrack > &  trackVector)

Definition at line 1598 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1599{
1600 int nDelete(0);
1601 if(trackVector.size()==0)return 0;
1602 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1603
1604 vector<HoughTrack>::iterator trkIT1 = trackVector.end();
1605 //vector<HoughTrack>::iterator trkIT1 = trackVector.rbegin();
1606 bool loop = true;
1607
1608 while(loop)
1609 {
1610 trkIT1--;
1611 if(trkIT1==trackVector.begin()) loop=false;
1612 if((trkIT1)->getFlag() == 1) continue;
1613 trkIT1->dropRedundentCgemXHits();
1614 //cout<<"hit size "<<hot1->size()<<endl;
1615 int nHits = trkIT1->getVecHitPnt().size();
1616 if(nHits<3){
1617 (trkIT1)->setFlag(1);
1618 trkIT1->releaseSelHits();
1619 continue;
1620 }
1621 int nShared = trkIT1->getNHitsShared();
1622
1623 //cout<<"shareHit = "<<shareHit<<endl;
1624 if(nShared>m_shareHitRate*nHits){
1625 trkIT1->setFlag(1);
1626 trkIT1->releaseSelHits();
1627 nDelete++;
1628 //cout<<"too many shared hits ("<<nShared<<"/"<<nHits
1629 //<<") in track "<<(trkIT1)->getTrkID()
1630 //<<endl;
1631 }
1632 }
1633 return nDelete;
1634}
bool moreHot(HoughTrack trk1, HoughTrack trk2)

Referenced by execute().

◆ checkHot() [2/2]

int HoughFinder::checkHot ( vector< HoughTrack > &  trackVector)

◆ checkSharedHits() [1/2]

void HoughFinder::checkSharedHits ( vector< HoughHit * > &  hitList)

Definition at line 2758 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2759{
2760 vector<HoughHit*>::iterator itHit = hitList.begin();
2761 for(; itHit!=hitList.end(); itHit++)
2762 {
2763 //vector<HoughTrack*> trkPntList = it->getTrkPntVec();
2764 //int nTrkSharing = trkPntList.size();
2765 vector<int> trkIdList = (*itHit)->getTrkID();
2766 int nTrkSharing = trkIdList.size();
2767 //cout<<"nTrkSharing = "<<nTrkSharing<<endl;
2768 if(nTrkSharing>1) {
2769 //cout<<"found one hit shared by "<<nTrkSharing<<" trks: layer "<<(*itHit)->getLayer()
2770 //<<", pos("<<(*itHit)->getHitPosition().x()<<","<<(*itHit)->getHitPosition().y()<<") with ";
2771 int trkId_minDelD=-1;
2772 double minResid = 99.0;
2773 //vector<int> trkIdList_toDel;
2774 vector<double> vecRes = (*itHit)->getVecResid();
2775 int nTrk = vecRes.size();
2776 //cout<<nTrk<<" residuals: ";
2777 for(int i=0; i<nTrk; i++)// find out the minimum residual
2778 {
2779 //cout<<" trk "<<trkIdList[i]<<":"<<vecRes[i];
2780 if(minResid>fabs(vecRes[i])) {
2781 minResid=fabs(vecRes[i]);
2782 trkId_minDelD=trkIdList[i];
2783 }
2784 }
2785 //cout<<endl;
2786 //cout<<" trkId_minDelD, minResid = "<<trkId_minDelD<<", "<<minResid<<endl;
2787 for(int i=0; i<nTrk; i++)// remove other associations
2788 {
2789 if(trkIdList[i]!=trkId_minDelD) {
2790 (*itHit)->dropTrkID(trkIdList[i]);// remove the trk from the hit
2791 vector<HoughTrack>::iterator itTrk = getHoughTrkIt(m_houghTrackList,trkIdList[i]);
2792 if(itTrk!=m_houghTrackList.end()) itTrk->dropHitPnt(&(*(*itHit)));// remove the hit from the trk
2793 }
2794 }
2795 trkIdList = (*itHit)->getTrkID();
2796 nTrkSharing = trkIdList.size();
2797 //cout<<"new nTrkSharing = "<<nTrkSharing<<endl;
2798 }// if(nTrkSharing>1)
2799 }// loop hits
2800}
vector< HoughTrack >::iterator getHoughTrkIt(vector< HoughTrack > &houghTrkList, int trkId)

Referenced by execute().

◆ checkSharedHits() [2/2]

void HoughFinder::checkSharedHits ( vector< HoughHit * > &  hitList)

◆ checkTrack() [1/2]

int HoughFinder::checkTrack ( vector< HoughTrack > &  trackVector)

Definition at line 1637 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1638{
1639 int nDelete(0);
1640 double dDr(999);
1641 double dPhi0(999);
1642 double dKappa(999);
1643 double dDz(999);
1644 double dTanl(999);
1645
1646 std::sort(trackVector.begin(),trackVector.end(),moreHot);
1647 for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
1648 if((trkIT1)->getFlag() == 1)continue;
1649 int charge = (trkIT1)->getCharge();
1650 double dr = (trkIT1)->dr();
1651 double phi0 = (trkIT1)->phi0();
1652 double kappa = (trkIT1)->kappa();
1653 double dz = (trkIT1)->dz();
1654 double tanl = (trkIT1)->tanl();
1655 for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
1656 if((trkIT2)->getFlag() == 1)continue;
1657 bool sameTrack(false);
1658 dDr = fabs(dr - (trkIT2)->dr());
1659 dPhi0 = fabs(phi0 - (trkIT2)->phi0());
1660 if((trkIT2)->getCharge() == charge){
1661 dKappa = fabs(kappa - (trkIT2)->kappa());
1662 dTanl = fabs(tanl - (trkIT2)->tanl());
1663 }else{
1664 dKappa = fabs(fabs(kappa) - fabs((trkIT2)->kappa()));
1665 dTanl = fabs(fabs(tanl) - fabs((trkIT2)->tanl()));
1666 }
1667 double r = fabs(dz)<fabs((trkIT2)->dz())?(trkIT1)->radius():(trkIT2)->radius();
1668 double k = dz<fabs((trkIT2)->dz())?tanl:(trkIT2)->tanl();
1669 dDz = fabs(dz - (trkIT2)->dz());
1670 int nCircle = fabs(dDz)/(2*M_PI*r*k);
1671 dDz = fabs(dDz - nCircle*(2*M_PI*r*k));
1672 if(dDz>M_PI*r*k)dDz = fabs(dDz - 2*M_PI*r*k);
1673 sameTrack = sameTrack&&(dDr<m_drCut);//FIXME
1674 sameTrack = sameTrack&&(dPhi0<m_phi0Cut);//FIXME
1675 sameTrack = sameTrack&&(dKappa<m_kappaCut);//FIXME
1676 sameTrack = sameTrack&&(dDz<m_dzCut);//FIXME
1677 sameTrack = sameTrack&&(dTanl<m_tanlCut);//FIXME
1678 if(sameTrack){
1679 //delete *trkIT2;
1680 //trackVector.erase(trkIT2);
1681 //trkIT2--;
1682 (trkIT2)->setFlag(1);
1683 nDelete++;
1684 }
1685 }
1686 }
1687 return nDelete;
1688}
#define M_PI
Definition: TConstant.h:4

◆ checkTrack() [2/2]

int HoughFinder::checkTrack ( vector< HoughTrack > &  trackVector)

◆ clearMemory() [1/2]

void HoughFinder::clearMemory ( )

Definition at line 2718 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2719{
2720 //cout<<"m_mdcHitCol.size()"<<m_mdcHitCol.size()<<endl;
2721 vector<MdcHit*>::iterator imdcHit = m_mdcHitCol.begin();
2722 for(;imdcHit != m_mdcHitCol.end();imdcHit++){
2723 delete (*imdcHit);
2724 }
2725
2726 for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin();trkIT!=m_houghTrackList.end();trkIT++){
2727 trkIT->clearMemory();
2728 //TrkRecoTrk* trkRecoTrk = trkIT->getTrkRecoTrk();
2729 //delete trkRecoTrk;
2730 }
2731}

Referenced by execute().

◆ clearMemory() [2/2]

void HoughFinder::clearMemory ( )

◆ clearTrack() [1/2]

void HoughFinder::clearTrack ( vector< HoughTrack > &  trackVector)

Definition at line 2705 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2706{
2707 /*
2708 vector<HoughTrack*>::iterator it = trackVector.begin();
2709 for(; it!=trackVector.end(); it++)
2710 {
2711 HoughTrack* trk = (*it);
2712 delete trk;
2713 }
2714 */
2715 trackVector.clear();
2716}

Referenced by execute().

◆ clearTrack() [2/2]

void HoughFinder::clearTrack ( vector< HoughTrack > &  trackVector)

◆ execute() [1/2]

StatusCode HoughFinder::execute ( )

Definition at line 281 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

282{
283 MsgStream log(msgSvc(), name());
284 log << MSG::INFO << "HoughFinder::execute()" << endreq;
285 cout.precision(6);
286
287
288 //cout<<"line "<<__LINE__<<endl;
289 // --- Get event header
290 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
291 if(!eventHeader){
292 log << MSG::ERROR << "Could not find Event Header" << endreq;
293 return StatusCode::FAILURE;
294 }
295 m_run = eventHeader->runNumber();
296 m_event = eventHeader->eventNumber();
297 m_totEvtProcessed++; if(m_totEvtProcessed%100==0) cout<<"HoughFinder: "<<m_totEvtProcessed<<" events processed."<<endl;
298
299 //cout<<"line "<<__LINE__<<endl;
300 // --- event start time
301 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
302 if(!evTimeCol){
303 log << MSG::ERROR << "Could not retrieve RecEsTimeCol" << endreq;
304 return StatusCode::FAILURE;
305 }
306 m_bunchT0 = 0;
307 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
308 if (iter_evt != evTimeCol->end()){
309 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
310 }
311
312 // --- clear, register and store reconstructed track
313 registerTrack(m_trackList_tds,m_hitList_tds);
314
315
316 //cout<<"line "<<__LINE__<<endl;
317 if(m_filter){
318 ifstream lowPt_Evt;
319 lowPt_Evt.open(m_evtFile.c_str());
320 vector<int> evtlist;
321 int evtNum;
322 while( lowPt_Evt >> evtNum) {
323 evtlist.push_back(evtNum);
324 }
325 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),m_event);
326 if(iter_evt == evtlist.end()){
327 setFilterPassed(false);
328 return StatusCode::SUCCESS;
329 }
330 else{
331 setFilterPassed(true);
332 }
333 }
334 //if(m_event<2000||m_event>3000){
335 //setFilterPassed(false);
336 //return StatusCode::SUCCESS;
337 //}
338
339 if(m_debug){
340 //cout<<endl;
341 cout<<"===================================================================================================="<<endl;
342 cout<<"run:"<<m_run<<" , event: "<<m_event<<endl;
343 //cout<<m_event<<" , ";
344 //cout<<endl;
345 }
346
347 //cout<<"line "<<__LINE__<<endl;
348 // --- prepare hits for reconstruction
349 int nHit = makeHoughHitList();
350
351 //cout<<"line "<<__LINE__<<endl;
352 int nMcHit;
353 if(m_run<0&&m_mcTruth){
354 int nMcHit = getMcHitCol();
355 //cout<<"nMcHit "<<nMcHit<<endl;
356if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
358if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
359 if(printFlag)cout<<endl;
360 if(printFlag)cout<<endl;
361 if(printFlag)cout<<endl;
362 //halfCircle(m_mcTrackCol,m_mcHitCol);
363 }
364
365 //cout<<"line "<<__LINE__<<endl;
366 searchCircle();
367 //cout<<"ntracks:"<<m_houghTrackList.size()<<endl;
368 //cout<<"line "<<__LINE__<<endl;
369 checkHot(m_houghTrackList);
370 //cout<<"line "<<__LINE__<<endl;
371 checkSharedHits(m_XHoughHitList);
372 //cout<<"line "<<__LINE__<<endl;
373 associateVHits();
374 //cout<<"line "<<__LINE__<<endl;
375if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
376 if(storeFlag==1)storeTracks(m_trackList_tds, m_hitList_tds, m_houghTrackList);
377 if(storeFlag==2)storeTrack(m_houghTrackList,m_trackList_tds,m_hitList_tds);
378if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
379 if(printFlag)cout<<endl;
380 if(printFlag)cout<<endl;
381 if(printFlag)cout<<endl;
382 //cout<<"line "<<__LINE__<<endl;
383
384 if(m_fillNTuple)dumpHit();
385 if(m_fillNTuple==1)dumpHoughTrack();
386 if(m_fillNTuple==1)dumpHoughEvent();
387 if(m_fillNTuple==2)dumpTdsTrack();
388 if(m_fillNTuple==2)dumpTdsEvent();
389 //cout<<"line "<<__LINE__<<endl;
390
391if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
393if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
394 if(printFlag)cout<<endl;
395 if(printFlag)cout<<endl;
396 if(printFlag)cout<<endl;
397
398if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
400if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
401 if(printFlag)cout<<endl;
402 if(printFlag)cout<<endl;
403 if(printFlag)cout<<endl;
404
405 clearTrack(m_houghTrackList);
406 //cout<<"line "<<__LINE__<<endl;
407
408 clearMemory();
409 //cout<<"line "<<__LINE__<<endl;
410 if(m_debug)cout<<endl;
411
412 return StatusCode::SUCCESS;
413}
int storeTrack(vector< HoughTrack > &trackVector, RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
int storeTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
void checkSharedHits(vector< HoughHit * > &hitList)
void clearTrack(vector< HoughTrack > &trackVector)
int checkHot(vector< HoughTrack > &trackVector)
StatusCode registerTrack(RecMdcTrackCol *&trackList_tds, RecMdcHitCol *&hitList_tds)

◆ execute() [2/2]

StatusCode HoughFinder::execute ( )

◆ fillHistogram() [1/6]

int HoughFinder::fillHistogram ( HoughHit hit,
TH2D *  hitMap,
HoughTrack trkCandi,
int  vote 
)

Definition at line 1002 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1003{
1004 int charge = trkCandi->getCharge();
1005 int nBin = hitMap->GetXaxis()->GetNbins();
1006 double xMin = hitMap->GetXaxis()->GetXmin();
1007 double xMax = hitMap->GetXaxis()->GetXmax();
1008 double yMin = hitMap->GetYaxis()->GetXmin();
1009 double yMax = hitMap->GetYaxis()->GetXmax();
1010 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
1011 double x = xMin + dx/2;
1012
1013 int flag = hit->getFlag();
1014 int nPoint(0);
1015 if(flag==0){
1016 bool firstHalf = trkCandi->judgeHalf(hit)==1;
1017 if(!firstHalf) return 0;
1018 double xHit = hit->getHitPosition().x();
1019 double yHit = hit->getHitPosition().y();
1020 double rHit = hit->getDriftDist();//drift distance
1021 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
1022
1023 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
1024 double X1 = 2*xHit/denominator;
1025 double Y1 = 2*yHit/denominator;
1026 double X2 = 2*xHit/denominator;
1027 double Y2 = 2*yHit/denominator;
1028 double R = 2*rHit/denominator;
1029
1030 while(x<xMax){
1031 double y1 = X1*cos(x) + Y1*sin(x) + R;
1032 double y2 = X2*cos(x) + Y2*sin(x) - R;
1033 double slope1 = -X1*sin(x) + Y1*cos(x);
1034 double slope2 = -X2*sin(x) + Y2*cos(x);
1035 double hitOnCircle1 = charge*slope1*y1;
1036 double hitOnCircle2 = charge*slope2*y2;
1037
1038 bool cut(true);
1039 cut = yMin<y1 && y1<yMax;
1040 //cut = cut && hitOnCircle1 <=0;
1041 if(cut){
1042 hitMap->Fill(x,y1);
1043 nPoint++;
1044 }
1045
1046 cut = yMin<y2 && y2<yMax;
1047 //cut = cut && hitOnCircle2 <= 0;
1048 if(cut){
1049 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
1050 //int bin = hitMap->FindBin(x,y2);
1051 //if(hitMap->GetBinContent(bin)>0)continue;
1052 hitMap->Fill(x,y2);
1053 nPoint++;
1054 }
1055 x += dx;
1056 }
1057 }else{
1058 //if(hit.getPairHit()==NULL)return nPoint;// FIXME
1059 //if(hit.getPairHit()->getHalfCircle()!=1)return nPoint;// FIXME
1060 int nTangency = trkCandi->calculateZ_S(hit);
1061 //cout<<nTangency<<endl;
1062 if(nTangency==0){
1063 //hit.setUse(-1);// why? FIXME
1064 }else{
1065 vector<HoughHit::S_Z> sz = hit->getSZ();
1066 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1067 for(;iter!=sz.end();iter++){
1068 double S = iter->first;
1069 double Z = iter->second;
1070 while(x<xMax){
1071 double y = -S*x + Z;
1072 //y = S*cos(x) + Z*sin(x);
1073 bool cut(true);
1074 cut = yMin<y && y<yMax;
1075 if(cut){
1076 hitMap->Fill(x,y);
1077 nPoint++;
1078 }
1079 x += dx;
1080 }
1081 }
1082 }
1083 }
1084 return nPoint;
1085}
Double_t x[10]
double sin(const BesAngle a)
double cos(const BesAngle a)
*********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
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27

◆ fillHistogram() [2/6]

int HoughFinder::fillHistogram ( HoughHit hit,
TH2D *  hitMap,
HoughTrack trkCandi,
int  vote 
)

◆ fillHistogram() [3/6]

int HoughFinder::fillHistogram ( HoughHit hit,
TH2D *  hitMap,
int  charge,
int  vote 
)

Definition at line 926 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

927{
928 int nBin = hitMap->GetXaxis()->GetNbins();
929 double xMin = hitMap->GetXaxis()->GetXmin();
930 double xMax = hitMap->GetXaxis()->GetXmax();
931 double yMin = hitMap->GetYaxis()->GetXmin();
932 double yMax = hitMap->GetYaxis()->GetXmax();
933 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
934 double x = xMin + dx/2;
935
936 int flag = hit->getFlag();
937 int nPoint(0);
938 if(flag==0){
939 double xHit = hit->getHitPosition().x();
940 double yHit = hit->getHitPosition().y();
941 double rHit = hit->getDriftDist();//drift distance
942 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
943
944 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
945 double X1 = 2*xHit/denominator;
946 double Y1 = 2*yHit/denominator;
947 double X2 = 2*xHit/denominator;
948 double Y2 = 2*yHit/denominator;
949 double R = 2*rHit/denominator;
950
951 while(x<xMax){
952 double y1 = X1*cos(x) + Y1*sin(x) + R;
953 double y2 = X2*cos(x) + Y2*sin(x) - R;
954 double slope1 = -X1*sin(x) + Y1*cos(x);
955 double slope2 = -X2*sin(x) + Y2*cos(x);
956 double hitOnCircle1 = charge*slope1*y1;
957 double hitOnCircle2 = charge*slope2*y2;
958
959 bool cut(true);
960 cut = yMin<y1 && y1<yMax;
961 cut = cut && hitOnCircle1 <=0;
962 if(cut){
963 hitMap->Fill(x,y1);
964 nPoint++;
965 }
966
967 cut = yMin<y2 && y2<yMax;
968 cut = cut && hitOnCircle2 <= 0;
969 if(cut){
970 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
971 //int bin = hitMap->FindBin(x,y2);
972 //if(hitMap->GetBinContent(bin)>0)continue;
973 hitMap->Fill(x,y2);
974 nPoint++;
975 }
976 x += dx;
977 }
978 }else{
979 vector<HoughHit::S_Z> sz = hit->getSZ();
980 vector<HoughHit::S_Z>::iterator iter = sz.begin();
981 for(;iter!=sz.end();iter++){
982 double S = iter->first;
983 double Z = iter->second;
984 double dy = -S*dx;
985 double y = -dy + Z;
986 while(x<xMax){
987 x += dx;
988 y += dy;
989 //y = S*cos(x) + Z*sin(x);
990 bool cut(true);
991 cut = yMin<y && y<yMax;
992 if(cut){
993 hitMap->Fill(x,y);
994 nPoint++;
995 }
996 }
997 }
998 }
999 return nPoint;
1000}

Referenced by fillHistogram().

◆ fillHistogram() [4/6]

int HoughFinder::fillHistogram ( HoughHit hit,
TH2D *  hitMap,
int  charge,
int  vote 
)

◆ fillHistogram() [5/6]

int HoughFinder::fillHistogram ( vector< HoughHit * > &  hitList,
TH2D *  hitMap,
int  charge,
int  vote,
HoughTrack trkCandi = NULL 
)

Definition at line 1087 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1088{
1089 int nHitsFilled = 0;
1090 std::vector<HoughHit*>::iterator iter = hitList.begin();
1091 for(; iter!= hitList.end(); iter++)
1092 {
1093 int used = (*iter)->getUse();
1094 //(*iter)->print();
1095 if(used==0) {
1096 int nPoint = 0;
1097 if(trkCandi!=NULL){
1098 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1099 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1100 nPoint=fillHistogram((*iter),hitMap,trkCandi,vote);
1101 }
1102 else{
1103 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1104 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1105 nPoint=fillHistogram(*iter,hitMap,charge,vote);
1106 }
1107 if(nPoint>0) nHitsFilled++;
1108 }
1109 }
1110 return nHitsFilled;
1111}
int fillHistogram(HoughHit *hit, TH2D *hitMap, int charge, int vote)

◆ fillHistogram() [6/6]

int HoughFinder::fillHistogram ( vector< HoughHit * > &  hitList,
TH2D *  hitMap,
int  charge,
int  vote,
HoughTrack trkCandi = NULL 
)

◆ finalize() [1/2]

StatusCode HoughFinder::finalize ( )

Definition at line 415 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

416{
417 MsgStream log(msgSvc(), name());
418 log << MSG::INFO << "HoughFinder::finalize()" << endreq;
419 return StatusCode::SUCCESS;
420}

◆ finalize() [2/2]

StatusCode HoughFinder::finalize ( )

◆ getHoughTrkIt() [1/2]

vector< HoughTrack >::iterator HoughFinder::getHoughTrkIt ( vector< HoughTrack > &  houghTrkList,
int  trkId 
)

Definition at line 2803 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2804{
2805 vector<HoughTrack>::iterator it = houghTrkList.begin();
2806 for(; it!=houghTrkList.end(); it++)
2807 {
2808 if(it->getTrkID()==trkId) break;
2809 }
2810 return it;
2811}

Referenced by checkSharedHits().

◆ getHoughTrkIt() [2/2]

vector< HoughTrack >::iterator HoughFinder::getHoughTrkIt ( vector< HoughTrack > &  houghTrkList,
int  trkId 
)

◆ getMcHitCol() [1/2]

int HoughFinder::getMcHitCol ( )

— put match RecCgemCluster into CemgMcHit

Definition at line 493 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

494{
495 MsgStream log(msgSvc(), name());
496 int nMcHit(0);
497 if(!(m_mcHitCol.empty()))m_mcHitCol.clear();
498
499 // --- get cgemMcHit
500 if(m_cgem){
501 SmartDataPtr<CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
502 if (!cgemMcHitCol){
503 log << MSG::ERROR << "Could not retrieve CgemMcHitCol" << endreq;
504 }else{
505 CgemMcHitCol::iterator cgemMcHitIt = cgemMcHitCol->begin();
506 for(;cgemMcHitIt != cgemMcHitCol->end();cgemMcHitIt++,nMcHit++){
507 const CgemMcHit* cgemMcHit = (*cgemMcHitIt);
508 HoughHit hit(cgemMcHit,m_bunchT0,nMcHit);
509 m_mcHitCol.push_back(hit);
510 HoughHit* phit = &(*m_mcHitCol.rbegin());
511
512 // --- match fired strip ID between CemgMcHit and RecCgemCluster
513 vector<int> xFireStripID;
514 vector<int> vFireStripID;
515 int sheet;
516 TArrayI identifier = phit->getCgemMcHit()->GetIdentifier();
517 for(int ii = 0; ii < identifier.GetSize(); ii++){
518 //const Identifier ident = Identifier(identifier.GetAt(ii));
519 const Identifier ident = Identifier(identifier.At(ii));
520 if(CgemID::is_xstrip(ident))xFireStripID.push_back(CgemID::strip(ident));
521 else vFireStripID.push_back(CgemID::strip(ident));
522 sheet = CgemID::sheet(ident);
523 //cout
524 //<< " hitID:" << phit->getHitID()
525 //<< " layerID=" << CgemID::layer(ident)
526 //<< " sheetID=" << CgemID::sheet(ident)
527 //<< " isXstrip=" << CgemID::is_xstrip(ident)
528 //<< " stripID=" << CgemID::strip(ident)
529 //<<endl;
530 }
531 sort(xFireStripID.begin(),xFireStripID.end());
532 sort(vFireStripID.begin(),vFireStripID.end());
533 //vector<int>::iterator uniqueX = unique(xFireStripID.begin(),xFireStrip ID.end());
534 //vector<int>::iterator uniqueV = unique(vFireStripID.begin(),vFireStrip ID.end());
535 //xFireStripID.erase(uniqueX,xFireStripID.end());
536 //vFireStripID.erase(uniqueV,vFireStripID.end());
537 //for(vector<int>::iterator i = vFireStripID.begin();i!=vFireStripID.end();i++)cout<<*i<<endl;
538 //cout<<endl;
539
540 bool findX(false), findV(false);
541 const RecCgemCluster* xCluster = NULL;
542 const RecCgemCluster* vCluster = NULL;
543 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
544 if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
545 const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
546 if(recCgemCluster->getlayerid() != phit->getLayer())continue;
547 if(!(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1))continue;
548 //if(recCgemCluster->getsheetid() != sheet)continue;
549 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
550 if(!findX){
551 findX = recCgemCluster->getflag() == 0
552 && xFireStripID.size() > 0
553 && recCgemCluster->getclusterflagb() == *(xFireStripID.begin())
554 && recCgemCluster->getclusterflage() == *(xFireStripID.rbegin());
555 if(findX){
556 xCluster = recCgemCluster;
557 hitIt->setPairHit(phit);
558 }
559 }
560 if(!findV){
561 findV = recCgemCluster->getflag() == 1
562 && vFireStripID.size() > 0
563 && recCgemCluster->getclusterflagb() == *(vFireStripID.begin())
564 && recCgemCluster->getclusterflage() == *(vFireStripID.rbegin());
565 if(findV){
566 vCluster = recCgemCluster;
567 hitIt->setPairHit(phit);
568 }
569 }
570
571 }
572
573 /// --- put match RecCgemCluster into CemgMcHit
574 if(findX&&findV){
575 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
576 if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
577 const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
578 if(recCgemCluster->getlayerid() != phit->getLayer())continue;
579 if(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1)continue;
580 //if(recCgemCluster->getsheetid() != sheet)continue;
581 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
582 if(recCgemCluster->getclusterflagb() == xCluster->getclusterid()
583 && recCgemCluster->getclusterflage() == vCluster->getclusterid()){
584 phit->setCgemCluster(recCgemCluster);
585 phit->setWire(recCgemCluster->getclusterid());
586 phit->setFlag(recCgemCluster->getflag());
587 phit->setPairHit(&(*hitIt));
588 hitIt->setPairHit(phit);
589 }
590 }
591 }
592 else if(findX&&!findV){
593 phit->setCgemCluster(xCluster);
594 phit->setWire(xCluster->getclusterid());
595 phit->setFlag(xCluster->getflag());
596 }
597 else if(!findX&&findV){
598 phit->setCgemCluster(vCluster);
599 phit->setWire(vCluster->getclusterid());
600 phit->setFlag(vCluster->getflag());
601 }
602 }
603 }
604 }
605
606 // --- get mdcMcHit
607 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
608 if(!mdcMcHitCol){
609 log << MSG::ERROR << "Could not retrieve MdcMcHitCol" << endreq;
610 }else{
611 MdcMcHitCol::iterator mdcMcHitIt = mdcMcHitCol->begin();
612 for(;mdcMcHitIt != mdcMcHitCol->end();mdcMcHitIt++,nMcHit++){
613 const MdcMcHit* mdcMcHit = (*mdcMcHitIt);
614 HoughHit hit(mdcMcHit,m_bunchT0,nMcHit);
615 m_mcHitCol.push_back(hit);
616 HoughHit* phit = &(*m_mcHitCol.rbegin());
617 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
618 if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
619 const MdcDigi* mdcDigi = hitIt->getDigi();
620 if(mdcDigi->identify() == phit->getMdcMcHit()->identify()){
621 phit->setDigi(mdcDigi);
622 phit->setPairHit(&(*hitIt));
623 hitIt->setPairHit(phit);
624 }
625 }
626 }
627 }
628
629 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
630 for(HitVector_Iterator mcHitIt = m_mcHitCol.begin(); mcHitIt != m_mcHitCol.end();mcHitIt++){
631 if(hitIt->getLayer()!=mcHitIt->getLayer())continue;
632 if(hitIt->getHitType() == HoughHit::CGEMHIT && mcHitIt->getHitType() == HoughHit::CGEMMCHIT){
633 if(mcHitIt->getCgemCluster()==NULL)continue;
634 int recClusterID(-1);
635 int mcClusterID(-2);
636 int recFlag = hitIt->getCgemCluster()->getflag();
637 int mcFlag = mcHitIt->getCgemCluster()->getflag();;
638 if(recFlag==mcFlag){
639 recClusterID = hitIt->getCgemCluster()->getclusterid();
640 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
641 }
642
643 if(recFlag==1&&mcFlag==2){
644 recClusterID = hitIt->getCgemCluster()->getclusterid();
645 mcClusterID = mcHitIt->getCgemCluster()->getclusterflage();
646 }
647
648 if(recFlag==0&&mcFlag==2){
649 recClusterID = hitIt->getCgemCluster()->getclusterid();
650 mcClusterID = mcHitIt->getCgemCluster()->getclusterflagb();
651 }
652
653 if(recFlag==2&&mcFlag==1){
654 recClusterID = hitIt->getCgemCluster()->getclusterflage();
655 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
656 }
657
658 if(recFlag==2&&mcFlag==0){
659 recClusterID = hitIt->getCgemCluster()->getclusterflagb();
660 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
661 }
662
663 if(recClusterID==mcClusterID){
664 hitIt->setPairHit(&(*mcHitIt));
665 hitIt->setCgemMcHit(mcHitIt->getCgemMcHit());
666 //mcHitIt->addPairHit(&(*hitIt));
667
668 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" "<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
669 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
670 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" ";
671 //cout<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
672 //cout<<recFlag<<" "<<mcFlag<<" ";
673 //cout<<recClusterID<<" "<<mcClusterID<<endl;
674
675 //cout<<hitIt->getHitPosition();
676 //cout<<endl;
677 //cout<<mcHitIt->getHitPosition();
678 //cout<<endl;
679 //cout<<endl;
680 }
681 }
682 else if(hitIt->getHitType() == HoughHit::MDCHIT && mcHitIt->getHitType() == HoughHit::MDCMCHIT){
683 if(mcHitIt->getDigi()==NULL)continue;
684 //if(hitIt->getLayer() == mcHitIt->getLayer() && hitIt->getWire() == mcHitIt->getWire())
685 if(hitIt->getDigi()->identify() == mcHitIt->getDigi()->identify())
686 {
687 hitIt->setPairHit(&(*mcHitIt));
688 hitIt->setMdcMcHit(mcHitIt->getMdcMcHit());
689 //mcHitIt->addPairHit(&(*hitIt));
690 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<endl;
691 //cout<<mcHitIt->getLayer()<<" "<<mcHitIt->getWire()<<endl;
692 //cout<<endl;
693 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
694 }
695 }
696 }
697 }
698 //for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
699 //if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
700 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<" ";
701 //if(hitIt->getPairHit()!=NULL){
702 //cout<<hitIt->getPairHit()->getHitID()<<" ";
703 //cout<<hitIt->getPairHit()->getLayer()<<" "<<hitIt->getPairHit()->getWire()<<" ";
704 //}
705 //cout<<endl;
706 //}
707
708 return nMcHit;
709}
vector< HoughHit >::iterator HitVector_Iterator
static int strip(const Identifier &id)
static int sheet(const Identifier &id)
static bool is_xstrip(const Identifier &id)
void setCgemCluster(const RecCgemCluster *cgemCluster)
virtual Identifier identify() const
Definition: RawData.cxx:15

Referenced by execute().

◆ getMcHitCol() [2/2]

int HoughFinder::getMcHitCol ( )

◆ getMcParticleCol() [1/2]

int HoughFinder::getMcParticleCol ( )

Definition at line 711 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

711 {
712 const int maxCell[43] = {40,44,48,56, 64,72,80,80,
713 76,76,88,88, 100,100,112,112, 128,128,140,140,
714 160,160,160,160, 176,176,176,176, 208,208,208,208, 240,240,240,240,
715 256,256,256,256, 288,288,288 };
716 m_mcTrackCol.clear();
717 MsgStream log(msgSvc(), name());
718 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
719 if (!mcParticleCol) {
720 log << MSG::ERROR << "Could not find McParticle" << endreq;
721 }else{
722 bool psipDecay(false);
723 McParticleCol::iterator iter_mc = mcParticleCol->begin();
724 for (;iter_mc != mcParticleCol->end(); iter_mc++){
725 int trackIndex = (*iter_mc)->trackIndex();
726 int pid = (*iter_mc)->particleProperty();
727 bool primaryParticle = (*iter_mc)->primaryParticle();
728 bool leafParticle = (*iter_mc)->leafParticle();
729 bool decayFromGenerator = (*iter_mc)->decayFromGenerator();
730 bool decayInFlight = (*iter_mc)->decayInFlight();
731 HepLorentzVector initialPosition = (*iter_mc)->initialPosition();
732 HepLorentzVector initialFourMomentum = (*iter_mc)->initialFourMomentum();
733 //HepLorentzVector finalPosition = (*iter_mc)->finalPosition();
734 //HepLorentzVector finalFourMomentum = (*iter_mc)->finalFourMomentum();
735 unsigned int statusFlags = (*iter_mc)->statusFlags();
736 //string process = (*iter_mc)->getProcess();
737 //if(primaryParticle)continue;
738 //if(!decayFromGenerator)continue;
739 //if(pid==100443)psipDecay = true;
740 //if(!psipDecay) continue;
741 //if(!(fabs(pid)==11||fabs(pid)==211))continue;
742
743 string name;
744 int charge(0);
745 if( pid == 0 ) {
746 //cout << "Wrong particle id" <<endl;
747 continue;
748 }else{
749 IPartPropSvc* p_PartPropSvc;
750 static const bool CREATEIFNOTTHERE(true);
751 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
752 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
753 log << MSG::ERROR << "Could not initialize Particle Properties Service" << endreq;
754 }
755 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
756 if( p_particleTable->particle(pid) ){
757 name = p_particleTable->particle(pid)->name();
758 charge = p_particleTable->particle(pid)->charge();
759 }else if( p_particleTable->particle(-pid) ){
760 name = "anti " + p_particleTable->particle(-pid)->name();
761 charge = (-1)*p_particleTable->particle(-pid)->charge();
762 }
763 }
764
765 //cout
766 //<<setw(8)<<trackIndex
767 //<<setw(12)<<name
768 //<<setw(8)<<pid
769 //<<setw(8)<<primaryParticle
770 //<<setw(8)<<leafParticle
771 //<<setw(8)<<decayFromGenerator
772 //<<setw(8)<<decayInFlight
773 //<<setw(8)<<statusFlags
774 //<<endl;
775
776 HoughTrack track(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
777 //HoughTrack* track = new HoughTrack(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
778 HepPoint3D pivot(0,0,0);
779 track.pivot(pivot);
780
781 int nHot(0);
782 vector<HoughHit*> hotList;
783 vector<HoughHit>& hitList = m_mcHitCol;
784 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
785 vector<int> trkID = hitIt->getTrkID();
786 if(trkID[0]==track.getTrkID()){
787 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
788 track.addHitPnt(&(*hitIt));
789 track.addVecStereoHitPnt(&(*hitIt));
790 hotList.push_back(&(*hitIt));
791 }
792 if(hitIt->getHitType()==HoughHit::MDCMCHIT){
793 if(hitIt->getFlag()==0)track.addHitPnt(&(*hitIt));
794 else track.addVecStereoHitPnt(&(*hitIt));
795 hotList.push_back(&(*hitIt));
796 }
797 nHot++;
798 }
799 }
800 if(nHot>0)m_mcTrackCol.push_back(track);
801
802 ///*
803 bool classifyHotByHotLayer(false);
804 bool classifyHotByTrackCenter(true);
805 if(classifyHotByHotLayer){
806 int lastHitLayer(0);
807 int deltaLayer(0);
808 int halfCircle(0);
809 HoughHit* lastLayerHit;
810 vector<HoughHit*> hitOnLayer;
811 int nMdcHot(0);
812 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
813 if((*hitIt)->getHitType()!=HoughHit::MDCHIT)continue;
814 if(nMdcHot==0){//the first hit on track
815 hitOnLayer.clear();
816 hitOnLayer.push_back(*hitIt);
817 lastLayerHit = *hitIt;
818 lastHitLayer = (*hitIt)->getLayer();
819 deltaLayer = 0;
820 halfCircle = 1;
821 }else{
822 if((*hitIt)->getLayer()==lastHitLayer){//find hits on the same layer
823 hitOnLayer.push_back(*hitIt);
824 lastHitLayer = (*hitIt)->getLayer();
825 }else{//another layer hit
826 if(deltaLayer*((*hitIt)->getLayer()-lastHitLayer)>=0){//not the turn layer
827 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
828 (*iter)->setHalfCircle(halfCircle);
829 lastLayerHit = *iter;
830 }
831 }else{//the turn layer
832 halfCircle++;
833 int maxWireGap(0);
834 int turnPoint(0);
835 int deltaWire(0);
836 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
837 deltaWire = fabs((*iter)->getWire()-lastLayerHit->getWire());
838 if(deltaWire>(maxCell[(*iter)->getLayer()])/2) deltaWire = maxCell[(*iter)->getLayer()] - deltaWire;
839 if(deltaWire>maxWireGap){
840 maxWireGap = deltaWire;
841 turnPoint = iter - hitOnLayer.begin();
842 }
843 lastLayerHit = *iter;
844 }
845 deltaWire = fabs((*hitIt)->getWire()-lastLayerHit->getWire());
846 if(deltaWire>(maxCell[(*hitIt)->getLayer()])/2) deltaWire = maxCell[(*hitIt)->getLayer()] - deltaWire;
847 if(deltaWire>maxWireGap){
848 maxWireGap = deltaWire;
849 turnPoint = hitOnLayer.size();
850 }
851
852 if(maxWireGap>2){
853 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
854 int shift = iter-hitOnLayer.begin();
855 if(shift<turnPoint)(*iter)->setHalfCircle(halfCircle-1);
856 else (*iter)->setHalfCircle(halfCircle);
857 }
858 }else{
859 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
860 int shift = iter-hitOnLayer.begin();
861 if(shift<hitOnLayer.size()/2)(*iter)->setHalfCircle(halfCircle-1);
862 else (*iter)->setHalfCircle(halfCircle);
863 }
864 }
865 }
866 hitOnLayer.clear();
867 hitOnLayer.push_back(*hitIt);
868 lastHitLayer = (*hitIt)->getLayer();
869 deltaLayer = (*hitIt)->getLayer()-lastHitLayer;
870 }
871
872 }
873 nMdcHot++;
874 }
875 }
876 //*/
877 ///*
878 if(classifyHotByTrackCenter){
879 int cgemHalf(0), mdcHalf(0);
880 int cgemSign(0), mdcSign(0);
881 int half(0), sign(0);
882 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
883 //(*hitIt)->print();
884 if((*hitIt)->getHitType() == HoughHit::CGEMMCHIT){
885 if(track.judgeHalf(*hitIt) != cgemSign){
886 cgemHalf++;
887 cgemSign = track.judgeHalf(*hitIt);
888 half = cgemHalf;
889 sign = cgemSign;
890 }
891 }
892 else if((*hitIt)->getHitType() == HoughHit::MDCMCHIT){
893 if(track.judgeHalf(*hitIt) != mdcSign){
894 mdcHalf++;
895 mdcSign = track.judgeHalf(*hitIt);
896 half = mdcHalf;
897 sign = mdcSign;
898 }
899 }
900 (*hitIt)->setHalfCircle(half);
901 }
902 }
903 if(hotList.size()>0){
904 if(printFlag)track.print();
905 //track.printHot();
906 //cout<<endl;
907 }
908 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++){
909 HoughHit* hitIt = *hotIt;
910 if(printFlag)hitIt->print();
911 //cout<<hitIt->getHitID();
912 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
913 //cout<<endl;
914 }
915 if(printFlag)if(hotList.size())cout<<endl;
916 }
917 sort(m_mcTrackCol.begin(),m_mcTrackCol.end(),largerPt);
918 //cout<<"N MCtrack: "<<m_mcTrackCol.size()<<endl;
919 for(vector<HoughTrack>::iterator trkIter=m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
920 //trkIter->print();
921 }
922 }
923 return m_mcTrackCol.size();
924}
bool largerPt(HoughTrack trk1, HoughTrack trk2)

Referenced by execute().

◆ getMcParticleCol() [2/2]

int HoughFinder::getMcParticleCol ( )

◆ initialize() [1/2]

StatusCode HoughFinder::initialize ( )

Definition at line 94 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

95{
96 MsgStream log(msgSvc(), name());
97 log << MSG::INFO << "HoughFinder::initialize()" << endreq;
98 StatusCode sc;
99
100 // --- RawData
101 IRawDataProviderSvc* irawDataProviderSvc;
102 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
103 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
104 if ( sc.isFailure() ){
105 log << MSG::ERROR << name()<<" Could not load RawDataProviderSvc!" << endreq;
106 return StatusCode::FAILURE;
107 }
108
109
110 // --- MDC Geometry
111 IMdcGeomSvc* imdcGeomSvc;
112 sc = service ("MdcGeomSvc", imdcGeomSvc);
113 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
114 if ( sc.isFailure() ){
115 log << MSG::ERROR << "Could not load MdcGeoSvc!" << endreq;
116 return StatusCode::FAILURE;
117 }
118
119 // --- MdcCalibFunSvc
120 IMdcCalibFunSvc* imdcCalibSvc;
121 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
122 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
123 if ( sc.isFailure() ){
124 log << MSG::ERROR << "Could not load MdcCalibFunSvc!" << endreq;
125 return StatusCode::FAILURE;
126 }
127
128 if(m_cgem){
129 // --- CGEM geometry
130 ICgemGeomSvc* icgemGeomSvc;
131 sc = service ("CgemGeomSvc", icgemGeomSvc);
132 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
133 if ( sc.isFailure() ){
134 log << MSG::ERROR << "Could not load CgemGeomSvc!" << endreq;
135 return StatusCode::FAILURE;
136 }
137
138 // --- CgemCalibFunSvc
139 ICgemCalibFunSvc* icgemCalibSvc;
140 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
141 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
142 if ( sc.isFailure() ){
143 log << MSG::ERROR << "Could not load CgemCalibFunSvc!" << endreq;
144 return StatusCode::FAILURE;
145 }
146
147 }
148
149 // --- MagneticFieldSvc
150 sc = service ("MagneticFieldSvc",m_pIMF);
151 if( sc.isFailure() ) {
152 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
153 return StatusCode::FAILURE;
154 }
155 m_bfield = new BField(m_pIMF);
156 //log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
157 m_trkContextEv = new TrkContextEv(m_bfield);
158
159 m_mdcDetector = new MdcDetector();
160 //read pdt
161 Pdt::readMCppTable(m_pdtFile);
162
163 // --- Mdc Print Service
164 //IMdcPrintSvc* imdcPrintSvc;
165 //sc = service ("MdcPrintSvc", imdcPrintSvc);
166 //m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
167 //if ( sc.isFailure() ){
168 //log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
169 //return StatusCode::FAILURE;
170 //}
171
172 // --- Bes Timer Service
173 //sc = service( "BesTimerSvc", m_timersvc);
174 //if( sc.isFailure() ) {
175 //log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
176 //return StatusCode::FAILURE;
177 //}
178 //m_timer_all = m_timersvc->addItem("Execution");
179 //m_timer_all->propName("nExecution");
180
181 HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
182 HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
183 HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
184 HoughHit::setCgemCalibFunSvc(m_cgemCalibFunSvc);
185 HoughHit::setMdcDetector(m_mdcDetector);
188
189 m_roughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",m_nBinTheta,0.,M_PI, m_nBinRho, -1.*m_rhoRange, m_rhoRange);
190 //m_tanlRange = 0.93/sqrt(1-0.93*0.93);//FIXME
191 m_tanlRange = 3.0;//FIXME
192 m_dzRange = 50;
193 m_roughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",m_nBinTanl,-1.*m_tanlRange,m_tanlRange,m_nBinDz,-1.*m_dzRange,m_dzRange);
194
195 const int nPt=12;
196 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};
197 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};
198 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};
199 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 };
200 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};
201 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 };
202 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};
203 m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
204 m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
205 m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
206 m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
207 m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
208 m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
209
210
211 int s(0),l(0),npt(0);
212 double momenta[100] = {0};
213 double cuts[2][43][100] = {0};
214 string cutFilePath = getenv("HOUGHTRANSALGROOT");
215 cutFilePath += "/share/cut.txt";
216 ifstream fin(cutFilePath.c_str(), ios::in);
217 cout<<"cutfiledir:"<<cutFilePath.c_str()<<endl;
218 //ifstream fin("/workfs/bes/huangzhen/workarea/CgemBoss/6.6.5.c/Reconstruction/HoughFinder/HoughFinder-00-00-13/share/cut.txt",ios::in);
219 if (!fin.good())
220 cout << "Error in " << __FILE__<<endl;
221 bool firstline(true);
222 string strcom;
223 while(std::getline(fin, strcom)){
224 if(firstline){
225 fin >> npt
226 >> momenta[0] >> momenta[1] >> momenta[2] >> momenta[3] >> momenta[4] >> momenta[5] >> momenta[6] >> momenta[7] >> momenta[8] >> momenta[9]
227 >> momenta[10] >> momenta[11] >> momenta[12] >> momenta[13] >> momenta[14] >> momenta[15] >> momenta[16] >> momenta[17] >> momenta[18] >> momenta[19]
228 >> momenta[20] >> momenta[21] >> momenta[22] >> momenta[23] >> momenta[24] >> momenta[25] >> momenta[26] >> momenta[27] >> momenta[28] >> momenta[29]
229 >> momenta[30] >> momenta[31] >> momenta[32] >> momenta[33] >> momenta[34] >> momenta[35] >> momenta[36] >> momenta[37] >> momenta[38] >> momenta[39]
230 >> momenta[40] >> momenta[41] >> momenta[42] >> momenta[43] >> momenta[44] >> momenta[45] >> momenta[46] >> momenta[47] >> momenta[48] >> momenta[49]
231 >> momenta[50] >> momenta[51] >> momenta[52] >> momenta[53] >> momenta[54] >> momenta[55] >> momenta[56] >> momenta[57] >> momenta[58] >> momenta[59]
232 >> momenta[60] >> momenta[61] >> momenta[62] >> momenta[63] >> momenta[64] >> momenta[65] >> momenta[66] >> momenta[67] >> momenta[68] >> momenta[69]
233 >> momenta[70] >> momenta[71] >> momenta[72] >> momenta[73] >> momenta[74] >> momenta[75] >> momenta[76] >> momenta[77] >> momenta[78] >> momenta[79]
234 >> momenta[80] >> momenta[81] >> momenta[82] >> momenta[83] >> momenta[84] >> momenta[85] >> momenta[86] >> momenta[87] >> momenta[88] >> momenta[89]
235 >> momenta[90] >> momenta[91] >> momenta[92] >> momenta[93] >> momenta[94] >> momenta[95] >> momenta[96] >> momenta[97] >> momenta[98] >> momenta[99]
236 ;
237 firstline = false;
238 continue;
239 }
240 fin >> s >> l;
241 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]
242 >> 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]
243 >> 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]
244 >> 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]
245 >> 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]
246 >> 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]
247 >> 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]
248 >> 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]
249 >> 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]
250 >> 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]
251 ;
252 if(fin.peek()<0)break;
253 }
254 for(int is=0;is<=1;is++){
255 for(int il=0;il<43;il++){
256 TGraph* gr = new TGraph();
257 int point(0);
258 //cout<<il<<" ";
259 for(int ipt=0;ipt<npt;ipt++){
260 //cout<<cuts[is][il][ipt]<<" ";
261 if(fabs(cuts[is][il][ipt])<1e-6)continue;
262 //cut[is][il]->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
263 gr->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
264 point++;
265 }
266 if(point>0)HoughTrack::m_cut[is][il] = gr;
267 //cout<<endl;
268 }
269 }
270
271
272 if(m_fillNTuple) sc = bookTuple();
273 if ( sc.isFailure() ){
274 log << MSG::FATAL << "Could not book Tuple !" << endreq;
275 return StatusCode::FAILURE;
276 }
277
278 return StatusCode::SUCCESS;
279}
TGraph * gr
XmlRpcServer s
Definition: HelloServer.cpp:11
static void setMdcGeomSvc(MdcGeomSvc *mdcGeomSvc)
static void setMdcCalibFunSvc(MdcCalibFunSvc *mdcCalibFunSvc)
static void setCgemCalibFunSvc(CgemCalibFunSvc *cgemCalibFunSvc)
static void setCgemGeomSvc(CgemGeomSvc *cgemGeomSvc)
static void setMdcDetector(const MdcDetector *mdcDetector)
static void readMCppTable(std::string filenm)

◆ initialize() [2/2]

StatusCode HoughFinder::initialize ( )

◆ makeHoughHitList() [1/2]

int HoughFinder::makeHoughHitList ( )

Definition at line 427 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

428{
429 MsgStream log(msgSvc(), name());
430 int nHit(0);
431 if(!(m_mdcHitCol.empty())) m_mdcHitCol.clear();
432 if(!(m_houghHitList.empty())) m_houghHitList.clear();
433 if(!(m_XHoughHitList.empty())) m_XHoughHitList.clear();
434 if(!(m_VHoughHitList.empty())) m_VHoughHitList.clear();
435
436 if(m_cgem){
437 // --- RecCgemCluster
438 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
439 if(!recCgemClusterCol){
440 log << MSG::ERROR << "Could not retrieve RecCgemClusterCol" << endreq;
441 }else{
442 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterCol->begin();
443 for (;cgemClusterIter != recCgemClusterCol->end(); cgemClusterIter++,nHit++){
444 const RecCgemCluster* recCgemCluster = (*cgemClusterIter);
445 //m_recCgemCluster = (*cgemClusterIter);
446 HoughHit hit(recCgemCluster, m_bunchT0, nHit);
447 m_houghHitList.push_back(hit);
448 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
449 }
450 }
451 }
452
453 uint32_t getDigiFlag = 0;
454 if(m_dropHot)getDigiFlag |= MdcRawDataProvider::b_dropHot;
455 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
456 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
457 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
458
459 MdcDigiVec::iterator mdcDigiIter = mdcDigiVec.begin();
460 for (;mdcDigiIter != mdcDigiVec.end(); mdcDigiIter++,nHit++){
461 const MdcDigi* mdcDigi = (*mdcDigiIter);
462 HoughHit hit(mdcDigi, m_bunchT0, nHit);
463 if(fabs(hit.driftTime())>1500.)continue;
464 MdcHit *mdcHit= new MdcHit(mdcDigi,m_mdcDetector);
465 hit.setMdcHit(mdcHit);
466
467 m_mdcHitCol.push_back(mdcHit);
468 m_houghHitList.push_back(hit);
469 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
470 }
471
472 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
473 int flag=hitIt->getFlag();
474 if(hitIt->getHitType()==HoughHit::CGEMHIT){
475 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
476 if(flag==2)m_VHoughHitList.push_back(&(*hitIt));
477 }
478 if(hitIt->getHitType()==HoughHit::MDCHIT){
479 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
480 else m_VHoughHitList.push_back(&(*hitIt));
481 }
482 }
483 //for(vector<HoughHit*>::iterator hitIt = m_XHoughHitList.begin(); hitIt != m_XHoughHitList.end(); hitIt++){
484 //(*hitIt)->print();
485 //}
486 //for(vector<HoughHit*>::iterator hitIt = m_VHoughHitList.begin(); hitIt != m_VHoughHitList.end(); hitIt++){
487 //(*hitIt)->print();
488 //}
489
490 return nHit;
491}

Referenced by execute().

◆ makeHoughHitList() [2/2]

int HoughFinder::makeHoughHitList ( )

◆ printHitList() [1/2]

void HoughFinder::printHitList ( )

Definition at line 2329 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

2330{
2331 //cout<<"nHit ==========> "<<m_houghHitList.size()<<endl;
2332 cout<<"========== truth =========="<<endl;
2333 for(HitVector_Iterator hitIt=m_mcHitCol.begin();hitIt!=m_mcHitCol.end();hitIt++){
2334 hitIt->print();
2335 }
2336 cout<<endl;
2337 cout<<"========== Digi =========="<<endl;
2338 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++)
2339 {
2340 hitIt->print();
2341 }
2342 cout<<endl;
2343 /*
2344 {
2345 vector<HoughHit>& hitList = m_mcHitCol;
2346 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2347 cout<<setw(4)<<hitIt->getHitID();
2348 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2349 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2350 vector<int> trkID = hitIt->getTrkID();
2351 cout<<"trkID:"<<setw(4)<<trkID[0];
2352 cout<<"half:"<<setw(4)<<hitIt->getHalfCircle();
2353 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
2354 if(hitIt->getCgemCluster()!=NULL){
2355 //cout<<hitIt->getFlag()<<" ";
2356 //if(hitIt->getFlag()==2)
2357 cout<<setw(5)<<hitIt->getCgemCluster()->getclusterid();
2358 cout<<"["<<setw(3)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(3)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2359 }
2360 //else cout<<"NULL ";
2361 }
2362 else{
2363 if(hitIt->getDigi()!=NULL){
2364 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2365 cout<<hitIt->getDigi()->identify();
2366 }
2367 //else cout<<"NULL ";
2368
2369 }
2370 cout<<endl;
2371 }
2372 }
2373 cout<<endl;
2374 //*/
2375 /*
2376 {
2377 vector<HoughHit>& hitList = m_houghHitList;
2378 cout<<hitList.size()<<endl;
2379 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2380 cout<<setw(4)<<hitIt->getHitID();
2381 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2382 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2383 cout<<setw(4)<<hitIt->getFlag();
2384 if(hitIt->getHitType()==HoughHit::CGEMHIT){
2385 cout<<setw(4)<<hitIt->getCgemCluster()->getclusterid();
2386 cout<<"["<<setw(4)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2387 }
2388 else{
2389 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2390 cout<<hitIt->getDigi()->identify()<<" ";
2391 //cout<<"("<<hitIt->getPairHit()->getLayer()<<", "<<hitIt->getPairHit()->getWire()<<") ";
2392 //if(hitIt->getPairHit()!=NULL)cout<<hitIt->getDriftDist()<<" "<<hitIt->getPairHit()->getDriftDist()<<" "<<hitIt->getDriftDist()-hitIt->getPairHit()->getDriftDist()<<endl;
2393 }
2394 if(hitIt->getPairHit()!=NULL)cout<<hitIt->getPairHit()->getHitID();
2395 //else cout<<"NULL";
2396 cout<<endl;
2397 }
2398 }
2399 cout<<endl;
2400 //*/
2401}

Referenced by execute().

◆ printHitList() [2/2]

void HoughFinder::printHitList ( )

◆ printTdsTrack() [1/2]

int HoughFinder::printTdsTrack ( )

Definition at line 3897 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

3898{
3899 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
3900 if(!recMdcTrackCol)return 0;
3901 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
3902 double dr = (*iter)->helix(0) ;
3903 double phi0 = (*iter)->helix(1) ;
3904 double kappa = (*iter)->helix(2) ;
3905 double dz = (*iter)->helix(3) ;
3906 double tanl = (*iter)->helix(4) ;
3907 cout<<setw(12)<<"dr:" <<setw(15)<<dr
3908 <<setw(12)<<"phi0:" <<setw(15)<<phi0
3909 <<setw(12)<<"kappa:" <<setw(15)<<kappa
3910 <<setw(12)<<"dz:" <<setw(15)<<dz
3911 <<setw(12)<<"tanl:" <<setw(15)<<tanl
3912 <<endl;
3913
3914 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
3915 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
3916 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
3917 //cout<<(*clusterIt)->;
3918 }
3919
3920 HitRefVec hitRefVec = (*iter)->getVecHits();
3921 HitRefVec::iterator hotIt = hitRefVec.begin();
3922 //cout<<"("<<"layer"<<","<<"wire"<<") "
3923 cout<<"("<<"l "<<","<<" w "<<") "
3924 //<<"Stat "
3925 //<<"FlagLR "
3926 <<"S "
3927 <<"F "
3928 <<"DriftLeft "
3929 <<" DriftRight "
3930 <<"ErrDrift_L "
3931 <<"ErrDrift_R "
3932 <<"ChisqAdd "
3933 //<<" Tdc "
3934 //<<" Adc "
3935 <<" DriftT "
3936 <<" Doca "
3937 <<" Entra "
3938 <<" Zhit "
3939 <<" FltLen "
3940 <<endl;
3941 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
3942 int layer = MdcID::layer((*hotIt)->getMdcId());
3943 int wire = MdcID::wire((*hotIt)->getMdcId());
3944 cout<<"("<<setw( 2)<<layer<<","<<setw( 3)<<wire<<") "
3945 <<setw( 3)<<(*hotIt)->getStat()
3946 <<setw( 3)<<(*hotIt)->getFlagLR()
3947 <<setw(12)<<(*hotIt)->getDriftDistLeft()
3948 <<setw(12)<<(*hotIt)->getDriftDistRight()
3949 <<setw(12)<<(*hotIt)->getErrDriftDistLeft()
3950 <<setw(12)<<(*hotIt)->getErrDriftDistRight()
3951 <<setw(12)<<(*hotIt)->getChisqAdd()
3952 //<<setw( 6)<<(*hotIt)->getTdc()
3953 //<<setw( 6)<<(*hotIt)->getAdc()
3954 <<setw(10)<<(*hotIt)->getDriftT()
3955 <<setw(12)<<(*hotIt)->getDoca()
3956 <<setw(12)<<(*hotIt)->getEntra()
3957 <<setw(10)<<(*hotIt)->getZhit()
3958 <<setw(8 )<<(*hotIt)->getFltLen()
3959 <<endl;
3960 }
3961 cout<<endl;
3962 }
3963 return recMdcTrackCol->size();
3964}
SmartRefVector< RecCgemCluster > ClusterRefVec
SmartRefVector< RecMdcHit > HitRefVec
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)

Referenced by execute().

◆ printTdsTrack() [2/2]

int HoughFinder::printTdsTrack ( )

◆ registerTrack() [1/2]

StatusCode HoughFinder::registerTrack ( RecMdcTrackCol *&  trackList_tds,
RecMdcHitCol *&  hitList_tds 
)

Definition at line 1748 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1749{
1750 MsgStream log(msgSvc(), name());
1751 StatusCode sc;
1752 IDataProviderSvc* eventSvc = NULL;
1753 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1754 if (eventSvc) {
1755 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1756 } else {
1757 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1758 return StatusCode::FAILURE;
1759 //return StatusCode::SUCCESS;
1760 }
1761 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1762
1763
1764 // --- register RecMdcTrack
1765 recMdcTrackCol = new RecMdcTrackCol;
1766 DataObject *aRecMdcTrackCol;
1767 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1768 if(aRecMdcTrackCol != NULL) {
1769 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1770 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1771 }
1772
1773 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1774 if(sc.isFailure()) {
1775 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1776 return StatusCode::FAILURE;
1777 }
1778 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1779
1780
1781
1782 recMdcHitCol = new RecMdcHitCol;
1783 DataObject *aRecMdcHitCol;
1784 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1785 if(aRecMdcHitCol != NULL) {
1786 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1787 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1788 }
1789
1790 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1791 if(sc.isFailure()) {
1792 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1793 return StatusCode::FAILURE;
1794 }
1795 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1796
1797
1798 return StatusCode::SUCCESS;
1799}//end of stroeTracks
ObjectVector< RecMdcHit > RecMdcHitCol
ObjectVector< RecMdcTrack > RecMdcTrackCol

Referenced by execute().

◆ registerTrack() [2/2]

StatusCode HoughFinder::registerTrack ( RecMdcTrackCol *&  trackList_tds,
RecMdcHitCol *&  hitList_tds 
)

◆ separateHoughHitByPhi() [1/2]

int HoughFinder::separateHoughHitByPhi ( )

◆ separateHoughHitByPhi() [2/2]

int HoughFinder::separateHoughHitByPhi ( )

◆ storeTrack() [1/2]

int HoughFinder::storeTrack ( vector< HoughTrack > &  trackVector,
RecMdcTrackCol *&  recMdcTrackCol,
RecMdcHitCol *&  recMdcHitCol 
)

Definition at line 1801 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1802{
1803 MdcTrackParams mdcTrackParams;
1804 //mdcTrackParams.lPrint=8;
1805 mdcTrackParams.lRemoveInActive=1;
1806 mdcTrackParams.maxGapLength=99;
1807 mdcTrackParams.nHitDeleted=99;
1808 MdcTrackList mdcTrackList(mdcTrackParams);
1809 mdcTrackList.setNoInner(true);
1810 vector<MdcTrack*> mdcTrackVector;
1811 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
1812 if((trkIT)->getFlag() == 1)continue;
1813 if((trkIT)->getFlag() == -1)continue;
1814 if((trkIT)->getFlag() == -2)continue;
1815 //(trkIT)->sortHot();
1816 MdcTrack* mdcTrack = new MdcTrack((trkIT)->getTrkRecoTrk());
1817 mdcTrackList.append(mdcTrack);
1818 mdcTrackVector.push_back(mdcTrack);
1819 }
1820 int nDeleted(0);
1821 mdcTrackList.setNoInner(true);
1822 if(m_checkHits)nDeleted = mdcTrackList.arbitrateHits();
1823 if(nDeleted>0)cout<<"nDeleted "<<nDeleted<<endl;
1824
1825 int trkID=0;
1826 for(int l=0;l<mdcTrackList.length();l++){
1827 mdcTrackList[l]->storeTrack(trkID,recMdcTrackCol,recMdcHitCol,4);
1828 trkID++;
1829 }
1830
1831 vector<MdcTrack*>::iterator iter = mdcTrackVector.begin();
1832 for(;iter!=mdcTrackVector.end();iter++){
1833 if(m_clearTrack)delete *iter;
1834 }
1835 return trkID;
1836}

Referenced by execute().

◆ storeTrack() [2/2]

int HoughFinder::storeTrack ( vector< HoughTrack > &  trackVector,
RecMdcTrackCol *&  recMdcTrackCol,
RecMdcHitCol *&  recMdcHitCol 
)

◆ storeTracks() [1/2]

int HoughFinder::storeTracks ( RecMdcTrackCol trackList,
RecMdcHitCol hitList,
vector< HoughTrack > &  trackVector 
)

Definition at line 1842 of file HoughTransAlg/HoughTransAlg-00-00-14/src/Hough.cxx.

1842 {
1843 int trackId(0);
1844 int tkStat =4;
1845 //cout<<"N Rectrack: "<<trackVector.size()<<endl;
1846 sort(trackVector.begin(),trackVector.end(),largerPt);
1847
1848 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
1849 if((trkIT)->getFlag() == 1)continue;
1850 if((trkIT)->getFlag() == -1)continue;
1851 if((trkIT)->getFlag() == -2)continue;
1852 if(trkIT->getTrkRecoTrk()==NULL)continue;
1853 //get the results of the fit to this track
1854 const TrkFit* fitresult = trkIT->getTrkRecoTrk()->fitResult();
1855 //check the fit worked
1856 if (fitresult == 0) continue;
1857 if(printFlag)trkIT->print();
1858 if(printFlag)trkIT->printHot();
1859 //cout<<endl;
1860
1861 //new a Rec. Track MdcTrack
1862 RecMdcTrack* recMdcTrack = new RecMdcTrack();
1863 const TrkHitList* aList = trkIT->getTrkRecoTrk()->hits();
1864 //some track info
1865 const BField& theField = trkIT->getTrkRecoTrk()->bField();
1866 double Bz = theField.bFieldZ();
1867 //Fit info
1868 int hitId = 0;
1869 int nHits=0;
1870 int nSt=0;
1871 //int nAct=0; int nSeg=0;
1872 //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
1873 //****************************
1874 //* active flag open this
1875 //****************************
1876 //* if (allHit>0){
1877 //* nHits = aList->nHit();//add inActive hit also
1878 //* }else{
1879 //* nHits = fitresult->nMdc();// = nActive()
1880 //* }
1881 //****************************
1882 //for 5.1.0 use all hits
1883 nHits = aList->nHit();
1884 //****************************
1885 // use active only
1886 // nHits = fitresult->nMdc();// = nActive()
1887 //****************************
1888
1889 int q = fitresult->charge();
1890 double chisq = fitresult->chisq();
1891 int nDof = fitresult->nDof();
1892 //track parameters at closest approach to beamline
1893 double fltLenPoca = 0.0;
1894 TrkExchangePar helix = fitresult->helix(fltLenPoca);
1895 //std::cout<< __FILE__ << " phi0 " << helix.phi0() << " omega " <<helix.omega()<<std::endl;
1896 double phi0 = helix.phi0();
1897 double tanDip = helix.tanDip();
1898 //double z0 = helix.z0();
1899 double d0 = helix.d0();
1900 //momenta and position
1901 //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
1902 HepPoint3D poca = fitresult->position(fltLenPoca);
1903
1904 double helixPar[5];
1905 //dro =-d0
1906 helixPar[0] = -d0;
1907 //phi0 = phi0 - pi/2 [0,2pi)
1908 double tphi0 = phi0 - Constants::pi/2.;
1909 if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
1910 helixPar[1] = tphi0;
1911 //kappa = q/pxy
1912 double pxy = fitresult->pt();
1913 if (pxy == 0.) helixPar[2] = 9999.;
1914 else helixPar[2] = q/fabs(pxy);
1915 if(pxy>9999.) helixPar[2] = 0.00001;
1916 //dz = z0
1917 helixPar[3] = helix.z0();
1918 //tanl
1919 helixPar[4] = tanDip;
1920 //error
1921 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
1922 HepSymMatrix mS(helix.params().num_row(),0);
1923 mS[0][0]=-1.;//dr0=-d0
1924 mS[1][1]=1.;
1925 mS[2][2]=-333.567/Bz;//pxy -> cpa
1926 mS[3][3]=1.;
1927 mS[4][4]=1.;
1928 HepSymMatrix mVy = helix.covariance().similarity(mS);
1929 double errorMat[15];
1930 int k = 0;
1931 for (int ie = 0 ; ie < 5 ; ie ++){
1932 for (int je = ie ; je < 5 ; je ++) {
1933 errorMat[k] = mVy[ie][je];
1934 k++;
1935 }
1936 }
1937 double p,px,py,pz;
1938 px = pxy * (-sin(helixPar[1]));
1939 py = pxy * cos(helixPar[1]);
1940 pz = pxy * helixPar[4];
1941 p = sqrt(pxy*pxy + pz*pz);
1942 //theta, phi
1943 double theta = acos(pz/p);
1944 double phi = atan2(py,px);
1945 recMdcTrack->setTrackId(trackId);
1946 recMdcTrack->setHelix(helixPar);
1947 recMdcTrack->setCharge(q);
1948 recMdcTrack->setPxy(pxy);
1949 recMdcTrack->setPx(px);
1950 recMdcTrack->setPy(py);
1951 recMdcTrack->setPz(pz);
1952 recMdcTrack->setP(p);
1953 recMdcTrack->setTheta(theta);
1954 recMdcTrack->setPhi(phi);
1955 recMdcTrack->setPoca(poca);
1956 recMdcTrack->setX(poca.x());//poca
1957 recMdcTrack->setY(poca.y());
1958 recMdcTrack->setZ(poca.z());
1959 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
1960 HepPoint3D pivot(0.,0.,0.);
1961 recMdcTrack->setPivot(pivot);
1962 recMdcTrack->setVX0(0.);//pivot
1963 recMdcTrack->setVY0(0.);
1964 recMdcTrack->setVZ0(0.);
1965 recMdcTrack->setError(errorMat);
1966 recMdcTrack->setError(mVy);
1967 recMdcTrack->setChi2(chisq);
1968 recMdcTrack->setStat(tkStat);
1969 recMdcTrack->setNhits(nHits);
1970 //-----------hit list-------------
1971 HitRefVec hitRefVec;
1972 ClusterRefVec clusterRefVec;
1973 vector<int> vecHits;
1974 map<int,int> clusterFitStat;
1975 //for fiterm
1976 int maxLayer = 0;
1977 //cout<<__FILE__ <<" "<<__LINE__ <<endl;
1978 for(TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
1979 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
1980 const MdcRecoHitOnTrack* recoHot;
1981 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
1982 int layerId = recoHot->mdcHit()->layernumber();
1983 int wireId = recoHot->mdcHit()->wirenumber();
1984 if(maxLayer < layerId)maxLayer = layerId;
1985 }
1986 }
1987 //cout<<"maxLayer:"<<maxLayer<<endl;
1988 int maxLayerId = -1;
1989 int minLayerId = 43;
1990 double fiTerm = 999.;
1991 const MdcRecoHitOnTrack* fiTermHot = NULL;
1992 TrkHitList::hot_iterator hot = aList->begin();
1993 for (;hot!=aList->end();hot++){
1994 //cout<<__FILE__ <<" "<<__LINE__ <<" "<< typeid(*hot).name()<<endl;
1995 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
1996 const MdcRecoHitOnTrack* recoHot;
1997 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
1998 //if (!recoHot->mdcHit()) return;
1999 RecMdcHit* recMdcHit = new RecMdcHit;
2000 //id
2001 recMdcHit->setId(hitId);
2002 //---------BESIII----------
2003 //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
2004 //phiWire - phiHit >0 flagLR=1 right driftright>0;
2005 //flag = 2 ambig;
2006 //-----Babar wireAmbig----
2007 //-1 -> right, 0 -> don't know, +1 -> left
2008 //+1 phiWire-phiHit<0 Left,
2009 //-1 phiWire-phiHit>0 Right,
2010 //0 don't know
2011 //ambig() w.r.t track direction
2012 //wireAmbig() w.r.t. wire location
2013 double hotWireAmbig = recoHot->wireAmbig();
2014 double driftDist = fabs(recoHot->drift());
2015 double sigma = recoHot->hitRms();
2016 double doca = fabs(recoHot->dcaToWire());
2017 int flagLR=2;
2018 if ( hotWireAmbig == 1){
2019 flagLR = 0; //left driftDist <0
2020 doca *= -1.;//2012-07-19
2021 }else if( hotWireAmbig == -1){
2022 flagLR = 1; //right driftDist >0
2023 }else if( hotWireAmbig == 0){
2024 flagLR = 2; //don't know
2025 }
2026 recMdcHit->setFlagLR(flagLR);
2027 recMdcHit->setDriftDistLeft((-1)*driftDist);
2028 recMdcHit->setErrDriftDistLeft(sigma);
2029 recMdcHit->setDriftDistRight(driftDist);
2030 recMdcHit->setErrDriftDistRight(sigma);
2031 //Mdc Id
2032 Identifier mdcId = recoHot->mdcHit()->digi()->identify();
2033 recMdcHit->setMdcId(mdcId);
2034 //corrected ADC
2035
2036 //contribution to chisquare
2037 //fill chisq
2038 double res=999.,rese=999.;
2039 if (recoHot->resid(res,rese,false)){
2040 }else{}
2041 double deltaChi=0;
2042 recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20
2043 recMdcHit->setChisqAdd( deltaChi * deltaChi );
2044 //set tracks containing this hit
2045 recMdcHit->setTrkId(trackId);
2046 //doca
2047 recMdcHit->setDoca(doca);//doca sign left<0
2048 //entrance angle
2049
2050 recMdcHit->setEntra(recoHot->entranceAngle());
2051 //z of hit
2052 HepPoint3D pos; Hep3Vector dir;
2053 double fltLen = recoHot->fltLen();
2054 recMdcHit->setFltLen(fltLen);
2055 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
2056 recMdcHit->setZhit(pos.z());
2057 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
2058 recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
2059 recMdcHit->setAdc(recoHot->mdcHit()->adcIndex());
2060 //drift time
2061 recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
2062 //for fiterm
2063 int layerId = recoHot->mdcHit()->layernumber();
2064 int wireId = recoHot->mdcHit()->wirenumber();
2065 if(layerId>(maxLayer - m_removeNOuterHits))continue;
2066 //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig "<<hotWireAmbig<<" y "<<pos.y()<<std::endl;
2067 // <<recoHot->entranceAngle()<<std::endl;
2068 if (layerId >= maxLayerId){
2069 maxLayerId = layerId;
2070 fiTermHot = recoHot;
2071 }
2072 if (layerId < minLayerId){
2073 minLayerId = layerId;
2074 }
2075 // status flag
2076 if (recoHot->isActive()) {
2077 recMdcHit->setStat(1);
2078 }else{
2079 recMdcHit->setStat(0);
2080 }
2081 // for 5.1.0 use all hits
2082 if (recoHot->layer()->view()) {
2083 ++nSt;
2084 }
2085 hitList->push_back(recMdcHit);
2086 SmartRef<RecMdcHit> refHit(recMdcHit);
2087 hitRefVec.push_back(refHit);
2088 vecHits.push_back(mdcId.get_value());
2089 ++hitId;
2090 }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2091 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2092 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2093 int clusterid = recCgemCluster->getclusterid();
2094 int stat = recoHot->isActive();
2095 //cout<<"stat:"<<stat<<endl;
2096 clusterFitStat[clusterid] = stat;
2097 //cout<<clusterid<<" "<<stat<<endl;
2098 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2099 clusterRefVec.push_back(recCgemCluster);
2100 }
2101 }
2102 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2103 //fi terminal angle
2104 if (fiTermHot!=NULL){
2105 fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
2106 }
2107 recMdcTrack->setFiTerm(fiTerm);
2108 // number of stereo hits contained
2109 recMdcTrack->setNdof(nDof);
2110 recMdcTrack->setNster(nSt);
2111 //recMdcTrack->setFirstLayer(maxLayerId);
2112 //recMdcTrack->setLastLayer(minLayerId);
2113 recMdcTrack->setFirstLayer(minLayerId);
2114 recMdcTrack->setLastLayer(maxLayerId);
2115 recMdcTrack->setVecHits(hitRefVec);
2116 //recMdcTrack->setClusterFitStat(clusterFitStat);
2117 //recMdcTrack->setVecClusters(clusterRefVec);
2118 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2119 recMdcTrack->setNcluster(clusterRefVec.size());
2120 trackList->push_back(recMdcTrack);
2121 trackId++;
2122 }
2123//cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
2124/*
2125 int trackId(0);
2126 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2127 //TrkErrCode trkErrCode = trkIT->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
2128 HoughTrack* houghTrack = &(*trkIT);
2129 if(houghTrack->getFlag() == 1)continue;
2130 if(houghTrack->getFlag() < 0)continue;
2131 //houghTrack->print();
2132 //houghTrack->printHot();
2133 trackId=trackList->size();
2134 //cout<<"start store trk "<<trackId;
2135 RecMdcTrack* recMdcTrack = new RecMdcTrack();
2136 //vector<HoughHit>* hotList = houghTrack->getHot();
2137 //vector<HoughHit*> hotList = houghTrack->getVecHitPnt();
2138 vector<HoughHit*> hotList = houghTrack->getHotList();
2139 //cout<<" with "<<hotList.size()<<" hits";//<<endl;
2140 double helixPar[5];
2141 helixPar[0] = houghTrack->dr();
2142 helixPar[1] = houghTrack->phi0();
2143 helixPar[2] = houghTrack->kappa();
2144 helixPar[3] = houghTrack->dz();
2145 helixPar[4] = houghTrack->tanl();
2146 int q = houghTrack->getCharge();
2147 double pxy = houghTrack->pt();
2148 double px = houghTrack->momentum(0).x();
2149 double py = houghTrack->momentum(0).y();
2150 double pz = houghTrack->momentum(0).z();
2151 double p = houghTrack->momentum(0).mag();
2152 double theta = houghTrack->direction(0).theta();
2153 double phi = houghTrack->direction(0).phi();
2154 HepPoint3D poca = houghTrack->x(0);
2155 HepPoint3D pivot = houghTrack->pivot();
2156 double r = poca.perp();
2157 HepSymMatrix Ea = houghTrack->Ea();
2158 //cout<<" phi, pxy, pz = "<<phi<<", "<<pxy<<", "<<pz<<endl;
2159 double errorMat[15];
2160 int k = 0;
2161 for (int ie = 0 ; ie < 5 ; ie ++){
2162 for (int je = ie ; je < 5 ; je ++){
2163 errorMat[k] = Ea[ie][je];
2164 k++;
2165 }
2166 }
2167 double chisq = houghTrack->getChi2();
2168 int tkStat = 4+houghTrack->getCircleFitStat();
2169 //if(tkStat==4) cout<<" track "<<trackId<<" has no good circle fit"<<endl;
2170 int nHits = hotList.size();
2171
2172 recMdcTrack->setTrackId(trackId);
2173 recMdcTrack->setHelix(helixPar);
2174 recMdcTrack->setCharge(q);
2175 recMdcTrack->setPxy(pxy);
2176 recMdcTrack->setPx(px);
2177 recMdcTrack->setPy(py);
2178 recMdcTrack->setPz(pz);
2179 recMdcTrack->setP(p);
2180 recMdcTrack->setTheta(theta);
2181 recMdcTrack->setPhi(phi);
2182 recMdcTrack->setPoca(poca);
2183 recMdcTrack->setX(poca.x());//poca
2184 recMdcTrack->setY(poca.y());
2185 recMdcTrack->setZ(poca.z());
2186 recMdcTrack->setR(r);
2187 recMdcTrack->setPivot(pivot);
2188 recMdcTrack->setVX0(0.);//pivot
2189 recMdcTrack->setVY0(0.);
2190 recMdcTrack->setVZ0(0.);
2191 recMdcTrack->setError(errorMat);
2192 recMdcTrack->setError(Ea);
2193 recMdcTrack->setChi2(chisq);
2194 recMdcTrack->setStat(tkStat);
2195 recMdcTrack->setNhits(nHits);
2196 //-----------hit list-------------
2197 int hitId = 0;
2198 int nSt(0);
2199 int maxLayerId = -1;
2200 int minLayerId = 43;
2201 double fiTerm = 999.;
2202 HoughHit fiTermHot;
2203
2204 HitRefVec hitRefVec;
2205 ClusterRefVec clusterRefVec;
2206 map<int,int> clusterFitStat;
2207 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2208 {
2209 //cout<<"handle hit "<<hitId<<endl;
2210 HoughHit* hot(*hotIt);
2211 //int layer = (**hotIt).getLayer();
2212 //int wire = (**hotIt).getWire();
2213 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2214 if(hot->getHitType()==HoughHit::MDCHIT){
2215 //cout<<"create a new RecMdcHit"<<endl;
2216 RecMdcHit* recMdcHit = new RecMdcHit;
2217
2218 recMdcHit->setId(hitId);
2219 recMdcHit->setTrkId(trackId);
2220 HepPoint3D hitPoint = hot->getHitPosition();
2221 double dPhi = houghTrack->dPhi(hitPoint);
2222 HepPoint3D position = houghTrack->x(dPhi);
2223 Hep3Vector direction = houghTrack->direction(dPhi);
2224 //cout<<"get 1"<<endl;
2225
2226 double ddl = -1*hot->getDriftDist();
2227 double ddr = hot->getDriftDist();
2228 //double erddl = hot->hitSigma();
2229 //double erddr = hot->hitSigma();
2230 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2231 int lr = 0;
2232 int stat = hot->getUse();
2233 stat = 1; // llwang
2234 //cout<<"hit stat = "<<stat<<endl;
2235 //cout<<"get 2"<<endl;
2236 Identifier mdcId = hot->getDigi()->identify();
2237 //cout<<"get 2.1"<<endl;
2238 double tdc = hot->getDigi()->getTimeChannel();
2239 //cout<<"get 2.2"<<endl;
2240 double adc = hot->getDigi()->getChargeChannel();
2241 //cout<<"get 2.3"<<endl;
2242 double driftT = hot->driftTime();
2243 //cout<<"get 2.4"<<endl;
2244 double doca = hot->getDriftDist() +hot->getResidual();
2245 //cout<<"get 2.5"<<endl;
2246 double entra = direction.phi()-position.phi();//FIXME
2247 //cout<<"get 2.6"<<endl;
2248 while(entra<-0.5*M_PI)entra+= M_PI;
2249 while(entra> 0.5*M_PI)entra -= M_PI;
2250 //cout<<"get 3"<<endl;
2251 double zhit = hitPoint.z();
2252 double fltLen = houghTrack->flightLength(hitPoint);
2253 //cout<<"after get "<<hitId<<endl;
2254
2255 recMdcHit->setDriftDistLeft(ddl);
2256 recMdcHit->setDriftDistRight(ddr);
2257 //recMdcHit->setErrDriftDistLeft(erddl);
2258 //recMdcHit->setErrDriftDistRight(erddr);
2259 recMdcHit->setChisqAdd(pChisq);
2260 recMdcHit->setFlagLR(lr);
2261 recMdcHit->setStat(stat);
2262 recMdcHit->setMdcId(mdcId);
2263 recMdcHit->setTdc(tdc);
2264 recMdcHit->setAdc(adc);
2265 recMdcHit->setDriftT(driftT);
2266 recMdcHit->setDoca(doca);
2267 recMdcHit->setEntra(entra);
2268 recMdcHit->setZhit(zhit);
2269 recMdcHit->setFltLen(fltLen);
2270 //cout<<"after set "<<hitId<<endl;
2271
2272 hitList->push_back(recMdcHit);
2273 SmartRef<RecMdcHit> refHit(recMdcHit);
2274 hitRefVec.push_back(refHit);
2275
2276 if(hot->getLayer()>maxLayerId){
2277 maxLayerId = hot->getLayer();
2278 fiTermHot = *hot;
2279 }
2280 if(hot->getLayer()<minLayerId){
2281 minLayerId = hot->getLayer();
2282 }
2283 if(hot->getUse()==1){
2284 recMdcHit->setStat(1);
2285 }else{
2286 recMdcHit->setStat(0);
2287 }
2288 if(hot->getFlag()!=0)++nSt;
2289 //cout<<"finish a DC hit "<<hitId<<endl;
2290 }else if(hot->getHitType()==HoughHit::CGEMHIT){
2291 //hot->print();
2292 //cout<<endl;
2293 const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
2294 int clusterid = hot->getCgemCluster()->getclusterid();
2295 //int stat = hot->getUse();
2296 int stat = 1;
2297 clusterFitStat[clusterid] = stat;
2298 clusterRefVec.push_back(recCgemCluster);
2299 }
2300 hitId++;
2301 }
2302
2303 //HepPoint3D point = fiTermHot.getHitPosition();
2304 //fiTerm = houghTrack->flightArc(point)/houghTrack->radius();
2305 //recMdcTrack->setFiTerm(fiTerm);
2306 recMdcTrack->setNdof(nHits-5);
2307 recMdcTrack->setNster(nSt);
2308 recMdcTrack->setFirstLayer(minLayerId);
2309 recMdcTrack->setLastLayer(maxLayerId);
2310 //cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2311 recMdcTrack->setVecHits(hitRefVec);
2312 //cout<<"before setVecClusters"<<endl;
2313 //recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2314 //recMdcTrack->setVecClusters(clusterRefVec);
2315 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2316 //cout<<"clusterRefVec.size = "<<clusterRefVec.size()<<endl;
2317 //cout<<"end setVecClusters"<<endl;
2318 recMdcTrack->setNcluster(clusterRefVec.size());
2319 //cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2320 trackList->push_back(recMdcTrack);
2321 //trackId++;
2322 //cout<<helixPar[2]<<endl;
2323 }
2324 return trackList->size();
2325//*/
2326 return trackList->size();
2327}
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
****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
void setError(double err[15])
void setHelix(double helix[5])
void setPoca(double poca[3])
int wireAmbig() const
double dcaToWire() const
double entranceAngle() const
const MdcLayer * layer() const
const MdcDigi * digi() const
double driftTime(double tof, double z) const
Definition: MdcHit.cxx:142
const MdcHit * mdcHit() const
void setVecClusters(ClusterRefVec vecclusters)
Definition: RecMdcTrack.cxx:86
void setVecHits(HitRefVec vechits)
Definition: RecMdcTrack.cxx:80
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
virtual TrkExchangePar helix(double fltL) const =0
double resid(bool exclude=false) const
TrkErrCode getFitStuff(HepVector &derivs, double &deltaChi) const
virtual double arrivalTime(double fltL) const
Definition: TrkRep.cxx:192

Referenced by execute().

◆ storeTracks() [2/2]

int HoughFinder::storeTracks ( RecMdcTrackCol trackList,
RecMdcHitCol hitList,
vector< HoughTrack > &  trackVector 
)

Member Data Documentation

◆ Layer

int HoughFinder::Layer

Definition at line 60 of file InstallArea/include/HoughTransAlg/HoughTransAlg/Hough.h.

Referenced by HoughFinder().

◆ m_clearTrack

int HoughFinder::m_clearTrack

◆ m_useCgemInGlobalFit

int HoughFinder::m_useCgemInGlobalFit

◆ printFlag

int HoughFinder::printFlag

◆ storeFlag

int HoughFinder::storeFlag

Definition at line 63 of file InstallArea/include/HoughTransAlg/HoughTransAlg/Hough.h.

Referenced by execute(), and HoughFinder().


The documentation for this class was generated from the following files: