#include <cmath>
#include <cstdlib>
#include <iostream>
#include <map>
#include <vector>
#include <fstream>
#include <cstring>

#ifdef SINGLE
#define double float
#define cos cosf
#define sin sinf
#define log logf
#define sqrt sqrtf
#define ceil ceilf
#define floor floorf
#endif
#define XACCL

#ifdef ROTASM
float asm_nq[4], asm_cos[2];
float asm_zsn=0.6, asm_zcs=0.8;
double asm_zi=0.2, asm_z2=2.0, asm_dbl;
char asm_sse[512];
char *asm_fmt="%f\n";
extern "C" void rot();
#endif
namespace ppc{
  using namespace std;

  /*
  int floor(double x){
    return lround(x-0.5);
  }

  int ceil(double x){
    return lround(x+0.5);
  }
  */

#ifndef ICETRAY
  const string dir("");
#else
  struct begin{ begin(){ cerr<<endl; } } b;
  const string dir(getenv("I3_BUILD")+string("/ppc/resources/"));
#endif

#include "ice.cxx"

  double xrnd(){
    const double a=1./RAND_MAX;
    double rnd=0;
    while(rnd==0) rnd=rand();
    return a*rnd;
  }

  struct KEY{
    int str, om;

    KEY(){
      str=-1; om=-1;
    }

    bool operator< (const KEY& rhs) const{
      return str < rhs.str || (str == rhs.str && om < rhs.om);
    }
  };

  struct STR{
    bool zf;
    int l, h;    // low (1) and high (60) DOM numbers
    double zn, zx, sp;   // tolerances, number of DOMs
    double zl, zh, st;   // min and max z, 1/step size
    double rl[3], rh[3]; // box around each string (+R)

    STR(){
      zf=true;
      sp=-1;
    }
  };

  map<int,STR> strs;

#ifdef XACCL
  static const int ncel=10;
  vector<int> cells[ncel][ncel];
  double rl[3], rh[3], rst[2];
#endif

  static const double zoff=-1948.07;
  static const double OVR=5;
  static const double R=OVR*0.178;

  struct DOM{
    double r[3];
  };
  map<KEY,DOM> doms;

#ifndef ICETRAY
  void load(){
    ifstream inFile;
    inFile.open((dir+"geo-f2k").c_str(), ifstream::in);

    if(!inFile.fail()){
      KEY key;
      DOM dom;
      string id;
      unsigned long long domid;
      while(inFile>>id>>hex>>domid>>dec>>dom.r[0]>>dom.r[1]>>dom.r[2]>>key.str>>key.om)
#else
  double angsens(double x){ // hole ice
    double s[11]={0.32752, 0.66696, 0.35389, -1.4107, -1.7738, 4.3302, 5.7151, -5.5262, -7.9000, 2.0155, 3.3965};
    double sum=s[0], y=1;
    for(int i=1; i<11; i++){ y*=x; sum+=s[i]*y; }
    return sum;
  }

  int hitnum=0;
  pair<int,unsigned long long> mcpids;
  I3MCHitSeriesMapPtr mchits(new I3MCHitSeriesMap);

  void load(const I3GeometryConstPtr& geo){
    strs=map<int,STR>();
    doms=map<KEY,DOM>();

    KEY key;
    DOM dom;

    for(I3OMGeoMap::const_iterator it=geo->omgeo.begin(); it!=geo->omgeo.end(); ++it){
      OMKey omkey=it->first;
      key.str=omkey.GetString();
      key.om=omkey.GetOM();
      const I3Position& pos = (it->second).position;
      dom.r[0]=pos.GetX()/I3Units::m;
      dom.r[1]=pos.GetY()/I3Units::m;
      dom.r[2]=pos.GetZ()/I3Units::m+zoff;
#endif
      if(key.str>0 && key.om>=1 && key.om<=60){
	doms[key]=dom;
	STR & s=strs[key.str];
	if(s.sp<0){
	  s.l=key.om; s.h=key.om;
	  for(int i=0; i<3; i++){
	    s.rl[i]=dom.r[i];
	    s.rh[i]=dom.r[i];
	  }
	  s.sp=0;
	}
	else{
	  if(key.om<s.l) s.l=key.om;
	  if(key.om>s.h) s.h=key.om;
	  for(int i=0; i<3; i++){
	    if(dom.r[i]<s.rl[i]) s.rl[i]=dom.r[i];
	    if(dom.r[i]>s.rh[i]) s.rh[i]=dom.r[i];
	  }
	  s.sp++;
	}
      }
#ifndef ICETRAY
      inFile.close();
#endif
    }

    for(map<int,STR>::iterator str=strs.begin(); str!=strs.end(); str++){
      STR & s=str->second;
      s.zl=s.rl[2]; s.zh=s.rh[2]; s.st=(s.h-s.l)/(s.zh-s.zl);
      for(int i=0; i<3; i++){ s.rl[i]-=R; s.rh[i]+=R; }
    }

    for(map<KEY,DOM>::iterator s=doms.begin(); s!=doms.end(); s++){
      STR & str=strs[s->first.str];

      double dz=s->second.r[2]-(str.zl+(str.h-s->first.om)/str.st);
      if(str.zf){
	str.zn=dz;
	str.zx=dz;
	str.zf=false;
      }
      else{
	if(dz<str.zn) str.zn=dz;
	if(dz>str.zx) str.zx=dz;
      }
    }

    for(map<int,STR>::iterator str=strs.begin(); str!=strs.end(); str++){
      STR & s=str->second;
      s.zn-=R; s.zx+=R;
    }

#ifdef XACCL
    bool first=true;
    for(map<int,STR>::const_iterator str=strs.begin(); str!=strs.end(); str++){
      const STR & s=str->second;
      if(first){
	for(int i=0; i<3; i++){
	  rl[i]=s.rl[i];
	  rh[i]=s.rh[i];
	}
	first=false;
      }
      else{
	for(int i=0; i<3; i++){
	  if(rl[i]>s.rl[i]) rl[i]=s.rl[i];
	  if(rh[i]<s.rh[i]) rh[i]=s.rh[i];
	}
      }
    }
    for(int i=0; i<3; i++){ rl[i]-=0.1; rh[i]+=0.1; if(i<2) rst[i]=ncel/(rh[i]-rl[i]); }

    for(int i=0; i<ncel; i++) for(int j=0; j<ncel; j++) cells[i][j].erase(cells[i][j].begin(), cells[i][j].end());
    for(map<int,STR>::const_iterator str=strs.begin(); str!=strs.end(); str++){
      const STR & s=str->second;
      int xl[2], xh[2];
      for(int i=0; i<2; i++){
	xl[i]=(int) floor((s.rl[i]-rl[i])*rst[i]);
	xh[i]=(int) floor((s.rh[i]-rl[i])*rst[i]);
      }
      for(int i=xl[0]; i<=xh[0]; i++) for(int j=xl[1]; j<=xh[1]; j++) cells[i][j].push_back(str->first);
    }

    if(false){
      for(int i=0; i<3; i++) cout<<i<<" "<<rl[i]<<" "<<rh[i]<<endl;
      for(int i=0; i<ncel; i++) for(int j=0; j<ncel; j++){
	cout<<i<<" "<<j;
	for(vector<int>::const_iterator s=cells[i][j].begin(); s!=cells[i][j].end(); s++) cout<<" "<<*s;
	cout<<endl;
      }
    }
#endif

    cerr<<doms.size()<<" DOMs loaded on "<<strs.size()<<" strings"<<endl;
  }

  struct WVD{
    int size;
    double *rx, *wv;

    double ppm, crf[4];

    WVD(){
      ppm=2107.84/(OVR*OVR);
      crf[0]=4.889*0.894; // em shower
      crf[1]=4.076*0.860; // hdr m/GeV
      crf[2]=1.1720;  // bare
      crf[3]=0.0324;  // muon

      size=0;
      ifstream inFile;
      inFile.open((dir+"wv.dat").c_str(), ifstream::in);

      if(!inFile.fail()){
	vector<double> x, y;
	double xa, ya, xo=0, yo=0;
	int num=0;
	bool flag=true;
	while(inFile>>xa>>ya){
	  if(( xa<0 || 1<xa ) || num==0 && xa!=0 || num>0 && ( xa<=xo || ya<=yo )){ flag=false; break; }
	  x.push_back(xa); y.push_back(ya);
	  xo=xa; yo=ya; num++;
	}
	if(xo!=1 || x.size()<2) flag=false;
	inFile.close();
	if(flag){
	  size=x.size();
	  rx = new double[size];
	  wv = new double[size];
	  for(int i=0; i<size; i++){ rx[i]=x[i]; wv[i]=y[i]*1.e-3; }
	  cerr<<"Loaded "<<size<<" wavelenth points from file wv.dat"<<endl;
	}
	else cerr<<"File wv.dat did not contain valid data"<<endl;
      }
      else cerr<<"Could not open file wv.dat"<<endl;
    }

    ~WVD(){
      if(size>0){
	delete rx;
	delete wv;
      }
    }

    double next(){
      double l=0.405;
      double xi=xrnd();
      if(size>1){
	int i=1; for(; i<size-1; i++) if(xi<rx[i]) break;
	l=((xi-rx[i-1])*wv[i]+(rx[i]-xi)*wv[i-1])/(rx[i]-rx[i-1]);
      }
      return l;
    }
  } w;

  struct SRC{
    double r[3], n[3], t;
    int type;

    SRC(){
#ifndef ICETRAY
      load();
#endif
      for(int i=0; i<3; i++){ r[i]=i<2?0:zoff; n[i]=i==0?1:0; }
    }

    void cascade(){
      type=1;
      m.init(w.next());
    }

    void advance(double a){
      t+=a/m.c;
      for(int i=0; i<3; i++) r[i]+=n[i]*a;
    }

    void init_flasher(int str, int om){
      type=0;
      t=0;
      for(map<KEY,DOM>::iterator s=doms.begin(); s!=doms.end(); s++) if(s->first.str==str && s->first.om==om){
	for(int i=0; i<3; i++) r[i]=s->second.r[i]; doms.erase(s); break;
      }
      m.init(0.405);
    }

    void flasher(){
      static const double rms=9.2*sqrt(3)*M_PI/180;
      double xi=xrnd();
      xi=rms*(2*xi-1);
      double nz=sin(xi), np=cos(xi);

      xi=xrnd();
      xi*=2*M_PI;
      n[0]=np*cos(xi);
      n[1]=np*sin(xi);
      n[2]=nz;
    }
  } x;

  struct photon{
    double r[3], n[3], t;
    double xx;

    photon(){
      t=x.t;
      for(int i=0; i<3; i++){
	r[i]=x.r[i]; n[i]=x.n[i];
      }

      if(x.type==1) rotate(m.coschr, m.sinchr);
      xx=1.e-5;
    }

    void rotate(double cs, double si){
      double zi=xrnd();
#ifdef ROTASM
      asm_nq[0]=0, asm_nq[1]=n[0], asm_nq[2]=n[1], asm_nq[3]=n[2];
      asm_zi=zi, asm_zsn=si, asm_zcs=cs;
      rot();
      n[0]=asm_nq[1], n[1]=asm_nq[2]; n[2]=asm_nq[3];
#else
      double u[3];
      {
	zi*=2*M_PI;
	double px=cos(zi), py=sin(zi);

	double p1x, p1y, p1z, p2x, p2y, p2z;
	{
	  double nx=n[0], ny=n[1], nz=n[2];
	  double nnx=nx*nx, nny=ny*ny, nnz=nz*nz;
	  int i=nny>nnx?nnz>nny?2:1:nnz>nnx?2:0;

	  if(i==0){
	    double ry=1/sqrt(nnx+nny), rz=1/sqrt(nnx+nnz);
	    p1x=-ny*ry; p1y=nx*ry; p1z=0;
	    p2x=-nz*rz; p2y=0; p2z=nx*rz;
	  }
	  else if(i==1){
	    double rz=1/sqrt(nny+nnz), rx=1/sqrt(nny+nnx);
	    p1y=-nz*rz; p1z=ny*rz; p1x=0;
	    p2y=-nx*rx; p2z=0; p2x=ny*rx;
	  }
	  else{
	    double rx=1/sqrt(nnz+nnx), ry=1/sqrt(nnz+nny);
	    p1z=-nx*rx; p1x=nz*rx; p1y=0;
	    p2z=-ny*ry; p2x=0; p2y=nz*ry;
	  }
	}
	{
	  double q1x, q1y, q1z, q2x, q2y, q2z;
	  double a1x=p1x-p2x; q1x=a1x; a1x*=a1x;
	  double a1y=p1y-p2y; q1y=a1y; a1y*=a1y;
	  double a1z=p1z-p2z; q1z=a1z; a1z*=a1z;

	  double a2x=p1x+p2x; q2x=a2x; a2x*=a2x;
	  double a2y=p1y+p2y; q2y=a2y; a2y*=a2y;
	  double a2z=p1z+p2z; q2z=a2z; a2z*=a2z;

	  double r1=1/sqrt(a1x+a1y+a1z), r2=1/sqrt(a2x+a2y+a2z);

	  p1x=q1x*r1;
	  p1y=q1y*r1;
	  p1z=q1z*r1;

	  p2x=q2x*r2;
	  p2y=q2y*r2;
	  p2z=q2z*r2;
	}

	u[0]=px*p1x+py*p2x;
	u[1]=px*p1y+py*p2y;
	u[2]=px*p1z+py*p2z;
      }

      {
	double r=0;
	for(int i=0; i<3; i++){
	  double a=n[i]*cs+u[i]*si;
	  n[i]=a; r+=a*a;
	}
	if(abs(1-r)>xx){
	  r=1/sqrt(r);
	  for(int i=0; i<3; i++) n[i]*=r;
	}
      }
#endif
    }

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

      double si=sqrt(1-xi*xi);
      rotate(xi, si);
    }

    void advance(double a){
      t+=a/m.cm;
      for(int i=0; i<3; i++) r[i]+=n[i]*a;
    }

    void propagate(){
      double TOT=-log(xrnd());
      while(true){
	double SCA=-log(xrnd());
	double sca=m.sX(SCA, r[2], n[2]);
	double tot=m.aX(TOT, r[2], n[2]);
	if(sca>tot) sca=tot;

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

	if(om.str>=0){
#ifndef ICETRAY
	  cout<<"HIT "<<om.str<<" "<<om.om<<"  "<<t<<" "<<n[2]<<endl;
#else
	  OMKey key(om.str, om.om);
	  I3MCHit hit(mcpids.second, mcpids.first);
	  hit.SetHitID(hitnum++);
	  hit.SetTime(t*I3Units::ns);
	  // hit.SetWeight(angsens(n[2]));
	  hit.SetWeight(1.0);
	  if(xrnd()<angsens(n[2])) (*mchits)[key].push_back(hit);
#endif
	  break;
	}

	if(abs(TOT)<xx) TOT=0;
	if(TOT==0) break;

	rotate();
      }
    }

    KEY sphere(double &sca){
      KEY imin;
      double hmin=sca;
      double dx[3], dn[3];
      for(int i=0; i<3; i++){
	double start=r[i];
	double end=start+sca*n[i];
	if(n[i]>0){
	  dn[i]=start;
	  dx[i]=end;
	}
	else{
	  dn[i]=end;
	  dx[i]=start;
	}
      }

#ifdef XACCL
      int xl[2], xh[2];
      for(int i=0; i<2; i++){
	xl[i]=max((int) floor((dn[i]-rl[i])*rst[i]), 0);
	xh[i]=min((int) floor((dx[i]-rl[i])*rst[i]), ncel-1);
      }
      if(dn[2]<rh[2] && dx[2]>rl[2]) for(int i=xl[0]; i<=xh[0]; i++) for(int j=xl[1]; j<=xh[1]; j++) for(vector<int>::const_iterator c=cells[i][j].begin(); c!=cells[i][j].end(); c++){
	int sn=*c;
	const STR & s=strs[sn];
#else
      for(map<int,STR>::const_iterator str=strs.begin(); str!=strs.end(); str++){
	int sn=str->first;
	const STR & s=str->second;
#endif
	bool pass=false;
	for(int i=0; i<3; i++) if(s.rh[i]<dn[i] || s.rl[i]>dx[i]){ pass=true; break; }
	if(pass) continue;

	int n1=(int) ceil (s.h-(dx[2]-s.zn-s.zl)*s.st);
	int n2=(int) floor(s.h-(dn[2]-s.zx-s.zl)*s.st);
	if(n1>s.h || n2<s.l) continue;

	for(int om=n1; om<=n2; om++){
	  KEY key;
	  key.str=sn;
	  key.om=om;
	  DOM dom=doms[key];

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

	  if(h>0 && h<sca) if(imin.str<0 || h<hmin){ imin=key; hmin=h; }
	}
      }

      if(imin.str>=0) sca=hmin;
      return imin;
    }
  };

  void flasher(int str, int dom, unsigned long long num){
    x.init_flasher(str, dom);
    for(unsigned long long i=0; i<num; i++){
      x.flasher();
      photon p;
      p.propagate();
    }
  }

  void cascade(double rx, double ry, double rz, double nx, double ny, double nz, double t, double E, int type){
    x.t=t;
    x.r[0]=rx; x.r[1]=ry; x.r[2]=rz;
    x.n[0]=nx; x.n[1]=ny; x.n[2]=nz;
    unsigned long long num=(unsigned long long) roundl(w.ppm*w.crf[type]*E);
    // cerr<<"cascade type "<<type<<", "<<num<<" photons"<<endl;
    for(unsigned long long i=0; i<num; i++){
      x.cascade();
      photon p;
      p.propagate();
    }
  }

  void muon(double rx, double ry, double rz, double nx, double ny, double nz, double t, double E, double dr){
    x.n[0]=nx; x.n[1]=ny; x.n[2]=nz;
    unsigned long long num=(unsigned long long) roundl(max((double) 0, w.ppm*dr*(w.crf[2]+w.crf[3]*log(E))));
    // cerr<<"muon, "<<num<<" photons"<<endl;
    for(unsigned long long i=0; i<num; i++){
      x.t=t;
      x.r[0]=rx; x.r[1]=ry; x.r[2]=rz;
      x.advance(dr*i/num);
      x.cascade();
      photon p;
      p.propagate();
    }
  }
}

#ifndef ICETRAY

using namespace ppc;

void f2k(){
  string in;
  while(getline(cin, in)){
    cout<<in<<endl;
    char name[20];
    int gens, igen;
    double x, y, z, th, ph, l, E, t;
    if(sscanf(in.c_str(), "TR %d %d %s %lf %lf %lf %lf %lf %lf %lf %lf", &gens, &igen, &name, &x, &y, &z, &th, &ph, &l, &E, &t)==11){
      z+=zoff; th=180-th; ph=ph<180?ph+180:ph-180;
      th*=M_PI/180; ph*=M_PI/180;
      double costh=cos(th), sinth=sin(th), cosph=cos(ph), sinph=sin(ph);
      double nx=sinth*cosph, ny=sinth*sinph, nz=costh;
      if(0==strcmp(name, "amu+") || 0==strcmp(name, "amu-") || 0==strcmp(name, "amu")) muon(x, y, z, nx, ny, nz, t, E, l);
      else if(0==strcmp(name, "delta") || 0==strcmp(name, "brems") || 0==strcmp(name, "epair")
	      || 0==strcmp(name, "e+") || 0==strcmp(name, "e-") || 0==strcmp(name, "e")) cascade(x, y, z, nx, ny, nz, t, E, 0);
      else if(0==strcmp(name, "munu") || 0==strcmp(name, "hadr")) cascade(x, y, z, nx, ny, nz, t, E, 1);
    }
  }
}

main(int arg_c, char *arg_a[]){
  if(arg_c<=1) cerr<<"Use: "<<arg_a[0]<<" [str] [om] [num] [seed]"<<endl;
  else if(arg_a[1][0]=='-'){
    m.init(-atof(arg_a[1]));
    for(int i=0; i<m.size; i++) cout<<m.hmin+m.dh*(0.5+i)<<" "<<m.abs[i]<<" "<<m.sca[i]<<endl;
  }
  else if(arg_c<=2){
    unsigned int seed=atoi(arg_a[1]);
    srand(1+seed);
    f2k();
  }
  else{
    int str=0, dom=0;
    unsigned int seed=0;
    unsigned long long num=1000000ULL;

    if(arg_c>1) str=atoi(arg_a[1]);
    if(arg_c>2) dom=atoi(arg_a[2]);
    if(arg_c>3) num=(unsigned long long) atof(arg_a[3]);
    if(arg_c>4) seed=atoi(arg_a[4]);

    srand(1+seed);
    flasher(str, dom, num);
  }
}

#endif