[ALICE3] Add inner- and outermost hit information to alice 3 tracks (#18260)

Author: jesgum

ALICE3/Core/DetLayer.h

@@ -108,25 +108,25 @@ struct DetLayer {
   TString mName;
 
   // position variables
-  float mR; // radius in centimeters
-  float mZ; // mZ dimension in centimeters
+  float mR{}; // radius in centimeters
+  float mZ{}; // mZ dimension in centimeters
 
   // material variables
-  float mX0;   // radiation length
-  float mXrho; // density
+  float mX0{};   // radiation length
+  float mXrho{}; // density
 
   // resolution variables for active layers
-  float mResRPhi; // RPhi resolution in centimeters
-  float mResZ;    // Z resolution in centimeters
+  float mResRPhi{}; // RPhi resolution in centimeters
+  float mResZ{};    // Z resolution in centimeters
 
   // efficiency
-  float mEff; // detection efficiency
+  float mEff{1.f}; // detection efficiency
 
   // dead regions in phi (in radians)
   TGraph* mDeadPhiRegions = nullptr;
 
   // layer type
-  int mType; // 0: undefined/inert, 1: silicon, 2: gas/tpc
+  int mType{0}; // 0: undefined/inert, 1: silicon, 2: gas/tpc
 };
 
 } // namespace o2::fastsim

ALICE3/Core/FastTracker.cxx

@@ -312,6 +312,8 @@ int FastTracker::fastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
   nIntercepts = 0;
   nSiliconPoints = 0;
   nGasPoints = 0;
+  trackIU = -1.f;
+  trackOU = -1.f;
   std::array<float, 3> posIni{}; // provision for != PV
   inputTrack.getXYZGlo(posIni);
   const float initialRadius = std::hypot(posIni[0], posIni[1]);
@@ -416,6 +418,9 @@ int FastTracker::fastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
       firstLayerReached = il;
     }
     lastLayerReached = il;
+    if (layers[il].isTOF()) {
+      continue; // TOF layers are not used for tracking, but are kept in range for the innermost/outermost update
+    }
     nIntercepts++;
   }
 
@@ -475,7 +480,7 @@ int FastTracker::fastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
       continue;
     }
 
-    if (!layers[il].isInert()) { // only update covm for tracker hits
+    if (!layers[il].isInert() && !layers[il].isTOF()) { // only update covm for tracker hits
       const o2::track::TrackParametrization<float>::dim2_t hitpoint = {xyz1[1], xyz1[2]};
       const o2::track::TrackParametrization<float>::dim3_t hitpointcov = {layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi(), 0.f, layers[il].getResolutionZ() * layers[il].getResolutionZ()};
 
@@ -510,7 +515,14 @@ int FastTracker::fastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
     }
 
     hits.push_back(thisHit);
-    if (!layers[il].isInert()) { // good hit probability calculation
+    if (!layers[il].isInert()) {
+      trackIU = std::hypot(spacePoint[0], spacePoint[1]);
+      if (trackOU < 0.f) {
+        trackOU = trackIU;
+      }
+    }
+
+    if (!layers[il].isInert() && !layers[il].isTOF()) { // good hit probability calculation
       float sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi());
       float sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ());
       goodHitProbability[il] = probGoodChiSqHit(layers[il].getRadius() * 100, sigYCmb * 100, sigZCmb * 100);

ALICE3/Core/FastTracker.h

@@ -58,7 +58,7 @@ class FastTracker
   /// \param phiEnd End angle of the dead region (in radians)
   void addDeadPhiRegionInLayer(const std::string& layerName, float phiStart, float phiEnd);
   [[nodiscard]] DetLayer getLayer(const int layer) const { return layers[layer]; }
-  [[nodiscard]] std::vector<DetLayer> getLayers() const { return layers; }
+  [[nodiscard]] const std::vector<DetLayer>& getLayers() const { return layers; }
   [[nodiscard]] int getLayerIndex(const std::string& name) const;
   [[nodiscard]] size_t getNLayers() const { return layers.size(); }
   [[nodiscard]] bool isLayerInert(const int layer) const { return layers[layer].isInert(); }
@@ -127,6 +127,8 @@ class FastTracker
   [[nodiscard]] int getNIntercepts() const { return nIntercepts; }
   [[nodiscard]] int getNSiliconPoints() const { return nSiliconPoints; }
   [[nodiscard]] int getNGasPoints() const { return nGasPoints; }
+  [[nodiscard]] float getTrackIU() const { return trackIU; }
+  [[nodiscard]] float getTrackOU() const { return trackOU; }
   [[nodiscard]] float getGoodHitProb(int layer) const
   {
     return (layer >= 0 && static_cast<size_t>(layer) < goodHitProbability.size()) ? goodHitProbability[layer] : 0.0f;
@@ -172,6 +174,8 @@ class FastTracker
   int nIntercepts = 0;    /// found in first outward propagation
   int nSiliconPoints = 0; /// silicon-based space points added to track
   int nGasPoints = 0;     /// tpc-based space points added to track
+  float trackIU = -1.f;   /// innermost update: radius of the innermost active-layer hit added to track
+  float trackOU = -1.f;   /// outermost update: radius of the outermost active-layer hit added to track
   std::vector<float> goodHitProbability;
 
   ClassDef(FastTracker, 1);

ALICE3/Core/GeometryContainer.cxx

@@ -19,7 +19,6 @@
 
 #include "Common/Core/TableHelper.h"
 
-#include <CCDB/BasicCCDBManager.h>
 #include <Framework/InitContext.h>
 #include <Framework/Logger.h>
 
@@ -41,15 +40,15 @@ void GeometryContainer::init(o2::framework::InitContext& initContext)
   }
   LOG(info) << "Size of detector configuration: " << detectorConfiguration.size();
 
-  bool cleanLutWhenLoaded;
-  const bool foundCleanLutWhenLoaded = common::core::getTaskOptionValue(initContext, "on-the-fly-detector-geometry-provider", "cleanLutWhenLoaded", cleanLutWhenLoaded, false);
+  bool cleanLut;
+  const bool foundCleanLutWhenLoaded = common::core::getTaskOptionValue(initContext, "on-the-fly-detector-geometry-provider", "cleanLutWhenLoaded", cleanLut, false);
   if (!foundCleanLutWhenLoaded) {
     LOG(fatal) << "Could not retrieve foundCleanLutWhenLoaded option from OnTheFlyDetectorGeometryProvider task.";
     return;
   }
-  setLutCleanupSetting(cleanLutWhenLoaded);
+  setLutCleanupSetting(cleanLut);
 
-  for (std::string& configFile : detectorConfiguration) {
+  for (const std::string& configFile : detectorConfiguration) {
     LOG(info) << "Detector geometry configuration file used: " << configFile;
     addEntry(configFile);
   }

ALICE3/Core/GeometryContainer.h

@@ -36,18 +36,18 @@ struct GeometryEntry {
   // Default constructor
   GeometryEntry() = default;
   explicit GeometryEntry(const std::string& filename, o2::ccdb::BasicCCDBManager* ccdb = nullptr)
+    : mFileName(ConfigurationParser::accessFile(filename, "./.ALICE3/Configuration/", ccdb))
   {
-    mFileName = ConfigurationParser::accessFile(filename, "./.ALICE3/Configuration/", ccdb);
     mConfigurations = ConfigurationParser::parseTEnvConfiguration(mFileName, mLayerNames);
     LOG(info) << "Loaded geometry configuration from file: " << mFileName << " with " << mLayerNames.size() << " layers.";
     if (mLayerNames.empty()) {
       LOG(warning) << "No layers found in geometry configuration file: " << filename;
     }
   }
 
-  std::map<std::string, std::map<std::string, std::string>> getConfigurations() const { return mConfigurations; }
+  const std::map<std::string, std::map<std::string, std::string>>& getConfigurations() const { return mConfigurations; }
   std::map<std::string, std::string> getConfiguration(const std::string& layerName) const;
-  std::vector<std::string> getLayerNames() const { return mLayerNames; }
+  const std::vector<std::string>& getLayerNames() const { return mLayerNames; }
   bool hasValue(const std::string& layerName, const std::string& key) const;
   std::string getValue(const std::string& layerName, const std::string& key, bool require = true) const;
   void setValue(const std::string& layerName, const std::string& key, const std::string& value) { mConfigurations[layerName][key] = value; }

ALICE3/DataModel/tracksAlice3.h

@@ -32,9 +32,13 @@ DECLARE_SOA_COLUMN(NTPCHits, nTPCHits, int);                //! number of tpc hi
 DECLARE_SOA_COLUMN(PdgCode, pdgCode, int);                  //! PDG code of the linked truth MC particle
 DECLARE_SOA_COLUMN(TrackType, trackType, int);              //! Type of the track
 DECLARE_SOA_COLUMN(IsPVContributor, isPVContributor, bool); //! Has track contributed to the PV fit?
+DECLARE_SOA_COLUMN(TrackIU, trackIU, float);                //! Innermost update: radius (cm) of the innermost hit assigned to the track
+DECLARE_SOA_COLUMN(TrackOU, trackOU, float);                //! Outermost update: radius (cm) of the outermost hit assigned to the track
 } // namespace track_alice3
 DECLARE_SOA_TABLE(TracksAlice3, "AOD", "TRACKSALICE3",
-                  track_alice3::IsReconstructed);
+                  track_alice3::IsReconstructed,
+                  track_alice3::TrackIU,
+                  track_alice3::TrackOU);
 using TrackAlice3 = TracksAlice3::iterator;
 
 DECLARE_SOA_TABLE(TracksAlice3Pdg, "AOD", "TRACKSALICE3PDG",

ALICE3/Macros/Configuration/a3geometry_v3.ini

@@ -92,15 +92,6 @@ B05.resZ: 0.001
 B05.eff: 1.00
 B05.type: 1
 
-iTOF.r: 19
-iTOF.z: 250
-iTOF.x0: 0.01*3.
-iTOF.xrho: 0.03
-iTOF.resRPhi: 0.001
-iTOF.resZ: 0.001
-iTOF.eff: 1.00
-iTOF.type: 0
-
 B06.r: 20.
 B06.z: 250
 B06.x0: 0.01
@@ -110,6 +101,15 @@ B06.resZ: 0.001
 B06.eff: 1.00
 B06.type: 1
 
+iTOF.r: 21
+iTOF.z: 250
+iTOF.x0: 0.01*3.
+iTOF.xrho: 0.03
+iTOF.resRPhi: 0.001
+iTOF.resZ: 0.001
+iTOF.eff: 1.00
+iTOF.type: 3
+
 B07.r: 30.
 B07.z: 250
 B07.x0: 0.01
@@ -146,6 +146,15 @@ B10.resZ: 0.001
 B10.eff: 1.00
 B10.type: 1
 
+oTOF.r: 92
+oTOF.z: 250
+oTOF.x0: 0.01*3.
+oTOF.xrho: 0.03
+oTOF.resRPhi: 0.001
+oTOF.resZ: 0.001
+oTOF.eff: 1.00
+oTOF.type: 3
+
 # Lookup tables
 global.lutEl: ccdb:/Users/j/jekarlss/LUTs/v3/b/el
 global.lutMu: ccdb:/Users/j/jekarlss/LUTs/v3/b/mu

ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

@@ -299,6 +299,8 @@ struct OnTheFlyTracker {
     int nSiliconHits = 0;         ///< number of silicon hits
     int nTPCHits = 0;             ///< number of TPC hits
     TrackType trackType = TrackType::kNone;
+    float trackIU = -1.f; ///< radius of the innermost hit (-1 if not available)
+    float trackOU = -1.f; ///< radius of the outermost hit (-1 if not available)
   };
 
   // Helper struct to pass cascade information
@@ -1828,7 +1830,7 @@ struct OnTheFlyTracker {
         tableTrackSelection(static_cast<uint8_t>(0), false, false, false, false, false, false);
         tableTrackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false);
       }
-      tableTracksAlice3(true);
+      tableTracksAlice3(true, trackParCov.trackIU, trackParCov.trackOU);
     }
   }
 
@@ -1949,12 +1951,14 @@ struct OnTheFlyTracker {
 
       bool reconstructed = true;
       int nTrkHits = 0;
+      std::array<float, 2> hitRadii = {-1.f, -1.f}; // {innermost, outermost}
       if (enablePrimarySmearing) {
         if (fastPrimaryTrackerSettings.fastTrackPrimaries && longLivedToBeHandled) {
           o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB);
           perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
           computeBremsstrahlungLoss(mcParticle, perfectTrackParCov);
           nTrkHits = fastTracker[icfg]->fastTrack(perfectTrackParCov, trackParCov, dNdEta);
+          hitRadii = {fastTracker[icfg]->getTrackIU(), fastTracker[icfg]->getTrackOU()};
           if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
             reconstructed = false;
           }
@@ -1965,6 +1969,7 @@ struct OnTheFlyTracker {
           const std::array<float, 3> decayVtx = decayer.generateDecayVertex(mcParticle, pdgDB);
           const float decayRadius2D = std::hypot(decayVtx[0], decayVtx[1]);
           nTrkHits = fastTracker[icfg]->fastTrack(perfectTrackParCov, trackParCov, dNdEta, decayRadius2D);
+          hitRadii = {fastTracker[icfg]->getTrackIU(), fastTracker[icfg]->getTrackOU()};
           if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
             reconstructed = false;
           }
@@ -2033,8 +2038,12 @@ struct OnTheFlyTracker {
       bool isPVContributor = trackType == TrackType::kRecoPrimary;
       if (reconstructed) {
         recoPrimaries.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), trackTime, TimeResolutionUs, isDecayDaughter, isPVContributor, false, 0, nTrkHits, trackType});
+        recoPrimaries.back().trackIU = hitRadii[0];
+        recoPrimaries.back().trackOU = hitRadii[1];
       } else {
         ghostPrimaries.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), trackTime, TimeResolutionUs, isDecayDaughter, isPVContributor, false, 0, nTrkHits, trackType});
+        ghostPrimaries.back().trackIU = hitRadii[0];
+        ghostPrimaries.back().trackOU = hitRadii[1];
       }
     }
 
@@ -2176,6 +2185,7 @@ struct OnTheFlyTracker {
 
       bool reconstructed = false;
       int nTrkHits = 0;
+      std::array<float, 2> hitRadii = {-1.f, -1.f}; // {innermost, outermost}
       const bool isSecondary = !otfParticle.isPrimary() && otfParticle.checkBit(o2::upgrade::DecayerBits::ProducedByDecayer) && otfParticle.isAlive();
       if (enablePrimarySmearing && longLivedToBeHandled && otfParticle.isPrimary()) {
         o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
@@ -2186,6 +2196,7 @@ struct OnTheFlyTracker {
         perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
         computeBremsstrahlungLoss(mcParticle, perfectTrackParCov);
         nTrkHits = fastTracker[icfg]->fastTrack(perfectTrackParCov, trackParCov, dNdEta, mcParticle.decayRadius());
+        hitRadii = {fastTracker[icfg]->getTrackIU(), fastTracker[icfg]->getTrackOU()};
         getHist<TH2>(histPath + "h2dGenShortLivedParticleRadius")->Fill(mcParticle.decayRadius(), perfectTrackParCov.getPt());
         getHist<TH2>(histPath + "h2dGenRadiusIniVsDecay")->Fill(std::hypot(perfectTrackParCov.getX(), perfectTrackParCov.getY()), mcParticle.decayRadius());
         if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
@@ -2203,6 +2214,7 @@ struct OnTheFlyTracker {
         computeBremsstrahlungLoss(mcParticle, perfectTrackParCov);
         perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
         nTrkHits = fastTracker[icfg]->fastTrack(perfectTrackParCov, trackParCov, dNdEta);
+        hitRadii = {fastTracker[icfg]->getTrackIU(), fastTracker[icfg]->getTrackOU()};
         reconstructed = nTrkHits >= fastTrackerSettings.minSiliconHits;
       }
 
@@ -2260,9 +2272,13 @@ struct OnTheFlyTracker {
 
       if (reconstructed) {
         tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, TimeResolutionUs, isSecondary, true, false, 0, nTrkHits, kRecoPrimary});
+        tracksAlice3.back().trackIU = hitRadii[0];
+        tracksAlice3.back().trackOU = hitRadii[1];
         getHist<TH1>(histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt());
       } else {
         ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, TimeResolutionUs, isSecondary, false, false, false, 0, nTrkHits, kGhostPrimary});
+        ghostTracksAlice3.back().trackIU = hitRadii[0];
+        ghostTracksAlice3.back().trackOU = hitRadii[1];
       }
     }
 

ALICE3/TableProducer/alice3TrackingTranslator.cxx

@@ -419,7 +419,7 @@ struct Alice3TrackingTranslator {
       // Convert tracks from ACTS to ALICE format
       const size_t nParticlesSim = fileParticlesSim.m_vx->size();
       const size_t nTracks = fileTracksummary.m_ePX_fit->size();
-      std::vector<size_t> idMCparticles;
+      // std::vector<size_t> idMCparticles; // fixme: unused
 
       // local index k within this event -> global AO2D index
       int firstIdxThisEvent = tableStoredMcParticles.lastIndex() + 1;
@@ -438,7 +438,7 @@ struct Alice3TrackingTranslator {
           fileParticlesSim.m_particle->at(iPart),
           fileParticlesSim.m_generation->at(iPart),
           fileParticlesSim.m_sub_particle->at(iPart));
-        barcodeToGlobalIdx[key] = firstIdxThisEvent + (int)iPart;
+        barcodeToGlobalIdx[key] = firstIdxThisEvent + static_cast<int>(iPart);
 
         if (iPart == 0) {
           tableMcCollisions(0,                                // mccollision::BCId,
@@ -486,7 +486,7 @@ struct Alice3TrackingTranslator {
         }
       }
       for (size_t iPart = 0; iPart < nParticlesSim; ++iPart) {
-        int globalIdx = firstIdxThisEvent + (int)iPart;
+        int globalIdx = firstIdxThisEvent + static_cast<int>(iPart);
         int motherIdx = -1;
         if (GlobalIdxToMotherIdx.find(globalIdx) != GlobalIdxToMotherIdx.end()) {
           motherIdx = GlobalIdxToMotherIdx[globalIdx];
@@ -666,7 +666,9 @@ struct Alice3TrackingTranslator {
         tableTracksDCACov(0.0f,  // sigmaDcaXY2
                           0.0f); // sigmaDcaZ2
         // Fill ALICE3 specific tables
-        tableTracksAlice3(true); // isReconstructed
+        tableTracksAlice3(true,  // isReconstructed
+                          -1.f,  // todo: track innermost update
+                          -1.f); // todo: track outermost update
         if (mcParticleIdx > 0)
           tableTracksAlice3Pdg(GlobalIdxToPDGCode[mcParticleIdx]); // PdgCode to the linked MC truth particle
         else