#include "hgeocoils.h"

typedef HVector<float> FVec;

int HGeoCoils::readDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia) {
  // uptil no not defined in database
  // it's not clear how it will look like
  cerr << "Uptil now he coils are not defined in database" << endl;
  return HFAILURE;
}


int HGeoCoils::readDet(HGeoInfo & geoInfo) {
  HGeoParam ele;
  HString volName(' ',120);

// reads Mdc-detector keepinvolume, needed for calculation of coil-keepinvolume
  if (readKeepIn(geoInfo)==0) return HFAILURE;
  keepinName.replace(2,'3');
  if (readKeepIn(geoInfo)==0) return HFAILURE;

  geoInfo.fin.seekg(0, ios::beg);
  geoInfo.fin.clear();
  while(!geoInfo.fin.eof()) {
    geoInfo.fin >> volName;
    if (volName.empty() || volName(0)=='/') geoInfo.fin.getline(volName,120);
    else {
      if (volName==detName) {
        ele.vol=FVec(0.0F,3);
        geoInfo.fin >> ele.vol[0][0] >> ele.vol[0][1];
                 // thickness of upper and lower section
        for(int i=1;i<8;i++) {geoInfo.fin >> ele.vol[i];}
                 // z, rmin, rmax for pgon
        append(detName,ele);
      }
      if (volName==eleName) {
        ele.vol=FVec(0.0F,3);
        int i=0;
        while (i<4) { geoInfo.fin >> ele.vol[i++] >> ele.vol[i++][0]; }
                 // center, radius
        while (i<6) { geoInfo.fin >> ele.vol[i][0] >> ele.vol[i++][1]; }
                 // delta radius inner and outer parts
                 // phi1 and phi2
        geoInfo.fin >> ele.vol[i];
                 // delta-radius keepinvolume and thickness/4 for front and
                 // other parts
        append(eleName,ele);
      }
    }
  }
  if (length()>2) return HSUCCESS;
  cerr << "cannot create Coils" << endl;
  return HFAILURE;
}

int HGeoCoils::createDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia){
// name convention
// 1st character C for coil
// 2nd character F M R for front middle rear section of coil
// 3rd character L R O M I for Left-Keepin Right-Keepin
//               outer middle inner shell of coil
// 4th character C for curved part and B T for bottom top trap
// CKIV keepinvolume for coil
//      with 4th character F R for inner front or rear part

  HGeoParam ele;
  int i,j,l, nrot[2], matno, ivolu;
  char kiName[]="CKIV";
  char tubsName[]="CFLC";
  char trapName[]="CMLB";
  HString kiMaterial="ALUMINIUM$";
  HString kfMaterial="HELIUM$";
  HString lMaterial="HELIUM$";
  HString oMaterial="ALUMINIUM$";
  HString mMaterial="COIL_MIXT$";
  HString iMaterial="COIL_ALU$";
  char vName[5], mother[5];
  float ckpar[19], par[19], tubsparam[5], trapparam[11], pos[2][3];
  float fac=10.F;

  int ns=geoInfo.getNs();
  float asec=360.F / ns;
  float as0=60.F;
  float raddeg=pShapes->shapePar.raddeg;
  float rasec=asec / raddeg;
  float b1=asec-90.F;
  float a2=2.F * asec;
  float b2=a2 - 90.F;
  if (b2<0.0F) b2=360.F+b2;
  float cosb2=cos(b2/raddeg);
  float sinb2=sin(b2/raddeg);

// keepinvolume
  ele=peek(detName);
  float fac1=cosb2 / fac;
  ckpar[0]=as0;
  ckpar[1]=asec;
  ckpar[2]=1;
  ckpar[3]=5;
  ckpar[4]=ele.vol[1][0] / fac;
  ckpar[5]=ele.vol[1][1] * fac1;
  ckpar[6]=ele.vol[1][2] * fac1;
  ckpar[7]=ele.vol[2][0] / fac;
  ckpar[8]=ele.vol[2][1] * fac1;
  ckpar[9]=ele.vol[2][2] * fac1;
  ckpar[10]=ele.vol[4][0] / fac;
  ckpar[11]=ele.vol[4][1] * fac1;
  ckpar[12]=ele.vol[4][2] * fac1;
  ckpar[13]=ele.vol[6][0] / fac;
  ckpar[14]=ele.vol[6][1] * fac1;
  ckpar[15]=ele.vol[6][2] * fac1;
  ckpar[16]=ele.vol[7][0] / fac;
  ckpar[17]=ele.vol[7][1] * fac1;
  ckpar[18]=ele.vol[7][2] * fac1;
  matno=geoMedia.createMedium(kiMaterial,geoInfo);
  if (!matno) {
    cout << "material " << kiMaterial << " for volume " << kiName <<
            " not defined " << endl;
    return HFAILURE;
  }
#ifdef ILSESHOW
  cout << "creating volume  " << kiName << '\n';
#endif
#ifdef WITHGEANT
  GSVOLU(kiName,"PGON",matno,ckpar,19,ivolu);
  HString kiMother="SEC1";
  HString sno;
  i=1;
  for(j=1;j<=ns;j++) {
    if (geoInfo.getnsSet(j)) {
      sno=j;
      kiMother.replace(3,sno);
      char* kiM=kiMother;
      GSPOS(kiName,i,kiM,0.0,0.0,0.0,0,"ONLY");
      i++;
    }
  }
#endif

  strcpy(vName,kiName);
  vName[3]='F';
  pos[0][0]=0.0F;
  pos[0][1]=ele.vol[0][0] /2.F / fac / sinb2;
  pos[0][2]=0.0F;
  par[0]=as0;
  par[1]=asec;
  par[2]=1;
  par[3]=5;
  par[4]=ckpar[4];
  par[5]=ckpar[5] - pos[0][1];
  par[6]=ckpar[6] - pos[0][1];
  par[7]=ckpar[7];
  par[8]=ckpar[8] - pos[0][1];
  par[9]=ckpar[9] - pos[0][1];
  par[10]=ele.vol[3][0] / fac;
  par[11]=ele.vol[3][1] * fac1 - pos[0][1];
  par[12]=ele.vol[3][2] * fac1 - pos[0][1];
  par[13]=ckpar[10];
  par[14]=(ele.vol[4][0] - ele.vol[5][0] +  ele.vol[5][2]) * fac1 - pos[0][1];
  par[15]=ckpar[12] - pos[0][1];
  par[16]=ele.vol[5][0] / fac;
  par[17]=ele.vol[5][2] * fac1 - pos[0][1];
  par[18]=par[17];
  matno=geoMedia.createMedium(kfMaterial,geoInfo);
  if (!matno) {
    cout << "material " << kfMaterial << " for volume " << kiName << 
            " not defined " << endl;
    return HFAILURE;
  }
#ifdef ILSESHOW
  cout << "creating volume  " << vName << '\n';
#endif
#ifdef WITHGEANT
  GSVOLU(vName,"PGON",matno,par,19,ivolu);
  GSPOS(vName,1,kiName,pos[0][0],pos[0][1],pos[0][2],0,"ONLY");
#endif

  vName[3]='R';
  pos[0][1]=ele.vol[0][1] /2.F / fac / sinb2;
  par[3]=4;
  par[4]=par[10];
  par[5]=ele.vol[3][1] * fac1 - pos[0][1];
  par[6]=par[5];
  par[7]=par[16];
  par[8]=ele.vol[5][1] * fac1 - pos[0][1];
  par[9]=ele.vol[5][2] * fac1 - pos[0][1];
  par[10]=ckpar[13];
  par[11]=ckpar[14] - pos[0][1];
  par[12]=ckpar[15] - pos[0][1];
  par[13]=ckpar[16];
  par[14]=ckpar[17] - pos[0][1];
  par[15]=ckpar[18] - pos[0][1];
#ifdef ILSESHOW
  cout << "creating volume  " << vName << '\n';
#endif
#ifdef WITHGEANT
  GSVOLU(vName,"PGON",matno,par,16,ivolu);
  GSPOS(vName,1,kiName,pos[0][0],pos[0][1],pos[0][2],0,"ONLY");
#endif

// curved part of coils
  ele=peek(eleName);
  tubsparam[0]=0.0F;
#ifdef WITHGEANT
  GSROTM(++pShapes->shapePar.geantRotNo,0.,0.,90.,asec,90.,b1);
  GSROTM(++pShapes->shapePar.geantRotNo,0.,0.,90.,a2,90.,b2);
#else
  ++pShapes->shapePar.geantRotNo;
#endif

  float r0;
  for(i=0;i<2;i++) {
    l=2-2*i;
    r0=ele.vol[l+1][0];
    if (i==0) {
      tubsparam[3]=ele.vol[5][0];
      tubsparam[4]=ele.vol[5][1];
    } else {
      tubsName[1]='R';
      tubsparam[3]=ele.vol[5][1] - 360.F;
      tubsparam[4]=ele.vol[5][0];
    }

    for(j=0;j<4;j++) {
      switch(j) {
        case 0:
          strcpy(vName,tubsName);
          strcpy(mother,kiName);
          matno=geoMedia.createMedium(lMaterial,geoInfo);
          if (!matno) {
            cout << "material " << lMaterial << " for volume " << vName <<
                    " not defined " << endl;
            return HFAILURE;
          }
          tubsparam[1]=(r0 + ele.vol[4][0] + ele.vol[4][1] + ele.vol[6][0]) / fac;
          if (i==0) tubsparam[2]=ele.vol[6][1] / fac;
          else tubsparam[2]=ele.vol[6][2] / fac;
          nrot[0]=pShapes->shapePar.geantRotNo - 1;
          nrot[1]=pShapes->shapePar.geantRotNo;
          pos[0][0]=(ele.vol[l][0] * cos(rasec) -
                     tubsparam[2] * fac * sin(rasec)) / fac;
          pos[0][1]=(ele.vol[l][0] * sin(rasec) +
                     tubsparam[2] * fac * cos(rasec)) / fac;
          pos[0][2]=ele.vol[l][2] / fac;
          pos[1][0]=-pos[0][0];
          pos[1][1]=pos[0][1];
          pos[1][2]=pos[0][2];
          break;
        case 1:
          strcpy(vName,tubsName);
          vName[2]='O';
          strcpy(mother,tubsName);
          matno=geoMedia.createMedium(oMaterial,geoInfo);
          if (!matno) {
            cout << "material " << oMaterial << " for volume " << vName <<
                    " not defined " << endl;
            return HFAILURE;
          }
          tubsparam[1]=(r0 + ele.vol[4][0] + ele.vol[4][1]) / fac;
          tubsparam[2]=ele.vol[2][1] / fac;
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=0.0F;
          pos[0][1]=0.0F;
          if (i==0) pos[0][2]=(ele.vol[6][1] - ele.vol[l][1]) / fac;
          else pos[0][2]=(ele.vol[6][2] - ele.vol[l][1]) / fac;
          pos[1][0]=0.0F;
          pos[1][1]=0.0F;
          pos[1][2]=-pos[0][2];
          break;
        case 2:
          strcpy(vName,tubsName);
          vName[2]='M';
          strcpy(mother,tubsName);
          mother[2]='O';
          matno=geoMedia.createMedium(mMaterial,geoInfo);
          if (!matno) {
            cout << "material " << mMaterial << " for volume " << vName <<
                    " not defined " << endl;
            return HFAILURE;
          }
          tubsparam[1]=(r0 + ele.vol[4][0]) / fac;
          tubsparam[2]=ele.vol[2][1] / fac;
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=0.0F;
          pos[0][1]=0.0F;
          pos[0][2]=0.0F;
          break;
        case 3:
          strcpy(vName,tubsName);
          vName[2]='I';
          strcpy(mother,tubsName);
          mother[2]='M';
          matno=geoMedia.createMedium(iMaterial,geoInfo);
          if (!matno) {
            cout << "material " << iMaterial << " for volume " << vName <<
                    " not defined " << endl;
            return HFAILURE;
          }
          tubsparam[1]=(r0 - ele.vol[4][0]) / fac;
          tubsparam[2]=ele.vol[2][1] / fac;
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=0.0F;
          pos[0][1]=0.0F;
          pos[0][2]=0.0F;
          break;
      }
#ifdef ILSESHOW
      cout << "creating volume  " << vName << '\n';
#endif
#ifdef WITHGEANT
      GSVOLU(vName,"TUBS",matno,tubsparam,5,ivolu);
      if (j==0) {
        vName[2]='R';
        GSVOLU(vName,"TUBS",matno,tubsparam,5,ivolu);
        vName[2]='L';
      }
      GSPOS(vName,1,mother,pos[0][0],pos[0][1],pos[0][2],nrot[0],"ONLY");
      if (j==0) {
        vName[2]='R';
        GSPOS(vName,1,mother,pos[1][0],pos[1][1],pos[1][2],nrot[1],"ONLY");
      }
      if (j==1) {
        mother[2]='R';
        GSPOS(vName,2,mother,pos[1][0],pos[1][1],pos[1][2],nrot[0],"ONLY");
      }
#endif
    }
  }


// trap part of coils

  float gamma;
  float alpha=(360.F - ele.vol[5][1] + ele.vol[5][0]) / 2.F;
  float sinalpha=sin(alpha/raddeg);
  float cosalpha=cos(alpha/raddeg);
  float beta=(90.F - alpha) / 2.F;
  float dist=sqrt(pow(ele.vol[2][0]-ele.vol[0][0],2) +
            pow(ele.vol[2][2]-ele.vol[0][2],2));
  trapparam[1]=0.0F;
  trapparam[2]=0.0F;
  trapparam[3]=dist * sinalpha / 2.F / fac;
  trapparam[7]=trapparam[3];
  for(i=0;i<2;i++) {        // i=0 for bottom i=1 for top trap
    if (i==0) {
      trapparam[6]=-beta;
      trapparam[10]=-beta;
      gamma=ele.vol[5][1] - 180.F;
#ifdef WITHGEANT
      GSROTM(++pShapes->shapePar.geantRotNo,135.,240.,45.,240.,90.,150.); // bottom left
      GSROTM(++pShapes->shapePar.geantRotNo,135.,300.,45.,300.,90.,210.); // bottom right
#else
      ++pShapes->shapePar.geantRotNo;
#endif
    } else {
      trapName[3]='T';
      trapparam[6]=beta;
      trapparam[10]=beta;
      gamma=ele.vol[5][0];
#ifdef WITHGEANT
      GSROTM(++pShapes->shapePar.geantRotNo,130.,240.,40.,240.,90.,150.); // top left
      GSROTM(++pShapes->shapePar.geantRotNo,130.,300.,40.,300.,90.,210.); // top right
#else
      ++pShapes->shapePar.geantRotNo;
#endif
    }

    for(j=0;j<4;j++) {
      switch(j) {
        case 0: {
          strcpy(vName,trapName);
          strcpy(mother,kiName);
          matno=geoMedia.createMedium(lMaterial,geoInfo);
          trapparam[0]=ele.vol[6][2] / fac;
          trapparam[4]=(ele.vol[3][0] + ele.vol[4][0] + ele.vol[4][1]
                        + ele.vol[6][0]) /2.F /fac;
          trapparam[5]=(ele.vol[1][0] + ele.vol[4][0] + ele.vol[4][1]
                        + ele.vol[6][0]) /2.F /fac;
          trapparam[8]=trapparam[4];
          trapparam[9]=trapparam[5];
          nrot[0]=pShapes->shapePar.geantRotNo - 1;
          nrot[1]=pShapes->shapePar.geantRotNo;
          float xi=(dist -
                  (trapparam[4] + trapparam[5]) * fac * cosalpha) / 2.F;
          float zi=(trapparam[4] + trapparam[5]) * fac * sinalpha / 2.F;
          if (i==1) zi=-zi;
          gamma=ele.vol[5][0] + 2.F * beta;
          float singamma=sin(gamma/raddeg);
          float cosgamma=cos(gamma/raddeg);
          float x=ele.vol[0][0] + zi * singamma + xi * cosgamma;
          float z=ele.vol[0][2] + zi * cosgamma - xi * singamma;
          pos[0][0]=(x * cos(rasec) - ele.vol[6][2] * sin(rasec)) / fac;
          pos[0][1]=(x * sin(rasec) + ele.vol[6][2] * cos(rasec)) / fac;
          pos[0][2]=z / fac;
          pos[1][0]=-pos[0][0];
          pos[1][1]=pos[0][1];
          pos[1][2]=pos[0][2];
          break;
        }
        case 1:
          strcpy(vName,trapName);
          vName[2]='O';
          strcpy(mother,trapName);
          matno=geoMedia.createMedium(oMaterial,geoInfo);
          trapparam[0]=ele.vol[2][1] / fac;
          trapparam[4]=(ele.vol[3][0] + ele.vol[4][0] +
                        ele.vol[4][1]) / 2.F / fac;
          trapparam[5]=(ele.vol[1][0] + ele.vol[4][0]
                        + ele.vol[4][1]) / 2.F / fac;
          trapparam[8]=trapparam[4];
          trapparam[9]=trapparam[5];
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=ele.vol[6][0] /2.F / fac;
          pos[0][1]=0.0F;
          pos[0][2]=(ele.vol[2][1] - ele.vol[6][2]) / fac;
          if (i==1) pos[0][0]=-ele.vol[6][0] / 2.F / fac;
          pos[1][0]=pos[0][0];
          pos[1][1]=pos[0][1];
          pos[1][2]=-pos[0][2];
          break;
        case 2:
          strcpy(vName,trapName);
          vName[2]='M';
          strcpy(mother,trapName);
          mother[2]='O';
          matno=geoMedia.createMedium(mMaterial,geoInfo);
          trapparam[0]=ele.vol[2][1] / fac;
          trapparam[4]=(ele.vol[3][0] + ele.vol[4][0]) / 2.F / fac;
          trapparam[5]=(ele.vol[1][0] + ele.vol[4][0]) / 2.F / fac;
          trapparam[8]=trapparam[4];
          trapparam[9]=trapparam[5];
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=ele.vol[4][1] / 2.F / fac;
          pos[0][1]=0.0F;
          pos[0][2]=0.0F;
          if (i==1) pos[0][0]=-ele.vol[4][1] / 2.F / fac;
          break;
        case 3:
          strcpy(vName,trapName);
          vName[2]='I';
          strcpy(mother,trapName);
          mother[2]='M';
          matno=geoMedia.createMedium(iMaterial,geoInfo);
          trapparam[0]=ele.vol[2][1] / fac;
          trapparam[4]=(ele.vol[3][0] - ele.vol[4][0]) /2.F / fac;
          trapparam[5]=(ele.vol[1][0] - ele.vol[4][0]) /2.F / fac;
          trapparam[8]=trapparam[4];
          trapparam[9]=trapparam[5];
          nrot[0]=0;
          nrot[1]=0;
          pos[0][0]=ele.vol[4][0] / fac;
          pos[0][1]=0.0F;
          pos[0][2]=0.0F;
          if (i==1) pos[0][0]=-ele.vol[4][0] / fac;
          break;
      }
#ifdef ILSESHOW
      cout << "creating volume  " << vName << '\n';
#endif
#ifdef WITHGEANT
      GSVOLU(vName,"TRAP",matno,trapparam,11,ivolu);
      if (j==0) {
        vName[2]='R';
        GSVOLU(vName,"TRAP",matno,trapparam,11,ivolu);
        vName[2]='L';
      }
      GSPOS(vName,1,mother,pos[0][0],pos[0][1],pos[0][2],nrot[0],"ONLY");
      if (j==0) {
        vName[2]='R';
        GSPOS(vName,1,mother,pos[1][0],pos[1][1],pos[1][2],nrot[1],"ONLY");
      }
      if (j==1) {
        mother[2]='R';
        GSPOS(vName,2,mother,pos[1][0],pos[1][1],pos[1][2],nrot[0],"ONLY");
      }
#endif
    }
  }
  return HSUCCESS;
}