#if defined(MKOW) || defined(XLIB)
float xrnd(){
  unsigned int rnd;
  do rnd=rand_r(&sv); while(rnd==0);
  const float a=1.0f/(1ll+RAND_MAX);
  return a*rnd;
}

float grnd(){  // gaussian distribution
  return sqrtf(-2*logf(xrnd()))*sinf(2*FPI*xrnd());
}
#endif

static const float ppm=2450.08;     // photons per meter
static const float rho=0.9216;      // density of ice [mwe]
static const float m0=0.105658389;  // muon rest mass [GeV]

photon p;

float yield(float E, int type){
  p.l=0;
#ifdef ANGW
  p.f=0;
#endif
#ifdef LONG
  float logE=logf(max(m0, E));
  const float Lrad=0.358/rho;
  if(type>0){  // hardonic shower
    p.a=1.49f+0.359f*logE, p.b=Lrad/0.772f;
  }
  else{  // em shower
    p.a=2.03f+0.604f*logE, p.b=Lrad/0.633f;
  }
#endif

  float nph;
#ifdef MKOW
  const float em=5.21*0.924f/rho;

  float f=1.0f;
  if(type>0){
    const float E0=0.399;
    const float m=0.130;
    const float f0=0.467;
    const float rms0=0.379;
    const float gamma=1.160;

    float e=max(10.0f, E);
    float F=1-powf(e/E0, -m)*(1-f0);
    float dF=F*rms0*powf(log10f(e), -gamma);
    do f=F+dF*grnd(); while(f<0 || 1<f);
  }
  nph=f*em;
#else
  const float em=4.889*0.894f/rho;  // em shower
  const float hr=4.076*0.860f/rho;  // hdr [m/GeV]
  nph=type>0?hr:em;
#endif
  return d.eff*ppm*nph*E;
}

float yield(float E, float dr){     // bare muon
  float logE=logf(max(m0, E));
#ifdef MKOW
  float extr=1+max(0.0f, 0.1880f+0.0206f*logE);
#else
  float extr=1+max(0.0f, 0.1720f+0.0324f*logE);
#endif
  float nph=dr>0?dr*extr:0;
  p.l=dr;
#ifdef ANGW
  p.f=1/extr;
#endif
#ifdef LONG
  p.a=0, p.b=0;
#endif
  return d.eff*ppm*nph;
}

void output();

#ifdef XLIB
struct mcid:pair<int,unsigned long long>{
  int frame;
};

struct ihit{
  ikey omkey;
  mcid track;
  float time;
} tmph;

deque<mcid> flnz;
vector<ihit> hitz;
#else
deque<string> flnz;
#endif

unsigned int flnb=0, flne=0;
#ifdef XCPU
unsigned int & flnd = flne;
#else
unsigned int flnd=0;
#endif

#ifndef RAND
int sign(float x){ return x<0?-1:x>0?1:0; }
#endif

int hcmp(const void *a, const void *b){
  hit & ah = * (hit *) a;
  hit & bh = * (hit *) b;
#ifdef RAND
  return (int) ( ah.n - bh.n );
#else
  return ah.n!=bh.n ? (int)(ah.n-bh.n) : ah.i!=bh.i ? (int)(ah.i-bh.i) : ah.t!=bh.t ? sign(ah.t-bh.t) : sign(ah.z-bh.z);
#endif
}

void print(){
#ifdef RAND
  if(flnb<flnd)
#endif
    qsort(q.hits, d.hidx, sizeof(hit), hcmp);

  for(unsigned int i=0; i<d.hidx; i++){
    hit & h = q.hits[i];

    for(; (int) ( flnb - h.n ) < 0; flnb++){
#ifdef XLIB
      tmph.track=flnz.front();
#else
      cout<<flnz.front()<<endl;
#endif
      flnz.pop_front();
    }

    name n=q.names[h.i];
#ifdef XLIB
    tmph.omkey.str=(int) n.str;
    tmph.omkey.dom=(int) n.dom;
    tmph.time=h.t;
    hitz.push_back(tmph);
#else
    printf("HIT %d %d %f %f\n", (int) n.str, (int) n.dom, h.t, h.z);
#endif
  }

  for(; (int) ( flnb - flnd ) < 0; flnb++){
#ifdef XLIB
    tmph.track=flnz.front();
#else
    cout<<flnz.front()<<endl;
#endif
    flnz.pop_front();
  }
}

unsigned long long f2ki, f2kn;

template <class T> void addp(float rx, float ry, float rz, float t, float E, T xt){
  p.r.w=t; p.r.x=rx; p.r.y=ry; p.r.z=rz;
  unsigned long long num=llroundf(yield(E, xt));
  for(f2kn+=num; f2ki<f2kn; f2ki+=ovr){
    q.pz[pn++]=p; if(pn==pmxo) output();
  }
}

#ifdef XLIB
template void addp(float, float, float, float, float, int);
template void addp(float, float, float, float, float, float);
#endif

void eini(){
  f2ki=0;
  f2kn=0;
}

void finc(){
  flne++;
  if((int) ( flne - flnb ) < 0){ cerr<<"Error: too many event segments"<<endl; exit(3); }
}

void eout(){
  output();
#ifndef XCPU
  output();
#endif
}

#ifdef XLIB
void efin(){
  hitz.erase(hitz.begin(), hitz.end());
}

void sett(float nx, float ny, float nz, pair<int,unsigned long long> id, int frame){
  mcid ID; ID.first=id.first; ID.second=id.second; ID.frame=frame;
  flnz.push_back(ID); finc();
  p.q=flne; p.n.x=nx; p.n.y=ny; p.n.z=nz;
}
#else

void f2k(){
  ini(0);

  string in;
  while(getline(cin, in)){
    flnz.push_back(in); finc();
    char name[32];
    int gens, igen;
    float x, y, z, th, ph, l, E, t;
    const char * str = "TR %d %d %31s %f %f %f %f %f %f %f %f";

    if(sscanf(in.c_str(), "%31s", name)==1){ if(strcmp(name, "EM")==0) eini(); }
    if(sscanf(in.c_str(), str, &gens, &igen, name, &x, &y, &z, &th, &ph, &l, &E, &t)==11){
      th=180-th; ph=ph<180?ph+180:ph-180;
      th*=FPI/180; ph*=FPI/180;
      float costh=cosf(th), sinth=sinf(th), cosph=cosf(ph), sinph=sinf(ph);
      p.q=flne; p.n.x=sinth*cosph; p.n.y=sinth*sinph; p.n.z=costh;
      if(0==strcmp(name, "amu+") || 0==strcmp(name, "amu-") || 0==strcmp(name, "amu")) addp(x, y, z, 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")) addp(x, y, z, t, E, 0);
      else if(0==strcmp(name, "munu") || 0==strcmp(name, "hadr")) addp(x, y, z, t, E, 1);
    }
  }
  eout();

  fin();
}
#endif