#define OFLA         // omit the flasher DOM
#define ROMB         // use rhomb cells aligned with the array

#define ASENS        // enable angular sensitivity
#define TILT         // enable tilted ice layers

#define MKOW         // use Marek Kowalski's photon yield parametrization
#define LONG         // simulate longitudinal cascade development
#define ANGW         // smear cherenkov cone due to shower development

#ifdef ASENS
#define ANUM 11      // number of coefficients in the angular sensitivity curve
#endif

#ifdef TILT
#define LMAX 6
#define LYRS 170
#endif

#ifdef XCPU
#define MLTP 0
#define WNUM 40
#define OVER 1
#define NTHR 4
#else
#define MLTP 30      // number of multiprocessors
#define WNUM 32      // optimized for 30 multiprocessors
#define OVER 10      // size of photon bunches along the muon track
#define NTHR 384     // NTHR*NBLK should not exceed the count of different random number multipliers
#endif

#define NPHO 1024    // maximum number of photons propagated by one thread
#define HNUM 1048576 // size of the output hit buffer, must hold hits from up to NPHO*NTHR*NBLK photons

#define MAXLYS 180   // maximum number of ice layers
#define MAXGEO 5200  // maximum number of OMs

#define OVR 5        // over-R: DOM radius "oversize" scaling factor

static const float xx=1.e-5;
static const float R=OVR*0.16510;
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 ikey{
  int str, dom;
};

struct OM:DOM,ikey{};
vector<OM> i3oms;

struct gless{
  bool operator()(const OM & lhs, const OM & rhs) const {
    return lhs.str == rhs.str ? lhs.dom < rhs.dom : lhs.str < rhs.str;
  }
};

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

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

#ifdef ROMB
static const float dir1=9.3, dir2=129.3;
static const unsigned int cx=21, cy=19, nstr=94;
#else
static const unsigned int cx=13, cy=13, nstr=94;
#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 line{
  short n, max;
  float x, y, r;
  float h, d;
  float dl, dh;
};

struct dats{
  union{
    struct{
      unsigned int hidx;
      float r[3];
    };
    struct{
      DOM * oms;
      name * names;
    };
  };
#ifndef XCPU
  unsigned int dt;  // kernel time clocks
  unsigned int ab;  // if TOT was abnormal
#endif
  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, step/2, 1/step, and min depth

  int cn[2];
  float cl[2], crst[2];

  unsigned char is[cx][cy];
  unsigned char ls[nstr];
  line sc[nstr];
  float rx;

  float fldr; // horizontal direction of the flasher led #1

#ifdef ASENS
  float eff;  // OM efficiency correction
  float mas;  // maximum angular sensitivity
  float s[ANUM]; // ang. sens. coefficients
#define EFF d.eff*d.mas
#else
#define EFF 1
#endif

#ifdef ROMB
  float cb[2][2];
#endif
#ifdef OFLA
  short fla;
#endif
#ifdef TILT
  int lnum, lpts, l50;
  float lmin, lrdz, r50;
  float lnx, lny;
  float lr[LMAX];
  float lp[LMAX][LYRS];
#endif

  hit * hits;
  unsigned int * rm;
  unsigned long long * rs;
  photon * pz;
  ices * w[WNUM];
} d;

unsigned char sname(int n){
  name & s = d.names[n];
  return s.str>78&&s.dom>10?s.str+10:s.str;
}

static const float zoff=1948.07;

struct ini{
  float ctr(line & s, int m){
#ifdef ROMB
    return d.cb[0][m]*s.x+d.cb[1][m]*s.y;
#else
    return m==0?s.x:s.y;
#endif
  }

  ini(){
    string dir("");
    {
      char * env = getenv("PPCTABLESDIR");
      if(env!=NULL) dir=string(env)+"/";
#ifdef XLIB
      if(env==NULL) dir=getenv("I3_BUILD")+string("/ppc/resources/");
      cerr<<endl<<"Configuring in "<<dir<<endl;
#endif
    }

#ifdef ASENS
    {
      d.eff=1, d.mas=1, d.s[0]=1;
      for(int i=1; i<ANUM; i++) d.s[i]=0;

      ifstream inFile((dir+"as.dat").c_str(), ifstream::in);
      if(!inFile.fail()){
	int n=0;
	float aux;
	if(inFile >> aux) d.eff=aux, n++;
	if(inFile >> aux) d.mas=aux, n++;
	for(int i=0; i<ANUM; i++) if(inFile >> aux) d.s[i]=aux, n++;
	if(n>0) cerr<<"Loaded "<<n<<" angsens coefficients"<<endl;
	else{ cerr<<"File as.dat did not contain valid data"<<endl; exit(1); }
      }
      else{ cerr<<"Could not open file as.dat"<<endl; exit(1); }
    }
#endif

#ifdef TILT
    {
      d.lnum=0;
      const float cv=fpi/180, th=225;
      d.lnx=cos(cv*th), d.lny=sin(cv*th);

      ifstream inFile((dir+"tilt.par").c_str(), ifstream::in);
      if(!inFile.fail()){
	int str;
	float aux;
	vector<float> lr;
	while(inFile >> str >> aux){ if(aux==0) d.l50=str; lr.push_back(aux); }
	inFile.close();

	int size=lr.size();
	if(size>LMAX){ cerr << "File tilt.par defines too many dust maps" << endl; exit(1); }
	for(int i=1; i<size; i++) if(lr[i]<lr[i-1]) { cerr << "Tilt map does not use increasing range order" << endl; exit(1); };
	for(int i=0; i<size; i++) d.lr[i]=lr[i];

	ifstream inFile((dir+"tilt.dat").c_str(), ifstream::in);
	if(!inFile.fail()){
	  d.lnum=size;
	  float depth, pts[d.lnum];
	  vector<float> ds, lp[d.lnum];
	  while(inFile >> depth){
	    int i=0;
	    while(inFile >> pts[i++]) if(i>=d.lnum) break;
	    if(i!=d.lnum) break;
	    if(pts[0]!=0){ cerr << "Format error in file tilt.dat: expecting 0" << endl; exit(1); }
	    ds.push_back(depth);
	    for(i=0; i<d.lnum; i++) lp[i].push_back(pts[i]);
	  }
	  inFile.close();

	  int size=ds.size();
	  if(size>LYRS){ cerr << "File tilt.dat defines too many map points" << endl; exit(1); }
	  for(int i=1; i<size; i++) if(ds[i]<ds[i-1]) { cerr << "Tilt map does not use increasing depth order" << endl; exit(1); };
	  for(int i=0; i<d.lnum; i++) for(int j=0; j<size; j++) d.lp[i][j]=lp[i][size-1-j];
	  d.lpts=size;

	  if(size<2) d.lnum=0;
	  else{
	    float lmin=ds[0], lmax=ds[size-1];
	    d.lmin=zoff-lmax; d.lrdz=(d.lpts-1)/(lmax-lmin);
	  }
	}

	if(d.lnum>0) cerr<<"Loaded "<<d.lnum<<"x"<<d.lpts<<" dust layer points"<<endl;
      }
    }
#endif

    { // initialize random numbers
      int size;
      vector<unsigned int> rx;

      ifstream inFile((dir+"rnd.txt").c_str(), ifstream::in);
      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, step/2, 1/step, and min depth

      bool flag=true;
      vector<float> wx, wy;

      {
	ifstream inFile((dir+"wv.dat").c_str(), ifstream::in);
	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;

      {
	ifstream inFile((dir+"icemodel.par").c_str(), 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); }
      }

      {
	ifstream inFile((dir+"icemodel.dat").c_str(), 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.5f)/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;
#ifdef OFLA
      d.fla=-1;
#endif

      {
	char * FLDR=getenv("FLDR");
	d.fldr=FLDR==NULL?-1:atof(FLDR);
      }
    }

    {
      d.hnum=HNUM;
    }

    {
      d.oms=NULL;
      d.names=NULL;
    }

#ifndef XLIB
    {
      ifstream inFile((dir+"geo-f2k").c_str(), ifstream::in);
      if(!inFile.fail()){
	OM om;
	string mbid;
	unsigned long long omid;
	while(inFile>>mbid>>hex>>omid>>dec>>om.r[0]>>om.r[1]>>om.r[2]>>om.str>>om.dom){
	  om.r[2]+=zoff;
	  i3oms.push_back(om);
	}
	inFile.close();
      }

      geo();
    }
#endif
  }

  void geo(){ // initialize geometry
    {
      {
	if(d.oms!=NULL) delete d.oms;
	if(d.names!=NULL) delete d.names;
      }

      vector<DOM> oms;
      vector<name> names;

      sort(i3oms.begin(), i3oms.end(), gless());
      for(vector<OM>::const_iterator i=i3oms.begin(); i!=i3oms.end(); ++i) if(i->str>0 && i->dom>=1 && i->dom<=60){
	name n;
	n.dom=i->dom;
	n.str=i->str;
	oms.push_back(*i);
	names.push_back(n);
      }

      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]; }

      d.gsize=gsize;
    }

    map<unsigned char, short> num;
    {
      map<unsigned char, float> l;
      map<unsigned char, line> sc;
      for(int n=0; n<d.gsize; n++){
	unsigned char str=sname(n);
	line & s = sc[str];
	DOM & om = d.oms[n];
	if(num.find(str)==num.end()){
	  l[str]=om.r[2];
	  s.h=om.r[2];
	  s.n=n;
	}
	else{
	  if(l[str]>om.r[2]) l[str]=om.r[2];
	  if(s.h<om.r[2]) s.h=om.r[2];
	  if(s.n>n) s.n=n;
	}
	num[str]++; s.x+=om.r[0], s.y+=om.r[1];
      }
      if(sc.size()>nstr){ cerr<<"Number of strings exceeds capacity of "<<nstr<<endl; exit(1); }

      for(map<unsigned char, short>::iterator i=num.begin(); i!=num.end(); ++i){
	unsigned char str=i->first, n=i->second;
	line & s = sc[str];
	float d=s.h-l[str];
	if(n<2 || d<=0){ cerr<<"Cannot estimate the spacing along string "<<str<<endl; exit(1); }
	s.x/=n, s.y/=n, s.r=0; s.d=(n-1)/d; s.dl=0; s.dh=0;
      }

      for(int n=0; n<d.gsize; n++){
	unsigned char str=sname(n);
	line & s = sc[str];
	DOM & om = d.oms[n];
#ifdef TILT
	if(str==d.l50) d.r50=d.lnx*s.x+d.lny*s.y;
#endif

	float dx=s.x-om.r[0], dy=s.y-om.r[1];
	float dr=dx*dx+dy*dy; if(dr>s.r) s.r=dr;

	float dz=om.r[2]-(s.h+(s.n-n)/s.d);
	if(s.dl>dz) s.dl=dz; if(s.dh<dz) s.dh=dz;
      }

      d.rx=0;
      int n=0;
      for(map<unsigned char, short>::iterator i=num.begin(); i!=num.end(); ++i, n++){
	unsigned char str=i->first;
	line & s = sc[str];
	s.max=i->second-1;
	i->second=n;
	s.r=R+sqrt(s.r);
	if(d.rx<s.r) d.rx=s.r;
	s.dl-=R, s.dh+=R;
	d.sc[n]=s;
      }
    }

    float sin12=0;
    {
#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;

      for(int i=0; i<2; i++) sin12+=bv[0][i]*bv[1][i];
#endif
      sin12=sqrt(1-sin12*sin12);
      d.rx/=sin12;
    }

    map<unsigned char, int> cells[cx][cy];
    {
      float cl[2]={0,0}, ch[2]={0,0}, crst[2];

      int n=0;
      for(map<unsigned char, short>::iterator i=num.begin(); i!=num.end(); ++i, n++){
	line & s = d.sc[i->second];
	for(int m=0; m<2; m++){
	  if(n==0 || ctr(s, m)<cl[m]) cl[m]=ctr(s, m);
	  if(n==0 || ctr(s, m)>ch[m]) ch[m]=ctr(s, m);
	}
      }

      d.cn[0]=cx;
      d.cn[1]=cy;

      for(int m=0; m<2; m++){
	float diff=ch[m]-cl[m];
#ifdef ROMB
	d.cn[m]=min(d.cn[m], max(2, 2*(int)lroundf(diff/125)+1));
#endif
	float s=R*(d.cn[m]-1);
	if(diff<2*s){
	  cerr<<"Error in cell["<<m<<"] initialization"<<endl;
	  float ave=(cl[m]+ch[m])/2;
	  cl[m]=ave-s; ch[m]=ave+s; diff=2*s;
	}
	crst[m]=(d.cn[m]-1)/diff;
      }

      bool flag=true;
      for(map<unsigned char, short>::iterator i=num.begin(); i!=num.end(); ++i){
	line & s = d.sc[i->second];
	int n[2];
	for(int m=0; m<2; m++){
	  n[m]=lroundf((ctr(s, m)-cl[m])*crst[m]);
	  if(n[m]<0 && n[m]>=d.cn[m]){ cerr<<"Error in cell initialization"<<endl; exit(1); }

	  float d1=fabsf(ctr(s, m)-(cl[m]+(n[m]-0.5f)/crst[m]));
	  float d2=fabsf(ctr(s, m)-(cl[m]+(n[m]+0.5f)/crst[m]));
	  float d=min(d1, d2)*sin12-s.r;
	  if(d<0){ flag=false; cerr<<"Warning: string "<<(int)i->first<<" too close to cell boundary"<<endl; }
	}

	cells[n[0]][n[1]][i->first]++;
      }
      if(flag) d.rx=0;

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

    {
      unsigned int pos=0;
      for(int i=0; i<d.cn[0]; i++) for(int j=0; j<d.cn[1]; j++){
	map<unsigned char, int> & c = cells[i][j];

	if(c.size()>0){
	  d.is[i][j]=pos;
	  for(map<unsigned char, int>::const_iterator n=c.begin(); n!=c.end(); ++n){
	    if(pos==nstr){ cerr<<"Number of string cells exceeds capacity of "<<nstr<<endl; exit(1); }
	    d.ls[pos++]=num[n->first];
	  }
	  d.ls[pos-1]|=0x80;
	}
	else d.is[i][j]=0x80;
      }
    }

    cerr<<"Loaded "<<d.gsize<<" DOMs ("<<d.cn[0]<<"x"<<d.cn[1]<<")"<<endl;
  }

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