From d083431f4db98baead5d30ce59694ce6143ccabe Mon Sep 17 00:00:00 2001
From: Nicolas Koehler <nicolas.koehler@cern.ch>
Date: Wed, 30 Sep 2020 10:40:09 +0000
Subject: [PATCH] Modifications to support combined NSW&BIS78 layouts

---
 .../AmdcDb/src/AmdcDbSvcMakerFromAmdc.cxx     |  28 +++-
 .../MuonAGDDDescription/MMDetectorHelper.h    |   7 +-
 .../MuonAGDDDescription/sTGCDetectorHelper.h  |   7 +-
 .../src/MMDetectorHelper.cxx                  |  15 +-
 .../src/sTGCDetectorHelper.cxx                |  16 +-
 .../MuonReadoutGeometry/MMReadoutElement.h    |   1 +
 .../MuonReadoutGeometry/MuonChannelDesign.h   |  16 +-
 .../MuonReadoutGeometry/sTgcReadoutElement.h  |  15 +-
 .../src/MMReadoutElement.cxx                  |  15 +-
 .../src/sTgcReadoutElement.cxx                |  48 ++++--
 .../MM_Digitization/src/MM_StripResponse.cxx  |   2 +-
 .../MuonGeoModel/MuonGeoModel/MultiLayer.h    |   5 +
 MuonSpectrometer/MuonGeoModel/src/Mdt.cxx     | 102 +++++++++++-
 .../MuonGeoModel/src/Micromegas.cxx           |  15 +-
 .../MuonGeoModel/src/MultiLayer.cxx           | 152 +++++++++++-------
 .../MuonGeoModel/src/MuonChamber.cxx          |  15 +-
 MuonSpectrometer/MuonGeoModel/src/sTGC.cxx    |  13 +-
 17 files changed, 355 insertions(+), 117 deletions(-)

diff --git a/MuonSpectrometer/Amdcsimrec/AmdcDb/src/AmdcDbSvcMakerFromAmdc.cxx b/MuonSpectrometer/Amdcsimrec/AmdcDb/src/AmdcDbSvcMakerFromAmdc.cxx
index 89285e9eacf..0db7b7f9184 100755
--- a/MuonSpectrometer/Amdcsimrec/AmdcDb/src/AmdcDbSvcMakerFromAmdc.cxx
+++ b/MuonSpectrometer/Amdcsimrec/AmdcDb/src/AmdcDbSvcMakerFromAmdc.cxx
@@ -7,6 +7,9 @@
 #include "AmdcDb/AmdcDbSvc.h"
 #include "AmdcDb/AmdcDbRecordset.h"
 #include "AmdcDb/AmdcDbRecord.h"
+#include "GaudiKernel/MsgStream.h"
+#include "AthenaKernel/getMessageSvc.h"
+
 #include <cmath>
  
 AmdcDbSvcMakerFromAmdc::AmdcDbSvcMakerFromAmdc(){
@@ -1339,7 +1342,30 @@ if ((StationNameHEAD[0] != 'T'
               DbVar = "VERS"     ; DbVarComment="VERSION"                                  ; iDbVal = m_version		                                                                     ; pAmdcDbRecord->addInt(DbVar,DbVarComment,iDbVal);
 
               DbVar = "DX"       ; DbVarComment="X RELATIVE POSITION OF THE SUB-CUT"       ; dDbVal = pAmdcsimrec->Cutdx (pAmdcsimrec->INOCUT(DB_JTYP,DB_INDX,DB_ICUT),KounterCutLines)/ 10. ; pAmdcDbRecord->addDouble(DbVar,DbVarComment,dDbVal,LocalEpsLengthCM);
-              DbVar = "DY"       ; DbVarComment="Y RELATIVE POSITION OF THE SUB-CUT"       ; dDbVal = pAmdcsimrec->Cutdy (pAmdcsimrec->INOCUT(DB_JTYP,DB_INDX,DB_ICUT),KounterCutLines)/ 10. ; 
+              DbVar = "DY";
+              DbVarComment="Y RELATIVE POSITION OF THE SUB-CUT";
+              int theDbIdx = DB_INDX;
+              // when amdb is mirroring chambers from the Aside to the Cside, it increments the station numbers, e.g. BIS7 becomes BIS17 etc.
+              // however, something inside the INOCUT function (https://acode-browser1.usatlas.bnl.gov/lxr/source/athena/MuonSpectrometer/Amdcsimrec/AmdcStand/src/readmdb.F90#0611)
+              // cannot translate these larger station numbers to the correct amdb numbers needed by the amdb Cut_dy function 
+              // (https://acode-browser1.usatlas.bnl.gov/lxr/source/athena/MuonSpectrometer/Amdcsimrec/AmdcStand/src/readmdb.F90#0971)
+              // this Cut_dy function only returns the correct H Line dy values in case the original amdb station number is given.
+              // This only happens for BI upgrades, BIS78 and BI1-6
+              if (TheStationName.find("BI")!=std::string::npos) {
+                if (DB_INDX>16) {
+                  MsgStream log(Athena::getMessageSvc(),"AmdcDbSvcMakerFromAmdc::ALIN");
+                  log<<MSG::WARNING<<"ATTENTION: Adjusting DB_INDX value for H Line Cutdy method call. This should only be done for BI upgrade chambers."<<endmsg;
+                }
+                if (DB_INDX==17) theDbIdx=7;
+                else if (DB_INDX==18) theDbIdx=13;
+                else if (DB_INDX==19) theDbIdx=9;
+                else if (DB_INDX==20) theDbIdx=10;
+                else if (DB_INDX==21) theDbIdx=8;
+                else if (DB_INDX==24) theDbIdx=11;
+                else if (DB_INDX==25) theDbIdx=12;
+                else if (DB_INDX==26) theDbIdx=14;
+              }
+              dDbVal = pAmdcsimrec->Cutdy(pAmdcsimrec->INOCUT(DB_JTYP,theDbIdx,DB_ICUT),KounterCutLines)/10.; 
 // JFL Thu Apr  3 14:47:44 CEST 2008: 
 //  found by Stefania Spagnolo: When computed for the -Z part there are rounding errors which giving non nul value to Dy while it should be null
 //  patch: put it a 0. if too small
diff --git a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/MMDetectorHelper.h b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/MMDetectorHelper.h
index 64a6f043540..28eea163974 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/MMDetectorHelper.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/MMDetectorHelper.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef MMDetectorHelper_H
@@ -12,6 +12,7 @@ class MMDetectorDescription;
 class AGDDDetectorPositioner;
 
 typedef std::map<std::string,MMDetectorDescription*> MicromegasMap;
+typedef std::map<std::string,MMDetectorDescription*> MicromegasMapSubType;
 typedef std::map<std::string,MMDetectorDescription*>::const_iterator MicromegasIterator;
 
 typedef std::pair<MMDetectorDescription*,AGDDDetectorPositioner*> AGDDPositionedDetector;
@@ -22,12 +23,14 @@ public:
 	MicromegasIterator MM_begin() {return m_MicromegasList.begin();}
 	MicromegasIterator MM_end()   {return m_MicromegasList.end();}
 	
-	MMDetectorDescription* Get_MMDetectorType(std::string type);
+	MMDetectorDescription* Get_MMDetectorType(const std::string& type);
+    MMDetectorDescription* Get_MMDetectorSubType(const std::string& type);
 	MMDetectorDescription* Get_MMDetector(char type,int ieta,int iphi,int layer=1,char side='A');
 	AGDDPositionedDetector Get_MMPositionedDetector(char type,int ieta,int iphi,int layer=1,char side='A');
 
 private:
 	MicromegasMap m_MicromegasList;
+    MicromegasMapSubType m_MicromegasListSubType;
 
 };
 
diff --git a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/sTGCDetectorHelper.h b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/sTGCDetectorHelper.h
index af50614a8a5..e03fb9b9e30 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/sTGCDetectorHelper.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/MuonAGDDDescription/sTGCDetectorHelper.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef sTGCDetectorHelper_H
@@ -12,6 +12,7 @@ class sTGCDetectorDescription;
 class AGDDDetectorPositioner;
 
 typedef std::map<std::string,sTGCDetectorDescription*> sTGCMap;
+typedef std::map<std::string,sTGCDetectorDescription*> sTGCMapSubType;
 typedef std::map<std::string,sTGCDetectorDescription*>::const_iterator sTGCIterator;
 
 typedef std::pair<sTGCDetectorDescription*,AGDDDetectorPositioner*> AGDDPositionedDetector;
@@ -23,11 +24,13 @@ public:
 	sTGCIterator sTGC_end()   {return m_sTGCList.end();}
 	
 	sTGCDetectorDescription* Get_sTGCDetector(char type,int ieta,int iphi,int layer=1,char side='A');
-	sTGCDetectorDescription* Get_sTGCDetectorType(std::string type);
+	sTGCDetectorDescription* Get_sTGCDetectorType(const std::string& type);
+	sTGCDetectorDescription* Get_sTGCDetectorSubType(const std::string& type);
 	AGDDPositionedDetector Get_sTGCPositionedDetector(char type,int ieta,int iphi,int layer=1,char side='A');
 	
 private:
 	sTGCMap m_sTGCList;
+	sTGCMapSubType m_sTGCListSubType;
 
 };
 
diff --git a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/MMDetectorHelper.cxx b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/MMDetectorHelper.cxx
index 83da1fede83..3596e70d6ad 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/MMDetectorHelper.cxx
+++ b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/MMDetectorHelper.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "MuonAGDDDescription/MMDetectorHelper.h"
@@ -18,8 +18,10 @@ MMDetectorHelper::MMDetectorHelper()
 	for ( auto vl_iter: vl)
 	{
 		MMDetectorDescription* st=dynamic_cast<MMDetectorDescription*>(vl_iter.second);
-		if (st) 
+		if (st) {
 			m_MicromegasList[vl_iter.first]=st;
+			m_MicromegasListSubType[vl_iter.second->subType()]=st;
+		}
 	}
 	
 }
@@ -96,8 +98,15 @@ AGDDPositionedDetector MMDetectorHelper::Get_MMPositionedDetector(char type,int
 	return p_mm;
 }
 
-MMDetectorDescription* MMDetectorHelper::Get_MMDetectorType(std::string type)
+MMDetectorDescription* MMDetectorHelper::Get_MMDetectorType(const std::string& type)
 {
 	if (m_MicromegasList.find(type) != m_MicromegasList.end()) return m_MicromegasList[type];
 	return nullptr;
 }
+
+MMDetectorDescription* MMDetectorHelper::Get_MMDetectorSubType(const std::string& type)
+{
+	if (m_MicromegasListSubType.find(type) != m_MicromegasListSubType.end()) return m_MicromegasListSubType[type];
+	return nullptr;
+}
+
diff --git a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/sTGCDetectorHelper.cxx b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/sTGCDetectorHelper.cxx
index 6be96f08d3f..5e51b3499bc 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/sTGCDetectorHelper.cxx
+++ b/MuonSpectrometer/MuonDetDescr/MuonAGDDDescription/src/sTGCDetectorHelper.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "MuonAGDDDescription/sTGCDetectorHelper.h"
@@ -18,9 +18,10 @@ sTGCDetectorHelper::sTGCDetectorHelper()
 	for ( auto vl_iter: vl)
 	{
 		sTGCDetectorDescription* st=dynamic_cast<sTGCDetectorDescription*>(vl_iter.second);
-		//AGDDMicromegas* st1=dynamic_cast<AGDDMicromegas*>(st);
-		if (st) 
+		if (st) {
 			m_sTGCList[vl_iter.first]=st;
+			m_sTGCListSubType[vl_iter.second->subType()]=st;
+		}
 	}
 	
 }
@@ -97,8 +98,15 @@ AGDDPositionedDetector sTGCDetectorHelper::Get_sTGCPositionedDetector(char type,
 	return p_sTGC;
 }
 
-sTGCDetectorDescription* sTGCDetectorHelper::Get_sTGCDetectorType(std::string type)
+sTGCDetectorDescription* sTGCDetectorHelper::Get_sTGCDetectorType(const std::string& type)
 {
 	if (m_sTGCList.find(type) != m_sTGCList.end()) return m_sTGCList[type];
 	return nullptr;
 }
+
+sTGCDetectorDescription* sTGCDetectorHelper::Get_sTGCDetectorSubType(const std::string& type)
+{
+	if (m_sTGCListSubType.find(type) != m_sTGCListSubType.end()) return m_sTGCListSubType[type];
+	return nullptr;
+}
+
diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
index 90f2f19df3a..3f33043c2b9 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MMReadoutElement.h
@@ -138,6 +138,7 @@ namespace MuonGM {
     double m_halfX;   // 0.5*radial_size
     double m_minHalfY; // 0.5*bottom length (active area)
     double m_maxHalfY; // 0.5*top length (active area)
+    double m_offset;
 
     double m_rots;
     double m_rotz;
diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonChannelDesign.h b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonChannelDesign.h
index a64d221ae21..d611d4f9f07 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonChannelDesign.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonChannelDesign.h
@@ -38,7 +38,6 @@ namespace MuonGM {
     double deadO; //this param is not used for MM
     double deadS; //this param is not used for MM
     double signY;
-    //Amg::Vector2D firstChannelPos;
     double firstPos; //the position of the first active strip
     double firstPitch; // Pitch of 1st strip or number of wires in 1st group
     double groupWidth; // Number of Wires per group
@@ -64,9 +63,6 @@ namespace MuonGM {
     double dlStereoTop; // length between the first eta and stereo
     double dlStereoBottom;
     int totalStrips; //total strips per MM module
-      
-    /** channel transform */
-    //HepGeom::Transform3D  channelTransform( int channel ) const;
 
     /** distance to readout */
     double distanceToReadout( const Amg::Vector2D& pos ) const;
@@ -100,7 +96,7 @@ namespace MuonGM {
     
     int chNum = channelNumber( pos );
 
-    if (chNum <0 ) return -1.;
+    if (chNum <1 ) return -1.;
     Amg::Vector2D chPos;
     if (!channelPosition( chNum, chPos) ) return -1.;
 
@@ -123,7 +119,7 @@ namespace MuonGM {
     // if channel number is out of bounds, get the nearest channel ( mostly for validation purposes )
     bool validMode = false;
     if (type==MuonChannelDesign::etaStrip && detType==MuonChannelDesign::DetType::MM) {
-        if( chNum < 0 || chNum > totalStrips ){
+        if( chNum < 1 || chNum > totalStrips ){
             chNum = channelNumber(pos);
             validMode = true;
         }
@@ -396,7 +392,7 @@ namespace MuonGM {
             
             else if (detType==MuonChannelDesign::DetType::MM)
             {
-                if( (st >= nMissedBottomEta) || (st < (totalStrips-nMissedTopEta))){
+                if( (st > nMissedBottomEta) || (st < (totalStrips-nMissedTopEta))){
                     
                     stLen = inputLength + 2*(0.5*(maxYSize-minYSize)*(st-nMissedBottomEta)*inputPitch/xSize);
                     return stLen;
@@ -415,16 +411,16 @@ namespace MuonGM {
                
 
                  
-             if( st >= nMissedBottomStereo && st < (nMissedBottomStereo+nRoutedBottom) )
+             if( st > nMissedBottomStereo && st < (nMissedBottomStereo+nRoutedBottom) )
              stLen = (minYPhiR + (st-nMissedBottomStereo)*inputPitch)/sin(sAngle);
                  
-                 if(st >= (nMissedBottomStereo+nRoutedBottom) )
+                 if(st > (nMissedBottomStereo+nRoutedBottom) )
                  stLen = (inputLength + 2*(0.5*(maxYSize-minYSize)*(st-nMissedBottomEta)*inputPitch/xSize) )/cos(sAngle);
                  
                  return stLen;
              }
              //length for routed strips is not defined
-             else if(st>=(totalStrips-(nMissedTopStereo+nRoutedTop)) && st < (totalStrips-nMissedTopStereo))
+             else if(st>(totalStrips-(nMissedTopStereo+nRoutedTop)) && st < (totalStrips-nMissedTopStereo))
              {
                  stLen = ( maxYPhi + (totalStrips-nMissedTopStereo-(st+1))*inputPitch)/sin(sAngle);
                  return stLen;
diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/sTgcReadoutElement.h b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/sTgcReadoutElement.h
index 12e869bd182..f1dfc00daa6 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/sTgcReadoutElement.h
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/sTgcReadoutElement.h
@@ -62,6 +62,8 @@ namespace MuonGM {
 	If the strip number is outside the range of valid strips, the function will return false */
     virtual bool stripPosition( const Identifier& id, Amg::Vector2D& pos ) const override;
 
+    bool stripGlobalPosition( const Identifier& id, Amg::Vector3D& gpos ) const;
+
     /** pad number corresponding to local position */
     int padNumber( const Amg::Vector2D& pos, const Identifier& id) const;
 
@@ -164,6 +166,7 @@ namespace MuonGM {
     int m_nlayers;
     
     int m_ml;
+    double m_offset;
 
     int m_sTGC_type;
 
@@ -221,14 +224,9 @@ namespace MuonGM {
   }
 
   inline int sTgcReadoutElement::layerHash( const Identifier& id ) const {
-    //return layerHash(manager()->stgcIdHelper()->gasGap(id));
     return surfaceHash(id);                                      // don't have a choice here : rewrite MuonClusterReadoutElement first  
   }
 
-  //inline int sTgcReadoutElement::layerHash( int gasGap ) const {
-  //  return gasGap-1;
-  //}
-
   inline int sTgcReadoutElement::boundaryHash( const Identifier& id ) const {
     int iphi = manager()->stgcIdHelper()->channelType(id)!=1 ? 1:0 ;      // wires and pads have locX oriented along phi
     if (std::abs(getStationEta())<3) iphi += 2*(manager()->stgcIdHelper()->gasGap(id)-1);
@@ -293,6 +291,13 @@ namespace MuonGM {
 
   }
 
+  inline bool sTgcReadoutElement::stripGlobalPosition( const Identifier& id, Amg::Vector3D& gpos ) const {
+    Amg::Vector2D lpos(0., 0.);
+    if (!stripPosition(id, lpos)) return false;
+    surface(id).localToGlobal(lpos, Amg::Vector3D(0., 0., 0.), gpos);
+    return true;
+  }
+
   inline bool sTgcReadoutElement::padPosition( const Identifier& id, Amg::Vector2D& pos ) const {
 
     const MuonPadDesign* design = getPadDesign(id);
diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/MMReadoutElement.cxx b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/MMReadoutElement.cxx
index 9bc2b1a9c4c..a715cfbd843 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/MMReadoutElement.cxx
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/MMReadoutElement.cxx
@@ -35,6 +35,7 @@ namespace MuonGM {
     m_rots = 0.;
     m_rotz = 0.;
     m_rott = 0.;
+    m_offset = 0.;
 
     m_hasALines = false;
     m_hasBLines = false;
@@ -168,6 +169,8 @@ namespace MuonGM {
       MMDetectorDescription* mm = aHelper.Get_MMDetector(sector_l, std::abs(getStationEta()), getStationPhi(), m_ml, side);
       MMReadoutParameters roParam = mm->GetReadoutParameters();
 
+      double ylFrame = mm->ylFrame();
+      double ysFrame = mm->ysFrame();
 
       m_halfX = roParam.activeH/2;    //0.5*radial length (active area)
       m_minHalfY = roParam.activeBottomLength/2;   //0.5*bottom length (active area)
@@ -199,6 +202,9 @@ namespace MuonGM {
       m_etaDesign[il].maxYPhi = roParam.maxYPhi;
       m_etaDesign[il].totalStrips = roParam.tStrips;
       m_etaDesign[il].sAngle = (roParam.stereoAngle).at(il);
+      m_etaDesign[il].ylFrame = ylFrame;
+      m_etaDesign[il].ysFrame = ysFrame;
+      m_offset = -0.5*(m_etaDesign[il].ylFrame - m_etaDesign[il].ysFrame);
       if (m_ml < 1 || m_ml > 2) {
         MsgStream log(Athena::getMessageSvc(),"MMReadoutElement");
         log << MSG::WARNING <<"MMReadoutElement -- Unexpected Multilayer: m_ml= " << m_ml <<endmsg;
@@ -208,7 +214,7 @@ namespace MuonGM {
 
         m_etaDesign[il].firstPos = -0.5*m_etaDesign[il].xSize + pitch;
         m_etaDesign[il].signY  = 1 ;
-        m_etaDesign[il].nch = ((int) std::round( (m_etaDesign[il].xSize/pitch))) + 1; // Total number of active strips
+        m_etaDesign[il].nch = ((int) std::round( (m_etaDesign[il].xSize/pitch))); // Total number of active strips
 	
       } else { // stereo layers
           
@@ -229,7 +235,7 @@ namespace MuonGM {
         double fPos = -0.5*m_etaDesign[il].xSize - low_swift + lm1_swift;
         double lPos = 0.5*m_etaDesign[il].xSize + up_swift;
 
-        m_etaDesign[il].nch = ((int)std::round( (lPos - fPos)/pitch )) + 1;
+        m_etaDesign[il].nch = ((int)std::round( (lPos - fPos)/pitch ));
         m_etaDesign[il].firstPos = ( -0.5*m_etaDesign[il].xSize + (m_etaDesign[il].nMissedBottomStereo - m_etaDesign[il].nMissedBottomEta)*pitch) + pitch;
 
       }
@@ -263,7 +269,7 @@ namespace MuonGM {
       // keep the tracking surface at the center of the gap
       // need to operate switch x<->z because of GeoTrd definition
       m_surfaceData->m_layerSurfaces.push_back( new Trk::PlaneSurface(*this, id) );
-      m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
+      m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*Amg::Translation3D(0.,0.,m_offset)*
 						 Amg::AngleAxis3D(-90.*CLHEP::deg,Amg::Vector3D(0.,1.,0.)) );
 
       // surface info (center, normal) 
@@ -298,8 +304,7 @@ namespace MuonGM {
   {
     int gg = manager()->mmIdHelper()->gasGap(id);
     
-    //    Amg::Vector3D  locP = (m_Xlg[gg-1].inverse())*locPos;
-    Amg::Vector3D  locP = (m_Xlg[gg-1])*locPos;
+    Amg::Vector3D  locP = (m_Xlg[gg-1])*Amg::Translation3D(0.,0.,m_offset)*locPos;
 #ifndef NDEBUG
     MsgStream log(Athena::getMessageSvc(),"MMReadoutElement");
     if (log.level()<=MSG::DEBUG) {
diff --git a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/sTgcReadoutElement.cxx b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/sTgcReadoutElement.cxx
index 9eb83dc2a66..c49edf66f2a 100644
--- a/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/sTgcReadoutElement.cxx
+++ b/MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry/src/sTgcReadoutElement.cxx
@@ -44,6 +44,7 @@ namespace MuonGM {
     m_hasALines = false;
     m_hasBLines = false;
     m_delta = Amg::Transform3D::Identity();
+    m_offset = 0.;
     m_ml = mL;
     
     // get the setting of the caching flag from the manager
@@ -240,6 +241,9 @@ namespace MuonGM {
       m_etaDesign[il].inputPitch = stgc->stripPitch();
       m_etaDesign[il].inputLength = m_etaDesign[il].minYSize;
       m_etaDesign[il].inputWidth = stgc->stripWidth();
+      // If the top and bottom frames are not the same widths, the active geometry is incorrectly positionned by half the difference
+      m_offset = 0.5*(m_etaDesign[il].ylFrame-m_etaDesign[il].ysFrame);
+
       if (!tech){
         throw std::runtime_error(Form("File: %s, Line: %d\nsTgcReadoutElement::initDesign() - Failed To get Technology for stgc element: %s", __FILE__, __LINE__, stgc->GetName().c_str()));
       }
@@ -396,9 +400,6 @@ namespace MuonGM {
         m_surfaceData->m_surfBounds.push_back( new Trk::RotatedTrapezoidBounds( m_halfX[layer], m_minHalfY[layer], m_maxHalfY[layer]));  // strips
         m_surfaceData->m_surfBounds.push_back( new Trk::TrapezoidBounds( m_PadminHalfY[layer], m_PadmaxHalfY[layer], m_PadhalfX[layer]));
       }
-      
-      // If the top and bottom frames are not the same widths, the active geometry is incorrectly positionned by half the difference
-      double offset = 0.5*(m_etaDesign[layer].ylFrame-m_etaDesign[layer].ysFrame);
 
       // identifier of the first channel - wire plane - locX along phi, locY max->min R
       Identifier id = manager()->stgcIdHelper()->channelID(getStationName(),getStationEta(),getStationPhi(),m_ml, layer+1, 2, 1);
@@ -407,12 +408,12 @@ namespace MuonGM {
       m_surfaceData->m_layerSurfaces.push_back( new Trk::PlaneSurface(*this, id) );
       if (m_sTGC_type == 1 || m_sTGC_type == 2)
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-						 Amg::Translation3D(0,0.,-offset)*
+						 Amg::Translation3D(0,0.,-m_offset)*
 						 Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.))*
 						 Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,0.,1.)) );
       else if (m_sTGC_type == 3) // if QL3, diamond. have to shift geometry to account for origin not being in center
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-					Amg::Translation3D(0,0.,-offset + m_PadhalfX[layer] - m_padDesign[layer].yCutout)*
+					Amg::Translation3D(0,0.,-m_offset + m_PadhalfX[layer] - m_padDesign[layer].yCutout)*
 					Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.))*
 					Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,0.,1.)) );
       else throw std::runtime_error(Form("File: %s, Line: %d\nsTgcReadoutElement::fillCache() - sTGC_type %d is not valid! Wire Geometry not Created!", __FILE__, __LINE__, m_sTGC_type));
@@ -436,12 +437,12 @@ namespace MuonGM {
 
       if (m_sTGC_type == 1 || m_sTGC_type == 2)
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-						 Amg::Translation3D(shift,0.,-offset)*
+						 Amg::Translation3D(shift,0.,-m_offset)*
 						 Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.)));
 
       else if (m_sTGC_type == 3) // if QL3, diamond. have to shift geometry to account for origin not being in center
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-					Amg::Translation3D(shift,0.,-offset + m_halfX[layer] - m_etaDesign[layer].yCutout)*
+					Amg::Translation3D(shift,0.,-m_offset + m_halfX[layer] - m_etaDesign[layer].yCutout)*
 					Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.)) );
       else std::runtime_error(Form("File: %s, Line: %d\nsTgcReadoutElement::fillCache() - sTGC_type %d is not valid! Strip Geometry not Created!", __FILE__, __LINE__, m_sTGC_type));
 
@@ -456,12 +457,12 @@ namespace MuonGM {
       m_surfaceData->m_layerSurfaces.push_back( new Trk::PlaneSurface(*this, id) );
       if (m_sTGC_type == 1 || m_sTGC_type == 2)
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-						 Amg::Translation3D(-shift,0.,-offset)*
+						 Amg::Translation3D(-shift,0.,-m_offset)*
 						 Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.))*
 						 Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,0.,1.)) );
       else if (m_sTGC_type == 3) // if QL3, diamond. have to shift geometry to account for origin not being in center
         m_surfaceData->m_layerTransforms.push_back(absTransform()*m_Xlg[layer]*
-					Amg::Translation3D(-shift,0.,-offset + m_PadhalfX[layer] - m_padDesign[layer].yCutout)*
+					Amg::Translation3D(-shift,0.,-m_offset + m_PadhalfX[layer] - m_padDesign[layer].yCutout)*
 					Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,1.,0.))*
 					Amg::AngleAxis3D(-90*CLHEP::deg,Amg::Vector3D(0.,0.,1.)) );
       else std::runtime_error(Form("File: %s, Line: %d\nsTgcReadoutElement::fillCache() - sTGC_type %d is not valid! Pad Geometry not Created!", __FILE__, __LINE__, m_sTGC_type));
@@ -500,7 +501,34 @@ namespace MuonGM {
   Amg::Vector3D sTgcReadoutElement::localToGlobalCoords(Amg::Vector3D locPos, Identifier id) const
   {
     int gg = manager()->stgcIdHelper()->gasGap(id);
-    Amg::Vector3D  locP = m_Xlg[gg-1]*locPos;
+    int channelType = manager()->stgcIdHelper()->channelType(id);
+
+    Amg::Vector3D locP(0,0,0);
+    // strip plane moved along normal, pad plane in the opposite direction
+    // We no longer want the readout elements to be seperated by the gas gas volume
+    // We place all 3 readouts at the center of the gas gap in z, with a 10 micron offset to seperate them
+    double shift = 0.01;
+    if((gg-1)%2) shift = -shift;
+
+    if(channelType == 0){ //if pad plane
+      if(m_sTGC_type == 1 || m_sTGC_type == 2)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(-shift,0.,-m_offset)*locPos;
+      else if(m_sTGC_type == 3)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(-shift,0.,-m_offset + m_PadhalfX[gg-1] - m_padDesign[gg-1].yCutout)*locPos;
+    }
+    else if( channelType == 1 ){ //if strip plane
+      if(m_sTGC_type == 1 || m_sTGC_type == 2)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(shift,0.,-m_offset)*locPos;
+      else if(m_sTGC_type == 3)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(shift,0.,-m_offset + m_halfX[gg-1] - m_etaDesign[gg-1].yCutout)*locPos;
+    }
+    else if(channelType == 2){ //if wire plane 
+      if(m_sTGC_type == 1 || m_sTGC_type == 2)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(0,0.,-m_offset)*locPos;
+      else if(m_sTGC_type == 3)
+        locP = m_Xlg[gg-1]*Amg::Translation3D(0,0.,-m_offset + m_PadhalfX[gg-1] - m_padDesign[gg-1].yCutout)*locPos;
+    }
+
 #ifndef NDEBUG
     MsgStream log(Athena::getMessageSvc(),"sTgcReadoutElement");
     if (log.level()<=MSG::DEBUG) {
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
index a1af25ef29a..663d14296df 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripResponse.cxx
@@ -107,7 +107,7 @@ void MM_StripResponse::calculateSummaries(float chargeThreshold) {
       // Last active strip numbrer is maxStripID-1
       if (stripVal < m_minstripID || stripVal > m_maxstripID-1) continue;
       // remove PCB gap strips
-      if (stripVal == 1023 || stripVal == 1024 || stripVal == 2047 || stripVal == 2048 || stripVal == 3071 || stripVal == 3072 || stripVal == 4095 || stripVal == 4096) continue;
+      if (stripVal == 1024 || stripVal == 1025 || stripVal == 2048 || stripVal == 2049 || stripVal == 3072 || stripVal == 3073 || stripVal == 4096 || stripVal == 4097) continue;
 			float stripChargeVal = stripCharge.second;
 			if(stripChargeVal < chargeThreshold) continue;
 			
diff --git a/MuonSpectrometer/MuonGeoModel/MuonGeoModel/MultiLayer.h b/MuonSpectrometer/MuonGeoModel/MuonGeoModel/MultiLayer.h
index 6615cdbd61f..6035b7d52b5 100755
--- a/MuonSpectrometer/MuonGeoModel/MuonGeoModel/MultiLayer.h
+++ b/MuonSpectrometer/MuonGeoModel/MuonGeoModel/MultiLayer.h
@@ -37,6 +37,11 @@ public:
    double cutoutTubeLength[5]; // tube length
    bool cutoutFullLength[5];   // true if this region is outside the cutout
    bool cutoutAtAngle;      // true if this station has cutouts at an angle; //EMS1,3 and BOS6
+
+   // the same but for several cutouts along the amdb x (GeoModel y)
+   // gives the (x1,x2) and (y1,y2) tuples of rectangles which are NOT cutout
+   std::vector<std::pair<double,double>> m_nonCutoutXSteps;
+   std::vector<std::pair<double,double>> m_nonCutoutYSteps;
     
 };
 } // namespace MuonGM
diff --git a/MuonSpectrometer/MuonGeoModel/src/Mdt.cxx b/MuonSpectrometer/MuonGeoModel/src/Mdt.cxx
index 9a477a8a478..d33adce59ca 100755
--- a/MuonSpectrometer/MuonGeoModel/src/Mdt.cxx
+++ b/MuonSpectrometer/MuonGeoModel/src/Mdt.cxx
@@ -2,7 +2,6 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-
 #include "MuonGeoModel/Mdt.h"
 #include "MuonGeoModel/MDT_Technology.h"
 #include "MuonGeoModel/MYSQL.h"
@@ -12,10 +11,19 @@
 #include "GeoModelKernel/GeoFullPhysVol.h"
 #include "GeoModelKernel/GeoShape.h"
 
+#include <TString.h> // for Form
+
+#include "float.h"
+
 #define verbose_mdt false
 
 namespace MuonGM {
 
+float round(const float toRound, const unsigned int decimals) {
+  unsigned int factor = std::pow(10,decimals);
+  return std::round(toRound*factor)/factor;
+}
+
 Mdt::Mdt(Component* ss, std::string lVName): DetectorElement(ss->name)
 {
   logVolName = lVName;
@@ -65,6 +73,23 @@ GeoFullPhysVol* Mdt::build(std::vector<Cutout*> vcutdef)
     if (verbose_mdt) {
       log << MSG::VERBOSE << " mdt cutouts are on " << endmsg;
     }
+
+    bool cutoutsVsX=false;
+    bool cutoutsVsY=false;
+    // first check whether there are several cutouts with different amdb x or y coordinates
+    float lastX(FLT_MAX);
+    float lastY(FLT_MAX);
+    for (int i = 0; i < Ncuts; ++i) {
+      if (lastX!=FLT_MAX && lastX*vcutdef[i]->dx<0) cutoutsVsX=true;
+      lastX = vcutdef[i]->dx;
+      if (lastY!=FLT_MAX && lastY*vcutdef[i]->dy<0) cutoutsVsY=true;
+      lastY = vcutdef[i]->dy;
+    }
+    if (cutoutsVsX&&cutoutsVsY) {
+      throw std::runtime_error(Form("File: %s, Line: %d\nMdt::build() - Found more than one cutout in amdb-x direction and more than one cutout in amdb-y direction, currently not supported", __FILE__, __LINE__));
+    }
+
+if (!cutoutsVsX) { // nominal case (used for BMS/BOG/BMG etc.)
     int cutoutNtubes[5];           // Number of tubes in sub-multilayer[i]
     bool cutoutFullLength[5];      // True if this region is outside the cutout
     double cutoutXtubes[5];        // Location of tube center within sub-ml[i] along x-amdb
@@ -95,7 +120,7 @@ GeoFullPhysVol* Mdt::build(std::vector<Cutout*> vcutdef)
     int cutLocationCode[3] = {0, 0, 0};
     for (int i = 0; i < Ncuts; i++) {
       if (vcutdef[i]->dy <= 0) cutLocationCode[i] = -1;
-      if (vcutdef[i]->dy + vcutdef[i]->lengthY >= top) cutLocationCode[i] = 1;
+      if (round(vcutdef[i]->dy + vcutdef[i]->lengthY,2) >= round(top,2)) cutLocationCode[i] = 1;
     }
 
     // Calculate quantities needed by multilayer
@@ -260,6 +285,79 @@ GeoFullPhysVol* Mdt::build(std::vector<Cutout*> vcutdef)
     }
 
     layer->cutoutAtAngle = cutAtAngle;
+  } else {
+    // there are several cutouts along the amdb-x coordinate
+
+    if (longWidth!=width) {
+      throw std::runtime_error(Form("File: %s, Line: %d\nMdt::build() - only support cutouts along amdb-x for rectangular chambers", __FILE__, __LINE__));
+    }
+
+    std::vector<std::pair<double,double>> nonCutoutXSteps;
+    std::vector<std::pair<double,double>> nonCutoutYSteps;
+
+    // Order cutouts by increasing dx
+    for (int i = 0; i < Ncuts; i++) {
+      for (int j = i+1; j < Ncuts; j++) {
+        if (vcutdef[j]->dx < vcutdef[i]->dx) {
+          Cutout* c = vcutdef[i];
+          vcutdef[i] = vcutdef[j];
+          vcutdef[j] = c;
+        }
+      }
+    }
+
+    // in amdb-coordinates
+    double xminChamber = round(-width/2,2);
+    double xmaxChamber = round(width/2,2);
+    double yminChamber = 0;
+    double ymaxChamber = round(length-tubePitch/2,2);
+
+    double latestXMax=xminChamber;
+
+    for (int i = 0; i < Ncuts; ++i) {
+      Cutout* c = vcutdef[i];
+      double lowerX = round(c->dx - c->widthXs/2,2);
+      double xmin = (lowerX <= xminChamber) ? xminChamber : lowerX;
+      if (xmin<latestXMax) {
+        throw std::runtime_error(Form("File: %s, Line: %d\nMdt::build() - cannot have cutouts along amdb-x which overlap in amdb-x", __FILE__, __LINE__));
+      }
+      if (i==0 && xmin>xminChamber) {
+        // we start with a full slice without cutout
+        nonCutoutXSteps.push_back(std::make_pair(xminChamber,lowerX));
+        nonCutoutYSteps.push_back(std::make_pair(yminChamber,ymaxChamber));
+      }
+      double upperX = round(c->dx + c->widthXs/2,2);
+      double xmax = (upperX >= xmaxChamber) ? xmaxChamber : upperX;
+
+      double ymin = round(c->dy+c->lengthY,2) < ymaxChamber ? c->dy+c->lengthY : 0;
+      double ymax = ymaxChamber <= round(c->dy+c->lengthY,2) ? c->dy : ymaxChamber;
+
+      if (latestXMax<xmin) {
+        // there is a full slice between latestXMax and xmin
+        nonCutoutXSteps.push_back(std::make_pair(latestXMax,xmin));
+        nonCutoutYSteps.push_back(std::make_pair(yminChamber,ymaxChamber));
+      }
+
+      nonCutoutXSteps.push_back(std::make_pair(xmin,xmax));
+      nonCutoutYSteps.push_back(std::make_pair(ymin,ymax));
+
+
+      if (i==Ncuts-1 && xmax<xmaxChamber) {
+        // we end with a full slice without cutout
+        nonCutoutXSteps.push_back(std::make_pair(xmax,xmaxChamber));
+        nonCutoutYSteps.push_back(std::make_pair(yminChamber,ymaxChamber));
+      }
+
+      latestXMax = xmax;
+    }
+
+    // Pass information to multilayer and MdtComponent
+    m_component->cutoutTubeXShift = 0;
+    layer->cutoutNsteps = nonCutoutXSteps.size();
+    layer->m_nonCutoutXSteps = nonCutoutXSteps;
+    layer->m_nonCutoutYSteps = nonCutoutYSteps;
+    layer->cutoutAtAngle = false;
+  }
 
     return layer->build();
 
diff --git a/MuonSpectrometer/MuonGeoModel/src/Micromegas.cxx b/MuonSpectrometer/MuonGeoModel/src/Micromegas.cxx
index 1e82edc57c1..e56450176cf 100755
--- a/MuonSpectrometer/MuonGeoModel/src/Micromegas.cxx
+++ b/MuonSpectrometer/MuonGeoModel/src/Micromegas.cxx
@@ -2,9 +2,11 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-
 #include "MuonGeoModel/Micromegas.h"
+
 #include "MuonAGDDDescription/MM_Technology.h"
+#include "MuonAGDDDescription/MMDetectorDescription.h"
+#include "MuonAGDDDescription/MMDetectorHelper.h"
 #include "AGDDKernel/AGDDDetectorStore.h"
 #include "MuonGeoModel/Station.h"
 #include "MuonGeoModel/MicromegasComponent.h"
@@ -21,7 +23,6 @@
 #include "GeoModelKernel/GeoSerialIdentifier.h"
 #include "GeoModelKernel/GeoIdentifierTag.h"
 #include "GeoModelKernel/GeoDefinitions.h"
-// for cutouts:
 #include "GeoModelKernel/GeoShapeSubtraction.h"
 #include "GeoModelKernel/GeoShapeIntersection.h"
 #include "GeoModelKernel/GeoShapeShift.h"
@@ -49,15 +50,17 @@ GeoFullPhysVol* Micromegas::build(int minimalgeo)
 GeoFullPhysVol* Micromegas::build(int minimalgeo, int , std::vector<Cutout*> )
 {
   AGDDDetectorStore* ds = AGDDDetectorStore::GetDetectorStore();
+  MMDetectorHelper mmHelper;
+  MMDetectorDescription* mm_descr = mmHelper.Get_MMDetectorSubType(m_component->subType);
+
   MM_Technology* t = (MM_Technology*) ds->GetTechnology(name);
   thickness = t->Thickness();
   double gasTck=t->gasThickness;
   double pcbTck=t->pcbThickness;
   double roTck=t->roThickness;
-  double f1=t->f1Thickness;
-  double f2=t->f2Thickness;
-  double f3=t->f3Thickness;
-
+  double f1=mm_descr->ylFrame();
+  double f2=mm_descr->ysFrame();
+  double f3=mm_descr->xFrame(); 
 
   minimalgeo=t->geoLevel;
 
diff --git a/MuonSpectrometer/MuonGeoModel/src/MultiLayer.cxx b/MuonSpectrometer/MuonGeoModel/src/MultiLayer.cxx
index 9e39419bd90..6dee99d7e56 100755
--- a/MuonSpectrometer/MuonGeoModel/src/MultiLayer.cxx
+++ b/MuonSpectrometer/MuonGeoModel/src/MultiLayer.cxx
@@ -31,6 +31,7 @@
 
 #include <vector>
 #include <cassert>
+
 using namespace GeoXF;
 
 #define verbose_multilayer false
@@ -39,13 +40,14 @@ namespace MuonGM {
 
 MultiLayer::MultiLayer(std::string n): DetectorElement(n),
   nrOfLayers(0), nrOfTubes(0), tubePitch(0.), width(0.), length(0.), thickness(0.),
-  mdtthickness(0.), longWidth(0.), nrOfSteps(0), cutoutNsteps(0), cutoutAtAngle(false)
-{
+  mdtthickness(0.), longWidth(0.), nrOfSteps(0), cutoutNsteps(0), cutoutAtAngle(false),
+  m_nonCutoutXSteps(),
+  m_nonCutoutYSteps() {
   MsgStream log(Athena::getMessageSvc(), "MultiLayer::MultiLayer");
 
   MYSQL* mysql = MYSQL::GetPointer();
   MDT* md = (MDT*)mysql->GetTechnology(name);
-  if (md != NULL) {
+  if (md) {
     nrOfLayers = md->numOfLayers;
     mdtthickness = md->totalThickness;
     tubePitch = md->pitch;
@@ -86,9 +88,9 @@ GeoFullPhysVol* MultiLayer::build()
         << endmsg;
   }
 
-  const GeoShape* slay = NULL;
-  const GeoShape* sfoamup = NULL;   // do them both together when there are cutouts
-  const GeoShape* sfoamlow = NULL;  // do them both together when there are cutouts
+  const GeoShape* slay = nullptr;
+  const GeoShape* sfoamup = nullptr;   // do them both together when there are cutouts
+  const GeoShape* sfoamlow = nullptr;  // do them both together when there are cutouts
   // so calculate other foam-related info first:
   double foamthicknesslow = yy[0] - tuberad;
   double foamthicknessup;
@@ -106,7 +108,9 @@ GeoFullPhysVol* MultiLayer::build()
     } else if (fabs(foamthicknessup - 10.00*Gaudi::Units::mm) < 0.1
                && logVolName.find("BMG") != std::string::npos ) {
       foamthicknessup = 10.00*Gaudi::Units::mm;
-    } else if ( logVolName == "BME1MDT09" || logVolName == "BME2MDT09" ) { //@@
+    } else if ( logVolName.find("MDT09") != std::string::npos || logVolName.find("MDT14") != std::string::npos ) {
+      // MDT09 is used for BME and BIS sMDTs, MDT14 is used for 2nd multilayer of BIS78
+      // All those sMDTs have no foam layer
       foamthicknesslow = 0.;
       foamthicknessup  = 0.;
     } else {
@@ -127,7 +131,9 @@ GeoFullPhysVol* MultiLayer::build()
     } else if (fabs(foamthicknesslow - 10.00*Gaudi::Units::mm) < 0.1
              && logVolName.find("BMG") != std::string::npos ) {
       foamthicknesslow = 10.00*Gaudi::Units::mm;
-    } else if ( logVolName == "BME1MDT09" || logVolName == "BME2MDT09" ) { //@@
+    } else if ( logVolName.find("MDT09") != std::string::npos || logVolName.find("MDT14") != std::string::npos ) {
+      // MDT09 is used for BME and BIS sMDTs, MDT14 is used for 2nd multilayer of BIS78
+      // All those sMDTs have no foam layer
       foamthicknesslow = 0.;
       foamthicknessup  = 0.;
     } else {
@@ -143,6 +149,38 @@ GeoFullPhysVol* MultiLayer::build()
     // No cutouts - layer is a simple box or trapezoid
     slay = new GeoTrd(mdtthickness/2, mdtthickness/2, width/2,
                       longWidth/2, length/2);
+  } else if (m_nonCutoutXSteps.size()) {
+    // there are cutouts sliced along amdb-x (slices in amdb-y)
+    double submlthick = mdtthickness/2.; // corresponds to half amdb-z (chamber height=constant)
+
+    // loop on the cutout slices and create rectangular boxes which are NOT cutout and add them together to an envelope
+    for (unsigned int i = 0; i < m_nonCutoutXSteps.size(); ++i) {
+
+      double submlwidth = (m_nonCutoutXSteps.at(i).second-m_nonCutoutXSteps.at(i).first)/2; // corresponds to half amdb-x
+      double submllength = (m_nonCutoutYSteps.at(i).second-m_nonCutoutYSteps.at(i).first)/2; // corresponds to half amdb-y
+
+      // need a transform between center multilayer and center of trapezoid (y is amdb-x, z is amdb-y)
+      double yPos = (m_nonCutoutXSteps.at(i).first+m_nonCutoutXSteps.at(i).second)/2;
+      double zPos = (m_nonCutoutYSteps.at(i).first+m_nonCutoutYSteps.at(i).second)/2-length/2;
+      GeoTrf::Transform3D submlpos = GeoTrf::Translate3D(0.,yPos,zPos);
+
+      const GeoTrd* tempSLay = nullptr;
+      const GeoShape* tempSLay1 = nullptr;
+
+      if (submlthick*submlwidth*submllength > 0) {
+        tempSLay = new GeoTrd(submlthick, submlthick, submlwidth, submlwidth, submllength);
+        tempSLay1 = &( (*tempSLay)<<submlpos);
+      } else {
+        log << MSG::WARNING << " problem with shape of temporary trapezoid in LogVolName "
+            << logVolName << " thick,width,length=" << submlthick << "," << submlwidth << "," << submllength << endmsg;
+      }
+
+      if (slay) {
+        slay = &(slay->add(*tempSLay1));
+      } else {
+        slay = tempSLay1;
+      }
+    } // Loop over cutout steps
   } else {
     // Layer to be built as a boolean of boxes and trapezoids to represent cutouts
     if (verbose_multilayer) {
@@ -176,12 +214,12 @@ GeoFullPhysVol* MultiLayer::build()
       double widthPos = cutoutXtubes[isub];
       GeoTrf::Transform3D submlpos = GeoTrf::Translate3D(0.,widthPos,lengthPos);
 
-      const GeoTrd* tempSLay = NULL;
-      const GeoShape* tempSLay1 = NULL;
-      const GeoTrd* tempSFoamup = NULL;
-      const GeoShape* tempSFoamup1 = NULL;
-      const GeoTrd* tempSFoamlow = NULL;
-      const GeoShape* tempSFoamlow1 = NULL;
+      const GeoTrd* tempSLay = nullptr;
+      const GeoShape* tempSLay1 = nullptr;
+      const GeoTrd* tempSFoamup = nullptr;
+      const GeoShape* tempSFoamup1 = nullptr;
+      const GeoTrd* tempSFoamlow = nullptr;
+      const GeoShape* tempSFoamlow1 = nullptr;
 
       if (submlthick*submlwidthl*submllength > 0) {
         if (verbose_multilayer) {
@@ -210,38 +248,40 @@ GeoFullPhysVol* MultiLayer::build()
             << endmsg;
       }
 
-      if (foamthicknessup > foamthicknesslow) {
-        if (foamthicknessup*submlwidths*submllength > 0) {
-          if (cutoutFullLength[isub]) {
-            tempSFoamup = new GeoTrd(foamthicknessup/2., foamthicknessup/2.,
-                                     submlwidths, submlwidthl, submllength);
+      if (!(foamthicknessup==0 && foamthicknesslow==0)) {
+        if (foamthicknessup > foamthicknesslow) {
+          if (foamthicknessup*submlwidths*submllength > 0) {
+            if (cutoutFullLength[isub]) {
+              tempSFoamup = new GeoTrd(foamthicknessup/2., foamthicknessup/2.,
+                                       submlwidths, submlwidthl, submllength);
+            } else {
+              tempSFoamup = new GeoTrd(foamthicknessup/2., foamthicknessup/2.,
+                                       cutoutTubeLength[isub]/2., cutoutTubeLength[isub]/2.,
+                                       submllength);
+            }
+
+            tempSFoamup1 = &( (*tempSFoamup)<<submlpos);
           } else {
-            tempSFoamup = new GeoTrd(foamthicknessup/2., foamthicknessup/2.,
-                                     cutoutTubeLength[isub]/2., cutoutTubeLength[isub]/2.,
-                                     submllength);
+            log << MSG::INFO << " problem with shape of upper foam trapezoid in LogVolName "
+                << logVolName << " thick,width,length " << foamthicknessup
+                << " " << submlwidths << " " << submllength << endmsg;
           }
-
-          tempSFoamup1 = &( (*tempSFoamup)<<submlpos);
         } else {
-          log << MSG::INFO << " problem with shape of upper foam trapezoid in LogVolName "
-              << logVolName << " thick,width,length " << foamthicknessup
-              << " " << submlwidths << " " << submllength << endmsg;
-        }
-      } else {
-        if (foamthicknesslow*submlwidths*submllength > 0) {
-          if (cutoutFullLength[isub]) {
-            tempSFoamlow = new GeoTrd(foamthicknesslow/2., foamthicknesslow/2.,
-                                      submlwidths, submlwidthl, submllength);
+          if (foamthicknesslow*submlwidths*submllength > 0) {
+            if (cutoutFullLength[isub]) {
+              tempSFoamlow = new GeoTrd(foamthicknesslow/2., foamthicknesslow/2.,
+                                        submlwidths, submlwidthl, submllength);
+            } else {
+              tempSFoamlow = new GeoTrd(foamthicknesslow/2., foamthicknesslow/2.,
+                                        cutoutTubeLength[isub]/2., cutoutTubeLength[isub]/2.,
+                                        submllength);
+            }
+            tempSFoamlow1 = &( (*tempSFoamlow)<<submlpos);
           } else {
-            tempSFoamlow = new GeoTrd(foamthicknesslow/2., foamthicknesslow/2.,
-                                      cutoutTubeLength[isub]/2., cutoutTubeLength[isub]/2.,
-                                      submllength);
+            log << MSG::INFO << " problem with shape of lower foam trapezoid in LogVolName "
+                << logVolName << " thick,width,length " << foamthicknesslow
+                << " " << submlwidths << " " << submllength << endmsg;
           }
-          tempSFoamlow1 = &( (*tempSFoamlow)<<submlpos);
-        } else {
-          log << MSG::INFO << " problem with shape of lower foam trapezoid in LogVolName "
-              << logVolName << " thick,width,length " << foamthicknesslow
-              << " " << submlwidths << " " << submllength << endmsg;
         }
       }
 
@@ -272,12 +312,12 @@ GeoFullPhysVol* MultiLayer::build()
 
   // Add/subtract cylinders at ends of layers 2 and 4 to accommodate the last MDT
 
-  const GeoShape* stube = NULL;
+  const GeoShape* stube = nullptr;
   double tL = longWidth/2.0 - (tubePitch/2.)*TrdDwoverL;
   stube = new GeoTube(0.0, tubePitch/2., tL);
   stube = & ( (*stube) << GeoTrf::RotateX3D(90.*Gaudi::Units::deg) );
-  const GeoShape* stubewithcut = NULL;
-  if (cutoutNsteps > 1) {
+  const GeoShape* stubewithcut = nullptr;
+  if (cutoutNsteps > 1 && !m_nonCutoutXSteps.size()) { // adaption of tube cuts only needed for cutouts along amdb-y
     double toptubelength = cutoutTubeLength[cutoutNsteps-1];
     if (cutoutFullLength[cutoutNsteps-1]) toptubelength = longWidth;
     stubewithcut = new GeoTube(0.0, tubePitch/2., toptubelength/2.0);
@@ -301,7 +341,7 @@ GeoFullPhysVol* MultiLayer::build()
       slay = &(slay->subtract( (*stube)<<GeoTrf::Translate3D(-mdtthickness/2.+yy[i],0.,-length/2.) ));
       // add tube at the end
       // distinguish stations with/without cutouts
-      if (cutoutNsteps == 1) {
+      if (cutoutNsteps == 1 || m_nonCutoutXSteps.size()) {
         // no cutouts
         if (verbose_multilayer) {
           log << MSG::VERBOSE
@@ -309,7 +349,7 @@ GeoFullPhysVol* MultiLayer::build()
               << " z = " << length/2. << endmsg;
         }
         slay = &(slay->add( (*stube)<<GeoTrf::Translate3D(-mdtthickness/2.+yy[i],0.,length/2.-tubePitch/2.) ));
-      } else {
+      } else { // adaption of tube cuts only needed for cutouts along amdb-y
         // there are cutouts
         if (verbose_multilayer) {
           log << MSG::VERBOSE
@@ -330,10 +370,10 @@ GeoFullPhysVol* MultiLayer::build()
   GeoFullPhysVol* play = new GeoFullPhysVol(llay);
 
   double foamposition = 0.;
-  const GeoShape* sfoam = NULL;
-  const GeoMaterial* mfoam = NULL;
-  GeoLogVol* lfoam = NULL;
-  GeoPhysVol* pfoam = NULL;
+  const GeoShape* sfoam = nullptr;
+  const GeoMaterial* mfoam = nullptr;
+  GeoLogVol* lfoam = nullptr;
+  GeoPhysVol* pfoam = nullptr;
 
   if (foamthicknesslow != 0) {
     foamposition = -(mdtthickness - foamthicknesslow)/2.;
@@ -359,14 +399,16 @@ GeoFullPhysVol* MultiLayer::build()
     mfoam = getMaterialManager()->getMaterial("muo::Foam");
     lfoam = new GeoLogVol("MultiLayerFoam", sfoam, mfoam);
 
-  } else if (logVolName == "BME1MDT09" || logVolName == "BME2MDT09") {
-    lfoam = 0;
+  } else if ( logVolName.find("MDT09") != std::string::npos || logVolName.find("MDT14") != std::string::npos ) {
+    // MDT09 is used for BME and BIS sMDTs, MDT14 is used for 2nd multilayer of BIS78
+    // All those sMDTs have no foam layer
+    lfoam = nullptr;
   } else {
     log << MSG::ERROR << " no foam thickeness, while it was expected " << endmsg;
     throw std::runtime_error("ATTENTION:  no foam");
   }
 
-  if ( logVolName != "BME1MDT09" && logVolName != "BME2MDT09" ) {  //@@
+  if ( lfoam ) {  //@@
     pfoam = new GeoPhysVol(lfoam);
     GeoTransform* xf = new GeoTransform (GeoTrf::TranslateX3D(foamposition));
     GeoNameTag* nt = new GeoNameTag(name+" MultiLayerFoam");
@@ -387,7 +429,7 @@ GeoFullPhysVol* MultiLayer::build()
   std::vector<int> Ntubes;
 
   // No cutouts
-  if (cutoutNsteps <= 1) {
+  if (cutoutNsteps <= 1 || m_nonCutoutXSteps.size()) { // adaption of tube cuts only needed for cutouts along amdb-y
     for (int j = 0; j < nrOfSteps; j++) {
       tube.length=width+j*diff/nrOfSteps;
 
@@ -621,7 +663,7 @@ GeoFullPhysVol* MultiLayer::build()
       }
     } // Loop over layers
 
-  } else if (cutoutNsteps > 1) {
+  } else if (cutoutNsteps > 1 && !m_nonCutoutXSteps.size()) { // adaption of tube cuts only needed for cutouts along amdb-y
     bool arrowpointoutwards=false;
     bool cutAtAngle = cutoutAtAngle;
     if (xx[1]-xx[0]>0.) {
diff --git a/MuonSpectrometer/MuonGeoModel/src/MuonChamber.cxx b/MuonSpectrometer/MuonGeoModel/src/MuonChamber.cxx
index ae08912e1a3..090ba7f9c2d 100755
--- a/MuonSpectrometer/MuonGeoModel/src/MuonChamber.cxx
+++ b/MuonSpectrometer/MuonGeoModel/src/MuonChamber.cxx
@@ -747,8 +747,8 @@ GeoVPhysVol* MuonChamber::build(
       std::string key=stName+techname;
 
       // for cutouts:
-      // MDT cutouts for BOS1,5, BMS7,14, (problem with BMS4,10),  EMS
-      bool mdtCutoutFlag = ((stname == "BOS" && std::abs(zi) == 6) || stname == "BMG" ||
+      // MDT cutouts for BOS1,5, BMS7,14, (problem with BMS4,10), EMS, BMG and BIS MDT14
+      bool mdtCutoutFlag = ((stname == "BOS" && std::abs(zi) == 6) || stname == "BMG" || techname=="MDT14" ||
                             (stname == "BMS" && (std::abs(zi) == 1 && fi == 3)) ||
                             (stname == "EMS" && (std::abs(zi) == 1 || std::abs(zi) == 3)));
       if (((manager->IncludeCutoutsFlag() &&  mdtCutoutFlag) ||
@@ -1396,17 +1396,18 @@ GeoVPhysVol* MuonChamber::build(
       ndbz[doubletR-1]++;
 
       float zdivision=100.;// point between doubletZ=1 and 2;
-      if (stname.find("BI")!=std::string::npos && std::abs(zi)==7) zdivision=400.;//BIS78 : RPC8 is smaller than other RPCs
 
-      // do the following check for all RPCs but the BI (station1-6) ones
-      if (stname.find("BI")==std::string::npos || std::abs(zi)==7) {
+      // the BIS RPCs are 3-gap RPCs mounted inside of the BIS sMDTs, 
+      // for BIS78, there is a second RPC doubletZ at amdb-y (MuonGeoModel-z)=144mm inside the station
+      if (stname.find("BI")!=std::string::npos) {
+        if (std::abs(stationPhi)>=7 && rp->posz>100) doubletZ=2;
+        else doubletZ = ndbz[doubletR-1];
+      } else {
         if (zi <= 0 && !is_mirrored) {
           if (zpos < -zdivision*Gaudi::Units::mm) doubletZ=2;
         } else {
           if (zpos > zdivision*Gaudi::Units::mm) doubletZ=2;
         }
-      } else {
-        doubletZ = ndbz[doubletR-1];
       }
 
       // BMS (BOG) RPCs can have |xpos|=950 (|xpos|=350)
diff --git a/MuonSpectrometer/MuonGeoModel/src/sTGC.cxx b/MuonSpectrometer/MuonGeoModel/src/sTGC.cxx
index 6741762fb40..b7b6a3a92fc 100755
--- a/MuonSpectrometer/MuonGeoModel/src/sTGC.cxx
+++ b/MuonSpectrometer/MuonGeoModel/src/sTGC.cxx
@@ -3,7 +3,10 @@
 */
 
 #include "MuonGeoModel/sTGC.h"
+
 #include "MuonAGDDDescription/sTGC_Technology.h"
+#include "MuonAGDDDescription/sTGCDetectorDescription.h"
+#include "MuonAGDDDescription/sTGCDetectorHelper.h"
 #include "AGDDKernel/AGDDDetectorStore.h"
 #include "MuonGeoModel/Station.h"
 #include "MuonGeoModel/MYSQL.h"
@@ -21,7 +24,6 @@
 #include "GeoModelKernel/GeoTransform.h"
 #include "GeoModelKernel/GeoSerialIdentifier.h"
 #include "GeoModelKernel/GeoIdentifierTag.h"
-// for cutouts:
 #include "GeoModelKernel/GeoShapeSubtraction.h"
 #include "GeoModelKernel/GeoShapeIntersection.h"
 #include "GeoModelKernel/GeoShapeShift.h"
@@ -54,14 +56,17 @@ GeoFullPhysVol* sTGC::build(int minimalgeo)
 GeoFullPhysVol* sTGC::build(int minimalgeo, int , std::vector<Cutout*> )
 {
   AGDDDetectorStore* ds = AGDDDetectorStore::GetDetectorStore();
+  sTGCDetectorHelper stgcHelper;
+  sTGCDetectorDescription* stgc_descr = stgcHelper.Get_sTGCDetectorSubType(m_component->subType);
 
   sTGC_Technology* t = (sTGC_Technology*) ds->GetTechnology(name);
   thickness = t->Thickness();
   double gasTck=t->gasThickness;
   double pcbTck=t->pcbThickness;
-  double f4=t->f4Thickness;
-  double f5=t->f5Thickness;
-  double f6=t->f6Thickness;
+  double f4=stgc_descr->ylFrame();
+  double f5=stgc_descr->ysFrame();
+  double f6=stgc_descr->xFrame();
+
   //Evaluate honeycomb thickness
   double chamberTck = gasTck+pcbTck; //Note: pcbTck is the xml value and is the combined thickness of 2 pcbs.
   double honeycombTck = (thickness - 4*chamberTck)/5;
-- 
GitLab