#include <MucRec2DRoad.h>
Inheritance diagram for MucRec2DRoad:
Public Member Functions | |
| void | AttachHit (MucRecHit *hit) |
| Attach the given hit to this road. | |
| void | AttachHit (MucRecHit *hit) |
| Attach the given hit to this road. | |
| void | AttachHitNoFit (MucRecHit *hit) |
| Attach the given hit to this road, but not fit this road. | |
| void | AttachHitNoFit (MucRecHit *hit) |
| Attach the given hit to this road, but not fit this road. | |
| int | Fit (const float &x, const float &y, const float &z, float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Fit (with weights) the hit center to a straight line. | |
| int | Fit (const float &x, const float &y, const float &z, float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Fit (with weights) the hit center to a straight line. | |
| int | GetDegreesOfFreedom () const |
| How many degrees of freedom in the trajectory fit? | |
| int | GetDegreesOfFreedom () const |
| How many degrees of freedom in the trajectory fit? | |
| MucRecHit * | GetHit (const int &gap) const |
| Get a pointer to the first found hit attached in a particular gap. | |
| MucRecHit * | GetHit (const int &gap) const |
| Get a pointer to the first found hit attached in a particular gap. | |
| float | GetHitDistance (const MucRecHit *hit) |
| Distance to a hit. | |
| float | GetHitDistance (const MucRecHit *hit) |
| Distance to a hit. | |
| vector< MucRecHit * > | GetHits () const |
| vector< MucRecHit * > | GetHits () const |
| vector< Identifier > | GetHitsID () const |
| Get indices of all hits in the road. | |
| vector< Identifier > | GetHitsID () const |
| Get indices of all hits in the road. | |
| int | GetHitsPerGap (const int &gap) const |
| How many hits per gap does this road contain? | |
| int | GetHitsPerGap (const int &gap) const |
| How many hits per gap does this road contain? | |
| int | GetIndex () const |
| A unique identifier for this road in the current event. | |
| int | GetIndex () const |
| A unique identifier for this road in the current event. | |
| float | GetIntercept () const |
| Intercept of trajectory. | |
| float | GetIntercept () const |
| Intercept of trajectory. | |
| int | GetLastGap () const |
| Which gap is the last one with hits attached to this road? | |
| int | GetLastGap () const |
| Which gap is the last one with hits attached to this road? | |
| int | GetMaxHitsPerGap () const |
| How many hits were attached in the gap with the most attached hits? | |
| int | GetMaxHitsPerGap () const |
| How many hits were attached in the gap with the most attached hits? | |
| int | GetNGapsWithHits () const |
| How many gaps provide hits attached to this road? | |
| int | GetNGapsWithHits () const |
| How many gaps provide hits attached to this road? | |
| int | GetNSharedHits (const MucRec2DRoad *road) const |
| How many hits do two roads share? | |
| int | GetNSharedHits (const MucRec2DRoad *road) const |
| How many hits do two roads share? | |
| int | GetOrient () const |
| In which orientation was this road found? | |
| int | GetOrient () const |
| In which orientation was this road found? | |
| int | GetPart () const |
| In which part was this road found? | |
| int | GetPart () const |
| In which part was this road found? | |
| void | GetPosSigma (float &possigma) const |
| void | GetPosSigma (float &possigma) const |
| bool | GetQuadFitOk () |
| bool | GetQuadFitOk () |
| float | GetReducedChiSquare () const |
| Chi-square parameter (per degree of freedom) of the trajectory fit. | |
| float | GetReducedChiSquare () const |
| Chi-square parameter (per degree of freedom) of the trajectory fit. | |
| float | GetSearchWindowSize (const int &gap) const |
| Determine the size of the search window in the given gap. | |
| float | GetSearchWindowSize (const int &gap) const |
| Determine the size of the search window in the given gap. | |
| int | GetSeg () const |
| In which segment was this road found? | |
| int | GetSeg () const |
| In which segment was this road found? | |
| void | GetSimpleFitParams (float &a, float &b, float &c, int &whichhalf, float &sigmaa, float &sigmab, float &sigmac, float &chisq, int &ndof) const |
| Get the parameters from the simple quad fit. | |
| void | GetSimpleFitParams (float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) const |
| Get the parameters from the simple fit. | |
| void | GetSimpleFitParams (float &a, float &b, float &c, int &whichhalf, float &sigmaa, float &sigmab, float &sigmac, float &chisq, int &ndof) const |
| Get the parameters from the simple quad fit. | |
| void | GetSimpleFitParams (float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) const |
| Get the parameters from the simple fit. | |
| float | GetSlope () const |
| Slope of trajectory. | |
| float | GetSlope () const |
| Slope of trajectory. | |
| int | GetTotalHits () const |
| How many hits in all does this road contain? | |
| int | GetTotalHits () const |
| How many hits in all does this road contain? | |
| void | GetVertexPos (float &x, float &y, float &z) const |
| Position of the vertex point. | |
| void | GetVertexPos (float &x, float &y, float &z) const |
| Position of the vertex point. | |
| bool | HasHit (MucRecHit *hit) const |
| Does the hit exist in the road . | |
| bool | HasHit (MucRecHit *hit) const |
| Does the hit exist in the road . | |
| bool | HasHitInGap (const int &gap) const |
| Does this road contain any hits in the given gap? | |
| bool | HasHitInGap (const int &gap) const |
| Does this road contain any hits in the given gap? | |
| MucRec2DRoad (const MucRec2DRoad &source) | |
| Copy constructor. | |
| MucRec2DRoad (const int &part=0, const int &seg=0, const int &orient=0, const float &xVertex=0.0, const float &yVertex=0.0, const float &zVertex=0.0, const int &fittingMethod=0) | |
| Constructor. | |
| MucRec2DRoad (const MucRec2DRoad &source) | |
| Copy constructor. | |
| MucRec2DRoad (const int &part=0, const int &seg=0, const int &orient=0, const float &xVertex=0.0, const float &yVertex=0.0, const float &zVertex=0.0, const int &fittingMethod=0) | |
| Constructor. | |
| MucRec2DRoad & | operator= (const MucRec2DRoad &orig) |
| Assignment constructor. | |
| MucRec2DRoad & | operator= (const MucRec2DRoad &orig) |
| Assignment constructor. | |
| void | PrintHitsInfo () const |
| Print Hits Infomation. | |
| void | PrintHitsInfo () const |
| Print Hits Infomation. | |
| void | Project (const int &gap, float &x, float &y, float &z, float &x2, float &y2, float &z2) |
| Where does the trajectory of this road intersect a specific gap? | |
| void | Project (const int &gap, float &x, float &y, float &z, float &x2, float &y2, float &z2) |
| Where does the trajectory of this road intersect a specific gap? | |
| void | SetIndex (const int &index) |
| Set the index for this road. | |
| void | SetIndex (const int &index) |
| Set the index for this road. | |
| void | SetMaxNSkippedGaps (const int &nGaps) |
| void | SetMaxNSkippedGaps (const int &nGaps) |
| int | SimpleFit (float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Calculate the best-fit straight line with "simple" weights. | |
| int | SimpleFit (float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Calculate the best-fit straight line with "simple" weights. | |
| int | SimpleFitNoCurrentGap (int currentgap, float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Calculate the best-fit straight line with "simple" weights. not use current gap!!! | |
| int | SimpleFitNoCurrentGap (int currentgap, float &slope, float &intercept, float &sigmaSlope, float &sigmaIntercept, float &chisq, int &ndof) |
| Calculate the best-fit straight line with "simple" weights. not use current gap!!! | |
| float | WeightFunc (const float &chisq, const float &distance) const |
| float | WeightFunc (const float &chisq, const float &distance) const |
| float | WindowFunc (const float &chisq, const float &distance) const |
| float | WindowFunc (const float &chisq, const float &distance) const |
| ~MucRec2DRoad () | |
| Destructor. | |
| ~MucRec2DRoad () | |
| Destructor. | |
Private Member Functions | |
| void | CountHits () |
| Calculate the hit counts and last gap quantities. | |
| void | CountHits () |
| Calculate the hit counts and last gap quantities. | |
Private Attributes | |
| float | m_Chi2 |
| int | m_DOF |
| bool | m_FitOK |
| int | m_fittingMethod |
| vector< float > | m_HitDistance |
| vector< float > | m_HitDistance |
| int | m_Index |
| int | m_LastGap |
| int | m_MaxHitsPerGap |
| int | m_MaxNSkippedGaps |
| int | m_Orient |
| int | m_Part |
| vector< MucRecHit * > | m_pHits |
| vector< MucRecHit * > | m_pHits |
| bool | m_QuadFitOK |
| int | m_Seg |
| float | m_SimpleChi2 |
| int | m_SimpleDOF |
| bool | m_SimpleFitOK |
| float | m_SimpleIntercept |
| float | m_SimpleInterceptSigma |
| float | m_SimplePosSigma |
| float | m_SimpleQuad_a |
| float | m_SimpleQuad_aSigma |
| float | m_SimpleQuad_b |
| float | m_SimpleQuad_bSigma |
| float | m_SimpleQuad_c |
| float | m_SimpleQuad_cSigma |
| int | m_SimpleQuad_whichhalf |
| float | m_SimpleSlope |
| float | m_SimpleSlopeSigma |
| int | m_TotalHits |
| Hep3Vector | m_VertexPos |
| Hep3Vector | m_VertexSigma |
|
||||||||||||||||||||||||||||||||
|
Constructor.
00035 : m_VertexPos(xVertex, yVertex, zVertex), 00036 m_VertexSigma(0.0, 0.0, 0.0), 00037 m_Part(part), m_Seg(seg), m_Orient(orient), 00038 m_Chi2(0.0), m_DOF(0), 00039 m_MaxHitsPerGap(0), m_LastGap(0), m_TotalHits(0), 00040 m_FitOK(false), m_QuadFitOK(false), 00041 m_HitDistance(5, float(muckInfinity)), 00042 m_pHits(0), m_fittingMethod(fittingMethod) 00043 { }
|
|
|
Copy constructor.
00073 : m_VertexPos(source.m_VertexPos), 00074 m_VertexSigma(source.m_VertexSigma), 00075 m_Index(source.m_Index), 00076 m_Part(source.m_Part), m_Seg(source.m_Seg), m_Orient(source.m_Orient), 00077 m_Chi2(source.m_Chi2), m_DOF(source.m_DOF), 00078 m_MaxHitsPerGap(source.m_MaxHitsPerGap), 00079 m_LastGap(source.m_LastGap), m_TotalHits(source.m_TotalHits), 00080 m_FitOK(source.m_FitOK), 00081 m_HitDistance(source.m_HitDistance), 00082 m_pHits(source.m_pHits), 00083 m_fittingMethod(source.m_fittingMethod) 00084 { }
|
|
|
Destructor.
00088 { }
|
|
||||||||||||||||||||||||||||||||
|
Constructor.
|
|
|
Copy constructor.
|
|
|
Destructor.
|
|
|
Attach the given hit to this road.
|
|
|
Attach the given hit to this road.
00101 {
00102 // cout << "MucRec2DRoad::AttachHit-I0 hit = " << hit << endl;
00103 if (!hit) {
00104 cout << "MucRec2DRoad::AttachHit-E1 null hit pointer!" << endl;
00105 return ;
00106 }
00107
00108 int gap = hit->Gap();
00109 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00110 // The gap number of the hit is out of range.
00111 cout << "MucRec2DRoad::AttachHit-W3 bad gap number = " << gap
00112 << endl;
00113 return;
00114 }
00115
00116 // Attach the hit to the road.
00117 //bool hitExist = false;
00118 for (int i = 0; i < (int)m_pHits.size(); i++) {
00119 if (m_pHits[i] == hit) return;
00120 }
00121 m_pHits.push_back(hit);
00122 //cout << "before Count " << m_pHits.size() << endl;
00123 // m_HitDistance[gap] = dX;
00124
00125 // Now recalculate the total number of hits and the max. number of
00126 // hits per gap.
00127 CountHits();
00128 //cout << "after Count " << m_pHits.size() << endl;
00129
00130 // Redo the "simple" least-squares fit to the positions ...
00131 float a, sa, b, sb, chi;
00132 int n;
00133 int status = SimpleFit(a, b, sa, sb, chi, n);
00134 if (status < 0) {
00135 //cout << "MucRec2DRoad::AttachHit-E5 simple fit fail status = " << status << endl;
00136 }
00137
00138 //cout << "Hit center = " << hit->GetCenterPos() << endl;
00139 //cout << "After attach hit, slope = " << a << " intercept = " << b << endl;
00140 }
|
|
|
Attach the given hit to this road, but not fit this road.
|
|
|
Attach the given hit to this road, but not fit this road.
00146 {
00147 // cout << "MucRec2DRoad::AttachHit-I0 hit = " << hit << endl;
00148 if (!hit) {
00149 cout << "MucRec2DRoad::AttachHit-E1 null hit pointer!" << endl;
00150 return ;
00151 }
00152
00153 int gap = hit->Gap();
00154 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00155 // The gap number of the hit is out of range.
00156 cout << "MucRec2DRoad::AttachHit-W3 bad gap number = " << gap
00157 << endl;
00158 return;
00159 }
00160
00161 // Attach the hit to the road.
00162 //bool hitExist = false;
00163 for (int i = 0; i < (int)m_pHits.size(); i++) {
00164 if (m_pHits[i] == hit) return;
00165 }
00166 m_pHits.push_back(hit);
00167 //cout << "before Count " << m_pHits.size() << endl;
00168 // m_HitDistance[gap] = dX;
00169
00170 // Now recalculate the total number of hits and the max. number of
00171 // hits per gap.
00172 CountHits();
00173 //cout << "after Count " << m_pHits.size() << endl;
00174
00175 }
|
|
|
Calculate the hit counts and last gap quantities.
|
|
|
Calculate the hit counts and last gap quantities.
01336 {
01337 m_MaxHitsPerGap = 0;
01338 m_TotalHits = 0;
01339 m_LastGap = 0;
01340
01341 vector<MucRecHit*>::const_iterator iHit;
01342 const int ngap = (int)MucID::getGapNum(m_Part);
01343 vector<int> HitsPerGap;
01344 for (int gap = 0; gap < ngap; gap++) {
01345 HitsPerGap.push_back(0);
01346 }
01347
01348 // Fill HitsPerGap
01349 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01350 if ( !(*iHit) ) {
01351 cout << "MucRec2DRoad::CountHits()-E1 null hit pointer !" << endl;
01352 return;
01353 }
01354 else {
01355 int gap = (*iHit)->Gap();
01356 HitsPerGap[gap]++;
01357 //cout << "gap " << gap << endl;
01358 }
01359 }
01360
01361 // Find the gap from where we should stop count,
01362 // namely in front of the gap there are more than
01363 // m_MaxNSkippedGaps gaps containing no hit.
01364
01365 int StopGap = ngap;
01366 /*
01367 for (int gap = m_MaxNSkippedGaps; gap < ngap; gap++) {
01368 int SubTotalHits = 0;
01369 for (int igap = gap-m_MaxNSkippedGaps; igap <= gap; igap++) {
01370 SubTotalHits += HitsPerGap[igap];
01371 }
01372 if (SubTotalHits == 0 ) {
01373 StopGap = gap - m_MaxNSkippedGaps;
01374 break;
01375 }
01376 }
01377 */
01378
01379 //cout << "StopGap " << StopGap << endl;
01380 // Now we might get correct counting on hits, last gap and MaxHitsPerGap
01381 for (int gap = 0; gap < StopGap; gap++) {
01382 m_TotalHits += HitsPerGap[gap];
01383 if(m_MaxHitsPerGap < HitsPerGap[gap]) m_MaxHitsPerGap = HitsPerGap[gap];
01384 if(HitsPerGap[gap] > 0) m_LastGap = gap;
01385 }
01386 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Fit (with weights) the hit center to a straight line.
|
|
||||||||||||||||||||||||||||||||||||||||
|
Fit (with weights) the hit center to a straight line.
00576 {
00577 int status;
00578
00579 // If the "simple" fit has not been done yet, do it now.
00580 if (!m_SimpleFitOK) {
00581 // Least-squares fit to the positions ...
00582 float a, sa, b, sb, chi;
00583 int n;
00584 status = SimpleFit(a, b, sa, sb, chi, n);
00585 if (status < 0) {
00586 //cout << "MucRec2DRoad::Fit-E2 simple fit fail status = "
00587 // << status << endl;
00588 return status;
00589 }
00590 }
00591
00592 // Assign to temporary arrays to be used in fit.
00593 float px[100];
00594 float py[100];
00595 float pw[100];
00596 int npoints = 0;
00597 float dx, dy, dr;
00598
00599
00600 // include vertex point when do the fancy fitting
00601 //if (m_Orient == kHORIZ) {
00602 // px[npoints] = m_VertexPos.y();
00603 // dx = px[npoints] - y;
00604 //} else {
00605 // px[npoints] = m_VertexPos.x();
00606 // dx = px[npoints] - x;
00607 //}
00608 //pz[npoints] = m_VertexPos.z();
00609 //dz = pz[npoints] - z;
00610 //dr = sqrt(dx*dx + dz*dz);
00611 //pw[npoints] = WeightFunc(m_SimpleChi2,dr);
00612 //npoints++;
00613
00614 vector<MucRecHit*>::const_iterator iHit;
00615
00616 // Determine the weights based on the chi-square of the fit
00617 // (the scatter of the points should be roughly related to the
00618 // momentum) and the distance from the center to the estimated
00619 // location.
00620
00621 //cout << m_pHits.size() << endl;
00622 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00623 if (*iHit) { // Check for a null pointer.
00624
00625 Hep3Vector center = (*iHit)->GetCenterPos();
00626 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00627
00628 if (m_Part == 1) { // change from 0 to 1 at 2006.11.30
00629 if ( m_Orient == 0) {
00630 px[npoints] = center.z();
00631 dx = px[npoints] - z;
00632 py[npoints] = sqrt(center.x()*center.x()
00633 + center.y()*center.y());
00634 if(m_Seg==2) py[npoints] = center.y(); //deal with seg2 seperately! because there is a hole here. 2006.11.9
00635 dy = py[npoints] - sqrt(x*x + y*y);
00636 }
00637 else {
00638 px[npoints] = center.x();
00639 dx = px[npoints] - x;
00640 py[npoints] = center.y();
00641 dy = py[npoints] - y;
00642 }
00643 }
00644 else {
00645 if ( m_Orient == 0) {
00646 px[npoints] = center.z();
00647 dx = px[npoints] - z;
00648 py[npoints] = center.y();
00649 dy = py[npoints] - y;
00650 }
00651 else {
00652 px[npoints] = center.z();
00653 dx = px[npoints] - z;
00654 py[npoints] = center.x();
00655 dy = py[npoints] - x;
00656 }
00657 }
00658
00659 dr = sqrt(dx*dx + dy*dy);
00660 pw[npoints] = WeightFunc(m_SimpleChi2, dr);
00661 //pw[npoints] = WeightFunc(m_SimpleChi2,dr);
00662
00663 // cout << " " << npoints << " "
00664 // << px[npoints] << " " << py[npoints] << " " << pw[npoints]
00665 // << " " << dx << " " << dy << " " << dr
00666 // << endl;
00667
00668 npoints++;
00669
00670 if(npoints > 99)
00671 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00672 return -1;
00673 }
00674
00675 }
00676 }
00677
00678 // Refit ...
00679 ndof = npoints - 2;
00680 if (npoints == 1) {
00681 if (m_Part == 1) { // change from 0 to 1 at 2006.11.30
00682 if ( m_Orient == 0) {
00683 px[npoints] = m_VertexPos.z();
00684 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00685 + m_VertexPos.y()*m_VertexPos.y());
00686 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00687 }
00688 else {
00689 px[npoints] = m_VertexPos.x();
00690 py[npoints] = m_VertexPos.y();
00691 }
00692 }
00693 else {
00694 if ( m_Orient == 0) {
00695 px[npoints] = m_VertexPos.z();
00696 py[npoints] = m_VertexPos.y();
00697 }
00698 else {
00699 px[npoints] = m_VertexPos.z();
00700 py[npoints] = m_VertexPos.x();
00701 }
00702 }
00703 pw[npoints] = 1.0;
00704 npoints++;
00705 }
00706 else {
00707 if (npoints == 0) {
00708 return -1;
00709 }
00710 }
00711
00712 MucRecLineFit fit;
00713 //cout << "npoints " << npoints << endl;
00714 status = fit.LineFit(px, py, pw, npoints,
00715 &slope, &intercept, &chisq,
00716 &sigmaSlope, &sigmaIntercept);
00717 m_DOF = ndof;
00718 m_Chi2 = chisq;
00719
00720 if (status < 0) {
00721 //cout << "MucRec2DRoad::Fit-E3 fit fail status = " << status << endl;
00722 }
00723
00724 return status;
00725 }
|
|
|
How many degrees of freedom in the trajectory fit?
|
|
|
How many degrees of freedom in the trajectory fit?
01059 { return m_SimpleDOF; }
|
|
|
Get a pointer to the first found hit attached in a particular gap.
|
|
|
Get a pointer to the first found hit attached in a particular gap.
01127 {
01128 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01129 cout << "MucRec2DRoad::Hit-E1 invalid gap = " << gap << endl;
01130 return 0;
01131 }
01132
01133 vector<MucRecHit*>::const_iterator iHit;
01134
01135 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01136 if (*iHit) { // Check for a null pointer.
01137 if ( (*iHit)->Gap() == gap) {
01138 return (*iHit);
01139 }
01140 }
01141 }
01142
01143 return 0L;
01144 }
|
|
|
Distance to a hit.
|
|
|
Distance to a hit.
01150 {
01151 // Use straight-line projection?
01152 if (!hit) {
01153 cout << "MucRec2DRoad::GetHitDistance-E1 null hit pointer!" << endl;
01154 return muckInfinity;
01155 }
01156
01157 int gap = hit->Gap();
01158 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01159 cout << "MucRec2DRoad::HitDistance-W2 bad gap number = " << gap << endl;
01160 return muckInfinity;
01161 }
01162
01163 if ( hit->GetGap()->Orient() != m_Orient) {
01164 // The orientation of the hit is different from that of the road.
01165 cout << "MucRec2DRoad::GetHitDistance-W3 "
01166 << " hit has wrong orientation = "
01167 << hit->GetGap()->Orient()
01168 << endl;
01169 return muckInfinity;
01170 }
01171
01172 HepPoint3D r = hit->GetCenterPos();
01173 HepPoint3D rl = hit->GetGap()->TransformToGap(r);
01174 // cout << "hit center " << r << endl;
01175 // cout << "hit center local " << rl << endl;
01176
01177 // Determine the projection point of the "simple" fit to the desired gap.
01178 // FIXME(?): need to use fit with fancy weights instead?
01179 float delta, delta1,delta2;
01180 float x0, y0, z0;
01181 float x2, y2, z2;
01182 //float sigmaX0, sigmaY0, sigmaZ0;
01183
01184 //cout << "y:x = " << m_SimpleSlope << endl;
01185
01186 Project(gap, x0, y0, z0, x2, y2, z2);
01187 // cout << " gap = " << gap << " x0 = " << x0
01188 // << " y0 = " << y0 << " z0 = " << z0
01189 // << endl;
01190
01191 if(x0==y0&&x0==z0&&x0==0) {x0+=-9999;y0+=-9999;z0+=-9999;}//cout<<"wrong intersection"<<endl;}
01192 if(x2==y2&&x2==z2&&x2==0) {x2+=-9999;y2+=-9999;z2+=-9999;}//cout<<"wrong intersection"<<endl;} //wrong intersection!!!
01193 HepPoint3D s = HepPoint3D(x0, y0, z0);
01194 HepPoint3D s_2 = HepPoint3D(x2, y2, z2);
01195 HepPoint3D sl = hit->GetGap()->TransformToGap(s);
01196 HepPoint3D s_2l = hit->GetGap()->TransformToGap(s_2);
01197 //cout << "project center " << s << endl;
01198
01199 // cout << "project center local sl= " << sl << endl;
01200 // cout << "project center local sl= " << s_2l << endl;
01201 // cout << "project center local rl= " << rl << endl;
01202 if ( m_Part == 0 ) {
01203 if ( m_Orient == 0 ) {
01204 delta1 = fabs( sl.y() - rl.y() );
01205 delta2= fabs( s_2l.y() - rl.y() );
01206 }
01207 else {
01208 delta1 = fabs( sl.x() - rl.x() );
01209 delta2= fabs( s_2l.x() - rl.x() );
01210 }
01211 }
01212 else {
01213 if ( m_Orient == 0 ) {
01214 delta1 = fabs( sl.y() - rl.y() );
01215 delta2= fabs( s_2l.y() - rl.y() );
01216 }
01217 else {
01218 delta1 = fabs( sl.x() - rl.x() );
01219 delta2= fabs( s_2l.x() - rl.x() );
01220 }
01221 }
01222 // if(which==0) {
01223 // if(delta1<delta2)delta = delta1;
01224 // else delta = delta2;
01225 // }
01226 //cout<<"which = "<<which<<" distance = "<<delta1<<" "<<delta2<<endl;
01227
01228 if(delta1 < delta2) return delta1;
01229 else return delta2;
01230 }
|
|
|
|
|
|
01435 {
01436 return m_pHits;
01437
01438 }
|
|
|
Get indices of all hits in the road.
|
|
|
Get indices of all hits in the road.
01414 {
01415 vector<Identifier> idCon;
01416
01417 vector<MucRecHit*>::const_iterator iHit;
01418 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01419 if (*iHit) { // Check for a null pointer.
01420 Identifier id = (*iHit)->GetID();
01421 idCon.push_back(id);
01422 /*
01423 cout << " MucRec2DRoad::HitIndices gap orientation twopack= "
01424 << (*iHit)->ChannelID().Plane() << " "
01425 << (*iHit)->ChannelID().Orient() << " "
01426 << (*iHit)->ChannelID().TwoPack() << endl;
01427 */
01428 }
01429 }
01430 return idCon;
01431 }
|
|
|
How many hits per gap does this road contain?
|
|
|
How many hits per gap does this road contain?
00973 {
00974 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00975 cout << "MucRec2DRoad::GetHitsPerGap-E1 invalid gap = " << gap << endl;
00976 return 0;
00977 }
00978
00979 vector<MucRecHit*>::const_iterator iHit;
00980 int hitsInGap = 0;
00981
00982 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00983
00984 if ( !(*iHit) ) {
00985 cout << "MucRec2DRoad::GetHitsPerGap()-E2 null hit pointer !" << endl;
00986 return 0;
00987 }
00988 else {
00989 if( gap == (*iHit)->Gap() ) hitsInGap += 1;
00990 }
00991 }
00992
00993 return hitsInGap;
00994 }
|
|
|
A unique identifier for this road in the current event.
|
|
|
A unique identifier for this road in the current event.
00905 {
00906 return m_Index;
00907 }
|
|
|
Intercept of trajectory.
|
|
|
Intercept of trajectory.
01055 { return m_SimpleIntercept; }
|
|
|
Which gap is the last one with hits attached to this road?
|
|
|
Which gap is the last one with hits attached to this road?
00934 { return m_LastGap; }
|
|
|
How many hits were attached in the gap with the most attached hits?
|
|
|
How many hits were attached in the gap with the most attached hits?
00938 { return m_MaxHitsPerGap; }
|
|
|
How many gaps provide hits attached to this road?
|
|
|
How many gaps provide hits attached to this road?
00947 {
00948 const int ngap = (int)MucID::getGapMax();
00949 int gap, count = 0;
00950 vector<int> firedGap;
00951 for ( gap = 0; gap < ngap; gap++) {
00952 firedGap.push_back(0);
00953 }
00954
00955 vector<MucRecHit*>::const_iterator iHit;
00956 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00957 if (*iHit) { // Check for a null pointer.
00958 gap = (*iHit)->Gap();
00959 firedGap[gap] = 1;
00960 }
00961 }
00962
00963 for ( gap = 0; gap < ngap; gap++) {
00964 count += firedGap[gap];
00965 }
00966
00967 return count;
00968 }
|
|
|
How many hits do two roads share?
|
|
|
How many hits do two roads share?
01022 {
01023 if (!road2) {
01024 return 0;
01025 }
01026
01027 int count = 0;
01028 vector<MucRecHit*>::const_iterator iHit1;
01029 vector<MucRecHit*>::const_iterator iHit2;
01030 MucRecHit *hit1, *hit2;
01031
01032 for( iHit1 = m_pHits.begin(); iHit1 != m_pHits.end(); iHit1++){
01033 for( iHit2 = road2->m_pHits.begin();
01034 iHit2 != road2->m_pHits.end(); iHit2++){
01035 hit1 = (*iHit1);
01036 hit2 = (*iHit2);
01037
01038 if ( (hit1 != 0) && (hit2 != 0) ) {
01039 if (hit1->GetID() == hit2->GetID()) {
01040 count++;
01041 }
01042 }
01043 }
01044 }
01045
01046 return count;
01047 }
|
|
|
In which orientation was this road found?
|
|
|
In which orientation was this road found?
00919 { return m_Orient; }
|
|
|
In which part was this road found?
|
|
|
In which part was this road found?
00911 { return m_Part; }
|
|
|
|
|
|
01095 {
01096 possigma = m_SimplePosSigma;
01097
01098 }
|
|
|
00160 {return m_QuadFitOK;}
|
|
|
00160 {return m_QuadFitOK;}
|
|
|
Chi-square parameter (per degree of freedom) of the trajectory fit.
|
|
|
Chi-square parameter (per degree of freedom) of the trajectory fit.
01064 {
01065 if ( (!m_SimpleFitOK) || (m_SimpleDOF < 0) ) {
01066 return -1.0;
01067 }
01068 else if (m_SimpleDOF == 0) {
01069 return 0.0;
01070 }
01071 else {
01072 return (m_SimpleChi2 / m_SimpleDOF);
01073 }
01074 }
|
|
|
Determine the size of the search window in the given gap.
|
|
|
Determine the size of the search window in the given gap.
01235 {
01236 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01237 cout << "MucRec2DRoad::GetSearchWindowSize-E1 invalid gap = " << gap << endl;
01238 return 0.0;
01239 }
01240
01241 // Determine the projection point of the "simple" fit to the last
01242 // gap and the desired gap.
01243 // FIXME? the "delta-x" variable is calculated as the scalar
01244 // difference between the positions obtained by projecting to the
01245 // last gap and to the desired gap.
01246
01247 MucGeoGeneral *geom = MucGeoGeneral::Instance();
01248 float x1, y1, z1, x2, y2, z2, dx, dy, dr;
01249 float sigmaX, sigmaY, sigmaZ;
01250
01251 if (m_Part == 0) {
01252 if (m_Orient == 0) {
01253 geom->FindIntersection(m_Part, 0, m_LastGap,
01254 1.0, m_SimpleSlope, 0.0,
01255 0.0, m_SimpleIntercept, 0.0,
01256 0.0, m_SimpleSlopeSigma, 0.0,
01257 0.0, m_SimpleInterceptSigma, 0.0,
01258 x1, y1, z1,
01259 sigmaX, sigmaY, sigmaZ);
01260 geom->FindIntersection(m_Part, 0, gap,
01261 1.0, m_SimpleSlope, 0.0,
01262 0.0, m_SimpleIntercept, 0.0,
01263 0.0, m_SimpleSlopeSigma, 0.0,
01264 0.0, m_SimpleInterceptSigma, 0.0,
01265 x2, y2, z2,
01266 sigmaX, sigmaY, sigmaZ);
01267 dx = z2 - z1;
01268 dy = sqrt(x2*x2 + y2*y2) - sqrt(x1*x1 + y1*y1);
01269 }
01270 else {
01271 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01272 m_SimpleSlope, 0.0, 1.0,
01273 m_SimpleIntercept, 0.0, 0.0,
01274 m_SimpleSlopeSigma, 0.0, 0.0,
01275 m_SimpleInterceptSigma, 0.0, 0.0,
01276 x1, y1, z1,
01277 sigmaX, sigmaY, sigmaZ);
01278 geom->FindIntersection(m_Part, m_Seg, gap,
01279 m_SimpleSlope, 0.0, 1.0,
01280 m_SimpleIntercept, 0.0, 0.0,
01281 m_SimpleSlopeSigma, 0.0, 0.0,
01282 m_SimpleInterceptSigma, 0.0, 0.0,
01283 x2, y2, z2,
01284 sigmaX, sigmaY, sigmaZ);
01285 dx = x2 - x1;
01286 dy = y2 - y1;
01287 }
01288 }
01289 else {
01290 if (m_Orient == 0) {
01291 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01292 0.0, m_SimpleSlope, 1.0,
01293 0.0, m_SimpleIntercept, 0.0,
01294 0.0, m_SimpleSlopeSigma, 0.0,
01295 0.0, m_SimpleInterceptSigma, 0.0,
01296 x1, y1, z1,
01297 sigmaX, sigmaY, sigmaZ);
01298 geom->FindIntersection(m_Part, m_Seg, gap,
01299 0.0, m_SimpleSlope, 1.0,
01300 0.0, m_SimpleIntercept, 0.0,
01301 0.0, m_SimpleSlopeSigma, 0.0,
01302 0.0, m_SimpleInterceptSigma, 0.0,
01303 x2, y2, z2,
01304 sigmaX, sigmaY, sigmaZ);
01305 dx = z2 - z1;
01306 dy = y2 - y1;
01307 }
01308 else {
01309 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01310 m_SimpleSlope, 0.0, 1.0,
01311 m_SimpleIntercept, 0.0, 0.0,
01312 m_SimpleSlopeSigma, 0.0, 0.0,
01313 m_SimpleInterceptSigma, 0.0, 0.0,
01314 x1, y1, z1,
01315 sigmaX, sigmaY, sigmaZ);
01316 geom->FindIntersection(m_Part, m_Seg, gap,
01317 m_SimpleSlope, 0.0, 1.0,
01318 m_SimpleIntercept, 0.0, 0.0,
01319 m_SimpleSlopeSigma, 0.0, 0.0,
01320 m_SimpleInterceptSigma, 0.0, 0.0,
01321 x2, y2, z2,
01322 sigmaX, sigmaY, sigmaZ);
01323 dx = z2 - z1;
01324 dy = x2 - x1;
01325 }
01326 }
01327
01328 dr = sqrt(dx*dx + dy*dy);
01329
01330 return WindowFunc(m_SimpleChi2,dr);
01331 }
|
|
|
In which segment was this road found?
|
|
|
In which segment was this road found?
00915 { return m_Seg; }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Get the parameters from the simple quad fit.
|
|
||||||||||||||||||||||||||||
|
Get the parameters from the simple fit.
|
|
||||||||||||||||||||||||||||||||||||||||
|
Get the parameters from the simple quad fit.
01106 {
01107 a = m_SimpleQuad_a;
01108 b = m_SimpleQuad_b;
01109 c = m_SimpleQuad_c;
01110 whichhalf = m_SimpleQuad_whichhalf;
01111
01112 sigmaa = m_SimpleQuad_aSigma;
01113 sigmab = m_SimpleQuad_bSigma;
01114 sigmac = m_SimpleQuad_cSigma;
01115
01116 chisq = m_SimpleChi2;
01117 ndof = m_SimpleDOF;
01118
01119 return;
01120 }
|
|
||||||||||||||||||||||||||||
|
Get the parameters from the simple fit.
01081 {
01082 slope = m_SimpleSlope;
01083 intercept = m_SimpleIntercept;
01084 sigmaSlope = m_SimpleSlopeSigma;
01085 sigmaIntercept = m_SimpleInterceptSigma;
01086 chisq = m_SimpleChi2;
01087 ndof = m_SimpleDOF;
01088
01089 return;
01090 }
|
|
|
Slope of trajectory.
|
|
|
Slope of trajectory.
01051 { return m_SimpleSlope; }
|
|
|
How many hits in all does this road contain?
|
|
|
How many hits in all does this road contain?
00942 { return m_TotalHits; }
|
|
||||||||||||||||
|
Position of the vertex point.
|
|
||||||||||||||||
|
Position of the vertex point.
00924 {
00925 x = m_VertexPos.x();
00926 y = m_VertexPos.y();
00927 z = m_VertexPos.z();
00928
00929 return;
00930 }
|
|
|
Does the hit exist in the road .
|
|
|
Does the hit exist in the road .
01391 {
01392
01393 vector<MucRecHit*>::const_iterator iHit;
01394 bool HitExist = false;
01395
01396 // Also, if a track hits overlap region of two panels, we avoid
01397 // to double count hits in two panels
01398
01399 Identifier id = hit->GetID();
01400
01401 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ ) {
01402 if ( *iHit ) { // Check for a null pointer.
01403 if ( (*iHit)->GetID() == id ) HitExist = true;
01404 }
01405 if (HitExist) break;
01406 }
01407
01408 return HitExist;
01409 }
|
|
|
Does this road contain any hits in the given gap?
|
|
|
Does this road contain any hits in the given gap?
00999 {
01000 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01001 cout << "MucRec2DRoad::HasHitInGap-E2 invalid gap = " << gap << endl;
01002 return false;
01003 }
01004
01005 bool found = false;
01006 vector<MucRecHit*>::const_iterator iHit;
01007
01008 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01009 if (*iHit) { // Check for a null pointer.
01010 if ( (*iHit)->Gap() == gap ) {
01011 found = true;
01012 }
01013 }
01014 }
01015
01016 return found;
01017 }
|
|
|
Assignment constructor.
|
|
|
Assignment constructor.
00048 {
00049 // Assignment operator.
00050 if ( this != &orig ) { // Watch out for self-assignment!
00051 m_VertexPos = orig.m_VertexPos;
00052 m_VertexSigma = orig.m_VertexSigma;
00053 m_Index = orig.m_Index;
00054 m_Part = orig.m_Part;
00055 m_Seg = orig.m_Seg;
00056 m_Orient = orig.m_Orient;
00057 m_Chi2 = orig.m_Chi2;
00058 m_DOF = orig.m_DOF;
00059 m_FitOK = orig.m_FitOK;
00060 m_MaxHitsPerGap = orig.m_MaxHitsPerGap;
00061 m_LastGap = orig.m_LastGap;
00062 m_TotalHits = orig.m_TotalHits;
00063 m_HitDistance = orig.m_HitDistance;
00064 m_pHits = orig.m_pHits;
00065 m_fittingMethod = orig.m_fittingMethod;
00066 }
00067
00068 return *this;
00069 }
|
|
|
Print Hits Infomation.
|
|
|
Print Hits Infomation.
01445 {
01446 vector<MucRecHit*>::const_iterator iHit;
01447 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01448 if (*iHit) { // Check for a null pointer.
01449 float xl, yl, zl;
01450 (*iHit)->GetStrip()->GetCenterPos(xl, yl, zl);
01451 HepPoint3D vl(xl, yl, zl);
01452 HepPoint3D vg = (*iHit)->GetGap()->TransformToGlobal(vl);
01453
01454 cout << " orient " << m_Orient
01455 << " part " << (*iHit)->Part()
01456 << " seg " << (*iHit)->Seg()
01457 << " gap " << (*iHit)->Gap()
01458 << " strip " << (*iHit)->Strip()
01459 << " pos (" << vg.x() << ", " << vg.y() << ", " << vg.z() << ")"
01460 << endl;
01461 }
01462 }
01463
01464 }
|
|
||||||||||||||||||||||||||||||||
|
Where does the trajectory of this road intersect a specific gap?
|
|
||||||||||||||||||||||||||||||||
|
Where does the trajectory of this road intersect a specific gap?
00731 {
00732 float sigmaX, sigmaY, sigmaZ;
00733
00734 x = 0.0; sigmaX = 0.0; x2 = 0.0;
00735 y = 0.0; sigmaY = 0.0; y2 = 0.0;
00736 z = 0.0; sigmaZ = 0.0; z2 = 0.0;
00737
00738 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00739 cout << "MucRec2DRoad::Project-E1 invalid gap = " << gap << endl;
00740 return;
00741 }
00742
00743 // Determine the projection point of the "simple" fit to the desired gap.
00744 float x0, y0, z0, sigmaX0, sigmaY0, sigmaZ0;
00745 float phi = m_Seg*0.25*kPi;
00746 MucGeoGeneral *geom = MucGeoGeneral::Instance();
00747
00748 if (m_Part == 1) {
00749 if (m_Orient == 0) {
00750 geom->FindIntersection(m_Part, m_Seg, gap,
00751 m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
00752 m_SimpleIntercept*cos(phi),m_SimpleIntercept*sin(phi), 0.0,
00753 m_SimpleSlopeSigma, 0.0, 0.0,
00754 m_SimpleInterceptSigma, 0.0, 0.0,
00755 x0, y0, z0,
00756 sigmaX0, sigmaY0, sigmaZ0);
00757
00758 if(m_SimpleSlope>10000){ //the line is in right center of this segment, we can not get intersection in y coordinate, in this case, m_SimpleIntercept is in z coordinate.
00759 geom->FindIntersection(m_Part, m_Seg, gap,
00760 m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
00761 0.0, 0.0, m_SimpleIntercept,
00762 m_SimpleSlopeSigma, 0.0, 0.0,
00763 m_SimpleInterceptSigma, 0.0, 0.0,
00764 x0, y0, z0,
00765 sigmaX0, sigmaY0, sigmaZ0);
00766
00767 }
00768
00769 }
00770 else {
00771 geom->FindIntersection(m_Part, m_Seg, gap,
00772 1.0, m_SimpleSlope, 0.0,
00773 0.0, m_SimpleIntercept, 0.0,
00774 0.0, m_SimpleSlopeSigma, 0.0,
00775 0.0, m_SimpleInterceptSigma, 0.0,
00776 x0, y0, z0,
00777 sigmaX0, sigmaY0, sigmaZ0);
00778 //cout<<"in MucRec2DRoad line Project xyz0 = "<<x0<<" "<<y0<<" "<<z0<<endl;
00779
00780 }
00781 }
00782 else {
00783 if (m_Orient == 0) {
00784 geom->FindIntersection(m_Part, m_Seg, gap,
00785 m_SimpleSlope, m_SimpleSlope, 1.0,
00786 0.0, m_SimpleIntercept, 0.0,
00787 0.0, m_SimpleSlopeSigma, 0.0,
00788 0.0, m_SimpleInterceptSigma, 0.0,
00789 x0, y0, z0,
00790 sigmaX0, sigmaY0, sigmaZ0);
00791 }
00792 else {
00793 geom->FindIntersection(m_Part, m_Seg, gap,
00794 m_SimpleSlope, m_SimpleSlope, 1.0,
00795 m_SimpleIntercept, 0.0, 0.0,
00796 m_SimpleSlopeSigma, 0.0, 0.0,
00797 m_SimpleInterceptSigma, 0.0, 0.0,
00798 x0, y0, z0,
00799 sigmaX0, sigmaY0, sigmaZ0);
00800 }
00801 //cout << " part " << m_Part
00802 // << " seg " << m_Seg
00803 // << " gap " << gap
00804 // << " orient " << m_Orient
00805 // << " slope = " << m_SimpleSlope
00806 // << endl;
00807 }
00808
00809 //cout << "In find intersection x0 = " << x0 << " y0 = " << y0 << " z0 = " << z0 << endl;
00810
00811 float a,b,sa,sb,chi2;
00812 int ndof;
00813
00814 int status = Fit(x0, y0, z0, a, b, sa, sb, chi2, ndof);
00815
00816 // m_FitOK = (status == 0) && (chi2<1000.0);
00817 // if (!fFitOK) {
00818 // cout << "MucRec2DRoad::Project-E2 fit fail status = "
00819 // << status << " npoints = " << npoints << " chi-sq = "
00820 // << chi2 << endl;
00821 // }
00822
00823 // // Assign to fields of TMuiRoad object.
00824 // m_DOF = npoints - 2;
00825 // m_Chi2 = chi2;
00826
00827 if (m_Part == 1) { //change from 0 to 1 at 2006.11.30
00828 if (m_Orient == 0) {
00829 geom->FindIntersection(m_Part, m_Seg, gap,
00830 a*cos(phi), a*sin(phi), 1.0,
00831 //a, 0.0, 1.0,
00832 b*cos(phi), b*sin(phi), 0.0,
00833 sa, 0.0, 0.0,
00834 sb, 0.0, 0.0,
00835 x, y, z,
00836 sigmaX, sigmaY, sigmaZ);
00837
00838 if(fabs(a)>10000){
00839 geom->FindIntersection(m_Part, m_Seg, gap,
00840 a*cos(phi), a*sin(phi), 1.0,
00841 0.0, 0.0, b,
00842 //a, 0.0, 1.0,
00843 //b, 0.0, 0.0,
00844 sa, 0.0, 0.0,
00845 sb, 0.0, 0.0,
00846 x, y, z,
00847 sigmaX, sigmaY, sigmaZ);
00848
00849 }
00850 }
00851 else {
00852 geom->FindIntersection(m_Part, m_Seg, gap,
00853 1.0, a, 0.0,
00854 0.0, b, 0.0,
00855 0.0, sa, 0.0,
00856 0.0, sb, 0.0,
00857 x, y, z,
00858 sigmaX, sigmaY, sigmaZ);
00859
00860 if(m_fittingMethod == 2 && m_QuadFitOK ){
00861 float a, b, c;
00862 float sa, sb, sc, chi2x; int ydof; int whichhalf;
00863
00864 GetSimpleFitParams(a, b, c, whichhalf, sa, sb, sc, chi2x, ydof);
00865 geom->FindIntersection(m_Part, m_Seg, gap,
00866 10E30, 0.0, m_SimpleIntercept, 0.0, //vy = infinite
00867 a, b, c, //y = a*x*x + b*x +c
00868 whichhalf,
00869 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00870 x, y, z, x2, y2, z2,
00871 sigmaX, sigmaY, sigmaZ);
00872
00873
00874 }
00875
00876 }
00877 }
00878 else {
00879 if (m_Orient == 0) {
00880 geom->FindIntersection(m_Part, m_Seg, gap,
00881 a, a, 1.0,
00882 0.0, b, 0.0,
00883 0.0, sa, 0.0,
00884 0.0, sb, 0.0,
00885 x, y, z,
00886 sigmaX, sigmaY, sigmaZ);
00887 }
00888 else {
00889 geom->FindIntersection(m_Part, m_Seg, gap,
00890 a, a, 1.0,
00891 b, 0.0, 0.0,
00892 sa, 0.0, 0.0,
00893 sb, 0.0, 0.0,
00894 x, y, z,
00895 sigmaX, sigmaY, sigmaZ);
00896 }
00897 }
00898
00899 return;
00900 }
|
|
|
Set the index for this road.
|
|
|
Set the index for this road.
00093 {
00094 if (index >= 0) m_Index = index;
00095 }
|
|
|
Max number of consecutive gaps allowed with no hits attached. This parameter affects the calculation of the last gap. |
|
|
Max number of consecutive gaps allowed with no hits attached. This parameter affects the calculation of the last gap. 00181 {
00182 m_MaxNSkippedGaps = numGaps;
00183 CountHits(); // recalculate the last gap and the hit counts.
00184 }
|
|
||||||||||||||||||||||||||||
|
Calculate the best-fit straight line with "simple" weights.
|
|
||||||||||||||||||||||||||||
|
Calculate the best-fit straight line with "simple" weights.
00194 {
00195 // Assign to temporary arrays to be used in fit.
00196 float px[100];
00197 float py[100];
00198 float pw[100];
00199 int npoints = 0;
00200
00201 vector<MucRecHit*>::const_iterator iHit;
00202
00203 float weight[100];
00204
00205 // for (int i = 0; i < m_pHits.size(); i++) {
00206 // cout<<"info: "<<m_pHits[i]->Seg()<<" "<<m_pHits[i]->Gap()<<" "<< m_pHits[i]->Strip()<<endl;
00207 // }
00208 for (int i = 0; i < m_pHits.size(); i++) {
00209 weight[i] = 1;
00210
00211 for(int j = 0; j < m_pHits.size(); j++){
00212
00213 if(j == i) continue;
00214 if(m_pHits[i]->Part() == m_pHits[j]->Part() &&
00215 m_pHits[i]->Seg() == m_pHits[j]->Seg() &&
00216 m_pHits[i]->Gap() == m_pHits[j]->Gap() )
00217 {
00218 int deltaStrip = fabs(m_pHits[i]->Strip()- m_pHits[j]->Strip());
00219
00220 //cout<<i<<" "<<m_pHits[i]->Seg()<<" "<<m_pHits[i]->Gap()<<" "<<m_pHits[i]->Strip()<<" - "<<m_pHits[j]->Strip()<<endl;
00221 if(deltaStrip == 0){
00222 cout<<"deltaStrip == 0 ? check it"<<endl;
00223 }
00224 else{
00225 weight[i] *= (deltaStrip+1) * (deltaStrip+1);
00226 }
00227
00228 }
00229 }
00230 }
00231
00232 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00233 if (*iHit) { // Check for a null pointer.
00234
00235 /*
00236 float localx, localy, localz;
00237 (*iHit)->GetStrip()->GetCenterPos(localx, localy, localz);
00238 if ( m_Orient == 0) {
00239 px[npoints] = localy;
00240 py[npoints] = localz;
00241 }
00242 else {
00243 px[npoints] = localx;
00244 py[npoints] = localz;
00245 }
00246 npoints++;
00247 }
00248 }
00249 */
00250
00251
00252 Hep3Vector center = (*iHit)->GetCenterPos();
00253 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00254 //Hep3Vector sigma(1.0, 1.0, 1.0);
00255
00256 if (m_Part == 1) {
00257 if ( m_Orient == 0) {
00258 px[npoints] = center.z();
00259 py[npoints] = sqrt(center.x()*center.x()
00260 + center.y()*center.y());
00261 if(m_Seg==2) py[npoints] = center.y(); //deal with seg2 seperately! because there is a hole here. 2006.11.9
00262
00263 pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
00264 }
00265 else {
00266 px[npoints] = center.x();
00267 py[npoints] = center.y();
00268 pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
00269 }
00270 }
00271 else {
00272 if ( m_Orient == 0) {
00273 px[npoints] = center.z();
00274 py[npoints] = center.y();
00275 pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
00276 }
00277 else {
00278 px[npoints] = center.z();
00279 py[npoints] = center.x();
00280 pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
00281 }
00282 }
00283
00284 //cout << " in MucRec2DRoad ID: " <<(*iHit)->Part()<<" "<<(*iHit)->Seg()<<" "<<(*iHit)->Gap()<<" "<< (*iHit)->Strip()<<" "<<px[npoints] << " " << py[npoints] << " " << pw[npoints] << endl;
00285 npoints++;
00286
00287 if(npoints > 99)
00288 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00289 return -1;
00290 }
00291 }
00292 }
00293 /*
00294 float px2[100];
00295 float py2[100];
00296 for (int i = 0; i < m_pHits.size(); i++) {
00297 int hitsInGap = 1;
00298 px2[i] = -999; py2[i] = -999;
00299 for(int j = 0; j < m_pHits.size(); j++){
00300
00301 if(j == i) continue;
00302 if(m_pHits[i]->Part() == m_pHits[j]->Part() &&
00303 m_pHits[i]->Seg() == m_pHits[j]->Seg() &&
00304 m_pHits[i]->Gap() == m_pHits[j]->Gap() )
00305 {
00306 hitsInGap++;
00307 px2[i] = (px[i]*(hitsInGap-1) + px[j])/hitsInGap;
00308 py2[i] = (py[i]*(hitsInGap-1) + py[j])/hitsInGap;
00309 cout<<hitsInGap<<" "<<px[i]<<" "<<px[j]<<" "<<px2[i]<<endl;
00310 }
00311 }
00312 }
00313
00314 for(int i = 0; i < m_pHits.size(); i++){
00315 if(px2[i] != -999&&py2[i] != -999){
00316 px[i] = px2[i]; py[i] = py2[i];
00317 }
00318 }
00319 */
00320
00321 ndof = npoints - 2;
00322 if (ndof < 0) return -1;
00323
00324 if (npoints == 1) {
00325 if (m_Part == 1) {
00326 if ( m_Orient == 0) {
00327 px[npoints] = m_VertexPos.z();
00328 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00329 + m_VertexPos.y()*m_VertexPos.y());
00330 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00331 }
00332 else {
00333 px[npoints] = m_VertexPos.x();
00334 py[npoints] = m_VertexPos.y();
00335 }
00336 }
00337 else {
00338 if ( m_Orient == 0) {
00339 px[npoints] = m_VertexPos.z();
00340 py[npoints] = m_VertexPos.y();
00341 }
00342 else {
00343 px[npoints] = m_VertexPos.z();
00344 py[npoints] = m_VertexPos.x();
00345 }
00346 }
00347 pw[npoints] = 1.0;
00348 npoints++;
00349 }
00350 else {
00351 if (npoints == 0 ) {
00352 return -1;
00353 }
00354 }
00355
00356 // Do the fits here.
00357 MucRecLineFit fit;
00358 int status = fit.LineFit(px, py, pw, npoints,
00359 &slope, &intercept, &chisq,
00360 &sigmaSlope, &sigmaIntercept);
00361
00362
00363 float tempslope, tempintercept,tempchisq, tempsigmaslope, sigmaPos;
00364 int status4 = fit.LineFit(px, py, pw,m_Part,m_Seg,m_Orient, npoints,
00365 &tempslope, &tempintercept, &tempchisq,
00366 &tempsigmaslope, &sigmaPos);
00367
00368 MucRecQuadFit quadfit;
00369 float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
00370 int whichhalf, status2;
00371
00372 if(m_fittingMethod == 2){
00373 status2 = quadfit.QuadFit(px, py, pw, npoints,
00374 &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
00375 &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
00376
00377 }
00378 //cout << " in MucRec2DRoad slope " << slope << " "<<intercept<<endl;
00379
00380 if (slope > 1.0e10 || slope < -1.0e10) {
00381 if (m_Seg > 4) slope *= -1.0; // to avoid wrong direction
00382 }
00383
00384 m_SimpleSlope = slope;
00385 m_SimpleSlopeSigma = sigmaSlope;
00386 m_SimpleIntercept = intercept;
00387 m_SimpleInterceptSigma = sigmaIntercept;
00388 m_SimplePosSigma = sigmaPos; //new 20071227
00389 m_SimpleChi2 = chisq;
00390 m_SimpleDOF = ndof;
00391 m_SimpleFitOK = (status == 0) && (chisq < 1000.0);
00392 m_QuadFitOK = (status2 == 1);
00393
00394 m_SimpleQuad_a = quad_a;
00395 m_SimpleQuad_b = quad_b;
00396 m_SimpleQuad_c = quad_c;
00397 m_SimpleQuad_whichhalf = whichhalf;
00398 m_SimpleQuad_aSigma = sigmaquad_a;
00399 m_SimpleQuad_bSigma = sigmaquad_b;
00400 m_SimpleQuad_cSigma = sigmaquad_c;
00401
00402 return status;
00403 }
|
|
||||||||||||||||||||||||||||||||
|
Calculate the best-fit straight line with "simple" weights. not use current gap!!!
|
|
||||||||||||||||||||||||||||||||
|
Calculate the best-fit straight line with "simple" weights. not use current gap!!!
00417 {
00418 // Assign to temporary arrays to be used in fit.
00419 float px[100];
00420 float py[100];
00421 float pw[100];
00422 int npoints = 0;
00423
00424 int notused = 0;
00425
00426 vector<MucRecHit*>::const_iterator iHit;
00427
00428 float pw2[100];
00429 for(int i = 0; i< m_pHits.size(); i++) {
00430 pw2[i] = 1;
00431 //----if( m_pHits[i]->Gap()==currentgap ) pw2[i] = 9999;
00432 }
00433
00434 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00435 if (*iHit) { // Check for a null pointer.
00436 if( (*iHit)->Gap()==currentgap ) continue;
00437
00438 Hep3Vector center = (*iHit)->GetCenterPos();
00439 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00440 //Hep3Vector sigma(1.0, 1.0, 1.0);
00441
00442 if (m_Part == 1) {
00443 if ( m_Orient == 0) {
00444 px[npoints] = center.z();
00445 py[npoints] = sqrt(center.x()*center.x()
00446 + center.y()*center.y());
00447 if(m_Seg==2) py[npoints] = center.y(); //deal with seg2 seperately! because there is a hole here. 2006.11.9
00448
00449 pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
00450 }
00451 else {
00452 px[npoints] = center.x();
00453 py[npoints] = center.y();
00454 pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
00455 }
00456 }
00457 else {
00458 if ( m_Orient == 0) {
00459 px[npoints] = center.z();
00460 py[npoints] = center.y();
00461 pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
00462 }
00463 else {
00464 px[npoints] = center.z();
00465 py[npoints] = center.x();
00466 pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
00467 }
00468 }
00469
00470 //cout << " in MucRec2DRoad ID: " <<(*iHit)->Part()<<" "<<(*iHit)->Seg()<<" "<<(*iHit)->Gap()<<" "<< (*iHit)->Strip()<<" "<<px[npoints] << " " << py[npoints] << " " << pw[npoints] << endl;
00471 //cout<<"process ngap: "<<currentgap<<" current: "<<(*iHit)->Gap()<<" x: "<<px[npoints]<<" y: "<<py[npoints]<<" weight: "<<pw[npoints]<<endl;
00472 npoints++;
00473
00474 if(npoints > 99)
00475 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00476 return -1;
00477 }
00478 }
00479 }
00480
00481 ndof = npoints - 2;
00482 if (ndof < 0) return -1;
00483
00484 if (npoints == 1) {
00485 if (m_Part == 1) {
00486 if ( m_Orient == 0) {
00487 px[npoints] = m_VertexPos.z();
00488 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00489 + m_VertexPos.y()*m_VertexPos.y());
00490 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00491 }
00492 else {
00493 px[npoints] = m_VertexPos.x();
00494 py[npoints] = m_VertexPos.y();
00495 }
00496 }
00497 else {
00498 if ( m_Orient == 0) {
00499 px[npoints] = m_VertexPos.z();
00500 py[npoints] = m_VertexPos.y();
00501 }
00502 else {
00503 px[npoints] = m_VertexPos.z();
00504 py[npoints] = m_VertexPos.x();
00505 }
00506 }
00507 pw[npoints] = 1.0;
00508 npoints++;
00509 }
00510 else {
00511 if (npoints == 0 ) {
00512 return -1;
00513 }
00514
00515 }
00516
00517 // Do the fits here.
00518 MucRecLineFit fit;
00519 int status = fit.LineFit(px, py, pw, npoints,
00520 &slope, &intercept, &chisq,
00521 &sigmaSlope, &sigmaIntercept);
00522
00523 MucRecQuadFit quadfit;
00524 float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
00525 int whichhalf, status2;
00526
00527 if(m_fittingMethod == 2){
00528 status2 = quadfit.QuadFit(px, py, pw, npoints,
00529 &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
00530 &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
00531
00532 }
00533 //cout << " in MucRec2DRoad slope " << slope << " "<<intercept<<endl;
00534
00535 if (slope > 1.0e10 || slope < -1.0e10) {
00536 if (m_Seg > 4) slope *= -1.0; // to avoid wrong direction
00537 }
00538
00540 m_SimpleSlope = slope;
00541 m_SimpleSlopeSigma = sigmaSlope;
00542 m_SimpleIntercept = intercept;
00543 m_SimpleInterceptSigma = sigmaIntercept;
00544 //m_SimplePosSigma = sigmaPos; //new 20071227
00545 m_SimpleChi2 = chisq;
00546 m_SimpleDOF = ndof;
00547 m_SimpleFitOK = (status == 0) && (chisq < 1000.0);
00548 m_QuadFitOK = (status2 == 1);
00549
00550 m_SimpleQuad_a = quad_a;
00551 m_SimpleQuad_b = quad_b;
00552 m_SimpleQuad_c = quad_c;
00553 m_SimpleQuad_whichhalf = whichhalf;
00554 m_SimpleQuad_aSigma = sigmaquad_a;
00555 m_SimpleQuad_bSigma = sigmaquad_b;
00556 m_SimpleQuad_cSigma = sigmaquad_c;
00557
00558
00559 return status;
00560 }
|
|
||||||||||||
|
|
|
||||||||||||
|
01469 {
01470 return 1.0;
01471 }
|
|
||||||||||||
|
|
|
||||||||||||
|
01476 {
01477 return 1.0;
01478 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1