#include <algorithm>
#include <iostream> 
#include <string>
#include <strstream>
#include <algorithm>
#include <cmath>

#include "DWUnits.hh"

#include "preXmlManager.hh"
#include "preXmlInteger.hh"
#include "preXmlAngle.hh"
#include "preXmlLength.hh"
#include "preXmlNumber.hh"
#include "preXmlNonNumber.hh"

#include "agddSection.hh"
#include "agddBox.hh"
#include "agddTrd.hh"
#include "agddTubs.hh"
#include "agddCons.hh"
#include "agddCompos.hh"
#include "agddUnion.hh"
#include "agddSubtraction.hh"
#include "agddMPosPhi.hh"
#include "agddPosXYZ.hh"
#include "agddAGDD.hh"
#include "agddNumberString.hh"

#include "agddUserGeometry.hh"

  // Some macros follow which save a bit of typing.
  // All they do is read values from the preXml file
  // and store them in constants and const vectors.
  // They make use of the fact that I have the same
  // names for my variables both in this program AND
  // in the preXml file.  You don't have to do it
  // this way if you don't want to.  By way of example,
  // I've typed the first "read" in it's expanded form.

#define DWSET(t,v)                                            \
  preXml:: t v ;                                              \
  useMy.fill( v , #v );                                       \
  cout << "I found "<< #v  <<" = "<< v << endl;

#define DWSETARRAY(t,v,n)                                     \
  vector< preXml:: t > v ( n ) ;                              \
  useMy.fill( v , #v ) ;                                      \
  {                                                           \
    for(int i=0; i< n ;i++) {                                 \
      cout << "I found "<< #v <<"["<<i<<"] = "<< v[i] <<endl; \
    }                                                         \
  };

// t : Type
// v : Name
// s : Structure
// n : Name

#define DWSTRUCTDEF(t,v)                                    \
  preXml::t v;

#define DWSTRUCTDEFARRAY(t,v,n)                             \
  vector< preXml:: t > v ;       

#define DWSTRUCTCFILL(t,v)                                    \
  useMy.fill( v , #v );                                       \
  cout << "I STRUCTC found " << #v  <<" = "<< v << endl;

#define DWSTRUCTCFILLARRAY(t,v,n)                            \
  v .resize( n );                                     \
  useMy.fill( v , #v );                                   \
  {                                                           \
    for(int i=0; i< n ;i++) {                             \
      cout << "I STRUCTC found " << #v <<"["<<i<<"] = "<< v [i] <<endl; \
    }                                                         \
  };

#define FORWARDINFOSTRUCT(mode)                                     \
    DWSTRUCT ## mode ## (Length,  dZWheel    )                      \
    DWSTRUCT ## mode ## (Length,  dZDisc     )                      \
    DWSTRUCT ## mode ## (Length,  rInt       )                      \
    DWSTRUCT ## mode ## (Length,  rInEnd     )                      \
    DWSTRUCT ## mode ## (Length,  rOuEnd     )                      \
    DWSTRUCT ## mode ## (Length,  zLoEnd     )                      \
    DWSTRUCT ## mode ## (Length,  zHiEnd     )                      \
    DWSTRUCT ## mode ## (Length,  thInsu     )                      \
    DWSTRUCT ## mode ## (Length,  thSF       )                      \
    DWSTRUCT ## mode ## (Length,  zTRTGap1   )                      \
    DWSTRUCT ## mode ## (Length,  zTRTGap2   )                      \
    DWSTRUCT ## mode ## (Length,  dZGap      )                      \
    DWSTRUCT ## mode ## (Length,  rBIP       )                      \
    DWSTRUCT ## mode ## (Length,  dZSC       )                      \
    DWSTRUCT ## mode ## (Length,  dZSub      )                      \
    DWSTRUCT ## mode ## (Length,  dZEle      )                      \
    DWSTRUCT ## mode ## (Length,  stagger    )                      \
    DWSTRUCT ## mode ## (Integer, radCabGran )                      \
    DWSTRUCT ## mode ## (Length,  rootAlArea )                      \
    DWSTRUCT ## mode ## (Integer, lastWheel  )                      \
                                                                    \
    DWSTRUCT ## mode ## (Integer, nRings )                          \
    DWSTRUCT ## mode ## ARRAY(Integer, nModul,   nRings )           \
    DWSTRUCT ## mode ## ARRAY(Angle,   dFMod,    nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  rISi,     nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  rOSi,     nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  rIEle,    nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  rOEle,    nRings )           \
    DWSTRUCT ## mode ## ARRAY(Angle,   stereo,   nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  dRCool,   nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  dQCool,   nRings )  	    \
    DWSTRUCT ## mode ## ARRAY(Length,  dZCool,   nRings )           \
    DWSTRUCT ## mode ## ARRAY(Length,  dQSub,    nRings )           \
    DWSTRUCT ## mode ## ARRAY(Integer, stagSign, nRings )           \
                                                                    \
    DWSTRUCT ## mode ## (Integer, nWheels )                         \
    DWSTRUCT ## mode ## ARRAY(Length,    zWheel,         nWheels )  \
    DWSTRUCT ## mode ## ARRAY(Angle,     dFWheel,        nWheels )  \
    DWSTRUCT ## mode ## ARRAY(Length,    rIDisc,         nWheels )  \
    DWSTRUCT ## mode ## ARRAY(Integer,   sgnSto,         nWheels )  \
    DWSTRUCT ## mode ## ARRAY(NonNumber, freeFormLiRing, nWheels )

#define DEFFORWARDINFOCLASS        \
      FORWARDINFOSTRUCT(DEF)
      
#define FILLFORWARDINFOCLASS  \
      FORWARDINFOSTRUCT(CFILL)
      
class ForwardInfo {
public:
 typedef vector<int> vectorOfIntegers;
  DEFFORWARDINFOCLASS
  vector<vectorOfIntegers> liRing;
public:
  ForwardInfo(preXml::Manager & useMy) {
    FILLFORWARDINFOCLASS
    {
	for (int i=0; i<nWheels; i++) {
	    const string freeForm = freeFormLiRing[i];
	    istrstream is(freeForm.c_str());
	    int tmp;
	    vector<int> theRingTypes;
	    while (is >> tmp) {
		theRingTypes.push_back(tmp-1);  // Note the -1.  I adjust so that 
		                                // the rings are numbered 0,1,2,... 
	    };
	    liRing.push_back(theRingTypes);
	};
    };
  }; // end of ForwardInfo constructor.

  inline const double alArea() const { return (rootAlArea*rootAlArea); };

  inline const bool electronicsIsOnOuterEdge(const int ring) const {
    return (rIEle[ring]>rISi[ring]);
  };
  
  // Average radius of electronics
  inline const double eAv(const int ring) const {
    return 0.5*(rOEle[ring]+rIEle[ring]);
  };
  
  inline const double ddzz(const int ring) const {
    return dZCool[ring]*dQCool[ring]*nModul[ring]/(eAv(ring)*2.0*DWpi);
  };
  
  inline const double dR(const int ring) const {
    return (rInt-rIEle[ring])/((double)radCabGran);
  };
  
  // Forward electronics length
  inline const double dQEle(const int ring) const {
    return tan((dFMod[ring]/2.0)/DWradian) * 
      ( 
       electronicsIsOnOuterEdge(ring) ?
       rIEle[ring] :
       rOEle[ring] 
       ); 
  };
  
  // Forward detector thickness
  inline const double tDet() const { return max(dZSC, dZEle); }

  inline const int noRing(const unsigned int wheel) const {
    return liRing[wheel].size();
  };
  
  // Stereo angle for tilted wafer in given forward ring and wheel.
  inline const double fi(const int ring, 
			 const int wheel) const {
    return sgnSto[wheel]*stereo[ring];
  };

  inline const double zCEnd()   const { return 0.5*(zHiEnd+zLoEnd); };
  inline const double zHEnd()   const { return 0.5*(zHiEnd-zLoEnd); };
  inline const double t2()      const { return 0.5*thInsu;          };
  inline const double zCentre() const { return zCEnd() - t2();      };
  inline const double zHLen()   const { return zHEnd() - t2();      };
};

agdd::Volume * generalZALP(const double rIn,
			   const double rOut,
			   const double length,
			   const string nameEnd) {
  agdd::Tubs * thisZALP = new agdd::Tubs;
  thisZALP->setName("sct_fwd_ZALP"+nameEnd);
  thisZALP->setTubs(rIn,
		    rOut,
		    length);
  thisZALP->setMaterial("Aluminium");
  return thisZALP;
};

// Cooling pipe of Al and water
agdd::Volume * generalZBIP(const double rIn,
			   const double rOut,
			   const double length,
			   const string nameEnd) {
  agdd::Tubs * thisZALP = new agdd::Tubs;
  thisZALP->setName("sct_fwd_ZBIP"+nameEnd);
  thisZALP->setTubs(rIn,
		    rOut,
		    length);
  thisZALP->setMaterial("Cooling");
  return thisZALP;
};

// Outer support frame
agdd::Volume * theZOSF(const double rMin,
		       const double rMax,
		       const double length) {
  agdd::Tubs * zosf = new agdd::Tubs;
  zosf->setName("sct_fwd_outerSupportFrame");
  zosf->setMaterial("Cfibre");
  zosf->setTubs(rMin, rMax, length);
  return zosf;
};

// ZOCP is the region occupied by the cooling pipes.
// Assume one cooling circuit per quadrant of each ring,
// i.e. eight pipes per ring.
agdd::Volume * theZOCP(const double rOSF,
		       /* const double rOCP, */
		       const ForwardInfo & fIS) {
  agdd::Compos * zocp = new agdd::Compos;
  zocp->setName("sct_fwd_coolingPipesZOCP");

  {
    const double area = 8.0 * DWpi * fIS.rBIP * fIS.rBIP;
    double rMin = rOSF + fIS.rBIP;
    for (int wheel=0; wheel<fIS.nWheels; wheel++) {
      const double wheelStartingAtOne = wheel+1;
      const double zGap = 
	(wheelStartingAtOne<=fIS.lastWheel) ? fIS.zTRTGap1 : fIS.zTRTGap2;
      const double halfLength=0.5*(abs(fIS.zWheel[wheel]-zGap));
      const double length = 2.0*halfLength;
      const double zOffset = 0.5*(fIS.zWheel[wheel]+zGap)-fIS.zCentre();
      const double thickness = fIS.noRing(wheel)*area/(DWtwoPi*rMin);
      const double rMax = rMin+thickness;
      agdd::PosXYZ * pos = new agdd::PosXYZ;
      {
	ostrstream os; os << "ForWheel" << wheel+1 
			  << "ForZOCP" << ends;
	pos->setVolume(generalZBIP(rMin, rMax, length, os.str()));
      };
      pos->setXYZ(0, 0, zOffset);

      zocp->addPos(pos);

      rMin+=thickness;
    };
  };
  return zocp;
};

agdd::Volume * theZOCA(const double rOCP,
		       const ForwardInfo & fIS) {
  agdd::Compos * zoca = new agdd::Compos;
  {
    ostrstream os; os << "sct_fwd_cablesZOCA" << ends;
    zoca->setName(os.str());
  };
  const double rStep = (fIS.rOuEnd-fIS.thInsu-rOCP)/(fIS.nWheels+2.0);
  for (int wheel=0; wheel<fIS.nWheels; wheel++) {
    const double wheelStartingAtOne = wheel+1;
    const double zGap = 
      (wheelStartingAtOne<=fIS.lastWheel) ? fIS.zTRTGap1 : fIS.zTRTGap2;
    const double halfLength = 0.5*abs(fIS.zWheel[wheel]-zGap);
    const double length = 2.0*halfLength;
    const double zOffset = 0.5*(fIS.zWheel[wheel]+zGap)-fIS.zCentre();
    int uMods=0;
    for (int n=0; n<fIS.noRing(wheel); n++) {
      const unsigned int ring=fIS.liRing[wheel][n]; // ring means `ringType'
      uMods += fIS.nModul[ring]; 
    };
    const double area = fIS.alArea()*uMods;
    const double rMin = rOCP+wheelStartingAtOne*rStep;
    const double thickness = area/(DWtwoPi*rMin);
    const double rMax = rMin+thickness;
    ostrstream os; os << "ForZOCAForWheel" << wheel+1 << ends;
    const string nameEnd = os.str();
    {
      agdd::PosXYZ * pos = new agdd::PosXYZ;
      pos->setVolume(generalZALP(rMin, rMax, length, nameEnd));
      pos->setXYZ(0, 0, zOffset);
      
      zoca->addPos(pos);
    };
  };
  return zoca;
};

// The radial space between the SCT wheels and the TRT:
agdd::Volume * theZOUT(const ForwardInfo & fIS) {
  agdd::Compos * zout = new agdd::Compos;
  zout->setName("sct_fwd_ZOUT");
  const double rMin       = fIS.rInt;
  // Not needed: const double rMax       = fIS.rOuEnd-fIS.thInsu;
  const double halfLength = fIS.zHLen()-2.0*fIS.t2();
  const double length     = 2.0*halfLength;
  const double rOSF       = fIS.rInt+fIS.thSF;
  const double rOCP       = (fIS.rOuEnd+rOSF-fIS.thInsu)/2.0;
  {
    agdd::PosXYZ * zosfPos = new agdd::PosXYZ;
    zosfPos->setVolume(theZOSF(rMin, rOSF, length));
    zout->addPos(zosfPos);
  };
  {
    agdd::PosXYZ * zocpPos = new agdd::PosXYZ;
    zocpPos->setVolume(theZOCP(rOSF, /* rOCP, */ fIS));
    zout->addPos(zocpPos);
  };
  {
    agdd::PosXYZ * zocaPos = new agdd::PosXYZ;
    zocaPos->setVolume(theZOCA(rOCP, fIS));
    zout->addPos(zocaPos);
  };
  return zout;
};

// Here is the sensitive silicon.
// Although actually built from two pieces of silicon, the
// sensitive material is described as just one piece here.
// We build it with it's length along X, and it's flat faces
// pointing along +/- Z.  The face at high X is the larger of
// the two unequal faces.  The origin is in the phys. centre.
agdd::Volume * ringSpecificZSEN(const int ring,
				const double thickness,
				const double innerWidth,
				const double outerWidth,
				const double length, 
				const ForwardInfo & fIS
				) {

  // The folowing precautions are now unnecessary, since this 
  // checking is now taken care of in ZDET.
  static bool firstTime=true;
  static vector < agdd::Volume * > senVec;
  if (firstTime) {
    firstTime=false;
    for (int i=0; i<fIS.nRings; i++) {
      senVec.push_back(NULL);
    };
  };

  if (!senVec[ring]) {
    // Must make this rZSEN
    agdd::Compos * thisZSEN = new agdd::Compos;
    
    {
      ostrstream os; os << "sct_fwd_betterOrientedSiliconForRing" << ring+1 << ends;
      thisZSEN->setName(os.str());
    }
    
    {
      agdd::Trd * thisSilicon = new agdd::Trd;
      thisSilicon->setTrd(thickness,  thickness,
			  innerWidth, outerWidth,
			  length);
      thisSilicon->setMaterial("silicon");
      {
	ostrstream os; os << "sct_fwd_siliconForRing" << ring+1 << ends;
	thisSilicon->setName(os.str());
      }    
      {
	agdd::PosXYZ * rot = new agdd::PosXYZ;
	rot->setRot(0, 90.0*DWdegree, 0);
	rot->setVolume(thisSilicon);
	thisZSEN->addPos(rot);
      }
    }
    senVec[ring]=thisZSEN;
  };
  return senVec[ring];
};

// Here is the electronics, oriented to be consistent with the ZSEN
agdd::Volume * ringSpecificZELE(const int ring,
		    const double thickness,
		    const double width,
		    const double length) {
  agdd::Box * zele = new agdd::Box;
  {
    ostrstream os; os << "sct_fwd_detectorElectronicsForRing" 
		   << ring+1 << ends;
    zele->setName(os.str());
  }
  zele->setBox(width, length, thickness);
  zele->setMaterial("hybridSilicon");
  return zele;
}; 

agdd::Volume * ringSpecificZSUB(const int ring,
				const ForwardInfo & fIS) {
  // A Substrate for forward modules.
  // The origin for this piece will be at the physics centre.
  static bool firstTime = true;
  static vector < agdd::Volume * > subVec;
  if (firstTime) {
    firstTime=false;
    for (int i=0; i<fIS.nRings; i++) {
      subVec.push_back(NULL);
    };
  };

  if (!subVec[ring]) {
    // We had better make this detector
    agdd::Box * box = new agdd::Box;
    {
      ostrstream os; os << "sct_fwd_substrateForRing" << ring+1 << ends;
      box->setName(os.str());
    };
    box->setMaterial("Beryllia");
    {
      const double length = 
	(fIS.electronicsIsOnOuterEdge(ring)) ?
	fIS.rOEle[ring]-fIS.rISi[ring] :
	fIS.rOSi[ring]-fIS.rIEle[ring] ;
      const double width = fIS.dQSub[ring];
      const double thickness = fIS.dZSub;
      box->setBox(length, width, thickness);
    };

    agdd::Compos * holder = new agdd::Compos;
    {
      ostrstream os; 
      os << "sct_fwd_holderForSubstrateForRing" 
	 << ring+1 << ends;
      holder->setName(os.str());
    };
    { 
      agdd::PosXYZ * pos    = new agdd::PosXYZ;
      pos->setVolume(box);
      const double elecHalfWidth=0.5*(fIS.rOEle[ring]-fIS.rIEle[ring]);
      if (fIS.electronicsIsOnOuterEdge(ring)) {
	pos->setXYZ(+elecHalfWidth, 0, 0);
      } else {
	pos->setXYZ(-elecHalfWidth, 0, 0);
      };
      holder->addPos(pos);
    };
    
    subVec[ring]=holder;
  }; // We made it!

  return subVec[ring];
};

agdd::Volume * ringSpecificZDET(const int ring,
				const ForwardInfo & fIS) {
  // ZMOD is a module.  They come in four type indexed by "ring".
  static bool firstTime = true;
  static vector < agdd::Volume * > detVec;
  if (firstTime) {
    firstTime=false;
    for (int i=0; i<fIS.nRings; i++) {
      detVec.push_back(NULL);
    };
  };

  if (!detVec[ring]) {
    // We had better make this detector
    const double siliconHalfThickness  =fIS.dZSC; // of semiconductor
    const double siliconOuterHalfWidth =tan((fIS.dFMod[ring]/2.0)/DWradian)*fIS.rOSi[ring];
    const double siliconInnerHalfWidth =tan((fIS.dFMod[ring]/2.0)/DWradian)*fIS.rISi[ring];
    const double siliconHalfLength     =(fIS.rOSi[ring]-fIS.rISi[ring])/2.0;
    
    const double siliconThickness  = 2.0*siliconHalfThickness;
    const double siliconOuterWidth = 2.0*siliconOuterHalfWidth;
    const double siliconInnerWidth = 2.0*siliconInnerHalfWidth;
    const double siliconLength     = 2.0*siliconHalfLength;
    
    const double elecHalfThickness = fIS.dZEle;
    const double elecHalfLength    = fIS.dQEle(ring);
    const double elecHalfWidth     = (fIS.rOEle[ring]-fIS.rIEle[ring])*0.5;
    
    const double elecThickness = 2.0 * elecHalfThickness;
    const double elecLength    = 2.0 * elecHalfLength;
    const double elecWidth     = 2.0 * elecHalfWidth;
    
    const double centreSeparation = siliconHalfLength+elecHalfWidth;
    
    agdd::Compos * thisZDET = new agdd::Compos;
    
    {
      ostrstream os; os << "sct_fwd_detectorForRing" << ring+1 << ends;	
      thisZDET->setName(os.str());
    }
    
    { 
      // Here comes the silicon:
      agdd::PosXYZ * siliconPos = new agdd::PosXYZ;
      siliconPos->setVolume(ringSpecificZSEN(ring, 
					     siliconThickness,
					     siliconInnerWidth,
					     siliconOuterWidth,
					     siliconLength,
					     fIS));
      thisZDET->addPos(siliconPos);
    }
    {
      // Here comes the electronics:
      agdd::PosXYZ * elecPos = new agdd::PosXYZ;
      elecPos->setVolume(ringSpecificZELE(ring,
					  elecThickness,
					  elecWidth,
					  elecLength)
			 );
      const double elecShift 
	= fIS.electronicsIsOnOuterEdge(ring) ? 
	+centreSeparation :
	-centreSeparation;
      
      elecPos->setXYZ(elecShift, 0, 0);
      thisZDET->addPos(elecPos);
    }
    
    detVec[ring]=thisZDET;
  } // We finished making it.
  return detVec[ring];
};  

agdd::Volume * ringSpecificZMOD(const int ring,
				const int wheel,
				const ForwardInfo & fIS) {
  // ZMOD is a module.  They come in four type indexed by "ring".
  
  agdd::Compos * thisZMOD = new agdd::Compos;
  
  {
    ostrstream os; os << "sct_fwd_moduleForRing" << ring+1 
		      << "Wheel" << wheel+1 << ends;	
    thisZMOD->setName(os.str());
  }
  
  { 
    const double vertOffset=fIS.tDet()+fIS.dZSub;
    agdd::Volume * thisRingSpecificZDET = ringSpecificZDET(ring, fIS);
    agdd::PosXYZ * detector1Pos = new agdd::PosXYZ;
    agdd::PosXYZ * detector2Pos = new agdd::PosXYZ;
    detector1Pos->setVolume(thisRingSpecificZDET);
    detector2Pos->setVolume(thisRingSpecificZDET);
    detector1Pos->setXYZ(0,0,+vertOffset);
    detector2Pos->setXYZ(0,0,-vertOffset);
    detector1Pos->setRot(0,0,0);
    detector2Pos->setRot(0,0,fIS.fi(ring, wheel));
    thisZMOD->addPos(detector1Pos);
    thisZMOD->addPos(detector2Pos);
  }

  {
    agdd::PosXYZ * substratePos = new agdd::PosXYZ;
    substratePos->setVolume(ringSpecificZSUB(ring, fIS));
    thisZMOD->addPos(substratePos);
  }

  return thisZMOD;
};  

agdd::Volume * ringSpecificZRII(const int ring,
				const int wheel,
				const double rIn,
				const double rOut,
				const ForwardInfo & fIS) {
  // ZRII is a ring of alternate modules.  
  // They come in up to ring*wheel (potentially) different types
  // indexed by "ring" and "wheel".
  
  // The following is no longer required now that ZRIIs are 
  // discovered to be wheel dependent
  /*
  static bool firstTime=true;
  static vector< agdd::Volume * > ringVec; // one for each type of ring
  if (firstTime) {
    firstTime=false;
    for (int i=0; i<fIS.nRings; i++) {
      ringVec.push_back(NULL);
    };
  };
  if (ringVec[ring]) {
    // I already made that ring of modules, so do nothing
  } else {
  */
  // I must make it
  agdd::Compos * thisZRII = new agdd::Compos;
  {
    ostrstream os; os << "sct_fwd_ZRIIForRing" << ring+1 
		      << "Wheel" << wheel+1 <<  ends;	
    thisZRII->setName(os.str());
  }
  {
    agdd::MPosPhi * mPos = new agdd::MPosPhi;
    thisZRII->addPos(mPos);
    { // Set N
      const unsigned int n = fIS.nModul[ring]; // Number of modules in a WHOLE ring
      const unsigned int halfN = n/2; // Number of modules in one of these ALTERNATE rings (every other module).
      if (2*halfN!=n) {
	cerr << "Error: Number of modules="<< n 
	     <<" in ring "<<ring+1<<" is not even!"<<endl;
      };
      assert(2*halfN==n);
      mPos->setN(halfN); 
    } // done setting N
    mPos->setRZ(0.5*(rIn+rOut), 0.0);
    mPos->setVolume(ringSpecificZMOD(ring, wheel, fIS)); 
  };
  /*
    ringVec[ring]=thisZRII;
    };
    return ringVec[ring];
  */
  return thisZRII;
};  

agdd::Volume * ringSpecificZBIPForTheZWSU(const int ring,
				       const ForwardInfo & fIS) {
  static bool firstTime = true;
  static vector < agdd::Volume * > zbipVec;
  if (firstTime) {
    firstTime=false;
    for (int i=0; i<fIS.nRings; i++) {
      zbipVec.push_back(NULL);
    };
  };
  // If we have not yet made this ZBIP, we must:
  if (!zbipVec[ring]) {
    const double rMin = fIS.eAv(ring)-fIS.rBIP;
    const double rMax = fIS.eAv(ring)+fIS.rBIP;
    const double halfLength = fIS.rBIP; // binary ice pipe?
    const double length = 2.0 * halfLength;

    ostrstream os; os << "ForZWSUForRing" << ring+1 << ends;
    const string nameEnd = os.str();

    zbipVec[ring]=generalZBIP(rMin, rMax, length, nameEnd);
  };
  // Now we can return it:
  return zbipVec[ring];
};

agdd::Volume * ringSpecificZALPForTheZWEI(const int ring,
					  const ForwardInfo & fIS) {
    static bool firstTime=true;
    static vector< agdd::Volume * > ringVec;
    if (firstTime) {
	firstTime=false;
	for (int i=0; i<fIS.nRings; i++) {
	    ringVec.push_back(NULL);
	};
    };
    if (ringVec[ring]) {
      // I already made that ring, so do nothing
    } else {
      // I must make it
      
      ostrstream os; os << ring+1<<"ForTheZWEI" << ends;	
      ringVec[ring]=generalZALP(fIS.eAv(ring)-0.5*fIS.dRCool[ring],
				fIS.eAv(ring)+0.5*fIS.dRCool[ring],
				2.0*fIS.ddzz(ring),
				os.str());
    };
    return ringVec[ring];
};

agdd::Volume * radialCableZALP(const int ring,
			       const ForwardInfo & fIS) {
    // Now we make a few ring-specific ZALPs which represent cables.
    static bool firstTime=true;
    static vector< agdd::Volume * > ringVec;
    if (firstTime) {
	firstTime=false;
	for (int i=0; i<fIS.nRings; i++) {
	    ringVec.push_back(NULL);
	};
    };
    if (ringVec[ring]) {
	// I already made that ring, so do nothing
    } else {
	// I must make it    
      //const double alArea  = fIS.rootAlArea*fIS.rootAlArea;
	agdd::Compos * thisRadialBit = new agdd::Compos;
	//radialCableZALP.push_back(thisRadialBit);
	{
	    ostrstream os; os << "sct_fwd_radialCableZALPForRing"
			      << ring+1 << ends;
	    thisRadialBit->setName(os.str());
	}
	for (int m=0; m<fIS.radCabGran; m++) {
	    const double rr=fIS.rIEle[ring]+((double)m)*fIS.dR(ring);
	    
	    {
	      ostrstream os; os << "ForRadialCableEffectRing" 
				<< ring+1 
				<< "Grain" << m+1 << "g" << ends;	
	      agdd::PosXYZ * thisZALPPos = new agdd::PosXYZ;
	      thisZALPPos->setVolume(generalZALP(rr,
						 rr+fIS.dR(ring),
						 2.0*(fIS.alArea()*fIS.nModul[ring]/(4.0*DWpi*rr)),
						 os.str())
				     );
	      
	      thisRadialBit->addPos(thisZALPPos);
	    }
	};
	ringVec[ring]=thisRadialBit;
	// That's all the ring speciific ZALPs which represent cables made.
    };
    return ringVec[ring];
};

// A ZWSU is a set of circular pipes which contact the 
// cooling blocks in a wheel.
agdd::Volume * wheelSpecificZWSU(const int wheel,
				 /* const vector< vector<int> > & liRing, */
				 const ForwardInfo & fIS) {
  agdd::Compos * zusu = new agdd::Compos;
  {
    ostrstream os; os << "sct_fwd_ZWSUForWheel" << wheel+1 << "g" << ends;
    zusu->setName(os.str());
  };
  for (int n=0; n<fIS.noRing(wheel); n++) {
    const unsigned int ring=fIS.liRing[wheel][n]; // ring means `ringType'
    agdd::PosXYZ * pos = new agdd::PosXYZ;
    pos->setVolume(ringSpecificZBIPForTheZWSU(ring,fIS));
    const double staggerSign = fIS.stagSign[ring];
    const double offset = (fIS.dZDisc-fIS.rBIP)*staggerSign;
    pos->setXYZ(0,0,offset);
    zusu->addPos(pos);
  };
  return zusu;
};

// One forward SCT wheel, known as a ZWEI in AGE.
agdd::Volume * newWheel(const int wheel, 
			/* const vector< vector<int> > & liRing, */
			const ForwardInfo & fIS ) {
    agdd::Compos * ZWEI = new agdd::Compos;
    { // Define the ZWEI
	{ // Name the ZWEI
	    ostrstream os;  os << "sct_fwd_ZWEIForWheel" << wheel+1 << ends; 
	    ZWEI->setName(os.str());
	}
	{ // Put in things only dependent on wheel first:
	  agdd::PosXYZ * zwsuPos = new agdd::PosXYZ;
	  zwsuPos -> setVolume(wheelSpecificZWSU(wheel, /* liRing,*/ fIS));
	  ZWEI->addPos(zwsuPos);
	};
	{ // put in ring specific things:
	  for (int n=0; n<fIS.noRing(wheel); 
	       n++) { // previously noRing()
	    const unsigned int ring=fIS.liRing[wheel][n]; // ring means `ringType'
	    
	    const double rIThi= 2.0*fIS.tDet()+fIS.dZSub;
	    const double staggerSign = fIS.stagSign[ring]; //(nIsEven)?-1.0:+1.0;
	    const double zSRing = staggerSign;

	    { // put in cables and cooling
	      
	      agdd::PosXYZ * ZALPPos = new agdd::PosXYZ;
	      ZALPPos->setVolume(ringSpecificZALPForTheZWEI(ring,fIS));
	      ZALPPos->setXYZ(0.0, 0.0, (fIS.ddzz(ring)+fIS.dZDisc)*(staggerSign));
	      ZWEI->addPos(ZALPPos);
	      const double xx = 
		staggerSign*(+2.0*rIThi //[ring]
			     +2.0*fIS.dZCool[ring]
			     +fIS.stagger
			     +fIS.dZDisc
			     +((double)(n+1))*0.1*DWcm);  
	      // I don't know what that 0.1*DWcm 
	      // is for, but it's in AGE.
	      
	      agdd::PosXYZ * radialCableZALPPos = new agdd::PosXYZ;
	      radialCableZALPPos->setVolume(radialCableZALP(ring,fIS));
	      radialCableZALPPos->setXYZ(0.0, 0.0, xx);
	      ZWEI->addPos(radialCableZALPPos);

	    } // That's all the cables and cooling
	    
	    { // Put in the detectors
	      
	      const bool electronicsOuter=
		fIS.electronicsIsOnOuterEdge(ring);
	      // (as opposed to on the inner edge)
	      
	      // AGE places its module origins at the centre of electronics AND silicon, whereas I put mine at the physics centreconst, consequently I do not use the following two lines, which would be right for AGE, but use the two lines after that instead.
	      //const double rIn =(electronicsIsOnOuterEdge ? fIS.rISi[ring] : fIS.rIEle[ring] );
	      //const double rOut=(electronicsIsOnOuterEdge ? fIS.rOEle[ring] : fIS.rOSi[ring] );
	      const double rIn = fIS.rISi[ring];
	      const double rOut= fIS.rOSi[ring];
	      const double dFDet=(electronicsOuter ?
				  fIS.dFMod[ring]/2.0 : atan2(fIS.dQEle(ring), fIS.rIEle[ring])*DWradian );
	      const double rR=0.5*(fIS.rOSi[ring]+fIS.rISi[ring]);
	      //const double eAv=fIS.eAv(ring);
	      /* const double fi=fIS.sgnSto[wheel]*fIS.stereo[ring];*/
	      const double xX=+rR*(1.0-cos(fIS.fi(ring,wheel)/DWradian));
	      const double yY=-rR*sin(fIS.fi(ring,wheel)/DWradian);

	      agdd::Volume * theZRII = ringSpecificZRII(ring, wheel, 
							rIn, rOut, 
							fIS);
	      for (unsigned int iStag=1; iStag<=3; ++++iStag) { 
		// First we put in one of the "alternate rings" (iStag=1) 
		// and then the other (iStag=3).

		const double fiFi=((iStag==1)?
				   0:
				   360.0*DWdegree/fIS.nModul[ring]);
		const double zPos=zSRing*(fIS.dZDisc
					  +2.0*fIS.dZCool[ring]
					  +rIThi
					  +(iStag-2.0)*fIS.stagger
					  );

		agdd::PosXYZ * ZRIIPos = new agdd::PosXYZ;
		ZRIIPos->setVolume(theZRII); 
		ZRIIPos->setXYZ(0.0, 0.0, zPos);
		ZRIIPos->setRot(0, 0, fiFi);
		ZWEI->addPos(ZRIIPos);

	      } 

	    } // That's done the detectors

	  } // end of loop over rings
	} // That's all the ring specific things
    } // Finished defining the ZWEI
    return ZWEI;
};

typedef agdd::UserGeometry::OutputStructure OutputStruct;
OutputStruct agdd::UserGeometry::defineUserGeometry() {
  preXml::Manager & useMy = *(_in.preXmlManager);
  agdd::UserGeometry::OutputStructure outStruct;
  
  // At the top of the preXml file, the dtd insists on the existence
  // of at least one "geometry" (and possibly more, nested if necessary).
  // The "geometries" are very much like directories.  They permit
  // the user to store more than one description of the detector in 
  // one preXml file.
  //
  // So, before we can proceed, we must tell the preXml::Manager which 
  // of these "geometries" (there may of course only be one!) we wish
  // to build the output-XML-file from.
  // This we do in the same was as "changing directories" on a file-
  // system.  When we get to the geometry we want, we tell the
  // preXml::Manager to setGeometry(), and that's all there is to it.

  // Here we choose the geometry to work on.
  useMy.geoManager().changeDir(string("G3Geometry"));
  useMy.setGeometry();
  
  //   Now that we've chosen the geometry, we have restricted ourselves
  // to a part of the preXml file that describes just one detector.
  //
  //   If we were confident enough (I'm not!) that nobody working on
  // a different sub-detector to us had created a name-clash for us,
  // then we could go right ahead and start reading data from the 
  // preXml file immediately.
  //   If we were to follow this "go right ahead" method, the 
  // preXml::Manager would spot any name clashes and tell us about them,
  // so nothing would go wrong. However, having to fixing the clash 
  // would still involve finding the person from the other sub-detector
  // responsible for the clash, etc, etc, etc.  This could take time. 
  //   A much better solution is for us to use the structure of the
  // "sections" (again used just like directories) which each geometry
  // in the preXml file is broken down into.
  //   If I "change directory" to  the section that is my sole  
  // responsibility, then if there is any problem, I only have to
  // grumble to myself.  Since I'm working on the SCT barrel, 
  // I'll go there:
  
  assert(useMy.secManager().changeDir(string("ATLAS")));
  assert(useMy.secManager().changeDir(string("InnerDetector")));
  assert(useMy.secManager().changeDir(string("SCT")));
  // I asserted these just to check that I succeeded.

  // Now we start the geometry proper.

  // Let's use this opportunity to define what 
  // we are trying to build.
  agdd::AGDD * theAGDD=new agdd::AGDD;
  outStruct.theBuildable=theAGDD;
  
  { // START OF BARREL
    assert(useMy.secManager().changeDir(string("Barrel")));

    // We're going to read those parameters which interest
    // us from the preXml file, and then produce the 
    // relevant output.
    // The xml file comes to us as the element called doc.

    // Since there are potential name clashes everywhere, let's
    // restrict ourselves to looking in only that part of the preXml file
    // which is relevant to the SCT barrel and this particular geometry.

    // -----------------------------------------------------
    // -------- Some SCT-Barrel specific quantites ---------
    // -------- (close to old AGE structure)       --------- 
    // -----------------------------------------------------

    // The first MYSET which follows would look like this if expanded:
    //
    // preXml::Length rMin;
    // useMy.fill(rMin, "rMin");
    // cout << "I found rMin = " << rMin << endl;

    DWSET(Length,  rMin)     //  overall min radius
    DWSET(Length,  rMax)     //  overall max radius
    DWSET(Length,  hLen)     //  overall half length
    DWSET(Integer, nLay)     //  Number of Layers
    DWSET(Length,  dYSect)   //  Width of a ski
    DWSET(Number,  rTckEle)  //  Thickness of the Electronics (in rad.lens.)
    DWSET(Number,  rTckLum)  //  THICKNESS OF THE LUMP (in % of r.l.)
    DWSET(Length,  wdthLum)  //  WIDTH OF THE LUMPED MATERIAL
    DWSET(Length,  zOvlp)    //  min Z OVERLAP +\- 2 sigma from origin
    //Not used. DWSET(Length, guard) // width of Guard Ring at crystal border
    DWSET(Length,  res)      //  width of polySi resistors at z ends of modules
    DWSET(Length,  elecZ)    //  The Board length
    DWSET(Length,  sep)      //  radial sep between modules adjacent in z
    DWSET(Length,  inn1)     //  INNER RADIUS OF CABLE 1
    DWSET(Length,  out1)     //  OUTER RADIUS OF CABLE 1
    DWSET(Length,  inn2)     //  INNER RADIUS OF CABLE 2
    DWSET(Length,  out2)     //  OUTER RADIUS OF CABLE 2
    DWSET(Number,  thInsu)   //  Insulator thickness (in % of r.l.)
  
    //----------------------------------------------------
    //---------- Some extra/redefined quantities ---------
    //----------------------------------------------------

    DWSET(Angle,   stereoInner)     // (as viewed from I.P.)
    DWSET(Angle,   stereoOuter)     // (as viewed from I.P.)
    DWSET(Length,  waferThickness)
    DWSET(Length,  waferWidth)   
    DWSET(Length,  waferLength)
    DWSET(Length,  sectGap)       
    DWSET(Length,  waferSep)     
    DWSET(Length,  radEle)   
    DWSET(Length,  radIso)  
    DWSET(Length,  radSup) 
    DWSET(Length,  radLum) 
    DWSET(Integer, nModulesPerSki)

    DWSETARRAY(Length, zModuleOffset, nModulesPerSki)
  
    //----------------------------------------------------
    //---------- Some layer dependant quantities ---------
    //----------------------------------------------------

    DWSETARRAY(Length,  layerRadius,  nLay)  
    DWSETARRAY(Number,  rTckSup,      nLay)  
    DWSETARRAY(Angle,   tiltOfLayer,  nLay)  
    DWSETARRAY(Integer, skisPerLayer, nLay)  
    DWSETARRAY(Length,  r0Sup,        nLay)  
    DWSETARRAY(Length,  flLen,        nLay)  
    DWSETARRAY(Number,  flTck,        nLay)  
    DWSETARRAY(Length,  cOLen,        nLay)  
    DWSETARRAY(Number,  cOTck,        nLay)  

    //----------------------------------------------------------
    //---------- Now let's derive some quantities.  ------------
    //---------- We don't HAVE to do this here;     ------------
    //---------- some params are better defined     ------------
    //---------- further on inside loops.           ------------
    //----------------------------------------------------------

    const double boardWidth      = dYSect-waferWidth;
    const double tckEle          = (rTckEle/DWfraction) * radEle;
    const double tckLum          = (rTckLum/DWfraction) * radLum;
    const double tckIso          = (thInsu /DWfraction) * radIso;
    const double sectThickness   = max((double) waferThickness, tckEle);
    const double sectWidth       = waferWidth+2.0*boardWidth;
    const double sectLength      = max(waferLength*2.0+waferSep, (double) elecZ);
    const double maxStereo       = max(abs(stereoInner), abs(stereoOuter));
    const double moduleThickness = 2.0*sectThickness+sectGap;
    const double moduleWidth     = sectWidth*cos(maxStereo/DWradian)+
                                 sectLength*sin(maxStereo/DWradian);
    const double moduleLength    = sectWidth*sin(maxStereo/DWradian)+
                                 sectLength*cos(maxStereo/DWradian);
    const double zCylEnd         = hLen - tckIso;
    const double dXSectB         = waferThickness/2.+sectGap/2.;
    const double dYSectB         = dYSect-waferWidth/2.+sectGap/2.;
    const double dZSectB         = sectLength/2.+sectGap/2.;
    const double stRA            = max(abs(stereoInner),abs(stereoOuter));
    const double dXTMou          = 2.*dXSectB+0.1*sep;
    const double dYTMou          = dYSectB*cos(stRA/DWradian)+
                                   dZSectB*sin(stRA/DWradian)+
                                   0.1*sep;
    const double dZTMou          = dYSectB*sin(stRA/DWradian)+
                                   dZSectB*cos(stRA/DWradian)+
                                   0.1*sep;

    vector<double> tckSup(nLay);
    vector<double> tckFl (nLay);
    vector<double> tckClo(nLay);
    for (int i=0; i<nLay; i++) {
	tckSup[i] = (rTckSup[i]/DWfraction)*radSup;
	tckFl [i] = (flTck  [i]/DWfraction)*radSup;
	tckClo[i] = (cOTck  [i]/DWfraction)*radSup;
    }

    // ----------------------------------------------------------
    // --------- Now here actually comes the geometry -----------
    // ----------------------------------------------------------

    // Let's define a section
    agdd::Section * barrelSection=new agdd::Section;;
    // Here are the section's parameters:
    barrelSection->setName("sct_barrel_section");
    barrelSection->setVersion("preliminary version of old geometry");
    barrelSection->setAuthor("Christopher Lester");
    // Better make a "top_volume".
    agdd::Compos * barrelCompos = new agdd::Compos();
    barrelCompos->setName("sct_barrel_barrel");
    barrelSection->setTopVol(barrelCompos);

    theAGDD->addSection(barrelSection); 

    // Here is a Wafer.
    agdd::Box * wafer = new agdd::Box;
    wafer->setName("sct_barrel_wafer");
    wafer->setBox(waferThickness, waferWidth, waferLength);
    wafer->setMaterial("silicon");

    // Here is a DoubleWafer.
    agdd::PosXYZ * topWafer = new agdd::PosXYZ;
    agdd::PosXYZ * botWafer = new agdd::PosXYZ;
    topWafer->setVolume(wafer);
    botWafer->setVolume(wafer);
    topWafer->setXYZ(0., 0., +(waferLength+waferSep)/2.0);
    botWafer->setXYZ(0., 0., -(waferLength+waferSep)/2.0);
    agdd::Compos * doubleWafer = new agdd::Compos;
    doubleWafer->setName("sct_barrel_doubleWafer");
    doubleWafer->addPos(topWafer);
    doubleWafer->addPos(botWafer);

    // Here is an ElectronicsBoard 
    agdd::Box * electronicsBoard = new agdd::Box;
    electronicsBoard->setName("sct_barrel_electronicsBoard");
    electronicsBoard->setBox(tckEle, boardWidth, elecZ);
    electronicsBoard->setMaterial("cheese");

    // Here is a SECT
    agdd::PosXYZ * waferPos = new agdd::PosXYZ;
    agdd::PosXYZ * boardPos = new agdd::PosXYZ;
    waferPos->setVolume(doubleWafer);
    boardPos->setVolume(electronicsBoard);
    boardPos->setXYZ(0., -(waferWidth+boardWidth)/2.0, 0.);
    agdd::Compos * sect = new agdd::Compos;
    sect->setName("sct_barrel_sect");
    sect->addPos(waferPos);
    sect->addPos(boardPos);

    // Here is an SLMP.
    // NOTE:
    // This implementation fixes a bug in the old G3 
    // definition which read something like:
    // df << DeeDoubleUBox(tckLum, wdthLum, hLen*2.);
    // which caused the lump volume to overlap the 
    // insulation volume.
    // I substitute this instead.
    agdd::Box * slmp = new agdd::Box;
    slmp->setName("sct_barrel_slmp");
    slmp->setBox(tckLum, wdthLum, zCylEnd*2.);
    slmp->setMaterial("cheese");

    // Here is a SCAB which is a representation of cabling.
    agdd::Cons * scab = new agdd::Cons;
    scab->setName("sct_barrel_scab");
    scab->setCons(inn1,out1, inn2,out2, zCylEnd);
    scab->setMaterial("cheese");

    // Here is a DoubleSCAB
    agdd::PosXYZ * posPos = new agdd::PosXYZ;
    agdd::PosXYZ * negPos = new agdd::PosXYZ;
    posPos->setVolume(scab);
    negPos->setVolume(scab);
    posPos->setXYZ(0., 0., +zCylEnd/2.);
    negPos->setXYZ(0., 0., -zCylEnd/2.);
    posPos->setRot(0.,   0.*DWdegree, 0.);
    negPos->setRot(0., 180.*DWdegree, 0.);
    agdd::Compos * doubleScab = new agdd::Compos;
    doubleScab->setName("sct_barrel_doubleScab");
    doubleScab->addPos(posPos);
    doubleScab->addPos(negPos);

    // Here are nLay Flanges
    vector<agdd::Tubs *> flange;
    {
	for (int i=0; i<nLay; i++) {
	    flange.push_back(new agdd::Tubs);
	    ostrstream o; o << "sct_barrel_flange" << i+1 << ends;
	    flange[i]->setName(o.str());
	    flange[i]->setTubs(r0Sup[i]-flLen[i], r0Sup[i], tckFl[i]);
	    flange[i]->setMaterial("cheese");
	}
    }

    // Here are nLay CloseOuts
    vector<agdd::Tubs *> closeOut;
    {
	for (int i=0; i<nLay; i++) {
	    closeOut.push_back(new agdd::Tubs);
	    ostrstream o; o << "sct_barrel_closeOut" << i+1 << ends;
	    closeOut[i]->setName(o.str());
	    closeOut[i]->setTubs(r0Sup[i]-tckClo[i], r0Sup[i], cOLen[i]);
	    closeOut[i]->setMaterial("cheese");
	}
    }

    // Here are nLay Supports
    vector<agdd::Tubs *> support;
    {
	for (int i=0; i<nLay; i++) {
	    support.push_back(new agdd::Tubs);
	    ostrstream o; o << "sct_barrel_support" << i+1 << ends;
	    support[i]->setName(o.str());
	    support[i]->setTubs(r0Sup[i]-tckSup[i],
				r0Sup[i],
				2.*(zCylEnd-cOLen[i]-tckFl[i]));
	    support[i]->setMaterial("cheese");
	}
    }

    // Here is the InsulationEnd
    agdd::Tubs * insulationEnd = new agdd::Tubs;
    insulationEnd->setName("sct_barrel_insulationEnd");
    insulationEnd->setTubs(rMin, rMax, tckIso);
    insulationEnd->setMaterial("cheese");

    // Here is the InsulationCyl
    agdd::Tubs * insulationCyl = new agdd::Tubs;
    insulationCyl->setName("sct_barrel_insulationCyl");
    insulationCyl->setTubs(rMax-tckIso, rMax, 2.*zCylEnd);
    insulationCyl->setMaterial("cheese");


    // Here are two BarrelModules, Type1 and Type2 (Positive/negative stereo)
    agdd::PosXYZ * posInnerSectPos  = new agdd::PosXYZ;
    agdd::PosXYZ * negInnerSectPos  = new agdd::PosXYZ;
    agdd::PosXYZ *    OuterSectPos1 = new agdd::PosXYZ;
    agdd::PosXYZ *    OuterSectPos2 = new agdd::PosXYZ;
    posInnerSectPos-> setVolume(sect);
    negInnerSectPos-> setVolume(sect);
    OuterSectPos1->setVolume(sect);
    OuterSectPos2->setVolume(sect);
    posInnerSectPos->setXYZ(-(sectThickness+sectGap)/2., 0., 0.);
    negInnerSectPos->setXYZ(-(sectThickness+sectGap)/2., 0., 0.);
    OuterSectPos1->setXYZ(+(sectThickness+sectGap)/2., 0., 0.);
    OuterSectPos2->setXYZ(+(sectThickness+sectGap)/2., 0., 0.);
    posInnerSectPos->setRot(+stereoInner/DWdegree, 0., 0.); // LESTER Check sign
    negInnerSectPos->setRot(-stereoInner/DWdegree, 0., 0.); // LESTER Check sign

    agdd::Compos * posModule = new agdd::Compos;
    agdd::Compos * negModule = new agdd::Compos;
    posModule->setName("sct_barrel_posBarrelModule");
    negModule->setName("sct_barrel_negBarrelModule");
    posModule->addPos(posInnerSectPos);
    negModule->addPos(negInnerSectPos);
    posModule->addPos(   OuterSectPos1);
    negModule->addPos(   OuterSectPos2);
    // That's the end of the barrel modules.

    // Here are the nLay types of ski.
    vector<agdd::Compos *> ski;
    {
	for (int i=0; i<nLay; i++) {
	    ostrstream o; o << "sct_barrel_skiType" << i+1 << ends;
	    ski.push_back(new agdd::Compos);
	    ski[i]->setName(o.str());

	    int sn; // Will oscillate with the radial staggering.
	    int j;  // Module number
	    for ((sn=-1, j=0); j<nModulesPerSki; (j++, sn=-sn)) {
		agdd::PosXYZ * modulePos = new agdd::PosXYZ;
		modulePos->setVolume( (i==0 || i==2) ? negModule : posModule );
		modulePos->setXYZ(sn*sep/2.0, 0., zModuleOffset[j]);
		ski[i]->addPos(modulePos);
	    }
	}
    }
    // There end all the skis.


    { // Here comes the Barrel. (This is a long one.)
	int i;     // Loop variable to represent the layer.
	int sn;    // Dummy variable indicating skiType.
  
	// First will position entities which
	// there are "one of" per layer.
  
	// Here starts the loop over the layers in the barrel.
	for ((sn=1,i=0); i<nLay; (i++,sn=1-sn)) {
	    // First a couple of temporary variables.      
	    double temp1;
	    double temp2;
    
	    // Now for some layer dependent quantities
	    // which will be useful for the positionings
	    // which follow.
	    const double angle  = 360. * DWdegree / skisPerLayer[i];
	    const double tilt   = tiltOfLayer[i];
	    const double dXSki  = sep/2.+dXTMou;	
	    temp1               = layerRadius[i]+dXSki*cos(tilt/DWradian)-dYTMou*sin(tilt/DWradian);
	    temp2               =                dXSki*sin(tilt/DWradian)+dYTMou*cos(tilt/DWradian);
	    const double rad2   = sqrt(temp1*temp1+temp2*temp2);
	    const double rad1a  = layerRadius[i]-(rad2-layerRadius[i]);
	    temp1               = rad1a;
	    temp2               = wdthLum/2.;
	    temp1               = temp1*temp1-temp2*temp2;
    
	    // Am about to take a square root so ...
	    assert(temp1>=0);
    
	    const double rad1   = sqrt(temp1)-max(tckLum, 2.*out2); 
	    const double rCab   = rad1 + out2;
	    const double yCab   = max(wdthLum/2., rCab*tan((angle/4.)/DWradian));	
	    const double wactrot = sectLength*sin(stRA/DWradian)+waferWidth*cos(stRA/DWradian);
	    const double delX    = layerRadius[i]+(-dXSki*cos(tilt/DWradian)+(wactrot/2.)*sin(tilt/DWradian));
	    const double delY    =                (+dXSki*sin(tilt/DWradian)+(wactrot/2.)*cos(tilt/DWradian));
	    const double alpha   = angle/2.-atan(delY/delX)*DWradian;
	    const double rLump = rad1 + tckLum/2.;
    
	    // OK.  That's all the variables we need to define.
	    // Let's get going with the placements.

	    // First we make a composition for the layer:
	    agdd::Compos * layer = new agdd::Compos;
	    layer->setName("sct_barrel_layer"+makeNumberString(i+1));
	    // and then we place it in the barrel:
	    agdd::PosXYZ * layerPos = new agdd::PosXYZ;
	    layerPos->setVolume(layer);
	    barrelCompos->addPos(layerPos);

	    // Position the Lump and Cables 
	    agdd::MPosPhi * slmpPos = new agdd::MPosPhi;
	    slmpPos->setVolume(slmp);
	    slmpPos->setN(skisPerLayer[i]);
	    slmpPos->setRZ(rLump, 0.);
	    layer->addPos(slmpPos);

	    agdd::MPosPhi * doubleScabPos = new agdd::MPosPhi;
	    doubleScabPos->setVolume(doubleScab);
	    doubleScabPos->setN(skisPerLayer[i]);
	    doubleScabPos->setRZ(sqrt(rCab*rCab+yCab*yCab), 0.);
	    doubleScabPos->setPhi0(atan(yCab/rCab)*DWradian); //LESTER Check sign
	    layer->addPos(doubleScabPos);
	    // That's finished positioning the cables.

	    // Here come the skis.
	    agdd::MPosPhi * skiPos = new agdd::MPosPhi;
	    skiPos->setVolume(ski[i]);
	    skiPos->setN(skisPerLayer[i]);
	    skiPos->setRZ(layerRadius[i], 0.);
	    skiPos->setRot(0., 0., +tiltOfLayer[i]); //LESTER Check sign
	    // The sign, above, was negative when I submitted this geometry
	    // for the first time, and all programs seemed to display
	    // the geometries correctly.
	    // Now (21/2/2000) running persint, I find that the modules
	    // have their electronics boards on the outside.
	    // I have changed the sign to positive to monitor the effect.
	    // Now things seem to look ok, but I don't understand which
	    // was wrong in the first place.  Me or my last viewer?
	    skiPos->setPhi0(alpha+angle);
	    layer->addPos(skiPos);
	    // That completes the skis.

	    // Position the CloseOuts, flanges and Supports
	    agdd::PosXYZ * posCloseOutPos = new agdd::PosXYZ;
	    agdd::PosXYZ * negCloseOutPos = new agdd::PosXYZ;
	    posCloseOutPos->setVolume(closeOut[i]);
	    negCloseOutPos->setVolume(closeOut[i]);
	    posCloseOutPos->setXYZ(0., 0., +(zCylEnd-tckFl[i]-cOLen[i]/2.));
	    negCloseOutPos->setXYZ(0., 0., -(zCylEnd-tckFl[i]-cOLen[i]/2.));
	    layer->addPos(posCloseOutPos);
	    layer->addPos(negCloseOutPos);

	    agdd::PosXYZ * posFlangePos = new agdd::PosXYZ;
	    agdd::PosXYZ * negFlangePos = new agdd::PosXYZ;
	    posFlangePos->setVolume(flange[i]);
	    negFlangePos->setVolume(flange[i]);
	    posFlangePos->setXYZ(0., 0., +(zCylEnd-tckFl[i]/2.));
	    negFlangePos->setXYZ(0., 0., -(zCylEnd-tckFl[i]/2.));
	    layer->addPos(posFlangePos);
	    layer->addPos(negFlangePos);

	    agdd::PosXYZ * supportPos = new agdd::PosXYZ;
	    supportPos->setVolume(support[i]);
	    layer->addPos(supportPos);
	    // That's all the CloseOuts, flanges and Supports done.
    
	}  // End of loop over layers.
  
	// Now come the layer independent placements.
  
	// The insulating stuff
	agdd::PosXYZ * insulationCylPos = new agdd::PosXYZ;
	insulationCylPos->setVolume(insulationCyl);
	barrelCompos->addPos(insulationCylPos);

	agdd::PosXYZ * posInsulationEndPos = new agdd::PosXYZ;
	agdd::PosXYZ * negInsulationEndPos = new agdd::PosXYZ;
	posInsulationEndPos->setVolume(insulationEnd);
	negInsulationEndPos->setVolume(insulationEnd);
	posInsulationEndPos->setXYZ(0., 0., +(hLen-tckIso/2.));
	negInsulationEndPos->setXYZ(0., 0., -(hLen-tckIso/2.));
	barrelCompos->addPos(posInsulationEndPos);
	barrelCompos->addPos(negInsulationEndPos);

    } // That was the end of the Barrel and skis.


    // should leave this scope as it was found
    assert(useMy.secManager().changeDir(string("..")));
  } // END OF BARREL

  { // START OF FORWARD
    assert(useMy.secManager().changeDir(string("Forward")));

    // We're going to read those parameters which interest
    // us from the preXml file, and then produce the 
    // relevant output.
    // The xml file comes to us as the element called doc.

    // Since there are potential name clashes everywhere, let's
    // restrict ourselves to looking in only that part of the preXml file
    // which is relevant to the SCT barrel and this particular geometry.

    // -----------------------------------------------------
    // -------- Some SCT-Barrel specific quantites ---------
    // -------- (close to old AGE structure)       --------- 
    // -----------------------------------------------------

    const ForwardInfo fIS(useMy);
    
    //----------------------------------------------------------
    //---------- Now let's derive some quantities.  ------------
    //---------- We don't HAVE to do this here;     ------------
    //---------- some params are better defined     ------------
    //---------- further on inside loops.           ------------
    //----------------------------------------------------------

    /*
      typedef vector<int> vectorOfIntegers;
      vector<vectorOfIntegers> liRing;
      {
      for (int i=0; i<fIS.nWheels; i++) {
      const string freeForm = fIS.freeFormLiRing[i];
      istrstream is(freeForm.c_str());
      int tmp;
      vector<int> theRingTypes;
      while (is >> tmp) {
      theRingTypes.push_back(tmp-1);  // Note the -1.  I adjust so that 
      // the rings are numbered 0,1,2,... 
      };
      liRing.push_back(theRingTypes);
      };
      }
    */
    
    /*
      const double zCEnd   = 0.5*(fIS.zHiEnd+fIS.zLoEnd);
      const double zHEnd   = 0.5*(fIS.zHiEnd-fIS.zLoEnd);
      const double t2      = 0.5*fIS.thInsu;
      const double zCentre = zCEnd - t2;
      const double zHLen   = zHEnd - t2;
    */
    
    // Some wheel/ring specific things:
    //vector<agdd::Tubs *> ringSpecificZALP;
    //vector<agdd::Compos *> radialCableZALP;
    //vector<double> rIThi; // one element for each ring
    //{ // Define rIThi
    //	for (int ring=0; ring<fIS.nRings; ring++) {
    //	    const double tDet=max(fIS.dZSC,fIS.dZEle);
    //	    rIThi.push_back(2.0*tDet+fIS.dZSub);
    //	};
    //}; 
    // end of wheel/ring specific things.

    // ----------------------------------------------------------
    // ----- Now here actually comes the forward geometry -------
    // ----------------------------------------------------------

    agdd::Section * forwardSection = new agdd::Section;
    forwardSection->setName("sct_fwd_section");
    forwardSection->setVersion("another preliminary version of old geometry");
    forwardSection->setAuthor("Christopher Lester");
    // Better make a "top_volume".
    agdd::Compos * forwardCompos = new agdd::Compos();
    forwardCompos->setName("sct_fwd_forward");
    forwardSection->setTopVol(forwardCompos);

    theAGDD->addSection(forwardSection); 
    
    // We will make both ends, and then rotate them into place.
    // For the moment, here is just one end.
    agdd::Compos * ZEND = new agdd::Compos;
    ZEND->setName("sct_fwd_ZEND");    

     // Here are the ZISOs (insulating surrounding forward sct)
    agdd::Tubs * ZISO2 = new agdd::Tubs;
    ZISO2->setName("sct_fwd_ZISO2InsulationTube");
    ZISO2->setTubs(fIS.rOuEnd-fIS.thInsu, fIS.rOuEnd, 2.0*(fIS.zHLen()-fIS.t2()));
    ZISO2->setMaterial("wool");
    agdd::Tubs * ZISO1 = new agdd::Tubs;
    ZISO1->setName("sct_fwd_ZISO1InsulationEndCap");
    ZISO1->setTubs(fIS.rInEnd, fIS.rOuEnd, 2.0*fIS.t2());
    ZISO1->setMaterial("wool");

   // Now we position the ZISOs in ZEND.
    agdd::PosXYZ * ZISO1Pos = new agdd::PosXYZ;
    ZISO1Pos->setVolume(ZISO1);
    ZISO1Pos->setXYZ(0.0, 0.0, fIS.zHLen());
    ZEND->addPos(ZISO1Pos);
    agdd::PosXYZ * ZISO2Pos = new agdd::PosXYZ;
    ZISO2Pos->setVolume(ZISO2);
    ZEND->addPos(ZISO2Pos);

    // Here follows the "inner part".
    {
	// Make the inner composition
	agdd::Compos * ZINN = new agdd::Compos;
	ZINN->setName("sct_fwd_ZINN");
	for (int wheel=0; wheel<fIS.nWheels; wheel++) {
	    agdd::Volume * ZWEI = newWheel(wheel, 
					   /* liRing, */
					   fIS);
	    { // Position the ZWEI
		agdd::PosXYZ * ZWEIPos = new agdd::PosXYZ;
		ZWEIPos->setXYZ(0.0, 0.0, fIS.zWheel[wheel]-fIS.zCentre());
		ZWEIPos->setRot(0.0, 0.0, fIS.dFWheel[wheel]);
		ZWEIPos->setVolume(ZWEI);
		ZINN->addPos(ZWEIPos);
	    } // Finished positioning the ZWEI
	}
	// now we've made the inner composition, we must position it in ZEND.
	agdd::PosXYZ * ZINNPos = new agdd::PosXYZ;
	ZINNPos->setVolume(ZINN);
	ZINNPos->setXYZ(0.0, 0.0, -fIS.t2());
	ZEND->addPos(ZINNPos);
    }
    // There ended the "inner part".

    { // Here starts the ZOUT - "outer part"
      agdd::PosXYZ * zoutPos = new agdd::PosXYZ;
      zoutPos->setVolume(theZOUT(fIS));
      ZEND->addPos(zoutPos);
    };
    
    // And here is the position:
    {
      agdd::PosXYZ * positiveZForwardPos = new agdd::PosXYZ;
      agdd::PosXYZ * negativeZForwardPos = new agdd::PosXYZ;
      
      positiveZForwardPos->setXYZ(0.0, 0.0, +fIS.zCEnd()); 
      negativeZForwardPos->setXYZ(0.0, 0.0, -fIS.zCEnd());
      
      positiveZForwardPos->setRot(0.0, 0.0, 0.0);
      negativeZForwardPos->setRot(0.0, 180.0*DWdegree, 0.0); // LESTER ?
      
      positiveZForwardPos->setVolume(ZEND);
      negativeZForwardPos->setVolume(ZEND);

      forwardCompos->addPos(positiveZForwardPos);
      forwardCompos->addPos(negativeZForwardPos);
    };
    
    // should leave this scope as it was found
    assert(useMy.secManager().changeDir(string("..")));
  } // END OF FORWARD

  return outStruct;

};