TRD: update of T0 value and possibility to disable shift in the tracklet transformer

Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx

@@ -186,7 +186,7 @@ void BarrelAlignmentSpec::init(InitContext& ic)
         mTRDTransformer->setApplyXOR();
       }
       auto prevShift = mTRDTransformer->isShiftApplied();
-      if (getenv("ALIEN_JDL_LPMPRODUCTIONTYPE") && std::strcmp(getenv("ALIEN_JDL_LPMPRODUCTIONTYPE"), "MC") == 0) {
+      if ((getenv("ALIEN_JDL_LPMPRODUCTIONTYPE") && std::strcmp(getenv("ALIEN_JDL_LPMPRODUCTIONTYPE"), "MC") == 0) || !ic.options().get<bool>("apply-shift")) {
         // apply artificial pad shift in case non-ideal alignment is used to compensate for shift in current alignment from real data
         mTRDTransformer->setApplyShift(false);
       }
@@ -365,6 +365,7 @@ DataProcessorSpec getBarrelAlignmentSpec(GTrackID::mask_t srcMP, GTrackID::mask_
   bool loadTPCCalib = false;
   Options opts{
     ConfigParamSpec{"apply-xor", o2::framework::VariantType::Bool, false, {"flip the 8-th bit of slope and position (for processing TRD CTFs from 2021 pilot beam)"}},
+    ConfigParamSpec{"apply-shift", o2::framework::VariantType::Bool, true, {"whether we apply the one pad shift for the TRD in data (never applied in MC)"}},
     ConfigParamSpec{"allow-afterburner-tracks", VariantType::Bool, false, {"allow using ITS-TPC afterburner tracks"}},
     ConfigParamSpec{"ignore-ccdb-alignment", VariantType::Bool, false, {"do not aplly CCDB alignment to ideal geometry"}},
     ConfigParamSpec{"initial-params-file", VariantType::String, "", {"initial parameters file"}},

Detectors/TRD/base/src/TrackletTransformer.cxx

@@ -61,8 +61,8 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane*
 
 float TrackletTransformer::calibrateX(double x) const
 {
-  // Hard-coded value will be replaced once t0 calibration is available from CCDB
-  float t0Correction = -0.279;
+  // Hard-coded value since t0 is rather stable as a function of time
+  float t0Correction = -0.230;
   return x += t0Correction;
 }
 
@@ -110,7 +110,6 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet, b
 
   float calibratedX = calibrateX(x);
 
-  // NOTE: Correction to y position based on x calibration NOT YET implemented. Need t0.
   if (trackingFrame) {
     std::array<float, 3> sectorSpacePoint = transformL2T(detector, std::array<double, 3>{calibratedX, y, z});
     LOG(debug) << "x: " << sectorSpacePoint[0] << " | "

Detectors/TRD/calibration/src/TrackBasedCalib.cxx

@@ -165,10 +165,10 @@ int TrackBasedCalib::filldEdx(gsl::span<const TrackTRD>& tracks, bool isTPCTRD)
 
       float localGainCorr = 1.;
       if (mLocalGain) {
-        localGainCorr = mLocalGain->getValue(trkltDet, tracklet.getPadCol(mApplyShift), tracklet.getPadRow());
+        localGainCorr = mLocalGain->getValue(trkltDet, tracklet.getPadCol(false), tracklet.getPadRow());
       }
       if (TMath::Abs(localGainCorr) < 0.0001f) {
-        LOGP(warn, "Invalid localGainCorr {} for det {}, pad col {}, pad row {}", localGainCorr, trkltDet, tracklet.getPadCol(mApplyShift), tracklet.getPadRow());
+        LOGP(warn, "Invalid localGainCorr {} for det {}, pad col {}, pad row {}", localGainCorr, trkltDet, tracklet.getPadCol(false), tracklet.getPadRow());
         continue;
       }
 

Detectors/TRD/qc/src/Tracking.cxx

@@ -85,6 +85,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
 
   // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets
   float tCorrPileUp = 0.;
+  float tCorrPileUp_weightedMean = 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
@@ -117,13 +118,15 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
     sumCorr += probBC * deltaBC;
     sumCorr2 += probBC * deltaBC * deltaBC;
     sumProb += probBC;
+    tCorrPileUp_weightedMean -= probBC * deltaBC;
     if (probBC > maxProb) {
       maxProb = probBC;
       tCorrPileUp = -deltaBC;
     }
   }
   if (sumProb > 1e-6) {
     tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp;
+    tCorrPileUp_weightedMean /= sumProb;
   }
 
   for (int iLayer = 0; iLayer < NLAYER; ++iLayer) {
@@ -158,7 +161,8 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
 
     // 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 yCorrPileUp = tCorrPileUp * slopeFactor;
+    float yCorrPileUp = tCorrPileUp_weightedMean * slopeFactor;
     float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
     float yPosCorrUp = mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp;
 
@@ -184,7 +188,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
     // \sqrt{ (dx/dx)^2 + (dy/dx)^2 + (dz/dx)^2}
     auto tphi = trk.getSnp() / std::sqrt((1.f - trk.getSnp()) * (1.f + trk.getSnp()));
     auto trackletLength = std::sqrt(1.f + tphi * tphi + trk.getTgl() * trk.getTgl());
-    auto cor = mLocalGain.getValue(tracklet.getHCID() / 2, tracklet.getPadCol(mApplyShift), tracklet.getPadRow()) * trackletLength;
+    auto cor = mLocalGain.getValue(tracklet.getHCID() / 2, tracklet.getPadCol(false), tracklet.getPadRow()) * trackletLength;
     float q0{tracklet.getQ0() / cor}, q1{tracklet.getQ1() / cor}, q2{tracklet.getQ2() / cor};
 
     // z-row merging
@@ -193,9 +197,9 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
         if (tracklet.getTrackletWord() == trklt.getTrackletWord()) { // skip original tracklet
           continue;
         }
-        if (std::abs(tracklet.getPadCol(mApplyShift) - trklt.getPadCol(mApplyShift)) <= 1 && std::abs(tracklet.getPadRow() - trklt.getPadRow()) == 1) {
+        if (std::abs(tracklet.getPadCol(false) - trklt.getPadCol(false)) <= 1 && std::abs(tracklet.getPadRow() - trklt.getPadRow()) == 1) {
           // Add charge information
-          auto cor = mLocalGain.getValue(trklt.getHCID() / 2, trklt.getPadCol(mApplyShift), trklt.getPadRow()) * trackletLength;
+          auto cor = mLocalGain.getValue(trklt.getHCID() / 2, trklt.getPadCol(false), trklt.getPadRow()) * trackletLength;
           q0 += trklt.getQ0() / cor;
           q1 += trklt.getQ1() / cor;
           q2 += trklt.getQ2() / cor;

Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h

@@ -38,10 +38,10 @@ struct TRDSimParams : public o2::conf::ConfigurableParamHelper<TRDSimParams> {
   bool doTR = true;                                       ///< switch for transition radiation
   SimParam::GasMixture gas = SimParam::GasMixture::Xenon; ///< the gas mixture in the TRD
   // TRF parameters
-  int trf = 0;                  // Sampled TRF function. 0: default TRF, 1: TRF described in TRF TDR, 2: No TRF, 3: Landau dist as TRF (parameters specified below)
+  int trf = 1;                  // Sampled TRF function. 0: default TRF, 1: TRF described in TRF TDR, 2: No TRF, 3: Landau dist as TRF (parameters specified below)
   float trf_landau_mu = 0.;     // Mu of the Landau distribution used to describe TRF
   float trf_landau_sigma = .03; // Sigma of the Landau distribution used to describe TRF
-  float t0 = -0.130;            // Time delay before start of amplification peak, in us
+  float t0 = -0.170;            // Time delay before start of amplification peak, in us
   O2ParamDef(TRDSimParams, "TRDSimParams");
 };
 

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

@@ -152,8 +152,7 @@ DataProcessorSpec getTRDGlobalTrackingQCSpec(o2::dataformats::GlobalTrackID::mas
     "trd-tracking-qc",
     dataRequest->inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<TRDGlobalTrackingQC>(dataRequest, ggRequest, isTPCavailable, src)},
-    Options{}};
+    AlgorithmSpec{adaptFromTask<TRDGlobalTrackingQC>(dataRequest, ggRequest, isTPCavailable, src)}};
 }
 
 template <typename T>

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

@@ -82,6 +82,7 @@ class TRDGlobalTracking : public o2::framework::Task
   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
+  int32_t mTCorrPileUp_weightedMean{0};                          ///< average 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;

Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx

@@ -542,6 +542,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
     float sumCorr = 0.f;
     float sumCorr2 = 0.f;
     float sumProb = 0.f;
+    mTCorrPileUp_weightedMean = 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()) {
@@ -565,13 +566,15 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
       sumCorr += probBC * deltaBC;
       sumCorr2 += probBC * deltaBC * deltaBC;
       sumProb += probBC;
+      mTCorrPileUp_weightedMean -= deltaBC * probBC;
       if (probBC > maxProb) {
         maxProb = probBC;
         mTCorrPileUp = -deltaBC;
       }
     }
     if (sumProb > 1e-6) {
       mTErrPileUp2 = sumCorr2 / sumProb - 2 * mTCorrPileUp * sumCorr / sumProb + mTCorrPileUp * mTCorrPileUp;
+      mTCorrPileUp_weightedMean /= sumProb;
     }
 
     nTrackletsAttached += trdTrack.getNtracklets();
@@ -926,7 +929,8 @@ 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 = mTCorrPileUp * slopeFactor;
+    // float yCorrPileUp = mTCorrPileUp * slopeFactor;
+    float yCorrPileUp = mTCorrPileUp_weightedMean * slopeFactor;
     float yAddErrPileUp2 = mTErrPileUp2 * slopeFactor * slopeFactor;
     float yPosCorrUp = mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp;
 

Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx

@@ -35,7 +35,8 @@ void TRDTrackletTransformerSpec::init(o2::framework::InitContext& ic)
   if (ic.options().get<bool>("apply-xor")) {
     mTransformer.setApplyXOR();
   }
-  if (getenv("ALIEN_JDL_LPMPRODUCTIONTYPE") && std::strcmp(getenv("ALIEN_JDL_LPMPRODUCTIONTYPE"), "MC") == 0) {
+  if ((getenv("ALIEN_JDL_LPMPRODUCTIONTYPE") && std::strcmp(getenv("ALIEN_JDL_LPMPRODUCTIONTYPE"), "MC") == 0) || ic.options().get<bool>("disable-shift")) {
+    LOGF(info, "TRD pad shift is not applied");
     // apply artificial pad shift in case non-ideal alignment is used to compensate for shift in current alignment from real data
     mTransformer.setApplyShift(false);
   }
@@ -174,7 +175,8 @@ o2::framework::DataProcessorSpec getTRDTrackletTransformerSpec(bool trigRecFilte
     outputs,
     AlgorithmSpec{adaptFromTask<TRDTrackletTransformerSpec>(dataRequest, ggRequest, trigRecFilterActive)},
     Options{
-      {"apply-xor", o2::framework::VariantType::Bool, false, {"flip the 8-th bit of slope and position (for processing CTFs from 2021 pilot beam)"}}}};
+      {"apply-xor", o2::framework::VariantType::Bool, false, {"flip the 8-th bit of slope and position (for processing CTFs from 2021 pilot beam)"}},
+      {"disable-shift", o2::framework::VariantType::Bool, false, {"whether we apply the one pad shift in data (always disabled in MC)"}}}};
 }
 
 } //end namespace trd

GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx

@@ -615,6 +615,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
           // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation.
           // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers
           float yCorrPileUp = 0.f;
+          float yCorrPileUp_weightedMean = 0.f;
           float yAddErrPileUp2 = 0.f;
           // additional zshiftTrk for TPC only tracks due to pile-up
           float zShiftTrkPileUp = 0.f;
@@ -632,6 +633,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
               sumCorr += probBC * slopeFactor * deltaBC;
               sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC;
               sumProb += probBC;
+              yCorrPileUp_weightedMean -= slopeFactor * deltaBC * probBC;
               if (probBC > maxProb) {
                 maxProb = probBC;
                 yCorrPileUp = -slopeFactor * deltaBC;
@@ -640,6 +642,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
             }
             if (sumProb > 1e-6f) {
               yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
+              yCorrPileUp_weightedMean /= sumProb;
             }
           }
 
@@ -649,7 +652,8 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
 
           // correction for mean z position of tracklet (is not the center of the pad if track eta != 0)
           float zPosCorr = spacePoints[trkltIdx].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl();
-          float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp;
+          // float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp;
+          float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp_weightedMean;
           zPosCorr -= zShiftTrk + zShiftTrkPileUp; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision
 
           // Correction of y position based on angular pull
@@ -808,6 +812,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
       // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation.
       // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers
       float yCorrPileUp = 0.f;
+      float yCorrPileUp_weightedMean = 0.f;
       float yAddErrPileUp2 = 0.f;
       float zShiftTrkPileUp = 0.f;
       if (nIdxBCMax - nIdxBCMin >= 2) {
@@ -824,6 +829,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
           sumCorr += probBC * slopeFactor * deltaBC;
           sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC;
           sumProb += probBC;
+          yCorrPileUp_weightedMean -= slopeFactor * deltaBC * probBC;
           if (probBC > maxProb) {
             maxProb = probBC;
             yCorrPileUp = -slopeFactor * deltaBC;
@@ -832,11 +838,13 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
         }
         if (sumProb > 1e-6f) {
           yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
+          yCorrPileUp_weightedMean /= sumProb;
         }
       }
 
       zPosCorrUp -= zShiftTrkPileUp;
-      yPosCorrUp += yCorrPileUp;
+      // yPosCorrUp += yCorrPileUp;
+      yPosCorrUp += yCorrPileUp_weightedMean;
 
       const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector();
       int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet];