#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
#define round roundf
#endif

#define XACCL

namespace ppc{
  using namespace std;

#ifdef SINGLE
  struct setprc{
    short store;
    setprc(){
      short byte=127, end;
      __asm__("fstcw %0\nfldcw %1\nfstcw %2"::"m" (store), "m" (byte), "m" (end));
      // cerr<<store<<" "<<byte<<" "<<end<<endl;
    }
    ~setprc(){
      short byte, end;
      __asm__("fstcw %0\nfldcw %1\nfstcw %2"::"m" (byte), "m" (store), "m" (end));
      // cerr<<byte<<" "<<store<<" "<<end<<endl;
    }
  } f;
#endif

  /*
  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(){
    double rnd;
    do rnd=rand(); while(rnd==0);
    const double a=1./(1ll+RAND_MAX);
    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 OVR=5;
  static const double R=OVR*0.16510;
  static const double zoff=-1948.07;

  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 rho, ppm, crf[4];

    WVD(){
      rho=0.9216;
      ppm=2450.08/(OVR*OVR);
      crf[0]=4.889*0.894/rho; // em shower
      crf[1]=4.076*0.860/rho; // 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 t, r[3], n[3];
    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{
#ifdef ROTASM
    volatile double t __attribute__ (( aligned(16) )), r[3] __attribute__ ((packed));
    volatile double q __attribute__ (( aligned(16) )), n[3] __attribute__ ((packed));
#else
    double t, r[3], n[3];
#endif
    double xx;

    photon(){
      t=x.t;
#ifdef ROTASM
      q=0;
#endif
      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(); zi*=2*M_PI;
      double px=cos(zi), py=sin(zi);
#ifdef ROTASM
      volatile __attribute__ (( aligned(16) )) float xcos[2]={px, py}, zcos[2]={cs, si};

      __asm__("	movaps %0,%%xmm0	\n"
	      "	movlps %1,%%xmm1	\n"
	      "	movlps %2,%%xmm5	\n"
	      ".intel_syntax noprefix	\n"
	      "	movaps xmm7,xmm0	\n"
	      "	mulps xmm7,xmm7		\n"
	      "	movaps xmm2,xmm7	\n"
	      "	shufps xmm2,xmm2,0xe5	\n"
	      "	movaps xmm3,xmm7	\n"
	      "	shufps xmm3,xmm3,0xe6	\n"
	      "	movaps xmm4,xmm7	\n"
	      "	shufps xmm4,xmm4,0xe7	\n"
	      "	comiss xmm2,xmm3	\n"
	      "	jc 1f			\n"
	      "	comiss xmm2,xmm4	\n"
	      "	jc 5f			\n"
	      "	jmp 3f			\n"
	      "1:			\n"
	      "	comiss xmm3,xmm4	\n"
	      "	jc 5f			\n"
	      "	jmp 4f			\n"
	      "3:			\n"
	      "	movaps xmm7,xmm2	\n"
	      "	shufps xmm7,xmm7,0x4e	\n"
	      "	addps xmm7,xmm2		\n"
	      "	rsqrtps xmm7,xmm7	\n"
	      "	movaps xmm3,xmm0	\n"
	      "	shufps xmm3,xmm3,0x10	\n"
	      "	movaps xmm4,xmm0	\n"
	      "	shufps xmm4,xmm4,0x08	\n"
	      "	subps xmm3,xmm4		\n"
	      "	movaps xmm6,xmm7	\n"
	      "	shufps xmm6,xmm6,0x00	\n"
	      "	mulps xmm6,xmm3		\n"
	      "	movaps xmm3,xmm0	\n"
	      "	shufps xmm3,xmm3,0x40	\n"
	      "	movaps xmm4,xmm0	\n"
	      "	shufps xmm4,xmm4,0x0c	\n"
	      "	subps xmm3,xmm4		\n"
	      "	shufps xmm7,xmm7,0x55	\n"
	      "	mulps xmm7,xmm3		\n"
	      "	jmp 2f			\n"
	      "4:			\n"
	      "	movaps xmm7,xmm3	\n"
	      "	shufps xmm7,xmm7,0x1b	\n"
	      "	addps xmm7,xmm3		\n"
	      "	rsqrtps xmm7,xmm7	\n"
	      "	movaps xmm4,xmm0	\n"
	      "	shufps xmm4,xmm4,0x80	\n"
	      "	movaps xmm2,xmm0	\n"
	      "	shufps xmm2,xmm2,0x30	\n"
	      "	subps xmm4,xmm2		\n"
	      "	movaps xmm6,xmm7	\n"
	      "	shufps xmm6,xmm6,0x00	\n"
	      "	mulps xmm6,xmm4		\n"
	      "	movaps xmm4,xmm0	\n"
	      "	shufps xmm4,xmm4,0x08	\n"
	      "	movaps xmm2,xmm0	\n"
	      "	shufps xmm2,xmm2,0x10	\n"
	      "	subps xmm4,xmm2		\n"
	      "	shufps xmm7,xmm7,0x55	\n"
	      "	mulps xmm7,xmm4		\n"
	      "	jmp 2f			\n"
	      "5:			\n"
	      "	movaps xmm7,xmm4	\n"
	      "	shufps xmm7,xmm7,0x4e	\n"
	      "	addps xmm7,xmm4		\n"
	      "	rsqrtps xmm7,xmm7	\n"
	      "	movaps xmm2,xmm0	\n"
	      "	shufps xmm2,xmm2,0x0c	\n"
	      "	movaps xmm3,xmm0	\n"
	      "	shufps xmm3,xmm3,0x40	\n"
	      "	subps xmm2,xmm3		\n"
	      "	movaps xmm6,xmm7	\n"
	      "	shufps xmm6,xmm6,0x55	\n"
	      "	mulps xmm6,xmm2		\n"
	      "	movaps xmm2,xmm0	\n"
	      "	shufps xmm2,xmm2,0x30	\n"
	      "	movaps xmm3,xmm0	\n"
	      "	shufps xmm3,xmm3,0x80	\n"
	      "	subps xmm2,xmm3		\n"
	      "	shufps xmm7,xmm7,0x00	\n"
	      "	mulps xmm7,xmm2		\n"
	      "2:			\n"
	      "	movaps xmm3,xmm6	\n"
	      "	subps xmm3,xmm7		\n"
	      "	addps xmm7,xmm6		\n"
	      "	movaps xmm4,xmm3	\n"
	      "	mulps xmm4,xmm4		\n"
	      "	movaps xmm6,xmm7	\n"
	      "	mulps xmm6,xmm6		\n"
	      "	movaps xmm2,xmm4	\n"
	      "	shufps xmm2,xmm6,0x11	\n"
	      "	shufps xmm4,xmm6,0xee	\n"
	      "	addps xmm4,xmm2		\n"
	      "	movaps xmm2,xmm4	\n"
	      "	shufps xmm2,xmm2,0xb1	\n"
	      "	addps xmm4,xmm2		\n"
	      "	rsqrtps xmm4,xmm4	\n"
	      "	movaps xmm6,xmm4	\n"
	      "	shufps xmm4,xmm4,0x00	\n"
	      "	shufps xmm6,xmm6,0xaa	\n"
	      "	mulps xmm3,xmm4		\n"
	      "	mulps xmm7,xmm6		\n"
	      ""
	      "	movaps xmm6,xmm1	\n"
	      "	shufps xmm6,xmm6,0x00	\n"
	      "	shufps xmm1,xmm1,0x55	\n"
	      "	mulps xmm6,xmm3		\n"
	      "	mulps xmm1,xmm7		\n"
	      "	addps xmm1,xmm6		\n"
	      ""
	      "	movaps xmm3,xmm5	\n"
	      "	shufps xmm3,xmm3,0x00	\n"
	      "	shufps xmm5,xmm5,0x55	\n"
	      "	mulps xmm1,xmm5		\n"
	      "	mulps xmm0,xmm3		\n"
	      "	addps xmm0,xmm1		\n"
	      "	movaps xmm7,xmm0	\n"
	      "	mulps xmm7,xmm7		\n"
	      "	shufps xmm7,xmm7,0x39	\n"
	      "	movss xmm5,xmm7		\n"
	      "	shufps xmm7,xmm7,0x39	\n"
	      "	movss xmm3,xmm7		\n"
	      "	shufps xmm7,xmm7,0x39	\n"
	      "	movss xmm1,xmm7		\n"
	      "	addss xmm5,xmm3		\n"
	      "	addss xmm5,xmm1		\n"
	      "	rsqrtss xmm5,xmm5	\n"
	      "	shufps xmm5,xmm5,0x00	\n"
	      "	mulps xmm0,xmm5		\n"
	      ".att_syntax prefix	\n"
	      "	movaps %%xmm0,%0	\n"
	      ""::"m" (q), "m" (xcos[0]), "m" (zcos[0]):);

      return;
#else
      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=1/sqrt(nq+n1*n1), r2=1/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=1/sqrt(r1); r2=1/sqrt(r2);

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

	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=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 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){
#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=llround(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=llround(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[32];
    int gens, igen;
    double x, y, z, th, ph, l, E, t;
#ifdef SINGLE
    const char * str = "TR %d %d %31s %f %f %f %f %f %f %f %f";
#else
    const char * str = "TR %d %d %31s %lf %lf %lf %lf %lf %lf %lf %lf";
#endif
    if(sscanf(in.c_str(), str, &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]));
    cerr<<"For l="<<m.wv<<"nm: np="<<m.np<<" ng="<<m.ng<<endl;
    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