Merge branch 'AliceO2Group:dev' into cholmcc_aod_to_hepmc

Author: Christian Holm Christensen

Common/MathUtils/CMakeLists.txt

@@ -16,6 +16,7 @@ o2_add_library(
           src/Chebyshev3D.cxx
           src/Chebyshev3DCalc.cxx
    src/SymMatrixSolver.cxx
+   src/Tsallis.cxx
   PUBLIC_LINK_LIBRARIES
     ROOT::Hist
     FairLogger::FairLogger
@@ -37,7 +38,8 @@ o2_target_root_dictionary(
           include/MathUtils/RandomRing.h
           include/MathUtils/Primitive2D.h
           include/MathUtils/SMatrixGPU.h
-          include/MathUtils/SymMatrixSolver.h)
+          include/MathUtils/SymMatrixSolver.h
+          include/MathUtils/Tsallis.h)
 
 o2_add_test(
   CachingTF1

Common/MathUtils/include/MathUtils/Tsallis.h

@@ -0,0 +1,48 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef ALICEO2_TSALLIS_H
+#define ALICEO2_TSALLIS_H
+
+namespace o2
+{
+
+namespace math_utils
+{
+
+struct Tsallis {
+  /// Tsallis/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 tsallisCharged(float pt, float mass, float sqrts);
+
+  /// Random downsampling trigger function using Tsallis/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 downsampleTsallisCharged(float pt, float factorPt, float sqrts, float& weight, float rnd, float mass = 0.13957);
+};
+
+} // namespace math_utils
+} // namespace o2
+
+#endif

Common/MathUtils/src/MathUtilsLinkDef.h

@@ -41,5 +41,6 @@
 #pragma link C++ class o2::math_utils::IntervalXYd_t + ;
 #pragma link C++ class o2::math_utils::Bracketf_t + ;
 #pragma link C++ class o2::math_utils::Bracketd_t + ;
+#pragma link C++ class o2::math_utils::Tsallis + ;
 
 #endif

Common/MathUtils/src/Tsallis.cxx

@@ -0,0 +1,39 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "MathUtils/Tsallis.h"
+#include <cmath>
+
+namespace o2::math_utils
+{
+
+float Tsallis::tsallisCharged(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;
+  return std::pow((1. + mt / p0), -n) * pt;
+}
+
+bool Tsallis::downsampleTsallisCharged(float pt, float factorPt, float sqrts, float& weight, float rnd, float mass)
+{
+  const float prob = tsallisCharged(pt, mass, sqrts);
+  const float probNorm = tsallisCharged(1., mass, sqrts);
+  weight = prob / probNorm;
+  return (rnd * (weight * pt * pt)) < factorPt;
+}
+
+} // namespace o2::math_utils

Common/Utils/include/CommonUtils/DebugStreamer.h

@@ -50,7 +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
+  sampleTsallis = 5,  ///< perform sampling on tsallis pdf
 };
 
 #if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
@@ -164,34 +164,12 @@ 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

@@ -107,29 +107,6 @@ bool o2::utils::DebugStreamer::checkStream(const StreamFlags streamFlag, const s
   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

DataFormats/common/include/CommonDataFormat/InteractionRecord.h

@@ -109,7 +109,7 @@ struct InteractionRecord {
 
   float differenceInBCMUS(const InteractionRecord& other) const
   {
-    // return difference in bunch-crossings in ms
+    // return difference in bunch-crossings in us
     return differenceInBC(other) * o2::constants::lhc::LHCBunchSpacingMUS;
   }
 

Detectors/AOD/src/AODProducerWorkflowSpec.cxx

@@ -319,7 +319,7 @@ void AODProducerWorkflowDPL::addToTracksExtraTable(TracksExtraCursorType& tracks
   // extra
   tracksExtraCursor(truncateFloatFraction(extraInfoHolder.tpcInnerParam, mTrack1Pt),
                     extraInfoHolder.flags,
-                    extraInfoHolder.itsClusterMap,
+                    extraInfoHolder.itsClusterSizes,
                     extraInfoHolder.tpcNClsFindable,
                     extraInfoHolder.tpcNClsFindableMinusFound,
                     extraInfoHolder.tpcNClsFindableMinusCrossedRows,
@@ -2370,14 +2370,14 @@ AODProducerWorkflowDPL::TrackExtraInfo AODProducerWorkflowDPL::processBarrelTrac
     int nClusters = itsTrack.getNClusters();
     float chi2 = itsTrack.getChi2();
     extraInfoHolder.itsChi2NCl = nClusters != 0 ? chi2 / (float)nClusters : 0;
-    extraInfoHolder.itsClusterMap = itsTrack.getPattern();
+    extraInfoHolder.itsClusterSizes = itsTrack.getClusterSizes();
     if (src == GIndex::ITS) { // standalone ITS track should set its time from the ROF
       const auto& rof = data.getITSTracksROFRecords()[mITSROFs[trackIndex.getIndex()]];
       double t = rof.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS + mITSROFrameHalfLengthNS;
       setTrackTime(t, mITSROFrameHalfLengthNS, false);
     }
   } else if (contributorsGID[GIndex::Source::ITSAB].isIndexSet()) { // this is an ITS-TPC afterburner contributor
-    extraInfoHolder.itsClusterMap = data.getITSABRefs()[contributorsGID[GIndex::Source::ITSAB].getIndex()].pattern;
+    extraInfoHolder.itsClusterSizes = data.getITSABRefs()[contributorsGID[GIndex::Source::ITSAB].getIndex()].getClusterSizes();
   }
   if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
     const auto& tpcOrig = data.getTPCTrack(contributorsGID[GIndex::TPC]);

Detectors/CTP/reconstruction/src/RawDataDecoder.cxx

@@ -384,19 +384,28 @@ int RawDataDecoder::shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digit
     } else if (lut == 2) { // L0
       shiftNew(dig.first, TFOrbit, inpmask, L0shift, 0, digitsMapShifted);
       if (dig.second.CTPClassMask.count()) {
-        LOG(error) << "Adding class mask without input ?";
+        // LOG(error) << "Adding class mask without input ?";
+        //  This is not needed as it can happen; Full checj done below - see next LOG(error)
         CTPDigit digi = {dig.first, 0, dig.second.CTPClassMask};
         digitsMapShifted[dig.first] = digi;
       }
     } else if (lut == 4) { // L1
       shiftNew(dig.first, TFOrbit, inpmask, L1shift, 1, digitsMapShifted);
+      if (dig.second.CTPClassMask.count()) {
+        CTPDigit digi = {dig.first, 0, dig.second.CTPClassMask};
+        digitsMapShifted[dig.first] = digi;
+      }
     } else if (lut == 6) { // L0 and L1
       shiftNew(dig.first, TFOrbit, inpmask, L0shift, 0, digitsMapShifted);
       shiftNew(dig.first, TFOrbit, inpmask, L1shift, 1, digitsMapShifted);
+      if (dig.second.CTPClassMask.count()) {
+        CTPDigit digi = {dig.first, 0, dig.second.CTPClassMask};
+        digitsMapShifted[dig.first] = digi;
+      }
     } else if (lut == 3) { // LM and L0
       shiftNew(dig.first, TFOrbit, inpmask, L0shift, 0, digitsMapShifted);
       CTPDigit digi = {dig.first, inpmask & (~L0MASKInputs), dig.second.CTPClassMask};
-      // LOG(info) << "LM-L0 present";
+      // if LM level do not need to add class as LM is not shifted;
       digitsMapShifted[dig.first] = digi;
     } else if (lut == 5) { // LM and L1
       shiftNew(dig.first, TFOrbit, inpmask, L1shift, 1, digitsMapShifted);
@@ -432,7 +441,7 @@ int RawDataDecoder::shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digit
     digits.push_back(dig.second);
   }
   if (nTwoI) { // Trigger class wo Input
-    LOG(error) << "LM:" << nLM << " L0:" << nL0 << " L1:" << nL1 << " TwI:" << nTwI << " Trigger cals wo inputTwoI:" << nTwoI;
+    LOG(error) << "LM:" << nLM << " L0:" << nL0 << " L1:" << nL1 << " TwI:" << nTwI << " Trigger classes wo input:" << nTwoI;
   }
   return 0;
 }

Detectors/Calibration/include/DetectorsCalibration/TimeSlotCalibration.h

@@ -489,6 +489,16 @@ bool TimeSlotCalibration<Container>::saveLastSlot()
   if (!updateSaveMetaData()) {
     return false;
   }
+
+  if (!mSaveDirectory.empty() && !std::filesystem::exists(mSaveDirectory)) {
+    std::filesystem::create_directories(mSaveDirectory);
+    if (!std::filesystem::exists(mSaveDirectory)) {
+      LOGP(fatal, "could not create output directory {}", mSaveDirectory);
+    } else {
+      LOGP(info, "created calibration directory {}", mSaveDirectory);
+    }
+  }
+
   auto pth = getSaveFilePath();
   auto pthTmp = pth + ".part";
   TFile flout(pthTmp.c_str(), "recreate");