CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughTrack Class Reference

#include <HoughTrack.h>

+ Inheritance diagram for HoughTrack:

Public Member Functions

 HoughTrack (int charge, double angle, double rho, double dAngle, double dRho, int trkID)
 
 HoughTrack (int charge, const HepPoint3D &position, const Hep3Vector &momentum, int trkID)
 
 HoughTrack (HepPoint3D &pivot, HepVector &a, int trkID)
 
 HoughTrack (const HoughTrack &other)
 

 
HoughTrackoperator= (const HoughTrack &other)
 
 HoughTrack ()
 
void setTrkID (int trkID)
 
void setFlag (int flag)
 
void setCharge (int charge)
 
void setAngle (double angle)
 
void setRho (double rho)
 
void setDAngle (double dAngle)
 
void setDRho (double dRho)
 
void setDTanl (double dTanl)
 
void setDDz (double dDz)
 
void setChi2 (double chi2)
 
void setDz (double dz)
 
void setTanl (double tanl)
 
int getTrkID () const
 
int getCharge () const
 
int getFlag () const
 
double getAngle () const
 
double getRho () const
 
double getDAngle () const
 
double getDRho () const
 
double getDTanl () const
 
double getDDz () const
 
double getChi2 () const
 
int getCircleFitStat () const
 
TrkRecoTrkgetTrkRecoTrk ()
 
double getDr () const
 
double getPhi0 () const
 
double getKappa () const
 
double getDz () const
 
double getTanl () const
 
int findXHot (vector< HoughHit * > &hitList, int charge)
 
void sortHot (vector< HoughHit * > &hotList)
 
vector< HoughHit * > getVecHitPnt ()
 
void addHitPnt (HoughHit *aHitPnt)
 
void dropHitPnt (HoughHit *aHitPnt)
 
void dropVHitPnt (HoughHit *aHitPnt)
 
int getNHitsShared ()
 
void dropRedundentCgemXHits ()
 
void dropRedundentCgemVHits ()
 
int getNhitFirstHalf ()
 
int getNhitSecondHalf ()
 
int getNhitUnusedFirstHalf ()
 
int getNhitUnusedSecondHalf ()
 
void resetNhitHalf ()
 
double driftDistRes (HoughHit *hit)
 
int judgeHalf (HoughHit *hit)
 
int judgeCharge (HoughHit *hit)
 
int judgeCharge (double xHit, double yHit)
 
TrkErrCode fitCircle (const MdcDetector *mdcDetector, TrkContextEv *trkContextEv, double bunchT0)
 
int fitCircle (DotsHelixFitter *fitter, double bunchT0)
 
int fitHelix (DotsHelixFitter *fitter, double bunchT0, RecCgemClusterCol::iterator recCgemClusterColBegin, double averageChi2cut=25)
 
int calculateZ_S (HoughHit *hit)
 
void updateHelix ()
 
void update (double angle, double rho)
 
void updateCirclePar (double dr, double phi0, double kappa)
 
void markUsedHot (vector< HoughHit * > &hitPntList, int use=1)
 
void markUsedHot (int use=1)
 
void print ()
 
void printHot ()
 
void clearHits ()
 
bool isAGoodCircle ()
 
int getNTried ()
 
void releaseSelHits ()
 
vector< double > getVecHitRes ()
 
vector< HoughHit * > getVecStereoHitPnt ()
 
void addVecStereoHitPnt (HoughHit *aHitPnt)
 
vector< double > getVecStereoHitRes ()
 
TrkErrCode fitHelix (const MdcDetector *mdcDetector, const BField *bField, double bunchT0, vector< HoughHit * > &hot, int Layer)
 
TrkErrCode fitHelix (const MdcDetector *mdcDetector, TrkContextEv *trkContextEv, double bunchT0, vector< MdcHit * > &mdcHitCol, vector< HoughHit * > &hot)
 
int findVHot (vector< HoughHit * > &hitList, int charge, int maxLayer, double cutFactor=1.0)
 
vector< HoughHit * > getHotList (int type=2)
 
HoughTrackgetMcTrack () const
 
void setMcTrack (HoughTrack *mcTrack)
 
void clearMemory ()
 
void setRecMdcHitVec (vector< RecMdcHit > &aRecMdcHitVec)
 
vector< RecMdcHit > * getRecMdcHitVec ()
 
void printRecMdcHitVec ()
 
- Public Member Functions inherited from Helix
 Helix ()
 
 Helix (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
 Constructor with pivot, helix parameter a, and its error matrix.
 
 Helix (const HepPoint3D &pivot, const HepVector &a)
 Constructor without error matrix.
 
 Helix (const HepPoint3D &position, const Hep3Vector &momentum, double charge)
 Constructor with position, momentum, and charge.
 
 Helix (const Helix &i)
 
virtual ~Helix ()
 Destructor.
 
const HepPoint3Dcenter (void) const
 returns position of helix center(z = 0.);
 
const HepPoint3Dpivot (void) const
 returns pivot position.
 
double radius (void) const
 returns radious of helix.
 
HepPoint3D x (double dPhi=0.) const
 returns position after rotating angle dPhi in phi direction.
 
double * x (double dPhi, double p[3]) const
 
HepPoint3D x (double dPhi, HepSymMatrix &Ex) const
 returns position and convariance matrix(Ex) after rotation.
 
Hep3Vector direction (double dPhi=0.) const
 returns direction vector after rotating angle dPhi in phi direction.
 
Hep3Vector momentum (double dPhi=0.) const
 returns momentum vector after rotating angle dPhi in phi direction.
 
Hep3Vector momentum (double dPhi, HepSymMatrix &Em) const
 returns momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass, HepSymMatrix &Em) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
HepLorentzVector momentum (double dPhi, double mass, HepPoint3D &x, HepSymMatrix &Emx) const
 returns 4momentum vector after rotating angle dPhi in phi direction.
 
double dr (void) const
 returns an element of parameters.
 
double phi0 (void) const
 
double kappa (void) const
 
double dz (void) const
 
double tanl (void) const
 
double curv (void) const
 
double sinPhi0 (void) const
 
double cosPhi0 (void) const
 
double alpha (void) const
 
const HepVector & a (void) const
 returns helix parameters.
 
const HepSymMatrix & Ea (void) const
 returns error matrix.
 
double pt (void) const
 
double cosTheta (void) const
 
const HepVector & a (const HepVector &newA)
 sets helix parameters.
 
const HepSymMatrix & Ea (const HepSymMatrix &newdA)
 sets helix paramters and error matrix.
 
const HepPoint3Dpivot (const HepPoint3D &newPivot)
 sets pivot position.
 
void set (const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
 sets helix pivot position, parameters, and error matrix.
 
void ignoreErrorMatrix (void)
 unsets error matrix. Error calculations will be ignored after this function call until an error matrix be set again. 0 matrix will be return as a return value for error matrix when you call functions which returns an error matrix.
 
double bFieldZ (double)
 sets/returns z componet of the magnetic field.
 
double bFieldZ (void) const
 
Helixoperator= (const Helix &)
 Copy operator.
 
HepMatrix delApDelA (const HepVector &ap) const
 
HepMatrix delXDelA (double phi) const
 
HepMatrix delMDelA (double phi) const
 
HepMatrix del4MDelA (double phi, double mass) const
 
HepMatrix del4MXDelA (double phi, double mass) const
 
double IntersectCylinder (double r) const
 
double flightArc (HepPoint3D &hit) const
 
double flightArc (double r) const
 
double flightLength (HepPoint3D &hit) const
 
double dPhi (HepPoint3D &hit) const
 

Static Public Attributes

static int m_clearTrack =1
 
static int m_useCgemInGlobalFit =3
 
static TGraph * m_cut [2][43] = {NULL}
 
- Static Public Attributes inherited from Helix
static const double ConstantAlpha = 333.564095
 Constant alpha for uniform field.
 

Additional Inherited Members

- Protected Attributes inherited from Helix
IMagneticFieldSvcm_pmgnIMF
 
double m_bField
 
double m_alpha
 

Detailed Description

Definition at line 22 of file HoughTrack.h.

Constructor & Destructor Documentation

◆ HoughTrack() [1/5]

HoughTrack::HoughTrack ( int  charge,
double  angle,
double  rho,
double  dAngle,
double  dRho,
int  trkID 
)

Definition at line 20 of file HoughTrack.cxx.

20 :Helix()
21{
22 bFieldZ(-bFieldZ());
23
24 m_trkID = trkID;
25 m_charge = (charge>0)?1:-1;
26 m_angle = angle;
27 m_rho = rho;
28 m_flag = 0;
29 m_dAngle = dAngle;
30 m_dRho = dRho;
31 m_dTanl = 0;
32 m_dDz = 0;
33 m_chi2 = 0;
34 m_circleFitStat = 0;
35
36 m_dr = 0;
37 m_phi0 = (rho>0)? angle:angle+M_PI;
38 m_phi0 = (m_charge<0)? m_phi0:m_phi0+M_PI;
39 if(m_phi0>2*M_PI)m_phi0-=2*M_PI;
40 if(m_phi0<0) m_phi0+=2*M_PI;
41 m_kappa = fabs(m_alpha*rho)*m_charge;
42 m_dz = 0.0;
43 m_tanl = 0.0;
44
45 HepPoint3D pivot(0,0,0);
46 HepVector a(5,0);
47 HepSymMatrix Ea(5,0);
48 a(2) = m_phi0;
49 a(3) = m_kappa;
50 set(pivot,a,Ea);
51
52 //Helix helix(pivot,a); helix.bFieldZ(-helix.bFieldZ());
53 //KalmanFit::Helix helix(pivot,a);
54 //m_helix = helix;
55 //m_hot.clear();
56 m_trkRecoTrk = NULL;
57
58 m_vecHitPnt.clear();
59 m_vecHitResidual.clear();
60 m_vecHitChi2.clear();
61
62 //m_dr = m_helix.dr();
63 //m_phi0 = m_helix.phi0();
64 //m_kappa = m_helix.kappa();
65 //m_dz = m_helix.dz();
66 //m_tanl = m_helix.tanl();
67
68 m_unused=0;
69 m_untried=0;
70 m_mcTrack = NULL;
71 m_cgemHitVector.clear();
72}
#define M_PI
Definition: TConstant.h:4
void set(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
sets helix pivot position, parameters, and error matrix.
const HepSymMatrix & Ea(void) const
returns error matrix.
const HepVector & a(void) const
returns helix parameters.
const HepPoint3D & pivot(void) const
returns pivot position.

◆ HoughTrack() [2/5]

HoughTrack::HoughTrack ( int  charge,
const HepPoint3D position,
const Hep3Vector &  momentum,
int  trkID 
)

Definition at line 74 of file HoughTrack.cxx.

74 :Helix(position,momentum,charge)
75{
76 bFieldZ(-bFieldZ());
77 //Helix helix(position, momentum, charge); helix.bFieldZ(-helix.bFieldZ());
78 //m_helix = helix;
79 m_trkID = trkID;
80 m_charge = charge;
81 //m_angle = m_helix.center().phi();
82 //m_rho = 1./m_helix.radius();
83 m_angle = center().phi();
84 m_rho = 1./radius();
85 m_flag = 0;
86 if(m_angle<0)m_angle=+2*M_PI;
87 m_circleFitStat = 0;
88
89 //m_dr = m_helix.dr();
90 //m_phi0 = m_helix.phi0();
91 //m_kappa = m_helix.kappa();
92 //m_dz = m_helix.dz();
93 //m_tanl = m_helix.tanl();
94 m_dr = dr();
95 m_phi0 = phi0();
96 m_kappa = kappa();
97 m_dz = dz();
98 m_tanl = tanl();
99
100 //m_hot.clear();
101 m_trkRecoTrk = NULL;
102
103 m_vecHitPnt.clear();
104 m_vecHitResidual.clear();
105 m_vecHitChi2.clear();
106
107 m_dAngle = 0;
108 m_dRho = 0;
109 m_dTanl = 0;
110 m_dDz = 0;
111 m_chi2 = 0;
112 m_mcTrack = NULL;
113 m_cgemHitVector.clear();
114}
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
double radius(void) const
returns radious of helix.
double dr(void) const
returns an element of parameters.

◆ HoughTrack() [3/5]

HoughTrack::HoughTrack ( HepPoint3D pivot,
HepVector &  a,
int  trkID 
)

Definition at line 116 of file HoughTrack.cxx.

116 :Helix(pivot,a)
117{
118 bFieldZ(-bFieldZ());
119 //Helix helix(pivot,a); helix.bFieldZ(-helix.bFieldZ());
120 //m_helix = helix;
121 m_trkID = trkID;
122 //m_charge = (m_helix.kappa()>0)? 1:-1;;
123 //m_angle = m_helix.center().phi();
124 //m_rho = 1./m_helix.radius();
125 m_charge = (kappa()>0)? 1:-1;;
126 m_angle = center().phi();
127 m_rho = 1./radius();
128 m_flag = 0;
129 m_circleFitStat = 0;
130 if(m_angle<0)m_angle=+2*M_PI;
131
132 //m_dr = m_helix.dr();
133 //m_phi0 = m_helix.phi0();
134 //m_kappa = m_helix.kappa();
135 //m_dz = m_helix.dz();
136 //m_tanl = m_helix.tanl();
137 m_dr = dr();
138 m_phi0 = phi0();
139 m_kappa = kappa();
140 m_dz = dz();
141 m_tanl = tanl();
142
143 //m_hot.clear();
144 m_trkRecoTrk = NULL;
145
146 m_vecHitPnt.clear();
147 m_vecHitResidual.clear();
148 m_vecHitChi2.clear();
149
150 m_dAngle = 0;
151 m_dRho = 0;
152 m_dTanl = 0;
153 m_dDz = 0;
154 m_chi2 = 0;
155 m_mcTrack = NULL;
156 m_cgemHitVector.clear();
157}

◆ HoughTrack() [4/5]

HoughTrack::HoughTrack ( const HoughTrack other)

Definition at line 222 of file HoughTrack.cxx.

222 :
223 Helix(other),
224 m_trkID(other.m_trkID),
225 m_charge(other.m_charge),
226 m_angle(other.m_angle),
227 m_rho(other.m_rho),
228 m_flag(other.m_flag),
229
230 m_dAngle(other.m_dAngle),
231 m_dRho(other.m_dRho),
232 m_dTanl(other.m_dTanl),
233 m_dDz(other.m_dTanl),
234
235 m_dr(other.m_dr),
236 m_phi0(other.m_phi0),
237 m_kappa(other.m_kappa),
238 m_dz(other.m_dz),
239 m_tanl(other.m_tanl),
240
241 m_chi2(other.m_chi2),
242 m_trkRecoTrk(other.m_trkRecoTrk),
243
244 m_circleFitStat(other.m_circleFitStat),
245 m_vecHitPnt(other.m_vecHitPnt),
246 m_vecHitResidual(other.m_vecHitResidual),
247 m_vecHitChi2(other.m_vecHitChi2),
248 m_vecStereoHitPnt(other.m_vecStereoHitPnt),
249 m_vecStereoHitRes(other.m_vecStereoHitRes),
250 m_mcTrack(other.m_mcTrack),
251 m_cgemHitVector(other.m_cgemHitVector),
252 m_vecRecMdcHit(other.m_vecRecMdcHit)
253{}
Index other(Index i, Index j)
Definition: EvtCyclic3.cc:118

◆ HoughTrack() [5/5]

HoughTrack::HoughTrack ( )

Definition at line 187 of file HoughTrack.cxx.

188{
189 m_trkID = -1;
190 m_charge = 0;
191 m_angle = 0;
192 m_rho = 1e-6;
193 m_flag = 0;
194 HepVector a(5,0);
195 HepPoint3D pivot(0,0,0);
196 bFieldZ(-bFieldZ());
197 m_circleFitStat = 0;
198
199 m_dr = dr();
200 m_phi0 = phi0();
201 m_kappa = kappa();
202 m_dz = dz();
203 m_tanl = tanl();
204
205 //m_hot.clear();
206 m_trkRecoTrk = NULL;
207
208 m_vecHitPnt.clear();
209 m_vecHitResidual.clear();
210 m_vecHitChi2.clear();
211
212 m_dAngle = 0;
213 m_dRho = 0;
214 m_dTanl = 0;
215 m_dDz = 0;
216 m_chi2 = 0;
217 m_mcTrack = NULL;
218 m_cgemHitVector.clear();
219}

Member Function Documentation

◆ addHitPnt()

void HoughTrack::addHitPnt ( HoughHit aHitPnt)
inline

Definition at line 84 of file HoughTrack.h.

84{m_vecHitPnt.push_back(aHitPnt);};// llwang 2018-12-26

Referenced by HoughFinder::getMcParticleCol().

◆ addVecStereoHitPnt()

void HoughTrack::addVecStereoHitPnt ( HoughHit aHitPnt)
inline

Definition at line 135 of file HoughTrack.h.

135{m_vecStereoHitPnt.push_back(aHitPnt);};

Referenced by HoughFinder::getMcParticleCol().

◆ calculateZ_S()

int HoughTrack::calculateZ_S ( HoughHit hit)

Definition at line 774 of file HoughTrack.cxx.

775{
776 int nTangency(0);
777 if(hit->getFlag()==0)return nTangency;
778 //if(hit->getPairHit()==NULL)return nTangency;
779 //if(hit->getPairHit()->getHalfCircle()!=1)return nTangency;
780 double s(0),z(0);
781 double xc = center().x();
782 double yc = center().y();
783 double rTrack = radius();// signed FIXME
784 if(hit->getHitType()==HoughHit::CGEMHIT){
785 if(hit->getFlag()!=2)return nTangency;
786 //for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++)
787 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++){
788 if((*hotIt)->getFlag()!=0)continue;
789 if((*hotIt)->getHitType()==HoughHit::MDCHIT)continue;
790 if((*hotIt)->getCgemCluster()->getclusterid() == hit->getCgemCluster()->getclusterflagb()){
791 HepPoint3D point = hit->getHitPosition();
792 s = flightArc(point);
793 z = point.z();
794 // double rHit = point.perp();
795 // s = flightArc(rHit);
796 // z = getZfrom ... // FIXME
797 pair<double,double> sz(s,z);
798 hit->addSZ(sz);
799 nTangency++;
800 //double s = flightArc(hit->getHitPosition());
801 }
802 }
803 }else{
804 hit->resetSZ();
805 double drift = hit->getDriftDist();
806 const MdcGeoWire* wire = hit->getMdcGeomSvc()->Wire(hit->getLayer(),hit->getWire());
807 HepPoint3D westPoint = wire->Forward();
808 HepPoint3D eastPoint = wire->Backward();
809 double xEast = eastPoint.x()/10.0;
810 double xWest = westPoint.x()/10.0;
811 double yEast = eastPoint.y()/10.0;
812 double yWest = westPoint.y()/10.0;
813 double zEast = eastPoint.z()/10.0;
814 double zWest = westPoint.z()/10.0;
815 //cout<<"wast:x,y,z: "<<xWest<<", "<<yWest<<", "<<zWest<<endl;
816 //cout<<"east:x,y,z: "<<xEast<<", "<<yEast<<", "<<zEast<<endl;
817 //cout<<"Xc, Yc, Rc: "<<xc<<", "<<yc<<", "<<rTrack<<endl;
818 double west2east = sqrt((xEast-xWest)*(xEast-xWest)+(yEast-yWest)*(yEast-yWest));
819
820 double slope = (yEast-yWest)/(xEast-xWest);
821 double intercept = (yWest-slope*xWest+yEast-slope*xEast)/2;
822 double a = 1+slope*slope;
823 double b = -2*(xc+slope*yc-slope*intercept);
824 double c1 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+drift)*(rTrack+drift);
825 double c2 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack-drift)*(rTrack-drift);
826 //double c1 = intercept*(2*yc-intercept)+drift*(drift+rTrack);
827 //double c2 = intercept*(2*yc-intercept)+drift*(drift-rTrack);
828 double delta1 = (b*b-4*a*c1);
829 double delta2 = (b*b-4*a*c2);
830 //cout<<"a,b,c1,c2: "<<a<<", "<<b<<", "<<c1<<", "<<c2<<endl;
831 //cout<<"delta: "<<delta1<<", "<<delta2<<endl;
832 if(delta1>=0){
833 double x1 = (-b+sqrt(delta1))/(2*a);
834 double x2 = (-b-sqrt(delta1))/(2*a);
835 double y1 = slope*x1+intercept;
836 double y2 = slope*x2+intercept;
837 if((x1-xWest)*(x1-xEast)<0){
838 HepPoint3D point(x1,y1,0);
839 s = flightArc(point);
840 //s = flightArc(hit->getHitPosition());
841 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
842 z = zWest + l/west2east*fabs((zEast-zWest));
843 pair<double,double> sz(s,z);
844 hit->addSZ(sz);
845 nTangency++;
846 HepPoint3D position(x1,y1,z);
847 //hit->addPosition(position);
848 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
849 }
850 if((x2-xWest)*(x2-xEast)<0){
851 HepPoint3D point(x2,y2,0);
852 s = flightArc(point);
853 //s = flightArc(hit->getHitPosition());
854 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
855 z = zWest + l/west2east*fabs((zEast-zWest));
856 pair<double,double> sz(s,z);
857 hit->addSZ(sz);
858 nTangency++;
859 HepPoint3D position(x2,y2,z);
860 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
861 //hit->addPosition(position);
862 }
863 }
864
865 if( delta2>=0){
866 double x1 = (-b+sqrt(delta2))/(2*a);
867 double x2 = (-b-sqrt(delta2))/(2*a);
868 double y1 = slope*x1+intercept;
869 double y2 = slope*x2+intercept;
870 if((x1-xWest)*(x1-xEast)<0){
871 HepPoint3D point(x1,y1,0);
872 s = flightArc(point);
873 //s = flightArc(hit->getHitPosition());
874 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
875 z = zWest + l/west2east*fabs((zEast-zWest));
876 pair<double,double> sz(s,z);
877 hit->addSZ(sz);
878 nTangency++;
879 HepPoint3D position(x1,y1,z);
880 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
881 //hit->addPosition(position);
882 }
883 if((x2-xWest)*(x2-xEast)<0){
884 HepPoint3D point(x2,y2,0);
885 s = flightArc(point);
886 //s = flightArc(hit->getHitPosition());
887 //double l = zWest + sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
888 //z = l/west2east*fabs((zEast-zWest));
889 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
890 z = zWest + l/west2east*fabs((zEast-zWest));
891 pair<double,double> sz(s,z);
892 hit->addSZ(sz);
893 nTangency++;
894 HepPoint3D position(x2,y2,z);
895 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
896 //hit->addPosition(position);
897 }
898 }
899 }
900 return nTangency;
901}
XmlRpcServer s
Definition: HelloServer.cpp:11
MdcGeomSvc * getMdcGeomSvc() const
Definition: HoughHit.h:65
const RecCgemCluster * getCgemCluster() const
Definition: HoughHit.h:41
double getDriftDist() const
Definition: HoughHit.h:54
void addSZ(S_Z sz)
Definition: HoughHit.h:106
HepPoint3D getHitPosition() const
Definition: HoughHit.h:51
HitType getHitType() const
Definition: HoughHit.h:40
int getFlag() const
Definition: HoughHit.h:47
int getLayer() const
Definition: HoughHit.h:45
int getWire() const
Definition: HoughHit.h:46
void resetSZ()
Definition: HoughHit.h:97
@ CGEMHIT
Definition: HoughHit.h:27
@ MDCHIT
Definition: HoughHit.h:27
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
int getclusterflagb(void) const
TCanvas * c1
Definition: tau_mode.c:75

Referenced by HoughFinder::fillHistogram(), and findVHot().

◆ clearHits()

void HoughTrack::clearHits ( )

Definition at line 1193 of file HoughTrack.cxx.

1194{
1195 //m_hot.clear();
1196 m_vecHitPnt.clear();
1197 m_vecHitResidual.clear();
1198 m_vecHitChi2.clear();
1199 m_map_lay_d.clear();
1200 m_map_lay_hit.clear();
1201 m_setLayer.clear();
1202 m_unused=0;
1203 m_untried=0;
1204 m_nHitFirstHalf=0;
1205 m_nHitSecondHalf=0;
1206 m_nHitUnused_FirstHalf=0;
1207 m_nHitUnused_SecondHalf=0;
1208
1209
1210}

◆ clearMemory()

void HoughTrack::clearMemory ( )

Definition at line 2814 of file HoughTrack.cxx.

2815{
2816 for(vector<CgemHitOnTrack*>::iterator iter = m_cgemHitVector.begin(); iter != m_cgemHitVector.end(); iter++){
2817 delete *iter;
2818 }
2819 if(m_clearTrack)delete m_trkRecoTrk;
2820}
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
static int m_clearTrack
Definition: HoughTrack.h:145

◆ driftDistRes()

double HoughTrack::driftDistRes ( HoughHit hit)

Definition at line 636 of file HoughTrack.cxx.

637{
638 double residual(9999.);
639 HepPoint3D hitPoint = hit->getHitPosition();
640 if(hit->getHitType() == HoughHit::MDCHIT || hit->getHitType() == HoughHit::MDCMCHIT){
641 //HepPoint3D circleCenter = m_helix.center();
642 HepPoint3D circleCenter = center();
643 //double distance = circleCenter.perp(hitPoint);
644 double distance = (circleCenter-hitPoint).perp();
645 //hit->print();
646 //cout<<distance<<endl
647 //<<sqrt((circleCenter-hitPoint).x()*(circleCenter-hitPoint).x()+(circleCenter-hitPoint).y()*(circleCenter-hitPoint).y())<<endl
648 //<<sqrt((circleCenter.x()-hitPoint.x())*(circleCenter.x()-hitPoint.x())+(circleCenter.y()-hitPoint.y())*(circleCenter.y()-hitPoint.y()))<<endl;
649 //double Rc = m_helix.radius();
650 double Rc = fabs(radius());
651 double driftDist(0);
652 //if(hit->getHitType() == HoughHit::MDCMCHIT) driftDist = 0;
653 if(hit->getHitType() == HoughHit::MDCHIT) driftDist = hit->getDriftDist();
654 //residual = fabs(distance - Rc) - driftDist;
655 residual = driftDist-fabs(distance - Rc);
656 //cout<<"driftDist, distance, Rc = "<<driftDist<<", "<<distance<<", "<<Rc<<endl;
657 //double distance = (m_helix.center()).perp(hit->getHitPosition());
658 //res = fabs(distance - m_helix.radius()) - hit->getDriftDist();
659 }else if(hit->getHitType() == HoughHit::CGEMHIT || hit->getHitType() == HoughHit::CGEMMCHIT){
660 double rCgem = hitPoint.perp();
661 //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
662 double phiTrkFlt = IntersectCylinder(rCgem);
663 phiTrkFlt=judgeHalf(hit)*phiTrkFlt;
664 //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
665 HepPoint3D crossPoint = x(phiTrkFlt);
666 double phiCrossPoint = crossPoint.phi();
667 double phiMeasure = hitPoint.phi();
668 residual = phiMeasure - phiCrossPoint;
669 while(residual<-M_PI)residual += 2*M_PI;
670 while(residual> M_PI)residual -= 2*M_PI;
671 residual = rCgem*residual;
672 }
673 hit->setResidual(residual);
674 return residual;
675}
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
@ CGEMMCHIT
Definition: HoughHit.h:27
@ MDCMCHIT
Definition: HoughHit.h:27
void setResidual(double residual)
Definition: HoughHit.h:90
int judgeHalf(HoughHit *hit)
Definition: HoughTrack.cxx:596

Referenced by findXHot().

◆ dropHitPnt()

void HoughTrack::dropHitPnt ( HoughHit aHitPnt)

Definition at line 1861 of file HoughTrack.cxx.

1862{
1863 vector<HoughHit*>::iterator it0 = m_vecHitPnt.begin();
1864 vector<HoughHit*>::iterator result = find(m_vecHitPnt.begin(), m_vecHitPnt.end(), aHitPnt);
1865 if(result!=m_vecHitPnt.end()) {
1866 //cout<<"remove hit "<<(*result)->getHitID();
1867 m_vecHitPnt.erase(result);
1868 vector<double>::iterator itRes = m_vecHitResidual.begin()+(result-it0);
1869 if(itRes!=m_vecHitResidual.end()) {
1870 //cout<<" with residual "<<*itRes<<" from trk "<<getTrkID()<<endl;
1871 m_vecHitResidual.erase(itRes);
1872 }
1873 }
1874}

Referenced by dropRedundentCgemXHits().

◆ dropRedundentCgemVHits()

void HoughTrack::dropRedundentCgemVHits ( )

Definition at line 1945 of file HoughTrack.cxx.

1946{
1947 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
1948 vector<HoughHit*> hitPnt_cgem[3];
1949 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
1950 double res_cgem_min[3]={999., 999., 999.};
1951 vector<HoughHit*>::iterator iter = m_vecStereoHitPnt.begin();
1952 for(; iter!=m_vecStereoHitPnt.end(); iter++)
1953 {
1954 int iLayer = (*iter)->getLayer();
1955 if(iLayer<3) {
1956 hitPnt_cgem[iLayer].push_back(*iter);
1957 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1958 double res = (*iter)->residual(this, positionOntrack, positionOnDetector);
1959 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
1960 res_cgem_min[iLayer] = res;
1961 hitPnt_cgem_resMin[iLayer] = *iter;
1962 }
1963 }
1964 }
1965 for(int iLayer=0; iLayer<3; iLayer++)
1966 {
1967 int nHits = hitPnt_cgem[iLayer].size();
1968 if(nHits>1)
1969 {
1970 HoughHit* hitToKeep = hitPnt_cgem_resMin[iLayer];
1971 for(int iHit=0; iHit<nHits; iHit++)
1972 {
1973 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
1974 if(aHoughHit==hitToKeep) continue;
1975 dropVHitPnt(aHoughHit);
1976 }
1977 }
1978 }
1979 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
1980}
void dropVHitPnt(HoughHit *aHitPnt)

◆ dropRedundentCgemXHits()

void HoughTrack::dropRedundentCgemXHits ( )

Definition at line 1886 of file HoughTrack.cxx.

1887{
1888 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
1889 vector<HoughHit*> hitPnt_cgem[3];
1890 HoughHit* hitPnt_cgem_keep[3]={NULL, NULL, NULL};
1891 double res_cgem_keep[3]={999., 999., 999.};
1892 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
1893 double res_cgem_min[3]={999., 999., 999.};
1894 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1895 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1896 vector<double>::iterator itRes = m_vecHitResidual.begin();
1897 for(; iter!=m_vecHitPnt.end(); iter++)
1898 {
1899 int iLayer = (*iter)->getLayer();
1900 if(iLayer<3) {
1901 hitPnt_cgem[iLayer].push_back(*iter);
1902 int idx = iter-iter0;
1903 double res = m_vecHitResidual[idx];
1904 vector<double> vecRes = (*iter)->getVecResid();
1905 int nTrk = vecRes.size();
1906 if(nTrk==1&&fabs(res)<fabs(res_cgem_keep[iLayer])) {
1907 res_cgem_keep[iLayer] = res;
1908 hitPnt_cgem_keep[iLayer] = *iter;
1909 }
1910 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
1911 res_cgem_min[iLayer] = res;
1912 hitPnt_cgem_resMin[iLayer] = *iter;
1913 }
1914 }
1915 }
1916 for(int iLayer=0; iLayer<3; iLayer++)
1917 {
1918 int nHits = hitPnt_cgem[iLayer].size();
1919 if(nHits>1)
1920 {
1921 HoughHit* hitToKeep;
1922 if(hitPnt_cgem_keep[iLayer]!=NULL) hitToKeep = hitPnt_cgem_keep[iLayer];
1923 else hitToKeep = hitPnt_cgem_resMin[iLayer];
1924 //cout<<"HoughTrack::dropRedundentCgemXHits: "<<endl;
1925 //cout<<" keep hit "<<hitToKeep->getHitID()
1926 //<<", layer "<<hitToKeep->getLayer()
1927 //<<" at ("<<hitToKeep->getHitPosition().x()
1928 //<<", "<<hitToKeep->getHitPosition().y()<<")"<<endl;
1929 for(int iHit=0; iHit<nHits; iHit++)
1930 {
1931 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
1932 if(aHoughHit==hitToKeep) continue;
1933 //cout<<" remove hit "<<aHoughHit->getHitID()
1934 //<<", layer "<<aHoughHit->getLayer()
1935 //<<" at ("<<aHoughHit->getHitPosition().x()
1936 //<<", "<<aHoughHit->getHitPosition().y()<<")"<<endl;
1937 dropHitPnt(aHoughHit);
1938 //aHoughHit->dropTrkID(m_trkID);
1939 }
1940 }
1941 }
1942 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
1943}
void dropHitPnt(HoughHit *aHitPnt)

◆ dropVHitPnt()

void HoughTrack::dropVHitPnt ( HoughHit aHitPnt)

Definition at line 1876 of file HoughTrack.cxx.

1877{
1878 vector<HoughHit*>::iterator it0 = m_vecStereoHitPnt.begin();
1879 vector<HoughHit*>::iterator result = find(m_vecStereoHitPnt.begin(), m_vecStereoHitPnt.end(), aHitPnt);
1880 if(result!=m_vecStereoHitPnt.end()) {
1881 //cout<<"remove hit "<<(*result)->getHitID();
1882 m_vecStereoHitPnt.erase(result);
1883 }
1884}

Referenced by dropRedundentCgemVHits().

◆ findVHot()

int HoughTrack::findVHot ( vector< HoughHit * > &  hitList,
int  charge,
int  maxLayer,
double  cutFactor = 1.0 
)

Definition at line 1996 of file HoughTrack.cxx.

1997{
1998 bool printFlag(false); //printFlag=true;
1999 int nHot(0);
2000 m_vecStereoHitPnt.clear();
2001 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end(); hitIt++){
2002 /*
2003 if(0 == flag){
2004 if((*hitIt)->getFlag() != 0) continue;
2005 if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
2006 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2007 //if((*hitIt)->getPairHit()->getHalfCircle()!=1)continue;//FIXME
2008 double cut1, cut2;
2009 int iLayer = (*hitIt)->getLayer();
2010 int used = (*hitIt)->getUse();
2011 XhitCutWindow(m_rho, iLayer, charge, cut1, cut2);
2012 double res = driftDistRes(*hitIt);
2013 map<int,double>::iterator it_map;
2014 it_map=m_map_lay_d.find(iLayer);
2015 if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
2016 {
2017 m_map_lay_d[iLayer]=res;
2018 m_map_lay_hit[iLayer]=(*hitIt);
2019 }
2020 if(res>cut1&&res<cut2)
2021 {
2022 //cout<<" selected!";
2023 //if(iLayer>3)
2024 m_vecHitPnt.push_back(*hitIt);
2025 m_vecHitResidual.push_back(res);
2026 //m_hot.push_back((*(*hitIt)));
2027 m_chi2 =+ res*res;
2028 m_setLayer.insert(iLayer);
2029 if(used==0) m_untried++;
2030 if(used==0||used==-1) m_unused++;
2031 if(judgeHalf(*hitIt)>0) {
2032 m_nHitFirstHalf++;
2033 if(used<=0) m_nHitUnused_FirstHalf++;
2034 }
2035 else {
2036 m_nHitSecondHalf++;
2037 if(used<=0) m_nHitUnused_SecondHalf++;
2038 }
2039 nHot++;
2040 }// within window
2041 //cout<<endl;
2042 }// X-view
2043 //*/
2044 //else {
2045 int layer = (*hitIt)->getLayer();
2046 int wire = (*hitIt)->getWire();
2047 if((*hitIt)->getFlag() == 0) continue;
2048 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
2049 if(calculateZ_S(*hitIt)==0){
2050 continue;
2051 }
2052 }
2053 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2054 //if((*hitIt)->getPairHit()->getHalfCircle()>1)continue;//FIXME
2055 double Pt = fabs(pt());
2056 double cut[2] = {0};
2057 if(layer<3){
2058 if(m_cut[0][layer+3]!=NULL)cut[0]=m_cut[0][layer+3]->Eval(Pt);
2059 if(m_cut[1][layer+3]!=NULL)cut[1]=m_cut[1][layer+3]->Eval(Pt);
2060 }else{
2061 if(m_cut[0][layer]!=NULL)cut[0]=m_cut[0][layer]->Eval(Pt);
2062 if(m_cut[1][layer]!=NULL)cut[1]=m_cut[1][layer]->Eval(Pt);
2063 }
2064 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
2065 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
2066 if(judgeCharge(positionOntrack.x(),positionOntrack.y())*charge<0) continue;// charge requirement
2067 bool isNewHot(true);
2068 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++){
2069 if((*hitIt)->getHitID() == (*hotIt)->getHitID()){
2070 isNewHot = false;
2071 break;
2072 }
2073 }
2074 if(!isNewHot){
2075 //cout<<" OK!";
2076 //cout<<endl;
2077 continue;
2078 }
2079 if(printFlag)(*hitIt)->print();
2080 if(printFlag)cout<<"res:"<<setw(12)<<res;
2081 if(printFlag)cout<<" win: "<<setw(5)<<cut[0]<<" ~ "<<setw(5)<<cut[1];
2082 double ext=cutFactor;
2083 if(ext*cut[0]<res&&res<ext*cut[1]){
2084 if((*hitIt)->getLayer()>=maxLayer)
2085 {
2086 //if(printFlag) cout<<;
2087 continue;
2088 }
2089 m_vecStereoHitPnt.push_back(*hitIt);
2090 if(printFlag)cout<<" selected!";
2091 nHot++;
2092 /*
2093 if(maxLayer<20){
2094 if((*hitIt)->getLayer()>=maxLayer)continue;
2095 //HepPoint3D szz(minS,zTrk,zHit);
2096 //(*hitIt)->setHitPosition(szz);
2097 m_vecStereoHitPnt.push_back(*hitIt);
2098 //m_vecStereoHitRes.push_back(minRes);
2099 nHot++;
2100 //cout<<" OK!";
2101 //(*hitIt)->print();
2102 }
2103 else{
2104 //if((*hitIt)->getLayer()!=maxLayer)continue;
2105 if((*hitIt)->getLayer()>=maxLayer)continue;
2106 //HepPoint3D szz(minS,zTrk,zHit);
2107 //(*hitIt)->setHitPosition(szz);
2108 m_vecStereoHitPnt.push_back(*hitIt);
2109 //m_vecStereoHitRes.push_back(minRes);
2110 nHot++;
2111 //cout<<" OK!";
2112 //(*hitIt)->print();
2113 }
2114 //*/
2115 }
2116 if(printFlag)cout<<endl;
2117 //if(layer<3)cout<<endl;
2118 //}// V-view
2119 }// loop hits
2120 return nHot;
2121}
*********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
int judgeCharge(HoughHit *hit)
Definition: HoughTrack.cxx:609
static TGraph * m_cut[2][43]
Definition: HoughTrack.h:147
int calculateZ_S(HoughHit *hit)
Definition: HoughTrack.cxx:774

◆ findXHot()

int HoughTrack::findXHot ( vector< HoughHit * > &  hitList,
int  charge 
)

Definition at line 481 of file HoughTrack.cxx.

482{
483
484 bool printFlag(false); //printFlag=true;
485 int nHot(0);
486 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end(); hitIt++){
487 //cout<<judgeCharge(*(*hitIt))*charge<<" ";
488 //if(0 == flag){
489 if((*hitIt)->getFlag() != 0) continue;
490 if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
491 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
492 //if((*hitIt)->getPairHit()->getHalfCircle()>1)continue;//FIXME
493 if(printFlag)(*hitIt)->print();
494 //double cut = 999.0;//FIXME
495 //double cut1, cut2;
496 double cut[2] = {0};
497 int iLayer = (*hitIt)->getLayer();
498 int used = (*hitIt)->getUse();
499 XhitCutWindow(m_rho, iLayer, charge, cut[0], cut[1]);
500 double res = driftDistRes(*hitIt);
501 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
502 double res2 = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
503 //if(fabs(res)<cut)
504 //cout<<"("<<(*hitIt)->getLayer()<<", "<<(*hitIt)->getWire()<<") ";
505 if(printFlag)cout<<"res:"<<setw(12)<<res;
506 //cout<<setw(12)<<res2;
507 if(printFlag)cout<<" win: "<<setw(5)<<cut[0]<<" ~ "<<setw(5)<<cut[1];
508 //cout<<endl;
509 map<int,double>::iterator it_map;
510 it_map=m_map_lay_d.find(iLayer);
511 if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
512 {
513 m_map_lay_d[iLayer]=res;
514 m_map_lay_hit[iLayer]=(*hitIt);
515 }
516 if(res>cut[0]&&res<cut[1])
517 {
518 if(printFlag)cout<<" selected!";
519 //if(iLayer>3)
520 m_vecHitPnt.push_back((*hitIt));
521 m_vecHitResidual.push_back(res);
522 //m_hot.push_back((*(*hitIt)));
523 //m_chi2 =+ res*res;
524 m_setLayer.insert(iLayer);
525 if(used==0) m_untried++;
526 if(used==0||used==-1) m_unused++;
527 if(judgeHalf(*hitIt)>0)
528 {
529 m_nHitFirstHalf++;
530 if(used<=0) m_nHitUnused_FirstHalf++;
531 }
532 else
533 {
534 m_nHitSecondHalf++;
535 if(used<=0) m_nHitUnused_SecondHalf++;
536 }
537 nHot++;
538 }// within window
539 if(printFlag)cout<<endl;
540 //}// X-view
541 /*
542 else {
543 if((*hitIt)->getFlag() == 0) continue;
544 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
545 //if((*hitIt)->getPairHit()->getHalfCircle()!=1)continue;//FIXME
546 int layer = (*hitIt)->getLayer();
547 int wire = (*hitIt)->getWire();
548 double Pt = fabs(pt());
549 double cut[2] = {0};
550 if(layer<3){
551 if(m_cut[0][layer+3]!=NULL)cut[0]=m_cut[0][layer+3]->Eval(Pt);
552 if(m_cut[1][layer+3]!=NULL)cut[1]=m_cut[1][layer+3]->Eval(Pt);
553 }else{
554 if(m_cut[0][layer]!=NULL)cut[0]=m_cut[0][layer]->Eval(Pt);
555 if(m_cut[1][layer]!=NULL)cut[1]=m_cut[1][layer]->Eval(Pt);
556 }
557 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
558 //cout<<"("<<setw(12)<<cut[0]<<" , "<<setw(13)<<cut[1]<<") ";
559 //cout<<Pt<<endl;
560 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
561 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
562 bool isNewHot(true);
563 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++){
564 if((*hitIt)->getHitID() == (*hotIt)->getHitID()){
565 isNewHot = false;
566 break;
567 }
568 }
569 if(!isNewHot){
570 //cout<<" OK!";
571 //cout<<endl;
572 continue;
573 }
574 if(cut[0]<res&&res<cut[1]){
575 //HepPoint3D szz(minS,zTrk,zHit);
576 //(*hitIt)->setHitPosition(szz);
577 m_vecStereoHitPnt.push_back((*hitIt));
578 //m_vecStereoHitRes.push_back(minRes);
579 //m_hot.push_back(*(*hitIt));
580 nHot++;
581 //cout<<" OK!";
582 //(*hitIt)->print();
583 }
584 }// V-view
585 // */
586 }// loop hits
587
588 m_maxGap=0;
589 m_nGap=XGapSize(m_setLayer,m_maxGap);
590 return nHot;
591}
double driftDistRes(HoughHit *hit)
Definition: HoughTrack.cxx:636

◆ fitCircle() [1/2]

TrkErrCode HoughTrack::fitCircle ( const MdcDetector mdcDetector,
TrkContextEv trkContextEv,
double  bunchT0 
)

Definition at line 2179 of file HoughTrack.cxx.

2180{
2181 bool debug=false; //debug=true;
2182 m_circleFitStat=0;
2183 double d0 = -m_dr;
2184 double phi0 = m_phi0+M_PI/2;
2185 while(phi0>M_PI)phi0-=2*M_PI;
2186 while(phi0<-M_PI)phi0+=2*M_PI;
2187 double omega = m_kappa/fabs(alpha());
2188 //double z0 = m_dz;
2189 //double tanl = m_tanl;
2190 double z0 = 0;
2191 double tanl = 0;
2192 //cout<<"hough params:"<<endl;
2193 //cout
2194 //<<setw(15)<<d0
2195 //<<setw(15)<<phi0
2196 //<<setw(15)<<omega
2197 //;
2198 //cout<<endl;
2199 TrkExchangePar circleTrk(d0,phi0,omega,z0,tanl);
2200 TrkCircleMaker circlefactory;
2201 double chisum =0.;
2202 TrkRecoTrk* trkRecoTrk = circlefactory.makeTrack(circleTrk, chisum, *trkContextEv, bunchT0);
2203 bool permitFlips = true;
2204 bool lPickHits = true;
2205 circlefactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2206 TrkHitList* trkHitList = trkRecoTrk->hits();
2207
2208 //---convert hits
2209 vector<MdcHit*> mdcHitVector;
2210 mdcHitVector.clear();
2211 vector<CgemHitOnTrack*> cgemHitVector;
2212 cgemHitVector.clear();
2213
2214 for(vector<HoughHit*>::iterator iter = m_vecHitPnt.begin(); iter != m_vecHitPnt.end(); iter++){
2215 HoughHit* hitIt = (*iter);
2216 if(hitIt->getFlag()!=0)continue;
2217 //hitIt->print();
2218 if(hitIt->getHitType()==HoughHit::MDCHIT || hitIt->getHitType()==HoughHit::MDCMCHIT){
2219 const MdcDigi* mdcDigi = hitIt->getDigi();
2220 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2221 mdcHitVector.push_back(mdcHit);
2222 mdcHit->setCalibSvc(hitIt->getMdcCalibFunSvc());
2223 mdcHit->setCountPropTime(true);
2224 //if(mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;//FIXME
2225 int newAmbig = 0;
2226 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2227 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2228 HepPoint3D midPoint = hitIt->getHitPosition();
2229 double fltLength = flightLength(midPoint);
2230 //double fltLength = flightLength(hitIt->getHitPosition());
2231 mdcHitOnTrack->setFltLen(fltLength);
2232 trkHitList->appendHot(mdcHitOnTrack);
2233 }else if(hitIt->getHitType()==HoughHit::CGEMHIT || hitIt->getHitType()==HoughHit::CGEMMCHIT){
2234 if(m_useCgemInGlobalFit<1)continue;
2235 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*(hitIt->getCgemCluster()),bunchT0*1.e9);
2236 cgemHitVector.push_back(cgemHitOnTrack);
2237 cgemHitOnTrack->setCgemGeomSvc(hitIt->getCgemGeomSvc());
2238 cgemHitOnTrack->setCgemCalibFunSvc(hitIt->getCgemCalibFunSvc());
2239 trkHitList->appendHot(cgemHitOnTrack);
2240 }
2241 }
2242 //---fitting
2243 TrkHitList* qhits = trkRecoTrk->hits();
2244 //cout<<"num of qhits to fit = "<<qhits->nHit()<<endl;
2245 //cout<<"before TrkHitList->fit()"<<endl;
2246 TrkErrCode err=qhits->fit();
2247 //cout<<"after TrkHitList->fit()"<<endl;
2248
2249 if(err.success()){
2250 //cout<<" circle fits well "<<endl;
2251 m_circleFitStat=1;
2252 const TrkFit* fitResult = trkRecoTrk->fitResult();
2253 TrkExchangePar par = fitResult->helix(0.);
2254 //if(fabs(m_kappa-par.omega()*fabs(alpha()))/fabs(m_kappa)>0.3)cout<<"fit divergent"<<endl;
2255 m_dr = -par.d0();
2256 m_phi0 = par.phi0()-M_PI/2;
2257 m_kappa = par.omega()*fabs(alpha());
2258 m_dz = par.z0();
2259 m_tanl = par.tanDip();
2260
2261 //cout<<"phi0 = "<<m_phi0<<endl;
2262 while(m_phi0>2.*M_PI) m_phi0-=2.*M_PI;
2263 while(m_phi0<0.) m_phi0+=2.*M_PI;
2264 //cout<<"phi0 = "<<m_phi0<<endl;
2265 //HepVector hepVector(5,0);
2266 //hepVector(1) = m_dr;
2267 //hepVector(2) = m_phi0;
2268 //hepVector(3) = m_kappa;
2269 //hepVector(4) = m_dz;
2270 //hepVector(5) = m_tanl;
2271 //a(hepVector);
2272 updateHelix();
2273 //cout<<"update Helix"<<endl;
2274
2275 //m_trkRecoTrk = trkRecoTrk;
2276 m_chi2 = fitResult->chisq();
2277 if(debug) cout<<"chi2 obtained "<<m_chi2<<endl;
2278 //cout<<"fitting params:"<<endl;
2279 //cout
2280 //<<setw(15)<<par.d0()
2281 //<<setw(15)<<par.phi0()
2282 //<<setw(15)<<par.omega()
2283 //;
2284 //cout<<endl;
2285
2286 //m_vecMdcDigiPnt.clear();
2287 //m_vecCgemClusterPnt.clear();
2288
2289 int nHotKept = 0;
2290 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2291 {
2292 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2293 {
2294 const MdcRecoHitOnTrack* recoHot;
2295 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2296 if(!(recoHot->isActive())) continue;
2297 nHotKept++;
2298 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2299 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2300 }
2301 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2302 {
2303 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2304 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2305 int clusterid = recCgemCluster->getclusterid();
2306 int stat = recoHot->isActive();
2307 if(stat==0) continue;//
2308 nHotKept++;
2309 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2310 }
2311 }
2312 if(debug) cout<<nHotKept<<" hits kept after TrkFitter"<<endl;
2313 }
2314 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2315 delete *iter;
2316 }
2317 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2318 delete *iter;
2319 }
2320 if(m_clearTrack)delete trkRecoTrk;
2321 return err;
2322}
void setCgemCalibFunSvc(const CgemCalibFunSvc *svc)
const RecCgemCluster * baseHit() const
void setCgemGeomSvc(const CgemGeomSvc *svc)
double flightLength(HepPoint3D &hit) const
const MdcDigi * getDigi() const
Definition: HoughHit.h:42
MdcCalibFunSvc * getMdcCalibFunSvc() const
Definition: HoughHit.h:66
CgemGeomSvc * getCgemGeomSvc() const
Definition: HoughHit.h:67
CgemCalibFunSvc * getCgemCalibFunSvc() const
Definition: HoughHit.h:68
static int m_useCgemInGlobalFit
Definition: HoughTrack.h:146
void updateHelix()
Definition: HoughTrack.cxx:923
Definition: MdcHit.h:44
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcHit.cxx:136
void setCountPropTime(const bool count)
Definition: MdcHit.h:86
int getclusterid(void) const
virtual double chisq() const =0
int success() const
Definition: TrkErrCode.h:62
double phi0() const
double omega() const
double z0() const
double d0() const
double tanDip() const
Definition: TrkFit.h:23
virtual TrkExchangePar helix(double fltL) const =0
TrkErrCode fit()
Definition: TrkHitList.cxx:59
hot_iterator begin() const
Definition: TrkHitList.h:45
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition: TrkHitList.cxx:91
hot_iterator end() const
Definition: TrkHitList.h:46
void setFltLen(double f)
Definition: TrkHitOnTrk.h:147
bool isActive() const
Definition: TrkHitOnTrk.h:200
TrkHitList * hits()
Definition: TrkRecoTrk.h:107
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const

◆ fitCircle() [2/2]

int HoughTrack::fitCircle ( DotsHelixFitter fitter,
double  bunchT0 
)

Definition at line 2324 of file HoughTrack.cxx.

2325{
2326 bool debug=false; //debug=true;
2327 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2328 ///*
2329 m_vecMdcDigiPnt.clear();
2330 m_vecCgemClusterPnt.clear();
2331 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2332 {
2333 if((*hotIt)->getFlag()!=0)continue;
2334 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2335 {
2336 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2337 m_vecMdcDigiPnt.push_back(aMdcDigi);
2338 }
2339 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2340 {
2341 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2342 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2343 }
2344 }
2345 //*/
2346
2347 // --- set
2348 HepPoint3D pivot(0,0,0);
2349 HepVector a(5,0);
2350 HepSymMatrix Ea(5,0);
2351
2352 // --- initial circle parameters from Hough transform
2353 double phi0_Hough = (m_rho>0)? m_angle:m_angle+M_PI;
2354 phi0_Hough = (m_charge<0)? phi0_Hough:phi0_Hough+M_PI;
2355 if(phi0_Hough>2*M_PI)phi0_Hough-=2*M_PI;
2356 if(phi0_Hough<0) phi0_Hough+=2*M_PI;
2357 double kappa_Hough = fabs(m_alpha*m_rho)*m_charge;
2358 a(2) = phi0_Hough;
2359 a(3) = kappa_Hough;
2360
2361 // --- from fit?
2362 //a(1) = m_dr;
2363 //a(2) = m_phi0;
2364 //a(3) = m_kappa;
2365 //a(5) = 0.00001;
2366
2367 KalmanFit::Helix ini_helix(pivot,a,Ea);
2368 if(debug) cout<<"circle ini_helix: "<<ini_helix.a()<<endl;
2369 fitter->setInitialHelix(ini_helix);
2370 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2371 fitter->setCgemClusters(m_vecCgemClusterPnt);
2372 if(debug) cout<<"finish DotsHelixFitter setting"<<endl;
2373
2374 // --- circle fit
2375 //fitter->fitCircleOnly();
2376 fitter->fitModeCircle();
2377 int fitFlag = 0;
2378 while(1)
2379 {
2380 fitFlag = fitter->calculateNewHelix();
2381 if(fitFlag==0)
2382 {
2383 HepVector a_new = fitter->getHelix();
2384 if(debug) cout<<"circle helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2385 if( fitter->getChi2() > 25*(fitter->getNActiveHits()) && fitter->deactiveHits(10.0))
2386 {
2387 continue;
2388 if(debug) cout<<"too large chi2, drop one worst DC hit!"<<endl;
2389 }
2390 else break;
2391 }
2392 else
2393 {
2394 if(debug) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2395 break;
2396 }
2397 }
2398 return fitFlag;
2399
2400}
void setInitialHelix(KalmanFit::Helix aHelix)
int deactiveHits(double chi_cut=10, int nMax=1)
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
HepVector getHelix()

◆ fitHelix() [1/3]

TrkErrCode HoughTrack::fitHelix ( const MdcDetector mdcDetector,
const BField bField,
double  bunchT0,
vector< HoughHit * > &  hot,
int  Layer 
)

Definition at line 2543 of file HoughTrack.cxx.

2544{
2545 // --- make TrkRecoTrk and trkHitList
2546 int nHot(0);
2547 double d0 = -m_dr;
2548 double phi0 = m_phi0+M_PI/2;
2549 double omega = m_kappa/fabs(alpha());
2550 double z0 = m_dz;
2551 double tanl = m_tanl;
2552 while(phi0>M_PI)phi0-=2*M_PI;
2553 while(phi0<-M_PI)phi0+=2*M_PI;
2554
2555 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2556 TrkHelixMaker helixfactory;
2557 double chisum =0.;
2558 long idnum(0);
2559 TrkRecoTrk* trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, bfield, bunchT0);
2560 bool permitFlips = true;
2561 bool lPickHits = true;
2562 helixfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2563 TrkHitList* trkHitList = trkRecoTrk->hits();
2564
2565 //--- convert hits and put them into trkHitList
2566 //cout<<maxLayer<<endl;
2567 vector<MdcHit*> mdcHitVector;
2568 mdcHitVector.clear();
2569 vector<CgemHitOnTrack*> cgemHitVector;
2570 cgemHitVector.clear();
2571
2572 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2573 //(*hitIt)->print();
2574 //cout<<endl;
2575 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2576 if((*hitIt)->getFlag()!=0&&(*hitIt)->getLayer()>maxLayer)continue;
2577 bool sameHit(false);
2578 const TrkHitList* aList = trkRecoTrk->hits();
2579 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2580 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2581 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2582 int layerId = recoHot->mdcHit()->layernumber();
2583 int wireId = recoHot->mdcHit()->wirenumber();
2584 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2585 }
2586 }
2587 if(sameHit)continue;
2588
2589 const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2590 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2591 mdcHitVector.push_back(mdcHit);
2592
2593 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2594 mdcHit->setCountPropTime(true);
2595 int newAmbig = 0;
2596 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2597 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2598 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2599 //double fltLength = flightLength(midPoint);
2600 double fltLength = flightArc(midPoint);
2601 //double fltLength = flightLength((*hitIt)->getHitPosition());
2602 mdcHitOnTrack->setFltLen(fltLength);
2603 trkHitList->appendHot(mdcHitOnTrack);
2604 nHot++;
2605 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2606 if(m_useCgemInGlobalFit<3)continue;
2607 bool sameHit(false);
2608 const TrkHitList* aList = trkRecoTrk->hits();
2609 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2610 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2611 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2612 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2613 int clusterid = recCgemCluster->getclusterid();
2614 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2615 }
2616 }
2617 if(sameHit)continue;
2618
2619 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2620 cgemHitVector.push_back(cgemHitOnTrack);
2621 if((*hitIt)->getFlag()==0)continue;
2622 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2623 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2624 trkHitList->appendHot(cgemHitOnTrack);
2625 nHot++;
2626 }
2627 }
2628 //cout<<trkHitList->nHit()<<endl;
2629 //---fitting
2630 TrkHitList* qhits = trkRecoTrk->hits();
2631 TrkErrCode err=qhits->fit();
2632 //err.print(std::cout);cout<<endl;cout<<endl;
2633
2634 //const TrkHitList* aList = trkRecoTrk->hits();
2635 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2636 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2637 //const MdcRecoHitOnTrack* recoHot;
2638 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2639 //int layerId = recoHot->mdcHit()->layernumber();
2640 //int wireId = recoHot->mdcHit()->wirenumber();
2641 //double res=999.,rese=999.;
2642 //if (recoHot->resid(res,rese,false)){
2643 //}else{}
2644 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2645 //}
2646 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2647 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2648 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2649 //int clusterid = recCgemCluster->getclusterid();
2650 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2651 //}
2652 //}
2653
2654 // --- if fit successes, update the helix parameters
2655 if(err.success()){
2656 const TrkFit* fitResult = trkRecoTrk->fitResult();
2657 TrkExchangePar par = fitResult->helix(0.);
2658 m_dr = -par.d0();
2659 m_phi0 = par.phi0()-M_PI/2;
2660 m_kappa = par.omega()*fabs(alpha());
2661 m_dz = par.z0();
2662 m_tanl = par.tanDip();
2663 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2664 while(m_phi0<0)m_phi0+=2*M_PI;
2665 updateHelix();
2666 m_chi2 = fitResult->chisq();
2667
2668 int nHotKept = 0;
2669 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2670 {
2671 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2672 {
2673 const MdcRecoHitOnTrack* recoHot;
2674 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2675 if(!(recoHot->isActive())) continue;
2676 nHotKept++;
2677 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2678 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2679 }
2680 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2681 {
2682 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2683 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2684 int clusterid = recCgemCluster->getclusterid();
2685 int stat = recoHot->isActive();
2686 if(stat==0) continue;//
2687 nHotKept++;
2688 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2689 }
2690 }
2691
2692 cout<<"chi2= "<<m_chi2
2693 <<", "<<nHotKept<<" hits kept after TrkFitter"
2694 <<", helix from TrkFitter = "<<a()
2695 <<endl;
2696 }
2697 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2698 delete *iter;
2699 }
2700 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2701 delete *iter;
2702 }
2703 if(m_clearTrack)delete trkRecoTrk;
2704 //cout<<"getId():"<<trkRecoTrk->id()<<endl;
2705 return err;
2706}
unsigned layernumber() const
Definition: MdcHit.h:61
unsigned wirenumber() const
Definition: MdcHit.h:62
const MdcHit * mdcHit() const

◆ fitHelix() [2/3]

TrkErrCode HoughTrack::fitHelix ( const MdcDetector mdcDetector,
TrkContextEv trkContextEv,
double  bunchT0,
vector< MdcHit * > &  mdcHitCol,
vector< HoughHit * > &  hot 
)

Definition at line 2402 of file HoughTrack.cxx.

2403{
2404 int nHot(0);
2405 double d0 = -m_dr;
2406 double phi0 = m_phi0+M_PI/2;
2407 double omega = m_kappa/fabs(alpha());
2408 double z0 = m_dz;
2409 double tanl = m_tanl;
2410 while(phi0>M_PI)phi0-=2*M_PI;
2411 while(phi0<-M_PI)phi0+=2*M_PI;
2412
2413 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2414 TrkHelixMaker helixfactory;
2415 double chisum =0.;
2416 m_trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, *trkContextEv, bunchT0);
2417 bool permitFlips = true;
2418 bool lPickHits = true;
2419 helixfactory.setFlipAndDrop(*m_trkRecoTrk, permitFlips, lPickHits);
2420 TrkHitList* trkHitList = m_trkRecoTrk->hits();
2421
2422 //--- convert hits
2423 //cout<<maxLayer<<endl;
2424 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2425 //(*hitIt)->print();
2426 //cout<<endl;
2427 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2428 bool sameHit(false);
2429 const TrkHitList* aList = m_trkRecoTrk->hits();
2430 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2431 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2432 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2433 int layerId = recoHot->mdcHit()->layernumber();
2434 int wireId = recoHot->mdcHit()->wirenumber();
2435 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2436 }
2437 }
2438 if(sameHit)continue;
2439
2440 //const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2441 //MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2442 MdcHit *mdcHit;
2443 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
2444 for(;imdcHit != mdcHitCol.end();imdcHit++){
2445 if((*imdcHit)->mdcId() == (*hitIt)->getDigi()->identify()){
2446 mdcHit = *imdcHit;
2447 break;
2448 }
2449 }
2450
2451 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2452 mdcHit->setCountPropTime(true);
2453 int newAmbig = 0;
2454 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2455 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2456 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2457 //double fltLength = flightLength(midPoint);
2458 double fltLength = flightArc(midPoint);
2459 //double fltLength = flightLength((*hitIt)->getHitPosition());
2460 mdcHitOnTrack->setFltLen(fltLength);
2461 trkHitList->appendHot(mdcHitOnTrack);
2462 nHot++;
2463 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2464 if(m_useCgemInGlobalFit<2)continue;
2465 bool sameHit(false);
2466 const TrkHitList* aList = m_trkRecoTrk->hits();
2467 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2468 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2469 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2470 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2471 int clusterid = recCgemCluster->getclusterid();
2472 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2473 }
2474 }
2475 if(sameHit)continue;
2476
2477 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2478 m_cgemHitVector.push_back(cgemHitOnTrack);
2479 if((*hitIt)->getFlag()==0)continue;
2480 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2481 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2482 trkHitList->appendHot(cgemHitOnTrack);
2483 nHot++;
2484 }
2485 }
2486 //cout<<trkHitList->nHit()<<endl;
2487 //---fitting
2488 TrkHitList* qhits = m_trkRecoTrk->hits();
2489 TrkErrCode err=qhits->fit();
2490 //err.print(std::cout);
2491 //cout<<endl;
2492 //cout<<endl;
2493
2494 //const TrkHitList* aList = m_trkRecoTrk->hits();
2495 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2496 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2497 //const MdcRecoHitOnTrack* recoHot;
2498 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2499 //int layerId = recoHot->mdcHit()->layernumber();
2500 //int wireId = recoHot->mdcHit()->wirenumber();
2501 //double res=999.,rese=999.;
2502 //if (recoHot->resid(res,rese,false)){
2503 //}else{}
2504 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2505 //}
2506 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2507 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2508 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2509 //int clusterid = recCgemCluster->getclusterid();
2510 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2511 //}
2512 //}
2513
2514 if(err.success()){
2515 const TrkFit* fitResult = m_trkRecoTrk->fitResult();
2516 TrkExchangePar par = fitResult->helix(0.);
2517 m_dr = -par.d0();
2518 m_phi0 = par.phi0()-M_PI/2;
2519 m_kappa = par.omega()*fabs(alpha());
2520 m_dz = par.z0();
2521 m_tanl = par.tanDip();
2522 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2523 while(m_phi0<0)m_phi0+=2*M_PI;
2524 //HepVector hepVector(5,0);
2525 //a(hepVector);
2526 updateHelix();
2527
2528 m_chi2 = fitResult->chisq();
2529 //cout<<"fitting params:"<<endl;
2530 //cout
2531 //<<setw(15)<<par.d0()
2532 //<<setw(15)<<par.phi0()
2533 //<<setw(15)<<par.omega()
2534 //<<setw(15)<<par.z0()
2535 //<<setw(15)<<par.tanDip()
2536 //;
2537 //cout<<endl;
2538 }
2539 //cout<<"getId()="<<m_trkRecoTrk->id()<<endl;
2540 return err;
2541}
TrkFundHit::hot_iterator begin() const
Definition: TrkFundHit.h:113

◆ fitHelix() [3/3]

int HoughTrack::fitHelix ( DotsHelixFitter fitter,
double  bunchT0,
RecCgemClusterCol::iterator  recCgemClusterColBegin,
double  averageChi2cut = 25 
)

Definition at line 2708 of file HoughTrack.cxx.

2709{
2710 bool debug=false; //debug=true;
2711 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2712 //cout<<"before fill"<<endl;
2713 m_vecMdcDigiPnt.clear();
2714 m_vecCgemClusterPnt.clear();
2715 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2716 {
2717 //if((*hotIt)->getFlag()!=0)continue;
2718 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2719 {
2720 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2721 m_vecMdcDigiPnt.push_back(aMdcDigi);
2722 }
2723 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2724 {
2725 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2726 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2727 }
2728 }
2729 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++)
2730 {
2731 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2732 {
2733 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2734 m_vecMdcDigiPnt.push_back(aMdcDigi);
2735 }
2736 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2737 {
2738 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2739 if(aRecCgemCluster->getflag()==2)
2740 {
2741 int clusterId = aRecCgemCluster->getclusterflage();
2742 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterColBegin+clusterId;
2743 //cout<<"V cluster ID = "<<setw(10)<<clusterId<<setw(10)<<(*cgemClusterIter)->getclusterid();
2744 if(clusterId==(*cgemClusterIter)->getclusterid())
2745 m_vecCgemClusterPnt.push_back(*cgemClusterIter);
2746 }
2747 }
2748 }
2749 //cout<<"end filling"<<endl;
2750
2751
2752 // --- set
2753 HepPoint3D pivot(0,0,0);
2754 HepVector a(5,0);
2755 HepSymMatrix aEa(5,0);
2756 a(1) = m_dr;
2757 a(2) = m_phi0;
2758 a(3) = m_kappa;
2759 a(4) = m_dz;
2760 a(5) = m_tanl;
2761 KalmanFit::Helix ini_helix(pivot,a,aEa);
2762 if(debug) cout<<"fitHelix ini_helix: "<<ini_helix.a()<<endl;
2763 fitter->setInitialHelix(ini_helix);
2764 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2765 fitter->setCgemClusters(m_vecCgemClusterPnt);
2766 int nMdcDigi = m_vecMdcDigiPnt.size();
2767
2768 // --- helix fit
2769 fitter->fitModeHelix();
2770 int fitFlag = 0;
2771 while(1)
2772 {
2773 fitFlag = fitter->calculateNewHelix();
2774 if(fitFlag==0)
2775 {
2776 HepVector a_new = fitter->getHelix();
2777 if(debug) {
2778 cout<<"helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2779 cout<<"nMdcDigi = "<<nMdcDigi<<endl;
2780 }
2781 //if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits())
2782 if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits(43,--nMdcDigi))
2783 {
2784 //cout<<"too large chi2, drop the worst DC hit!"<<endl;
2785 if(debug) cout<<"too large chi2, drop the outmost DC hit!"<<endl;
2786 continue;
2787 }
2788 else
2789 {
2790 if(debug) cout<<"stop droping hit! "<<endl;
2791 break;
2792 }
2793 }
2794 else
2795 {
2796 if(debug) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2797 break;
2798 }
2799 }
2800 if(fitFlag==0) //update dr ... and Ea
2801 {
2802 HepVector a = fitter->getHelix();
2803 m_dr = a[0];
2804 m_phi0 = a[1];
2805 m_kappa = a[2];
2806 m_dz = a[3];
2807 m_tanl = a[4];
2808 m_chi2 = fitter->getChi2();
2809 Ea(fitter->getEa());
2810 }
2811 return fitFlag;
2812}
const HepSymMatrix & getEa()
int getflag(void) const
int getclusterflage(void) const

◆ getAngle()

double HoughTrack::getAngle ( ) const
inline

Definition at line 52 of file HoughTrack.h.

52{return m_angle;}

◆ getCharge()

int HoughTrack::getCharge ( ) const
inline

Definition at line 50 of file HoughTrack.h.

50{return m_charge;}

Referenced by HoughFinder::fillHistogram().

◆ getChi2()

double HoughTrack::getChi2 ( ) const
inline

Definition at line 58 of file HoughTrack.h.

58{return m_chi2;}

Referenced by print().

◆ getCircleFitStat()

int HoughTrack::getCircleFitStat ( ) const
inline

Definition at line 59 of file HoughTrack.h.

59{return m_circleFitStat;};

◆ getDAngle()

double HoughTrack::getDAngle ( ) const
inline

Definition at line 54 of file HoughTrack.h.

54{return m_dAngle;}

◆ getDDz()

double HoughTrack::getDDz ( ) const
inline

Definition at line 57 of file HoughTrack.h.

57{return m_dDz;}

◆ getDr()

double HoughTrack::getDr ( ) const
inline

Definition at line 63 of file HoughTrack.h.

63{return m_dr;}

◆ getDRho()

double HoughTrack::getDRho ( ) const
inline

Definition at line 55 of file HoughTrack.h.

55{return m_dRho;}

◆ getDTanl()

double HoughTrack::getDTanl ( ) const
inline

Definition at line 56 of file HoughTrack.h.

56{return m_dTanl;}

◆ getDz()

double HoughTrack::getDz ( ) const
inline

Definition at line 66 of file HoughTrack.h.

66{return m_dz;}

◆ getFlag()

int HoughTrack::getFlag ( ) const
inline

Definition at line 51 of file HoughTrack.h.

51{return m_flag;}

◆ getHotList()

vector< HoughHit * > HoughTrack::getHotList ( int  type = 2)

Definition at line 2127 of file HoughTrack.cxx.

2128{
2129 vector<HoughHit*> hotList;
2130 vector<HoughHit*> axialHot = getVecHitPnt();
2131 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
2132
2133 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2134 if(type==1)continue;
2135 if(type==2&&(**hotIt).getHitType()==HoughHit::CGEMHIT)continue;
2136 hotList.push_back(*hotIt);
2137 }
2138 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2139 if(type>0)hotList.push_back(*hotIt);
2140 }
2141
2142 /*
2143 // --- CGEM XV cluster
2144 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2145 if((**hotIt).getLayer()<3)hotList.push_back(*hotIt);
2146 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2147 }
2148
2149 // --- MDC hits, layer 9-20
2150 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2151 //if(m_fitFlag==1&&(**hotIt).getLayer()<3)hot.push_back(**hotIt);
2152 if((**hotIt).getLayer()>3&&(**hotIt).getLayer()<20)hotList.push_back(*hotIt);
2153 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2154 }
2155
2156 // --- MDC hits, layer 21-36
2157 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2158 if((**hotIt).getLayer()>=20&&(**hotIt).getLayer()<36)hotList.push_back(*hotIt);
2159 //int layer = (*hotIt)->getLayer();
2160 //int wire = (*hotIt)->getWire();
2161 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2162 //cout<<(**hotIt).residual(&(*trkIT));
2163 //cout<<endl;
2164 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2165 }
2166
2167 // --- MDC hits, layer 37-43
2168 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2169 if((**hotIt).getLayer()>=36)hotList.push_back(*hotIt);
2170 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2171 }
2172 //cout<<"nHot: "<<hot.size()<<endl;
2173 // */
2174
2175 sortHot(hotList);
2176 return hotList;
2177}
void sortHot(vector< HoughHit * > &hotList)
Definition: HoughTrack.cxx:918
vector< HoughHit * > getVecHitPnt()
Definition: HoughTrack.h:83
vector< HoughHit * > getVecStereoHitPnt()
Definition: HoughTrack.h:134

Referenced by moreHot(), print(), and printHot().

◆ getKappa()

double HoughTrack::getKappa ( ) const
inline

Definition at line 65 of file HoughTrack.h.

65{return m_kappa;}

◆ getMcTrack()

HoughTrack * HoughTrack::getMcTrack ( ) const
inline

Definition at line 141 of file HoughTrack.h.

141{return m_mcTrack;}

◆ getNhitFirstHalf()

int HoughTrack::getNhitFirstHalf ( )
inline

Definition at line 92 of file HoughTrack.h.

92{return m_nHitFirstHalf;};

◆ getNhitSecondHalf()

int HoughTrack::getNhitSecondHalf ( )
inline

Definition at line 93 of file HoughTrack.h.

93{return m_nHitSecondHalf;};

◆ getNHitsShared()

int HoughTrack::getNHitsShared ( )

Definition at line 1983 of file HoughTrack.cxx.

1984{
1985 int nShared = 0;
1986 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1987 for(; iter!=m_vecHitPnt.end(); iter++)
1988 {
1989 vector<double> vecRes = (*iter)->getVecResid();
1990 int nTrk = vecRes.size();
1991 if(nTrk>1) nShared++;
1992 }
1993 return nShared;
1994}

◆ getNhitUnusedFirstHalf()

int HoughTrack::getNhitUnusedFirstHalf ( )
inline

Definition at line 94 of file HoughTrack.h.

94{return m_nHitUnused_FirstHalf;};

◆ getNhitUnusedSecondHalf()

int HoughTrack::getNhitUnusedSecondHalf ( )
inline

Definition at line 95 of file HoughTrack.h.

95{return m_nHitUnused_SecondHalf;};

◆ getNTried()

int HoughTrack::getNTried ( )
inline

Definition at line 129 of file HoughTrack.h.

129{return m_untried;};// llwang 2018-12-23

◆ getPhi0()

double HoughTrack::getPhi0 ( ) const
inline

Definition at line 64 of file HoughTrack.h.

64{return m_phi0;}

◆ getRecMdcHitVec()

vector< RecMdcHit > * HoughTrack::getRecMdcHitVec ( )
inline

Definition at line 151 of file HoughTrack.h.

151{return &m_vecRecMdcHit;};

◆ getRho()

double HoughTrack::getRho ( ) const
inline

Definition at line 53 of file HoughTrack.h.

53{return m_rho;}

◆ getTanl()

double HoughTrack::getTanl ( ) const
inline

Definition at line 67 of file HoughTrack.h.

67{return m_tanl;}

◆ getTrkID()

int HoughTrack::getTrkID ( ) const
inline

Definition at line 49 of file HoughTrack.h.

49{return m_trkID;}

Referenced by HoughFinder::getMcParticleCol(), and markUsedHot().

◆ getTrkRecoTrk()

TrkRecoTrk * HoughTrack::getTrkRecoTrk ( )
inline

Definition at line 62 of file HoughTrack.h.

62{return m_trkRecoTrk;}

◆ getVecHitPnt()

vector< HoughHit * > HoughTrack::getVecHitPnt ( )
inline

Definition at line 83 of file HoughTrack.h.

83{return m_vecHitPnt;}; // llwang 2018-12-26

Referenced by getHotList(), and print().

◆ getVecHitRes()

vector< double > HoughTrack::getVecHitRes ( )
inline

Definition at line 133 of file HoughTrack.h.

133{return m_vecHitResidual;}

◆ getVecStereoHitPnt()

vector< HoughHit * > HoughTrack::getVecStereoHitPnt ( )
inline

Definition at line 134 of file HoughTrack.h.

134{return m_vecStereoHitPnt;}; // llwang 2018-12-26

Referenced by getHotList(), and print().

◆ getVecStereoHitRes()

vector< double > HoughTrack::getVecStereoHitRes ( )
inline

Definition at line 136 of file HoughTrack.h.

136{return m_vecStereoHitRes;}

◆ isAGoodCircle()

bool HoughTrack::isAGoodCircle ( )

Definition at line 1507 of file HoughTrack.cxx.

1508{
1509 bool debug = false; //debug=true;
1510 int NtotLayer = m_setLayer.size();
1511
1512 std::set<int>::iterator it = m_setLayer.begin();
1513 int lastLayer = *it;
1514 // check N_CGEM_layer, N_ODC1_nearStereo, N_ODC2_nearStereo
1515 int nCgemLayer(0), N_ODC1_nearStereo(0), N_ODC2_nearStereo(0);
1516 for(; it!=m_setLayer.end(); it++)
1517 {
1518 if((*it)<3) nCgemLayer++;
1519 else if((*it)>=16&&(*it)<=19) N_ODC1_nearStereo++;
1520 else if((*it)>=36&&(*it)<=39) N_ODC2_nearStereo++;
1521 }
1522 bool enoughVHits;
1523 //enoughVHits = nCgemLayer==3 ;
1524 enoughVHits = nCgemLayer>=2
1525 || (nCgemLayer>0 && (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2))
1526 //|| (N_ODC1_nearStereo>=2 && N_ODC2_nearStereo>=2 && NtotLayer>=8) ;
1527 || ( (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2) && NtotLayer>=8) ; // modified in 2020-08
1528
1529 // gap fraction
1530 double r_gap = 1.0*m_nGap/(m_nGap+NtotLayer);
1531 bool smallGap = r_gap<0.5;
1532 if(debug)
1533 {
1534 if(!enoughVHits) cout<<"not enough V hits, ";
1535 if(!smallGap) cout<<"gap too large 50%, "<<endl;
1536 if(m_unused<3) cout<<"N_unused = "<<m_unused<<"<3"<<endl;
1537 if(m_untried<3) cout<<"N_untried = "<<m_untried<<"<3"<<endl;
1538 }
1539
1540 // --- check conditions
1541 bool good = true;
1542 good=good && m_unused>=2; // 3 -> 2, modified 2020-10
1543 good=good && m_untried>=2;// 3 -> 2, modified 2020-10
1544 good=good && (nCgemLayer==3||NtotLayer>3);
1545 good=good && enoughVHits;
1546 good=good && smallGap;
1547 //good=good && m_maxGap<4;
1548 return good;
1549
1550}

◆ judgeCharge() [1/2]

int HoughTrack::judgeCharge ( double  xHit,
double  yHit 
)

Definition at line 622 of file HoughTrack.cxx.

623{
624 //xHit = hit->getHitPosition().x();
625 //yHit = hit->getHitPosition().y();
626 //double xCenter = m_helix.center().x();
627 //double yCenter = m_helix.center().y();
628 double xCenter = center().x();
629 double yCenter = center().y();
630 double leftOrRight = xHit*yCenter - xCenter*yHit;
631 if(leftOrRight>0) return 1;
632 else return -1;
633}

◆ judgeCharge() [2/2]

int HoughTrack::judgeCharge ( HoughHit hit)

Definition at line 609 of file HoughTrack.cxx.

610{
611 double xHit = hit->getHitPosition().x();
612 double yHit = hit->getHitPosition().y();
613 //double xCenter = m_helix.center().x();
614 //double yCenter = m_helix.center().y();
615 double xCenter = center().x();
616 double yCenter = center().y();
617 double leftOrRight = xHit*yCenter - xCenter*yHit;
618 if(leftOrRight>0) return 1;
619 else return -1;
620}

Referenced by findVHot(), and findXHot().

◆ judgeHalf()

int HoughTrack::judgeHalf ( HoughHit hit)

Definition at line 596 of file HoughTrack.cxx.

597{
598 double xHit = hit->getHitPosition().x();
599 double yHit = hit->getHitPosition().y();
600 //double xCenter = m_helix.center().x();
601 //double yCenter = m_helix.center().y();
602 double xCenter = center().x();
603 double yCenter = center().y();
604 if(xHit*yCenter > xCenter*yHit)return m_charge;
605 else if(xHit*yCenter < xCenter*yHit)return -m_charge;
606 else return 0;
607}

Referenced by driftDistRes(), HoughFinder::fillHistogram(), findXHot(), HoughFinder::getMcParticleCol(), and HoughHit::residual().

◆ markUsedHot() [1/2]

void HoughTrack::markUsedHot ( int  use = 1)

Definition at line 1484 of file HoughTrack.cxx.

1485{
1486 vector<double>::iterator it0_res = m_vecHitResidual.begin();
1487 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1488 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1489 for(; iter!=m_vecHitPnt.end(); iter++)
1490 {
1491 int useOld = (*iter)->getUse();
1492 if(use==-1&&useOld>0) continue;
1493 (*iter)->setUse(use);
1494 if(use==1) {
1495 //cout<<"HoughTrack::markUsedHot(use) add trk ID "<<getTrkID()<<endl;
1496 //(*iter)->addTrkPnt(this);
1497 (*iter)->addTrkID(getTrkID());
1498 vector<double>::iterator it_res = it0_res+(iter-iter0);
1499 //cout<<"add residual "<<*it_res<<endl;
1500 (*iter)->addResid(*it_res);
1501 //vector<double> vecRes = (*iter)->getVecResid();
1502 //cout<<"VecResid.size = "<<vecRes.size()<<endl;
1503 }
1504 }
1505}
int getTrkID() const
Definition: HoughTrack.h:49

◆ markUsedHot() [2/2]

void HoughTrack::markUsedHot ( vector< HoughHit * > &  hitPntList,
int  use = 1 
)

Definition at line 1467 of file HoughTrack.cxx.

1468{
1469 vector<HoughHit*>::iterator iter = hitPntList.begin();
1470 for(; iter!=hitPntList.end(); iter++)
1471 {
1472 int useOld = (*iter)->getUse();
1473 if(use==-1&&useOld>0) continue;
1474 (*iter)->setUse(use);
1475
1476 if(use==1) {
1477 //cout<<"HoughTrack::markUsedHot add trk ID "<<getTrkID()<<endl;
1478 //(*iter)->addTrkPnt(this);
1479 (*iter)->addTrkID(getTrkID());
1480 }
1481 }
1482}

◆ operator=()

HoughTrack & HoughTrack::operator= ( const HoughTrack other)

Definition at line 255 of file HoughTrack.cxx.

256{
257 if (this == & other) return * this;
258
259 Helix::operator=(other);
260 m_trkID = other.m_trkID;
261 m_charge = other.m_charge;
262 m_angle = other.m_angle;
263 m_rho = other.m_rho;
264 m_flag = other.m_flag;
265
266 m_dAngle = other.m_dAngle;
267 m_dRho = other.m_dRho;
268 m_dTanl = other.m_dTanl;
269 m_dDz = other.m_dDz;
270
271 m_dr = other.m_dr;
272 m_phi0 = other.m_phi0;
273 m_kappa = other.m_kappa;
274 m_dz = other.m_dz;
275 m_tanl = other.m_tanl;
276
277 m_chi2 = other.m_chi2;
278 m_trkRecoTrk = other.m_trkRecoTrk;
279
280 m_circleFitStat=other.m_circleFitStat;
281 m_vecHitPnt = other.m_vecHitPnt;
282 m_vecHitResidual = other.m_vecHitResidual;
283 m_vecHitChi2 = other.m_vecHitChi2;
284
285 m_vecStereoHitPnt = other.m_vecStereoHitPnt;
286 m_vecStereoHitRes = other.m_vecStereoHitRes;
287 m_mcTrack = other.m_mcTrack;
288 m_cgemHitVector = other.m_cgemHitVector;
289 m_vecRecMdcHit = other.m_vecRecMdcHit;
290 return * this;
291}
Helix & operator=(const Helix &)
Copy operator.

◆ print()

void HoughTrack::print ( )

Definition at line 962 of file HoughTrack.cxx.

963{
964 cout
965 <<setw(12)<<"trkID:" <<setw(15)<<m_trkID
966 <<setw(12)<<"flag:" <<setw(15)<<m_flag
967 //<<setw(12)<<"q:" <<setw(15)<<m_charge
968 <<setw(12)<<"pt:" <<setw(15)<<pt()
969 <<setw(12)<<"angle:" <<setw(15)<<m_angle
970 <<setw(12)<<"rho:" <<setw(15)<<m_rho
971 <<endl
972 <<setw(12)<<"dr:" <<setw(15)<<dr()
973 <<setw(12)<<"phi0:" <<setw(15)<<phi0()
974 <<setw(12)<<"kappa:" <<setw(15)<<kappa()
975 <<setw(12)<<"dz:" <<setw(15)<<dz()
976 <<setw(12)<<"tanl:" <<setw(15)<<tanl()
977 <<setw(12)<<"chi2:" <<setw(15)<<getChi2()
978 <<endl
979 //<<setw(12)<<"dr:" <<setw(15)<<m_dr
980 //<<setw(12)<<"phi0:" <<setw(15)<<m_phi0
981 //<<setw(12)<<"kappa:" <<setw(15)<<m_kappa
982 //<<setw(12)<<"dz:" <<setw(15)<<m_dz
983 //<<setw(12)<<"tanl:" <<setw(15)<<m_tanl
984 //<<endl
985 //<<setw(12)<<"dAngle:" <<setw(15)<<m_dAngle
986 //<<setw(12)<<"dRho:" <<setw(15)<<m_dRho
987 //<<setw(12)<<"dTanl:" <<setw(15)<<m_dTanl
988 //<<setw(12)<<"dDz:" <<setw(15)<<m_dDz
989 //<<setw(12)<<"chi2:" <<setw(15)<<m_chi2
990 //<<setw(12)<<":"<<setw(15)<<
991 //<<setw(12)<<":"<<setw(15)<<
992 //<<setw(12)<<":"<<setw(15)<<
993 ;
994 int nHot = getHotList().size();
995 int nAxialHot = 0;
996 int nStereoHot = 0;
997 int nXCluster = 0;
998 int nXVCluster = 0;
999 vector<HoughHit*> hotList;
1000 vector<HoughHit*> axialHot = getVecHitPnt();
1001 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
1002 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
1003 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXCluster++;
1004 if((**hotIt).getHitType()==HoughHit::MDCHIT)nAxialHot++;
1005 }
1006 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
1007 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXVCluster++;
1008 if((**hotIt).getHitType()==HoughHit::MDCHIT)nStereoHot++;
1009 }
1010 cout
1011 <<setw(12)<<"nHot:" <<setw(15)<<nHot
1012 <<setw(12)<<"nAxialHot0:" <<setw(15)<<nAxialHot
1013 <<setw(12)<<"nStereoHot0:" <<setw(15)<<nStereoHot
1014 <<setw(12)<<"nXCluster0:" <<setw(15)<<nXCluster
1015 <<setw(12)<<"nXVCluster0:" <<setw(15)<<nXVCluster
1016 <<endl;
1017 //cout<<endl;
1018}
double getChi2() const
Definition: HoughTrack.h:58
vector< HoughHit * > getHotList(int type=2)

Referenced by HoughFinder::getMcParticleCol().

◆ printHot()

void HoughTrack::printHot ( )

Definition at line 1020 of file HoughTrack.cxx.

1021{
1022 vector<HoughHit*> hitList = getHotList();
1023 //vector<HoughHit*> hitList = getVecStereoHitPnt();
1024 for(vector<HoughHit*>::iterator hotIt = hitList.begin(); hotIt != hitList.end();hotIt++){
1025 HoughHit* hitIt = *hotIt;
1026 hitIt->print();
1027 //if(hitIt->getLayer()>3)continue;
1028 //cout<<hitIt->getHitID();
1029 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
1030 //if(hitIt->getHitType()==HoughHit::CGEMMCHIT||hitIt->getHitType()==HoughHit::MDCMCHIT)cout<<" halfCircle:"<<hitIt->getHalfCircle();
1031 //else if(hitIt->getPairHit()!=NULL)cout<<" halfCircle:"<<hitIt->getPairHit()->getHalfCircle();
1032 //HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1033 //double res = hitIt->residual(this, positionOntrack, positionOnDetector);
1034 //cout<<" res:"<<setw(10)<<res<<" ";
1035 //cout<<endl;
1036 }
1037 cout<<endl;
1038 /*
1039 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1040 cout<<setw(4)<<(*hitIt)->getHitID();
1041 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1042 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1043 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1044 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1045 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1046 vector<int> trkID = (*hitIt)->getTrkID();
1047 cout<<"TrkID:"<<setw(4)<<trkID[0];
1048 }
1049 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1050 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1051 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1052 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1053 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1054 cout<<"res:"<<setw(10)<<res<<" ";
1055 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1056 }
1057 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1058 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1059 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1060 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1061 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1062 cout<<"res:"<<setw(10)<<res<<" ";
1063 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1064 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1065 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1066 }
1067 if((*hitIt)->getPairHit()!=NULL){
1068 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1069 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1070 }
1071 //HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1072 //cout<<flightLength(hitPoint);
1073 //else cout<<"NULL";
1074 cout<<endl;
1075 }
1076
1077 //vector<HoughHit*> hitList = getVecHitPnt();
1078 hitList = getVecStereoHitPnt();
1079 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1080 cout<<setw(4)<<(*hitIt)->getHitID();
1081 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1082 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1083 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1084 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1085 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1086 vector<int> trkID = (*hitIt)->getTrkID();
1087 cout<<"TrkID:"<<setw(4)<<trkID[0];
1088 }
1089 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1090 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1091 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1092 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1093 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1094 cout<<"res:"<<setw(10)<<res<<" ";
1095 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1096 }
1097 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1098 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1099 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1100 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1101 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1102 cout<<"res:"<<setw(10)<<res<<" ";
1103 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1104 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1105 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1106 }
1107 if((*hitIt)->getPairHit()!=NULL){
1108 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1109 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1110 }
1111//HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1112//cout<<flightLength(hitPoint);
1113//else cout<<"NULL";
1114cout<<endl;
1115}
1116//*/
1117}
void print()
Definition: HoughHit.cxx:327

◆ printRecMdcHitVec()

void HoughTrack::printRecMdcHitVec ( )

Definition at line 1119 of file HoughTrack.cxx.

1120{
1121 cout<<"HoughTrack::printRecMdcHitVec(): "<<endl;
1122 vector<RecMdcHit>::iterator iter_recMdcHit = m_vecRecMdcHit.begin();
1123 for(; iter_recMdcHit!=m_vecRecMdcHit.end(); iter_recMdcHit++)
1124 {
1125 Identifier mdcid = iter_recMdcHit->getMdcId();
1126 int layer = MdcID::layer(mdcid);
1127 int wire = MdcID::wire(mdcid);
1128 cout<<" hit at layer "<<layer<<" wire "<<wire<<endl;
1129 }
1130}
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
static int wire(const Identifier &id)
Definition: MdcID.cxx:54

◆ releaseSelHits()

void HoughTrack::releaseSelHits ( )

Definition at line 1852 of file HoughTrack.cxx.

1853{
1854 vector<HoughHit*>::iterator it = m_vecHitPnt.begin();
1855 for(; it!=m_vecHitPnt.end(); it++)
1856 {
1857 (*it)->dropTrkID(m_trkID);// erase trk id and residual from the Hough hit
1858 }
1859}

◆ resetNhitHalf()

void HoughTrack::resetNhitHalf ( )
inline

Definition at line 96 of file HoughTrack.h.

96{m_nHitFirstHalf=0; m_nHitSecondHalf=0; m_nHitUnused_FirstHalf=0; m_nHitUnused_SecondHalf=0;};

◆ setAngle()

void HoughTrack::setAngle ( double  angle)
inline

Definition at line 37 of file HoughTrack.h.

37{m_angle = angle;}

◆ setCharge()

void HoughTrack::setCharge ( int  charge)
inline

Definition at line 36 of file HoughTrack.h.

36{m_charge = charge;}

◆ setChi2()

void HoughTrack::setChi2 ( double  chi2)
inline

Definition at line 43 of file HoughTrack.h.

43{m_chi2 = chi2;}

◆ setDAngle()

void HoughTrack::setDAngle ( double  dAngle)
inline

Definition at line 39 of file HoughTrack.h.

39{m_dAngle = dAngle;}

◆ setDDz()

void HoughTrack::setDDz ( double  dDz)
inline

Definition at line 42 of file HoughTrack.h.

42{m_dDz = dDz;}

◆ setDRho()

void HoughTrack::setDRho ( double  dRho)
inline

Definition at line 40 of file HoughTrack.h.

40{m_dRho = dRho;}

◆ setDTanl()

void HoughTrack::setDTanl ( double  dTanl)
inline

Definition at line 41 of file HoughTrack.h.

41{m_dTanl = dTanl;}

◆ setDz()

void HoughTrack::setDz ( double  dz)
inline

Definition at line 46 of file HoughTrack.h.

46{m_dz = dz;}

◆ setFlag()

void HoughTrack::setFlag ( int  flag)
inline

Definition at line 35 of file HoughTrack.h.

35{m_flag = flag;}

◆ setMcTrack()

void HoughTrack::setMcTrack ( HoughTrack mcTrack)
inline

Definition at line 142 of file HoughTrack.h.

142{m_mcTrack = mcTrack;}

◆ setRecMdcHitVec()

void HoughTrack::setRecMdcHitVec ( vector< RecMdcHit > &  aRecMdcHitVec)
inline

Definition at line 150 of file HoughTrack.h.

150{ m_vecRecMdcHit=aRecMdcHitVec; };

◆ setRho()

void HoughTrack::setRho ( double  rho)
inline

Definition at line 38 of file HoughTrack.h.

38{m_rho = rho;}

◆ setTanl()

void HoughTrack::setTanl ( double  tanl)
inline

Definition at line 47 of file HoughTrack.h.

47{m_tanl = tanl;}

◆ setTrkID()

void HoughTrack::setTrkID ( int  trkID)
inline

Definition at line 34 of file HoughTrack.h.

34{m_trkID = trkID;}

◆ sortHot()

void HoughTrack::sortHot ( vector< HoughHit * > &  hotList)

Definition at line 918 of file HoughTrack.cxx.

919{
920 std::sort(hotList.begin(),hotList.end(),sortByLayer);
921}
bool sortByLayer(HoughHit *hit1, HoughHit *hit2)
Definition: HoughTrack.cxx:913

Referenced by getHotList().

◆ update()

void HoughTrack::update ( double  angle,
double  rho 
)

Definition at line 935 of file HoughTrack.cxx.

936{
937 m_angle = angle;
938 m_rho = rho;
939 m_phi0 = (rho>0)? angle:angle+M_PI;
940 m_phi0 = (m_charge<0)? m_phi0:m_phi0+M_PI;
941 if(m_phi0>2*M_PI)m_phi0-=2*M_PI;
942 if(m_phi0<0) m_phi0+=2*M_PI;
943 m_kappa = fabs(m_alpha*rho)*m_charge;
944 HepPoint3D pivot(0,0,0);
945 HepVector a(5,0);
946 HepSymMatrix Ea(5,0);
947 a(2) = m_phi0;
948 a(3) = m_kappa;
949 set(pivot,a,Ea);
950}

◆ updateCirclePar()

void HoughTrack::updateCirclePar ( double  dr,
double  phi0,
double  kappa 
)

Definition at line 953 of file HoughTrack.cxx.

954{
955 m_dr = dr;
956 m_phi0=phi0;
957 m_kappa=kappa;
958 updateHelix();
959
960}

◆ updateHelix()

void HoughTrack::updateHelix ( )

Definition at line 923 of file HoughTrack.cxx.

924{
925 pivot(HepPoint3D(0,0,0));
926 HepVector hepVector(5,0);
927 hepVector(1) = m_dr;
928 hepVector(2) = m_phi0;
929 hepVector(3) = m_kappa;
930 hepVector(4) = m_dz;
931 hepVector(5) = m_tanl;
932 a(hepVector);
933}
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37

Referenced by fitCircle(), fitHelix(), and updateCirclePar().

Member Data Documentation

◆ m_clearTrack

int HoughTrack::m_clearTrack =1
static

Definition at line 145 of file HoughTrack.h.

Referenced by clearMemory(), fitCircle(), fitHelix(), and HoughFinder::initialize().

◆ m_cut

TGraph * HoughTrack::m_cut = {NULL}
static

Definition at line 147 of file HoughTrack.h.

Referenced by findVHot(), and HoughFinder::initialize().

◆ m_useCgemInGlobalFit

int HoughTrack::m_useCgemInGlobalFit =3
static

Definition at line 146 of file HoughTrack.h.

Referenced by fitCircle(), fitHelix(), and HoughFinder::initialize().


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