static const float xx=1.e-5;
static const float zoff=1948.07;

static const float OVR=5;
static const float R=OVR*0.178;

static const float ppm=2107.84/(OVR*OVR);
static const float crf[4]={ ppm*4.889*0.894, // em shower
			    ppm*4.076*0.860, // hdr m/GeV
			    ppm*1.1720,      // bare
			    ppm*0.0324 };    // muon

static const int wnum=40;
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 [];
} * iced[wnum];


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

  name(){
    f=1.0;
  }
};

struct DOM{
  name t;
  float r[3];

  DOM(){
    for(int i=0; i<3; i++) r[i]=0;
  }
};

vector<DOM> oms;

#ifdef ACCL2
  static const int cx=10, cy=10, cz=10;
  int cn[3]={cx, cy, cz};
  vector<int> cells[cx][cy][cz];
  float cl[3], ch[3], crst[3];
#endif

struct dats{
  unsigned int rsize, ml;
  union{
    unsigned long long da;
    struct{
      unsigned int ax;
      unsigned int dx;
    };
  };

  unsigned int * rm;
  unsigned long long * rs;

  float xrnd(){
    union{
      unsigned int tmp;
      float xrd;
    };
    do{
      da = ml * (unsigned long long) ax + dx;
      tmp = ax >> 9;
    } while(tmp==0);
    tmp |= 0x3f800000;
    return xrd-1;
  }

  void setr(unsigned int s){
    int idx=s%rsize;
    ml=rm[idx];
    da=rs[idx];
  }

  void endr(unsigned int s){
    int idx=s%rsize;
    rs[idx]=da;
  }

  int size;                   // size of kurt table
  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
  ices * w;

  // get distance for overburden x
  float aX(float x, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float ax, ai, ah;

    if(nz>0){
      ah=h+hdh;
      ai=(nz*x-(ah-z)*w->z[i].abs)*rdh;
      for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->z[++j].abs;
    }
    else{
      ah=h-hdh;
      ai=(nz*x-(ah-z)*w->z[i].abs)*rdh;
      for(j=i; j>0 && ai<0; ah-=dh) ai+=w->z[--j].abs;
    }

    if(i==j) ax=x/w->z[i].abs;
    else ax=(ai*dh/w->z[j].abs+ah-z)/nz;

    return ax;
  }

  // get distance for overburden x
  float sX(float x, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float sx, ai, ah;

    if(nz>0){
      ah=h+hdh;
      ai=(nz*x-(ah-z)*w->z[i].sca)*rdh;
      for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->z[++j].sca;
    }
    else{
      ah=h-hdh;
      ai=(nz*x-(ah-z)*w->z[i].sca)*rdh;
      for(j=i; j>0 && ai<0; ah-=dh) ai+=w->z[--j].sca;
    }

    if(i==j) sx=x/w->z[i].sca;
    else sx=(ai*dh/w->z[j].sca+ah-z)/nz;

    return sx;
  }

  // get overburden for distance a
  float xA(float a, float z, float nz){
    int i=lroundf((z-hmin)*rdh), j;
    if(i<0) i=0; else if(i>=size) i=size-1;

    float h=hmin+i*dh;
    float y=z+nz*a;

    j=lroundf((y-hmin)*rdh);
    if(j<0) j=0; else if(j>=size) j=size-1;

    float xa;

    if(i==j) xa=a*w->z[i].abs; 
    else{
      float g=hmin+j*dh;
      if(nz>0){
	xa=(y-(g-hdh))*w->z[j].abs;
	while(--j>i) xa+=hdh*w->z[j].abs;
	xa+=(h+hdh-z)*w->z[i].abs;
      }
      else{
	xa=(y-(g+hdh))*w->z[j].abs;
	while(++j<i) xa-=hdh*w->z[j].abs;
	xa+=(h-hdh-z)*w->z[i].abs;
      }
      xa/=nz;
    }

    return xa;
  }

  float t, r[3], n[3];
  int type;
  float xt, xr[3], xn[3];

  void photon(){
    t=xt;
    for(int i=0; i<3; i++){
      r[i]=xr[i]; n[i]=xn[i];
    }

    if(type==1) rotate(w->coschr, w->sinchr);
    else{
      static const float rms=9.2*sqrtf(3)*M_PI/180;
      float xi=xrnd();
      xi=rms*(2*xi-1);
      float nz=sinf(xi), np=cosf(xi);

      xi=xrnd();
      xi*=2*M_PI;
      n[0]=np*cosf(xi);
      n[1]=np*sinf(xi);
      n[2]=nz;
    }
  }

  void rotate(float cs, float si){
    float u[3];
    {
      float p1[3], p2[3];
      {
	int i0=0;
	float n0=n[0];
	for(int i=1; i<3; i++) if(fabsf(n[i])>fabsf(n0)){ i0=i; n0=n[i]; }
	int i1=(i0+1)%3, i2=(i0+2)%3;
	float n1=n[i1], n2=n[i2], nq=n0*n0;
	float r1=1/sqrtf(nq+n1*n1), r2=1/sqrtf(nq+n2*n2);

	p1[i0]=-n1*r1; p1[i1]=n0*r1; p1[i2]=0;
	p2[i0]=-n2*r2; p2[i1]=0; p2[i2]=n0*r2;
      }
      {
	float q1[3], q2[3], r1=0, r2=0;
	for(int i=0; i<3; i++){
	  float a1=p1[i]-p2[i]; q1[i]=a1; r1+=a1*a1;
	  float a2=p1[i]+p2[i]; q2[i]=a2; r2+=a2*a2;
	}

	r1=1/sqrtf(r1); r2=1/sqrtf(r2);

	for(int i=0; i<3; i++){
	  p1[i]=q1[i]*r1;
	  p2[i]=q2[i]*r2;
	}
      }

      {
	float zi=xrnd(); zi*=2*M_PI;
	float px=cosf(zi), py=sinf(zi);
	for(int i=0; i<3; i++) u[i]=px*p1[i]+py*p2[i];
      }
    }

    {
      float r=0;
      for(int i=0; i<3; i++){
	float a=n[i]*cs+u[i]*si;
	n[i]=a; r+=a*a;
      }
      {
	r=1/sqrtf(r);
	for(int i=0; i<3; i++) n[i]*=r;
      }
    }
  }

  void rotate(){
    float xi=xrnd();
    xi=2*xi-1;
    if(g!=0){
      float g2=g*g;
      float ga=(1-g2)/(1+g*xi);
      xi=(1+g2-ga*ga)/(2*g);
    }
    if(xi>1) xi=1; else if(xi<-1) xi=-1;

    float si=sqrtf(1-xi*xi);
    rotate(xi, si);
  }

  void advance(float a){
    t+=a*w->ocm;
    for(int i=0; i<3; i++) r[i]+=n[i]*a;
  }

  void propagate(){
    photon();
    float TOT=-logf(xrnd());
    while(true){
      float SCA=-logf(xrnd());
      float sca=sX(SCA, r[2], n[2]);
      float tot=aX(TOT, r[2], n[2]);
      if(tot<sca) sca=tot;

      name om=sphere(sca);
      TOT-=xA(sca, r[2], n[2]);
      advance(sca);

      if(om.k!=0){
	cout<<"HIT "<<(int)om.str<<" "<<(int)om.dom<<" "<<t<<" "<<n[2]<<endl;
	break;
      }

      if(TOT<xx) break;
      rotate();
    }
  }

  name sphere(float &sca){
    name imin;

#ifdef ACCL0
    float dx[3], dn[3];
    for(int m=0; m<3; m++){
      float ini=r[m];
      float fin=ini+sca*n[m];
      if(n[m]>0){
	dn[m]=ini;
	dx[m]=fin;
      }
      else{
	dn[m]=fin;
	dx[m]=ini;
      }
      dn[m]-=R;
      dx[m]+=R;
    }
#endif

#ifdef ACCL2
    int xl[3], xh[3];
    for(int m=0; m<3; m++){
      xl[m]=min(max((int) lroundf((dn[m]-cl[m])*crst[m]), 0), cn[m]);
      xh[m]=max(min((int) lroundf((dx[m]-cl[m])*crst[m]), cn[m]-1), -1);
    }
    for(int i=xl[0]; i<=xh[0]; i++) for(int j=xl[1]; j<=xh[1]; j++) for(int k=xl[2]; k<=xh[2]; k++) for(int l=0; l<cells[i][j][k].size(); l++){
      const DOM & om=oms[cells[i][j][k][l]];
#else
    for(int l=0; l<oms.size(); l++){
      const DOM & om=oms[l];
#endif
#ifdef ACCL1
      bool pass=false;
      for(int i=0; i<3; i++) if(om.r[i]<dn[i] || om.r[i]>dx[i]){ pass=true; break; }
      if(pass) continue;
#endif

      {
	float b=0, c=0;
	for(int i=0; i<3; i++){
	  float dr=om.r[i]-r[i];
	  b+=n[i]*dr; c+=dr*dr;
	}
	c-=R*R;
	float D=b*b-c, h=-1;
	if(D>0) h=b-sqrtf(D);

	if(h>0 && h<sca){ imin=om.t; sca=h; }
      }
    }

    return imin;
  }
} d;


struct ini{
  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); }
      }
      else{ cerr<<"Could not open file rnd.txt"<<endl; exit(1); }
      size=rx.size();
      cerr<<"Loaded "<<size<<" random multipliers"<<endl;

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

      d.ax=362436069;
      d.dx=1234567;

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

    { // initialize ice parameters
      ices *w;
      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;

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

      int nfla;
      float wfla=405, fmin;
      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);

	w = (ices *) new float[2*size+sizeof(ices)/sizeof(float)];
	iced[n] = w;

	d.g=g;
	d.size=size;
	d.dh=dh;
	d.hdh=hdh;
	d.rdh=rdh;
	d.hmin=hmin;

	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->ocm=ng/c;
	w->coschr=1/np; w->sinchr=sqrt(1-w->coschr*w->coschr);
      }
      d.w=iced[nfla];
    }

    { // initialize geometry
      ifstream inFile;
      inFile.open("geo-f2k", ifstream::in);

      if(!inFile.fail()){
	DOM dom;
	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;
	  dom.t.dom=om;
	  dom.t.str=str;
	  oms.push_back(dom);
	}
	inFile.close();
      }

#ifdef ACCL2
      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<oms.size(); n++) for(int m=0; m<3; m++){
	if(n==0 || oms[n].r[m]<cl[m]) cl[m]=oms[n].r[m];
	if(n==0 || oms[n].r[m]>ch[m]) ch[m]=oms[n].r[m];
      }
      for(int m=0; m<3; m++){
	float diff=ch[m]-cl[m];
	if(diff<R){ cerr<<"Error in cell initialization"<<endl; exit(1); }
	crst[m]=(cn[m]-1)/diff;
      }

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

      cerr<<"Loaded "<<oms.size()<<" DOMs"<<endl;
    }

  }

  void destroy(){   // empties booked arrays and tables
    delete d.rm;
    delete d.rs;
    for(int n=0; n<wnum; n++) delete iced[n];
  }

  ~ini(){
    destroy();
  }

  void next(float rnd){
    d.w=iced[lroundf(wnum*rnd-0.5)];
  }
} m;