Merge branch 'AliceO2Group:dev' into new_geometry

Author: pkurash

Common/MathUtils/include/MathUtils/detail/trigonometric.h

@@ -275,12 +275,11 @@ GPUhdi() constexpr T fastATan2(T y, T x)
   // Average inaccuracy: 0.00048
   // Max inaccuracy: 0.00084
   // Speed: 6.2 times faster than atan2f()
-  constexpr T Pi = 3.1415926535897932384626433832795;
-
   auto atan = [](T a) -> T {
     // returns the arctan for the angular range [-Pi/4, Pi/4]
     // the polynomial coefficients are taken from:
     // https://stackoverflow.com/questions/42537957/fast-accurate-atan-arctan-approximation-algorithm
+    constexpr T Pi = 3.1415926535897932384626433832795;
     constexpr T A = 0.0776509570923569;
     constexpr T B = -0.287434475393028;
     constexpr T C = ((Pi / 4) - A - B);
@@ -290,6 +289,7 @@ GPUhdi() constexpr T fastATan2(T y, T x)
 
   auto atan2P = [atan](T yy, T xx) -> T {
     // fast atan2(yy,xx) for the angular range [0,+Pi]
+    constexpr T Pi = 3.1415926535897932384626433832795;
     constexpr T Pi025 = 1 * Pi / 4;
     constexpr T Pi075 = 3 * Pi / 4;
     const T x1 = xx + yy; //  point p1 (x1,y1) = (xx,yy) - Pi/4

DataFormats/Reconstruction/include/ReconstructionDataFormats/GlobalTrackID.h

@@ -165,7 +165,7 @@ GPUconstexpr() GlobalTrackID::mask_t sMasks[GlobalTrackID::NSources] = ///< dete
     GlobalTrackID::mask_t(math_utils::bit2Mask(GlobalTrackID::MCHMID)) // Temporary ordering
 };
 
-#ifndef __OPENCL__
+#if !defined(__OPENCL__) && !defined(__METAL__)
 GPUconstexpr() GlobalTrackID::mask_t sTrackSources =
   sMasks[GlobalTrackID::ITS] | sMasks[GlobalTrackID::TPC] | sMasks[GlobalTrackID::MFT] | sMasks[GlobalTrackID::MCH] | sMasks[GlobalTrackID::ITSTPC] | sMasks[GlobalTrackID::TPCTOF] |
   sMasks[GlobalTrackID::TPCTRD] | sMasks[GlobalTrackID::MFTMCH] | sMasks[GlobalTrackID::ITSTPCTRD] | sMasks[GlobalTrackID::ITSTPCTOF] | sMasks[GlobalTrackID::TPCTRDTOF] |
@@ -177,7 +177,7 @@ GPUconstexpr() GlobalTrackID::mask_t sTrackSources =
 GPUdi() constexpr GlobalTrackID::DetID::mask_t GlobalTrackID::getSourceDetectorsMask(int i) { return globaltrackid_internal::SourceDetectorsMasks[i]; }
 GPUdi() constexpr GlobalTrackID::mask_t GlobalTrackID::getSourceMask(int s) { return globaltrackid_internal::sMasks[s]; }
 
-#ifndef __OPENCL__
+#if !defined(__OPENCL__) && !defined(__METAL__)
 GPUdi() bool GlobalTrackID::isTrackSource(int s)
 {
   return globaltrackid_internal::sTrackSources[s];

Detectors/Base/include/DetectorsBase/GRPGeomHelper.h

@@ -91,7 +91,8 @@ struct GRPGeomRequest {
   enum GeomRequest { None,
                      Aligned,
                      Ideal,
-                     Alignments };
+                     Alignments,
+                     AlignmentsAndAligned };
 
   bool askAggregateRunInfo = false;
   bool askGRPECS = false;

Detectors/Base/src/GRPGeomHelper.cxx

@@ -39,14 +39,14 @@ namespace o2d = o2::dataformats;
 GRPGeomRequest::GRPGeomRequest(bool orbitResetTime, bool GRPECS, bool GRPLHCIF, bool GRPMagField, bool askMatLUT, GeomRequest geom, std::vector<o2::framework::InputSpec>& inputs, bool askOnce, bool needPropD, std::string detMaskString)
   : askGRPECS(GRPECS), askGRPLHCIF(GRPLHCIF), askGRPMagField(GRPMagField), askMatLUT(askMatLUT), askTime(orbitResetTime), askOnceAllButField(askOnce), needPropagatorD(needPropD)
 {
-  if (geom == Aligned) {
+  if (geom == Aligned || geom == AlignmentsAndAligned) {
     askGeomAlign = true;
     addInput({"geomAlg", "GLO", "GEOMALIGN", 0, Lifetime::Condition, ccdbParamSpec("GLO/Config/GeometryAligned")}, inputs);
   } else if (geom == Ideal || geom == Alignments) {
     askGeomIdeal = true;
     addInput({"geomIdeal", "GLO", "GEOMIDEAL", 0, Lifetime::Condition, ccdbParamSpec("GLO/Config/Geometry")}, inputs);
   }
-  if (geom == Alignments) {
+  if (geom == Alignments || geom == AlignmentsAndAligned) {
     askAlignments = true;
     o2::detectors::DetID::mask_t algDetMask = DetID::getMask(detMaskString);
     for (auto id = DetID::First; id <= DetID::Last; id++) {

Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt

@@ -9,7 +9,7 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
-#add_compile_options(-O0 -g -fPIC)
+# add_compile_options(-O0 -g -fPIC)
 
 o2_add_library(GlobalTrackingStudy
                TARGETVARNAME targetName

Detectors/GlobalTrackingWorkflow/study/src/CheckResidSpec.cxx

@@ -26,6 +26,7 @@
 #include "DataFormatsGlobalTracking/RecoContainer.h"
 #include "DataFormatsITSMFT/TrkClusRef.h"
 #include "DataFormatsITSMFT/ClustersPerLayer.h"
+#include "ITSMFTReconstruction/ChipMappingITS.h"
 #include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h"
 #include "ReconstructionDataFormats/TrackTPCITS.h"
 #include "ReconstructionDataFormats/GlobalTrackID.h"
@@ -106,6 +107,7 @@ class CheckResidSpec final : public Task
   void postProcessHistos();
   void drawHistos();
 
+  o2::itsmft::ChipMappingITS mChipMap;
   o2::globaltracking::RecoContainer* mRecoData = nullptr;
   int mNThreads = 1;
   bool mMeanVertexUpdated = false;
@@ -632,6 +634,7 @@ bool CheckResidSpec::refitITStrack(o2::track::TrackParCov& track, GTrackID gid)
 void CheckResidSpec::fillHistos(const o2::checkresid::Track& trc)
 {
   const auto& params = CheckResidConfig::Instance();
+  auto geom = o2::its::GeometryTGeo::Instance();
   int np = trc.points.size();
   auto pt = trc.track.getPt();
   if (pt < params.minPt || pt > params.maxPt) {
@@ -662,6 +665,21 @@ void CheckResidSpec::fillHistos(const o2::checkresid::Track& trc)
       mHMan->getHisto2F(il * 10 + 1 * 100 + 5 + 2000)->Fill(pt, pull);
       mHMan->getHisto2F(il * 10 + 1 * 100 + 5 + 3000)->Fill(trc.track.getTgl(), pull);
     }
+    if (pnt.lr >= 0 && pnt.sens >= 0 && geom->getLayer(pnt.sens) == pnt.lr) {
+      const int halfBarrel = geom->getHalfBarrel(pnt.sens);
+      const int chip = geom->getFirstChipIndex(pnt.lr) + geom->getChipIdInLayer(pnt.sens);
+      const int chipHistoID = (2 * il + halfBarrel) * 100 + 15000;
+      if (halfBarrel >= 0 && halfBarrel < 2) {
+        mHMan->getHisto2F(chipHistoID)->Fill(chip, pnt.dy);
+        mHMan->getHisto2F(chipHistoID + 10)->Fill(chip, pnt.dz);
+        if (pnt.sig2y > 0) {
+          mHMan->getHisto2F(chipHistoID + 5)->Fill(chip, pnt.dy / std::sqrt(pnt.sig2y));
+        }
+        if (pnt.sig2z > 0) {
+          mHMan->getHisto2F(chipHistoID + 15)->Fill(chip, pnt.dz / std::sqrt(pnt.sig2z));
+        }
+      }
+    }
   }
   //--------------
   if (trc.trIBOut.getX() > 1 && std::abs(trc.trIBOut.getX() - trc.trOBInw.getX()) < 0.1) {
@@ -726,6 +744,18 @@ void CheckResidSpec::bookHistos()
       mHMan->addHisto(htgl2, il * 10 + iyz * 100 + 3000);
       auto htgl2p = new TH2F(fmt::format("{}_{}_{}_pull", dname, lrName, "tgl").c_str(), fmt::format("pull {}_{{{}}} vs {};tg#lambda; pull{}", dtit, lrName, "tg#lambda", dtit).c_str(), params.nBinsTgl, -params.maxTgl, params.maxTgl, params.nBinsRes, -params.maxPull, params.maxPull);
       mHMan->addHisto(htgl2p, il * 10 + iyz * 100 + 5 + 3000);
+
+      if (il > 0) {
+        const int ilr = il - 1;
+        const int nChips = mChipMap.getNChipsOnLayer(ilr) / 2;
+        for (int ib = 0; ib < 2; ++ib) { // chips per halfbarrel
+          const int startChip = ib * nChips + mChipMap.getFirstChipsOnLayer(ilr);
+          auto h2 = new TH2F(fmt::format("{}_{}_{}", dname, lrName, ib).c_str(), fmt::format("{} Lr{} HB{};chip;{}", dname, lrName, ib, dtit).c_str(), nChips, startChip - 0.5, startChip - 0.5 + nChips, params.nBinsRes, -params.maxDYZ[il], params.maxDYZ[il]);
+          mHMan->addHisto(h2, (2 * il + ib) * 100 + iyz * 10 + 15000);
+          auto h2p = new TH2F(fmt::format("{}_{}_{}_pull", dname, lrName, ib).c_str(), fmt::format("{} Lr{} HB{} pull;chip;pull{}", dname, lrName, ib, dtit).c_str(), nChips, startChip - 0.5, startChip - 0.5 + nChips, params.nBinsRes, -params.maxPull, params.maxPull);
+          mHMan->addHisto(h2p, (2 * il + ib) * 100 + iyz * 10 + 5 + 15000);
+        }
+      }
     }
   }
 
@@ -805,14 +835,23 @@ void CheckResidSpec::postProcessHistos()
         }
       }
     }
-    histm->write();
+    for (int il = 1; il < 8; il++) {
+      for (int ib = 0; ib < 2; ib++) {
+        for (int iyz = 0; iyz < 2; iyz++) {
+          for (int iht = 0; iht < 2; iht++) {
+            const int offsV = iht == 0 ? 0 : 5;
+            fitSlices((2 * il + ib) * 100 + iyz * 10 + offsV + 15000);
+          }
+        }
+      }
+    }
   }
   delete gs;
 }
 
 void CheckResidSpec::drawHistos()
 {
-  gROOT->SetBatch(true);
+  // gROOT->SetBatch(true);
   gStyle->SetTitleX(0.2);
   gStyle->SetTitleY(0.88);
   gStyle->SetTitleW(0.25);
@@ -928,6 +967,8 @@ void CheckResidSpec::drawHistos()
   }
   TCanvas cly("cly", "", 600, 800), clz("clz", "", 600, 800), clpar("clpar", "", 600, 800);
   TCanvas czly("czly", "", 600, 800), czlz("czlz", "", 600, 800), czlpar("czlpar", "", 600, 800);
+  TCanvas cchip("cchip", "", 600, 800);
+  constexpr std::array<const char*, 4> chipTitles{"#DeltaY", "pull #DeltaY", "#DeltaZ", "pull #DeltaZ"};
   auto AddLabel = [](const char* txt, float x = 0.1, float y = 0.9, int color = kBlack, float size = 0.04) {
     TLatex* lt = new TLatex(x, y, txt);
     lt->SetNDC();
@@ -1062,6 +1103,94 @@ void CheckResidSpec::drawHistos()
   AddLabel("IB-OB tracks params differences at R = 12 cm", 0.2, 0.8);
   czlpar.Print(Form("%s_hman.pdf", params.outname.c_str()));
 
+  auto getHMinMax = [this](int id, float marginFrac, float maxDev, int roundToMicrons) -> std::tuple<float, float, float> {
+    float maxV = std::numeric_limits<float>::min(), minV = std::numeric_limits<float>::max();
+    const int nh = this->mHManV.size();
+    for (int j = 0; j < nh; j++) {
+      if (auto histo = this->mHManV[j]->getHisto(id); histo != nullptr) {
+        const int nb = histo->GetNbinsX();
+        for (int ib = 1; ib <= nb; ib++) {
+          float v = histo->GetBinContent(ib);
+          maxV = std::max(maxV, v);
+          minV = std::min(minV, v);
+        }
+      }
+    }
+    if (maxV < minV) {
+      return {0.f, 1.f, 0.1f};
+    }
+    auto delta = maxV - minV;
+    maxV += delta * marginFrac;
+    minV -= delta * marginFrac;
+    if (roundToMicrons > 0) {
+      maxV = ((int(std::ceil(maxV * 1e4)) / roundToMicrons + 1) * roundToMicrons) * 1e-4;
+      minV = ((int(std::floor(minV * 1e4)) / roundToMicrons - 1) * roundToMicrons) * 1e-4;
+    }
+    minV = std::max(minV, -maxDev);
+    maxV = std::min(maxV, maxDev);
+    return {minV, maxV, std::max(std::abs(minV), std::abs(maxV))};
+  };
+  auto collectSpread = [](TH1* hin, TH1* hout) {
+    const int nb = hin->GetNbinsX();
+    for (int ib = 1; ib <= nb; ib++) {
+      if (hin->GetBinError(ib) != 0) {
+        hout->Fill(hin->GetBinContent(ib));
+      }
+    }
+    hout->SetLineColor(hin->GetLineColor());
+    hout->SetMarkerColor(hin->GetMarkerColor());
+    hout->SetMarkerStyle(hin->GetMarkerStyle());
+  };
+  int nh = this->mHManV.size();
+  for (int il = 1; il < 8; ++il) {
+    for (int ib = 0; ib < 2; ib++) {
+      cchip.Clear();
+      cchip.Divide(1, 4);
+      for (int iyz = 0; iyz < 2; iyz++) {
+        for (int iv = 0; iv < 2; iv++) { // mean or sigma
+          const float yRange = iyz == 0 ? params.resMMLrY[il] : params.resMMLrZ[il];
+          bool drawn = false;
+          const int id = 15000 + iyz * 10 + (2 * il + ib) * 100 + 1 + iv;
+          auto [minV, maxV, absMax] = getHMinMax(id, 0.15, (iv == 0 ? 2 : 4) * yRange, 10);
+          TH1* firstAuto = nullptr;
+          float maxAuto = 0.f;
+          for (int j = 0; j < nh; j++) {
+            auto hman = this->mHManV[j].get();
+            auto histo = hman->getHisto(id);
+            if (histo) {
+              histo->SetMaximum(absMax);
+              histo->SetMinimum(iv ? 0 : -absMax);
+              if (il > 3) {
+                histo->SetMarkerSize(0.3);
+                histo->SetLineStyle(3);
+              }
+              cchip.cd(2 * iv + 1);
+              histo->Draw(drawn ? "same" : "");
+              if (!drawn) {
+                gPad->SetGrid();
+              }
+              auto hmeans = new TH1F(Form("%s_spread", histo->GetName()), Form("%s spread", histo->GetTitle()), params.nBinsRes, iv == 0 ? minV : std::max(0.f, minV), maxV);
+              collectSpread(histo, hmeans);
+              hman->addHisto(hmeans, id + 2);
+              cchip.cd(2 * iv + 2);
+              if (!drawn) {
+                firstAuto = hmeans;
+              }
+              maxAuto = std::max(maxAuto, (float)hmeans->GetMaximum());
+              hmeans->Draw(drawn ? "same" : "");
+              drawn = true;
+            }
+          }
+          if (firstAuto) {
+            firstAuto->SetMaximum(maxAuto * 1.25);
+            gPad->SetGrid();
+            gPad->Modified();
+          }
+        }
+        cchip.Print(Form("%s_hman.pdf", params.outname.c_str()));
+      }
+    }
+  }
   cly.Print(Form("%s_hman.pdf]", params.outname.c_str()));
 }
 
@@ -1074,6 +1203,13 @@ void CheckResidSpec::endOfStream(EndOfStreamContext& ec)
   if (mDraw) {
     drawHistos();
   }
+  if (auto nh = mHManV.size(); nh > 0) {
+    for (int ihm = 0; ihm < nh; ihm++) {
+      if (mHManV[ihm]) {
+        mHManV[ihm].get()->write();
+      }
+    }
+  }
 }
 
 void CheckResidSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualWriterSpec.h

@@ -22,7 +22,7 @@ namespace tpc
 {
 
 /// create a processor spec
-framework::DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput);
+framework::DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput, bool writeTrackDataMC);
 
 } // namespace tpc
 } // namespace o2

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCUnbinnedResidualReaderSpec.h

@@ -30,26 +30,28 @@ namespace tpc
 class TPCUnbinnedResidualReader : public o2::framework::Task
 {
  public:
-  TPCUnbinnedResidualReader(bool trkInput) : mTrackInput(trkInput){};
+  TPCUnbinnedResidualReader(bool trkInput, bool mcInput) : mTrackInput(trkInput), mMCInput(mcInput) {};
   ~TPCUnbinnedResidualReader() override = default;
   void init(o2::framework::InitContext& ic) final;
   void run(o2::framework::ProcessingContext& pc) final;
 
  private:
   void connectTree();
   bool mTrackInput = false;
+  bool mMCInput = false;
   std::unique_ptr<TFile> mFileIn;
   std::unique_ptr<TTree> mTreeIn;
   std::string mInFileName;
   std::string mInTreeName;
   std::vector<UnbinnedResid> mUnbinnedResid, *mUnbinnedResidPtr = &mUnbinnedResid;
   std::vector<DetInfoResid> mDetInfoUnbRes, *mDetInfoUnbResPtr = &mDetInfoUnbRes;
   std::vector<TrackData> mTrackData, *mTrackDataPtr = &mTrackData;
+  std::vector<TrackDataMC> mTrackDataMC, *mTrackDataMCPtr = &mTrackDataMC;
   std::vector<TrackDataCompact> mTrackDataCompact, *mTrackDataCompactPtr = &mTrackDataCompact;
 };
 
-/// read unbinned TPC residuals and reference tracks from a root file
-framework::DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput);
+/// read unbinned TPC residuals and reference tracks (with mcInput also their MC truth) from a root file
+framework::DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput, bool mcInput);
 
 } // namespace tpc
 } // namespace o2

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx

@@ -29,7 +29,7 @@ namespace tpc
 template <typename T>
 using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition<T>;
 
-DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput)
+DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput, bool writeTrackDataMC)
 {
   bool writeUnfiltered = SpacePointsCalibConfParam::Instance().writeUnfiltered;
   return MakeRootTreeWriterSpec("tpc-residuals-writer",
@@ -39,6 +39,7 @@ DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput
                                 BranchDefinition<std::vector<DetInfoResid>>{InputSpec{"detInfo", "GLO", "DETINFORES"}, "detInfo"},
                                 BranchDefinition<std::vector<TrackDataCompact>>{InputSpec{"trackRefs", "GLO", "TRKREFS"}, "trackRefs"},
                                 BranchDefinition<std::vector<TrackData>>{InputSpec{"tracks", "GLO", "TRKDATA"}, "tracks", (writeTrackData ? 1 : 0)},
+                                BranchDefinition<std::vector<TrackDataMC>>{InputSpec{"tracksMC", "GLO", "TRKDATAMC"}, "tracksMC", (writeTrackDataMC ? 1 : 0)},
                                 BranchDefinition<std::vector<TrackDataExtended>>{InputSpec{"trackExt", "GLO", "TRKDATAEXT"}, "trackExt", (debugOutput ? 1 : 0)})();
 }
 

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCUnbinnedResidualReaderSpec.cxx

@@ -52,6 +52,12 @@ void TPCUnbinnedResidualReader::connectTree()
   if (mTrackInput) {
     mTreeIn->SetBranchAddress("tracks", &mTrackDataPtr);
   }
+  if (mMCInput) {
+    if (!mTreeIn->GetBranch("tracksMC")) {
+      LOG(fatal) << "No tracksMC branch found in the unbinned residuals tree, the interpolation workflow must run with --enable-mc";
+    }
+    mTreeIn->SetBranchAddress("tracksMC", &mTrackDataMCPtr);
+  }
   LOG(info) << "Loaded tree from " << mInFileName << " with " << mTreeIn->GetEntries() << " entries";
 }
 
@@ -68,14 +74,17 @@ void TPCUnbinnedResidualReader::run(ProcessingContext& pc)
     LOG(info) << "Pushing " << mTrackData.size() << " reference tracks for these residuals";
     pc.outputs().snapshot(Output{"GLO", "TRKDATA", 0}, mTrackData);
   }
+  if (mMCInput) {
+    pc.outputs().snapshot(Output{"GLO", "TRKDATAMC", 0}, mTrackDataMC);
+  }
 
   if (mTreeIn->GetReadEntry() + 1 >= mTreeIn->GetEntries()) {
     pc.services().get<ControlService>().endOfStream();
     pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
   }
 }
 
-DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput)
+DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput, bool mcInput)
 {
   std::vector<OutputSpec> outputs;
   outputs.emplace_back("GLO", "UNBINNEDRES", 0, Lifetime::Timeframe);
@@ -84,12 +93,15 @@ DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput)
   if (trkInput) {
     outputs.emplace_back("GLO", "TRKDATA", 0, Lifetime::Timeframe);
   }
+  if (mcInput) {
+    outputs.emplace_back("GLO", "TRKDATAMC", 0, Lifetime::Timeframe);
+  }
 
   return DataProcessorSpec{
     "TPCUnbinnedResidualReader",
     Inputs{},
     outputs,
-    AlgorithmSpec{adaptFromTask<TPCUnbinnedResidualReader>(trkInput)},
+    AlgorithmSpec{adaptFromTask<TPCUnbinnedResidualReader>(trkInput, mcInput)},
     Options{
       {"input-filename", VariantType::String, "o2residuals_tpc.root", {"Name of the input file"}},
       {"input-dir", VariantType::String, "none", {"Input directory"}},