CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Attributes | List of all members
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)
 

 
HoughTrack & operator= (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
 
TrkRecoTrk * getTrkRecoTrk ()
 
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 ()
 
int 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, double chiCut=7., int debug=0)
 
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 (int type=2)
 
void updateLayerStat ()
 
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)
 
HoughTrack * getMcTrack () const
 
void setMcTrack (HoughTrack *mcTrack)
 
void clearMemory ()
 
void setRecMdcHitVec (vector< RecMdcHit > &aRecMdcHitVec)
 
vector< RecMdcHit > * getRecMdcHitVec ()
 
void printRecMdcHitVec ()
 
bool nearStereoHits ()
 
void clearXVhitsCgem ()
 
int findHelixInCgem ()
 
- 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 HepPoint3D & center (void) const
 returns position of helix center(z = 0.);
 
const HepPoint3D & pivot (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 HepPoint3D & pivot (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
 
Helix & operator= (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
IMagneticFieldSvc * m_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 22 of file HoughTrack.cxx.

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

◆ HoughTrack() [2/5]

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

Definition at line 78 of file HoughTrack.cxx.

78 :Helix(position,momentum,charge)
79{
80 bFieldZ(-bFieldZ());
81 //Helix helix(position, momentum, charge); helix.bFieldZ(-helix.bFieldZ());
82 //m_helix = helix;
83 m_trkID = trkID;
84 m_charge = charge;
85 //m_angle = m_helix.center().phi();
86 //m_rho = 1./m_helix.radius();
87 m_angle = center().phi();
88 m_rho = 1./radius();
89 m_flag = 0;
90 if(m_angle<0)m_angle=+2*M_PI;
91 m_circleFitStat = 0;
92
93 //m_dr = m_helix.dr();
94 //m_phi0 = m_helix.phi0();
95 //m_kappa = m_helix.kappa();
96 //m_dz = m_helix.dz();
97 //m_tanl = m_helix.tanl();
98 m_dr = dr();
99 m_phi0 = phi0();
100 m_kappa = kappa();
101 m_dz = dz();
102 m_tanl = tanl();
103
104 //m_hot.clear();
105 m_trkRecoTrk = NULL;
106
107 m_vecHitPnt.clear();
108 m_vecHitResidual.clear();
109 m_vecHitChi2.clear();
110
111 m_dAngle = 0;
112 m_dRho = 0;
113 m_dTanl = 0;
114 m_dDz = 0;
115 m_chi2 = 0;
116 m_mcTrack = NULL;
117 m_cgemHitVector.clear();
118 m_nearStereoHits=false;
119}
momentum
**********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
Definition BesBdkRc/BesBdkRc-00-00-02/src/fortran/HepEvt.h:23
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:212
Helix::radius
double radius(void) const
returns radious of helix.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:224
Helix::dr
double dr(void) const
returns an element of parameters.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:236

◆ HoughTrack() [3/5]

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

Definition at line 121 of file HoughTrack.cxx.

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

◆ HoughTrack() [4/5]

HoughTrack::HoughTrack ( const HoughTrack & other)

Definition at line 229 of file HoughTrack.cxx.

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

◆ HoughTrack() [5/5]

HoughTrack::HoughTrack ( )

Definition at line 193 of file HoughTrack.cxx.

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

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 136 of file HoughTrack.h.

136{m_vecStereoHitPnt.push_back(aHitPnt);};

Referenced by HoughFinder::getMcParticleCol().

◆ calculateZ_S()

int HoughTrack::calculateZ_S ( HoughHit * hit)

Definition at line 783 of file HoughTrack.cxx.

784{
785 int nTangency(0);
786 if(hit->getFlag()==0)return nTangency;
787 //if(hit->getPairHit()==NULL)return nTangency;
788 //if(hit->getPairHit()->getHalfCircle()!=1)return nTangency;
789 double s(0),z(0);
790 double xc = center().x();
791 double yc = center().y();
792 double rTrack = radius();// signed FIXME
793 if(hit->getHitType()==HoughHit::CGEMHIT){
794 if(hit->getFlag()!=2)return nTangency;
795 hit->resetSZ();
796 //for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++)
797 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++){
798 if((*hotIt)->getFlag()!=0)continue;
799 if((*hotIt)->getHitType()==HoughHit::MDCHIT)continue;
800 if((*hotIt)->getUse()!=1) continue;
801 vector<int> vec_trkId = (*hotIt)->getTrkID();
802 vector<int>::iterator found_it = find(vec_trkId.begin(), vec_trkId.end(), getTrkID());
803 if(found_it==vec_trkId.end()) continue;
804 if((*hotIt)->getCgemCluster()->getclusterid() == hit->getCgemCluster()->getclusterflagb()){
805 HepPoint3D point = hit->getHitPosition();
806 s = flightArc(point);
807 z = point.z();
808 // double rHit = point.perp();
809 // s = flightArc(rHit);
810 // z = getZfrom ... // FIXME
811 pair<double,double> sz(s,z);
812 hit->addSZ(sz);
813 nTangency++;
814 //double s = flightArc(hit->getHitPosition());
815 int layer=hit->getLayer();
816 m_XVhits_cgem[layer].push_back(hit);// for self CGEM V-cluster selection
817 }
818 }
819 }else{
820 hit->resetSZ();
821 const MdcGeoWire* wire = hit->getMdcGeomSvc()->Wire(hit->getLayer(),hit->getWire());
822 int sz_version=1;
823
824 if(sz_version==0) {
825 double drift = hit->getDriftDist();
826 HepPoint3D westPoint = wire->Forward();
827 HepPoint3D eastPoint = wire->Backward();
828 double xEast = eastPoint.x()/10.0;
829 double xWest = westPoint.x()/10.0;
830 double yEast = eastPoint.y()/10.0;
831 double yWest = westPoint.y()/10.0;
832 double zEast = eastPoint.z()/10.0;
833 double zWest = westPoint.z()/10.0;
834 //cout<<"wast:x,y,z: "<<xWest<<", "<<yWest<<", "<<zWest<<endl;
835 //cout<<"east:x,y,z: "<<xEast<<", "<<yEast<<", "<<zEast<<endl;
836 //cout<<"Xc, Yc, Rc: "<<xc<<", "<<yc<<", "<<rTrack<<endl;
837 double west2east = sqrt((xEast-xWest)*(xEast-xWest)+(yEast-yWest)*(yEast-yWest));
838
839 double slope = (yEast-yWest)/(xEast-xWest);
840 double intercept = (yWest-slope*xWest+yEast-slope*xEast)/2;
841 double a = 1+slope*slope;
842 double b = -2*(xc+slope*yc-slope*intercept);
843 double c1 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+drift)*(rTrack+drift);
844 double c2 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack-drift)*(rTrack-drift);
845 //double c1 = intercept*(2*yc-intercept)+drift*(drift+rTrack);
846 //double c2 = intercept*(2*yc-intercept)+drift*(drift-rTrack);
847 double delta1 = (b*b-4*a*c1);
848 double delta2 = (b*b-4*a*c2);
849 //cout<<"a,b,c1,c2: "<<a<<", "<<b<<", "<<c1<<", "<<c2<<endl;
850 //cout<<"delta: "<<delta1<<", "<<delta2<<endl;
851 if(delta1>=0){
852 double x1 = (-b+sqrt(delta1))/(2*a);
853 double x2 = (-b-sqrt(delta1))/(2*a);
854 double y1 = slope*x1+intercept;
855 double y2 = slope*x2+intercept;
856 if((x1-xWest)*(x1-xEast)<0){
857 HepPoint3D point(x1,y1,0);
858 s = flightArc(point);
859 //s = flightArc(hit->getHitPosition());
860 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
861 z = zWest + l/west2east*fabs((zEast-zWest));
862 pair<double,double> sz(s,z);
863 hit->addSZ(sz);
864 nTangency++;
865 HepPoint3D position(x1,y1,z);
866 //hit->addPosition(position);
867 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
868 }
869 if((x2-xWest)*(x2-xEast)<0){
870 HepPoint3D point(x2,y2,0);
871 s = flightArc(point);
872 //s = flightArc(hit->getHitPosition());
873 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
874 z = zWest + l/west2east*fabs((zEast-zWest));
875 pair<double,double> sz(s,z);
876 hit->addSZ(sz);
877 nTangency++;
878 HepPoint3D position(x2,y2,z);
879 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
880 //hit->addPosition(position);
881 }
882 }
883
884 if( delta2>=0){
885 double x1 = (-b+sqrt(delta2))/(2*a);
886 double x2 = (-b-sqrt(delta2))/(2*a);
887 double y1 = slope*x1+intercept;
888 double y2 = slope*x2+intercept;
889 if((x1-xWest)*(x1-xEast)<0){
890 HepPoint3D point(x1,y1,0);
891 s = flightArc(point);
892 //s = flightArc(hit->getHitPosition());
893 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
894 z = zWest + l/west2east*fabs((zEast-zWest));
895 pair<double,double> sz(s,z);
896 hit->addSZ(sz);
897 nTangency++;
898 HepPoint3D position(x1,y1,z);
899 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
900 //hit->addPosition(position);
901 }
902 if((x2-xWest)*(x2-xEast)<0){
903 HepPoint3D point(x2,y2,0);
904 s = flightArc(point);
905 //s = flightArc(hit->getHitPosition());
906 //double l = zWest + sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
907 //z = l/west2east*fabs((zEast-zWest));
908 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
909 z = zWest + l/west2east*fabs((zEast-zWest));
910 pair<double,double> sz(s,z);
911 hit->addSZ(sz);
912 nTangency++;
913 HepPoint3D position(x2,y2,z);
914 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
915 //hit->addPosition(position);
916 }
917 }
918 // if(nTangency==0) { cout<<"delta1, delta2, " }
919 }
920 else if(sz_version==1) {
921 //vector<pair<double, double> > vec_sz;
922
923 // --- drift distance
924 double r_drift = hit->getDriftDist();
925
926 // --- wire line segment
927 HepPoint3D Xa=(wire->Forward())*0.1;
928 HepPoint3D Xb=(wire->Backward())*0.1;
929 HepPoint3D Xdelta = Xb - Xa;
930
931 // function
932 // A*lambda **2 + B*lambda + C ==0
933 double A = pow(Xdelta.x(),2 ) + pow(Xdelta.y(),2);
934 double B = 2 * Xdelta.x() * (Xa.x() - xc) + 2 * Xdelta.y() * (Xa.y() - yc);
935 double C_part = pow(Xa.x() - xc, 2) + pow(Xa.y() - yc, 2);
936 double C1 = C_part - pow(rTrack + r_drift,2);
937 double C2 = C_part - pow(rTrack - r_drift, 2);
938
939 // for each of the two C
940
941 //if(debug) cout<<"layer "<<myLayer<<": ";
942 for (int _ = 0; _ < 2; _++) {
943 double C = _ ? C2 : C1;
944 double Delta = B*B - 4 * A * C;
945
946 if (Delta < 0) continue;// require Delta >=0
947
948 for (int sign = -1; sign < 2; sign += 2) {// if Delta >0, we have 2 roots
949 double lambda = (-B + sign * sqrt(Delta)) / 2 / A;
950 if (lambda > 1.2 or lambda < -0.2) continue; // require local; +/-20% for safety
951
952 HepPoint3D Xhit = Xa + lambda * Xdelta;
953 double s = flightArc(Xhit);
954 //if(debug) cout<<"("<<s<<","<<Xhit.z()<<") ";
955 pair<double, double> sz(s, Xhit.z());
956 //vec_sz.push_back(sz);
957 hit->addSZ(sz);
958 nTangency++;
959 if (Delta == 0) break; // if Delta==0, run only once.
960 }
961 }
962 }// if(sz_version==1)
963 }// if MDC hits
964 return nTangency;
965}
c1
TCanvas * c1
Definition CgemLUTWriter.cxx:75
s
XmlRpcServer s
Definition HelloServer.cpp:11
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
_
#define _(A, B)
Definition cfortran.h:44
HoughHit::getMdcGeomSvc
MdcGeomSvc * getMdcGeomSvc() const
Definition HoughHit.h:65
HoughHit::getCgemCluster
const RecCgemCluster * getCgemCluster() const
Definition HoughHit.h:41
HoughHit::getDriftDist
double getDriftDist() const
Definition HoughHit.h:54
HoughHit::addSZ
void addSZ(S_Z sz)
Definition HoughHit.h:106
HoughHit::getHitPosition
HepPoint3D getHitPosition() const
Definition HoughHit.h:51
HoughHit::getHitType
HitType getHitType() const
Definition HoughHit.h:40
HoughHit::getFlag
int getFlag() const
Definition HoughHit.h:47
HoughHit::getLayer
int getLayer() const
Definition HoughHit.h:45
HoughHit::getWire
int getWire() const
Definition HoughHit.h:46
HoughHit::resetSZ
void resetSZ()
Definition HoughHit.h:97
HoughHit::CGEMHIT
@ CGEMHIT
Definition HoughHit.h:27
HoughTrack::getTrkID
int getTrkID() const
Definition HoughTrack.h:49
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:770
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition RecCgemCluster.h:73

Referenced by findVHot().

◆ clearHits()

void HoughTrack::clearHits ( )

Definition at line 1274 of file HoughTrack.cxx.

1275{
1276 //m_hot.clear();
1277 m_vecHitPnt.clear();
1278 m_vecHitResidual.clear();
1279 m_vecHitChi2.clear();
1280 m_map_lay_d.clear();
1281 m_map_lay_hit.clear();
1282 m_setLayer.clear();
1283 m_unused=0;
1284 m_untried=0;
1285 m_nHitFirstHalf=0;
1286 m_nHitSecondHalf=0;
1287 m_nHitUnused_FirstHalf=0;
1288 m_nHitUnused_SecondHalf=0;
1289
1290
1291}

◆ clearMemory()

void HoughTrack::clearMemory ( )

Definition at line 3051 of file HoughTrack.cxx.

3052{
3053 for(vector<CgemHitOnTrack*>::iterator iter = m_cgemHitVector.begin(); iter != m_cgemHitVector.end(); iter++){
3054 delete *iter;
3055 }
3056 if(m_clearTrack)delete m_trkRecoTrk;
3057}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
HoughTrack::m_clearTrack
static int m_clearTrack
Definition HoughTrack.h:146

◆ clearXVhitsCgem()

void HoughTrack::clearXVhitsCgem ( )

Definition at line 3059 of file HoughTrack.cxx.

3060{
3061 for(int i=0; i<3; i++) m_XVhits_cgem[i].clear();
3062}

◆ driftDistRes()

double HoughTrack::driftDistRes ( HoughHit * hit)

Definition at line 645 of file HoughTrack.cxx.

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

Referenced by findXHot().

◆ dropHitPnt()

void HoughTrack::dropHitPnt ( HoughHit * aHitPnt)

Definition at line 1942 of file HoughTrack.cxx.

1943{
1944 vector<HoughHit*>::iterator it0 = m_vecHitPnt.begin();
1945 vector<HoughHit*>::iterator result = find(m_vecHitPnt.begin(), m_vecHitPnt.end(), aHitPnt);
1946 if(result!=m_vecHitPnt.end()) {
1947 //cout<<"remove hit "<<(*result)->getHitID();
1948 m_vecHitPnt.erase(result);
1949 vector<double>::iterator itRes = m_vecHitResidual.begin()+(result-it0);
1950 if(itRes!=m_vecHitResidual.end()) {
1951 //cout<<" with residual "<<*itRes<<" from trk "<<getTrkID()<<endl;
1952 m_vecHitResidual.erase(itRes);
1953 }
1954 }
1955}

Referenced by dropRedundentCgemXHits().

◆ dropRedundentCgemVHits()

void HoughTrack::dropRedundentCgemVHits ( )

Definition at line 2032 of file HoughTrack.cxx.

2033{
2034 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
2035 vector<HoughHit*> hitPnt_cgem[3];
2036 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
2037 double res_cgem_min[3]={999., 999., 999.};
2038 vector<HoughHit*>::iterator iter = m_vecStereoHitPnt.begin();
2039 for(; iter!=m_vecStereoHitPnt.end(); iter++)
2040 {
2041 int iLayer = (*iter)->getLayer();
2042 if(iLayer<3) {
2043 hitPnt_cgem[iLayer].push_back(*iter);
2044 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
2045 double res = (*iter)->residual(this, positionOntrack, positionOnDetector);
2046 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
2047 res_cgem_min[iLayer] = res;
2048 hitPnt_cgem_resMin[iLayer] = *iter;
2049 }
2050 }
2051 }
2052 for(int iLayer=0; iLayer<3; iLayer++)
2053 {
2054 int nHits = hitPnt_cgem[iLayer].size();
2055 if(nHits>1)
2056 {
2057 HoughHit* hitToKeep = hitPnt_cgem_resMin[iLayer];
2058 for(int iHit=0; iHit<nHits; iHit++)
2059 {
2060 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
2061 if(aHoughHit==hitToKeep) continue;
2062 dropVHitPnt(aHoughHit);
2063 }
2064 }
2065 }
2066 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
2067}
HoughTrack::dropVHitPnt
void dropVHitPnt(HoughHit *aHitPnt)
Definition HoughTrack.cxx:1957

◆ dropRedundentCgemXHits()

int HoughTrack::dropRedundentCgemXHits ( )

Definition at line 1967 of file HoughTrack.cxx.

1968{
1969 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
1970 vector<HoughHit*> hitPnt_cgem[3];
1971 HoughHit* hitPnt_cgem_keep[3]={NULL, NULL, NULL};
1972 double res_cgem_keep[3]={999., 999., 999.};
1973 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
1974 double res_cgem_min[3]={999., 999., 999.};
1975 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1976 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1977 vector<double>::iterator itRes = m_vecHitResidual.begin();
1978 for(; iter!=m_vecHitPnt.end(); iter++)
1979 {
1980 int iLayer = (*iter)->getLayer();
1981 if(iLayer<3) {
1982 hitPnt_cgem[iLayer].push_back(*iter);
1983 int idx = iter-iter0;
1984 double res = m_vecHitResidual[idx];
1985 vector<double> vecRes = (*iter)->getVecResid();
1986 int nTrk = vecRes.size();
1987 if(nTrk==1&&fabs(res)<fabs(res_cgem_keep[iLayer])) {
1988 res_cgem_keep[iLayer] = res;
1989 hitPnt_cgem_keep[iLayer] = *iter;
1990 }
1991 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
1992 res_cgem_min[iLayer] = res;
1993 hitPnt_cgem_resMin[iLayer] = *iter;
1994 }
1995 }
1996 }
1997 int nDrop=0;
1998 for(int iLayer=0; iLayer<3; iLayer++)
1999 {
2000 int nHits = hitPnt_cgem[iLayer].size();
2001 if(nHits>1)
2002 {
2003 HoughHit* hitToKeep;
2004 if(hitPnt_cgem_keep[iLayer]!=NULL) hitToKeep = hitPnt_cgem_keep[iLayer];
2005 else hitToKeep = hitPnt_cgem_resMin[iLayer];
2006 //cout<<"HoughTrack::dropRedundentCgemXHits: "<<endl;
2007 //cout<<" keep hit "<<hitToKeep->getHitID()
2008 //<<", layer "<<hitToKeep->getLayer()
2009 //<<" at ("<<hitToKeep->getHitPosition().x()
2010 //<<", "<<hitToKeep->getHitPosition().y()<<")"<<endl;
2011 for(int iHit=0; iHit<nHits; iHit++)
2012 {
2013 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
2014 if(aHoughHit==hitToKeep) continue;
2015 //cout<<" remove hit "<<aHoughHit->getHitID()
2016 //<<", layer "<<aHoughHit->getLayer()
2017 //<<" at ("<<aHoughHit->getHitPosition().x()
2018 //<<", "<<aHoughHit->getHitPosition().y()<<")"<<endl;
2019 dropHitPnt(aHoughHit);
2020 nDrop++;
2021 //aHoughHit->dropTrkID(m_trkID);
2022 int used = aHoughHit->getUse();
2023 if(used==0) m_untried--;
2024 if(used==0||used==-1) m_unused--;// 0: unused, 1: used, -1: tried, but not used
2025 }
2026 }
2027 }
2028 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
2029 return nDrop;
2030}
HoughHit::getUse
int getUse() const
Definition HoughHit.h:58
HoughTrack::dropHitPnt
void dropHitPnt(HoughHit *aHitPnt)
Definition HoughTrack.cxx:1942

◆ dropVHitPnt()

void HoughTrack::dropVHitPnt ( HoughHit * aHitPnt)

Definition at line 1957 of file HoughTrack.cxx.

1958{
1959 vector<HoughHit*>::iterator it0 = m_vecStereoHitPnt.begin();
1960 vector<HoughHit*>::iterator result = find(m_vecStereoHitPnt.begin(), m_vecStereoHitPnt.end(), aHitPnt);
1961 if(result!=m_vecStereoHitPnt.end()) {
1962 //cout<<"remove hit "<<(*result)->getHitID();
1963 m_vecStereoHitPnt.erase(result);
1964 }
1965}

Referenced by dropRedundentCgemVHits().

◆ findHelixInCgem()

int HoughTrack::findHelixInCgem ( )

Definition at line 3064 of file HoughTrack.cxx.

3065{
3066}

◆ findVHot()

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

Definition at line 2083 of file HoughTrack.cxx.

2084{
2085 bool printFlag(false); //printFlag=true;
2086 int nHot(0);
2087 m_vecStereoHitPnt.clear();
2088 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end(); hitIt++){
2089 /*
2090 if(0 == flag){
2091 if((*hitIt)->getFlag() != 0) continue;
2092 if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
2093 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2094 //if((*hitIt)->getPairHit()->getHalfCircle()!=1)continue;//FIXME
2095 double cut1, cut2;
2096 int iLayer = (*hitIt)->getLayer();
2097 int used = (*hitIt)->getUse();
2098 XhitCutWindow(m_rho, iLayer, charge, cut1, cut2);
2099 double res = driftDistRes(*hitIt);
2100 map<int,double>::iterator it_map;
2101 it_map=m_map_lay_d.find(iLayer);
2102 if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
2103 {
2104 m_map_lay_d[iLayer]=res;
2105 m_map_lay_hit[iLayer]=(*hitIt);
2106 }
2107 if(res>cut1&&res<cut2)
2108 {
2109 //cout<<" selected!";
2110 //if(iLayer>3)
2111 m_vecHitPnt.push_back(*hitIt);
2112 m_vecHitResidual.push_back(res);
2113 //m_hot.push_back((*(*hitIt)));
2114 m_chi2 =+ res*res;
2115 m_setLayer.insert(iLayer);
2116 if(used==0) m_untried++;
2117 if(used==0||used==-1) m_unused++;
2118 if(judgeHalf(*hitIt)>0) {
2119 m_nHitFirstHalf++;
2120 if(used<=0) m_nHitUnused_FirstHalf++;
2121 }
2122 else {
2123 m_nHitSecondHalf++;
2124 if(used<=0) m_nHitUnused_SecondHalf++;
2125 }
2126 nHot++;
2127 }// within window
2128 //cout<<endl;
2129 }// X-view
2130 // */
2131 //else {
2132 int layer = (*hitIt)->getLayer();
2133 int wire = (*hitIt)->getWire();
2134 if((*hitIt)->getFlag() == 0) continue;
2135 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
2136 if(calculateZ_S(*hitIt)==0){
2137 continue;
2138 }
2139 }
2140 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2141 //if((*hitIt)->getPairHit()->getHalfCircle()>1)continue;//FIXME
2142 double Pt = fabs(pt());
2143 double cut[2] = {0};
2144 if(layer<3){
2145 if(m_cut[0][layer+3]!=NULL)cut[0]=m_cut[0][layer+3]->Eval(Pt);
2146 if(m_cut[1][layer+3]!=NULL)cut[1]=m_cut[1][layer+3]->Eval(Pt);
2147 }else{
2148 if(m_cut[0][layer]!=NULL)cut[0]=m_cut[0][layer]->Eval(Pt);
2149 if(m_cut[1][layer]!=NULL)cut[1]=m_cut[1][layer]->Eval(Pt);
2150 }
2151 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
2152 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
2153 if(judgeCharge(positionOntrack.x(),positionOntrack.y())*charge<0) continue;// charge requirement
2154 bool isNewHot(true);
2155 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++){
2156 if((*hitIt)->getHitID() == (*hotIt)->getHitID()){
2157 isNewHot = false;
2158 break;
2159 }
2160 }
2161 if(!isNewHot){
2162 //cout<<" OK!";
2163 //cout<<endl;
2164 continue;
2165 }
2166 if(printFlag)(*hitIt)->print();
2167 if(printFlag)cout<<"res:"<<setw(12)<<res;
2168 if(printFlag)cout<<" win: "<<setw(5)<<cut[0]<<" ~ "<<setw(5)<<cut[1];
2169 double ext=cutFactor;
2170 if(ext*cut[0]<res&&res<ext*cut[1]){
2171 if((*hitIt)->getLayer()>=maxLayer)
2172 {
2173 //if(printFlag) cout<<;
2174 continue;
2175 }
2176 m_vecStereoHitPnt.push_back(*hitIt);
2177 if(printFlag)cout<<" selected!";
2178 nHot++;
2179 /*
2180 if(maxLayer<20){
2181 if((*hitIt)->getLayer()>=maxLayer)continue;
2182 //HepPoint3D szz(minS,zTrk,zHit);
2183 //(*hitIt)->setHitPosition(szz);
2184 m_vecStereoHitPnt.push_back(*hitIt);
2185 //m_vecStereoHitRes.push_back(minRes);
2186 nHot++;
2187 //cout<<" OK!";
2188 //(*hitIt)->print();
2189 }
2190 else{
2191 //if((*hitIt)->getLayer()!=maxLayer)continue;
2192 if((*hitIt)->getLayer()>=maxLayer)continue;
2193 //HepPoint3D szz(minS,zTrk,zHit);
2194 //(*hitIt)->setHitPosition(szz);
2195 m_vecStereoHitPnt.push_back(*hitIt);
2196 //m_vecStereoHitRes.push_back(minRes);
2197 nHot++;
2198 //cout<<" OK!";
2199 //(*hitIt)->print();
2200 }
2201 // */
2202 }
2203 if(printFlag)cout<<endl;
2204 //if(layer<3)cout<<endl;
2205 //}// V-view
2206 }// loop hits
2207 return nHot;
2208}
cut
*********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
HoughTrack::judgeCharge
int judgeCharge(HoughHit *hit)
Definition HoughTrack.cxx:618
HoughTrack::m_cut
static TGraph * m_cut[2][43]
Definition HoughTrack.h:18
HoughTrack::calculateZ_S
int calculateZ_S(HoughHit *hit)
Definition HoughTrack.cxx:783

◆ findXHot()

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

Definition at line 490 of file HoughTrack.cxx.

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

◆ fitCircle() [1/2]

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

Definition at line 2266 of file HoughTrack.cxx.

2267{
2268 bool debug=false; //debug=true;
2269 m_circleFitStat=0;
2270 double d0 = -m_dr;
2271 double phi0 = m_phi0+M_PI/2;
2272 while(phi0>M_PI)phi0-=2*M_PI;
2273 while(phi0<-M_PI)phi0+=2*M_PI;
2274 double omega = m_kappa/fabs(alpha());
2275 //double z0 = m_dz;
2276 //double tanl = m_tanl;
2277 double z0 = 0;
2278 double tanl = 0;
2279 //cout<<"hough params:"<<endl;
2280 //cout
2281 //<<setw(15)<<d0
2282 //<<setw(15)<<phi0
2283 //<<setw(15)<<omega
2284 //;
2285 //cout<<endl;
2286 TrkExchangePar circleTrk(d0,phi0,omega,z0,tanl);
2287 TrkCircleMaker circlefactory;
2288 double chisum =0.;
2289 TrkRecoTrk* trkRecoTrk = circlefactory.makeTrack(circleTrk, chisum, *trkContextEv, bunchT0);
2290 bool permitFlips = true;
2291 bool lPickHits = true;
2292 circlefactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2293 TrkHitList* trkHitList = trkRecoTrk->hits();
2294
2295 //---convert hits
2296 vector<MdcHit*> mdcHitVector;
2297 mdcHitVector.clear();
2298 vector<CgemHitOnTrack*> cgemHitVector;
2299 cgemHitVector.clear();
2300
2301 for(vector<HoughHit*>::iterator iter = m_vecHitPnt.begin(); iter != m_vecHitPnt.end(); iter++){
2302 HoughHit* hitIt = (*iter);
2303 if(hitIt->getFlag()!=0)continue;
2304 //hitIt->print();
2305 if(hitIt->getHitType()==HoughHit::MDCHIT || hitIt->getHitType()==HoughHit::MDCMCHIT){
2306 const MdcDigi* mdcDigi = hitIt->getDigi();
2307 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2308 mdcHitVector.push_back(mdcHit);
2309 mdcHit->setCalibSvc(hitIt->getMdcCalibFunSvc());
2310 mdcHit->setCountPropTime(true);
2311 //if(mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;//FIXME
2312 int newAmbig = 0;
2313 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2314 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2315 HepPoint3D midPoint = hitIt->getHitPosition();
2316 double fltLength = flightLength(midPoint);
2317 //double fltLength = flightLength(hitIt->getHitPosition());
2318 mdcHitOnTrack->setFltLen(fltLength);
2319 trkHitList->appendHot(mdcHitOnTrack);
2320 }else if(hitIt->getHitType()==HoughHit::CGEMHIT || hitIt->getHitType()==HoughHit::CGEMMCHIT){
2321 if(m_useCgemInGlobalFit<1)continue;
2322 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*(hitIt->getCgemCluster()),bunchT0*1.e9);
2323 cgemHitVector.push_back(cgemHitOnTrack);
2324 cgemHitOnTrack->setCgemGeomSvc(hitIt->getCgemGeomSvc());
2325 cgemHitOnTrack->setCgemCalibFunSvc(hitIt->getCgemCalibFunSvc());
2326 trkHitList->appendHot(cgemHitOnTrack);
2327 }
2328 }
2329 //---fitting
2330 TrkHitList* qhits = trkRecoTrk->hits();
2331 //cout<<"num of qhits to fit = "<<qhits->nHit()<<endl;
2332 //cout<<"before TrkHitList->fit()"<<endl;
2333 TrkErrCode err=qhits->fit();
2334 //cout<<"after TrkHitList->fit()"<<endl;
2335
2336 if(err.success()){
2337 //cout<<" circle fits well "<<endl;
2338 m_circleFitStat=1;
2339 const TrkFit* fitResult = trkRecoTrk->fitResult();
2340 TrkExchangePar par = fitResult->helix(0.);
2341 //if(fabs(m_kappa-par.omega()*fabs(alpha()))/fabs(m_kappa)>0.3)cout<<"fit divergent"<<endl;
2342 m_dr = -par.d0();
2343 m_phi0 = par.phi0()-M_PI/2;
2344 m_kappa = par.omega()*fabs(alpha());
2345 m_dz = par.z0();
2346 m_tanl = par.tanDip();
2347
2348 //cout<<"phi0 = "<<m_phi0<<endl;
2349 while(m_phi0>2.*M_PI) m_phi0-=2.*M_PI;
2350 while(m_phi0<0.) m_phi0+=2.*M_PI;
2351 //cout<<"phi0 = "<<m_phi0<<endl;
2352 //HepVector hepVector(5,0);
2353 //hepVector(1) = m_dr;
2354 //hepVector(2) = m_phi0;
2355 //hepVector(3) = m_kappa;
2356 //hepVector(4) = m_dz;
2357 //hepVector(5) = m_tanl;
2358 //a(hepVector);
2359 updateHelix();
2360 //cout<<"update Helix"<<endl;
2361
2362 //m_trkRecoTrk = trkRecoTrk;
2363 m_chi2 = fitResult->chisq();
2364 if(debug) cout<<"chi2 obtained "<<m_chi2<<endl;
2365 //cout<<"fitting params:"<<endl;
2366 //cout
2367 //<<setw(15)<<par.d0()
2368 //<<setw(15)<<par.phi0()
2369 //<<setw(15)<<par.omega()
2370 //;
2371 //cout<<endl;
2372
2373 //m_vecMdcDigiPnt.clear();
2374 //m_vecCgemClusterPnt.clear();
2375
2376 int nHotKept = 0;
2377 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2378 {
2379 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2380 {
2381 const MdcRecoHitOnTrack* recoHot;
2382 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2383 if(!(recoHot->isActive())) continue;
2384 nHotKept++;
2385 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2386 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2387 }
2388 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2389 {
2390 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2391 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2392 int clusterid = recCgemCluster->getclusterid();
2393 int stat = recoHot->isActive();
2394 if(stat==0) continue;//
2395 nHotKept++;
2396 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2397 }
2398 }
2399 if(debug) cout<<nHotKept<<" hits kept after TrkFitter"<<endl;
2400 }
2401 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2402 delete *iter;
2403 }
2404 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2405 delete *iter;
2406 }
2407 if(m_clearTrack)delete trkRecoTrk;
2408 return err;
2409}
CgemHitOnTrack::setCgemCalibFunSvc
void setCgemCalibFunSvc(const CgemCalibFunSvc *svc)
Definition CgemHitOnTrack.h:115
CgemHitOnTrack::baseHit
const RecCgemCluster * baseHit() const
Definition CgemHitOnTrack.h:119
CgemHitOnTrack::setCgemGeomSvc
void setCgemGeomSvc(const CgemGeomSvc *svc)
Definition CgemHitOnTrack.h:114
HoughHit::getDigi
const MdcDigi * getDigi() const
Definition HoughHit.h:42
HoughHit::getMdcCalibFunSvc
MdcCalibFunSvc * getMdcCalibFunSvc() const
Definition HoughHit.h:66
HoughHit::getCgemGeomSvc
CgemGeomSvc * getCgemGeomSvc() const
Definition HoughHit.h:67
HoughHit::getCgemCalibFunSvc
CgemCalibFunSvc * getCgemCalibFunSvc() const
Definition HoughHit.h:68
HoughTrack::m_useCgemInGlobalFit
static int m_useCgemInGlobalFit
Definition HoughTrack.h:147
HoughTrack::updateHelix
void updateHelix()
Definition HoughTrack.cxx:987
MdcHit::setCalibSvc
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition MdcHit.cxx:136
MdcHit::setCountPropTime
void setCountPropTime(const bool count)
Definition MdcHit.h:86
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition RecCgemCluster.h:54
TrkAbsFit::chisq
virtual double chisq() const =0
TrkErrCode::success
int success() const
Definition TrkErrCode.h:62
TrkExchangePar::phi0
double phi0() const
Definition TrkExchangePar.h:43
TrkExchangePar::omega
double omega() const
Definition TrkExchangePar.h:44
TrkExchangePar::z0
double z0() const
Definition TrkExchangePar.h:45
TrkExchangePar::d0
double d0() const
Definition TrkExchangePar.h:42
TrkExchangePar::tanDip
double tanDip() const
Definition TrkExchangePar.h:46
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitList::fit
TrkErrCode fit()
Definition TrkHitList.cxx:59
TrkHitList::begin
hot_iterator begin() const
Definition TrkHitList.h:45
TrkHitList::appendHot
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition TrkHitList.cxx:91
TrkHitList::end
hot_iterator end() const
Definition TrkHitList.h:46
TrkHitOnTrk::setFltLen
void setFltLen(double f)
Definition TrkHitOnTrk.h:147
TrkHitOnTrk::isActive
bool isActive() const
Definition TrkHitOnTrk.h:200
TrkRecoTrk::hits
TrkHitList * hits()
Definition TrkRecoTrk.h:107
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition TrkRecoTrk.cxx:387
TrkSimpleMaker::setFlipAndDrop
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkSimpleMaker::makeTrack
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const

◆ fitCircle() [2/2]

int HoughTrack::fitCircle ( DotsHelixFitter * fitter,
double bunchT0,
double chiCut = 7.,
int debug = 0 )

Definition at line 2411 of file HoughTrack.cxx.

2412{
2413 //int debug=false; //debug=true;
2414 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2415 ///*
2416 m_vecMdcDigiPnt.clear();
2417 m_vecCgemClusterPnt.clear();
2418 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2419 {
2420 if((*hotIt)->getFlag()!=0)continue;
2421 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2422 {
2423 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2424 m_vecMdcDigiPnt.push_back(aMdcDigi);
2425 }
2426 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2427 {
2428 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2429 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2430 }
2431 }
2432 //*/
2433
2434 // --- set
2435 HepPoint3D pivot(0,0,0);
2436 HepVector a(5,0);
2437 HepSymMatrix Ea(5,0);
2438
2439 // --- initial circle parameters from Hough transform
2440 double phi0_Hough = (m_rho>0)? m_angle:m_angle+M_PI;
2441 phi0_Hough = (m_charge<0)? phi0_Hough:phi0_Hough+M_PI;
2442 if(phi0_Hough>2*M_PI)phi0_Hough-=2*M_PI;
2443 if(phi0_Hough<0) phi0_Hough+=2*M_PI;
2444 double kappa_Hough = fabs(m_alpha*m_rho)*m_charge;
2445 a(2) = phi0_Hough;
2446 a(3) = kappa_Hough;
2447
2448 // --- from fit?
2449 //a(1) = m_dr;
2450 //a(2) = m_phi0;
2451 //a(3) = m_kappa;
2452 //a(5) = 0.00001;
2453
2454 KalmanFit::Helix ini_helix(pivot,a,Ea);
2455 if(debug==1) cout<<"circle ini_helix: "<<ini_helix.a()<<endl;
2456 fitter->setInitialHelix(ini_helix);
2457 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2458 fitter->setCgemClusters(m_vecCgemClusterPnt);
2459 if(debug==1) cout<<"finish DotsHelixFitter setting"<<endl;
2460 if(debug==5) cout<<"initial pt = "<<ini_helix.pt()<<endl;
2461 int fitFlag = 0;
2462
2463 BesTimer timer("circle_fit");
2464
2465 // --- circle Taubin fit without initial helix
2466 //clock_t begin = clock();
2467 timer.start();
2468 vector<double> a_Taubin;
2469 while(1)
2470 {
2471 a_Taubin = fitter->calculateCirclePar_Taubin(false);
2472 if(a_Taubin[3]<0) break;
2473 if( fitter->deactiveHits(chiCut) )
2474 continue;
2475 else break;
2476 }
2477 //clock_t end = clock();
2478 //double time_Taubin_0 = (double)(end-begin) / CLOCKS_PER_SEC;
2479 timer.stop();
2480 double time_Taubin_0 = timer.elapsed();
2481 if(debug==5) cout<<setw(45)<<"circle Taubin fit without initial helix:"<<" pt = "<<1.0/a_Taubin[2]<<", dr="<<a_Taubin[0]<<", phi0="<<a_Taubin[1]<<", nHitKept="<<fitter->getNActiveHits()<<", chi2="<<fitter->getChi2()<<", converged="<<a_Taubin[3]<<", time="<<time_Taubin_0<<" ms"<<endl;
2482 int NHitsKept_0=fitter->getNActiveHits();
2483 double chi2_0=fitter->getChi2();
2484
2485 // --- circle Taubin fit with Hough initial helix
2486 //begin = clock();
2487 timer.start();
2488 vector<double> a_Taubin_2;
2489 fitter->resetChi2FitFlag();
2490 fitter->setInitialHelix(ini_helix);
2491 while(1)
2492 {
2493 a_Taubin_2 = fitter->calculateCirclePar_Taubin(true);
2494 if(a_Taubin[3]<0) break;
2495 if( fitter->deactiveHits(chiCut) )
2496 continue;
2497 else break;
2498 }
2499 //end = clock();
2500 //double time_Taubin_1 = (double)(end-begin) / CLOCKS_PER_SEC;
2501 timer.stop();
2502 double time_Taubin_1 = timer.elapsed();
2503 if(debug==5) cout<<setw(45)<<"circle Taubin fit with Hough initial helix:"<<" pt = "<<1.0/a_Taubin_2[2]<<", dr="<<a_Taubin_2[0]<<", phi0="<<a_Taubin_2[1]<<", nHitKept="<<fitter->getNActiveHits()<<", chi2="<<fitter->getChi2()<<", converged="<<a_Taubin_2[3]<<", time="<<time_Taubin_1<<" ms"<<endl;
2504
2505 // --- select better Taubin result
2506 int selTaubin=0;
2507 vector<double> a_Taubin_better=a_Taubin_2;
2508 if(a_Taubin_2[3]>0) {
2509 selTaubin=2;
2510 if(a_Taubin[3]>0) {
2511 if(a_Taubin[3]>0&&NHitsKept_0==fitter->getNActiveHits()) {
2512 if(chi2_0<fitter->getChi2()) {
2513 a_Taubin_better=a_Taubin;
2514 selTaubin=1;
2515 }
2516 }
2517 else if(NHitsKept_0>fitter->getNActiveHits())
2518 {
2519 //if(fitter->getChi2()-chi2_0<(fitter->getNActiveHits()-NHitsKept_0)*)
2520 a_Taubin_better=a_Taubin;
2521 selTaubin=1;
2522 }
2523 }
2524 }
2525 else if(a_Taubin[3]>0) {
2526 a_Taubin_better=a_Taubin;
2527 selTaubin=1;
2528 }
2529 a(1)=a_Taubin_better[0];
2530 a(2)=a_Taubin_better[1];
2531 a(3)=a_Taubin_better[2];
2532 KalmanFit::Helix ini_helix_Taubin(pivot,a,Ea);
2533 if(debug==5)cout<<"better Taubin "<<selTaubin<<endl;
2534
2535 // --- Taubin fit with MDC axial hits only (for MDC stereo hits selection) // to be done, FIXME
2536
2537 // --- circle fit with Hough helix
2538 //begin = clock();
2539 timer.start();
2540 fitter->setInitialHelix(ini_helix);
2541 //fitter->fitCircleOnly();
2542 fitter->resetChi2FitFlag();
2543 fitter->fitModeCircle();
2544 while(1)
2545 {
2546 fitFlag = fitter->calculateNewHelix();
2547 if(fitFlag==0)
2548 {
2549 if(debug==1) {
2550 HepVector a_new = fitter->getHelix();
2551 cout<<"circle helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2552 }
2553 //if( fitter->getChi2() > chiCut*(fitter->getNActiveHits()) && fitter->deactiveHits(chiCut))
2554 if( fitter->deactiveHits(chiCut) )
2555 {
2556 continue;
2557 if(debug==1) cout<<"too large chi2, drop one worst DC hit!"<<endl;
2558 }
2559 else break;
2560 }
2561 else
2562 {
2563 if(debug==1) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2564 break;
2565 }
2566 }
2567 HepVector a_new = fitter->getHelix();
2568 //end = clock();
2569 //double time_LS = (double)(end-begin) / CLOCKS_PER_SEC;
2570 timer.stop();
2571 double time_LS = timer.elapsed();
2572 if(debug==5) {
2573 cout<<setw(45)<<"circle fit with Hough initial helix:"<<" pt = "<<1.0/a_new[2]<<", dr="<<a_new[0]<<", phi0="<<a_new[1]<<", nHitKept="<<fitter->getNActiveHits()<<", chi2="<<fitter->getChi2()<<", time="<<time_LS<<" ms";
2574 if(fitFlag!=0) cout<<" (failed) ";
2575 cout<<endl;
2576 }
2577
2578 // --- circle fit with Taubin helix
2579 timer.start();
2580 if(selTaubin>0) {
2581 fitter->setInitialHelix(ini_helix_Taubin);
2582 fitter->resetChi2FitFlag();
2583 fitter->fitModeCircle();
2584 while(1)
2585 {
2586 fitFlag = fitter->calculateNewHelix();
2587 if(fitFlag==0)
2588 {
2589 if(debug==1) {
2590 HepVector a_new = fitter->getHelix();
2591 cout<<"circle helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2592 }
2593 //if( fitter->getChi2() > chiCut*(fitter->getNActiveHits()) && fitter->deactiveHits(chiCut))
2594 if( fitter->deactiveHits(chiCut) )
2595 {
2596 continue;
2597 if(debug==1) cout<<"too large chi2, drop one worst DC hit!"<<endl;
2598 }
2599 else break;
2600 }
2601 else
2602 {
2603 if(debug==1) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2604 break;
2605 }
2606 }
2607 HepVector a_new_2 = fitter->getHelix();
2608 timer.stop();
2609 double time_LS_2 = timer.elapsed();
2610 if(debug==5) {
2611 cout<<setw(45)<<"circle fit with Taubin initial helix:"<<" pt = "<<1.0/a_new_2[2]<<", dr="<<a_new_2[0]<<", phi0="<<a_new_2[1]<<", nHitKept="<<fitter->getNActiveHits()<<", chi2="<<fitter->getChi2()<<", time="<<time_LS_2<<" ms";
2612 if(fitFlag!=0) cout<<" (failed) ";
2613 cout<<endl;
2614 }
2615 }
2616
2617 return fitFlag;
2618
2619}
DotsHelixFitter::setInitialHelix
void setInitialHelix(KalmanFit::Helix aHelix)
Definition DotsHelixFitter.cxx:881
DotsHelixFitter::setDChits
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
Definition DotsHelixFitter.cxx:160
DotsHelixFitter::setCgemClusters
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
Definition DotsHelixFitter.cxx:214
DotsHelixFitter::deactiveHits
int deactiveHits(double chi_cut=10, int nMax=1, int layerMin=0, int layerMax=50)
Definition DotsHelixFitter.cxx:736
DotsHelixFitter::getHelix
HepVector getHelix()
Definition DotsHelixFitter.h:59
DotsHelixFitter::calculateCirclePar_Taubin
vector< double > calculateCirclePar_Taubin(bool useIniHelix=false, int layerMin=-1, int layerMax=50, bool useExtraPoints=false, int debug=0)
Definition DotsHelixFitter.cxx:1425
HoughTrack::getChi2
double getChi2() const
Definition HoughTrack.h:58

◆ fitHelix() [1/3]

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

Definition at line 2762 of file HoughTrack.cxx.

2763{
2764 // --- make TrkRecoTrk and trkHitList
2765 int nHot(0);
2766 double d0 = -m_dr;
2767 double phi0 = m_phi0+M_PI/2;
2768 double omega = m_kappa/fabs(alpha());
2769 double z0 = m_dz;
2770 double tanl = m_tanl;
2771 while(phi0>M_PI)phi0-=2*M_PI;
2772 while(phi0<-M_PI)phi0+=2*M_PI;
2773
2774 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2775 TrkHelixMaker helixfactory;
2776 double chisum =0.;
2777 long idnum(0);
2778 TrkRecoTrk* trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, bfield, bunchT0);
2779 bool permitFlips = true;
2780 bool lPickHits = true;
2781 helixfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2782 TrkHitList* trkHitList = trkRecoTrk->hits();
2783
2784 //--- convert hits and put them into trkHitList
2785 //cout<<maxLayer<<endl;
2786 vector<MdcHit*> mdcHitVector;
2787 mdcHitVector.clear();
2788 vector<CgemHitOnTrack*> cgemHitVector;
2789 cgemHitVector.clear();
2790
2791 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2792 //(*hitIt)->print();
2793 //cout<<endl;
2794 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2795 if((*hitIt)->getFlag()!=0&&(*hitIt)->getLayer()>maxLayer)continue;
2796 bool sameHit(false);
2797 const TrkHitList* aList = trkRecoTrk->hits();
2798 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2799 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2800 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2801 int layerId = recoHot->mdcHit()->layernumber();
2802 int wireId = recoHot->mdcHit()->wirenumber();
2803 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2804 }
2805 }
2806 if(sameHit)continue;
2807
2808 const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2809 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2810 mdcHitVector.push_back(mdcHit);
2811
2812 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2813 mdcHit->setCountPropTime(true);
2814 int newAmbig = 0;
2815 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2816 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2817 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2818 //double fltLength = flightLength(midPoint);
2819 double fltLength = flightArc(midPoint);
2820 //double fltLength = flightLength((*hitIt)->getHitPosition());
2821 mdcHitOnTrack->setFltLen(fltLength);
2822 trkHitList->appendHot(mdcHitOnTrack);
2823 nHot++;
2824 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2825 if(m_useCgemInGlobalFit<3)continue;
2826 bool sameHit(false);
2827 const TrkHitList* aList = trkRecoTrk->hits();
2828 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2829 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2830 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2831 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2832 int clusterid = recCgemCluster->getclusterid();
2833 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2834 }
2835 }
2836 if(sameHit)continue;
2837
2838 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2839 cgemHitVector.push_back(cgemHitOnTrack);
2840 if((*hitIt)->getFlag()==0)continue;
2841 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2842 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2843 trkHitList->appendHot(cgemHitOnTrack);
2844 nHot++;
2845 }
2846 }
2847 //cout<<trkHitList->nHit()<<endl;
2848 //---fitting
2849 TrkHitList* qhits = trkRecoTrk->hits();
2850 TrkErrCode err=qhits->fit();
2851 //err.print(std::cout);cout<<endl;cout<<endl;
2852
2853 //const TrkHitList* aList = trkRecoTrk->hits();
2854 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2855 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2856 //const MdcRecoHitOnTrack* recoHot;
2857 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2858 //int layerId = recoHot->mdcHit()->layernumber();
2859 //int wireId = recoHot->mdcHit()->wirenumber();
2860 //double res=999.,rese=999.;
2861 //if (recoHot->resid(res,rese,false)){
2862 //}else{}
2863 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2864 //}
2865 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2866 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2867 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2868 //int clusterid = recCgemCluster->getclusterid();
2869 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2870 //}
2871 //}
2872
2873 // --- if fit successes, update the helix parameters
2874 if(err.success()){
2875 const TrkFit* fitResult = trkRecoTrk->fitResult();
2876 TrkExchangePar par = fitResult->helix(0.);
2877 m_dr = -par.d0();
2878 m_phi0 = par.phi0()-M_PI/2;
2879 m_kappa = par.omega()*fabs(alpha());
2880 m_dz = par.z0();
2881 m_tanl = par.tanDip();
2882 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2883 while(m_phi0<0)m_phi0+=2*M_PI;
2884 updateHelix();
2885 m_chi2 = fitResult->chisq();
2886
2887 int nHotKept = 0;
2888 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2889 {
2890 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2891 {
2892 const MdcRecoHitOnTrack* recoHot;
2893 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2894 if(!(recoHot->isActive())) continue;
2895 nHotKept++;
2896 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2897 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2898 }
2899 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2900 {
2901 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2902 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2903 int clusterid = recCgemCluster->getclusterid();
2904 int stat = recoHot->isActive();
2905 if(stat==0) continue;//
2906 nHotKept++;
2907 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2908 }
2909 }
2910
2911 cout<<"chi2= "<<m_chi2
2912 <<", "<<nHotKept<<" hits kept after TrkFitter"
2913 <<", helix from TrkFitter = "<<a()
2914 <<endl;
2915 }
2916 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2917 delete *iter;
2918 }
2919 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2920 delete *iter;
2921 }
2922 if(m_clearTrack)delete trkRecoTrk;
2923 //cout<<"getId():"<<trkRecoTrk->id()<<endl;
2924 return err;
2925}
MdcHit::layernumber
unsigned layernumber() const
Definition MdcHit.h:61
MdcHit::wirenumber
unsigned wirenumber() const
Definition MdcHit.h:62
MdcRecoHitOnTrack::mdcHit
const MdcHit * mdcHit() const
Definition MdcRecoHitOnTrack.cxx:38

◆ fitHelix() [2/3]

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

Definition at line 2621 of file HoughTrack.cxx.

2622{
2623 int nHot(0);
2624 double d0 = -m_dr;
2625 double phi0 = m_phi0+M_PI/2;
2626 double omega = m_kappa/fabs(alpha());
2627 double z0 = m_dz;
2628 double tanl = m_tanl;
2629 while(phi0>M_PI)phi0-=2*M_PI;
2630 while(phi0<-M_PI)phi0+=2*M_PI;
2631
2632 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2633 TrkHelixMaker helixfactory;
2634 double chisum =0.;
2635 m_trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, *trkContextEv, bunchT0);
2636 bool permitFlips = true;
2637 bool lPickHits = true;
2638 helixfactory.setFlipAndDrop(*m_trkRecoTrk, permitFlips, lPickHits);
2639 TrkHitList* trkHitList = m_trkRecoTrk->hits();
2640
2641 //--- convert hits
2642 //cout<<maxLayer<<endl;
2643 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2644 //(*hitIt)->print();
2645 //cout<<endl;
2646 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2647 bool sameHit(false);
2648 const TrkHitList* aList = m_trkRecoTrk->hits();
2649 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2650 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2651 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2652 int layerId = recoHot->mdcHit()->layernumber();
2653 int wireId = recoHot->mdcHit()->wirenumber();
2654 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2655 }
2656 }
2657 if(sameHit)continue;
2658
2659 //const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2660 //MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2661 MdcHit *mdcHit;
2662 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
2663 for(;imdcHit != mdcHitCol.end();imdcHit++){
2664 if((*imdcHit)->mdcId() == (*hitIt)->getDigi()->identify()){
2665 mdcHit = *imdcHit;
2666 break;
2667 }
2668 }
2669
2670 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2671 mdcHit->setCountPropTime(true);
2672 int newAmbig = 0;
2673 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2674 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2675 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2676 //double fltLength = flightLength(midPoint);
2677 double fltLength = flightArc(midPoint);
2678 //double fltLength = flightLength((*hitIt)->getHitPosition());
2679 mdcHitOnTrack->setFltLen(fltLength);
2680 trkHitList->appendHot(mdcHitOnTrack);
2681 nHot++;
2682 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2683 if(m_useCgemInGlobalFit<2)continue;
2684 bool sameHit(false);
2685 const TrkHitList* aList = m_trkRecoTrk->hits();
2686 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2687 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2688 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2689 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2690 int clusterid = recCgemCluster->getclusterid();
2691 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2692 }
2693 }
2694 if(sameHit)continue;
2695
2696 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2697 m_cgemHitVector.push_back(cgemHitOnTrack);
2698 if((*hitIt)->getFlag()==0)continue;
2699 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2700 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2701 trkHitList->appendHot(cgemHitOnTrack);
2702 nHot++;
2703 }
2704 }
2705 //cout<<trkHitList->nHit()<<endl;
2706 //---fitting
2707 TrkHitList* qhits = m_trkRecoTrk->hits();
2708 TrkErrCode err=qhits->fit();
2709 //err.print(std::cout);
2710 //cout<<endl;
2711 //cout<<endl;
2712
2713 //const TrkHitList* aList = m_trkRecoTrk->hits();
2714 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2715 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2716 //const MdcRecoHitOnTrack* recoHot;
2717 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2718 //int layerId = recoHot->mdcHit()->layernumber();
2719 //int wireId = recoHot->mdcHit()->wirenumber();
2720 //double res=999.,rese=999.;
2721 //if (recoHot->resid(res,rese,false)){
2722 //}else{}
2723 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2724 //}
2725 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2726 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2727 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2728 //int clusterid = recCgemCluster->getclusterid();
2729 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2730 //}
2731 //}
2732
2733 if(err.success()){
2734 const TrkFit* fitResult = m_trkRecoTrk->fitResult();
2735 TrkExchangePar par = fitResult->helix(0.);
2736 m_dr = -par.d0();
2737 m_phi0 = par.phi0()-M_PI/2;
2738 m_kappa = par.omega()*fabs(alpha());
2739 m_dz = par.z0();
2740 m_tanl = par.tanDip();
2741 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2742 while(m_phi0<0)m_phi0+=2*M_PI;
2743 //HepVector hepVector(5,0);
2744 //a(hepVector);
2745 updateHelix();
2746
2747 m_chi2 = fitResult->chisq();
2748 //cout<<"fitting params:"<<endl;
2749 //cout
2750 //<<setw(15)<<par.d0()
2751 //<<setw(15)<<par.phi0()
2752 //<<setw(15)<<par.omega()
2753 //<<setw(15)<<par.z0()
2754 //<<setw(15)<<par.tanDip()
2755 //;
2756 //cout<<endl;
2757 }
2758 //cout<<"getId()="<<m_trkRecoTrk->id()<<endl;
2759 return err;
2760}

◆ fitHelix() [3/3]

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

Definition at line 2927 of file HoughTrack.cxx.

2928{
2929 bool debug=false; //debug=true;
2930 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2931 //cout<<"before fill"<<endl;
2932 m_vecMdcDigiPnt.clear();
2933 m_vecCgemClusterPnt.clear();
2934 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2935 {
2936 //if((*hotIt)->getFlag()!=0)continue;
2937 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2938 {
2939 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2940 m_vecMdcDigiPnt.push_back(aMdcDigi);
2941 }
2942 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2943 {
2944 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2945 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2946 }
2947 }
2948 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++)
2949 {
2950 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2951 {
2952 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2953 m_vecMdcDigiPnt.push_back(aMdcDigi);
2954 }
2955 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2956 {
2957 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2958 if(aRecCgemCluster->getflag()==2)
2959 {
2960 int clusterId = aRecCgemCluster->getclusterflage();
2961 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterColBegin+clusterId;
2962 //cout<<"V cluster ID = "<<setw(10)<<clusterId<<setw(10)<<(*cgemClusterIter)->getclusterid();
2963 if(clusterId==(*cgemClusterIter)->getclusterid())
2964 m_vecCgemClusterPnt.push_back(*cgemClusterIter);
2965 }
2966 }
2967 }
2968 //cout<<"end filling"<<endl;
2969
2970
2971 // --- set
2972 HepPoint3D pivot(0,0,0);
2973 HepVector a(5,0);
2974 HepSymMatrix aEa(5,0);
2975 a(1) = m_dr;
2976 a(2) = m_phi0;
2977 a(3) = m_kappa;
2978 a(4) = m_dz;
2979 a(5) = m_tanl;
2980 KalmanFit::Helix ini_helix(pivot,a,aEa);
2981 if(debug) cout<<"fitHelix ini_helix: "<<ini_helix.a()<<endl;
2982 fitter->setInitialHelix(ini_helix);
2983 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2984 fitter->setCgemClusters(m_vecCgemClusterPnt);
2985 int nMdcDigi = m_vecMdcDigiPnt.size();
2986
2987 // --- helix fit
2988 fitter->fitModeHelix();
2989 int fitFlag = 0;
2990 while(1)
2991 {
2992 fitFlag = fitter->calculateNewHelix();
2993 if(fitFlag==0)
2994 {
2995 HepVector a_new = fitter->getHelix();
2996 if(debug) {
2997 cout<<"helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2998 cout<<"nMdcDigi = "<<nMdcDigi<<endl;
2999 }
3000 //if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits())
3001 //if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits(43,--nMdcDigi))
3002 if( fitter->deactiveHits(sqrt(averageChi2cut),1) )
3003 {
3004 //if(debug) cout<<"too large chi2, drop the outmost DC hit!"<<endl;
3005 if(debug) cout<<"too large chi2, drop the worst DC hit!"<<endl;
3006 continue;
3007 }
3008 else
3009 {
3010 if(debug) cout<<"stop droping hit! "<<endl;
3011 break;
3012 }
3013 }
3014 else
3015 {
3016 if(debug) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
3017 break;
3018 }
3019 }// fit iterations
3020
3021 //int nReActiTot=0;
3022 //while(1)
3023 //{
3024 // if(fitFlag!=0) break;
3025 // int nReActi = fitter.activeHits(chiCut);
3026 // if(nReActi==0) break;
3027 // nReActiTot+=nReActi;
3028 // fitFlag = myDotsHelixFitter.calculateNewHelix();
3029 // if( fitter->deactiveHits(averageChi2cut,nReActi)<nReActi ) {
3030 // continue
3031 // }
3032 // else break;
3033 // nIter++;
3034 // //if(nIter>100) break;
3035 //}
3036
3037 if(fitFlag==0) //update dr ... and Ea
3038 {
3039 HepVector a = fitter->getHelix();
3040 m_dr = a[0];
3041 m_phi0 = a[1];
3042 m_kappa = a[2];
3043 m_dz = a[3];
3044 m_tanl = a[4];
3045 m_chi2 = fitter->getChi2();
3046 Ea(fitter->getEa());
3047 }
3048 return fitFlag;
3049}
DotsHelixFitter::getEa
const HepSymMatrix & getEa()
Definition DotsHelixFitter.h:60
RecCgemCluster::getflag
int getflag(void) const
Definition RecCgemCluster.h:58
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition RecCgemCluster.h:74

◆ 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;}

◆ getChi2()

double HoughTrack::getChi2 ( ) const
inline

Definition at line 58 of file HoughTrack.h.

58{return m_chi2;}

Referenced by fitCircle(), and 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 2214 of file HoughTrack.cxx.

2215{
2216 vector<HoughHit*> hotList;
2217 vector<HoughHit*> axialHot = getVecHitPnt();
2218 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
2219
2220 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2221 if(type==1)continue;
2222 if(type==2&&(**hotIt).getHitType()==HoughHit::CGEMHIT)continue;
2223 hotList.push_back(*hotIt);
2224 }
2225 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2226 if(type>0)hotList.push_back(*hotIt);
2227 }
2228
2229 /*
2230 // --- CGEM XV cluster
2231 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2232 if((**hotIt).getLayer()<3)hotList.push_back(*hotIt);
2233//cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2234}
2235
2236// --- MDC hits, layer 9-20
2237for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2238//if(m_fitFlag==1&&(**hotIt).getLayer()<3)hot.push_back(**hotIt);
2239if((**hotIt).getLayer()>3&&(**hotIt).getLayer()<20)hotList.push_back(*hotIt);
2240//cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2241}
2242
2243// --- MDC hits, layer 21-36
2244for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2245if((**hotIt).getLayer()>=20&&(**hotIt).getLayer()<36)hotList.push_back(*hotIt);
2246//int layer = (*hotIt)->getLayer();
2247//int wire = (*hotIt)->getWire();
2248//cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2249//cout<<(**hotIt).residual(&(*trkIT));
2250//cout<<endl;
2251//cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2252}
2253
2254// --- MDC hits, layer 37-43
2255for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2256if((**hotIt).getLayer()>=36)hotList.push_back(*hotIt);
2257//cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2258}
2259//cout<<"nHot: "<<hot.size()<<endl;
2260// */
2261
2262sortHot(hotList);
2263return hotList;
2264}
HoughTrack::sortHot
void sortHot(vector< HoughHit * > &hotList)
Definition HoughTrack.cxx:982
HoughTrack::getVecHitPnt
vector< HoughHit * > getVecHitPnt()
Definition HoughTrack.h:83
HoughTrack::getVecStereoHitPnt
vector< HoughHit * > getVecStereoHitPnt()
Definition HoughTrack.h:135

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 142 of file HoughTrack.h.

142{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 2070 of file HoughTrack.cxx.

2071{
2072 int nShared = 0;
2073 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
2074 for(; iter!=m_vecHitPnt.end(); iter++)
2075 {
2076 vector<double> vecRes = (*iter)->getVecResid();
2077 int nTrk = vecRes.size();
2078 if(nTrk>1) nShared++;
2079 }
2080 return nShared;
2081}

◆ 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 130 of file HoughTrack.h.

130{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 152 of file HoughTrack.h.

152{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 calculateZ_S(), HoughFinder::getMcParticleCol(), markUsedHot(), 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 134 of file HoughTrack.h.

134{return m_vecHitResidual;}

◆ getVecStereoHitPnt()

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

Definition at line 135 of file HoughTrack.h.

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

Referenced by getHotList(), and print().

◆ getVecStereoHitRes()

vector< double > HoughTrack::getVecStereoHitRes ( )
inline

Definition at line 137 of file HoughTrack.h.

137{return m_vecStereoHitRes;}

◆ isAGoodCircle()

bool HoughTrack::isAGoodCircle ( )

Definition at line 1588 of file HoughTrack.cxx.

1589{
1590 bool debug = false; //debug=true;
1591 int NtotLayer = m_setLayer.size();
1592
1593 std::set<int>::iterator it = m_setLayer.begin();
1594 int lastLayer = *it;
1595 // check N_CGEM_layer, N_ODC1_nearStereo, N_ODC2_nearStereo
1596 int nCgemLayer(0), N_ODC1_nearStereo(0), N_ODC2_nearStereo(0);
1597 for(; it!=m_setLayer.end(); it++)
1598 {
1599 if((*it)<3) nCgemLayer++;
1600 else if((*it)>=16&&(*it)<=19) N_ODC1_nearStereo++;
1601 else if((*it)>=36&&(*it)<=39) N_ODC2_nearStereo++;
1602 }
1603 bool enoughVHits;
1604 //enoughVHits = nCgemLayer==3 ;
1605 enoughVHits = nCgemLayer>=2
1606 || (nCgemLayer>0 && (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2))
1607 //|| (N_ODC1_nearStereo>=2 && N_ODC2_nearStereo>=2 && NtotLayer>=8) ;
1608 || ( (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2) && NtotLayer>=8) ; // modified in 2020-08
1609
1610 // gap fraction
1611 double r_gap = 1.0*m_nGap/(m_nGap+NtotLayer);
1612 bool smallGap = r_gap<0.5;
1613 if(debug)
1614 {
1615 if(!enoughVHits) cout<<"not enough V hits, ";
1616 if(!smallGap) cout<<"gap too large 50%, "<<endl;
1617 if(m_unused<3) cout<<"N_unused = "<<m_unused<<"<3"<<endl;
1618 if(m_untried<3) cout<<"N_untried = "<<m_untried<<"<3"<<endl;
1619 }
1620
1621 // --- check conditions
1622 bool good = true;
1623 good=good && m_unused>=2; // 3 -> 2, modified 2020-10
1624 good=good && m_untried>=2;// 3 -> 2, modified 2020-10
1625 good=good && (nCgemLayer==3||NtotLayer>3);
1626 good=good && enoughVHits;
1627 good=good && smallGap;
1628 //good=good && m_maxGap<4;
1629 return good;
1630
1631}

◆ judgeCharge() [1/2]

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

Definition at line 631 of file HoughTrack.cxx.

632{
633 //xHit = hit->getHitPosition().x();
634 //yHit = hit->getHitPosition().y();
635 //double xCenter = m_helix.center().x();
636 //double yCenter = m_helix.center().y();
637 double xCenter = center().x();
638 double yCenter = center().y();
639 double leftOrRight = xHit*yCenter - xCenter*yHit;
640 if(leftOrRight>0) return 1;
641 else return -1;
642}

◆ judgeCharge() [2/2]

int HoughTrack::judgeCharge ( HoughHit * hit)

Definition at line 618 of file HoughTrack.cxx.

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

Referenced by findVHot(), and findXHot().

◆ judgeHalf()

int HoughTrack::judgeHalf ( HoughHit * hit)

Definition at line 605 of file HoughTrack.cxx.

606{
607 double xHit = hit->getHitPosition().x();
608 double yHit = hit->getHitPosition().y();
609 //double xCenter = m_helix.center().x();
610 //double yCenter = m_helix.center().y();
611 double xCenter = center().x();
612 double yCenter = center().y();
613 if(xHit*yCenter > xCenter*yHit)return m_charge;
614 else if(xHit*yCenter < xCenter*yHit)return -m_charge;
615 else return 0;
616}

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

◆ markUsedHot() [1/2]

void HoughTrack::markUsedHot ( int use = 1)

Definition at line 1565 of file HoughTrack.cxx.

1566{
1567 vector<double>::iterator it0_res = m_vecHitResidual.begin();
1568 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1569 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1570 for(; iter!=m_vecHitPnt.end(); iter++)
1571 {
1572 int useOld = (*iter)->getUse();
1573 if(use==-1&&useOld>0) continue;
1574 (*iter)->setUse(use);
1575 if(use==1) {
1576 //cout<<"HoughTrack::markUsedHot(use) add trk ID "<<getTrkID()<<endl;
1577 //(*iter)->addTrkPnt(this);
1578 (*iter)->addTrkID(getTrkID());
1579 vector<double>::iterator it_res = it0_res+(iter-iter0);
1580 //cout<<"add residual "<<*it_res<<endl;
1581 (*iter)->addResid(*it_res);
1582 //vector<double> vecRes = (*iter)->getVecResid();
1583 //cout<<"VecResid.size = "<<vecRes.size()<<endl;
1584 }
1585 }
1586}

◆ markUsedHot() [2/2]

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

Definition at line 1548 of file HoughTrack.cxx.

1549{
1550 vector<HoughHit*>::iterator iter = hitPntList.begin();
1551 for(; iter!=hitPntList.end(); iter++)
1552 {
1553 int useOld = (*iter)->getUse();
1554 if(use==-1&&useOld>0) continue;
1555 (*iter)->setUse(use);
1556
1557 if(use==1) {
1558 //cout<<"HoughTrack::markUsedHot add trk ID "<<getTrkID()<<endl;
1559 //(*iter)->addTrkPnt(this);
1560 (*iter)->addTrkID(getTrkID());
1561 }
1562 }
1563}

◆ nearStereoHits()

bool HoughTrack::nearStereoHits ( )
inline

Definition at line 155 of file HoughTrack.h.

155{return m_nearStereoHits;};

◆ operator=()

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

Definition at line 263 of file HoughTrack.cxx.

264{
265 if (this == & other) return * this;
266
267 Helix::operator=(other);
268 m_trkID = other.m_trkID;
269 m_charge = other.m_charge;
270 m_angle = other.m_angle;
271 m_rho = other.m_rho;
272 m_flag = other.m_flag;
273
274 m_dAngle = other.m_dAngle;
275 m_dRho = other.m_dRho;
276 m_dTanl = other.m_dTanl;
277 m_dDz = other.m_dDz;
278
279 m_dr = other.m_dr;
280 m_phi0 = other.m_phi0;
281 m_kappa = other.m_kappa;
282 m_dz = other.m_dz;
283 m_tanl = other.m_tanl;
284
285 m_chi2 = other.m_chi2;
286 m_trkRecoTrk = other.m_trkRecoTrk;
287
288 m_circleFitStat=other.m_circleFitStat;
289 m_vecHitPnt = other.m_vecHitPnt;
290 m_vecHitResidual = other.m_vecHitResidual;
291 m_vecHitChi2 = other.m_vecHitChi2;
292
293 m_vecStereoHitPnt = other.m_vecStereoHitPnt;
294 m_vecStereoHitRes = other.m_vecStereoHitRes;
295 m_mcTrack = other.m_mcTrack;
296 m_cgemHitVector = other.m_cgemHitVector;
297 m_vecRecMdcHit = other.m_vecRecMdcHit;
298 m_nearStereoHits = other.m_nearStereoHits;
299 return * this;
300}
Helix::operator=
Helix & operator=(const Helix &)
Copy operator.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:403

◆ print()

void HoughTrack::print ( )

Definition at line 1026 of file HoughTrack.cxx.

1027{
1028 cout
1029 <<setw(12)<<"trkID:" <<setw(15)<<m_trkID
1030 <<setw(12)<<"flag:" <<setw(15)<<m_flag
1031 //<<setw(12)<<"q:" <<setw(15)<<m_charge
1032 <<setw(12)<<"pt:" <<setw(15)<<pt()
1033 <<setw(12)<<"angle:" <<setw(15)<<m_angle
1034 <<setw(12)<<"rho:" <<setw(15)<<m_rho
1035 <<endl
1036 <<setw(12)<<"dr:" <<setw(15)<<dr()
1037 <<setw(12)<<"phi0:" <<setw(15)<<phi0()
1038 <<setw(12)<<"kappa:" <<setw(15)<<kappa()
1039 <<setw(12)<<"dz:" <<setw(15)<<dz()
1040 <<setw(12)<<"tanl:" <<setw(15)<<tanl()
1041 <<setw(12)<<"chi2:" <<setw(15)<<getChi2()
1042 <<endl
1043 //<<setw(12)<<"dr:" <<setw(15)<<m_dr
1044 //<<setw(12)<<"phi0:" <<setw(15)<<m_phi0
1045 //<<setw(12)<<"kappa:" <<setw(15)<<m_kappa
1046 //<<setw(12)<<"dz:" <<setw(15)<<m_dz
1047 //<<setw(12)<<"tanl:" <<setw(15)<<m_tanl
1048 //<<endl
1049 //<<setw(12)<<"dAngle:" <<setw(15)<<m_dAngle
1050 //<<setw(12)<<"dRho:" <<setw(15)<<m_dRho
1051 //<<setw(12)<<"dTanl:" <<setw(15)<<m_dTanl
1052 //<<setw(12)<<"dDz:" <<setw(15)<<m_dDz
1053 //<<setw(12)<<"chi2:" <<setw(15)<<m_chi2
1054 //<<setw(12)<<":"<<setw(15)<<
1055 //<<setw(12)<<":"<<setw(15)<<
1056 //<<setw(12)<<":"<<setw(15)<<
1057 ;
1058 int nHot = getHotList().size();
1059 int nAxialHot = 0;
1060 int nStereoHot = 0;
1061 int nXCluster = 0;
1062 int nXVCluster = 0;
1063 vector<HoughHit*> hotList;
1064 vector<HoughHit*> axialHot = getVecHitPnt();
1065 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
1066 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
1067 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXCluster++;
1068 if((**hotIt).getHitType()==HoughHit::MDCHIT)nAxialHot++;
1069 }
1070 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
1071 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXVCluster++;
1072 if((**hotIt).getHitType()==HoughHit::MDCHIT)nStereoHot++;
1073 }
1074 cout
1075 <<setw(12)<<"nHot:" <<setw(15)<<nHot
1076 <<setw(12)<<"nAxialHot0:" <<setw(15)<<nAxialHot
1077 <<setw(12)<<"nStereoHot0:" <<setw(15)<<nStereoHot
1078 <<setw(12)<<"nXCluster0:" <<setw(15)<<nXCluster
1079 <<setw(12)<<"nXVCluster0:" <<setw(15)<<nXVCluster
1080 <<endl;
1081 //cout<<endl;
1082}
HoughTrack::getHotList
vector< HoughHit * > getHotList(int type=2)
Definition HoughTrack.cxx:2214

Referenced by HoughFinder::getMcParticleCol().

◆ printHot()

void HoughTrack::printHot ( int type = 2)

Definition at line 1084 of file HoughTrack.cxx.

1085{
1086 vector<HoughHit*> hitList = getHotList(type);
1087 //vector<HoughHit*> hitList = getVecStereoHitPnt();
1088 for(vector<HoughHit*>::iterator hotIt = hitList.begin(); hotIt != hitList.end();hotIt++){
1089 HoughHit* hitIt = *hotIt;
1090 hitIt->print();
1091 //if(hitIt->getLayer()>3)continue;
1092 //cout<<hitIt->getHitID();
1093 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
1094 //if(hitIt->getHitType()==HoughHit::CGEMMCHIT||hitIt->getHitType()==HoughHit::MDCMCHIT)cout<<" halfCircle:"<<hitIt->getHalfCircle();
1095 //else if(hitIt->getPairHit()!=NULL)cout<<" halfCircle:"<<hitIt->getPairHit()->getHalfCircle();
1096 //HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1097 //double res = hitIt->residual(this, positionOntrack, positionOnDetector);
1098 //cout<<" res:"<<setw(10)<<res<<" ";
1099 //cout<<endl;
1100 }
1101 cout<<endl;
1102 /*
1103 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1104 cout<<setw(4)<<(*hitIt)->getHitID();
1105 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1106 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1107 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1108 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1109 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1110 vector<int> trkID = (*hitIt)->getTrkID();
1111 cout<<"TrkID:"<<setw(4)<<trkID[0];
1112 }
1113 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1114 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1115 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1116 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1117 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1118 cout<<"res:"<<setw(10)<<res<<" ";
1119 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1120 }
1121 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1122 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1123 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1124 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1125 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1126 cout<<"res:"<<setw(10)<<res<<" ";
1127 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1128 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1129 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1130 }
1131 if((*hitIt)->getPairHit()!=NULL){
1132 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1133 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1134 }
1135 //HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1136 //cout<<flightLength(hitPoint);
1137 //else cout<<"NULL";
1138 cout<<endl;
1139 }
1140
1141 //vector<HoughHit*> hitList = getVecHitPnt();
1142 hitList = getVecStereoHitPnt();
1143 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1144 cout<<setw(4)<<(*hitIt)->getHitID();
1145 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1146 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1147 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1148 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1149 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1150 vector<int> trkID = (*hitIt)->getTrkID();
1151 cout<<"TrkID:"<<setw(4)<<trkID[0];
1152 }
1153 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1154 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1155 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1156 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1157 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1158 cout<<"res:"<<setw(10)<<res<<" ";
1159 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1160 }
1161 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1162 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1163 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1164 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1165 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1166 cout<<"res:"<<setw(10)<<res<<" ";
1167 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1168 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1169 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1170 }
1171 if((*hitIt)->getPairHit()!=NULL){
1172 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1173 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1174 }
1175//HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1176//cout<<flightLength(hitPoint);
1177//else cout<<"NULL";
1178cout<<endl;
1179}
1180// */
1181}
HoughHit::print
void print()
Definition HoughHit.cxx:327

◆ printRecMdcHitVec()

void HoughTrack::printRecMdcHitVec ( )

Definition at line 1200 of file HoughTrack.cxx.

1201{
1202 cout<<"HoughTrack::printRecMdcHitVec(): "<<endl;
1203 vector<RecMdcHit>::iterator iter_recMdcHit = m_vecRecMdcHit.begin();
1204 for(; iter_recMdcHit!=m_vecRecMdcHit.end(); iter_recMdcHit++)
1205 {
1206 Identifier mdcid = iter_recMdcHit->getMdcId();
1207 int layer = MdcID::layer(mdcid);
1208 int wire = MdcID::wire(mdcid);
1209 cout<<" hit at layer "<<layer<<" wire "<<wire<<endl;
1210 }
1211}
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:54

◆ releaseSelHits()

void HoughTrack::releaseSelHits ( )

Definition at line 1933 of file HoughTrack.cxx.

1934{
1935 vector<HoughHit*>::iterator it = m_vecHitPnt.begin();
1936 for(; it!=m_vecHitPnt.end(); it++)
1937 {
1938 (*it)->dropTrkID(m_trkID);// erase trk id and residual from the Hough hit
1939 }
1940}

◆ 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 143 of file HoughTrack.h.

143{m_mcTrack = mcTrack;}

◆ setRecMdcHitVec()

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

Definition at line 151 of file HoughTrack.h.

151{ 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 982 of file HoughTrack.cxx.

983{
984 std::sort(hotList.begin(),hotList.end(),sortByLayer);
985}
sortByLayer
bool sortByLayer(HoughHit *hit1, HoughHit *hit2)
Definition HoughTrack.cxx:977

Referenced by getHotList().

◆ update()

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

Definition at line 999 of file HoughTrack.cxx.

1000{
1001 m_angle = angle;
1002 m_rho = rho;
1003 m_phi0 = (rho>0)? angle:angle+M_PI;
1004 m_phi0 = (m_charge<0)? m_phi0:m_phi0+M_PI;
1005 if(m_phi0>2*M_PI)m_phi0-=2*M_PI;
1006 if(m_phi0<0) m_phi0+=2*M_PI;
1007 m_kappa = fabs(m_alpha*rho)*m_charge;
1008 HepPoint3D pivot(0,0,0);
1009 HepVector a(5,0);
1010 HepSymMatrix Ea(5,0);
1011 a(2) = m_phi0;
1012 a(3) = m_kappa;
1013 set(pivot,a,Ea);
1014}

◆ updateCirclePar()

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

Definition at line 1017 of file HoughTrack.cxx.

1018{
1019 m_dr = dr;
1020 m_phi0=phi0;
1021 m_kappa=kappa;
1022 updateHelix();
1023
1024}

◆ updateHelix()

void HoughTrack::updateHelix ( )

Definition at line 987 of file HoughTrack.cxx.

988{
989 pivot(HepPoint3D(0,0,0));
990 HepVector hepVector(5,0);
991 hepVector(1) = m_dr;
992 hepVector(2) = m_phi0;
993 hepVector(3) = m_kappa;
994 hepVector(4) = m_dz;
995 hepVector(5) = m_tanl;
996 a(hepVector);
997}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37

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

◆ updateLayerStat()

void HoughTrack::updateLayerStat ( )

Definition at line 1183 of file HoughTrack.cxx.

1184{
1185 m_setLayer.clear();
1186 m_nearStereoHits=false;
1187 int nCgemLayer(0), N_ODC1_nearStereo(0), N_ODC2_nearStereo(0);
1188 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
1189 {
1190 int iLayer=(*hotIt)->getLayer();
1191 m_setLayer.insert(iLayer);
1192 if(iLayer<3) nCgemLayer++;
1193 else if(iLayer>=16&&iLayer<=19) N_ODC1_nearStereo++;
1194 else if(iLayer>=36&&iLayer<=39) N_ODC2_nearStereo++;
1195 }
1196 if(N_ODC1_nearStereo>=2||N_ODC2_nearStereo>=2) m_nearStereoHits=true;
1197}

Member Data Documentation

◆ m_clearTrack

int HoughTrack::m_clearTrack =1
static

Definition at line 146 of file HoughTrack.h.

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

◆ m_cut

TGraph * HoughTrack::m_cut = {NULL}
static

Definition at line 18 of file HoughTrack.h.

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

◆ m_useCgemInGlobalFit

int HoughTrack::m_useCgemInGlobalFit =3
static

Definition at line 147 of file HoughTrack.h.

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

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