Fill alternative TPC dedx according to dEdxClusterRejectionFlagMaskAlt

DataFormats/Detectors/TPC/include/DataFormatsTPC/TrackTPC.h

@@ -128,6 +128,9 @@ class TrackTPC : public o2::track::TrackParCov
   GPUd() const dEdxInfo& getdEdx() const { return mdEdx; }
   GPUd() void setdEdx(const dEdxInfo& v) { mdEdx = v; }
 
+  GPUd() const dEdxInfo& getdEdxAlt() const { return mdEdxAlt; }
+  GPUd() void setdEdxAlt(const dEdxInfo& v) { mdEdxAlt = v; }
+
  private:
   float mTime0 = 0.f;                 ///< Assumed time of the vertex that created the track in TPC time bins, 0 for triggered data
   float mDeltaTFwd = 0;               ///< max possible increment to mTime0
@@ -136,9 +139,10 @@ class TrackTPC : public o2::track::TrackParCov
   float mChi2 = 0.f;                  // Chi2 of the track
   o2::track::TrackParCov mOuterParam; // Track parameters at outer end of TPC.
   dEdxInfo mdEdx;                     // dEdx Information
+  dEdxInfo mdEdxAlt;                  // dEdx alternative Information
   ClusRef mClustersReference;         // reference to externale cluster indices
 
-  ClassDefNV(TrackTPC, 4);
+  ClassDefNV(TrackTPC, 5);
 };
 
 } // namespace tpc

Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h

@@ -416,6 +416,7 @@ class AODProducerWorkflowDPL : public Task
     int16_t tpcdcaR{};
     int16_t tpcdcaZ{};
     uint8_t tpcClusterByteMask{};
+    uint8_t tpcSignalTodEdxNorm{};
     uint8_t tpcdEdxMax0R{};
     uint8_t tpcdEdxMax1R{};
     uint8_t tpcdEdxMax2R{};

Detectors/AOD/src/AODProducerWorkflowSpec.cxx

@@ -365,6 +365,7 @@ void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACurs
     trackQAInfoHolder.tpcdcaR,
     trackQAInfoHolder.tpcdcaZ,
     trackQAInfoHolder.tpcClusterByteMask,
+    trackQAInfoHolder.tpcSignalTodEdxNorm,
     trackQAInfoHolder.tpcdEdxMax0R,
     trackQAInfoHolder.tpcdEdxMax1R,
     trackQAInfoHolder.tpcdEdxMax2R,
@@ -2542,14 +2543,15 @@ AODProducerWorkflowDPL::TrackExtraInfo AODProducerWorkflowDPL::processBarrelTrac
   if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
     const auto& tpcOrig = data.getTPCTrack(contributorsGID[GIndex::TPC]);
     const auto& tpcClData = mTPCCounters[contributorsGID[GIndex::TPC]];
+    const auto& dEdx = tpcOrig.getdEdx().dEdxTotTPC > 0 ? tpcOrig.getdEdx() : tpcOrig.getdEdxAlt();
     extraInfoHolder.tpcInnerParam = tpcOrig.getP() / tpcOrig.getAbsCharge();
     extraInfoHolder.tpcChi2NCl = tpcOrig.getNClusters() ? tpcOrig.getChi2() / tpcOrig.getNClusters() : 0;
-    extraInfoHolder.tpcSignal = tpcOrig.getdEdx().dEdxTotTPC;
+    extraInfoHolder.tpcSignal = dEdx.dEdxTotTPC;
     extraInfoHolder.tpcNClsFindable = tpcOrig.getNClusters();
     extraInfoHolder.tpcNClsFindableMinusFound = tpcOrig.getNClusters() - tpcClData.found;
     extraInfoHolder.tpcNClsFindableMinusCrossedRows = tpcOrig.getNClusters() - tpcClData.crossed;
     extraInfoHolder.tpcNClsShared = tpcClData.shared;
-    uint32_t clsUsedForPID = tpcOrig.getdEdx().NHitsIROC + tpcOrig.getdEdx().NHitsOROC1 + tpcOrig.getdEdx().NHitsOROC2 + tpcOrig.getdEdx().NHitsOROC3;
+    uint32_t clsUsedForPID = dEdx.NHitsIROC + dEdx.NHitsOROC1 + dEdx.NHitsOROC2 + dEdx.NHitsOROC3;
     extraInfoHolder.tpcNClsFindableMinusPID = tpcOrig.getNClusters() - clsUsedForPID;
     if (src == GIndex::TPC) { // standalone TPC track should set its time from their timebins range
       if (needBCSlice) {
@@ -2625,16 +2627,19 @@ AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int
     }
     trackQAHolder.tpcTime0 = tpcOrig.getTime0();
     trackQAHolder.tpcClusterByteMask = byteMask;
-    const float dEdxNorm = (tpcOrig.getdEdx().dEdxTotTPC > 0) ? 100. / tpcOrig.getdEdx().dEdxTotTPC : 0;
-    trackQAHolder.tpcdEdxMax0R = uint8_t(tpcOrig.getdEdx().dEdxMaxIROC * dEdxNorm);
-    trackQAHolder.tpcdEdxMax1R = uint8_t(tpcOrig.getdEdx().dEdxMaxOROC1 * dEdxNorm);
-    trackQAHolder.tpcdEdxMax2R = uint8_t(tpcOrig.getdEdx().dEdxMaxOROC2 * dEdxNorm);
-    trackQAHolder.tpcdEdxMax3R = uint8_t(tpcOrig.getdEdx().dEdxMaxOROC3 * dEdxNorm);
+    const auto& dEdxInfoAlt = tpcOrig.getdEdxAlt(); // tpcOrig.getdEdx()
+    const float dEdxNorm = (dEdxInfoAlt.dEdxTotTPC > 0) ? 100. / dEdxInfoAlt.dEdxTotTPC : 0;
+    float aodTPCSignal = tpcOrig.getdEdx().dEdxTotTPC > 0 ? tpcOrig.getdEdx().dEdxTotTPC : dEdxInfoAlt.dEdxTotTPC; // what was stored in the AOD
+    trackQAHolder.tpcSignalTodEdxNorm = uint8_t(aodTPCSignal * dEdxNorm);
+    trackQAHolder.tpcdEdxMax0R = uint8_t(dEdxInfoAlt.dEdxMaxIROC * dEdxNorm);
+    trackQAHolder.tpcdEdxMax1R = uint8_t(dEdxInfoAlt.dEdxMaxOROC1 * dEdxNorm);
+    trackQAHolder.tpcdEdxMax2R = uint8_t(dEdxInfoAlt.dEdxMaxOROC2 * dEdxNorm);
+    trackQAHolder.tpcdEdxMax3R = uint8_t(dEdxInfoAlt.dEdxMaxOROC3 * dEdxNorm);
     //
-    trackQAHolder.tpcdEdxTot0R = uint8_t(tpcOrig.getdEdx().dEdxTotIROC * dEdxNorm);
-    trackQAHolder.tpcdEdxTot1R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC1 * dEdxNorm);
-    trackQAHolder.tpcdEdxTot2R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC2 * dEdxNorm);
-    trackQAHolder.tpcdEdxTot3R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC3 * dEdxNorm);
+    trackQAHolder.tpcdEdxTot0R = uint8_t(dEdxInfoAlt.dEdxTotIROC * dEdxNorm);
+    trackQAHolder.tpcdEdxTot1R = uint8_t(dEdxInfoAlt.dEdxTotOROC1 * dEdxNorm);
+    trackQAHolder.tpcdEdxTot2R = uint8_t(dEdxInfoAlt.dEdxTotOROC2 * dEdxNorm);
+    trackQAHolder.tpcdEdxTot3R = uint8_t(dEdxInfoAlt.dEdxTotOROC3 * dEdxNorm);
     ///
     float scaleTOF{0};
     auto contributorsGIDA = data.getSingleDetectorRefs(trackIndex);
@@ -2708,6 +2713,7 @@ AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int
                        << "trackQAHolder.tpcdcaR=" << trackQAHolder.tpcdcaR
                        << "trackQAHolder.tpcdcaZ=" << trackQAHolder.tpcdcaZ
                        << "trackQAHolder.tpcdcaClusterByteMask=" << trackQAHolder.tpcClusterByteMask
+                       << "trackQAHolder.tpcSignalTodEdxNorm=" << trackQAHolder.tpcSignalTodEdxNorm
                        << "trackQAHolder.tpcdEdxMax0R=" << trackQAHolder.tpcdEdxMax0R
                        << "trackQAHolder.tpcdEdxMax1R=" << trackQAHolder.tpcdEdxMax1R
                        << "trackQAHolder.tpcdEdxMax2R=" << trackQAHolder.tpcdEdxMax2R

Framework/Core/include/Framework/AnalysisDataModel.h

@@ -683,7 +683,8 @@ DECLARE_SOA_COLUMN(TPCTime0, tpcTime0, float);                            //! tp
 DECLARE_SOA_COLUMN(TPCDCAR, tpcdcaR, int16_t);                            //! tpc only DCAr
 DECLARE_SOA_COLUMN(TPCDCAZ, tpcdcaZ, int16_t);                            //! tpc only DCAz
 DECLARE_SOA_COLUMN(TPCClusterByteMask, tpcClusterByteMask, uint8_t);      //! tracklet bitmask - track defining 8 tracklets (152=8*19 rows) bit set if nCluster>thr (default 5)
-DECLARE_SOA_COLUMN(TPCdEdxMax0R, tpcdEdxMax0R, uint8_t);                  //! TPC dEdxQMax -ROC0/dEdx
+DECLARE_SOA_COLUMN(TPCSignalTodEdxNorm, tpcSignalTodEdxNorm, uint8_t);    //! 100*TrackExtra.TPCSignal/TrackTPC.mdEdxAlt
+DECLARE_SOA_COLUMN(TPCdEdxMax0R, tpcdEdxMax0R, uint8_t);                  //! TPC dEdxQMax -ROC0/dEdx from TrackTPC.mdEdxAlt
 DECLARE_SOA_COLUMN(TPCdEdxMax1R, tpcdEdxMax1R, uint8_t);                  //! TPC dEdxQMax -ROC1/dEdx
 DECLARE_SOA_COLUMN(TPCdEdxMax2R, tpcdEdxMax2R, uint8_t);                  //! TPC dEdxQMax -ROC2/dEdx
 DECLARE_SOA_COLUMN(TPCdEdxMax3R, tpcdEdxMax3R, uint8_t);                  //! TPC dEdxQMax -ROC3/dEdx
@@ -736,7 +737,17 @@ DECLARE_SOA_TABLE_VERSIONED(TracksQA_002, "AOD", "TRACKQA", 2, //! trackQA infor
                             trackqa::IsDummy<trackqa::DeltaRefContParamY, trackqa::DeltaRefContParamZ, trackqa::DeltaRefContParamSnp, trackqa::DeltaRefContParamTgl, trackqa::DeltaRefContParamQ2Pt,
                                              trackqa::DeltaRefGloParamY, trackqa::DeltaRefGloParamZ, trackqa::DeltaRefGloParamSnp, trackqa::DeltaRefGloParamTgl, trackqa::DeltaRefGloParamQ2Pt>);
 
-using TracksQAVersion = TracksQA_002;
+DECLARE_SOA_TABLE_VERSIONED(TracksQA_003, "AOD", "TRACKQA", 3, //! trackQA information - version 3 - including alternative dedx normalization
+                            o2::soa::Index<>, trackqa::TrackId, trackqa::TPCTime0, trackqa::TPCDCAR, trackqa::TPCDCAZ, trackqa::TPCClusterByteMask, trackqa::TPCSignalTodEdxNorm,
+                            trackqa::TPCdEdxMax0R, trackqa::TPCdEdxMax1R, trackqa::TPCdEdxMax2R, trackqa::TPCdEdxMax3R,
+                            trackqa::TPCdEdxTot0R, trackqa::TPCdEdxTot1R, trackqa::TPCdEdxTot2R, trackqa::TPCdEdxTot3R,
+                            trackqa::DeltaRefContParamY, trackqa::DeltaRefContParamZ, trackqa::DeltaRefContParamSnp, trackqa::DeltaRefContParamTgl, trackqa::DeltaRefContParamQ2Pt,
+                            trackqa::DeltaRefGloParamY, trackqa::DeltaRefGloParamZ, trackqa::DeltaRefGloParamSnp, trackqa::DeltaRefGloParamTgl, trackqa::DeltaRefGloParamQ2Pt,
+                            trackqa::DeltaTOFdX, trackqa::DeltaTOFdZ,
+                            trackqa::IsDummy<trackqa::DeltaRefContParamY, trackqa::DeltaRefContParamZ, trackqa::DeltaRefContParamSnp, trackqa::DeltaRefContParamTgl, trackqa::DeltaRefContParamQ2Pt,
+                                             trackqa::DeltaRefGloParamY, trackqa::DeltaRefGloParamZ, trackqa::DeltaRefGloParamSnp, trackqa::DeltaRefGloParamTgl, trackqa::DeltaRefGloParamQ2Pt>);
+
+using TracksQAVersion = TracksQA_003;
 using TracksQA = TracksQAVersion::iterator;
 
 namespace fwdtrack

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -158,6 +158,7 @@ AddOptionRTC(enablePID, int8_t, 1, "", 0, "Enable PID response")
 AddOptionRTC(PID_useNsigma, int8_t, 1, "", 0, "Use nSigma instead of absolute distance in PID response")
 AddOptionRTC(adddEdxSubThresholdClusters, int8_t, 1, "", 0, "Add sub threshold clusters in TPC dEdx computation")
 AddOptionRTC(dEdxClusterRejectionFlagMask, int8_t, o2::gpu::GPUTPCGMMergedTrackHit::flagEdge, "", 0, "OR mask of TPC flags that will reject the cluster in dEdx")
+AddOptionRTC(dEdxClusterRejectionFlagMaskAlt, int8_t, o2::gpu::GPUTPCGMMergedTrackHit::flagEdge, "", 0, "OR mask of TPC flags that will reject the cluster in alternative dEdx")
 AddOptionRTC(rejectEdgeClustersInSeeding, int8_t, 0, "", 0, "Reject edge clusters based on uncorrected track Y during seeding")
 AddOptionRTC(rejectEdgeClustersInTrackFit, int8_t, 0, "", 0, "Reject edge clusters based on uncorrected track Y during track fit")
 AddOptionArray(PID_remap, int8_t, 9, (0, 1, 2, 3, 4, 5, 6, 7, 8), "", 0, "Remap Ipid to PID_reamp[Ipid] (no remap if<0)") // BUG: CUDA cannot yet hand AddOptionArrayRTC

GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx

@@ -299,6 +299,7 @@ void* GPUTPCGMMerger::SetPointersOutput(void* mem)
   computePointerWithAlignment(mem, mOutputTracks, mNMaxTracks);
   if (mRec->GetParam().dodEdxDownscaled) {
     computePointerWithAlignment(mem, mOutputTracksdEdx, mNMaxTracks);
+    computePointerWithAlignment(mem, mOutputTracksdEdxAlt, mNMaxTracks);
   }
   computePointerWithAlignment(mem, mClusters, mNMaxOutputTrackClusters);
   if (mRec->GetParam().par.earlyTpcTransform) {

GPU/GPUTracking/Merger/GPUTPCGMMerger.h

@@ -108,6 +108,8 @@ class GPUTPCGMMerger : public GPUProcessor
   GPUhdi() GPUTPCGMMergedTrack* OutputTracks() { return mOutputTracks; }
   GPUhdi() const GPUdEdxInfo* OutputTracksdEdx() const { return mOutputTracksdEdx; }
   GPUhdi() GPUdEdxInfo* OutputTracksdEdx() { return mOutputTracksdEdx; }
+  GPUhdi() const GPUdEdxInfo* OutputTracksdEdxAlt() const { return mOutputTracksdEdxAlt; }
+  GPUhdi() GPUdEdxInfo* OutputTracksdEdxAlt() { return mOutputTracksdEdxAlt; }
   GPUhdi() uint32_t NClusters() const { return mNClusters; }
   GPUhdi() uint32_t NMaxClusters() const { return mNMaxClusters; }
   GPUhdi() uint32_t NMaxTracks() const { return mNMaxTracks; }
@@ -261,6 +263,7 @@ class GPUTPCGMMerger : public GPUProcessor
   int32_t mNClusters = 0;                           // Total number of incoming clusters (from sector tracks)
   GPUTPCGMMergedTrack* mOutputTracks = nullptr;     //* array of output merged tracks
   GPUdEdxInfo* mOutputTracksdEdx = nullptr;         //* dEdx information
+  GPUdEdxInfo* mOutputTracksdEdxAlt = nullptr;      //* dEdx alternative information
   GPUTPCGMSectorTrack* mSectorTrackInfos = nullptr; //* additional information for sector tracks
   int32_t* mSectorTrackInfoIndex = nullptr;
   GPUTPCGMMergedTrackHit* mClusters = nullptr;

GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx

@@ -114,6 +114,7 @@ GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::output>(int32_t nBlocks
   constexpr float MinDelta = 0.1f;
   const GPUTPCGMMergedTrack* tracks = merger.OutputTracks();
   GPUdEdxInfo* tracksdEdx = merger.OutputTracksdEdx();
+  GPUdEdxInfo* tracksdEdxAlt = merger.OutputTracksdEdxAlt();
   const int32_t nTracks = merger.NOutputTracksTPCO2();
   const GPUTPCGMMergedTrackHit* trackClusters = merger.Clusters();
   constexpr uint8_t flagsReject = getFlagsReject();
@@ -146,6 +147,7 @@ GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::output>(int32_t nBlocks
     auto& outerPar = tracks[i].OuterParam();
     if (merger.Param().par.dodEdx && merger.Param().dodEdxDownscaled) {
       oTrack.setdEdx(tracksdEdx[i]);
+      oTrack.setdEdxAlt(tracksdEdxAlt[i]);
     }
 
     auto snpOut = outerPar.P[2];

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -60,7 +60,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
 
   const GPUParam& GPUrestrict() param = merger->Param();
 
-  GPUdEdx dEdx;
+  GPUdEdx dEdx, dEdxAlt;
   GPUTPCGMPropagator prop;
   gputpcgmmergertypes::InterpolationErrors interpolation;
   prop.SetMaterialTPC();
@@ -220,6 +220,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
         dodEdx = AttachClustersPropagate(merger, cluster.sector, lastRow, cluster.row, iTrk, cluster.leg == clusters[maxN - 1].leg, prop, inFlyDirection, GPUCA_MAX_SIN_PHI, dodEdx);
         if (dodEdx) {
           dEdx.fillSubThreshold(lastRow - wayDirection);
+          dEdxAlt.fillSubThreshold(lastRow - wayDirection);
         }
       }
 
@@ -366,25 +367,33 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
           CADEBUG(printf("Reinit linearization\n"));
           prop.SetTrack(this, prop.GetAlpha());
         }
-        if (param.par.dodEdx && param.dodEdxDownscaled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg && (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMask) == 0) { // TODO: Costimize flag to remove, and option to remove double-clusters
-          float qtot = 0, qmax = 0, pad = 0, relTime = 0;
-          const int32_t clusterCount = (ihit - ihitMergeFirst) * wayDirection + 1;
-          for (int32_t iTmp = ihitMergeFirst; iTmp != ihit + wayDirection; iTmp += wayDirection) {
-            if (merger->GetConstantMem()->ioPtrs.clustersNative == nullptr) {
-              qtot += clustersXYZ[ihit].amp;
-            } else {
-              const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num];
-              qtot += cl.qTot;
-              qmax = CAMath::Max<float>(qmax, cl.qMax);
-              pad += cl.getPad();
-              relTime += cl.getTime();
+        if (param.par.dodEdx && param.dodEdxDownscaled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg) { // TODO: Costimize flag to remove, and option to remove double-clusters
+          bool acc = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMask) == 0, accAlt = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMaskAlt) == 0;
+          if (acc || accAlt) {
+            float qtot = 0, qmax = 0, pad = 0, relTime = 0;
+            const int32_t clusterCount = (ihit - ihitMergeFirst) * wayDirection + 1;
+            for (int32_t iTmp = ihitMergeFirst; iTmp != ihit + wayDirection; iTmp += wayDirection) {
+              if (merger->GetConstantMem()->ioPtrs.clustersNative == nullptr) {
+                qtot += clustersXYZ[ihit].amp;
+              } else {
+                const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num];
+                qtot += cl.qTot;
+                qmax = CAMath::Max<float>(qmax, cl.qMax);
+                pad += cl.getPad();
+                relTime += cl.getTime();
+              }
+            }
+            qtot /= clusterCount; // TODO: Weighted Average
+            pad /= clusterCount;
+            relTime /= clusterCount;
+            relTime = relTime - CAMath::Round(relTime);
+            if (acc) {
+              dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
+            }
+            if (accAlt) {
+              dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
             }
           }
-          qtot /= clusterCount; // TODO: Weighted Average
-          pad /= clusterCount;
-          relTime /= clusterCount;
-          relTime = relTime - CAMath::Round(relTime);
-          dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
         }
       } else if (retVal >= GPUTPCGMPropagator::updateErrorClusterRejected) { // cluster far away form the track
         if (allowModification) {
@@ -419,6 +428,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
 
   if (param.par.dodEdx && param.dodEdxDownscaled) {
     dEdx.computedEdx(merger->OutputTracksdEdx()[iTrk], param);
+    dEdxAlt.computedEdx(merger->OutputTracksdEdxAlt()[iTrk], param);
   }
   Alpha = prop.GetAlpha();
   MoveToReference(prop, param, Alpha);