#define ACCL1
#define ACCL2
#define ACCL3
#define ACCL4
#define ACCL5

#ifdef ACCL1
#define ACCL0
#endif
#ifdef ACCL2
#define ACCL0
#define xOFLA        // omit the flasher DOM
#define xROMB        // use rhomb cells aligned with the array
#endif

#define ASENS        // enable hole ice sensitivity
#define xTILT        // enable tilted ice layers

#define MLTP 30      // number of multiprocessors
#ifdef ACCL5
#define WNUM 34      // optimized for 30 multiprocessors
#else
#define WNUM MLTP    // cannot be divisible by 13
#endif
#define MAXLYS 180   // maximum number of ice layers
#define MAXGEO 5200  // maximum number of OMs

#define OVER 10      // size of photon bunches along the muon track
#define NPHO 1042    // maximum number of photons propagated by one thread
#define NTHR 320     // NTHR*NBLK should not exceed the count of different random number multipliers
#define HNUM 1048576 // size of the output hit buffer, must hold hits from up to NPHO*NTHR*NBLK photons

static const float xx=1.e-5;
static const float OVR=5;
static const float R=OVR*0.178;
static const float fpi=3.141592653589793f;

union name{
  unsigned short k;
  struct{
    unsigned char dom;
    unsigned char str;
    unsigned short one;
  };
  float f;

  name(){
    f=1.0;
  }
};

struct DOM{
  float r[3];
};

struct hit{
  unsigned int i;
  float t;
  unsigned int n;
  float z;
};

struct photon{
  int q;
  float n[3];
  float t, r[3];
};

#ifdef ACCL2
static const int cx=10, cy=10, cz=20;
#ifdef ROMB
static const float dir1=9.3, dir2=129.3;
#endif
#endif

struct ices{
  float wvl;                  // wavelength of this block
  float ocm;                  // 1 / speed of light in medium
  float coschr, sinchr;       // cos and sin of the cherenkov angle
  struct{
    float abs;                // absorption
    float sca;                // scattering
  } z [MAXLYS];
};

struct dats{
  union{
    struct{
      int hidx;
      float r[3];
    };
    struct{
      DOM * oms;
      name * names;
    };
  };
  int type;   // 0=cascade/1=flasher/2=flasher 45/3=laser up/4=laser down
  int nfla;   // wv. bin of flasher
  int hnum;   // size of hits buffer
  int size;   // size of kurt table
  int rsize;  // count of multipliers
  int gsize;  // count of initialized OMs

  float g;                    // <cos(scattering angle)>
  float ocv;                  // 1 / speed of light in vacuum
  float dh, hdh, rdh, hmin;   // step, 1/step, min and max depth

#ifdef ACCL2
  int cn[3];
  float cl[3], crst[3];
  short cidx[cx][cy][cz];
  short carr[MAXGEO];
#ifdef ROMB
  float cb[2][2];
#endif
#endif
  hit * hits;
  unsigned int * rm;
  unsigned long long * rs;
  photon * pz;
  ices * w[WNUM];
} d;

static const float zoff=1948.07;

struct ini{
#ifdef ACCL2
  float ctr(float r[], int m){
#ifdef ROMB
    if(m<2) return d.cb[0][m]*r[0]+d.cb[1][m]*r[1];
    else return r[2];
#else
    return r[m];
#endif
  }
#endif

  ini(){

    { // initialize random numbers
      int size;

      ifstream inFile;
      inFile.open("rnd.txt", ifstream::in);
      vector<unsigned int> rx;

      if(!inFile.fail()){
	unsigned int a;
	unsigned long long b,c;
	while(inFile>>a>>b>>c) rx.push_back(a);
	if(rx.size()<1){ cerr<<"File rnd.txt did not contain valid data"<<endl; exit(1); }
	inFile.close();
      }
      else{ cerr<<"Could not open file rnd.txt"<<endl; exit(1); }
      size=rx.size();
      cerr<<"Loaded "<<size<<" random multipliers"<<endl;

      d.rm = new unsigned int[size];
      d.rs = new unsigned long long[size];

      union{
	unsigned long long da;
	struct{
	  unsigned int ax;
	  unsigned int dx;
	};
      } s;

      s.ax=362436069;
      s.dx=1234567;

      for(int i=0; i<size; i++){
	d.rsize=size;
	d.rm[i]=rx[i];
	d.rs[i]=s.da;
      }
    }

    { // initialize ice parameters
      int size;                   // size of kurt table
      float g;                    // <cos(scattering angle)>
      float dh, hdh, rdh, hmin;   // step, 1/step, min and max depth

      ifstream inFile;

      bool flag=true;
      inFile.open("wv.dat", ifstream::in);
      vector<float> wx, wy;

      if(!inFile.fail()){
	float xa, ya, xo=0, yo=0;
	int num=0;
	while(inFile>>xa>>ya){
	  if(( xa<0 || 1<xa ) || num==0 && xa!=0 || num>0 && ( xa<=xo || ya<=yo )){ flag=false; break; }
	  wx.push_back(xa); wy.push_back(ya);
	  xo=xa; yo=ya; num++;
	}
	if(xo!=1 || wx.size()<2) flag=false;
	inFile.close();
	if(flag){ cerr<<"Loaded "<<wx.size()<<" wavelenth points"<<endl; }
	else{ cerr<<"File wv.dat did not contain valid data"<<endl; exit(1); }
      }
      else{ cerr<<"Could not open file wv.dat"<<endl; exit(1); }

      flag=true;
      float A, B, D, E, a, k;
      float Ae, Be, De, Ee, ae, ke;
      vector<float> dp, be, ba, td;

      inFile.open("icemodel.par", ifstream::in);

      if(flag=!inFile.fail()){
	if(flag) flag=(inFile >> a >> ae);
	if(flag) flag=(inFile >> k >> ke);
	if(flag) flag=(inFile >> A >> Ae);
	if(flag) flag=(inFile >> B >> Be);
	if(flag) flag=(inFile >> D >> De);
	if(flag) flag=(inFile >> E >> Ee);
	if(!flag) cerr << "File icemodel.par found, but is corrupt" << endl;
	inFile.close();
	if(!flag) exit(1);
      }
      else{ cerr << "File icemodel.par was not found" << endl; exit(1); }

      inFile.open("icemodel.dat", ifstream::in);
      if(!inFile.fail()){
	size=0;
	float dpa, bea, baa, tda;
	while(inFile >> dpa >> bea >> baa >> tda){
	  dp.push_back(dpa);
	  be.push_back(bea);
	  ba.push_back(baa);
	  td.push_back(tda);
	  size++;
	}
	inFile.close();
	if(size<2){ cerr << "File icemodel.dat found, but is corrupt" << endl; exit(1); }
      }
      else{ cerr << "File icemodel.dat was not found" << endl; exit(1); }

      dh=dp[1]-dp[0];
      if(dh<=0){ cerr << "Ice table does not use increasing depth spacing" << endl; exit(1); }

      for(int i=0; i<size; i++) if(i>0) if(fabsf(dp[i]-dp[i-1]-dh)>dh*xx){
	cerr << "Ice table does not use uniform depth spacing" << endl; exit(1);
      }
      cerr<<"Loaded "<<size<<" ice layers"<<endl;

      if(size>MAXLYS){
	cerr<<"Error: too many layers ("<<size<<"), truncating to "<<MAXLYS<<endl;
	size=MAXLYS;
      }

      hdh=dh/2; rdh=1/dh; hmin=zoff-dp[size-1];
      g=0.8;
      {
	d.g=g;
	d.size=size;
	d.dh=dh;
	d.hdh=hdh;
	d.rdh=rdh;
	d.hmin=hmin;
      }

      int nfla;
      float wfla=405.f, fmin;
      {
	char * WFLA=getenv("WFLA");
	float aux=WFLA==NULL?0:atof(WFLA);
	if(aux>0) wfla=aux;
      }

      for(int n=0; n<WNUM; n++){
	float p=(n+0.5)/WNUM;
	int m=0; while(wx[++m]<p);
	float wva=(wy[m-1]*(wx[m]-p)+wy[m]*(p-wx[m-1]))/(wx[m]-wx[m-1]);

	float faux=fabsf(wva-wfla);
	if(n==0 || faux<fmin){ nfla=n; fmin=faux; }

	float wv0=400;
	float l_a=pow(wva/wv0, -a);
	float l_k=pow(wva, -k);
	float ABl=A*exp(-B/wva);

	ices * w = new ices;
	d.w[n] = w;

	for(int i=0; i<size; i++){
	  int j=size-1-i;
	  w->z[i].sca=be[j]*l_a/(1-g); w->z[i].abs=(D*ba[j]+E)*l_k+ABl*(1+0.01*td[j]);
	  if(w->z[i].sca<=0 || w->z[i].abs<=0){ cerr << "Invalid value of ice parameter, cannot proceed" << endl; exit(1); }
	}

	float wv=wva*1.e-3;
	float wv2=wv*wv;
	float wv3=wv*wv2;
	float wv4=wv*wv3;
	float np=1.55749-1.57988*wv+3.99993*wv2-4.68271*wv3+2.09354*wv4;
	float ng=np*(1+0.227106-0.954648*wv+1.42568*wv2-0.711832*wv3);
	float c=0.299792458; d.ocv=1/c; w->wvl=wva; w->ocm=ng/c;
	w->coschr=1/np; w->sinchr=sqrt(1-w->coschr*w->coschr);
      }
      d.nfla=nfla;
    }

    { // initialize geometry
      ifstream inFile;
      inFile.open("geo-f2k", ifstream::in);
      vector<DOM> oms;
      vector<name> names;

      if(!inFile.fail()){
	DOM dom;
	name n;
	string mbid;
	int om, str;
	unsigned long long omid;
	while(inFile>>mbid>>hex>>omid>>dec>>dom.r[0]>>dom.r[1]>>dom.r[2]>>str>>om) if(str>0 && om>=1 && om<=60){
	  dom.r[2]+=zoff;
	  n.dom=om;
	  n.str=str;
	  oms.push_back(dom);
	  names.push_back(n);
	}
	inFile.close();
      }

      int gsize = oms.size();
      if(gsize>MAXGEO){
	cerr<<"Error: too many OMs ("<<gsize<<"), truncating to "<<MAXGEO<<endl;
	gsize=MAXGEO;
      }

      d.oms = new DOM[gsize];
      d.names = new name[gsize];
      for(int n=0; n<gsize; n++){ d.oms[n]=oms[n]; d.names[n]=names[n]; }

#ifdef ACCL2
      d.cn[0]=cx;
      d.cn[1]=cy;
      d.cn[2]=cz;

#ifdef ROMB
      const float cv=fpi/180;
      float bv[2][2];
      bv[0][0]=cos(cv*dir1);
      bv[0][1]=sin(cv*dir1);
      bv[1][0]=cos(cv*dir2);
      bv[1][1]=sin(cv*dir2);

      float det=bv[0][0]*bv[1][1]-bv[0][1]*bv[1][0];
      d.cb[0][0]=bv[1][1]/det;
      d.cb[0][1]=-bv[0][1]/det;
      d.cb[1][0]=-bv[1][0]/det;
      d.cb[1][1]=bv[0][0]/det;
#endif

      float cl[3]={0,0,0}, ch[3]={0,0,0}, crst[3];

      vector<int> cells[cx][cy][cz];
      for(int i=0; i<cx; i++) for(int j=0; j<cy; j++) for(int k=0; k<cz; k++)
	cells[i][j][k].erase(cells[i][j][k].begin(), cells[i][j][k].end());

      for(int n=0; n<gsize; n++) for(int m=0; m<3; m++){
	if(n==0 || ctr(oms[n].r, m)<cl[m]) cl[m]=ctr(oms[n].r, m);
	if(n==0 || ctr(oms[n].r, m)>ch[m]) ch[m]=ctr(oms[n].r, m);
      }
      for(int m=0; m<3; m++){
	float diff=ch[m]-cl[m], s=R*(d.cn[m]-1);
	if(diff<2*s){
	  cerr<<"Error in cell["<<m<<"] initialization"<<endl;
	  double ave=(cl[m]+ch[m])/2;
	  cl[m]=ave-s; ch[m]=ave+s; diff=2*s;
	}
	crst[m]=(d.cn[m]-1)/diff;
      }

      for(int n=0; n<gsize; n++){
	int i[3];
	for(int m=0; m<3; m++){
	  i[m]=lroundf((ctr(oms[n].r, m)-cl[m])*crst[m]);
	  if(i[m]<0 && i[m]>=d.cn[m]){ cerr<<"Error in cell initialization"<<endl; exit(1); }
	}
	cells[i[0]][i[1]][i[2]].push_back(n);
      }

      for(int m=0; m<3; m++){ d.cl[m]=cl[m]; d.crst[m]=crst[m]; }

      int m=0;
      for(int i=0; i<cx; i++) for(int j=0; j<cy; j++) for(int k=0; k<cz; k++){
	int size=cells[i][j][k].size();
	d.cidx[i][j][k]=size>0?m:-1;
	for(int l=0; l<size; l++, m++) d.carr[m]=cells[i][j][k][l];
	if(size>0) d.carr[m-1]|=0x8000;
      }
#endif
      d.gsize=gsize;

      cerr<<"Loaded "<<gsize<<" DOMs"<<endl;
    }

    {
      d.hnum=HNUM;
      d.hits = new hit[d.hnum];
    }
  }

  ~ini(){   // empties booked arrays and tables
    delete d.rm;
    delete d.rs;
    for(int n=0; n<WNUM; n++) delete d.w[n];
    delete d.oms;
    delete d.names;
    delete d.hits;
  }
} m;