static const unsigned int ovr=OVR*OVR*OVER;
static const float ppm=2450.08*EFF;
static const float m0=0.105658389;  // muon rest mass [GeV]

float maxe=1.0f;

struct light{
#if defined(MKOW) || defined(LONG) || defined(ANGW)
#define XRND
  float xrnd(){
    float rnd;
    do rnd=rand(); while(rnd==0);
    const float a=1.0f/(1ll+RAND_MAX);
    return a*rnd;
  }
#endif

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

#ifdef LONG
  float mrnd(float k){  // gamma distribution
    float x;
    if(k<1){  // Weibull algorithm
      float c=1/k;
      float d=(1-k)*powf(k, k/(1-k));
      float z, e;
      do{
	z=-logf(xrnd());
	e=-logf(xrnd());
	x=powf(z, c);
      } while(z+e<d+x);
    }
    else{  // Cheng's algorithm
      float b=k-logf(4);
      float l=sqrtf(2*k-1);
      float c=1+logf(4.5);
      float u, v, y, z, r;
      do{
	u=xrnd(); v=xrnd();
	y=logf(v/(1-v))/l;
	x=k*expf(y);
	z=u*v*v;
	r=b+(k+l)*y-x;
      } while(r<4.5*z-c && r<logf(z));
    }
    return x;
  }

  float a, b;  // parameters of longitudinal shower profile

  float dr(){
    const float Lrad=0.358/0.917;
    return Lrad*mrnd(a)/b;
  }
#endif

#ifdef ANGW
  float r;  // fraction of light from muon alone (without cascades)

  void rset(unsigned int & x){
    if(r==0) x|=0x80000000;
    else{ if(r<xrnd()) x|=0x80000000; else x&=0x7FFFFFFF; }
  }
#endif

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

    float nph;
#ifdef MKOW
    const float em=5.21;

    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=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;  // em shower
    const float hr=4.076*0.860f;  // hdr [m/GeV]
    nph=type>0?hr:em;
#endif
    return maxe*ppm*nph*E;
  }

  float yield(float E, float l){  // bare muon
    float logE=logf(max(m0, E));
    float extr=max(0.0f, 0.1720f+0.0324f*logE);
    float nph=l>0?l*(1+extr):0;
#ifdef ANGW
    r=1/(1+extr);
#endif
    return maxe*ppm*nph;
  }
} l;

photon p;
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

bool rcmp(unsigned int a, unsigned int b){
  return ! ( a - b < 0x40000000 );
}

void rinc(unsigned int & x){
  if( ++x == 0x80000000 ) x=0;
}

int hcmp(const void *a, const void *b){
  return rcmp(((hit *) b)->n, ((hit *) a)->n);
}

void print(){
  if(flnb<flnd) qsort(d.hits, eh->hidx, sizeof(hit), hcmp);

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

    for(; rcmp(flnb, h.n); rinc(flnb)){
#ifdef XLIB
      tmph.track=flnz.front();
#else
      cout<<flnz.front()<<endl;
#endif
      flnz.pop_front();
    }

    name n=d.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(; rcmp(flnb, flnd); rinc(flnb)){
#ifdef XLIB
    tmph.track=flnz.front();
#else
    cout<<flnz.front()<<endl;
#endif
    flnz.pop_front();
  }
}

unsigned long long f2ki, f2kn;

void cascade(float rx, float ry, float rz, float t, float E, int type){
#ifndef LONG
  p.t=t; p.r[0]=rx; p.r[1]=ry; p.r[2]=rz;
#endif
  unsigned long long num=llroundf(l.yield(E, type));
  f2kn+=num;
  for(; f2ki<f2kn; f2ki+=ovr){
#ifdef LONG
    p.t=t; p.r[0]=rx; p.r[1]=ry; p.r[2]=rz;
    float a=l.dr(); p.t+=a*d.ocv;
    for(int i=0; i<3; i++) p.r[i]+=p.n[i]*a;
#endif
#ifdef ANGW
    l.rset(p.q);
#endif
    d.pz[pn++]=p; if(pn==pmxo) output();
  }
}

void muon(float rx, float ry, float rz, float t, float E, float dr){
  unsigned long long num=llroundf(l.yield(E, dr));
  f2kn+=num;
  for(unsigned long long i=0; f2ki<f2kn; f2ki+=ovr, i+=ovr){
    p.t=t; p.r[0]=rx; p.r[1]=ry; p.r[2]=rz;
    float a=(dr*i)/num; p.t+=a*d.ocv;
    for(int i=0; i<3; i++) p.r[i]+=p.n[i]*a;
#ifdef ANGW
    l.rset(p.q);
#endif
    d.pz[pn++]=p; if(pn==pmxo) output();
  }
}

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

void finc(){
  rinc(flne);
  if(rcmp(flne, flnb)){ 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[0]=nx; p.n[1]=ny; p.n[2]=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[0]=sinth*cosph; p.n[1]=sinth*sinph; p.n[2]=costh;
      if(0==strcmp(name, "amu+") || 0==strcmp(name, "amu-") || 0==strcmp(name, "amu")) muon(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")) cascade(x, y, z, t, E, 0);
      else if(0==strcmp(name, "munu") || 0==strcmp(name, "hadr")) cascade(x, y, z, t, E, 1);
    }
  }
  eout();

  fin();
}
#endif