Refactor track quality assessment and trigger logic - fix TRD extra information

Updated track quality assessment logic to conditionally retain only dEdx information and adjusted related calculations for tracklet merging based on trigger records.

Author: KangkanGoswami

Detectors/AOD/src/AODProducerWorkflowSpec.cxx

@@ -365,10 +365,10 @@ void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACurs
 {
   tracksQACursor(
     trackQAInfoHolder.trackID,
-    truncateFloatFraction(trackQAInfoHolder.tpcTime0, mTPCTime0),
+    mTrackQCRetainOnlydEdx ? 0.0f : truncateFloatFraction(trackQAInfoHolder.tpcTime0, mTPCTime0),
     truncateFloatFraction(trackQAInfoHolder.tpcdEdxNorm, mTrackSignal),
-    trackQAInfoHolder.tpcdcaR,
-    trackQAInfoHolder.tpcdcaZ,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int16_t>::min() : trackQAInfoHolder.tpcdcaR,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int16_t>::min() : trackQAInfoHolder.tpcdcaZ,
     trackQAInfoHolder.tpcClusterByteMask,
     trackQAInfoHolder.tpcdEdxMax0R,
     trackQAInfoHolder.tpcdEdxMax1R,
@@ -378,18 +378,18 @@ void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACurs
     trackQAInfoHolder.tpcdEdxTot1R,
     trackQAInfoHolder.tpcdEdxTot2R,
     trackQAInfoHolder.tpcdEdxTot3R,
-    trackQAInfoHolder.dRefContY,
-    trackQAInfoHolder.dRefContZ,
-    trackQAInfoHolder.dRefContSnp,
-    trackQAInfoHolder.dRefContTgl,
-    trackQAInfoHolder.dRefContQ2Pt,
-    trackQAInfoHolder.dRefGloY,
-    trackQAInfoHolder.dRefGloZ,
-    trackQAInfoHolder.dRefGloSnp,
-    trackQAInfoHolder.dRefGloTgl,
-    trackQAInfoHolder.dRefGloQ2Pt,
-    trackQAInfoHolder.dTofdX,
-    trackQAInfoHolder.dTofdZ);
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefContY,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefContZ,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefContSnp,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefContTgl,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefContQ2Pt,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefGloY,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefGloZ,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefGloSnp,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefGloTgl,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dRefGloQ2Pt,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dTofdX,
+    mTrackQCRetainOnlydEdx ? std::numeric_limits<int8_t>::min() : trackQAInfoHolder.dTofdZ);
 }
 
 template <typename TRDsExtraCursorType>
@@ -433,7 +433,7 @@ void AODProducerWorkflowDPL::addToTRDsExtra(const o2::globaltracking::RecoContai
 
     auto tphi = trkC.getSnp() / std::sqrt((1.f - trkC.getSnp()) * (1.f + trkC.getSnp()));
     auto trackletLength = std::sqrt(1.f + tphi * tphi + trkC.getTgl() * trkC.getTgl());
-    float cor = mTRDLocalGain->getValue(tracklet.getHCID() / 2, tracklet.getPadCol(), tracklet.getPadRow()) * trackletLength;
+    float cor = mTRDLocalGain->getValue(tracklet.getHCID() / 2, tracklet.getPadCol(), tracklet.getPadRow()) * mTRDGainCalib->getMPVdEdx(tracklet.getDetector()) / o2::trd::constants::MPVDEDXDEFAULT * trackletLength;
     q0s[iLay] = tracklet.getQ0();
     q1s[iLay] = tracklet.getQ1();
     q2s[iLay] = tracklet.getQ2();
@@ -443,14 +443,56 @@ void AODProducerWorkflowDPL::addToTRDsExtra(const o2::globaltracking::RecoContai
     ttgls[iLay] = trkC.getTgl();
     tphis[iLay] = tphi;
 
-    // z-row merging
+    // z-row merging, we want to merge only with tracklets from the same trigger record
     if (trk.getIsCrossingNeighbor(iLay) && trk.getHasNeighbor()) {
-      for (const auto& trklt : trklets) {
-        if (tracklet.getTrackletWord() == trklt.getTrackletWord()) {
+      // find the trigger the tracklet belongs to
+      auto trigsTRD = recoData.getTRDTriggerRecords();
+      size_t trdSelID = -1;
+
+      const auto& trig = trigsTRD[mCurrentTRDTrigID];
+      bool foundTRDTrigger = false;
+      // first check current trigger
+      if (trkltId >= trig.getFirstTracklet() && trkltId < trig.getFirstTracklet() + trig.getNumberOfTracklets()) {
+        trdSelID = mCurrentTRDTrigID;
+        foundTRDTrigger = true;
+      } else {
+        // then check next trigger
+        if (mCurrentTRDTrigID < trigsTRD.size() - 1) {
+          const auto& trig = trigsTRD[mCurrentTRDTrigID + 1];
+          if (trkltId >= trig.getFirstTracklet() && trkltId < trig.getFirstTracklet() + trig.getNumberOfTracklets()) {
+            trdSelID = mCurrentTRDTrigID + 1;
+            foundTRDTrigger = true;
+          }
+        }
+      }
+
+      size_t low = 0, up = trigsTRD.size() - 1;
+
+      // otherwise binary search
+      while (low <= up && !foundTRDTrigger) {
+        trdSelID = low + std::floor((up - low) / 2);
+        const auto& trig = trigsTRD[trdSelID];
+        if (trig.getFirstTracklet() > trkltId) {
+          up = trdSelID - 1;
+        } else {
+          if (trig.getFirstTracklet() + trig.getNumberOfTracklets() <= trkltId) {
+            low = trdSelID + 1;
+          } else {
+            foundTRDTrigger = true;
+          }
+        }
+      }
+      //-------------------
+      mCurrentTRDTrigID = trdSelID;
+      const auto& trigSel = trigsTRD[trdSelID];
+
+      // loop on other tracklets from the same trigger record
+      for (const auto& trklt : trklets.subspan(trigSel.getFirstTracklet(), trigSel.getNumberOfTracklets())) {
+        if (tracklet.getTrackletWord() == trklt.getTrackletWord() || tracklet.getDetector() != trklt.getDetector()) {
           continue;
         }
         if (std::abs(tracklet.getPadCol() - trklt.getPadCol()) <= 1 && std::abs(tracklet.getPadRow() - trklt.getPadRow()) == 1) {
-          cor = mTRDLocalGain->getValue(trklt.getHCID() / 2, trklt.getPadCol(), trklt.getPadRow()) * trackletLength;
+          cor = mTRDLocalGain->getValue(trklt.getHCID() / 2, trklt.getPadCol(), trklt.getPadRow()) * mTRDGainCalib->getMPVdEdx(tracklet.getDetector()) / o2::trd::constants::MPVDEDXDEFAULT * trackletLength;
           q0s[iLay] += trklt.getQ0();
           q1s[iLay] += trklt.getQ1();
           q2s[iLay] += trklt.getQ2();
@@ -574,7 +616,7 @@ void AODProducerWorkflowDPL::fillTrackTablesPerCollision(int collisionID,
 
           float weight = 0;
           static std::uniform_real_distribution<> distr(0., 1.);
-          bool writeQAData = o2::math_utils::Tsallis::downsampleTsallisCharged(data.getTrackParam(trackIndex).getPt(), mTrackQCFraction, mSqrtS, weight, distr(mGenerator));
+          bool writeQAData = o2::math_utils::Tsallis::downsampleTsallisCharged(data.getTrackParam(trackIndex).getPt(), mTrackQCFraction, mSqrtS, weight, distr(mGenerator)) || ((src != GIndex::TPC || mGIDUsedBySVtx.find(trackIndex) != mGIDUsedBySVtx.end() || mGIDUsedByStr.find(trackIndex) != mGIDUsedByStr.end()) && mTrackQCKeepGlobalTracks);
           auto extraInfoHolder = processBarrelTrack(collisionID, collisionBC, trackIndex, data, bcsMap);
 
           if (writeQAData) {
@@ -1140,7 +1182,8 @@ void AODProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader&
                              source == 0, // background
                              mMcParticleW,
                              mMcParticleMom,
-                             mMcParticlePos);
+                             mMcParticlePos,
+                             mUseSigFiltMC);
 
     mcReader.releaseTracksForSourceAndEvent(source, event);
   }
@@ -1795,6 +1838,8 @@ void AODProducerWorkflowDPL::init(InitContext& ic)
       LOGP(warn, "Specified non-default empty streamer mask!");
     }
   }
+  mTrackQCKeepGlobalTracks = ic.options().get<bool>("trackqc-keepglobaltracks");
+  mTrackQCRetainOnlydEdx = ic.options().get<bool>("trackqc-retainonlydedx");
   mTrackQCFraction = ic.options().get<float>("trackqc-fraction");
   mTrackQCNTrCut = ic.options().get<int64_t>("trackqc-NTrCut");
   mTrackQCDCAxy = ic.options().get<float>("trackqc-tpc-dca");
@@ -2215,9 +2260,7 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc)
                          0,
                          sourceID);
         }
-        if (sourceID != 0 || !mUseSigFiltMC) {
-          mcColToEvSrc.emplace_back(std::vector<int>{iCol, sourceID, eventID}); // point background and injected signal events to one collision
-        }
+        mcColToEvSrc.emplace_back(std::vector<int>{iCol, sourceID, eventID}); // point background and injected signal events to one collision
       }
     }
   }
@@ -2382,6 +2425,7 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc)
     }
   }
 
+  mCurrentTRDTrigID = 0; // reinitialize index for TRD trigger record search
   // filling unassigned tracks first
   // so that all unassigned tracks are stored in the beginning of the table together
   auto& trackRef = primVer2TRefs.back(); // references to unassigned tracks are at the end
@@ -2390,6 +2434,7 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc)
                               ambigTracksCursor, mftTracksCursor, mftTracksCovCursor, ambigMFTTracksCursor,
                               fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, fwdTrkClsCursor, bcsMap);
 
+  mCurrentTRDTrigID = 0; // reinitialize index for TRD trigger record search
   // filling collisions and tracks into tables
   collisionID = 0;
   collisionsCursor.reserve(primVertices.size());
@@ -2735,6 +2780,12 @@ AODProducerWorkflowDPL::TrackExtraInfo AODProducerWorkflowDPL::processBarrelTrac
     if (tpcOrig.getdEdx().dEdxTotTPC == 0) {
       extraInfoHolder.flags |= o2::aod::track::TPCdEdxAlt;
     }
+    if (tpcOrig.hasASideClusters()) {
+      extraInfoHolder.flags |= o2::aod::track::TPCSideA;
+    }
+    if (tpcOrig.hasCSideClusters()) {
+      extraInfoHolder.flags |= o2::aod::track::TPCSideC;
+    }
     extraInfoHolder.tpcInnerParam = tpcOrig.getP() / tpcOrig.getAbsCharge();
     extraInfoHolder.tpcChi2NCl = tpcOrig.getNClusters() ? tpcOrig.getChi2() / tpcOrig.getNClusters() : 0;
     extraInfoHolder.tpcSignal = dEdx.dEdxTotTPC;
@@ -3096,6 +3147,7 @@ void AODProducerWorkflowDPL::updateTimeDependentParams(ProcessingContext& pc)
     if (mEnableTRDextra) {
       mTRDLocalGain = pc.inputs().get<o2::trd::LocalGainFactor*>("trdlocalgainfactors").get();
       mTRDNoiseMap = pc.inputs().get<o2::trd::NoiseStatusMCM*>("trdnoisemap").get();
+      mTRDGainCalib = pc.inputs().get<o2::trd::CalGain*>("trdgaincalib").get(); // time dependent gain
     }
   }
   if (mPropTracks) {
@@ -3368,11 +3420,8 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo
     OutputForTable<Collisions>::spec(),
     OutputForTable<Decay3Bodys>::spec(),
     OutputForTable<FDDs>::spec(),
-    OutputForTable<FDDsExtra>::spec(),
     OutputForTable<FT0s>::spec(),
-    OutputForTable<FT0sExtra>::spec(),
     OutputForTable<FV0As>::spec(),
-    OutputForTable<FV0AsExtra>::spec(),
     OutputForTable<StoredFwdTracks>::spec(),
     OutputForTable<StoredFwdTracksCov>::spec(),
     OutputForTable<StoredMFTTracks>::spec(),
@@ -3410,6 +3459,7 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo
     outputs.push_back(OutputForTable<TRDsExtra>::spec());
     dataRequest->inputs.emplace_back("trdlocalgainfactors", "TRD", "LOCALGAINFACTORS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/LocalGainFactor"));
     dataRequest->inputs.emplace_back("trdnoisemap", "TRD", "NOISEMAP", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/NoiseMapMCM"));
+    dataRequest->inputs.emplace_back("trdgaincalib", "TRD", "CALGAIN", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/CalGain"));
   }
 
   if (useMC) {
@@ -3454,6 +3504,8 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo
       ConfigParamSpec{"hepmc-update", VariantType::String, "always", {"When to update HepMC Aux tables: always - force update, never - never update, all - if all keys are present, any - when any key is present (not valid yet)"}},
       ConfigParamSpec{"propagate-muons", VariantType::Bool, false, {"Propagate muons to IP"}},
       ConfigParamSpec{"thin-tracks", VariantType::Bool, false, {"Produce thinned track tables"}},
+      ConfigParamSpec{"trackqc-keepglobaltracks", VariantType::Bool, false, {"Always keep TrackQA for global tracks"}},
+      ConfigParamSpec{"trackqc-retainonlydedx", VariantType::Bool, false, {"Keep only dEdx information, zero out everything else"}},
       ConfigParamSpec{"trackqc-fraction", VariantType::Float, float(0.1), {"Fraction of tracks to QC"}},
       ConfigParamSpec{"trackqc-NTrCut", VariantType::Int64, 4L, {"Minimal length of the track - in amount of tracklets"}},
       ConfigParamSpec{"trackqc-tpc-dca", VariantType::Float, 3.f, {"Keep TPC standalone track with this DCAxy to the PV"}},