#include "geantdef.h"
#include "hgeantinput.h"
#include "hgeorich.h"
#include <ctype.h>
#include "hnamedlistfc.h"
#include "utility.h"

#ifdef WITHGEANT
#include "cfortran.h"
#include "geant321.h"
#include "geafuncdec.h"
#endif

typedef HFloatVec FVec;
typedef HVector<double> DVec;

// Common Blocks for target informations to be availible
// also in the GEANT Fortran routines
#ifdef WITHGEANT
#define SEGMAXNO 20
struct GEANTTARGET_DEF {
  int nSegments;
  int targetMednum;
  float size[SEGMAXNO][3];
  float position[SEGMAXNO][3];
  char shape[5];
};
#define GEATARGET COMMON_BLOCK(GEATARGET,geatarget)
COMMON_BLOCK_DEF(GEANTTARGET_DEF,GEATARGET);
#endif

// ************************************************************************
// ----------------------  readDet  ---------------------------------------
// ************************************************************************

int HGeoRich::readDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia) {
#ifdef WITHORACLE
  // get object index of target mothervolume
  // no daughters of this volume should be read when reading the Rich
  HString mName("RTAM");
  int breakLoop=geoInfo.geoDb->getObjIndex(mName,1);
  noOfKeepIn=readSingleModule(geoInfo,geoMedia,breakLoop);
  if (noOfKeepIn<=0) return HFAILURE;

// read the target
  HString s("target");
  if (geoInfo.connectInput(s)==HSUCCESS) {
    int oType=getTargetVersion(geoInfo);  //target version
    if (oType<=1) return HFAILURE;
    // get object index of target mothervolume with desired target version
    // and read all daughters
    int iMo=geoInfo.geoDb->getObjIndex(mName,oType);
    int nD=0;
    if (iMo>0) nD=readDaughters(geoInfo,iMo,geoMedia,0);
    if (nD==0) {
      cout<<"*************************************************"<<endl;
      cout<<"******    ATTENTION: No target defined!    ******"<<endl;
      cout<<"******  Cannot read target from database!  ******"<<endl;
      cout<<"*************************************************"<<endl;
      return HFAILURE;
    }    
  } else {
    cout<<"*************************************************"<<endl;
    cout<<"******    ATTENTION: No target defined!    ******"<<endl;
    cout<<"*************************************************"<<endl;
  }
  if (length()>0) return HSUCCESS;
#endif
  return HFAILURE;
}

// ------------------------------------------------------------------------

int HGeoRich::readDet(HGeoInfo & geoInfo) {
  noOfKeepIn=readModule(geoInfo);
  if (noOfKeepIn<=0) return HFAILURE;
  geoInfo.closeFile();

// read the target
  HString s("target");
  if (geoInfo.connectInput(s)==HSUCCESS) {
    if (readTarget(geoInfo)==HFAILURE) return HFAILURE;
  } else {
    cout<<"*************************************************"<<endl;
    cout<<"******    ATTENTION: No target defined!    ******"<<endl;
    cout<<"*************************************************"<<endl;
  }
  return HSUCCESS;
}


// ************************************************************************
// ----------------------  readTarget  ------------------------------------
// ************************************************************************
int HGeoRich::readTarget(HGeoInfo & geoInfo){
  HGeoParam ele;
  HString volName(' ',120);
  int lName, rc;
  ele.inout=1;
  cout << "reading target" << endl;
  while(!geoInfo.fin.eof()) {
    geoInfo.fin >> volName;
    if (volName.empty() || volName(0)=='/') geoInfo.fin.getline(volName,120);
    else {
      if (volName[0]=='T') {
        lName=volName.length();
        geoInfo.fin >> ele.mother.assign(10);
        if (lName<=4) rc=readParam(geoInfo,ele);
        else rc=readCopy(geoInfo,ele);
        if (rc==HSUCCESS) append(volName,ele);
        else return HFAILURE;
      }
    }
  }
  return HSUCCESS;
}


// ************************************************************************
// ----------------------  createDet  -------------------------------------
// ************************************************************************

int HGeoRich::createDet(HGeoInfo & geoInfo, HGeoMedia & geoMedia){
  int matNo, copyNo, copyId, ivolu, lName, nSeg=0;
  HGeoParam ele;
  HGeoParam mo;
  HNamedList<int> copyList;
  HString volName, copyName;
  int lDetName=detName.length();
  int lEleName=eleName.length();
  reset();
  do {
    volName=name();
    if (volName(0,lDetName)==detName || volName(0,lEleName)==eleName 
          || volName[0]=='T') {
      lName=volName.length();
      ele=element();
      if (ele.inout==1) {
        if (volName(0,lDetName)==detName) {
          mo.shape="PGON";
          mo.geaPos=0.0F;
          mo.geaRot=0.0F;
          mo.geaRot[0]=1.0F;
          mo.geaRot[4]=1.0F;
          mo.geaRot[8]=1.0F;
        }
        else mo=peek(ele.mother);
        HGeoShapePar & shapePar=pShapes->shapePar;
        copyName=volName(0,4);
        copyNo=1;
        copyId=1;
        if (lName>4) {
          if (copyList.find(copyName)) copyNo=++copyList.element(0);
          else copyList.append(copyName,copyNo);
          HString s=volName(4,lName-4);
          copyId=int(s);
        }
        char* cN=copyName;
        char* eS=ele.shape;
        if (copyNo==1) {
          matNo=geoMedia.createMedium(ele.material,geoInfo);
          if (!matNo) {
            cerr << "material for " << volName << " not defined" << endl;
            return HFAILURE;
          }
          if (pShapes->calcShape(ele)==HSUCCESS) {
#ifdef WITHGEANT
            GSVOLU(cN,eS,matNo,&(shapePar.volPar[0]),shapePar.noPar,ivolu);
#endif
          }
          else return HFAILURE;
        }
        if (pShapes->posShape(ele,mo,shapePar)==HSUCCESS) {
          HString tmp=ele.mother(0,4);
          char* moN=tmp;
#ifdef WITHGEANT
          char cTmp[]="ONLY"; 
          GSPOS(cN,copyId,moN,shapePar.pos[0],shapePar.pos[1],
                shapePar.pos[2],shapePar.nRot,cTmp);
#endif
          replace(volName,ele);
        }
        else return HFAILURE;
#ifdef ILSESHOW
        cout << volName << "  matNo: " << matNo << '\n';
        pShapes->showShapePar();
#endif
        if (copyNo==1 && geoMedia.getCurrentSensFlag()==1) {
          if (createHits(volName)==HFAILURE) {
            cerr << "cannot create hits for " << volName << endl;
            return HFAILURE;
          }
        }

#ifdef WITHGEANT
        // fill Common block for GEANT Fortran routines
        if (copyName=="TARG") {
          if (copyNo==1) {
            GEATARGET.targetMednum=matNo;
            strcpy(GEATARGET.shape,ele.shape);
          }
          if (ele.shape=="BOX " || ele.shape=="TUBE"  || ele.shape=="ELTU") {
            if (ele.shape=="TUBE") GEATARGET.size[nSeg][0]=shapePar.volPar[1]*10.F;
            else GEATARGET.size[nSeg][0]=shapePar.volPar[0]*10.F;
            GEATARGET.size[nSeg][1]=shapePar.volPar[1]*10.F;   // in mm
            GEATARGET.size[nSeg][2]=shapePar.volPar[2]*10.F;
          } else {
            cerr << "wrong shape for target" << endl;
            return HFAILURE;
          }
          FVec tPos(0.F,3);
          calcTargetLabPos(ele,tPos);
          GEATARGET.position[nSeg][0]=tPos[0];  // in mm
          GEATARGET.position[nSeg][1]=tPos[1];
          GEATARGET.position[nSeg][2]=tPos[2];
          nSeg++;
        }
#endif
      }
    }
  } while(next(1));
#ifdef WITHGEANT
  GEATARGET.nSegments=nSeg;

#ifdef ILSESHOW
  cout<<"*********** target **************"<<endl;
  cout<<GEATARGET.nSegments<<"  "<<GEATARGET.targetMednum<<"  "
      <<GEATARGET.shape<<endl;  
  for(int m=0;m<nSeg;m++) {
    cout<<GEATARGET.size[m][0]<<"  "<<GEATARGET.size[m][1]<<"  "
        <<GEATARGET.size[m][2]<<"  "
        <<GEATARGET.position[m][0]<<"  "<<GEATARGET.position[m][1]<<"  "
        <<GEATARGET.position[m][2]<<endl;
  }
  cout<<"*********************************"<<endl;
#endif

#endif
  return HSUCCESS;
}


// ************************************************************************
// ----------------------  createHits  ------------------------------------
// ************************************************************************

int HGeoRich::createHits(HString & volName) {
  extern HGeantInput geantInput;
  int idtype, iset, idet;
  HString copyName=volName(0,4);
  idtype=1;
  if (hitsRead<1) {
    int i=0;
    int n=geantInput.getHitNoFiles();
    while (i<=n) {
      hitsFileName=geantInput.getHitFileName(i);
      if (hitsFileName.find("rich")>=0) break;
      i++;
    }
    if (i>n) {
      cerr << "No file for hit definition of Rich found" << endl;
      return HFAILURE;
    }
    if (readHits(hitsFileName)==HFAILURE) return HFAILURE;
  }
  if (restoreTree(volName)==HFAILURE) return HFAILURE;
#ifdef WITHGEANT
  char* cS=chset;
  char* cN=copyName;
  GSDET(cS,cN,nv,chnmsv,nbitsv,idtype,1000,1000,iset,idet);
  GSDETH(cS,cN,nh,chnamh,nbitsh,orig,fact);
#endif
#ifdef ILSESHOW
  showHits(volName,idtype);
#endif
  return HSUCCESS;
}


// ************************************************************************
// ----------------------  writeToFile  -----------------------------------
// ************************************************************************

int HGeoRich::writeToFile(HString & modName, HGeoInfo & geoInfo) {
  ofstream fout;
  HString volName, copyName;
  HGeoParam ele;

  HString mName[2]={modName,"target"};
  char eName[2]={eleName[0],'T'};
  for (int k=0;k<2;k++) {
    connectOutputFile(mName[k],geoInfo,fout);
    reset();
    do {
      volName=name();
      if (volName[0]==eName[k]) {
        ele=element();
        fout << volName << '\n';
        if (ele.inout==0) fout << ele.inout << '\n';
        fout << ele.mother << '\n';
        fout << ele.shape << '\n';
        fout << ele.material << '\n';
        fout.precision(3);
        if (pShapes->writePoints(fout,ele)==HFAILURE) {
          fout.close();
          return HFAILURE;
        }
        fout<<ele.trans[0]<<"  "<<ele.trans[1]<<"  "<<ele.trans[2]<<'\n';
        fout.precision(6);
        for (int i=0;i<9;i++) {fout<<"  "<<ele.rot[i];}
        fout <<'\n'<<
             "// ***********************************************************"<<'\n';
      }
    } while(next(1));
    fout.close();
  }
  return HSUCCESS;
}


// ************************************************************************
// ----------------------  getTargetVersion  ------------------------------
// ************************************************************************

int HGeoRich::getTargetVersion(HGeoInfo & geoInfo) {
  extern HGeantInput geantInput;
  int vi=0;
#ifdef WITHORACLE
  HString version=geantInput.getTargetVersion();
  if (version.length()==0) vi=geoInfo.geoDb->showTargetVersions();
  else {
    int i=int(version[0]);
    if (isdigit(i)) vi=int(version);
    else vi=geoInfo.geoDb->getTargetVersion(version);
  }
#endif
  if (vi<=1) {
    cerr << "cannot read target version " << vi << " from database" << endl;
  }
  return vi;
}


// ************************************************************************
// ----------------------  calcTargetLabPos  ------------------------------
// ************************************************************************

void HGeoRich::calcTargetLabPos(HGeoParam & ele, FVec & targPos) {
  DVec t(0.,3), tg(0.,3);
  for(int i=3;i--;tg[i]=(double)ele.geaPos[i]);
  if (ele.labRot.norm()<.1) {
    // when reading from file lab transformation is not defined
    HString mName(ele.mother);
    HGeoParam mother;
    mother=peek(mName);
    DVec r(0.0,9);

    // lab transformation of mother coordinate system
    // only for first target if all have the same mother
    if (mother.labRot.norm()<.1) {
      HGeoParam mo;
      mo=mother;
      HGeoParam gm;
      DVec gt(0.,3), gr(0.0,9);
      for(int i=3;i--;mother.labTrans[i]=(double)mo.trans[i]);
      for(int i=9;i--;mother.labRot[i]=(double)mo.rot[i]);
      while (mo.mother.find("CAVE")<0) {
        t=mother.labTrans;
        r=mother.labRot;
        gm=peek(mo.mother);
        for(int i=3;i--;gt[i]=(double)gm.trans[i]);
        for(int i=9;i--;gr[i]=(double)gm.rot[i]);
        multiplyMatrices(gr,r,mother.labRot);
        rotateVector(gr,t,mother.labTrans);
        mother.labTrans+=gt;
        mo=gm;
      }
      replace(mName,mother);
    }

    // lab transformation of target coordinate system
    for(int i=3;i--;t[i]=(double)ele.trans[i]);
    for(int i=9;i--;r[i]=(double)ele.rot[i]);
    multiplyMatrices(mother.labRot,r,ele.labRot);
    rotateVector(mother.labRot,t,ele.labTrans);
    ele.labTrans+=mother.labTrans;
  }
  rotateVector(ele.labRot,tg,t);
  t+=ele.labTrans;
  for(int i=3;i--;targPos[i]=(float)Round(t[i],3));
  return;
}