#include "hgeotof.h"

typedef HVector<float> FVec;

int HGeoTof::readDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia) {
  // reads shower-detector keepinvolumeneeded for calculation of tof-keepinvolume
  int n=readKeepIn(geoInfo,geoMedia);
  if (n==0) return HFAILURE;
  for(int i=1;i<=numModules;i++) {
    changeName(detName,1,i);
    changeName(eleName,1,i);
    noOfKeepIn[i]=readModule(geoInfo,geoMedia);
  }
  if (length()>0) return HSUCCESS;
  else return HFAILURE;
}


int HGeoTof::readDet(HGeoInfo & geoInfo) {
  int n=readKeepIn(geoInfo);    // reads shower-detector keepinvolume
                                // needed for calculation of tof-keepinvolume
  if (n==0) return HFAILURE;
  for(int i=1;i<=numModules;i++) {
    changeName(detName,1,i);
    changeName(eleName,1,i);
    noOfKeepIn[i]=readModule(geoInfo);
    geoInfo.fin.seekg(0, ios::beg);
    geoInfo.fin.clear();
  }
  if (length()>0) return HSUCCESS;
  return HFAILURE;
}


int HGeoTof::createDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia){
  reset();
  if (createKeepIn(geoInfo,geoMedia)==HFAILURE) return HFAILURE;
  keepinName="";
  for(int i=1;i<=numModules;i++) {
    reset();
    changeName(detName,1,i);
    changeName(eleName,1,i);
    if(noOfKeepIn[i]>0) createModule(geoInfo,geoMedia);
  }
  return HSUCCESS;
}


int HGeoTof::createKeepIn(HGeoInfo & geoInfo, HGeoMedia & geoMedia){
  HGeoParam ele;
  int matNo, npar, ivolu, pointno, no=0;
  float tkpar[22], ms, as, mt, at;
  float delta=1.F;
  float deltay1=300.F;
  float deltay2=350.F;
  float fac=10.0F;
  FVec labpoint(0.0F,3), s1(0.0F,3), s2(0.0F,3);
  FVec t1(0.0F,3), t2(0.0F,3), t3(0.0F,3), t4(0.0F,3);
  HString sno;
  HString tk1Name="TCM1";
  HString tk2Name="TMM1";
  HString tk3Name="TOM1";
  HString moName="SEC1";
  HString tofName=detName + "1";
  HString showerName=keepinName + "1";
  char* tN;

  HString material="VAKUUM$";
  matNo=geoMedia.createMedium(material,geoInfo);
  if (matNo==0) {
    cerr << "material for tof-keepinvolume not defined" << endl;
    return HFAILURE;
  }

  int ns=geoInfo.getNs();
  tkpar[0]=60.0F;
  tkpar[1]=360.0F / ns;
  tkpar[2]=1;

  for(int i=1;i<=ns;i++) {
    if(geoInfo.getnsSet(i)) {
      sno=i;
      tk1Name.replace(3,sno);
      tk2Name.replace(3,sno);
      tk3Name.replace(3,sno);
      moName.replace(3,sno);
      tofName.replace(4,sno);
      showerName.replace(3,sno);

      int in[3]={0,0,0};
      int n=1;
      do {
        changeName(tofName,1,n);
        ele=element(tofName);
        if (ele.inout) in[0]=1;
        n++;
      } while(in[0]==0 && n<=14);
      n=15;
      do {
        changeName(tofName,1,n);
        ele=element(tofName);
        if (ele.inout) in[1]=1;
        n++;
      } while(in[1]==0 && n<=18);
      n=19;
      do {
        changeName(tofName,1,n);
        ele=element(tofName);
        if (ele.inout) in[2]=1;
        n++;
      } while(in[2]==0 && n<=22);

// keepinvolume for tof 1 - 14
      for(int j=0;j<3;j++) {
        switch (j) {
          case 0:                        // keepinvolume for tof 1 - 14
            tN=tk1Name;
            npar=16;
            changeName(tofName,1,1);
            ele=element(tofName);
            transCatiaToLab(ele,0,t1);
            changeName(tofName,1,14);
            ele=element(tofName);
            transCatiaToLab(ele,1,t2);
            ele=element(showerName);
            transCatiaToLab(ele,0,s1);
            transCatiaToLab(ele,1,s2);
            ms=(s2[1] - s1[1]) / (s2[2] - s1[2]);
            as=s1[1] - ms * s1[2];
            mt=(t2[1] - t1[1]) / (t2[2] - t1[2]);
            at=t1[1] - mt * t1[2];
            tkpar[3]=4;
            tkpar[4]=t2[2] / fac;
            tkpar[5]=t2[1] / fac;
            tkpar[6]=tkpar[5];
            tkpar[7]=as / (t2[1] / t2[2] - ms) / fac;
            tkpar[8]=(mt * tkpar[7] * fac + at) / fac;
            tkpar[9]=(ms * tkpar[7] * fac + as) / fac;
            tkpar[10]=t1[2] / fac;
            tkpar[11]=t1[1] / fac;
            tkpar[12]=(ms * tkpar[10] * fac + as) / fac;
            tkpar[13]=as / (t1[1] / t1[2] - ms) / fac;
            tkpar[14]=(ms * tkpar[13] * fac + as) / fac;
            tkpar[15]=tkpar[14];
            break;

          case 1:                        // keepinvolume for tof 15 - 18
            tN=tk2Name;
            npar=22;
            tkpar[3]=6;
            tkpar[16]=tkpar[4];
            tkpar[17]=tkpar[5];
            tkpar[19]=tkpar[7];
            tkpar[20]=tkpar[9];
            tkpar[21]=tkpar[9];
            changeName(tofName,1,15);
            ele=element(tofName);
            transCatiaToLab(ele,0,t3);
            transCatiaToLab(ele,4,t4);
            t1=t3 - (t4 - t3) * delta;
            changeName(tofName,1,18);
            ele=element(tofName);
            transCatiaToLab(ele,1,t2);
            mt=(t2[1] - t1[1]) / (t2[2] - t1[2]);
            at=t1[1] - mt * t1[2];
            ms=(t2[1] + deltay1 - tkpar[20] * fac) / (t2[2] - tkpar[19] * fac);
            as=tkpar[20] * fac - ms * tkpar[19] * fac;
            tkpar[4]=t2[2] / fac;
            tkpar[5]=t2[1] / fac;
            tkpar[6]=tkpar[5];
            tkpar[7]=as / (t2[1] / t2[2] - ms) / fac;
            tkpar[8]=(mt * tkpar[7] * fac + at) / fac;
            tkpar[9]=(ms * tkpar[7] * fac + as) / fac;
            tkpar[10]=t1[2] / fac;
            tkpar[11]=t1[1] / fac;
            tkpar[12]=(ms * tkpar[10] * fac + as) / fac;
            tkpar[13]=t4[2] / fac;
            tkpar[14]=t4[1] / fac;
            tkpar[15]=(ms * tkpar[13] * fac + as) / fac;
            tkpar[18]=(ms * tkpar[16] * fac + as) / fac;
            break;

          case 2:                        // keepinvolume for tof 19 - 22
            tN=tk3Name;
            npar=22;
            tkpar[3]=6;
            tkpar[16]=tkpar[4];
            tkpar[17]=tkpar[5];
            tkpar[19]=tkpar[7];
            tkpar[20]=tkpar[9];
            tkpar[21]=tkpar[9];
            changeName(tofName,1,19);
            ele=element(tofName);
            transCatiaToLab(ele,0,t1);
            transCatiaToLab(ele,4,t4);
            changeName(tofName,1,22);
            ele=element(tofName);
            transCatiaToLab(ele,1,t2);
            mt=(t2[1] - t1[1]) / (t2[2] - t1[2]);
            at=t1[1] - mt * t1[2];
            ms=(t2[1] + deltay2 - tkpar[20] * fac) / (t2[2] - tkpar[19] * fac);
            as=tkpar[20] * fac - ms * tkpar[19] * fac;
            tkpar[4]=t2[2] / fac;
            tkpar[5]=t2[1] / fac;
            tkpar[6]=tkpar[5];
            tkpar[7]=as / (t2[1] / t2[2] - ms) / fac;
            tkpar[8]=(mt * tkpar[7] * fac + at) / fac;
            tkpar[9]=(ms * tkpar[7] * fac + as) / fac;
            tkpar[10]=t1[2] / fac;
            tkpar[11]=t1[1] / fac;
            tkpar[12]=(ms * tkpar[10] * fac + as) / fac;
            tkpar[13]=t4[2] / fac;
            tkpar[14]=t4[1] / fac;
            tkpar[15]=(ms * tkpar[13] * fac + as) / fac;
            tkpar[18]=(ms * tkpar[16] * fac + as) / fac;
            break;
        } 
        if (in[j]) {
#ifdef WITHGEANT
          char* mN=moName;
          GSVOLU(tN,"PGON",matNo,tkpar,npar,ivolu);
          GSPOS(tN,1,mN,0.0,0.0,0.0,0,"ONLY");
#endif
#ifdef ILSESHOW
          cout << tN << "  matNo: " << matNo << '\n';
          for(int l=0;l<npar;l++) {cout << tkpar[l] << '\n';}
#endif
        }
      }
    }
  }
  return HSUCCESS;
}


void HGeoTof::changeName(HString & name, int firstPos, int no) {
  int j=(int)no/10;
  HString s1no=j;
  HString s2no=no-j*10;
  name.replace(firstPos,s1no);
  name.replace(++firstPos,s2no);
  return;
}


void HGeoTof::transCatiaToLab(HGeoParam & ele, int pno, FVec & p) {
  float ctheta=ele.rot[8];
  float stheta=sqrt(1.F-ctheta*ctheta);
  float dz=(float)ele.trans.norm() + ele.vol[pno][2];
  p[1]=dz * stheta + ele.vol[pno][1] * ctheta;
  p[2]=dz * ctheta - ele.vol[pno][1] * stheta;
  return;
}


int HGeoTof::createHits(HString & copyName) {
  int l1, l2, idtype, iset, idet;
  char setName[]="META";
  int hitNh=20;
  char hitNamesh[20][5]={"X   ","Y   ","Z   ","E1  ","TOF ",
                         "E2  ","PART","SECT","ITRA","ETRA",
                         "IDET","E3  ","E4  ","E5  ","E6  ",
                         "INWV","Px  ","Py  ","Pz  ","E   "};
  int hitNbitsh[20]={30,30,30,30,30,30,8,8,12,30,
                     12,30,30,30,30,8,30,30,30,30};
  float hitOrig[20]={340.F,340.F,340.F,0.F,0.F,0.F,0.F,0.F,0.F,0.F,
                     0.F,0.F,20.F,0.F,0.F,0.F,5.F,5.F,5.F,5.F};
  float hitFact[20]={1.e4F,1.e4F,1.e4F,1.e8F,1.e14F,1.e5F,1.F,1.F,1.F,1.e8F,
                     1.F,1.e6F,1.e6F,1.e6F,1.e6F,1.F,1.e6F,1.e6F,1.e6F,1.e6F};
  l1=(int)copyName(1,1);
  l2=(int)copyName(2,1);
  idtype=100+10*l1+l2;
#ifdef WITHGEANT
  char* cN=copyName;
  GSDETV(setName,cN,idtype,1000,1000,iset,idet);
  GSDETH(setName,cN,hitNh,hitNamesh,hitNbitsh,hitOrig,hitFact);
#endif
#ifdef ILSESHOW
  cout << "detector: " << copyName << "   sensitive volume: " << setName <<
          "   idtype: " << idtype << '\n';
#endif
  return HSUCCESS;
}