static const unsigned int ovr=OVR*OVR*OVER;
static const float ppm=2107.84*EFF*MAS;
static const float crf[4]={ ppm*4.889*0.894, // em shower
			    ppm*4.076*0.860, // hdr m/GeV
			    ppm*1.1720,      // bare
			    ppm*0.0324 };    // muon
#ifdef MKOW
static const float mkw[6]={ ppm*5.21,        // em shower
			    0.399,           // E0
			    0.130,           // m
			    0.467,           // f0
			    0.379,           // RMS0
			    1.160 };         // gamma

float xrnd(){
  float rnd;
  do rnd=rand(); while(rnd==0);
  const float a=1.0f/(1ll+RAND_MAX);
  return a*rnd;
}

float grnd(){
  return sqrtf(-2*logf(xrnd()))*sinf(2*fpi*xrnd());
}
#endif

float yield(float E, int type){
#ifdef MKOW
  float f=1.0f;
  if(type>0){
    float e=max(10.0f, E);
    float F=1+powf(e/mkw[1], -mkw[2])*(mkw[3]-1);
    float dF=mkw[4]*pow(log10f(e), -mkw[5]);
    do f=F+dF*grnd(); while(f<0 || 1<f);
  }
  return f*mkw[0]*E;
#else
  return crf[type]*E;
#endif
}

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;

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

void print(){
  if(flnb<flne) qsort(d.hits, e.hidx, sizeof(hit), hcmp);

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

    for(; flnb<h.n; 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(; flnb<flne; flnb++){
#ifdef XLIB
    tmph.track=flnz.front();
#else
    cout<<flnz.front()<<endl;
#endif
    flnz.pop_front();
  }
  p.q=0; flnb=0; flne=0;
}

unsigned long long f2ki, f2kn;

void cascade(float rx, float ry, float rz, float t, float E, int type){
  p.t=t; p.r[0]=rx; p.r[1]=ry; p.r[2]=rz;
  unsigned long long num=llroundf(yield(E, type));
  f2kn+=num;
  for(; f2ki<f2kn; f2ki+=ovr){
    pz[pn++]=p; if(pn==pmax/OVER) output();
  }
}

void muon(float rx, float ry, float rz, float t, float E, float dr){
  unsigned long long num=llroundf(max(0.0f, dr*(crf[2]+crf[3]*logf(E))));
  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;
    pz[pn++]=p; if(pn==pmax/OVER) output();
  }
}

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

#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); flne++;
  p.q=flne; p.n[0]=nx; p.n[1]=ny; p.n[2]=nz;
}

void igeo(){
  int type=e.type;
  fin(); m.geo();
  e=d; ini(type);
}
#else

void f2k(){
  ini(0);

  string in;
  while(getline(cin, in)){
    flnz.push_back(in); flne++;
    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);
    }
  }
  output();

  fin();
}
#endif