Remove erroneous chi2 calculation causing FPEs

Author: martenole

GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx

@@ -731,7 +731,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
             float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp);
             // GPUInfo("layer %i: chi2 = %f", iLayer, chi2);
             if (chi2 < mMaxChi2 && CAMath::Abs(GetAngularPull(mSpacePoints[trkltIdx].mDy, trkWork->getSnp())) < 4) {
-              Hypothesis hypo(trkWork->GetNlayers(), iCandidate, trkltIdx, trkWork->GetChi2() + chi2, GetPredictedChi2(trkltPosTmpYZ, trkltCovTmp, trkWork->getPar(), trkWork->getCov()));
+              Hypothesis hypo(trkWork->GetNlayers(), iCandidate, trkltIdx, trkWork->GetChi2() + chi2);
               InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset);
             } // end tracklet chi2 < mMaxChi2
           }   // end tracklet in window
@@ -769,7 +769,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
         My_Float covReal[3] = {0.};
         RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(mTracklets[realTrkltId].GetZbin()), covReal);
         mDebug->SetChi2Real(prop->getPredictedChi2(yzPosReal, covReal), iLayer);
-        mDebug->SetChi2YZPhiReal(GetPredictedChi2(yzPosReal, covReal, trkWork->getPar(), trkWork->getCov()), iLayer);
         mDebug->SetRawTrackletPositionReal(mSpacePoints[realTrkltId].mR, mSpacePoints[realTrkltId].mX, iLayer);
         mDebug->SetCorrectedTrackletPositionReal(yzPosReal, iLayer);
         mDebug->SetTrackletPropertiesReal(mTracklets[realTrkltId].GetDetector(), iLayer);
@@ -778,7 +777,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
 #endif
     //
     mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->GetChi2(), iLayer); // only meaningful for ONE candidate!!!
-    mDebug->SetChi2YZPhiUpdate(mHypothesis[0 + hypothesisIdxOffset].mChi2YZPhi, iLayer);      // only meaningful for ONE candidate!!!
     mDebug->SetRoad(roadY, roadZ, iLayer);                                                    // only meaningful for ONE candidate
     bool wasTrackStored = false;
     // --------------------------------------------------------------------------------
@@ -989,40 +987,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
   return true;
 }
 
-template <class TRDTRK, class PROP>
-GPUd() float GPUTRDTracker_t<TRDTRK, PROP>::GetPredictedChi2(const My_Float* pTRD, const My_Float* covTRD, const My_Float* pTrk, const My_Float* covTrk) const
-{
-  // FIXME: This function cannot yet be used in production (maybe it is not necessary at all?)
-  // conversion from tracklet deflection to sin(phi) needs to be parametrized
-  // predict chi2 for update of track with pTrk and covTrk with TRD space point with pTRD and covTRD
-  // taking into account y, z and azimuthal angle of the track and tracklet
-  float deltaY = pTrk[0] - pTRD[0];
-  float deltaZ = pTrk[1] - pTRD[1];
-  float deltaS = pTrk[2] - pTRD[2]; // FIXME: pTRD[2] is not defined, needs conversion dy -> sin(phi)
-
-  // add errors for track and space point, assume no correlation in y-sin(phi) and z-sin(phi) for space point
-  float sigmaZ2 = covTrk[2] + covTRD[2];
-  float sigmaS2 = covTrk[5] + GetAngularResolution(pTrk[2]); // FIXME: convert angular resolution from dy to sin(phi) for the TRD space point
-  float sigmaY2 = covTrk[0] + covTRD[0];
-  float sigmaZS = covTrk[4];
-  float sigmaYS = covTrk[3];
-  float sigmaYZ = covTrk[1] + covTRD[1];
-  // inverse of the covariance matrix
-  float c11 = sigmaZ2 * sigmaS2 - sigmaZS * sigmaZS;
-  float c21 = sigmaZS * sigmaYS - sigmaYZ * sigmaS2;
-  float c22 = sigmaY2 * sigmaS2 - sigmaYS * sigmaYS;
-  float c31 = sigmaYZ * sigmaZS - sigmaZ2 * sigmaYS;
-  float c32 = sigmaYZ * sigmaYS - sigmaY2 * sigmaZS;
-  float c33 = sigmaY2 * sigmaZ2 - sigmaYZ * sigmaYZ;
-  // determinant
-  float det = sigmaY2 * sigmaZ2 * sigmaS2 + 2 * sigmaYZ * sigmaZS * sigmaYS - sigmaYS * sigmaYS * sigmaZ2 - sigmaZS * sigmaZS * sigmaY2 - sigmaS2 * sigmaYZ * sigmaYZ;
-  if (CAMath::Abs(det) < 1.e-10f) {
-    printf("Determinant too small: %f\n", det);
-    det = 1.e-10f;
-  }
-  det = 1.f / det;
-  return (c11 * deltaY * deltaY + 2.f * c21 * deltaY * deltaZ + 2.f * c31 * deltaY * deltaS + c22 * deltaZ * deltaZ + 2.f * c32 * deltaZ * deltaS + c33 * deltaS * deltaS) * det;
-}
 
 template <class TRDTRK, class PROP>
 GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::InsertHypothesis(Hypothesis hypo, int& nCurrHypothesis, int idxOffset)

GPU/GPUTracking/TRDTracking/GPUTRDTracker.h

@@ -95,11 +95,10 @@ class GPUTRDTracker_t : public GPUProcessor
     int mCandidateId; // to which track candidate the hypothesis belongs
     int mTrackletId;  // tracklet index to be used for update
     float mChi2;      // predicted chi2 for given space point
-    float mChi2YZPhi; // not yet ready (see GetPredictedChi2 method in cxx file)
 
     GPUd() float GetReducedChi2() { return mLayers > 0 ? mChi2 / mLayers : mChi2; }
     GPUd() Hypothesis() : mLayers(0), mCandidateId(-1), mTrackletId(-1), mChi2(9999.f) {}
-    GPUd() Hypothesis(int layers, int candidateId, int trackletId, float chi2, float chi2YZPhi = -1.f) : mLayers(layers), mCandidateId(candidateId), mTrackletId(trackletId), mChi2(chi2), mChi2YZPhi(chi2YZPhi) {}
+    GPUd() Hypothesis(int layers, int candidateId, int trackletId, float chi2, float chi2YZPhi = -1.f) : mLayers(layers), mCandidateId(candidateId), mTrackletId(trackletId), mChi2(chi2) {}
   };
 
   short MemoryPermanent() const { return mMemoryPermanent; }
@@ -149,7 +148,6 @@ class GPUTRDTracker_t : public GPUProcessor
   GPUd() void DoTrackingThread(int iTrk, int threadId = 0);
   GPUd() bool CalculateSpacePoints(int iCollision = 0);
   GPUd() bool FollowProlongation(PROP* prop, TRDTRK* t, int threadId, int collisionId);
-  GPUd() float GetPredictedChi2(const My_Float* pTRD, const My_Float* covTRD, const My_Float* pTrk, const My_Float* covTrk) const;
   GPUd() int GetDetectorNumber(const float zPos, const float alpha, const int layer) const;
   GPUd() bool AdjustSector(PROP* prop, TRDTRK* t, const int layer) const;
   GPUd() int GetSector(float alpha) const;

GPU/GPUTracking/TRDTracking/GPUTRDTrackerDebug.h

@@ -95,9 +95,7 @@ class GPUTRDTrackerDebug
     fTrackYReal.ResizeTo(6);
     fTrackZReal.ResizeTo(6);
     fTrackSecReal.ResizeTo(6);
-    fChi2YZPhiUpdate.ResizeTo(6);
     fChi2Update.ResizeTo(6);
-    fChi2YZPhiReal.ResizeTo(6);
     fChi2Real.ResizeTo(6);
     fNmatchesAvail.ResizeTo(6);
     fFindable.ResizeTo(6);
@@ -151,9 +149,7 @@ class GPUTRDTrackerDebug
     fTrackYReal.Zero();
     fTrackZReal.Zero();
     fTrackSecReal.Zero();
-    fChi2YZPhiUpdate.Zero();
     fChi2Update.Zero();
-    fChi2YZPhiReal.Zero();
     fChi2Real.Zero();
     fNmatchesAvail.Zero();
     fFindable.Zero();
@@ -290,8 +286,6 @@ class GPUTRDTrackerDebug
   // update information
   void SetChi2Update(float chi2, int ly) { fChi2Update(ly) = chi2; }
   void SetChi2Real(float chi2, int ly) { fChi2Real(ly) = chi2; }
-  void SetChi2YZPhiUpdate(float chi2, int ly) { fChi2YZPhiUpdate(ly) = chi2; }
-  void SetChi2YZPhiReal(float chi2, int ly) { fChi2YZPhiReal(ly) = chi2; }
 
   // other infos
   void SetRoad(float roadY, float roadZ, int ly)
@@ -370,8 +364,6 @@ class GPUTRDTrackerDebug
       "trackletDetReal.=" << &fTrackletDetReal <<          // detector number for matching or related tracklet if available, otherwise -1
       "chi2Update.=" << &fChi2Update <<                    // chi2 for update
       "chi2Real.=" << &fChi2Real <<                        // chi2 for first tracklet w/ matching MC label
-      "chi2YZPhiUpdate.=" << &fChi2YZPhiUpdate <<          // chi2 for update taking into account full tracklet information (y, z, dY aka sin(phi))
-      "chi2YZPhiReal.=" << &fChi2YZPhiReal <<              // chi2 for first tracklet w/ matching MC label taking into account full tracklet information (y, z, dY aka sin(phi))
       "chi2Total=" << fChi2 <<                             // total chi2 for track
       "nLayers=" << fNlayers <<                            // number of layers in which track was findable
       "nTracklets=" << fNtrklts <<                         // number of attached tracklets
@@ -456,8 +448,6 @@ class GPUTRDTrackerDebug
   TVectorF fTrackletDetReal;
   TVectorF fChi2Update;
   TVectorF fChi2Real;
-  TVectorF fChi2YZPhiUpdate;
-  TVectorF fChi2YZPhiReal;
   TVectorF fRoadY;
   TVectorF fRoadZ;
   TVectorF fFindable;
@@ -512,7 +502,6 @@ class GPUTRDTrackerDebug
   GPUd() void SetChi2Update(float chi2, int ly) {}
   GPUd() void SetChi2Real(float chi2, int ly) {}
   GPUd() void SetChi2YZPhiUpdate(float chi2, int ly) {}
-  GPUd() void SetChi2YZPhiReal(float chi2, int ly) {}
 
   // other infos
   GPUd() void SetRoad(float roadY, float roadZ, int ly) {}