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

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

namespace ppc{
  using namespace std;
#include "ice.cxx"

  double xrnd(){
    double rnd=0;
    while(rnd==0){
      rnd=rand();
      rnd/=RAND_MAX;
    }
    return 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;
    double zn, zx, sp;
    double rl[3], rh[3];

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

  map<int,STR> strs;

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

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

  void load(){
    ifstream inFile;
    inFile.open("geo-f2k", 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) 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++;
	}
      }
      inFile.close();
    }

    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.rl[2]+(s->first.om-str.h)*(str.rh[2]-str.rl[2])/(str.l-str.h));
      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;
      }
    }

    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("wv.dat", 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(){
      load();
      for(int i=0; i<3; i++){ r[i]=i<2?0:-1950; 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 u[3];
      {
	double p1[3], p2[3];
	{
	  int i0=0;
	  double n0=n[0];
	  for(int i=1; i<3; i++) if(abs(n[i])>abs(n0)){ i0=i; n0=n[i]; }
	  int i1=(i0+1)%3, i2=(i0+2)%3;
	  double n1=n[i1], n2=n[i2], nq=n0*n0;
	  double r1=sqrt(nq+n1*n1), r2=sqrt(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;
	}
	{
	  double q1[3], q2[3], r1=0, r2=0;
	  for(int i=0; i<3; i++){
	    double a1=p1[i]-p2[i]; q1[i]=a1; r1+=a1*a1;
	    double a2=p1[i]+p2[i]; q2[i]=a2; r2+=a2*a2;
	  }

	  r1=sqrt(r1); r2=sqrt(r2);

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

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

      {
	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=sqrt(r);
	  for(int i=0; i<3; i++) n[i]/=r;
	}
      }
    }

    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 tot=m.aX(TOT, r[2], n[2]);
	double SCA=-log(xrnd());
	double sca=m.sX(SCA, 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){
	  cout<<"HIT "<<om.str<<" "<<om.om<<"  "<<t<<" "<<n[2]<<endl;
	  break;
	}

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

	rotate();
      }
    }

    KEY sphere(double &sca){
      KEY imin;
      double hmin=sca;
      for(map<int,STR>::const_iterator str=strs.begin(); str!=strs.end(); str++){
	bool pass=false;
	const STR & s=str->second;
	for(int i=0; i<3; i++){
	  if(n[i]>0){
	    if(s.rh[i]+R<r[i] || s.rl[i]-R>r[i]+sca*n[i]){
	      pass=true;
	      break;
	    }
	  }
	  else{
	    if(s.rl[i]-R>r[i] || s.rh[i]+R<r[i]+sca*n[i]){
	      pass=true;
	      break;
	    }
	  }
	}
	if(pass) continue;

	double min, max;
	double st=(s.h-s.l)/(s.rh[2]-s.rl[2]);
	if(n[2]>0){
	  min=s.h-(r[2]+sca*n[2]-(s.zn-R)-s.rl[2])*st;
	  max=s.h-(r[2]-(s.zx+R)-s.rl[2])*st;
	}
	else{
	  min=s.h-(r[2]-(s.zn-R)-s.rl[2])*st;
	  max=s.h-(r[2]+sca*n[2]-(s.zx+R)-s.rl[2])*st;
	}
	int n1=(int)floor(min);
	int n2=(int)ceil(max);
	if(n1>s.h | n2<s.l) continue;

	for(int om=n1; om<=n2; om++){
	  KEY key;
	  key.str=str->first;
	  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-=1950; 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