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

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,20 @@ 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,28 +139,52 @@ 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);
   }
 
   for (auto& layer : mLayers) {
     layer.createLayer(vMID);
   }
+
 }
 
 void Detector::Reset()
 {
   if (!o2::utils::ShmManager::Instance().isOperational()) {
     mHits->clear();
   }
+  mTrackData.mHitStarted = false;
 }
 
 bool Detector::ProcessHits(FairVolume* vol)
@@ -159,14 +197,15 @@ 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;
   if (fMC->IsTrackEntering()) {
@@ -213,14 +252,14 @@ 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) { 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());

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),
+                   int nstaves,
+                   float zOffset,
+                   int nModulesZ,
+                   float staveWidth,
+                   int nBars) : mName(layerName),
                                   mRadius(rInn),
                                   mLength(length),
+                                  mZOffset(zOffset),
+                                  mStaveWidth(staveWidth),
                                   mNumber(layerNumber),
-                                  mNStaves(nstaves)
+                                  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,7 +64,8 @@ MIDLayer::Stave::Stave(std::string staveName,
                        float staveLength,
                        float staveWidth,
                        float staveThickness,
-                       int nModulesZ) : mName(staveName),
+                       int nModulesZ,
+                       int nBars) : mName(staveName),
                                         mRadDistance(radDistance),
                                         mRotAngle(rotAngle),
                                         mLength(staveLength),
@@ -64,17 +75,20 @@ MIDLayer::Stave::Stave(std::string staveName,
                                         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);

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

@@ -26,12 +26,15 @@ namespace passive
 
 enum MagnetLayout : int {
   AluminiumStabilizer = 0,
-  CopperStabilizer = 1
+  CopperStabilizer = 1,
+  WindingPack = 2,
+  SuperconductingMagnet = 3
 };
 
 enum DetLayout : int {
   StandardRadius = 0,
-  ReducedRadius = 1
+  ReducedRadius = 1,
+  SteppedAbsorber = 2
 };
 
 struct Alice3PassiveBaseParam : public o2::conf::ConfigurableParamHelper<Alice3PassiveBaseParam> {

Detectors/Upgrades/ALICE3/Passive/src/Absorber.cxx

@@ -130,8 +130,9 @@ void Alice3Absorber::ConstructGeometry()
     LOG(fatal) << "Could not find the barrel volume while constructing absorber geometry";
   }
 
-  TGeoPcon* absorings = new TGeoPcon(0., 360., 18);
   auto& passiveBaseParam = Alice3PassiveBaseParam::Instance();
+  int nSections = (passiveBaseParam.mDetLayout == o2::passive::DetLayout::SteppedAbsorber) ? 6 : 18;
+  TGeoPcon* absorings = new TGeoPcon(0., 360., nSections);
   switch (passiveBaseParam.mDetLayout) {
     case o2::passive::DetLayout::StandardRadius:
       absorings->DefineSection(0, 500, 236, 274);
@@ -173,6 +174,16 @@ void Alice3Absorber::ConstructGeometry()
       absorings->DefineSection(16, -400, 201, 239);
       absorings->DefineSection(17, -500, 201, 239);
       break;
+    case o2::passive::DetLayout::SteppedAbsorber:
+      // Geometria 6 (Ian/tesis): Rext=290 constante, escalon en Rmin.
+      // Externas 45 cm (Rmin=245), central 70 cm (Rmin=220). Ref: Ian DetectorConstruction.cc abs_thickness={45,70,45}
+      absorings->DefineSection(0, -500, 245, 290);
+      absorings->DefineSection(1, -300, 245, 290);
+      absorings->DefineSection(2, -300, 220, 290);
+      absorings->DefineSection(3,  300, 220, 290);
+      absorings->DefineSection(4,  300, 245, 290);
+      absorings->DefineSection(5,  500, 245, 290);
+      break;
     default:
       LOG(fatal) << "Unknown detector layout " << passiveBaseParam.mDetLayout;
       break;