Merge branch 'master' into mc-sampling-fixes

Author: EmilGorm

ALICE3/Core/CMakeLists.txt

@@ -20,6 +20,7 @@ o2physics_add_library(ALICE3Core
 o2physics_target_root_dictionary(ALICE3Core
                       HEADERS TrackUtilities.h
                               FlatLutEntry.h
+                              OTFParticle.h
                               GeometryContainer.h
                       LINKDEF ALICE3CoreLinkDef.h)
 

ALICE3/Macros/drawFastTracker.C

@@ -9,6 +9,10 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+/// \file   drawFastTracker.C
+/// \author Nicolò Jacazio nicolo.jacazio@cern.ch
+/// \brief  Draw FastTracker performace
+
 #include "ALICE3/Core/FastTracker.h"
 #include "ALICE3/Core/TrackUtilities.h"
 
@@ -26,6 +30,7 @@
 #include <TH1.h>
 #include <TLatex.h>
 #include <TLorentzVector.h>
+#include <TPDGCode.h>
 #include <TParticlePDG.h>
 #include <TString.h>
 
@@ -34,67 +39,64 @@
 #include <cstddef>
 #include <vector>
 
-void drawFastTracker(float magneticField = 5.f, // in units of kGauss
-                     const int nch = 100,       // number of charged particles per unit rapidity
-                     const int pdg = 211)       // PDG code of the particle to track
+void drawFastTracker(const float magneticField = 20.f, // in units of kGauss
+                     const int nch = 100,              // number of charged particles per unit rapidity
+                     const int pdg = PDG_t::kPiPlus)   // PDG code of the particle to track
 {
   TDatabasePDG* db = TDatabasePDG::Instance();
-  TParticlePDG* p = 0;
+  TParticlePDG* p = nullptr;
   p = db->GetParticle(pdg);
   if (!p) {
     LOG(fatal) << "Particle with PDG code " << pdg << " not found in TDatabasePDG";
     return;
   }
-  const float mass = p->Mass();      // particle mass in GeV/c^2
-  const float q = p->Charge() / 3.0; // charge in e
+  const float mass = p->Mass();                    // particle mass in GeV/c^2
+  const int q = static_cast<int>(p->Charge() / 3); // charge in e
 
   o2::parameters::GRPMagField grpmag;
   grpmag.setFieldUniformity(true);
   grpmag.setL3Current(30000.f * (magneticField / 5.0f));
-  auto field = grpmag.getNominalL3Field();
+  // auto field = grpmag.getNominalL3Field();
   o2::base::Propagator::initFieldFromGRP(&grpmag);
 
   fair::Logger::SetVerbosity(fair::Verbosity::verylow);
   o2::fastsim::FastTracker fastTracker = o2::fastsim::FastTracker();
-  if (0) {
-    fastTracker.SetApplyEffCorrection(false);
-    Double_t x0IB = 0.001;
-    Double_t x0OB = 0.01;
-    Double_t xrhoIB = 2.3292e-02; // 100 mum Si
-    Double_t xrhoOB = 2.3292e-01; // 1000 mum Si
-    Double_t resRPhiIB = 0.00025;
-    Double_t resZIB = 0.00025;
-    Double_t resRPhiOB = 0.00100;
-    Double_t resZOB = 0.00100;
-    Double_t eff = 0.98;
-    // fastTracker.AddLayer("vertex", 0.0, 250, 0, 0);                // dummy vertex for matrix calculation
-    fastTracker.AddLayer("bpipe0", 0.48, 250, 0.00042, 2.772e-02); // 150 mum Be
-    fastTracker.AddLayer("B00", 0.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1);
-    fastTracker.AddLayer("B01", 1.20, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1);
-    fastTracker.AddLayer("B02", 2.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1);
-    fastTracker.AddLayer("bpipe1", 3.7, 250, 0.0014, 9.24e-02); // 500 mum Be
-    fastTracker.AddLayer("B03", 3.75, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B04", 7.00, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B05", 12.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B06", 20.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B07", 30.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B08", 45.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B09", 60.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B10", 80.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-    fastTracker.AddLayer("B11", 100., 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1);
-  } else {
-    std::vector<float> pixelRes{0.025, 0.025, 0.01, 0.01};
-    // fastTracker.AddSiliconALICE3v4(pixelRes); // FIXME
-  }
-
-  fastTracker.Print();
-  fastTracker.SetMagneticField(magneticField);
+  fastTracker.setApplyEffCorrection(false);
+  Double_t x0IB = 0.001;
+  Double_t x0OB = 0.01;
+  Double_t xrhoIB = 2.3292e-02; // 100 mum Si
+  Double_t xrhoOB = 2.3292e-01; // 1000 mum Si
+  Double_t resRPhiIB = 0.00025;
+  Double_t resZIB = 0.00025;
+  Double_t resRPhiOB = 0.00100;
+  Double_t resZOB = 0.00100;
+  Double_t trackerEff = 0.98;
+  // fastTracker.AddLayer("vertex", 0.0, 250, 0, 0);                // dummy vertex for matrix calculation
+
+  // N.B. Tracker configuration outdated and is only meant to be used as a template
+  fastTracker.addLayer("bpipe0", 0.48, 250, 0.00042, 2.772e-02); // 150 mum Be
+  fastTracker.addLayer("B00", 0.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, trackerEff, 1);
+  fastTracker.addLayer("B01", 1.20, 250, x0IB, xrhoIB, resRPhiIB, resZIB, trackerEff, 1);
+  fastTracker.addLayer("B02", 2.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, trackerEff, 1);
+  fastTracker.addLayer("bpipe1", 3.7, 250, 0.0014, 9.24e-02); // 500 mum Be
+  fastTracker.addLayer("B03", 3.75, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B04", 7.00, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B05", 12.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B06", 20.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B07", 30.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B08", 45.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B09", 60.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B10", 80.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+  fastTracker.addLayer("B11", 100., 250, x0OB, xrhoOB, resRPhiOB, resZOB, trackerEff, 1);
+
+  fastTracker.print();
+  fastTracker.setMagneticField(magneticField);
 
   TAxis ptBinning(1000, 0., 10);
-  TGraph* gPt = new TGraph();
+  auto gPt = new TGraph();
   gPt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/c)");
   gPt->GetYaxis()->SetTitle("Efficiency");
-  TEfficiency* hEfficiency = new TEfficiency("hEfficiency", ";#it{p}_{T} (GeV/c);Efficiency", ptBinning.GetNbins(), ptBinning.GetBinLowEdge(1), ptBinning.GetBinUpEdge(ptBinning.GetNbins()));
+  auto hEfficiency = new TEfficiency("hEfficiency", ";#it{p}_{T} (GeV/c);Efficiency", ptBinning.GetNbins(), ptBinning.GetBinLowEdge(1), ptBinning.GetBinUpEdge(ptBinning.GetNbins()));
   TH1F* hFastTrackerQA = new TH1F("hFastTrackerQA", ";#it{p}_{T} (GeV/c;Tracking code", ptBinning.GetNbins(), ptBinning.GetBinLowEdge(1), ptBinning.GetBinUpEdge(ptBinning.GetNbins()));
 
   TLorentzVector tlv;
@@ -120,11 +122,12 @@ void drawFastTracker(float magneticField = 5.f, // in units of kGauss
       continue;
     }
 
-    for (int trial = 0; trial < 200; trial++) {
-      hEfficiency->Fill(fastTracker.FastTrack(trkIn, trkOut, nch) > 0, pt);
+    const int nTrials = 200;
+    for (int trial = 0; trial < nTrials; trial++) {
+      hEfficiency->Fill(fastTracker.fastTrack(trkIn, trkOut, nch) > 0, pt);
     }
     int status = 4;
-    status = fastTracker.FastTrack(trkIn, trkOut, nch);
+    status = fastTracker.fastTrack(trkIn, trkOut, nch);
     hFastTrackerQA->Fill(pt, status);
     if (status < 0) {
       LOG(debug) << " --- fatSolve: FastTrack failed with status " << status << " --- ";
@@ -134,12 +137,12 @@ void drawFastTracker(float magneticField = 5.f, // in units of kGauss
     }
     // define the efficiency
     float eff = 1.;
-    for (size_t l = 1; l < fastTracker.GetNLayers(); ++l) {
-      if (fastTracker.IsLayerInert(l)) {
+    for (size_t l = 1; l < fastTracker.getNLayers(); ++l) {
+      if (fastTracker.isLayerInert(l)) {
         continue; // skip inert layers
       }
       float igoodhit = 0.f;
-      igoodhit = fastTracker.GetGoodHitProb(l);
+      igoodhit = fastTracker.getGoodHitProb(l);
       if (igoodhit <= 0.) {
         continue;
       }

ALICE3/Macros/testFastTracker.C

@@ -14,20 +14,34 @@
 /// \brief  Test the FastTracker functionality
 
 #include "ALICE3/Core/FastTracker.h"
+#include "ALICE3/Core/GeometryContainer.h"
 
+#include <CCDB/BasicCCDBManager.h>
 #include <Framework/Logger.h>
 
 #include <string>
 
-void testFastTracker(const std::string& geometryFile = "a3geo.ini")
+void testFastTracker(const std::string& geometryFile = "Configuration/a3geo.ini")
 {
 
   fair::Logger::SetConsoleSeverity(fair::Severity::debug);
 
-  // auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
-  // ccdb.setURL("http://alice-ccdb.cern.ch");
-  o2::fastsim::FastTracker fastTracker;
-  // fastTracker.AddGenericDetector(geometryFile); // FIXME
-  // fastTracker.AddGenericDetector(geometryFile, &ccdb);
-  fastTracker.Print();
+  auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
+  ccdb.setURL("http://alice-ccdb.cern.ch");
+
+  o2::fastsim::GeometryContainer geometryContainer;
+  geometryContainer.setCcdbManager(&ccdb);
+  geometryContainer.addEntry(geometryFile);
+  const int nGeometries = geometryContainer.getNumberOfConfigurations();
+  const float magneticField = geometryContainer.getFloatValue(0, "global", "magneticfield");
+
+  for (int icfg = 0; icfg < nGeometries; ++icfg) {
+    const o2::fastsim::GeometryEntry& geometry = geometryContainer.getEntry(icfg);
+    o2::fastsim::FastTracker fastTracker;
+
+    fastTracker.setMagneticField(magneticField);
+    fastTracker.addGenericDetector(geometry, &ccdb);
+
+    fastTracker.print();
+  }
 }

Common/Core/TPCVDriftManager.h

@@ -37,10 +37,18 @@ class TPCVDriftManager
     mCCDB = ccdb;
   }
 
+  // Use the TPC side flags (with legacy-data fallback) instead of the tgl-sign-based correction.
+  // Off by default so that existing analyses see no change in results until this is explicitly
+  // enabled for testing. Tasks using TPCVDriftManager can expose this via their own Configurable.
+  void setUseSideBasedCorrection(bool value) noexcept
+  {
+    mUseSideBasedCorrection = value;
+  }
+
   void update(uint64_t timestamp) noexcept
   {
     // Keep the object we already have if it is still valid for this timestamp.
-    // firstTime/lastTime are the first and last timestamps of the TFs the correction
+    // firstTime/lastTime are the first and last timestamps of the TF the correction
     // was derived from, so the validity range is closed on both ends.
     if (mVD != nullptr && timestamp >= static_cast<uint64_t>(mVD->firstTime) && timestamp <= static_cast<uint64_t>(mVD->lastTime)) {
       return;
@@ -136,7 +144,7 @@ class TPCVDriftManager
     float dTime = tTB - trackExtra.trackTime();
     float dDrift = dTime * mTPCVDriftNS;
     float dDriftErr = tTBErr * mTPCVDriftNS;
-    if (dDriftErr < 0.f || dDrift > 250.f) { // we cannot move a track outside the drift volume
+    if (dDriftErr < 0.f || dDrift > mMaxDriftCm) { // we cannot move a track outside the drift volume
       if (mOutside < mWarningLimit) {
         LOGP(warn, "Skipping correction outside of tpc volume with dDrift={} +- {}", dDrift, dDriftErr);
         const auto trackBC = trackExtra.template collision_as<Collisions>().template foundBC_as<BCs>().globalBC();
@@ -152,7 +160,54 @@ class TPCVDriftManager
     }
 
     // impose new Z coordinate
-    track.setZ(track.getZ() + ((track.getTgl() < 0.) ? -dDrift : dDrift));
+    float zShift = 0.f;
+    if (!mUseSideBasedCorrection) {
+      // Legacy behaviour (default): infer the side from tgl alone.
+      zShift = (track.getTgl() < 0.f) ? -dDrift : dDrift;
+    } else {
+      const auto sides = (trackExtra.flags() & (o2::aod::track::TrackFlags::TPCSideA | o2::aod::track::TrackFlags::TPCSideC));
+      if (sides == o2::aod::track::TrackFlags::TPCSideA) {
+        zShift = dDrift;
+      } else if (sides == o2::aod::track::TrackFlags::TPCSideC) {
+        zShift = -dDrift;
+      } else if (sides == 0) {
+        // Fallback for datasets produced before the TPC side flags were introduced (Feb. 2026).
+        o2::aod::track::extensions::TPCTimeErrEncoding tEnc;
+        tEnc.encoding.timeErr = trackExtra.trackTimeRes();
+        const float dFwd = tEnc.getDeltaTFwd();
+        const float dBwd = tEnc.getDeltaTBwd();
+        // Equal, small forward/backward margins mean the track is bounded on both ends,
+        // i.e. it crosses the CE: it cannot be moved and is already corrected elsewhere.
+        const bool crossesCE = (dFwd == dBwd) && (dFwd < mMaxCECrossingDeltaTNS);
+        if (!crossesCE) {
+          const bool zPositive = track.getZ() > 0.f;
+          const bool tglPositive = track.getTgl() > 0.f;
+          int side = 0; // +1 = A, -1 = C, 0 = undetermined -> leave uncorrected
+          if (zPositive == tglPositive) {
+            // Consistent sign: the track converges to Z=0 at the beamline by construction.
+            side = tglPositive ? 1 : -1;
+          } else if (dBwd == 0.f && dFwd > 0.f) {
+            // Bounded backward at the CE with room forward: no clusters on the opposite
+            // side, so trust the measured Z rather than the (here inverted) tgl.
+            side = zPositive ? 1 : -1;
+          } else if (dBwd > 0.f) {
+            // Large tgl track bounded at the readout side instead: trust tgl.
+            side = tglPositive ? 1 : -1;
+          }
+          // else: degenerate case, track touches both CE and readout -> cannot be deduced/moved.
+          // (in practice unreachable here: dFwd==dBwd==0 would already satisfy crossesCE above,
+          // since dFwd/dBwd are always >= 0; kept explicit to mirror the reference logic 1:1.)
+
+          if (side > 0) {
+            zShift = dDrift;
+          } else if (side < 0) {
+            zShift = -dDrift;
+          }
+        }
+      }
+      // else: track has clusters on both sides (crossed the CE) and is already corrected elsewhere
+    }
+    track.setZ(track.getZ() + zShift);
     if constexpr (std::is_base_of_v<o2::track::TrackParCov, Track>) {
       track.setCov(track.getSigmaZ2() + dDriftErr * dDriftErr, o2::track::kSigZ2);
     }
@@ -169,6 +224,7 @@ class TPCVDriftManager
 
  private:
   bool mValid{false};
+  bool mUseSideBasedCorrection{false}; // off by default: preserves legacy tgl-sign behaviour
   // Factors
   float mTPCVDriftNS{0.f}; // drift velocity in cm/ns
 
@@ -177,6 +233,8 @@ class TPCVDriftManager
   o2::ccdb::BasicCCDBManager* mCCDB{};  // reference to initialized ccdb manager
 
   static constexpr unsigned int mWarningLimit{10};
+  static constexpr float mMaxDriftCm{250.f};             // TPC drift volume half-length in cm
+  static constexpr float mMaxCECrossingDeltaTNS{1000.f}; // ~1 us, ballpark forward/backward time margin of a CE-crossing track
 
   // Counters
   unsigned int mCalls{0};       // total number of calls