[Common] adding median time for occupancy estimation (#14202)

Author: altsybee

Common/DataModel/EventSelection.h

@@ -65,6 +65,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(FoundFDD, foundFDD, int, FDDs, "_foundFDD");  //!
 DECLARE_SOA_INDEX_COLUMN_FULL(FoundZDC, foundZDC, int, Zdcs, "_foundZDC");  //! ZDC entry index in ZDCs table (-1 if doesn't exist)
 DECLARE_SOA_COLUMN(NumTracksInTimeRange, trackOccupancyInTimeRange, int);   //! Occupancy in specified time interval by a number of tracks from nearby collisions // o2-linter: disable=name/o2-column
 DECLARE_SOA_COLUMN(SumAmpFT0CInTimeRange, ft0cOccupancyInTimeRange, float); //! Occupancy in specified time interval by a sum of FT0C amplitudes from nearby collisions // o2-linter: disable=name/o2-column
+DECLARE_SOA_COLUMN(OccupancyMedianTime, occupancyMedianTime, float);        //! Median time for the occupancy in specified time interval (the delta-time wrt to the collision at which occupancy crosses 1/2 of its full value)
 } // namespace evsel
 
 // bc-joinable event selection decisions
@@ -85,7 +86,8 @@ DECLARE_SOA_TABLE(EvSels, "AOD", "EVSEL", //!
                   evsel::FoundFDDId,
                   evsel::FoundZDCId,
                   evsel::NumTracksInTimeRange,
-                  evsel::SumAmpFT0CInTimeRange);
+                  evsel::SumAmpFT0CInTimeRange,
+                  evsel::OccupancyMedianTime);
 using EvSel = EvSels::iterator;
 } // namespace o2::aod
 

Common/TableProducer/eventSelection.cxx

@@ -694,7 +694,7 @@ struct EventSelectionTask {
       }
     }
 
-    evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, 0, 0);
+    evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, 0, 0, 0);
   }
   PROCESS_SWITCH(EventSelectionTask, processRun2, "Process Run2 event selection", true);
 
@@ -750,7 +750,7 @@ struct EventSelectionTask {
         int32_t foundFDD = bc.foundFDDId();
         int32_t foundZDC = bc.foundZDCId();
         uint32_t rct = 0;
-        evsel(bc.alias_raw(), bc.selection_raw(), rct, kFALSE, kFALSE, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, -1, -1);
+        evsel(bc.alias_raw(), bc.selection_raw(), rct, kFALSE, kFALSE, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, -1, -1, -1);
       }
       return;
     }
@@ -1166,7 +1166,7 @@ struct EventSelectionTask {
       }
 
       evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC,
-            vNumTracksITS567inFullTimeWin[colIndex], vSumAmpFT0CinFullTimeWin[colIndex]);
+            vNumTracksITS567inFullTimeWin[colIndex], vSumAmpFT0CinFullTimeWin[colIndex], 0);
     }
   }
 

Common/Tools/EventSelectionModule.h

@@ -96,7 +96,7 @@ struct bcselConfigurables : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<int> confTimeFrameStartBorderMargin{"TimeFrameStartBorderMargin", -1, "Number of bcs to cut at the start of the Time Frame. Take from CCDB if -1"}; // o2-linter: disable=name/configurable (temporary fix)
   o2::framework::Configurable<int> confTimeFrameEndBorderMargin{"TimeFrameEndBorderMargin", -1, "Number of bcs to cut at the end of the Time Frame. Take from CCDB if -1"};       // o2-linter: disable=name/configurable (temporary fix)
   o2::framework::Configurable<bool> confCheckRunDurationLimits{"checkRunDurationLimits", false, "Check if the BCs are within the run duration limits"};                           // o2-linter: disable=name/configurable (temporary fix)
-  o2::framework::Configurable<std::vector<int>> maxInactiveChipsPerLayer{"maxInactiveChipsPerLayer", {8, 8, 8, 111, 111, 195, 195}, "Maximum allowed number of inactive ITS chips per layer"};
+  o2::framework::Configurable<std::vector<int>> confMaxInactiveChipsPerLayer{"maxInactiveChipsPerLayer", {8, 8, 8, 111, 111, 195, 195}, "Maximum allowed number of inactive ITS chips per layer"};
   o2::framework::Configurable<int> confNumberOfOrbitsPerTF{"NumberOfOrbitsPerTF", -1, "Number of orbits per Time Frame. Take from CCDB if -1"}; // o2-linter: disable=name/configurable (temporary fix)
 };
 
@@ -121,6 +121,8 @@ struct evselConfigurables : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<bool> confUseWeightsForOccupancyVariable{"UseWeightsForOccupancyEstimator", 1, "Use or not the delta-time weights for the occupancy estimator"};                                   // o2-linter: disable=name/configurable (temporary fix)
   o2::framework::Configurable<int> confNumberOfOrbitsPerTF{"NumberOfOrbitsPerTF", -1, "Number of orbits per Time Frame. Take from CCDB if -1"};                                                                  // o2-linter: disable=name/configurable (temporary fix)
 
+  // o2::framework::Configurable<std::vector<float>> confTimeIntervalsForGranularOccupancy{"timeIntervalsForGranularOccupancy", {-100, -60, -30, -10, 0, 10, 30, 60, 100}, "Delta-time intervals (wrt given collision) for which we store occupancy"};
+
   // configurables for light-ion event selection
   o2::framework::Configurable<float> confLightIonsNsigmaOnVzDiff{"VzDiffNsigma", 3.0, "+/- nSigma on vZ difference by FT0 and by tracks"};                  // o2-linter: disable=name/configurable (temporary fix)
   o2::framework::Configurable<float> confLightIonsMarginVzDiff{"VzDiffMargin", 0.2, "margin for +/- nSigma cut on vZ difference by FT0 and by tracks, cm"}; // o2-linter: disable=name/configurable (temporary fix)
@@ -561,14 +563,14 @@ class BcSelectionModule
         LOGP(debug, "orbit: {}, previous orbit: {}, next orbit: {} ", orbit, prevOrbitForInactiveChips, nextOrbitForInactiveChips);
         LOGP(debug, "next inactive chips: {} {} {} {} {} {} {}", vNextInactiveChips[0], vNextInactiveChips[1], vNextInactiveChips[2], vNextInactiveChips[3], vNextInactiveChips[4], vNextInactiveChips[5], vNextInactiveChips[6]);
         LOGP(debug, "prev inactive chips: {} {} {} {} {} {} {}", vPrevInactiveChips[0], vPrevInactiveChips[1], vPrevInactiveChips[2], vPrevInactiveChips[3], vPrevInactiveChips[4], vPrevInactiveChips[5], vPrevInactiveChips[6]);
-        isGoodITSLayer3 = vPrevInactiveChips[3] <= bcselOpts.maxInactiveChipsPerLayer->at(3) && vNextInactiveChips[3] <= bcselOpts.maxInactiveChipsPerLayer->at(3);
+        isGoodITSLayer3 = vPrevInactiveChips[3] <= bcselOpts.confMaxInactiveChipsPerLayer->at(3) && vNextInactiveChips[3] <= bcselOpts.confMaxInactiveChipsPerLayer->at(3);
         isGoodITSLayer0123 = true;
         for (int i = 0; i < 4; i++) { // o2-linter: disable=magic-number (counting first 4 ITS layers)
-          isGoodITSLayer0123 &= vPrevInactiveChips[i] <= bcselOpts.maxInactiveChipsPerLayer->at(i) && vNextInactiveChips[i] <= bcselOpts.maxInactiveChipsPerLayer->at(i);
+          isGoodITSLayer0123 &= vPrevInactiveChips[i] <= bcselOpts.confMaxInactiveChipsPerLayer->at(i) && vNextInactiveChips[i] <= bcselOpts.confMaxInactiveChipsPerLayer->at(i);
         }
         isGoodITSLayersAll = true;
         for (int i = 0; i < o2::itsmft::ChipMappingITS::NLayers; i++) {
-          isGoodITSLayersAll &= vPrevInactiveChips[i] <= bcselOpts.maxInactiveChipsPerLayer->at(i) && vNextInactiveChips[i] <= bcselOpts.maxInactiveChipsPerLayer->at(i);
+          isGoodITSLayersAll &= vPrevInactiveChips[i] <= bcselOpts.confMaxInactiveChipsPerLayer->at(i) && vNextInactiveChips[i] <= bcselOpts.confMaxInactiveChipsPerLayer->at(i);
         }
       }
 
@@ -875,7 +877,7 @@ class EventSelectionModule
         }
       }
 
-      evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, 0, 0);
+      evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, 0, 0, 0);
     }
   } // end processRun2
 
@@ -924,7 +926,7 @@ class EventSelectionModule
         int32_t foundFDD = bcselbuffer[bc.globalIndex()].foundFDDId;
         int32_t foundZDC = bcselbuffer[bc.globalIndex()].foundZDCId;
         uint32_t rct = 0;
-        evsel(bcselbuffer[bc.globalIndex()].alias, bcselbuffer[bc.globalIndex()].selection, rct, kFALSE, kFALSE, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, -1, -1);
+        evsel(bcselbuffer[bc.globalIndex()].alias, bcselbuffer[bc.globalIndex()].selection, rct, kFALSE, kFALSE, foundBC, foundFT0, foundFV0, foundFDD, foundZDC, -1, -1, -1);
       }
       return;
     }
@@ -1138,9 +1140,10 @@ class EventSelectionModule
 
     // pre-loop for occupancy calculation
     std::vector<bool> vIsFullInfoForOccupancy(cols.size(), 0);               // info for occupancy in +/- windows is available (i.e. a given coll is not too close to the TF borders)
-    std::vector<bool> vIsCollAtROFborder(cols.size(), 0);                    // collision is close to ITS ROF border
     std::vector<bool> vIsCollRejectedByTFborderCut(cols.size(), 0);          // helper vector with
     std::vector<bool> vCanHaveAssocCollsWithinLastDriftTime(cols.size(), 0); // to see if for some collisions in the occupancy calc (that are close to TF border) we will switch to FT0C based occupancy estimation
+    std::vector<int> vProxyForCollNtracks(cols.size(), 0);                   // n tracks, or - if a collision is close to ITS ROF border - scaled FT0C mult
+    std::vector<float> vMedianTimeForOccupancy(cols.size(), 0);              // median time for the occupancy in a defined time window
 
     const float timeWinOccupancyCalcMinNS = evselOpts.confTimeIntervalForOccupancyCalculationMin * 1e3; // ns
     const float timeWinOccupancyCalcMaxNS = evselOpts.confTimeIntervalForOccupancyCalculationMax * 1e3; // ns
@@ -1153,10 +1156,9 @@ class EventSelectionModule
 
       int32_t foundBC = vFoundBCindex[colIndex];
       auto bcselEntry = bcselbuffer[foundBC];
-      // check if we are close to ROF or TF borders => N tracks are not reliable, but FT0 can be used for occupancy estimation
-      if (!bitcheck64(bcselEntry.selection, aod::evsel::kNoITSROFrameBorder)) {
-        vIsCollAtROFborder[colIndex] = true;
-      }
+      // check if we are close to ITS ROF borders => N ITS tracks is not reliable, and FT0C ampl can be used for occupancy estimation
+      // denominator for vAmpFT0CperColl is the approximate conversion factor b/n FT0C ampl and number of PV tracks after cuts
+      vProxyForCollNtracks[colIndex] = bitcheck64(bcselEntry.selection, aod::evsel::kNoITSROFrameBorder) ? vTracksITS567perColl[colIndex] : vAmpFT0CperColl[colIndex] / 10.;
 
       if (!bitcheck64(bcselEntry.selection, aod::evsel::kNoTimeFrameBorder)) {
         vIsCollRejectedByTFborderCut[colIndex] = true;
@@ -1239,6 +1241,7 @@ class EventSelectionModule
       // ### for occupancy in time windows
       std::vector<int> vAssocToThisCol;
       std::vector<float> vCollsTimeDeltaWrtGivenColl;
+      std::vector<float> vProxyNtracksAssocColls; // mult of associated collisions, for the median time calc
       // find all collisions in time window before the current one
       minColIndex = colIndex - 1;
       while (minColIndex >= 0) {
@@ -1253,6 +1256,7 @@ class EventSelectionModule
           break;
         vAssocToThisCol.push_back(minColIndex);
         vCollsTimeDeltaWrtGivenColl.push_back(dt);
+        vProxyNtracksAssocColls.push_back(vProxyForCollNtracks[minColIndex]);
         minColIndex--;
       }
       // find all collisions in time window after the current one
@@ -1267,10 +1271,31 @@ class EventSelectionModule
           break;
         vAssocToThisCol.push_back(maxColIndex);
         vCollsTimeDeltaWrtGivenColl.push_back(dt);
+        vProxyNtracksAssocColls.push_back(vProxyForCollNtracks[maxColIndex]);
         maxColIndex++;
       }
       vCollsInTimeWin.push_back(vAssocToThisCol);
       vTimeDeltaForColls.push_back(vCollsTimeDeltaWrtGivenColl);
+
+      // calculation of the median time for the occupancy in a given time window
+      std::vector<std::pair<float, float>> pairsDeltaTimeMult;
+      for (size_t iCol = 0; iCol < vCollsTimeDeltaWrtGivenColl.size(); iCol++)
+        pairsDeltaTimeMult.emplace_back(vCollsTimeDeltaWrtGivenColl[iCol], vProxyNtracksAssocColls[iCol]);
+      std::sort(pairsDeltaTimeMult.begin(), pairsDeltaTimeMult.end()); // sorts by first element by default
+      int proxyTotalMultInTimeWin = std::accumulate(vProxyNtracksAssocColls.begin(), vProxyNtracksAssocColls.end(), 0);
+
+      float sumMult = 0.0;
+      for (size_t iCol = 0; iCol < vCollsTimeDeltaWrtGivenColl.size(); iCol++) {
+        sumMult += pairsDeltaTimeMult[iCol].second;
+        if (sumMult > proxyTotalMultInTimeWin / 2.0) {
+          vMedianTimeForOccupancy[colIndex] = pairsDeltaTimeMult[iCol].first / 1e3; // ns -> us
+          break;
+        }
+      }
+      for (size_t iCol = 0; iCol < vCollsTimeDeltaWrtGivenColl.size(); iCol++) {
+        LOGP(debug, "dt={} mult={}", pairsDeltaTimeMult[iCol].first, pairsDeltaTimeMult[iCol].second);
+      }
+      LOGP(debug, "   --> median time = {}", vMedianTimeForOccupancy[colIndex]);
     }
 
     // perform the occupancy calculation per ITS ROF and also in the pre-defined time window
@@ -1330,15 +1355,11 @@ class EventSelectionModule
       for (uint32_t iCol = 0; iCol < vAssocToThisCol.size(); iCol++) {
         int thisColIndex = vAssocToThisCol[iCol];
         float dt = vCollsTimeDeltaWrtGivenColl[iCol] / 1e3; // ns -> us
-
-        // check if we are close to ITS ROF borders => N ITS tracks is not reliable, and FT0C ampl can be used for occupancy estimation
-        // denominator for vAmpFT0CperColl is the approximate conversion factor b/n FT0C ampl and number of PV tracks after cuts
-        int nItsTracksAssocColl = !vIsCollAtROFborder[thisColIndex] ? vTracksITS567perColl[thisColIndex] : vAmpFT0CperColl[thisColIndex] / 10.;
         // counting tracks from other collisions in fixed time windows
         if (std::fabs(dt) < evselOpts.confTimeRangeVetoOnCollNarrow)
-          nITS567tracksForVetoNarrow += nItsTracksAssocColl;
+          nITS567tracksForVetoNarrow += vProxyForCollNtracks[thisColIndex];
         if (std::fabs(dt) < evselOpts.confTimeRangeVetoOnCollStrict)
-          nITS567tracksForVetoStrict += nItsTracksAssocColl;
+          nITS567tracksForVetoStrict += vProxyForCollNtracks[thisColIndex];
 
         // veto on high-mult collisions nearby, where artificial structures in the dt-occupancy plots are observed
         if (dt > -4.0 && dt < 2.0 && vAmpFT0CperColl[thisColIndex] > evselOpts.confFT0CamplCutVetoOnCollInTimeRange) { // dt in us // o2-linter: disable=magic-number
@@ -1359,7 +1380,7 @@ class EventSelectionModule
             // weighted occupancy
             wOccup = calcWeightForOccupancy(dt);
           }
-          nITS567tracksInFullTimeWindow += wOccup * nItsTracksAssocColl;
+          nITS567tracksInFullTimeWindow += wOccup * vProxyForCollNtracks[thisColIndex];
           sumAmpFT0CInFullTimeWindow += wOccup * vAmpFT0CperColl[thisColIndex];
         }
       }
@@ -1410,7 +1431,7 @@ class EventSelectionModule
       vNoCollInTimeRangeNarrow[colIndex] = (nITS567tracksForVetoNarrow == 0);
       vNoCollInTimeRangeStrict[colIndex] = (nITS567tracksForVetoStrict == 0);
       vNoHighMultCollInTimeRange[colIndex] = (nCollsWithFT0CAboveVetoStandard == 0 && nITS567tracksForVetoNarrow == 0);
-    }
+    } // end of the occupancy calculation
 
     for (const auto& col : cols) {
       int32_t colIndex = col.globalIndex();
@@ -1490,7 +1511,7 @@ class EventSelectionModule
       }
 
       evsel(alias, selection, rct, sel7, sel8, foundBC, foundFT0, foundFV0, foundFDD, foundZDC,
-            vNumTracksITS567inFullTimeWin[colIndex], vSumAmpFT0CinFullTimeWin[colIndex]);
+            vNumTracksITS567inFullTimeWin[colIndex], vSumAmpFT0CinFullTimeWin[colIndex], vMedianTimeForOccupancy[colIndex]);
     }
   } // end processRun3
 }; // end EventSelectionModule