[ALICE3] MI3: Add superconducting magnet/cryostat geometry to ALICE3 MID simulation (#15585)

* [MID] Add superconducting magnet/cryostat geometry to simulation

* Apply clang-format

* [MI3] Add superconducting magnet/cryostat geometry to ALICE3 MID simulation

* [ALICE3] Magnet: add WindingPack (NbTi+Cu+Al) coil material option

* [ALICE3] Magnet: improve WindingPack material documentation

* ICN-UNAM: IanMagnet, IanAbsorber, ICNStepped MID layout — 0 overlaps verified

* ICN-UNAM: IanAbsorber — Geometria 6 (Antonio): cara externa plana Rmax=290, escalon hacia adentro

* ICN-UNAM: cleanup diagnostic instrumentation and translate comments to English

* ICN-UNAM: fix hit production — simple sensor names for VMC resolution and reset mHitStarted on Reset()

* ICN-UNAM: fix physLay in ProcessHits — extract layer index from sensor name; remove ITS TrackReference boilerplate; guard against invalid physLay

* ICN-UNAM: document CurrentVolOffID limitation; full sensor location available via sensor name pattern

* ICN-UNAM: IanAbsorber — correct external Rmin to 245 cm (45 cm thickness per Ian reference code)

* ICN-UNAM: rename enum values — IanMagnet->SuperconductingMagnet, IanAbsorber->SteppedAbsorber, ICNStepped->SteppedLayout

* ICN-UNAM: fix layer 1 half-length to 525 cm (Ian/SD Table 16); revert absorber to Geometry 6 (Rext=290 const, step in Rmin)

* ICN-UNAM: fix ALICE3Field.C — B1 in numerator, Rc=170 cm (R_out_coil per Ian; confirmed by A. Ortiz Aug 2026)

* Refactor field function and improve readability

Refactored field function to use static constexpr for constants and improved readability. Updated condition checks to use std::abs.

* Refactor magnetic field calculation and conditions

Refactor magnetic field calculation by defining variables for Rc, R1, R2, B1, and B2. Update conditions for magnetic field assignment based on radius and beam start position.

* Fix formatting of B1 assignment in ALICE3Field.C

* Refactor member initializer list formatting

* Refactor magnet and detector layout enums

* Refactor absorber geometry initialization logic

* Refactor magnet layout handling in Magnet.cxx

* Refactor constants and improve code structure

Refactor constants and improve readability in Detector.cxx.

---------

Co-authored-by: Nicolò Jacazio <njacazio@users.noreply.github.com>

Author: joseesquivel-hub

Detectors/Upgrades/ALICE3/MID/base/include/MI3Base/MI3BaseParam.h

@@ -26,7 +26,8 @@ namespace mi3
 
 enum MIDLayout : int {
   StandardRadius = 0,
-  ReducedRadius = 1
+  ReducedRadius = 1,
+  SteppedLayout = 2
 };
 
 struct MIDBaseParam : public o2::conf::ConfigurableParamHelper<MIDBaseParam> {

Detectors/Upgrades/ALICE3/MID/simulation/include/MI3Simulation/MIDLayer.h

@@ -92,7 +92,8 @@ class MIDLayer
           float staveLength = 500.f,
           float staveWidth = 50.f,
           float staveThickness = 0.5f,
-          int nModulesZ = 10);
+          int nModulesZ = 10,
+          int nBars = -1);
     void createStave(TGeoVolume* motherVolume);
 
    private:
@@ -110,16 +111,20 @@ class MIDLayer
 
  public:
   MIDLayer() = default;
-  MIDLayer(int layerNumber, std::string layerName, float rInn, float length, int nstaves = 16);
+  MIDLayer(int layerNumber, std::string layerName, float rInn, float length, int nstaves = 16, float zOffset = 0.f, int nModulesZ = 10, float staveWidth = -1.f, int nBars = -1);
   void createLayer(TGeoVolume* motherVolume);
 
  private:
   std::string mName;
   std::vector<Stave> mStaves;
   float mRadius;
   float mLength;
+  float mZOffset;
+  float mStaveWidth;
   int mNumber;
   int mNStaves;
+  int mNModulesZ;
+  int mNBars;
 };
 } // namespace o2::mi3
 

Detectors/Upgrades/ALICE3/MID/simulation/src/Detector.cxx

@@ -20,6 +20,7 @@
 #include "DetectorsBase/Stack.h"
 #include "ITSMFTSimulation/Hit.h"
 #include "MI3Simulation/Detector.h"
+#include <set>
 #include "MI3Base/MI3BaseParam.h"
 
 using o2::itsmft::Hit;
@@ -92,7 +93,19 @@ void Detector::InitializeO2Detector()
 {
   LOG(info) << "Initialize MID O2Detector";
   mGeometryTGeo = GeometryTGeo::Instance();
-  // defineSensitiveVolumes();
+  // Register sensitive volumes
+  TObjArray* allVols = gGeoManager->GetListOfVolumes();
+  TString sensorPattern = GeometryTGeo::getMIDSensorPattern();
+  std::set<TGeoVolume*> registered;
+  for (int i = 0; i < allVols->GetEntries(); i++) {
+    TGeoVolume* v = (TGeoVolume*)allVols->At(i);
+    TString vname = v->GetName();
+    if (vname.Contains(sensorPattern) && registered.find(v) == registered.end()) {
+      AddSensitiveVolume(v);
+      registered.insert(v);
+    }
+  }
+  LOGP(info, "Total MI3 sensitive volumes registered: {}", registered.size());
 }
 
 void Detector::EndOfEvent() { Reset(); }
@@ -125,14 +138,36 @@ void Detector::createGeometry()
   vMID->SetTitle(vstrng);
 
   // Build the MID
-  mLayers.resize(2);
   auto& midParam = MIDBaseParam::Instance();
   const bool standardRadius = (midParam.mLayout == o2::mi3::MIDLayout::StandardRadius);
 
   if (standardRadius) {
+    mLayers.resize(2);
     mLayers[0] = MIDLayer(0, GeometryTGeo::composeSymNameLayer(0), 301.f, 500.f);
-    mLayers[1] = MIDLayer(1, GeometryTGeo::composeSymNameLayer(1), 311.f, 520.f); // arbitrarily reduced to get multiple of 5.2f
+    mLayers[1] = MIDLayer(1, GeometryTGeo::composeSymNameLayer(1), 311.f, 525.f); // 10 modules x 52.5 cm = 525 cm — matches Ian ref. code and SD Table 16 (10.5 m)
+  } else if (midParam.mLayout == o2::mi3::MIDLayout::SteppedLayout) {
+    // Ian Perez Garcia design (ICN-UNAM) — tesis §3.4.7 Geometria 8
+    // 11 cm gap from absorber outer face to MID layer
+    // mLayer index is flat 0-5: even = physical layer 0, odd = physical layer 1
+    // Module step: layer0=99.8cm (2x49.9), layer1=104cm (2x52=2xsumWidth)
+    // Central segment: Rmax_abso=290 -> Layer0=301, Layer1=311, nMod=6, semi-dz=299.4/312 at Z=0
+    // External segments: Rmax_abso=265 -> Layer0=276, Layer1=286, nMod=2, semi-dz=99.8/104 at Z=+-400
+    constexpr float kAbsGap = 11.f;
+    constexpr float kPitch = 10.f;
+    constexpr float kRCen0 = 290.f + kAbsGap; // 301 cm
+    constexpr float kRCen1 = kRCen0 + kPitch; // 311 cm
+    constexpr float kRExt0 = 265.f + kAbsGap; // 276 cm
+    constexpr float kRExt1 = kRExt0 + kPitch; // 286 cm
+    mLayers.resize(6);
+    // length = semi-length = nModulesZ x step (layer0: step=49.9cm, layer1: step=52cm)
+    mLayers[0] = MIDLayer(0, "MIDLayer0_central", kRCen0, 299.4f, 16, 0.f, 6);              // 6 modules x 49.9 cm step
+    mLayers[1] = MIDLayer(1, "MIDLayer1_central", kRCen1, 312.f, 16, 0.f, 6);               // 6 modules x 52 cm step
+    mLayers[2] = MIDLayer(2, "MIDLayer0_forward", kRExt0, 99.8f, 16, +400.f, 2, -1.f, 21);  // 2 modules x 49.9 cm step, nBars=21 for R=276 cm
+    mLayers[3] = MIDLayer(3, "MIDLayer1_forward", kRExt1, 104.f, 16, +405.f, 2);            // 2 modules x 52 cm step, +5 cm offset to clear absorber transition
+    mLayers[4] = MIDLayer(4, "MIDLayer0_backward", kRExt0, 99.8f, 16, -400.f, 2, -1.f, 21); // 2 modules x 49.9 cm step, nBars=21 for R=276 cm
+    mLayers[5] = MIDLayer(5, "MIDLayer1_backward", kRExt1, 104.f, 16, -405.f, 2);           // 2 modules x 52 cm step, -5 cm offset to clear absorber transition
   } else {
+    mLayers.resize(2);
     mLayers[0] = MIDLayer(0, GeometryTGeo::composeSymNameLayer(0), 266.f, 500.f);
     mLayers[1] = MIDLayer(1, GeometryTGeo::composeSymNameLayer(1), 276.f, 520.f);
   }
@@ -147,6 +182,7 @@ void Detector::Reset()
   if (!o2::utils::ShmManager::Instance().isOperational()) {
     mHits->clear();
   }
+  mTrackData.mHitStarted = false;
 }
 
 bool Detector::ProcessHits(FairVolume* vol)
@@ -159,13 +195,16 @@ bool Detector::ProcessHits(FairVolume* vol)
   int lay = vol->getVolumeId();
   int volID = vol->getMCid();
 
-  // Is it needed to keep a track reference when the outer ITS volume is encountered?
+  // TrackReference block removed: ITS boilerplate whose condition (lay == 0
+  // against a TGeo volume ID) never fired. No MID reconstruction consumes
+  // MID track references at present.
   auto stack = (o2::data::Stack*)fMC->GetStack();
-  if (fMC->IsTrackExiting() && (lay == 0)) {
-    o2::TrackReference tr(*fMC, GetDetId());
-    tr.setTrackID(stack->GetCurrentTrackNumber());
-    tr.setUserId(lay);
-    stack->addTrackReference(tr);
+  // Extract physical layer index (0 or 1) from sensor name: MIDSensor_L<lay>_S...
+  int physLay = -1;
+  const char* volName = fMC->CurrentVolName();
+  sscanf(volName, "MIDSensor_L%d", &physLay);
+  if (physLay >= 0) {
+    physLay = physLay % 2;
   }
   bool startHit = false, stopHit = false;
   unsigned char status = 0;
@@ -213,14 +252,17 @@ bool Detector::ProcessHits(FairVolume* vol)
   if (stopHit) {
     TLorentzVector positionStop;
     fMC->TrackPosition(positionStop);
-    // Retrieve the indices with the volume path
-    int stave(0), halfstave(0), chipinmodule(0), module;
-    fMC->CurrentVolOffID(1, chipinmodule);
-    fMC->CurrentVolOffID(2, module);
-    fMC->CurrentVolOffID(3, halfstave);
-    fMC->CurrentVolOffID(4, stave);
-
-    Hit* p = addHit(stack->GetCurrentTrackNumber(), lay, mTrackData.mPositionStart.Vect(), positionStop.Vect(),
+    // CurrentVolOffID(1..4) yields copy numbers of module/halfstave/stave ancestors.
+    // With TGeoVolumeAssembly nodes these are always 0 except the stave level.
+    // Full sensor location (layer, stave, module, bar) is encoded in the sensor
+    // name (MIDSensor_L<l>_S<s>_M<m>_B<b>) and can be decoded with sscanf if needed.
+    // Left as future work for hit digitization.
+
+    if (physLay < 0) {
+      LOGP(warn, "MID sensor name {} did not match expected pattern, cannot extract physical layer index", volName);
+      return false;
+    } // guard: sensor name did not match expected pattern
+    Hit* p = addHit(stack->GetCurrentTrackNumber(), physLay, mTrackData.mPositionStart.Vect(), positionStop.Vect(),
                     mTrackData.mMomentumStart.Vect(), mTrackData.mMomentumStart.E(), positionStop.T(),
                     mTrackData.mEnergyLoss, mTrackData.mTrkStatusStart, status);
     // p->SetTotalEnergy(vmc->Etot());
@@ -241,4 +283,4 @@ o2::itsmft::Hit* Detector::addHit(int trackID, int detID, const TVector3& startP
   return &(mHits->back());
 }
 } // namespace o2::mi3
-ClassImp(o2::mi3::Detector);
+ClassImp(o2::mi3::Detector);

Detectors/Upgrades/ALICE3/MID/simulation/src/MIDLayer.cxx

@@ -26,23 +26,33 @@ MIDLayer::MIDLayer(int layerNumber,
                    std::string layerName,
                    float rInn,
                    float length,
-                   int nstaves) : mName(layerName),
-                                  mRadius(rInn),
-                                  mLength(length),
-                                  mNumber(layerNumber),
-                                  mNStaves(nstaves)
+                   int nstaves,
+                   float zOffset,
+                   int nModulesZ,
+                   float staveWidth,
+                   int nBars) : mName(layerName),
+                                mRadius(rInn),
+                                mLength(length),
+                                mZOffset(zOffset),
+                                mStaveWidth(staveWidth),
+                                mNumber(layerNumber),
+                                mNStaves(nstaves),
+                                mNModulesZ(nModulesZ),
+                                mNBars(nBars)
 {
   mStaves.reserve(nstaves);
-  LOGP(debug, "Constructing MIDLayer: {} with inner radius: {}, length: {} cm and {} staves", mName, mRadius, mLength, mNStaves);
+  LOGP(debug, "Constructing MIDLayer: {} with inner radius: {}, length: {} cm, {} staves and {} modules/stave", mName, mRadius, mLength, mNStaves, mNModulesZ);
   for (int iStave = 0; iStave < mNStaves; ++iStave) {
     mStaves.emplace_back(GeometryTGeo::composeSymNameStave(layerNumber, iStave),
                          mRadius,
                          TMath::TwoPi() / (float)nstaves * iStave,
                          mNumber,
                          iStave,
                          mLength,
-                         !layerNumber ? 59.8f : 61.75f,
-                         0.5f);
+                         !(layerNumber % 2) ? 59.8f : 61.75f,
+                         0.5f,
+                         mNModulesZ,
+                         mNBars);
   }
 }
 
@@ -54,27 +64,31 @@ MIDLayer::Stave::Stave(std::string staveName,
                        float staveLength,
                        float staveWidth,
                        float staveThickness,
-                       int nModulesZ) : mName(staveName),
-                                        mRadDistance(radDistance),
-                                        mRotAngle(rotAngle),
-                                        mLength(staveLength),
-                                        mWidth(staveWidth),
-                                        mThickness(staveThickness),
-                                        mLayer(layer),
-                                        mNumber(number),
-                                        mNModulesZ(nModulesZ)
+                       int nModulesZ,
+                       int nBars) : mName(staveName),
+                                    mRadDistance(radDistance),
+                                    mRotAngle(rotAngle),
+                                    mLength(staveLength),
+                                    mWidth(staveWidth),
+                                    mThickness(staveThickness),
+                                    mLayer(layer),
+                                    mNumber(number),
+                                    mNModulesZ(nModulesZ)
 {
+  // nBars=-1 uses default calibrated for standard radii
+  int effNBars = (nBars < 0) ? (!(mLayer % 2) ? 23 : 20) : nBars;
+  float moduleOffset = -effNBars * 5.2f / 2.f; // 5.2 = 2*barWidth + barSpacing
   // Staves are ideal shapes made of air including the modules, for now.
-  LOGP(debug, "\t\tConstructing MIDStave: {} layer: {} at angle {}", mName, mLayer, mRotAngle * TMath::RadToDeg());
+  LOGP(debug, "\t\tConstructing MIDStave: {} layer: {} at angle {} nBars={}", mName, mLayer, mRotAngle * TMath::RadToDeg(), effNBars);
   mModules.reserve(nModulesZ);
   for (int iModule = 0; iModule < mNModulesZ; ++iModule) {
     mModules.emplace_back(GeometryTGeo::composeSymNameModule(mLayer, mNumber, iModule),
                           mLayer,
                           mNumber,
                           iModule,
-                          !mLayer ? 23 : 20,
+                          effNBars,
                           -staveLength,
-                          !mLayer ? 49.9f : 61.75f);
+                          !(mLayer % 2) ? 49.9f : 61.75f);
   }
 }
 
@@ -106,8 +120,8 @@ MIDLayer::Stave::Module::Module(std::string moduleName,
                           mStave,
                           mNumber,
                           iBar,
-                          !mLayer ? -59.8f : -52.f,  // offset
-                          !mLayer ? 49.9f : 61.75f); // sensor length
+                          -mNBars * 5.2f / 2.f,            // moduleOffset derived from nBars
+                          !(mLayer % 2) ? 49.9f : 61.75f); // sensor length
   }
 }
 
@@ -136,9 +150,13 @@ MIDLayer::Stave::Module::Sensor::Sensor(std::string sensorName,
 
 void MIDLayer::createLayer(TGeoVolume* motherVolume)
 {
-  LOGP(debug, "Creating MIDLayer: {}", mName);
+  LOGP(debug, "Creating MIDLayer: {} at zOffset={} cm", mName, mZOffset);
   TGeoVolumeAssembly* layerVolume = new TGeoVolumeAssembly(mName.c_str());
-  motherVolume->AddNode(layerVolume, 0);
+  if (mZOffset != 0.f) {
+    motherVolume->AddNode(layerVolume, 0, new TGeoTranslation(0, 0, mZOffset));
+  } else {
+    motherVolume->AddNode(layerVolume, 0);
+  }
   for (auto& stave : mStaves) {
     stave.createStave(layerVolume);
   }
@@ -172,7 +190,7 @@ void MIDLayer::Stave::Module::createModule(TGeoVolume* motherVolume)
     sensor.createSensor(moduleVolume);
   }
   TGeoCombiTrans* modTrans = nullptr;
-  if (!mLayer) {
+  if (!(mLayer % 2)) {
     modTrans = new TGeoCombiTrans(0, 0, mZOffset + mNumber * 2 * mBarLength + mBarLength, nullptr);
   } else {
     modTrans = new TGeoCombiTrans(0, 0, mZOffset + mNumber * 2 * sumWidth + sumWidth, nullptr);
@@ -184,17 +202,19 @@ void MIDLayer::Stave::Module::Sensor::createSensor(TGeoVolume* motherVolume)
 {
   LOGP(debug, "\t\t\t\tCreating MIDSensor: {}", mName);
   TGeoBBox* sensor = nullptr;
-  if (!mLayer) {
+  if (!(mLayer % 2)) {
     sensor = new TGeoBBox(mName.c_str(), mWidth, mThickness, mLength);
   } else {
     sensor = new TGeoBBox(mName.c_str(), mLength, mThickness, mWidth);
   }
   auto* polyMed = gGeoManager->GetMedium("MI3_POLYSTYRENE");
-  TGeoVolume* sensorVolume = new TGeoVolume(mName.c_str(), sensor, polyMed);
+  // Simple unique name without slashes so gMC->VolId() resolves correctly during stepping
+  auto volName = Form("MIDSensor_L%d_S%d_M%d_B%d", mLayer, mStave, mNumber, mNumber);
+  TGeoVolume* sensorVolume = new TGeoVolume(volName, sensor, polyMed);
   sensorVolume->SetVisibility(true);
   auto totWidth = mWidth + mSpacing / 2;
   TGeoTranslation* sensorTrans = nullptr;
-  if (!mLayer) {
+  if (!(mLayer % 2)) {
     sensorTrans = new TGeoTranslation(mModuleOffset + 2 * totWidth * mNumber + totWidth, 0, 0);
     sensorVolume->SetLineColor(kAzure + 4);
     sensorVolume->SetTransparency(50);
@@ -205,4 +225,4 @@ void MIDLayer::Stave::Module::Sensor::createSensor(TGeoVolume* motherVolume)
   }
   motherVolume->AddNode(sensorVolume, 0, sensorTrans);
 }
-} // namespace o2::mi3
+} // namespace o2::mi3

Detectors/Upgrades/ALICE3/Passive/include/Alice3DetectorsPassive/PassiveBaseParam.h

@@ -24,21 +24,24 @@ namespace passive
 // ** Parameters for Passive base configuration
 // **
 
-enum MagnetLayout : int {
-  AluminiumStabilizer = 0,
-  CopperStabilizer = 1
+enum MagnetType : int {
+  AluminiumStabilizer = 0,  // Using Aluminium stabilizer for the magnet
+  CopperStabilizer = 1,     // Using Copper stabilizer for the magnet
+  WindingPack = 2,          // Using Winding Pack for the magnet
+  SuperconductingMagnet = 3 // Using Superconducting magnet (NbTi+Cu+Al) for the magnet
 };
 
-enum DetLayout : int {
-  StandardRadius = 0,
-  ReducedRadius = 1
+enum MagnetAbsorberLayout : int {
+  StandardRadius = 0, // Using standard radius for the magnet and absorber
+  ReducedRadius = 1,  // Using reduced radius for the magnet and absorber
+  SteppedAbsorber = 2 // Using stepped absorber for the magnet and absorber
 };
 
 struct Alice3PassiveBaseParam : public o2::conf::ConfigurableParamHelper<Alice3PassiveBaseParam> {
   // Geometry Builder parameters
 
-  int mLayout = MagnetLayout::AluminiumStabilizer;
-  int mDetLayout = DetLayout::StandardRadius;
+  MagnetType mMagType = MagnetType::AluminiumStabilizer;                                   // Magnet type: as in MagnetType enum
+  MagnetAbsorberLayout mMagAbsLayout = o2::passive::MagnetAbsorberLayout::SteppedAbsorber; // Magnet and absorber layout: as in MagnetAbsorberLayout enum
 
   O2ParamDef(Alice3PassiveBaseParam, "Alice3PassiveBase");
 };