Debug streamer: Add downsampling by tsalis (#12112)

* Debug streamer: Add downsampling by tsalis

* Adding time series for vertex position

- adding time series for DCA to vertex for ITS-TPC tracks

* increase ClassDef for TimeSeriesITSTPC

* fix typo

Author: matthias-kleiner

Common/Utils/include/CommonUtils/DebugStreamer.h

@@ -50,6 +50,7 @@ enum SamplingTypes {
   sampleID = 2,       ///< sample every n IDs (per example track)
   sampleIDGlobal = 3, ///< in case different streamers have access to the same IDs use this gloabl ID
   sampleWeights = 4,  ///< perform sampling on weights, defined where the streamer is called
+  sampleTsalis = 5,   ///< perform sampling on tsalis pdf
 };
 
 #if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
@@ -163,10 +164,34 @@ class DebugStreamer
   /// \return returns integer index for given streamer flag
   static int getIndex(const StreamFlags streamFlag);
 
+  /// Random downsampling trigger function using Tsalis/Hagedorn spectra fit (sqrt(s) = 62.4 GeV to 13 TeV) as in https://iopscience.iop.org/article/10.1088/2399-6528/aab00f/pdf
+  /// \return flat q/pt trigger
+  /// \param pt pat of particle
+  /// \param factorPt defines the sampling
+  /// \param sqrts centre of mass energy
+  /// \param weight weight which is internally calculated
+  /// \param rnd random value between (0->1) used to check for sampling
+  /// \param mass particles mass (use pion if not known)
+  static bool downsampleTsalisCharged(float pt, float factorPt, float sqrts, float& weight, float rnd, float mass = 0.13957);
+
+  /// get random value between min and max
+  static float getRandom(float min = 0, float max = 1);
+
  private:
   using StreamersPerFlag = tbb::concurrent_unordered_map<size_t, std::unique_ptr<o2::utils::TreeStreamRedirector>>;
   StreamersPerFlag mTreeStreamer; ///< streamer which is used for the debugging
 
+  /// Tsalis/Hagedorn function describing charged pt spectra (m s = 62.4 GeV to 13 TeV) as in https://iopscience.iop.org/article/10.1088/2399-6528/aab00f/pdf
+  /// https://github.com/alisw/AliPhysics/blob/523f2dc8b45d913e9b7fda9b27e746819cbe5b09/PWGPP/AliAnalysisTaskFilteredTree.h#L145
+  /// \param pt     - transverse momentum
+  /// \param mass   - mass of particle
+  /// \param sqrts  - centre of mass energy
+  /// \return       - invariant yields of the charged particle *pt
+  ///    n(sqrts)= a + b/sqrt(s)                             - formula 6
+  ///    T(sqrts)= c + d/sqrt(s)                             - formula 7
+  ///    a = 6.81 ± 0.06       and b = 59.24 ± 3.53 GeV      - for charged particles page 3
+  ///    c = 0.082 ± 0.002 GeV and d = 0.151 ± 0.048 (GeV)   - for charged particles page 4
+  static float tsalisCharged(float pt, float mass, float sqrts);
 #else
 
   // empty implementation of the class for GPU or when the debug streamer is not build for CPU

Common/Utils/src/DebugStreamer.cxx

@@ -65,27 +65,22 @@ bool o2::utils::DebugStreamer::checkStream(const StreamFlags streamFlag, const s
   // check sampling frequency
   const auto sampling = getSamplingTypeFrequency(streamFlag);
   if (sampling.first != SamplingTypes::sampleAll) {
-    // init random number generator for each thread
-    static thread_local std::mt19937 generator(std::random_device{}());
-    std::uniform_real_distribution<> distr(0, 1);
-
     auto sampleTrack = [&]() {
       if (samplingID == -1) {
         LOGP(fatal, "Sampling type sampleID not supported for stream flag {}", streamFlag);
       }
-      std::uniform_real_distribution<> distr(0, 1);
       // sample on samplingID (e.g. track level)
       static thread_local std::unordered_map<StreamFlags, std::pair<size_t, bool>> idMap;
       // in case of first call samplingID in idMap is 0 and always false and first ID rejected
       if (idMap[streamFlag].first != samplingID) {
-        idMap[streamFlag] = std::pair<size_t, bool>{samplingID, (distr(generator) < sampling.second)};
+        idMap[streamFlag] = std::pair<size_t, bool>{samplingID, (getRandom() < sampling.second)};
       }
       return idMap[streamFlag].second;
     };
 
     if (sampling.first == SamplingTypes::sampleRandom) {
       // just sample randomly
-      return (distr(generator) < sampling.second);
+      return (getRandom() < sampling.second);
     } else if (sampling.first == SamplingTypes::sampleID) {
       return sampleTrack();
     } else if (sampling.first == SamplingTypes::sampleIDGlobal) {
@@ -106,12 +101,44 @@ bool o2::utils::DebugStreamer::checkStream(const StreamFlags streamFlag, const s
       }
     } else if (sampling.first == SamplingTypes::sampleWeights) {
       // sample with weight
-      return (weight * distr(generator) < sampling.second);
+      return (weight * getRandom() < sampling.second);
     }
   }
   return true;
 }
 
+float o2::utils::DebugStreamer::tsalisCharged(float pt, float mass, float sqrts)
+{
+  const float a = 6.81;
+  const float b = 59.24;
+  const float c = 0.082;
+  const float d = 0.151;
+  const float mt = std::sqrt(mass * mass + pt * pt);
+  const float n = a + b / sqrts;
+  const float T = c + d / sqrts;
+  const float p0 = n * T;
+  const float result = std::pow((1. + mt / p0), -n) * pt;
+  return result;
+}
+
+bool o2::utils::DebugStreamer::downsampleTsalisCharged(float pt, float factorPt, float sqrts, float& weight, float rnd, float mass)
+{
+  const float prob = tsalisCharged(pt, mass, sqrts);
+  const float probNorm = tsalisCharged(1., mass, sqrts);
+  weight = prob / probNorm;
+  const bool isSampled = (rnd * (weight * pt * pt)) < factorPt;
+  return isSampled;
+}
+
+float o2::utils::DebugStreamer::getRandom(float min, float max)
+{
+  // init random number generator for each thread
+  static thread_local std::mt19937 generator(std::random_device{}());
+  std::uniform_real_distribution<> distr(min, max);
+  const float rnd = distr(generator);
+  return rnd;
+}
+
 int o2::utils::DebugStreamer::getIndex(const StreamFlags streamFlag)
 {
   // see: https://stackoverflow.com/a/71539401

Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h

@@ -362,11 +362,28 @@ struct ITSTPC_Matching {
 };
 
 struct TimeSeriesITSTPC {
-  TimeSeries mTSTPC;                  ///< TPC standalone DCAs
-  TimeSeries mTSITSTPC;               ///< ITS-TPC standalone DCAs
-  ITSTPC_Matching mITSTPCAll;         ///< ITS-TPC matching efficiency for ITS standalone + afterburner
-  ITSTPC_Matching mITSTPCStandalone;  ///< ITS-TPC matching efficiency for ITS standalone
-  ITSTPC_Matching mITSTPCAfterburner; ///< ITS-TPC matchin efficiency  fir ITS afterburner
+  TimeSeries mTSTPC;                             ///< TPC standalone DCAs
+  TimeSeries mTSITSTPC;                          ///< ITS-TPC standalone DCAs
+  ITSTPC_Matching mITSTPCAll;                    ///< ITS-TPC matching efficiency for ITS standalone + afterburner
+  ITSTPC_Matching mITSTPCStandalone;             ///< ITS-TPC matching efficiency for ITS standalone
+  ITSTPC_Matching mITSTPCAfterburner;            ///< ITS-TPC matchin efficiency  fir ITS afterburner
+  std::vector<float> nPrimVertices;              ///< number of primary vertices
+  std::vector<float> nVertexContributors_Median; ///< number of primary vertices
+  std::vector<float> nVertexContributors_RMS;    ///< number of primary vertices
+  std::vector<float> vertexX_Median;             ///< vertex x position
+  std::vector<float> vertexY_Median;             ///< vertex y position
+  std::vector<float> vertexZ_Median;             ///< vertex z position
+  std::vector<float> vertexX_RMS;                ///< vertex x RMS
+  std::vector<float> vertexY_RMS;                ///< vertex y RMS
+  std::vector<float> vertexZ_RMS;                ///< vertex z RMS
+  std::vector<float> mDCAr_comb_A_Median;        ///< DCAr for ITS-TPC track - A-side
+  std::vector<float> mDCAz_comb_A_Median;        ///< DCAz for ITS-TPC track - A-side
+  std::vector<float> mDCAr_comb_A_RMS;           ///< DCAr RMS for ITS-TPC track - A-side
+  std::vector<float> mDCAz_comb_A_RMS;           ///< DCAz RMS for ITS-TPC track - A-side
+  std::vector<float> mDCAr_comb_C_Median;        ///< DCAr for ITS-TPC track - C-side
+  std::vector<float> mDCAz_comb_C_Median;        ///< DCAz for ITS-TPC track - C-side
+  std::vector<float> mDCAr_comb_C_RMS;           ///< DCAr RMS for ITS-TPC track - C-side
+  std::vector<float> mDCAz_comb_C_RMS;           ///< DCAz RMS for ITS-TPC track - C-side
 
   // dump this object to a tree
   void dumpToTree(const char* outFileName, const int nHBFPerTF = 32);
@@ -390,6 +407,8 @@ struct TimeSeriesITSTPC {
   bool areSameSize() const { return (mTSTPC.areSameSize() && mTSITSTPC.areSameSize()); } ///< check if stored currents have same number of entries
   bool isEmpty() const { return mTSTPC.isEmpty(); }                                      ///< check if values are empty
   size_t getEntries() const { return mTSTPC.getEntries(); }                              ///< \return returns number of values stored
+  void fill(const std::vector<float>& vecFrom, std::vector<float>& vecTo, const unsigned int posIndex) { std::copy(vecFrom.begin(), vecFrom.end(), vecTo.begin() + posIndex); }
+  void insert(std::vector<float>& vec, const std::vector<float>& vecTmp) { vec.insert(vec.begin(), vecTmp.begin(), vecTmp.end()); }
 
   /// acummulate integrated currents at given index
   /// \param posIndex index where data will be copied to
@@ -401,6 +420,25 @@ struct TimeSeriesITSTPC {
     mITSTPCAll.fill(posIndex, data.mITSTPCAll);
     mITSTPCStandalone.fill(posIndex, data.mITSTPCStandalone);
     mITSTPCAfterburner.fill(posIndex, data.mITSTPCAfterburner);
+    fill(data.mDCAr_comb_A_Median, mDCAr_comb_A_Median, posIndex);
+    fill(data.mDCAz_comb_A_Median, mDCAz_comb_A_Median, posIndex);
+    fill(data.mDCAr_comb_A_RMS, mDCAr_comb_A_RMS, posIndex);
+    fill(data.mDCAz_comb_A_RMS, mDCAz_comb_A_RMS, posIndex);
+    fill(data.mDCAr_comb_C_Median, mDCAr_comb_C_Median, posIndex);
+    fill(data.mDCAz_comb_C_Median, mDCAz_comb_C_Median, posIndex);
+    fill(data.mDCAr_comb_C_RMS, mDCAr_comb_C_RMS, posIndex);
+    fill(data.mDCAz_comb_C_RMS, mDCAz_comb_C_RMS, posIndex);
+
+    const int iTF = posIndex / mTSTPC.getNBins();
+    nPrimVertices[iTF] = data.nPrimVertices.front();
+    nVertexContributors_Median[iTF] = data.nVertexContributors_Median.front();
+    nVertexContributors_RMS[iTF] = data.nVertexContributors_RMS.front();
+    vertexX_Median[iTF] = data.vertexX_Median.front();
+    vertexY_Median[iTF] = data.vertexY_Median.front();
+    vertexZ_Median[iTF] = data.vertexZ_Median.front();
+    vertexX_RMS[iTF] = data.vertexX_RMS.front();
+    vertexY_RMS[iTF] = data.vertexY_RMS.front();
+    vertexZ_RMS[iTF] = data.vertexZ_RMS.front();
   }
 
   /// \param nDummyValues number of empty values which are inserted at the beginning of the accumulated integrated currents
@@ -411,6 +449,27 @@ struct TimeSeriesITSTPC {
     mITSTPCAll.insert(nDummyValues);
     mITSTPCStandalone.insert(nDummyValues);
     mITSTPCAfterburner.insert(nDummyValues);
+    std::vector<float> vecTmp(nDummyValues, 0);
+    insert(mDCAr_comb_A_Median, vecTmp);
+    insert(mDCAz_comb_A_Median, vecTmp);
+    insert(mDCAr_comb_A_RMS, vecTmp);
+    insert(mDCAz_comb_A_RMS, vecTmp);
+    insert(mDCAr_comb_C_Median, vecTmp);
+    insert(mDCAz_comb_C_Median, vecTmp);
+    insert(mDCAr_comb_C_RMS, vecTmp);
+    insert(mDCAz_comb_C_RMS, vecTmp);
+
+    const int nDummyValuesVtx = nDummyValues / mTSTPC.getNBins();
+    std::vector<float> vecTmpVtx(nDummyValuesVtx, 0);
+    insert(nPrimVertices, vecTmp);
+    insert(nVertexContributors_Median, vecTmp);
+    insert(nVertexContributors_RMS, vecTmp);
+    insert(vertexX_Median, vecTmp);
+    insert(vertexY_Median, vecTmp);
+    insert(vertexZ_Median, vecTmp);
+    insert(vertexX_RMS, vecTmp);
+    insert(vertexY_RMS, vecTmp);
+    insert(vertexZ_RMS, vecTmp);
   }
 
   /// resize buffer for accumulated currents
@@ -421,9 +480,28 @@ struct TimeSeriesITSTPC {
     mITSTPCAll.resize(nTotal);
     mITSTPCStandalone.resize(nTotal);
     mITSTPCAfterburner.resize(nTotal);
-  }
-
-  ClassDefNV(TimeSeriesITSTPC, 1);
+    mDCAr_comb_A_Median.resize(nTotal);
+    mDCAz_comb_A_Median.resize(nTotal);
+    mDCAr_comb_A_RMS.resize(nTotal);
+    mDCAz_comb_A_RMS.resize(nTotal);
+    mDCAr_comb_C_Median.resize(nTotal);
+    mDCAz_comb_C_Median.resize(nTotal);
+    mDCAr_comb_C_RMS.resize(nTotal);
+    mDCAz_comb_C_RMS.resize(nTotal);
+
+    const int nTotalVtx = nTotal / mTSTPC.getNBins();
+    nPrimVertices.resize(nTotalVtx);
+    nVertexContributors_Median.resize(nTotalVtx);
+    nVertexContributors_RMS.resize(nTotalVtx);
+    vertexX_Median.resize(nTotalVtx);
+    vertexY_Median.resize(nTotalVtx);
+    vertexZ_Median.resize(nTotalVtx);
+    vertexX_RMS.resize(nTotalVtx);
+    vertexY_RMS.resize(nTotalVtx);
+    vertexZ_RMS.resize(nTotalVtx);
+  }
+
+  ClassDefNV(TimeSeriesITSTPC, 2);
 };
 
 } // end namespace tpc