some fixes and parameter updates

Author: glegras

Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h

@@ -26,8 +26,10 @@ namespace trd
 struct TRDCalibParams : public o2::conf::ConfigurableParamHelper<TRDCalibParams> {
   unsigned int nTrackletsMin = 5;  ///< minimum amount of tracklets
   unsigned int nTrackletsMinLoose = 4; ///< minimum amount of tracklets if two layers with a large lever arm both have a hit
-  unsigned int chi2RedMax = 6;     ///< maximum reduced chi2 acceptable for track quality
-  size_t minEntriesChamber = 200;  ///< minimum number of entries per chamber to fit single time slot
+  unsigned int chi2RedMax = 6;      ///< maximum reduced chi2 acceptable for track quality
+  float minPtCalib = 1.;           ///< min pt for vd and ExB calib
+  bool rejectTPCTRD = true;         ///< reject TPC-TRD tracks for vd ExB calib
+  size_t minEntriesChamber = 200;   ///< minimum number of entries per chamber to fit single time slot
   size_t minEntriesTotal = 400'000; ///< minimum total required for meaningful fits
 
   // For gain calibration

Detectors/TRD/calibration/src/TrackBasedCalib.cxx

@@ -22,6 +22,7 @@
 #include "TRDBase/Geometry.h"
 #include "TRDBase/PadPlane.h"
 #include "CommonUtils/NameConf.h"
+#include "CommonConstants/GeomConstants.h"
 #include "DataFormatsTPC/TrackTPC.h"
 #include "ReconstructionDataFormats/TrackTPCITS.h"
 #include <fairlogger/Logger.h>
@@ -77,8 +78,9 @@ void TrackBasedCalib::calculateAngResHistos()
 
   LOGF(info, "As input tracks are available: %lu ITS-TPC-TRD tracks and %lu TPC-TRD tracks", mTracksInITSTPCTRD.size(), mTracksInTPCTRD.size());
 
+  auto& params = TRDCalibParams::Instance();
   int nTracksSuccessITSTPCTRD = doTrdOnlyTrackFits(mTracksInITSTPCTRD);
-  int nTracksSuccessTPCTRD = doTrdOnlyTrackFits(mTracksInTPCTRD);
+  int nTracksSuccessTPCTRD = params.rejectTPCTRD ? 0 : doTrdOnlyTrackFits(mTracksInTPCTRD);
 
   LOGF(info, "Successfully processed %i tracks (%i from ITS-TPC-TRD and %i from TPC-TRD) and collected %lu angular residuals",
        nTracksSuccessITSTPCTRD + nTracksSuccessTPCTRD, nTracksSuccessITSTPCTRD, nTracksSuccessTPCTRD, mAngResHistos.getNEntries());
@@ -195,6 +197,11 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
         continue;
       }
     }
+    if (trkIn.getPt() < params.minPtCalib) {
+      // we reject low pt tracks which might suffer from multiple scattering (giving lower quality of the TRD-only fit)
+      continue;
+    }
+
     auto trkWork = trkIn; // input is const, so we need to create a copy
     bool trackFailed = false;
 
@@ -205,6 +212,12 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       LOG(alarm) << "Track with invalid parameters found: " << trkWork.getRefGlobalTrackId();
       continue;
     }
+    
+    // reject tracks which cross sectors within TRD (if the extrapolation from the outer TRD to the outer TPC leads to a change in sector or close to the sector edges with 5 cm  margin), which have larger uncertainties and probably more fakes
+    float yOuterTPC = trkIn.getOuterParam().getYAt(o2::constants::geom::XTPCOuterRef, bz);
+    if (std::fabs(yOuterTPC) > o2::constants::geom::XTPCOuterRef * tan(M_PI / 18.) - 5.) {
+      continue;
+    }
 
     // first inward propagation (TRD track fit)
     int currLayer = NLAYER;
@@ -274,6 +287,7 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trkWork.getZ()) < padLength))) {
         tiltCorrUp = 0.f;
       }
+
       // use uncalibrated dy because online calibration does not work otherwise
       float trkltDy = mTrackletsRaw[trkltId].getUncalibratedDy(30.f / o2::trd::constants::VDRIFTDEFAULT) + tiltCorrUp;
       float trkltAngle = o2::math_utils::atan(trkltDy / Geometry::cdrHght()) * TMath::RadToDeg();

Detectors/TRD/qc/src/Tracking.cxx

@@ -165,13 +165,6 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
     float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp())) / std::sqrt(mRecoParam.getDyRes(trk.getSnp(), 0));
     // Correction of y position based on angular pull
     float corrPull = - angularPull * mRecoParam.getCorrYDy(trk.getSnp());
-    // in the tails with very large angle difference, the correlation becomes flat
-    if (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp()) > 0.6) {
-      corrPull = - 0.6 / mRecoParam.getDyRes(trk.getSnp(), 0) * mRecoParam.getCorrYDy(trk.getSnp());
-    }
-    if (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp()) < -0.6) {
-      corrPull = 0.6 / mRecoParam.getDyRes(trk.getSnp(), 0) * mRecoParam.getCorrYDy(trk.getSnp());
-    }
     yPosCorrUp += corrPull;
 
     std::array<float, 2> trkltPosUp{yPosCorrUp, zPosCorrUp};

Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx

@@ -934,13 +934,6 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
     // Correction of y position based on angular pull
     if (mRec->GetParam().rec.trd.useAngularPull == 3 || mRec->GetParam().rec.trd.useAngularPull == 4) {
       float corrPull = - angularPull * mRecoParam.getCorrYDy(trkParam->getSnp());
-      // in the tails with very large angle difference, the correlation becomes flat
-      if (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp()) > 0.6) {
-        corrPull = - 0.6 / mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber) * mRecoParam.getCorrYDy(trkParam->getSnp());
-      }
-      if (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp()) < -0.6) {
-        corrPull = 0.6 / mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber) * mRecoParam.getCorrYDy(trkParam->getSnp());
-      }
       yPosCorrUp += corrPull;
     }
 

GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx

@@ -45,18 +45,18 @@ void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec)
   if (CAMath::Abs(CAMath::Abs(bz) - 2) < 0.1) {
     if (bz > 0) {
       // magnetic field +0.2 T
-      mRPhiC2 = 4.55e-2f;
+      mRPhiC2 = 0.098f;
     } else {
       // magnetic field -0.2 T
-      mRPhiC2 = 4.55e-2f;
+      mRPhiC2 = 0.098f;
     }
   } else if (CAMath::Abs(CAMath::Abs(bz) - 5) < 0.1) {
     if (bz > 0) {
       // magnetic field +0.5 T
-      mRPhiC2 = 0.0961f;
+      mRPhiC2 = 0.058f;
     } else {
       // magnetic field -0.5 T
-      mRPhiC2 = 0.1156f;
+      mRPhiC2 = 0.072f;
     }
   } else {
     LOGP(warning, "No error parameterization available for Bz= {}. Keeping default value (sigma_y = const. = 1cm)", bz);

GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h

@@ -70,14 +70,14 @@ class GPUTRDRecoParam
   float mRPhiATgp{1.f}; ///< parameterization for tracklet position resolution
   float mRPhiC2{0.f};   ///< parameterization for tracklet position resolution
   // angle
-  float mDyA2{1.225e-3f}; ///< parameterization for tracklet angular resolution
-  float mDyC2{0.f};       ///< parameterization for tracklet angular resolution
+  float mDyA2{6e-3f}; ///< parameterization for tracklet angular resolution
+  float mDyC2{0.3f};       ///< parameterization for tracklet angular resolution
   // variation in y when dy variates by one sigma (= cov / sigma_dy = corr * sigma_y) for tracklet (valid within ~[-0.6,0.6] of dy_trklt-dy_trk)
-  float mCorrYDyA{7.8e-4f};
-  float mCorrYDyC{2.2e-2f};
+  float mCorrYDyA{9e-4f};
+  float mCorrYDyC{1.5e-2f};
   // error parametrization vs angular pull (pol2)
-  float mPullA{6.8e-3f};
-  float mPullB{0.049f};
+  float mPullA{2.8e-3f};
+  float mPullB{0.040f};
   // error parametrization of y position vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy))
   float mOccA{3.3e-4f};
   // error parametrization for dy vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy))

GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx

@@ -496,8 +496,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
   const int32_t nMaxChambersToSearch = 4;
 
   mDebug->SetGeneralInfo(mNEvents, mNTracks, iTrk, t->getPt());
-  
-//LOGF(info, "start id: %d errTrkY: %f", (int)trkWork->getRefGlobalTrackIdRaw(), CAMath::Sqrt(trkWork->getSigmaY2()));
+
   for (int32_t iLayer = 0; iLayer < kNLayers; ++iLayer) {
     nCurrHypothesis = 0;
     bool isOK = false; // if at least one candidate could be propagated or the track was stopped this becomes true
@@ -535,15 +534,15 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
         }
         continue;
       }
-//LOGF(info, "prop id: %d, layer: %d errTrkY: %f", (int)trkWork->getRefGlobalTrackIdRaw(), iLayer, CAMath::Sqrt(trkWork->getSigmaY2()));
+
       // rotate track in new sector in case of sector crossing
       if (!AdjustSector(prop, trkWork)) {
         if (ENABLE_INFO) {
           GPUInfo("Adjusting sector failed for track %i candidate %i in layer %i", iTrk, iCandidate, iLayer);
         }
         continue;
       }
-//LOGF(info, "adj id: %d, layer: %d errTrkY: %f", (int)trkWork->getRefGlobalTrackIdRaw(), iLayer, CAMath::Sqrt(trkWork->getSigmaY2()));
+
       // check if track is findable
       if (IsGeoFindable(trkWork, iLayer, prop->getAlpha(), zShiftTrk)) {
         trkWork->setIsFindable(iLayer);
@@ -594,7 +593,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
           }
         }
 
-        //LOGF(info, "id: %d, layer: %d errTrkY: %f", (int)trkWork->getRefGlobalTrackIdRaw(), iLayer, CAMath::Sqrt(trkWork->getSigmaY2()));
         // first propagate track to x of tracklet
         for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) {
           if (CAMath::Abs(trkWork->getY() - spacePoints[trkltIdx].getY()) > roadY || CAMath::Abs(trkWork->getZ() + zShiftTrk - spacePoints[trkltIdx].getZ()) > roadZ) {
@@ -658,13 +656,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
           // Correction of y position based on angular pull
           if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
             float corrPull = - angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
-            // in the tails with very large angle difference, the correlation becomes flat
-            if (spacePoints[trkltIdx].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp()) > 0.6) {
-              corrPull = - 0.6 / mRecoParam->getDyRes(trkWork->getSnp(), nTrackletsChamber) * mRecoParam->getCorrYDy(trkWork->getSnp());
-            }
-            if (spacePoints[trkltIdx].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp()) < -0.6) {
-              corrPull = 0.6 / mRecoParam->getDyRes(trkWork->getSnp(), nTrackletsChamber) * mRecoParam->getCorrYDy(trkWork->getSnp());
-            }
             yPosCorr += corrPull;
           }
 
@@ -678,8 +669,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
             RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), (Param().rec.trd.useAngularPull == 2 || Param().rec.trd.useAngularPull == 4 ? angularPull : 0.f), nTrackletsChamber, trkltCovTmp);
             trkltCovTmp[0] += yAddErrPileUp2;
             float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp);
-            //if (currDet/30 != sector0 && chi2 < 12) LOGF(info, "track matching : id: %d layer %d trkPt: %f trkltY: %f trkY: %f %f errYtrklt: %f errYtrk: %f diff: %f chi2: %f alpha: %f X: %f", (int)trkWork->getRefGlobalTrackIdRaw(), iLayer, 1./trkWork->getQ2Pt(), yPosCorr, trkWork->getY(), projY, CAMath::Sqrt(trkltCovTmp[0]), CAMath::Sqrt(trkWork->getSigmaY2()), yPosCorr - projY, chi2, trkWork->getAlpha(), trkWork->getX());
-                      
+
             if (Param().rec.trd.addDeflectionInChi2 >= 1 && (trkWork->getSnp() < 1.f - 1e-6f) && (trkWork->getSnp() > -1.f + 1e-6f)) {
               // we add the slope in the chi2 calculation
               float trkltCovTmpWithDy[6] = {trkltCovTmp[0], trkltCovTmp[1], trkltCovTmp[2], 0.f, 0.f, 0.f};
@@ -709,15 +699,13 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
                   // In this case we take into account the full likelihood, so we replace (deltaDy/sigmaDy)^2 by -2*ln(likelihood), which is the same in the default Gaussian case
                   double likelihood = mRecoParam->getDyLikelihood(trkWork->getSnp(), spacePoints[trkltIdx].getDy() + dyTiltCorr, nTrackletsChamber);
                   if (likelihood < 1e-6f) continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation 
-                  deltaDy = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaDy2);
+                  deltaDy = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaDy2) * (deltaDy > 0.f ? 1.f : -1.f);
                 }
                 if (Param().rec.trd.addDeflectionInChi2 == 3) {
                   // We do the same for deltaZ
-                  //LOGF(info, "before z new: %f %f", deltaZ, deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy);
                   double likelihood = mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2()));
                   if (likelihood < 1e-6f) continue;
-                  deltaZ = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaZ2);
-                  //LOGF(info, "after z new: %f %f likelihood: %f %f", deltaZ, deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy, mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2())), -2.f * CAMath::Log(mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2()))));
+                  deltaZ = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaZ2) * (deltaZ > 0.f ? 1.f : -1.f);
                 }
                 chi2 = deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy;
               }
@@ -855,13 +843,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
       // Correction of y position based on angular pull
       if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
         float corrPull = - angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
-        // in the tails with very large angle difference, the correlation becomes flat
-        if (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp()) > 0.6) {
-          corrPull = - 0.6 / mRecoParam->getDyRes(trkWork->getSnp(), nTrackletsChamber) * mRecoParam->getCorrYDy(trkWork->getSnp());
-        }
-        if (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp()) < -0.6) {
-          corrPull = 0.6 / mRecoParam->getDyRes(trkWork->getSnp(), nTrackletsChamber) * mRecoParam->getCorrYDy(trkWork->getSnp());
-        }
         yPosCorrUp += corrPull;
       }