create_stsgeo_v24c.C

/* Copyright (C) 2012-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: David Emschermann [committer] */

/******************************************************************************
 ** Creation of STS geometry in ROOT format (TGeo).
 **
 ** @file create_stsgeo_v24b.C
 ** @author David Emschermann <d.emschermann@gsi.de>
 ** @author Volker Friese <v.friese@gsi.de>
 ** @since 15 June 2012
 ** @date 09.05.2014
 ** @author Tomas Balog <T.Balog@gsi.de>
 **
 ** 2024-04-29 - MS - v24c: based on v24b, shift STS in Z to position measured in the cave
 **                         as it was seen that the station 0 module position was not correct 
 **			    and the measurement was done again in the CAVE. The position of
 **                         the Back-side of the STS box was expected to be at 57.8 cm. 
 **			    We opened the side wall of the sts-box and did the measurement
 **                         inside the box for the placement of the C-Frames. 
 **			    The thickness of the back-side wall is 0.5 cm which means from
 **                         inside the box (57.8-0.5=57.3) the last C-Frame edge is at 7 cm
 **                         away (50.3). The edge of the front wall from out-side is at 22.8 cm
 **                         (35 cm away from the backside wall). And the first module box wall
 **                         from the wall is 1.6 cm away.  
 ** 2024-03-27 - DE - v24b: based on v24a, shift STS in z to position measured in the cave
 ** 2024-03-26 - DE - v24a: add preliminary version of mSTS v24a for the 2024_03_22 gold setup
 ** 2024-02-06 - DE - v24a: position the STS sensors as measured in the cave
 ** 2023-06-29 - DE - v23b: add a 6x6 cm2 module at z=16.5 cm, keep STS box around 2x2 and 3x3 (from v22e)
 ** 2023-06-27 - DE - v23a: add a 6x6 cm2 module at z=14.0 cm (nominal), keep STS box around 2x2 and 3x3 (from v22e)
 ** 2023-02-09 - DE - v22f: decrease distance between statios 0 and 1 from 140 mm to 135 mm, to reflect CAD
 ** 2023-02-09 - DE - v22f: increase gap between units of a station gkLadderGapZ from 10 mm to 12 mm, to reflect CAD
 ** 2023-02-09 - DE - v22f: shift ladder U2 downwards in y by 2.2 mm
 ** 2023-02-08 - DE - v22f: introduce ladder type 12 (based on 10) for U2, but with original gkSectorOverlapY of 4.6 mm
 ** 2022-09-13 - DE - v22e: correct gkSectorOverlapY to 9.0 mm for the ladder type 10 (with 12 cm sensors) 
 ** 2022-03-18 - DE - v22d: add position offset for mSTS 2022_03_22_iron as surveyed in the cave
 ** 2022-02-01 - DE - v22b: add position offset for mSTS 07_2021 as surveyed in the cave
 ** 2022-01-25 - DE - v22a: add a aluminium box around v20e
 ** 2020-05-16 - DE - v20e: swap ladders in unit 3 - 3x 6x6 on beam side
 ** 2020-03-21 - DE - v20d: build mSTS for March 2020 - 1 ladder - shift -3 cm in x
 ** 2020-03-11 - DE - v20c: build mSTS for March 2021 - 5 ladders
 ** 2020-03-11 - DE - v20b: build mSTS for May   2020 - 2 ladders
 ** 2020-03-11 - DE - v20a: build mSTS for March 2020 - 1 ladder
 **
 ** 2019-12-04 - DE - v19d: build 2nd mSTS v19d station from one 6x6 and one 6x12 cm x cm sensors
 ** 2019-12-04 - DE - v19c: build 1st mSTS v19c station at nominal position
 ** 2019-08-12 - DE - v19a: mSTS as built in March 2019 - the mSTS FEBs are in the bottom
 ** 2019-03-15 - DE - v18n: mSTS as built in March 2019 - downstream ladder of station 0 at position of station 1
 ** 2018-01-18 - DE - v18g: set overlaps in X and Y according to mSTS CAD model
 **
 ** v18f: flip orientation of carbon ladders for primary beam left of mSTS, change z-positions to 30 and 45 cm
 ** v18e: 2 stations, derived from v15b, 1st 2x2, 2nd of 3x3 sensor array, sensor size 6x6 cm2 with carbon ladder supports
 ** v18d: 2 stations of 3x3 sensor array, sensor size 6x6 cm2 with carbon ladder supports
 ** v18c: (causes segfault due to divide) 2 stations of 3x3 sensor array, sensor size 6x6 cm2 with carbon ladder supports
 ** v18b: 2 stations of 4x4 sensor array, sensor size 6x6 cm2
 ** v18a: 2 stations of 3x3 sensor array, sensor size 6x6 cm2
 **
 ** v15b: introduce modified carbon ladders from v13z
 ** v15a: with flipped ladder orientation for stations 0,2,4,6 to match CAD design
 **
 ** TODO:
 **
 ** DONE:
 ** v15b - use carbon macaroni as ladder support
 ** v15b - introduce a small gap between lowest sensor and carbon ladder
 ** v15b - build small cones for the first 2 stations
 ** v15b - within a station the ladders of adjacent units should not touch eachother - set gkLadderGapZ to 10 mm
 ** v15b - for all ladders set an even number of ladder elements 
 ** v15b - z offset of cones to ladders should not be 0.3 by default, but 0.26
 ** v15b - within a station the ladders should be aligned in z, defined either by the unit or the ladder with most sensors
 ** v15b - get rid of cone overlap in stations 7 and 8 - done by adapting rHole size
 **
 ** The geometry hierarachy is:
 **
 ** 1. Sensors  (see function CreateSensors)
 **    The sensors are the active volumes and the lowest geometry level.
 **    They are built as TGeoVolumes, shape box, material silicon.
 **    x size is determined by strip pitch 58 mu and 1024 strips 
 **    plus guard ring of 1.3 mm at each border -> 6.1992 cm.
 **    Sensor type 1 is half of that (3.0792 cm).
 **    y size is determined by strip length (2.2 / 4.2 / 6.3 cm) plus
 **    guard ring of 1.3 mm at top and bottom -> 2.46 / 4.46 / 6.46 cm.
 **    z size is a parameter, to be set by gkSensorThickness.
 **
 ** 2. Sectors  (see function CreateSectors)
 **    Sectors consist of several chained sensors. These are arranged
 **    vertically on top of each other with a gap to be set by
 **    gkChainGapY. Sectors are constructed as TGeoVolumeAssembly.
 **    The sectors are auxiliary volumes used for proper placement
 **    of the sensor(s) in the module. They do not show up in the
 **    final geometry.
 **
 ** 3. Modules (see function ConstructModule)
 **    A module is a readout unit, consisting of one sensor or
 **    a chain of sensors (see sector) and a cable.
 **    The cable extends from the top of the sector vertically to the
 **    top of the halfladder the module is placed in. The cable and module
 **    volume thus depend on the vertical position of the sector in 
 **    the halfladder. The cables consist of silicon with a thickness to be
 **    set by gkCableThickness.
 **    Modules are constructed as TGeoVolume, shape box, medium gStsMedium.
 **    The module construction can be switched off (gkConstructCables)
 **    to reproduce older geometries.
 **
 ** 4. Halfladders (see function ConstructHalfLadder)
 **    A halfladder is a vertical assembly of several modules. The modules
 **    are placed vertically such that their sectors overlap by 
 **    gkSectorOverlapY. They are displaced in z direction to allow for the 
 **    overlap in y by gkSectorGapZ.
 **    The horizontal placement of modules in the halfladder can be choosen
 **    to left aligned or right aligned, which only matters if sensors of
 **    different x size are involved.
 **    Halfladders are constructed as TGeoVolumeAssembly.
 **
 ** 5. Ladders (see function CreateLadders and ConstructLadder)
 **    A ladder is a vertical assembly of two halfladders, and is such the
 **    vertical building block of a station. The second (bottom) half ladder
 **    is rotated upside down. The vertical arrangement is such that the
 **    inner sectors of the two halfladders have the overlap gkSectorOverlapY
 **    (function CreateLadder) or that there is a vertical gap for the beam
 **    hole (function CreateLadderWithGap).
 **    Ladders are constructed as TGeoVolumeAssembly.
 **   
 ** 6. Stations (see function ConstructStation)
 **    A station represents one layer of the STS geometry: one measurement
 **    at (approximately) a given z position. It consist of several ladders
 **    arranged horizontally to cover the acceptance.
 **    The ladders are arranged such that there is a horizontal overlap
 **    between neighbouring ladders (gkLadderOverLapX) and a vertical gap
 **    to allow for this overlap (gkLadderGapZ). Each second ladder is
 **    rotated around its y axis to face away from or into the beam.
 **    Stations are constructed as TGeoVolumes, shape box minus tube (for
 **    the beam hole), material gStsMedium.
 **
 ** 7. STS
 **    The STS is a volume hosting the entire detectors system. It consists
 **    of several stations located at different z positions.
 **    The STS is constructed as TGeoVolume, shape box minus cone (for the
 **    beam pipe), material gStsMedium. The size of the box is computed to
 **    enclose all stations.
 *****************************************************************************/


// Remark: With the proper steering variables, this should exactly reproduce
// the geometry version v11b of A. Kotynia's described in the ASCII format.
// The only exception is a minimal difference in the z position of the
// sectors/sensors. This is because of ladder types 2 and 4 containing the half
// sensors around the beam hole (stations 1,2 and 3). In v11b, the two ladders
// covering the beam hole cannot be transformed into each other by rotations,
// but only by a reflection. This means they are constructionally different.
// To avoid introducing another two ladder types, the difference in z position
// was accepted.


// Differences to v12:
// gkChainGap reduced from 1 mm to 0
// gkCableThickness increased from 100 mum to 200 mum (2 cables per module)
// gkSectorOverlapY reduced from 3 mm to 2.4 mm
// New sensor types 05 and 06
// New sector types 07 and 08
// Re-definiton of ladders (17 types instead of 8)
// Re-definiton of station from new ladders


#include "TGeoCompositeShape.h"
#include "TGeoCone.h"
#include "TGeoManager.h"
#include "TGeoPara.h"
#include "TGeoTube.h"

#include <iomanip>
#include <iostream>


// -------------   Installation offset      -----------------------------------

// position in x and y agrees within 2 mm to CbmRoot geometry (01.02.2022)
const Double_t gOffX = 0.0;  // Offset from nominal position
const Double_t gOffY = 0.0;  // Offset from nominal position
// 2023
// const Double_t gOffZ = 5.0;  // Offset from nominal position
// 2024
const Double_t gOffZ = 4.55;  // Offset from nominal position

const Bool_t IncludeBox = true; // false;  // true, if STS box is plotted

// -------------   Steering variables       -----------------------------------

// ---> Horizontal width of sensors [cm]
const Double_t gkSensorSizeX = 6.2092;

// ---> Thickness of sensors [cm]
const Double_t gkSensorThickness = 0.03;

// ---> Vertical gap between chained sensors [cm]
const Double_t gkChainGapY = 0.00;

// ---> Thickness of cables [cm]
const Double_t gkCableThickness = 0.02;

// ---> Horizontal overlap of neighbouring ladders [cm]
const Double_t gkLadderOverlapX = 0.25;  // delta X - Oleg CAD 14/05/2020

// ---> Vertical overlap of neighbouring sectors in a ladder [cm]
const Double_t gkSectorOverlapY = 0.46;  // delta Y - Oleg CAD 14/05/2020
const Double_t gkSectorOverlapYL10 = 0.90;  // for ladders with 124 mm sensors - Oleg CAD 13/09/2022

// ---> Gap in z between neighbouring sectors in a ladder [cm]
const Double_t gkSectorGapZ = 0.17;  // gap + thickness = pitch // delta Z pitch = 0.20 - Oleg CAD 14/05/2020

// ---> Gap in z between neighbouring ladders [cm]
const Double_t gkLadderGapZ = 1.20;  // DEJH -> 0.90 / 0.50
// v22e const Double_t gkLadderGapZ = 1.00;  // DEJH -> 0.90 / 0.50

// ---> Gap in z between lowest sector to carbon support structure [cm]
const Double_t gkSectorGapZFrame = 0.10;

// ---> Switch to construct / not to construct readout cables
const Bool_t gkConstructCables = kTRUE;

// ---> Switch to construct / not to construct frames
const Bool_t gkConstructCones       = kFALSE;  // kFALSE;  // switch this false for v15a
const Bool_t gkConstructFrames      = kTRUE;   // kFALSE;  // switch this false for v15a
const Bool_t gkConstructSmallFrames = kTRUE;   // kFALSE;
const Bool_t gkCylindricalFrames    = kTRUE;   // kFALSE;

// ---> Size of the frame
const Double_t gkFrameThickness     = 0.2;
const Double_t gkThinFrameThickness = 0.05;
const Double_t gkFrameStep          = 4.0;  // size of frame cell along y direction

const Double_t gkCylinderDiaInner =
  0.07;  // properties of cylindrical carbon supports, see CBM-STS Integration Meeting (10 Jul 2015)
const Double_t gkCylinderDiaOuter =
  0.15;  // properties of cylindrical carbon supports, see CBM-STS Integration Meeting (10 Jul 2015)

// ----------------------------------------------------------------------------


// --------------   Parameters of beam pipe in the STS region    --------------
// ---> Needed to compute stations and STS such as to avoid overlaps
const Double_t gkPipeZ1 = 22.0;
const Double_t gkPipeR1 = 1.8;
const Double_t gkPipeZ2 = 50.0;
const Double_t gkPipeR2 = 1.8;
const Double_t gkPipeZ3 = 125.0;
const Double_t gkPipeR3 = 5.5;

//DE const Double_t gkPipeZ1 =  27.0;
//DE const Double_t gkPipeR1 =   1.05;
//DE const Double_t gkPipeZ2 = 160.0;
//DE const Double_t gkPipeR2 =   3.25;
// ----------------------------------------------------------------------------


// -------------   Other global variables   -----------------------------------
// ---> STS medium (for every volume except silicon)
TGeoMedium* gStsMedium = NULL;  // will be set later
// ---> TGeoManager (too lazy to write out 'Manager' all the time
TGeoManager* gGeoMan = NULL;  // will be set later
// ----------------------------------------------------------------------------

Int_t CreateSubplates();
Int_t CreatePlates();
Int_t CreateSensors();
Int_t CreateSectors();
Int_t CreateLadders();
void CheckVolume(TGeoVolume* volume);
void CheckVolume(TGeoVolume* volume, fstream& file);
TGeoVolume* ConstructFrameElement(const TString& name, TGeoVolume* frameBoxVol, Double_t x);
TGeoVolume* ConstructSmallCone(Double_t coneDz);
TGeoVolume* ConstructBigCone(Double_t coneDz);
TGeoVolume* ConstructHalfLadder(Int_t ladderid, const TString& name, Int_t nSectors, Int_t* sectorTypes, char align);
TGeoVolume* ConstructLadder(Int_t LadderIndex, TGeoVolume* halfLadderU, TGeoVolume* halfLadderD, Double_t shiftZ);
TGeoVolume* ConstructLadderWithGap(Int_t LadderIndex, TGeoVolume* halfLadderU, TGeoVolume* halfLadderD, Double_t gapY);
TGeoVolume* ConstructStation(Int_t iStation, Int_t nLadders, Int_t* ladderTypes, Double_t rHole);

// ============================================================================
// ======                         Main function                           =====
// ============================================================================

void create_stsgeo_v24c(const char* geoTag = "v24c_mcbm")
{

  // -------   Geometry file name (output)   ----------------------------------
  TString geoFileName = "sts_";
  geoFileName         = geoFileName + geoTag + ".geo.root?reproducible";
  // --------------------------------------------------------------------------


  // -------   Open info file   -----------------------------------------------
  TString infoFileName = geoFileName;
  infoFileName.ReplaceAll("root", "info");
  fstream infoFile;
  infoFile.open(infoFileName.Data(), fstream::out);
  infoFile << "STS geometry created with create_stsgeo_v24c.C" << endl << endl;
  infoFile << "Global variables: " << endl;
  infoFile << "Sensor thickness = " << gkSensorThickness << " cm" << endl;
  infoFile << "Vertical gap in sensor chain = " << gkChainGapY << " cm" << endl;
  infoFile << "Vertical overlap of sensors = " << gkSectorOverlapY << " cm" << endl;
  infoFile << "Gap in z between neighbour sensors = " << gkSectorGapZ << " cm" << endl;
  infoFile << "Horizontal overlap of sensors = " << gkLadderOverlapX << " cm" << endl;
  infoFile << "Gap in z between neighbour ladders = " << gkLadderGapZ << " cm" << endl;
  if (gkConstructCables) infoFile << "Cable thickness = " << gkCableThickness << " cm" << endl;
  else
    infoFile << "No cables" << endl;
  infoFile << endl;
  infoFile << "Beam pipe: R1 = " << gkPipeR1 << " cm at z = " << gkPipeZ1 << " cm" << endl;
  infoFile << "Beam pipe: R2 = " << gkPipeR2 << " cm at z = " << gkPipeZ2 << " cm" << endl;
  infoFile << "Beam pipe: R3 = " << gkPipeR3 << " cm at z = " << gkPipeZ3 << " cm" << endl;
  // --------------------------------------------------------------------------

  // Load FairRunSim to ensure the correct unit system
  FairRunSim* sim = new FairRunSim();

  // -------   Load media from media file   -----------------------------------
  FairGeoLoader* geoLoad    = new FairGeoLoader("TGeo", "FairGeoLoader");
  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  TString geoPath           = gSystem->Getenv("VMCWORKDIR");
  TString medFile           = geoPath + "/geometry/media.geo";
  geoFace->setMediaFile(medFile);
  geoFace->readMedia();
  gGeoMan = gGeoManager;
  // --------------------------------------------------------------------------


  // -----------------   Get and create the required media    -----------------
  FairGeoMedia* geoMedia   = geoFace->getMedia();
  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();

  // ---> air
  FairGeoMedium* mAir = geoMedia->getMedium("air");
  if (!mAir) Fatal("Main", "FairMedium air not found");
  geoBuild->createMedium(mAir);
  TGeoMedium* air = gGeoMan->GetMedium("air");
  if (!air) Fatal("Main", "Medium air not found");

  // ---> silicon
  FairGeoMedium* mSilicon = geoMedia->getMedium("silicon");
  if (!mSilicon) Fatal("Main", "FairMedium silicon not found");
  geoBuild->createMedium(mSilicon);
  TGeoMedium* silicon = gGeoMan->GetMedium("silicon");
  if (!silicon) Fatal("Main", "Medium silicon not found");

  // ---> carbon
  FairGeoMedium* mCarbon = geoMedia->getMedium("carbon");
  if (!mCarbon) Fatal("Main", "FairMedium carbon not found");
  geoBuild->createMedium(mCarbon);
  TGeoMedium* carbon = gGeoMan->GetMedium("carbon");
  if (!carbon) Fatal("Main", "Medium carbon not found");

  // ---> aluminium
  FairGeoMedium* mAluminium = geoMedia->getMedium("aluminium");
  if (!mAluminium) Fatal("Main", "FairMedium aluminium not found");
  geoBuild->createMedium(mAluminium);
  TGeoMedium* aluminium = gGeoMan->GetMedium("aluminium");
  if (!aluminium) Fatal("Main", "Medium aluminium not found");

  // ---> STScable
  FairGeoMedium* mSTScable = geoMedia->getMedium("STScable");
  if (!mSTScable) Fatal("Main", "FairMedium STScable not found");
  geoBuild->createMedium(mSTScable);
  TGeoMedium* STScable = gGeoMan->GetMedium("STScable");
  if (!STScable) Fatal("Main", "Medium STScable not found");

  // ---
  gStsMedium = air;
  // --------------------------------------------------------------------------


  // --------------   Create geometry and top volume  -------------------------
  gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
  gGeoMan->SetName("STSgeom");
  TGeoVolume* top = new TGeoVolumeAssembly("TOP");
  gGeoMan->SetTopVolume(top);
  // --------------------------------------------------------------------------


  // --------------   Create media   ------------------------------------------
  /*
  cout << endl;
  cout << "===> Creating media....";
  cout << CreateMedia();
  cout << " media created" << endl;
  TList* media = gGeoMan->GetListOfMedia();
  for (Int_t iMedium = 0; iMedium < media->GetSize(); iMedium++ ) {
    cout << "Medium " << iMedium << ": " 
	 << ((TGeoMedium*) media->At(iMedium))->GetName() << endl;
  }
  gStsMedium = gGeoMan->GetMedium("air");
  if ( ! gStsMedium ) Fatal("Main", "medium sts_air not found");
  */
  // --------------------------------------------------------------------------

  // ---------------   Create sensors   ---------------------------------------
  cout << endl << endl;
  cout << "===> Creating sensors...." << endl << endl;
  infoFile << endl << "Sensors: " << endl;
  Int_t nSensors = CreateSensors();
  for (Int_t iSensor = 1; iSensor <= nSensors; iSensor++) {
    TString name       = Form("Sensor%02d", iSensor);
    TGeoVolume* sensor = gGeoMan->GetVolume(name);

    // add color to sensors
    if (iSensor == 1) sensor->SetLineColor(kYellow);
    if (iSensor == 2) sensor->SetLineColor(kRed);
    if (iSensor == 3) sensor->SetLineColor(kBlue);
    if (iSensor == 4) sensor->SetLineColor(kAzure + 7);
    if (iSensor == 5) sensor->SetLineColor(kGreen);
    if (iSensor == 6) sensor->SetLineColor(kYellow);

    CheckVolume(sensor);
    CheckVolume(sensor, infoFile);
  }
  // --------------------------------------------------------------------------


  // ----------------   Create sectors   --------------------------------------
  cout << endl << endl;
  cout << "===> Creating sectors...." << endl;
  infoFile << endl << "Sectors: " << endl;
  Int_t nSectors = CreateSectors();
  for (Int_t iSector = 1; iSector <= nSectors; iSector++) {
    cout << endl;
    TString name       = Form("Sector%02d", iSector);
    TGeoVolume* sector = gGeoMan->GetVolume(name);
    CheckVolume(sector);
    CheckVolume(sector, infoFile);
  }
  // --------------------------------------------------------------------------


  // ----------------   Create ladders   --------------------------------------
  TString name = "";
  cout << endl << endl;
  cout << "===> Creating ladders...." << endl;
  infoFile << endl << "Ladders:" << endl;
  Int_t nLadders = CreateLadders();
  for (Int_t iLadder = 1; iLadder <= nLadders; iLadder++) {
    cout << endl;
    name               = Form("Ladder%02d", iLadder);
    TGeoVolume* ladder = gGeoMan->GetVolume(name);
    CheckVolume(ladder);
    CheckVolume(ladder, infoFile);
    CheckVolume(ladder->GetNode(0)->GetVolume(), infoFile);
  }
  // --------------------------------------------------------------------------


  // ----------------   Create cone   -----------------------------------------
  Double_t coneDz            = 1.64;
  TGeoVolume* coneSmallVolum = ConstructSmallCone(coneDz);
  if (!coneSmallVolum) Fatal("ConstructSmallCone", "Volume Cone not found");
  TGeoVolume* coneBigVolum = ConstructBigCone(coneDz);
  if (!coneBigVolum) Fatal("ConstructBigCone", "Volume Cone not found");
  // --------------------------------------------------------------------------


  // ----------------   Create stations   -------------------------------------
  //  Float_t statPos[8] = {30., 40., 50., 60., 70., 80., 90., 100.};
  //DE23  Float_t statPos[8] = {13.75, 28.25, 41.75, 50., 60., 70., 80., 90.};
  //  Float_t statPos[8] = {14., 28., 42., 50., 60., 70., 80., 90.};
  // 2023
  //  Float_t statPos[8] = {11.5, 28., 42., 50., 60., 70., 80., 90.};
  // 2024
  //
  // from cave (Jens)
  //  hier die Sensorpositionen in Z-Richtung für das mSTS Upgrade.
  //  Target to Station 0: dZ=170   165
  //  Station 0 to 1:      dZ=125   290
  //  Station 1 to 2:      dZ=135   425
  //
  // 2024
  
  //Float_t statPos[8] = {16.5, 29.0, 42.5, 50., 60., 70., 80., 90.};
  Float_t statPos[8] = {12.065, 26.50, 40.0, 50., 60., 70., 80., 90.};
  //
  //  Float_t statPos[8] = {30., 45., 50., 60., 70., 80., 90., 100.};

  cout << endl << endl;
  cout << "===> Creating stations...." << endl;
  infoFile << endl << "Stations: ";
  nLadders = 0;
  Int_t ladderTypes[20];
  Double_t statZ             = 0.;
  Double_t rHole             = 0.;
  TGeoBBox* statShape        = NULL;
  TGeoTranslation* statTrans = NULL;


  // --- Station 01: 8 ladders, type 3 2 2 1 1 2 2 3
  cout << endl;
  statZ                 = 30.;
  rHole                 = 2.0;
  nLadders              = 1;
  ladderTypes[0]        = 13;
  TGeoVolume* station01 = ConstructStation(0, nLadders, ladderTypes, rHole);

  CheckVolume(station01);
  CheckVolume(station01, infoFile);
  infoFile << "Position z = " << statPos[0] << endl;


  // --- Station 02: 12 ladders, type 4 3 3 2 2 1 1 2 2 3 3 4
  cout << endl;
  statZ                 = 40.;
  rHole                 = 2.0;
  nLadders              = 3;
  nLadders              = 2;
  ladderTypes[0]        = 9;
  ladderTypes[1]        = 9;
// v22e  ladderTypes[0]        = 10;
// v22e  ladderTypes[1]        = 10;
// v22e  ladderTypes[2]        = 11;  // triple ladder
  TGeoVolume* station02 = ConstructStation(1, nLadders, ladderTypes, rHole);

  CheckVolume(station02);
  CheckVolume(station02, infoFile);
  infoFile << "Position z = " << statPos[1] << endl;


  // --- Station 03: 12 ladders, type 8 7 6 6 6 5 5 6 6 6 7 8
  cout << endl;
  statZ = 50.;
  rHole = 2.9;
  nLadders              = 3;
  ladderTypes[0]        = 10;
  ladderTypes[1]        = 12;
  ladderTypes[2]        = 11;  // triple ladder
  TGeoVolume* station03 = ConstructStation(2, nLadders, ladderTypes, rHole);

  CheckVolume(station03);
  CheckVolume(station03, infoFile);
  infoFile << "Position z = " << statPos[2] << endl;


  //  // --- Station 04: 14 ladders, type 9 8 7 6 6 6 5 5 6 6 7 8 9
  //  cout << endl;
  //  statZ = 60.;
  //  rHole = 2.9;
  //  nLadders = 14;
  //  ladderTypes[0]  = 15;  // 42;  // 9;
  //  ladderTypes[1]  = 14;  // 34;  // 8;
  //  ladderTypes[2]  = 13;  // 33;  // 7;
  //  ladderTypes[3]  = 12;  // 32;  // 6;
  //  ladderTypes[4]  = 12;  // 32;  // 6;
  //  ladderTypes[5]  = 12;  // 32;  // 6;
  //  ladderTypes[6]  =  4;  // 41;  // 5;
  //  ladderTypes[7]  =  4;  // 41;  // 5;
  //  ladderTypes[8]  = 12;  // 32;  // 6;
  //  ladderTypes[9]  = 12;  // 32;  // 6;
  //  ladderTypes[10] = 12;  // 32;  // 6;
  //  ladderTypes[11] = 13;  // 33;  // 7;
  //  ladderTypes[12] = 14;  // 34;  // 8;
  //  ladderTypes[13] = 15;  // 42;  // 9;
  //  TGeoVolume* station04 = ConstructStation(3, nLadders, ladderTypes, rHole);
  //
  //  if (gkConstructCones) {
  //    // upstream
  //    TGeoRotation* coneRot41 = new TGeoRotation;
  //    coneRot41->RotateZ(-90);
  //    coneRot41->RotateY(180);
  //    //    TGeoCombiTrans* conePosRot41 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.3-gkLadderGapZ/2., coneRot41);
  //    TGeoCombiTrans* conePosRot41 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.285-gkLadderGapZ/2., coneRot41);
  //    station04->AddNode(coneBigVolum, 1, conePosRot41);
  //
  //    // downstream
  //    TGeoRotation* coneRot42 = new TGeoRotation;
  //    coneRot42->RotateZ(-90);
  //    //    TGeoCombiTrans* conePosRot42 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.3+gkLadderGapZ/2., coneRot42);
  //    TGeoCombiTrans* conePosRot42 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.285+gkLadderGapZ/2., coneRot42);
  //    station04->AddNode(coneBigVolum, 2, conePosRot42);
  //
  //    station04->GetShape()->ComputeBBox();
  //  }
  //
  //  CheckVolume(station04);
  //  CheckVolume(station04, infoFile);
  //  infoFile << "Position z = " << statPos[3] << endl;
  //
  //
  //  // --- Station 05: 14 ladders, type 14 13 12 12 11 11 10 10 11 11 12 12 13 14
  //  cout << endl;
  //  statZ = 70.;
  //  rHole = 3.7;
  //  nLadders = 14;
  //  ladderTypes[0]  = 19;  //  55;  // 14;
  //  ladderTypes[1]  = 18;  //  54;  // 13;
  //  ladderTypes[2]  = 17;  //  53;  // 12;
  //  ladderTypes[3]  = 17;  //  53;  // 12;
  //  ladderTypes[4]  = 16;  //  52;  // 11;
  //  ladderTypes[5]  = 16;  //  52;  // 11;
  //  ladderTypes[6]  =  5;  //  51;  // 23;   // 10;
  //  ladderTypes[7]  =  5;  //  51;  // 23;   // 10;
  //  ladderTypes[8]  = 16;  //  52;  // 11;
  //  ladderTypes[9]  = 16;  //  52;  // 11;
  //  ladderTypes[10] = 17;  //  53;  // 12;
  //  ladderTypes[11] = 17;  //  53;  // 12;
  //  ladderTypes[12] = 18;  //  54;  // 13;
  //  ladderTypes[13] = 19;  //  55;  // 14;
  //  TGeoVolume* station05 = ConstructStation(4, nLadders, ladderTypes, rHole);
  //
  //  if (gkConstructCones) {
  //    // upstream
  //    TGeoRotation* coneRot51 = new TGeoRotation;
  //    coneRot51->RotateZ(90);
  //    coneRot51->RotateY(180);
  //    //    TGeoCombiTrans* conePosRot51 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.3-gkLadderGapZ/2., coneRot51);
  //    TGeoCombiTrans* conePosRot51 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.285-gkLadderGapZ/2., coneRot51);
  //    station05->AddNode(coneBigVolum, 1, conePosRot51);
  //
  //    // downstream
  //    TGeoRotation* coneRot52 = new TGeoRotation;
  //    coneRot52->RotateZ(90);
  //    //    TGeoCombiTrans* conePosRot52 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.3+gkLadderGapZ/2., coneRot52);
  //    TGeoCombiTrans* conePosRot52 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.285+gkLadderGapZ/2., coneRot52);
  //    station05->AddNode(coneBigVolum, 2, conePosRot52);
  //
  //    station05->GetShape()->ComputeBBox();
  //  }
  //
  //  CheckVolume(station05);
  //  CheckVolume(station05, infoFile);
  //  infoFile << "Position z = " << statPos[4] << endl;
  //
  //
  //  // --- Station 06: 14 ladders, type 14 13 12 12 11 11 10 10 11 11 12 12 13 14
  //  cout << endl;
  //  statZ = 80.;
  //  rHole = 3.7;
  //  nLadders = 14;
  //  ladderTypes[0]  = 19;  // 55;  // 14;
  //  ladderTypes[1]  = 18;  // 54;  // 13;
  //  ladderTypes[2]  = 17;  // 53;  // 12;
  //  ladderTypes[3]  = 17;  // 53;  // 12;
  //  ladderTypes[4]  = 16;  // 52;  // 11;
  //  ladderTypes[5]  = 16;  // 52;  // 11;
  //  ladderTypes[6]  =  6;  // 61;  // 10;
  //  ladderTypes[7]  =  6;  // 61;  // 10;
  //  ladderTypes[8]  = 16;  // 52;  // 11;
  //  ladderTypes[9]  = 16;  // 52;  // 11;
  //  ladderTypes[10] = 17;  // 53;  // 12;
  //  ladderTypes[11] = 17;  // 53;  // 12;
  //  ladderTypes[12] = 18;  // 54;  // 13;
  //  ladderTypes[13] = 19;  // 55;  // 14;
  //  TGeoVolume* station06 = ConstructStation(5, nLadders, ladderTypes, rHole);
  //
  //  if (gkConstructCones) {
  //    // upstream
  //    TGeoRotation* coneRot61 = new TGeoRotation;
  //    coneRot61->RotateZ(-90);
  //    coneRot61->RotateY(180);
  //    //    TGeoCombiTrans* conePosRot61 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.3-gkLadderGapZ/2., coneRot61);
  //    TGeoCombiTrans* conePosRot61 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.285-gkLadderGapZ/2., coneRot61);
  //    station06->AddNode(coneBigVolum, 1, conePosRot61);
  //
  //    // downstream
  //    TGeoRotation* coneRot62 = new TGeoRotation;
  //    coneRot62->RotateZ(-90);
  //    //    TGeoCombiTrans* conePosRot62 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.3+gkLadderGapZ/2., coneRot62);
  //    TGeoCombiTrans* conePosRot62 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.285+gkLadderGapZ/2., coneRot62);
  //    station06->AddNode(coneBigVolum, 2, conePosRot62);
  //
  //    station06->GetShape()->ComputeBBox();
  //  }
  //
  //  CheckVolume(station06);
  //  CheckVolume(station06, infoFile);
  //  infoFile << "Position z = " << statPos[5] << endl;
  //
  //
  //  // --- Station 07: 16 ladders, type 14 13 17 17 16 16 16 15 15 16 16 16 17 17 13 14
  //  cout << endl;
  //  statZ = 90.;
  //  rHole = 4.2;
  //  nLadders = 16;
  //  ladderTypes[0]  = 21;  // 73;  // 17;
  //  ladderTypes[1]  = 19;  // 55;  // 14;
  //  ladderTypes[2]  = 18;  // 54;  // 13;
  //  ladderTypes[3]  = 20;  // 72;  // 16;
  //  ladderTypes[4]  = 20;  // 72;  // 16;
  //  ladderTypes[5]  = 20;  // 72;  // 16;
  //  ladderTypes[6]  = 20;  // 72;  // 16;
  //  ladderTypes[7]  =  7;  // 71;  // 15;
  //  ladderTypes[8]  =  7;  // 71;  // 15;
  //  ladderTypes[9]  = 20;  // 72;  // 16;
  //  ladderTypes[10] = 20;  // 72;  // 16;
  //  ladderTypes[11] = 20;  // 72;  // 16;
  //  ladderTypes[12] = 20;  // 72;  // 16;
  //  ladderTypes[13] = 18;  // 54;  // 13;
  //  ladderTypes[14] = 19;  // 55;  // 14;
  //  ladderTypes[15] = 21;  // 73;  // 17;
  //  TGeoVolume* station07 = ConstructStation(6, nLadders, ladderTypes, rHole);
  //
  //  if (gkConstructCones) {
  //    // upstream
  //    TGeoRotation* coneRot71 = new TGeoRotation;
  //    coneRot71->RotateZ(90);
  //    coneRot71->RotateY(180);
  //    //    TGeoCombiTrans* conePosRot71 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.3-gkLadderGapZ/2., coneRot71);
  //    TGeoCombiTrans* conePosRot71 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.285-gkLadderGapZ/2., coneRot71);
  //    station07->AddNode(coneBigVolum, 1, conePosRot71);
  //
  //    // downstream
  //    TGeoRotation* coneRot72 = new TGeoRotation;
  //    coneRot72->RotateZ(90);
  //    //    TGeoCombiTrans* conePosRot72 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.3+gkLadderGapZ/2., coneRot72);
  //    TGeoCombiTrans* conePosRot72 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.285+gkLadderGapZ/2., coneRot72);
  //    station07->AddNode(coneBigVolum, 2, conePosRot72);
  //
  //    station07->GetShape()->ComputeBBox();
  //  }
  //
  //  CheckVolume(station07);
  //  CheckVolume(station07, infoFile);
  //  infoFile << "Position z = " << statPos[6] << endl;
  //
  //
  //  // --- Station 08: 16 ladders, type 14 13 17 17 16 16 16 15 15 16 16 16 17 17 13 14
  //  cout << endl;
  //  statZ = 100.;
  //  rHole = 4.2;
  //  nLadders = 16;
  //  ladderTypes[0]  = 19;  // 55;  // 14;
  //  ladderTypes[1]  = 17;  // 53;  // 12;
  //  ladderTypes[2]  = 23;  // 83;  // 20;
  //  ladderTypes[3]  = 22;  // 82;  // 19;
  //  ladderTypes[4]  = 22;  // 82;  // 19;
  //  ladderTypes[5]  = 22;  // 82;  // 19;
  //  ladderTypes[6]  = 22;  // 82;  // 19;
  //  ladderTypes[7]  =  8;  // 81;  // 18;
  //  ladderTypes[8]  =  8;  // 81;  // 18;
  //  ladderTypes[9]  = 22;  // 82;  // 19;
  //  ladderTypes[10] = 22;  // 82;  // 19;
  //  ladderTypes[11] = 22;  // 82;  // 19;
  //  ladderTypes[12] = 22;  // 82;  // 19;
  //  ladderTypes[13] = 23;  // 83;  // 20;
  //  ladderTypes[14] = 17;  // 53;  // 12;
  //  ladderTypes[15] = 19;  // 55;  // 14;
  //  TGeoVolume* station08 = ConstructStation(7, nLadders, ladderTypes, rHole);
  //
  //  if (gkConstructCones) {
  //    // upstream
  //    TGeoRotation* coneRot81 = new TGeoRotation;
  //    coneRot81->RotateZ(-90);
  //    coneRot81->RotateY(180);
  //    //    TGeoCombiTrans* conePosRot81 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.3-gkLadderGapZ/2., coneRot81);
  //    TGeoCombiTrans* conePosRot81 = new TGeoCombiTrans(name+"conePosRot2", 0., 0., -coneDz-0.285-gkLadderGapZ/2., coneRot81);
  //    station08->AddNode(coneBigVolum, 1, conePosRot81);
  //
  //    // downstream
  //    TGeoRotation* coneRot82 = new TGeoRotation;
  //    coneRot82->RotateZ(-90);
  //    //    TGeoCombiTrans* conePosRot82 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.3+gkLadderGapZ/2., coneRot82);
  //    TGeoCombiTrans* conePosRot82 = new TGeoCombiTrans(name+"conePosRot1", 0., 0., coneDz+0.285+gkLadderGapZ/2., coneRot82);
  //    station08->AddNode(coneBigVolum, 2, conePosRot82);
  //
  //    station08->GetShape()->ComputeBBox();
  //  }
  //
  //  CheckVolume(station08);
  //  CheckVolume(station08, infoFile);
  //  infoFile << "Position z = " << statPos[7] << endl;
  // --------------------------------------------------------------------------


  // ---------------   Create subplates   ---------------------------------------
  cout << endl << endl;
  cout << "===> Creating subplates...." << endl << endl;
  infoFile << endl << "Subplates: " << endl;
  Int_t nSubplates = CreateSubplates();
  for (Int_t iSubplate = 1; iSubplate <= nSubplates; iSubplate++) {
    TString name         = Form("Subplate%02d", iSubplate);
    TGeoVolume* subplate = gGeoMan->GetVolume(name);

    // add color to plates
    subplate->SetTransparency(60);
    if (iSubplate == 1) subplate->SetLineColor(kOrange);  // kOrange);
    if (iSubplate == 2) subplate->SetLineColor(kOrange);  // kRed);
    if (iSubplate == 3) subplate->SetLineColor(kOrange);  // kGreen);

    if (iSubplate == 4) subplate->SetLineColor(kOrange);  // kRed);
    if (iSubplate == 5) subplate->SetLineColor(kOrange);  // kOrange);
    if (iSubplate == 6) subplate->SetLineColor(kOrange);  // kCyan);

    if (iSubplate == 7) subplate->SetLineColor(kOrange);  // kRed);

    if (iSubplate == 8) subplate->SetLineColor(kOrange);   // kGreen);
    if (iSubplate == 9) subplate->SetLineColor(kOrange);   // kGreen);
    if (iSubplate == 10) subplate->SetLineColor(kOrange);  // kGreen);

    //    CheckVolume(plate);
    //    CheckVolume(plate, infoFile);
  }
  // --------------------------------------------------------------------------

  // ---------------   Create plates   ---------------------------------------
  cout << endl << endl;
  cout << "===> Creating plates...." << endl << endl;
  infoFile << endl << "Plates: " << endl;
  Int_t nPlates = CreatePlates();
  for (Int_t iPlate = 1; iPlate <= nPlates; iPlate++) {
    TString name      = Form("Plate%02d", iPlate);
    TGeoVolume* plate = gGeoMan->GetVolume(name);
    //    CheckVolume(plate);
    //    CheckVolume(plate, infoFile);
  }
  // --------------------------------------------------------------------------


  // ---------------   Create STS volume   ------------------------------------
  cout << endl << endl;
  cout << "===> Creating STS...." << endl;

  TString stsName = "sts_";
  stsName += geoTag;

  // --- Determine size of STS box
  Double_t stsX      = 0.;
  Double_t stsY      = 0.;
  Double_t stsZ      = 0.;
  Double_t stsBorder = 2 * 5.;  // 5 cm space for carbon ladders on each side
  Int_t nStation     = 3;       // set number of stations
  for (Int_t iStation = 1; iStation <= nStation; iStation++) {
    TString statName    = Form("Station%02d", iStation);
    TGeoVolume* station = gGeoMan->GetVolume(statName);
    TGeoBBox* shape     = (TGeoBBox*) station->GetShape();
    stsX                = TMath::Max(stsX, 2. * shape->GetDX());
    stsY                = TMath::Max(stsY, 2. * shape->GetDY());
    cout << "Station " << iStation << ":  Y " << stsY << endl;
  }
  // --- Some border around the stations
  stsX += stsBorder;
  stsY += stsBorder;
  stsZ             = (statPos[2] - statPos[1]) + stsBorder;
  Double_t stsPosZ = 0.5 * (statPos[2] + statPos[1]);
//  stsZ             = (statPos[1] - statPos[0]) + stsBorder;
//  Double_t stsPosZ = 0.5 * (statPos[1] + statPos[0]);

  // --- Create box  around the stations
  TGeoBBox* stsBox = new TGeoBBox("stsBox", stsX / 2., stsY / 2., stsZ / 2.);
  cout << "size of STS box: x " << stsX << " - y " << stsY << " - z " << stsZ << endl;

  //  // --- Create cone hosting the beam pipe
  //  // --- One straight section with constant radius followed by a cone
  //  Double_t z1 = statPos[0] - 0.5 * stsBorder;  // start of STS box
  //  Double_t z2 = gkPipeZ2;
  //  Double_t z3 = statPos[1] + 0.5 * stsBorder;  // end of STS box
  //  Double_t r1 = BeamPipeRadius(z1);
  //  Double_t r2 = BeamPipeRadius(z2);
  //  Double_t r3 = BeamPipeRadius(z3);
  //  r1 += 0.01;    // safety margin
  //  r2 += 0.01;    // safety margin
  //  r3 += 0.01;    // safety margin
  //
  //  cout << endl;
  //  cout << z1 << "  " << r1 << endl;
  //  cout << z2 << "  " << r2 << endl;
  //  cout << z3 << "  " << r3 << endl;
  //
  //  cout << endl;
  //  cout << "station1 :  " << BeamPipeRadius(statPos[0]) << endl;
  //  cout << "station2 :  " << BeamPipeRadius(statPos[1]) << endl;
  //  cout << "station3 :  " << BeamPipeRadius(statPos[2]) << endl;
  //  cout << "station4 :  " << BeamPipeRadius(statPos[3]) << endl;
  //  cout << "station5 :  " << BeamPipeRadius(statPos[4]) << endl;
  //  cout << "station6 :  " << BeamPipeRadius(statPos[5]) << endl;
  //  cout << "station7 :  " << BeamPipeRadius(statPos[6]) << endl;
  //  cout << "station8 :  " << BeamPipeRadius(statPos[7]) << endl;
  //
  //  //  TGeoPcon* cutout = new TGeoPcon("stsCone", 0., 360., 3); // 2.*TMath::Pi(), 3);
  //  //  cutout->DefineSection(0, z1, 0., r1);
  //  //  cutout->DefineSection(1, z2, 0., r2);
  //  //  cutout->DefineSection(2, z3, 0., r3);
  //  new TGeoTrd2("stsCone1", r1, r2, r1, r2, (z2-z1)/2.+.1);  // add .1 in z length for a clean cutout
  //  TGeoTranslation *trans1 = new TGeoTranslation("trans1", 0., 0., -(z3-z1)/2.+(z2-z1)/2.);
  //  trans1->RegisterYourself();
  //  new TGeoTrd2("stsCone2", r2, r3, r2, r3, (z3-z2)/2.+.1);  // add .1 in z length for a clean cutout
  //  TGeoTranslation *trans2 = new TGeoTranslation("trans2", 0., 0., +(z3-z1)/2.-(z3-z2)/2.);
  //  trans2->RegisterYourself();

  //DE   Double_t z1 = statPos[0] - 0.5 * stsBorder;  // start of STS box
  //DE   Double_t z2 = statPos[7] + 0.5 * stsBorder;  // end of STS box
  //DE   Double_t slope = (gkPipeR2 - gkPipeR1) / (gkPipeZ2 - gkPipeZ1);
  //DE   Double_t r1 = gkPipeR1 + slope * (z1 - gkPipeZ1); // at start of STS
  //DE   Double_t r2 = gkPipeR1 + slope * (z2 - gkPipeZ1); // at end of STS
  //DE   r1 += 0.1;    // safety margin
  //DE   r2 += 0.1;    // safety margin
  //DE   //  new TGeoCone("stsCone", stsZ/2., 0., r1, 0., r2);
  //DE   new TGeoTrd2("stsCone", r1, r2, r1, r2, stsZ/2.);

  // --- Create STS volume
  //  TGeoShape* stsShape = new TGeoCompositeShape("stsShape",
  //                                               "stsBox-stsCone1:trans1-stsCone2:trans2");
  //  TGeoVolume* sts = new TGeoVolume(stsName.Data(), stsShape, gStsMedium);
  //  TGeoVolume* sts = new TGeoVolume(stsName.Data(), stsBox, gStsMedium);
  TGeoVolumeAssembly* sts = new TGeoVolumeAssembly(stsName.Data());  // do not produce keeping volumes

  // --- Place stations in the STS
  for (Int_t iStation = 1; iStation <= nStation; iStation++) {
    TString statName       = Form("Station%02d", iStation);
    TGeoVolume* station    = gGeoMan->GetVolume(statName);
    Double_t posZ          = statPos[iStation - 1] - stsPosZ;
    TGeoTranslation* trans = new TGeoTranslation(0., 0., posZ);  // standard
    sts->AddNode(station, iStation, trans);
    sts->GetShape()->ComputeBBox();
  }


  // plates
  if (IncludeBox)
    for (Int_t iPlate = 1; iPlate <= nPlates; iPlate++) {
      TString platName  = Form("Plate%02d", iPlate);
      TGeoVolume* plate = gGeoMan->GetVolume(platName);
      TGeoTranslation* boxdz01 = new TGeoTranslation("ty01", 0, 0, 2.50); // v23a - shift to 2x2 and 3x3 stations
      sts->AddNode(plate, iPlate, boxdz01);
  //      sts->AddNode(plate, iPlate);
      sts->GetShape()->ComputeBBox();
    }

  cout << endl;
  CheckVolume(sts);
  // --------------------------------------------------------------------------


  // ---------------   Finish   -----------------------------------------------
  TGeoTranslation* stsTrans = new TGeoTranslation(gOffX, gOffY, stsPosZ + gOffZ);
  top->AddNode(sts, 1, stsTrans);
  top->GetShape()->ComputeBBox();
  cout << endl << endl;
  CheckVolume(top);
  cout << endl << endl;
  gGeoMan->CloseGeometry();
  gGeoMan->CheckOverlaps(0.0001);
  gGeoMan->PrintOverlaps();
  gGeoMan->Test();

  TFile* geoFile = new TFile(geoFileName, "RECREATE");
  sts->Export(geoFileName);
  cout << endl;
  cout << "Geometry " << sts->GetName() << " exported to " << geoFileName << endl;
  geoFile->Close();

  geoFile = new TFile(geoFileName, "UPDATE");
  stsTrans->Write();
  geoFile->Close();

  TString geoFileName_ = "sts_";
  geoFileName_         = geoFileName_ + geoTag + "_geo.root?reproducible";

  geoFile = new TFile(geoFileName_, "RECREATE");
  gGeoMan->Write();  // use this is you want GeoManager format in the output
  geoFile->Close();

  top->Draw("ogl");
  gGeoManager->SetVisLevel(6);

  infoFile.close();

  // create medialist for this geometry
  TString createmedialist = gSystem->Getenv("VMCWORKDIR");
  createmedialist += "/macro/geometry/create_medialist.C";
  std::cout << "Loading macro " << createmedialist << std::endl;
  gROOT->LoadMacro(createmedialist);
  gROOT->ProcessLine("create_medialist(\"\", false)");
}
// ============================================================================
// ======                   End of main function                          =====
// ============================================================================


// ****************************************************************************
// *****      Definition of media, sensors, sectors and ladders           *****
// *****                                                                  *****
// *****     Decoupled from main function for better readability          *****
// ****************************************************************************


/** ===========================================================================
 ** Create media
 **
 ** Currently created: air, active silicon, passive silion
 **
 ** Not used for the time being
 **/
Int_t CreateMedia()
{

  Int_t nMedia     = 0;
  Double_t density = 0.;

  // --- Material air
  density             = 1.205e-3;  // [g/cm^3]
  TGeoMixture* matAir = new TGeoMixture("sts_air", 3, density);
  matAir->AddElement(14.0067, 7, 0.755);  // Nitrogen
  matAir->AddElement(15.999, 8, 0.231);   // Oxygen
  matAir->AddElement(39.948, 18, 0.014);  // Argon

  // --- Material silicon
  density             = 2.33;  // [g/cm^3]