#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 "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"

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")));
  // I asserted these just to check that I succeeded.
  
  // Now we start the geometry proper.
  
  // 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");
  barrelSection->setTopVol(barrelCompos);

  agdd::AGDD * theAGDD=new agdd::AGDD;
  theAGDD->addSection(barrelSection); 

  outStruct.theBuildable=theAGDD;

{{ // This is where the old SCT.geo.hh started

  // 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 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; \
    }                                                         \
  };


assert(useMy.secManager().changeDir(string("SCT")));
assert(useMy.secManager().changeDir(string("Barrel")));

  // -----------------------------------------------------
  // -------- 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 -----------
  // ----------------------------------------------------------

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

// 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_doubleWafer");
doubleWafer->addPos(topWafer);
doubleWafer->addPos(botWafer);

// Here is an ElectronicsBoard 
agdd::Box * electronicsBoard = new agdd::Box;
electronicsBoard->setName("sct_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_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_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_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_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_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_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_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_insulationEnd");
insulationEnd->setTubs(rMin, rMax, tckIso);
insulationEnd->setMaterial("cheese");

// Here is the InsulationCyl
agdd::Tubs * insulationCyl = new agdd::Tubs;
insulationCyl->setName("sct_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_posBarrelModule");
negModule->setName("sct_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_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_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.

}} // This is where the old SCT.geo.hh ended

return outStruct;

};