Methods for layer-dependend mat.LUT rescaling

See the macro O2/Detectors/Base/test/rescaleLUT.C as an example, which writes a rescaled LUT to the <original_name>_rescaled.root file.
The rescaling parameters are provided as a set of RescRange structs: radial range + scaling coefficient.
It uses a method o2::base::MatLayerCylSet::scaleLayersByR(rmin,rmax, factor, bool scaleX2X0=true, bool scaleRho=true).
All LUT layers overlapping with rmin:rmax range will be rescaled. Alternatively, one can use directly the method o2::base::MatLayerCylSet::scaleLayersByID(lrIDmin,lrIDmax, factor, bool scaleX2X0=true, bool scaleRho=true)
to scale a set of layers [lrIDmin:lrIDmax].

Author: shahor02

Detectors/Base/CMakeLists.txt

@@ -87,6 +87,7 @@ endif()
 
 install(FILES test/buildMatBudLUT.C
               test/extractLUTLayers.C
+              test/rescaleLUT.C
               DESTINATION share/macro/)
 
 o2_add_test_root_macro(test/buildMatBudLUT.C
@@ -96,3 +97,7 @@ o2_add_test_root_macro(test/buildMatBudLUT.C
 o2_add_test_root_macro(test/extractLUTLayers.C
                        PUBLIC_LINK_LIBRARIES O2::DetectorsBase
                        LABELS detectorsbase)
+
+o2_add_test_root_macro(test/rescaleLUT.C
+                       PUBLIC_LINK_LIBRARIES O2::DetectorsBase
+                       LABELS detectorsbase)

Detectors/Base/include/DetectorsBase/MatLayerCyl.h

@@ -93,10 +93,12 @@ class MatLayerCyl : public o2::gpu::FlatObject
   GPUd() const MatCell& getCell(int iphiSlice, int iz) const { return mCells[getCellID(iphiSlice, iz)]; }
 
 #ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
-  GPUd() MatCell& getCellPhiBin(int iphi, int iz)
+  MatCell& getCellPhiBin(int iphi, int iz)
   {
     return mCells[getCellIDPhiBin(iphi, iz)];
   }
+
+  void scale(float factor, bool _x2x0 = true, bool _rho = true);
 #endif
 
   // ---------------------- Z slice manipulation

Detectors/Base/include/DetectorsBase/MatLayerCylSet.h

@@ -98,6 +98,10 @@ class MatLayerCylSet : public o2::gpu::FlatObject
     // get material budget traversed on the line between point0 and point1
     return getMatBudget(point0.X(), point0.Y(), point0.Z(), point1.X(), point1.Y(), point1.Z());
   }
+
+  void scaleLayersByID(int lrFrom, int lrTo, float factor, bool _x2x0 = true, bool _rho = true);
+  void scaleLayersByR(float rFrom, float rTo, float factor, bool _x2x0 = true, bool _rho = true);
+
 #endif // !GPUCA_ALIGPUCODE
   GPUd() MatBudget getMatBudget(float x0, float y0, float z0, float x1, float y1, float z1) const;
 

Detectors/Base/src/MatLayerCyl.cxx

@@ -319,9 +319,24 @@ void MatLayerCyl::flatten(char* newPtr)
   mConstructionMask = Constructed;
 }
 
+//________________________________________________________________________________
+void MatLayerCyl::scale(float factor, bool _x2x0, bool _rho)
+{
+  LOGP(info, "Scaling layer {:.3f}<r<{:.3f} by {:.3f}", getRMin(), getRMax(), factor);
+  for (int i = 0; i < mNPhiSlices * mNZBins; i++) {
+    if (_x2x0) {
+      mCells[i].meanX2X0 *= factor;
+    }
+    if (_rho) {
+      mCells[i].meanRho *= factor;
+    }
+  }
+}
+
 #endif // ! GPUCA_ALIGPUCODE
 
 #ifndef GPUCA_GPUCODE
+
 //________________________________________________________________________________
 void MatLayerCyl::fixPointers(char* oldPtr, char* newPtr)
 {

Detectors/Base/src/MatLayerCylSet.cxx

@@ -256,6 +256,33 @@ void MatLayerCylSet::print(bool data) const
          float(getFlatBufferSize()) / 1024 / 1024);
 }
 
+//________________________________________________________________________________
+void MatLayerCylSet::scaleLayersByID(int lrFrom, int lrTo, float factor, bool _x2x0, bool _rho)
+{
+  lrFrom = std::max(0, std::min(lrFrom, get()->mNLayers - 1));
+  lrTo = std::max(0, std::min(lrTo, get()->mNLayers - 1));
+  int dir = lrFrom >= lrTo ? -1 : 1;
+  lrTo += dir;
+  for (int i = lrFrom; i != lrTo; i += dir) {
+    get()->mLayers[i].scale(factor, _x2x0, _rho);
+  }
+}
+
+//________________________________________________________________________________
+void MatLayerCylSet::scaleLayersByR(float rFrom, float rTo, float factor, bool _x2x0, bool _rho)
+{
+  if (rFrom > rTo) {
+    std::swap(rFrom, rTo);
+  }
+  Ray ray(std::max(getRMin(), rFrom), 0., 0., std::min(getRMax(), rTo), 0., 0.);
+  short lmin, lmax;
+  if (!getLayersRange(ray, lmin, lmax)) {
+    LOGP(warn, "No layers found for {} < r < {}", rFrom, rTo);
+    return;
+  }
+  scaleLayersByID(lmin, lmax, factor, _x2x0, _rho);
+}
+
 #endif //!GPUCA_ALIGPUCODE
 
 #ifndef GPUCA_GPUCODE

Detectors/Base/test/rescaleLUT.C

@@ -0,0 +1,82 @@
+// 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.
+
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "DetectorsBase/MatLayerCylSet.h"
+#include "Framework/Logger.h"
+#include "CCDB/BasicCCDBManager.h"
+#include <regex>
+#endif
+
+// Macro to extract layers covering selected radial range into the separate LUT file.
+
+void rescaleLUT(o2::base::MatLayerCylSet* src, const std::string& outName)
+{
+  struct RescRange {
+    float rMin, rMax, factor;
+  };
+  std::vector<RescRange> task = {
+    // put here radial ranges in increasing order with corresponding factors to rescale
+    {0.1f, 6.f, 1.05},  // e.g. rescale layers covering 0.1<r<6 by factor 1.05
+    {30.f, 40.f, 1.15}, // e.g. rescale layers covering 30.f<r<40.f by factor 1.15
+  };
+
+  // check if there are no overlaps in ranges, to avoid double rescaling
+  for (size_t il = 1; il < task.size(); il++) {
+    short lmax, lmin;
+    float rmin = task[il - 1].rMax, rmax = task[il].rMin;
+    if (rmin > rmax) {
+      LOGP(error, "rMax={:.2f} of range {} is larger then rMin={:.2f} of range {}, must be in increasing order", rmin, il - 1, rmax, il);
+      return;
+    }
+    o2::base::Ray ray(std::max(src->getRMin(), rmin), 0., 0., std::min(src->getRMax(), rmax), 0., 0.);
+    if (!src->getLayersRange(ray, lmin, lmax)) {
+      LOGP(error, "No layers found for {:.2f} < r < {:.2f}", rmin, rmax);
+      return;
+    }
+    if (lmin == lmax) {
+      LOGP(error, "rMax={:.2f} of range {} and rMin={:.2f} of range {}, correspond to the same slice {} with {:.2f}<r<{:.2f}",
+           rmin, il - 1, rmax, il, lmin, src->getLayer(lmin).getRMin(), src->getLayer(lmin).getRMax());
+      return;
+    }
+  }
+
+  for (size_t il = 0; il < task.size(); il++) {
+    src->scaleLayersByR(task[il].rMin, task[il].rMax, task[il].factor);
+  }
+  if (outName.size()) {
+    src->writeToFile(outName);
+  }
+}
+
+void rescaleLUT(const std::string& fname)
+{
+  auto src = o2::base::MatLayerCylSet::loadFromFile(fname);
+  if (!src) {
+    LOGP(error, "failed to open source LUT from {}", fname);
+    return;
+  }
+  auto fnameOut = std::regex_replace(fname, std::regex(R"(.root)"), "_rescaled.root");
+  rescaleLUT(src, fnameOut);
+}
+
+void rescaleLUT(long timestamp = -1)
+{
+  auto& mg = o2::ccdb::BasicCCDBManager::instance();
+  mg.setTimestamp(timestamp);
+  auto src = o2::base::MatLayerCylSet::rectifyPtrFromFile(mg.get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
+  if (!src) {
+    LOGP(error, "failed to open load LUT from CCDB for timestamp {}", timestamp);
+    return;
+  }
+  auto fnameOut = fmt::format("matbudLUT_ts{}_rescaled.root", timestamp);
+  rescaleLUT(src, fnameOut);
+}