TRD: using slope to correct y position and reject more fakes (#15799)

* add possibility to account for exponential tail of tracklet slope in chi2

* TRD: using slope to correct y position and reject more fakes

* some fixes and parameter updates

* clang format

* clang format

* small fix

* test to find origin of the build errors

* clang format

* Please consider the following formatting changes

---------

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: glegras

Detectors/TRD/base/include/TRDBase/PadPlane.h

@@ -120,7 +120,41 @@ class PadPlane
 
   GPUd() int getPadRowNumberROC(double z) const;
   GPUd() double getPadRow(double z) const;
-  GPUd() int getPadColNumber(double rphi) const;
+
+  GPUd() int getPadColNumber(double rphi) const
+  {
+    //
+    // Finds the pad column number for a given rphi-position
+    //
+
+    int col = 0;
+    int nabove = 0;
+    int nbelow = 0;
+    int middle = 0;
+
+    if ((rphi < getCol0()) || (rphi > getColEnd())) {
+      col = -1;
+
+    } else {
+      nabove = mNcols;
+      nbelow = 0;
+      while (nabove - nbelow > 1) {
+        middle = (nabove + nbelow) / 2;
+        if (rphi == mPadCol[middle]) {
+          col = middle;
+        }
+        if (rphi > mPadCol[middle]) {
+          nbelow = middle;
+        } else {
+          nabove = middle;
+        }
+      }
+      col = nbelow;
+    }
+
+    return col;
+  }
+
   GPUd() double getPad(double y, double z) const;
 
   GPUd() double getTiltOffset(int row, double rowOffset) const

Detectors/TRD/base/src/PadPlane.cxx

@@ -76,41 +76,6 @@ int PadPlane::getPadRowNumberROC(double z) const
   return row;
 }
 
-//_____________________________________________________________________________
-int PadPlane::getPadColNumber(double rphi) const
-{
-  //
-  // Finds the pad column number for a given rphi-position
-  //
-
-  int col = 0;
-  int nabove = 0;
-  int nbelow = 0;
-  int middle = 0;
-
-  if ((rphi < getCol0()) || (rphi > getColEnd())) {
-    col = -1;
-
-  } else {
-    nabove = mNcols;
-    nbelow = 0;
-    while (nabove - nbelow > 1) {
-      middle = (nabove + nbelow) / 2;
-      if (rphi == mPadCol[middle]) {
-        col = middle;
-      }
-      if (rphi > mPadCol[middle]) {
-        nbelow = middle;
-      } else {
-        nabove = middle;
-      }
-    }
-    col = nbelow;
-  }
-
-  return col;
-}
-
 void PadPlane::setNcols(int n)
 {
   if (n > MAXCOLS) {

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;
 
@@ -206,6 +213,12 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       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;
     for (int iLayer = NLAYER - 1; iLayer >= 0; --iLayer) {
@@ -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

@@ -160,10 +160,14 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
     float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f;
     float yCorrPileUp = tCorrPileUp * slopeFactor;
     float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
+    float yPosCorrUp = mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp;
 
     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());
+    yPosCorrUp += corrPull;
 
-    std::array<float, 2> trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp};
+    std::array<float, 2> trkltPosUp{yPosCorrUp, zPosCorrUp};
     std::array<float, 3> trkltCovUp;
     mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp, angularPull, 0);
     trkltCovUp[0] += yAddErrPileUp2;

Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h

@@ -78,6 +78,8 @@ class TRDGlobalTracking : public o2::framework::Task
   float mTPCTBinMUSInv{1.f / mTPCTBinMUS};                       ///< inverse width of a TPC time bin in 1/us
   float mTPCVdrift{2.58f};                                       ///< TPC drift velocity (for shifting TPC tracks along Z)
   float mTPCTDriftOffset{0.f};                                   ///< TPC drift time additive offset
+  int32_t mTCorrPileUp{0};                                       ///< most probable correction in number of BCs due to pile-up in TRD
+  float mTErrPileUp2{0.f};                                       ///< error on correction in number of BCs due to pile-up in TRD
   std::shared_ptr<o2::globaltracking::DataRequest> mDataRequest; ///< seeding input (TPC-only, ITS-TPC or both)
   std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
   o2::tpc::VDriftHelper mTPCVDriftHelper{};

Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx

@@ -153,6 +153,9 @@ void TRDGlobalTracking::updateTimeDependentParams(ProcessingContext& pc)
       mBase->init(pc);
       mBase->setLocalGainFactors(pc.inputs().get<o2::trd::LocalGainFactor*>("localgainfactors").get());
     }
+
+    pc.inputs().get<std::array<int, constants::MAXCHAMBER>*>("chamberstatus"); // called to trigger finaliseCCDB
+    // pc.inputs().get<o2::trd::PadStatus*>("padstatus");                         // called to trigger finaliseCCDB
   }
 
   const auto& trackTune = TrackTuneParams::Instance();
@@ -198,6 +201,34 @@ void TRDGlobalTracking::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
     return;
   }
 #endif
+  if (matcher == ConcreteDataMatcher("TRD", "CHAMBERSTATUS", 0)) {
+    LOG(info) << "chamber status object updated";
+    const std::array<int, constants::MAXCHAMBER>* chamberStatus = (const std::array<int, constants::MAXCHAMBER>*)obj;
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
+      if ((*chamberStatus)[iDet] == 3) {
+        mTracker->SetChamberStatus(iDet, false); // chamber is good
+      } else {
+        mTracker->SetChamberStatus(iDet, true); // chamber is bad
+      }
+    }
+    return;
+  }
+  /*if (matcher == ConcreteDataMatcher("TRD", "PADSTATUS", 0)) {
+    LOG(info) << "pad status object updated";
+    const o2::trd::PadStatus* padStatus = (const o2::trd::PadStatus*)obj;
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
+      for (int iCol = 0; iCol < constants::NCOLUMN; iCol++) {
+        for (int iRow = 0; iRow < ((iDet % 30) / 6 == 2 ? constants::NROWC0 : constants::NROWC1); iRow++) {
+          if (padStatus->isMasked(iDet, iCol, iRow) || padStatus->isNotConnected(iDet, iCol, iRow)) {
+            mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, true); // pad is masked
+          } else {
+            mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, false); // pad is not masked
+          }
+        }
+      }
+    }
+    return;
+  }*/
 }
 
 void TRDGlobalTracking::fillMCTruthInfo(const TrackTRD& trk, o2::MCCompLabel lblSeed, std::vector<o2::MCCompLabel>& lblContainerTrd, std::vector<o2::MCCompLabel>& lblContainerMatch, const o2::dataformats::MCTruthContainer<o2::MCCompLabel>* trkltLabels) const
@@ -500,6 +531,45 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
     if (trdTrack.getChi2() / trdTrack.getNtracklets() > mTracker->Param().rec.trd.maxChi2Red) {
       continue;
     }
+
+    // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets
+    float maxProb = 0.f;
+    // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
+    float sumCorr = 0.f;
+    float sumCorr2 = 0.f;
+    float sumProb = 0.f;
+    for (int iBC = 0; iBC < mTriggeredBCFT0.size(); iBC++) {
+      int deltaBC = roundf(mTriggeredBCFT0[iBC] - mChainTracking->mIOPtrs.trdTriggerTimes[trdTrack.getCollisionId()] / o2::constants::lhc::LHCBunchSpacingMUS);
+      if (deltaBC <= mRecoParam.getPileUpRangeBefore() || deltaBC >= mRecoParam.getPileUpRangeAfter()) {
+        continue;
+      }
+      // collect the charges
+      std::array<int, 6> q0;
+      std::array<int, 6> q1;
+      for (int iLy = 0; iLy < NLAYER; iLy++) {
+        int trkltId = trdTrack.getTrackletIndex(iLy);
+        if (trkltId < 0) {
+          q0[iLy] = -1;
+          q1[iLy] = -1;
+        } else {
+          q0[iLy] = mTrackletsRaw[trkltId].getQ0();
+          q1[iLy] = mTrackletsRaw[trkltId].getQ1();
+        }
+      }
+      // get pile-up probability
+      float probBC = mRecoParam.getPileUpProbTrack(deltaBC, q0, q1);
+      sumCorr += probBC * deltaBC;
+      sumCorr2 += probBC * deltaBC * deltaBC;
+      sumProb += probBC;
+      if (probBC > maxProb) {
+        maxProb = probBC;
+        mTCorrPileUp = -deltaBC;
+      }
+    }
+    if (sumProb > 1e-6) {
+      mTErrPileUp2 = sumCorr2 / sumProb - 2 * mTCorrPileUp * sumCorr / sumProb + mTCorrPileUp * mTCorrPileUp;
+    }
+
     nTrackletsAttached += trdTrack.getNtracklets();
     auto trackGID = trdTrack.getRefGlobalTrackId();
     if (trackGID.includesDet(GTrackID::Source::ITS)) {
@@ -529,7 +599,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
       } else {
         tracksOutTPC.back().setPileUpDistance(mTracker->Param().rec.trd.pileupBwdNBC, mTracker->Param().rec.trd.pileupFwdNBC);
       }
-      if (!refitTPCTRDTrack(tracksOutTPC.back(), mChainTracking->mIOPtrs.trdTriggerTimes[trdTrack.getCollisionId()], &inputTracks) || std::isnan(tracksOutTPC.back().getSnp())) {
+      if (!refitTPCTRDTrack(tracksOutTPC.back(), mChainTracking->mIOPtrs.trdTriggerTimes[trdTrack.getCollisionId()] - mTCorrPileUp * o2::constants::lhc::LHCBunchSpacingMUS, &inputTracks) || std::isnan(tracksOutTPC.back().getSnp())) {
         tracksOutTPC.pop_back();
         ++nTracksFailedTPCTRDRefit;
         continue;
@@ -734,7 +804,8 @@ bool TRDGlobalTracking::refitTPCTRDTrack(TrackTRD& trk, float timeTRD, o2::globa
     return false;
   }
   if (pileUpOn) { // account pileup time uncertainty in Z errors
-    timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS();
+    // timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS();
+    timeZErr = mTPCVdrift * mTPCVdrift * mTErrPileUp2;
     outerParam.updateCov(timeZErr, o2::track::CovLabels::kSigZ2);
   }
   if (!refitTRDTrack(trk, chi2Out, false, true)) {
@@ -818,46 +889,6 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
     }
   }
 
-  // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets
-  float tCorrPileUp = 0.;
-  float tErrPileUp2 = 0;
-  float maxProb = 0.f;
-  // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
-  float sumCorr = 0.f;
-  float sumCorr2 = 0.f;
-  float sumProb = 0.f;
-  for (int iBC = 0; iBC < mTriggeredBCFT0.size(); iBC++) {
-    int deltaBC = roundf(mTriggeredBCFT0[iBC] - mChainTracking->mIOPtrs.trdTriggerTimes[trk.getCollisionId()] / o2::constants::lhc::LHCBunchSpacingMUS);
-    if (deltaBC <= mRecoParam.getPileUpRangeBefore() || deltaBC >= mRecoParam.getPileUpRangeAfter()) {
-      continue;
-    }
-    // collect the charges
-    std::array<int, 6> q0;
-    std::array<int, 6> q1;
-    for (int iLy = 0; iLy < NLAYER; iLy++) {
-      int trkltId = trk.getTrackletIndex(iLy);
-      if (trkltId < 0) {
-        q0[iLy] = -1;
-        q1[iLy] = -1;
-      } else {
-        q0[iLy] = mTrackletsRaw[trkltId].getQ0();
-        q1[iLy] = mTrackletsRaw[trkltId].getQ1();
-      }
-    }
-    // get pile-up probability
-    float probBC = mRecoParam.getPileUpProbTrack(deltaBC, q0, q1);
-    sumCorr += probBC * deltaBC;
-    sumCorr2 += probBC * deltaBC * deltaBC;
-    sumProb += probBC;
-    if (probBC > maxProb) {
-      maxProb = probBC;
-      tCorrPileUp = -deltaBC;
-    }
-  }
-  if (sumProb > 1e-6) {
-    tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp;
-  }
-
   if (inwards) {
     // reset covariance to something big for inwards refit
     trkParam->resetCovariance(100);
@@ -891,13 +922,20 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
 
     // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin
     float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f;
-    float yCorrPileUp = tCorrPileUp * slopeFactor;
-    float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
+    float yCorrPileUp = mTCorrPileUp * slopeFactor;
+    float yAddErrPileUp2 = mTErrPileUp2 * slopeFactor * slopeFactor;
+    float yPosCorrUp = mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp;
 
     int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet);
     float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber));
 
-    std::array<float, 2> trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp};
+    // 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());
+      yPosCorrUp += corrPull;
+    }
+
+    std::array<float, 2> trkltPosUp{yPosCorrUp, zPosCorrUp};
     std::array<float, 3> trkltCovUp;
     mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp, (mRec->GetParam().rec.trd.useAngularPull != 0 ? angularPull : 0.), nTrackletsChamber);
     trkltCovUp[0] += yAddErrPileUp2;
@@ -1007,6 +1045,10 @@ DataProcessorSpec getTRDGlobalTrackingSpec(bool useMC, GTrackID::mask_t src, boo
     inputs.emplace_back("localgainfactors", "TRD", "LOCALGAINFACTORS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/LocalGainFactor"));
   }
 
+  // request list of bad chambers and masked pads to estimate better the number of findable tracklets
+  inputs.emplace_back("chamberstatus", "TRD", "CHAMBERSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/DCSDPsFedChamberStatus"));
+  // inputs.emplace_back("padstatus", "TRD", "PADSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/PadStatus"));
+
   if (GTrackID::includesSource(GTrackID::Source::ITSTPC, src)) {
     outputs.emplace_back(o2::header::gDataOriginTRD, "MATCH_ITSTPC", 0, Lifetime::Timeframe);
     outputs.emplace_back(o2::header::gDataOriginTRD, "TRGREC_ITSTPC", 0, Lifetime::Timeframe);

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);