#include <stdio.h>

#include "ini.h"
#include "pro.cu"

extern dats d;
dats e=d; // copy of "d" for device
dats *ed; // pointer to structure on device

static const unsigned long long pmax=(unsigned long long)NPHO*NBLK*NTHR;
unsigned long long pn;
photon * pz;

unsigned int mode=0, fsim=0;
float deviceTime=0;

void checkError(cudaError result){
  if(result!=cudaSuccess){
    fprintf(stderr, "Error: %s\n", cudaGetErrorString(result));
    exit(1);
  }
}

void ini(int type){
  {
    pn=0; fsim=type; e.type=type; e.hidx=0;
    for(int m=0; m<3; m++) e.r[m]=0;
  }

  {
    unsigned int size=d.rsize;
    if(size<NBLK*NTHR) fprintf(stderr, "Error: not enough multipliers: only have %d (need %d)!\n", size, NBLK*NTHR);
  }

  {
    unsigned int size=d.rsize*sizeof(unsigned int);
    checkError(cudaMalloc((void**) &e.rm, size));
    checkError(cudaMemcpy(e.rm, d.rm, size, cudaMemcpyHostToDevice));
  }

  {
    unsigned int size=d.rsize*sizeof(unsigned long long);
    checkError(cudaMalloc((void**) &e.rs, size));
    checkError(cudaMemcpy(e.rs, d.rs, size, cudaMemcpyHostToDevice));
  }

  for(int i=0; i<WNUM; i++){
    unsigned int size=sizeof(ices);
    checkError(cudaMalloc((void**) &e.w[i], size));
    checkError(cudaMemcpy(e.w[i], d.w[i], size, cudaMemcpyHostToDevice));
  }

  {
    unsigned int size=d.hnum*sizeof(hit);
    checkError(cudaMalloc((void**) &e.hits, size));
  }

  {
    pz = new photon[pmax];
    unsigned int size=pmax*sizeof(photon);
    checkError(cudaMalloc((void**) &e.pz, size));
  }

  {
    unsigned int size=d.gsize;
    checkError(cudaMemcpyToSymbol(oms, d.oms, size*sizeof(DOM)));
  }

  {
    checkError(cudaMalloc((void**) &ed, sizeof(dats)));
    checkError(cudaMemcpy(ed, &e, sizeof(dats), cudaMemcpyHostToDevice));
  }
}

void fin(){
  checkError(cudaFree(e.rm));
  checkError(cudaFree(e.rs));
  for(int i=0; i<WNUM; i++) checkError(cudaFree(e.w[i]));
  checkError(cudaFree(e.hits));
  checkError(cudaFree(e.pz));
  checkError(cudaFree(ed));
  delete pz;
}

void listDevices(){
  int deviceCount, driver, runtime;
  cudaGetDeviceCount(&deviceCount);
  cudaDriverGetVersion(&driver);
  cudaRuntimeGetVersion(&runtime);
  fprintf(stderr, "\nFound %d devices, driver %d, runtime %d\n", deviceCount, driver, runtime);
  for(int device=0; device<deviceCount; ++device){
    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, device);
    fprintf(stderr, "%d(%d.%d): %s %g GHz G(%lu) S(%lu) C(%lu) R(%d) W(%d)\n"
	    "\tl%d o%d c%d h%d i%d m%d a%lu M(%lu) T(%d: %d,%d,%d) G(%d,%d,%d)\n",
	    device, prop.major, prop.minor, prop.name, prop.clockRate/1.e6,
	    (unsigned long)prop.totalGlobalMem, (unsigned long)prop.sharedMemPerBlock,
	    (unsigned long)prop.totalConstMem, prop.regsPerBlock, prop.warpSize,
	    prop.kernelExecTimeoutEnabled, prop.deviceOverlap, prop.computeMode,
	    prop.canMapHostMemory, prop.integrated, prop.multiProcessorCount,
	    (unsigned long)prop.textureAlignment,
	    (unsigned long)prop.memPitch, prop.maxThreadsPerBlock,
	    prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2],
	    prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
  }
}

void kernel(unsigned long long num){
  {
    float dt;
    cudaEvent_t evt1, evt2;

    checkError(cudaEventCreateWithFlags(&evt1, cudaEventBlockingSync));
    checkError(cudaEventCreateWithFlags(&evt2, cudaEventBlockingSync));

    checkError(cudaEventRecord(evt1, NULL));
    checkError(cudaEventSynchronize(evt1));

    propagate<<< fsim>0&&mode>0?WNUM:NBLK, NTHR >>>(ed, num);
    checkError(cudaGetLastError());

    checkError(cudaEventRecord(evt2, NULL));
    checkError(cudaEventSynchronize(evt2));

    checkError(cudaEventElapsedTime(&dt, evt1, evt2));
    deviceTime+=dt;

    checkError(cudaEventDestroy(evt1));
    checkError(cudaEventDestroy(evt2));

    checkError(cudaThreadSynchronize());
  }

  {
    checkError(cudaMemcpy(&e, ed, sizeof(int), cudaMemcpyDeviceToHost));
    unsigned int size=e.hidx*sizeof(hit);
    checkError(cudaMemcpy(d.hits, e.hits, size, cudaMemcpyDeviceToHost));
  }

  if(fsim==0||mode==0) fprintf(stderr, "photons: %llu  hits: %u\n", num, e.hidx);

  for(int i=0; i<e.hidx; i++){
    name n=d.names[d.hits[i].i];
    printf("HIT %d %d %f\n", (int) n.str, (int) n.dom, d.hits[i].t);
  }

  if(e.hidx>=e.hnum) fprintf(stderr, "Error: data buffer overflow occurred!\n");
}

void flasher(int str, int dom, unsigned long long num){
  ini(0);

  for(long long i=num; i>0; i-=pmax){
    e.hidx=0;
    for(int n=0; n<d.gsize; n++) if(d.names[n].str==str && d.names[n].dom==dom){
      for(int m=0; m<3; m++) e.r[m]=d.oms[n].r[m]; break;
    }
    checkError(cudaMemcpy(ed, &e, 4*sizeof(int), cudaMemcpyHostToDevice));

    if(mode==0) kernel(min(i, pmax));
    else{ kernel(min(i, num)); break; }
  }

  fin();
}

void output(){
  {
    e.hidx=0;
    checkError(cudaMemcpy(ed, &e, sizeof(int), cudaMemcpyHostToDevice));
  }

  {
    unsigned int size=pn*sizeof(photon);
    checkError(cudaMemcpy(e.pz, pz, size, cudaMemcpyHostToDevice));
  }

  kernel(pn);

  pn=0;
}

static const float ppm=2107.84/(OVR*OVR);
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

photon p;

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(crf[type]*E);
  for(unsigned long long i=0; i<num; i++){
    pz[pn++]=p; if(pn==pmax) output();
  }
}

void muon(float rx, float ry, float rz, float t, float E, float dr){
  unsigned long long num=llroundf(max((float) 0, dr*(crf[2]+crf[3]*logf(E))));
  for(unsigned long long i=0; i<num; i++){
    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) output();
  }
}

void f2k(){
  ini(1);

  char line[512];
  while(scanf("%511[^\n]%*1[\n]", line)!=-1){
    char name[32];
    if(mode>0){
      sscanf(line, "%31s", name);
      if(0==strcmp(name, "EE")) output();
    }
    printf("%s\n", line);
    int gens, igen;
    float x, y, z, th, ph, l, E, t;
    char * str = "TR %d %d %31s %f %f %f %f %f %f %f %f";
    if(sscanf(line, str, &gens, &igen, &name, &x, &y, &z, &th, &ph, &l, &E, &t)==11){
      th=180-th; ph=ph<180?ph+180:ph-180;
      th*=M_PI/180; ph*=M_PI/180;
      float costh=cosf(th), sinth=sinf(th), cosph=cosf(ph), sinph=sinf(ph);
      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();
}

int main(int arg_c, char *arg_a[]){
  cudaSetDeviceFlags(cudaDeviceBlockingSync);
  if(arg_c<=1){
    listDevices();
    fprintf(stderr, "\nUse: %s [device] [mode] (f2k muons)\n"
	    "     %s [str] [om] [num] [device] [mode] (flasher)\n", arg_a[0], arg_a[0]);
  }
  else if(arg_c<=3){
    int device=0;
    if(arg_c>1) device=atoi(arg_a[1]);
    if(arg_c>2) mode=atoi(arg_a[2]);
    checkError(cudaSetDevice(device));
    fprintf(stderr, "Processing f2k muons from stdin on device %d\n", device);
    f2k();
  }
  else{
    int str=0, dom=0, device=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) device=atoi(arg_a[4]);
    if(arg_c>5) mode=atoi(arg_a[5]);
    checkError(cudaSetDevice(device));
    fprintf(stderr, "Running flasher simulation on device %d\n", device);
    flasher(str, dom, num);
  }

  fprintf(stderr, "\nDevice time: %2.1f [ms]\n", deviceTime);
  checkError(cudaThreadExit());
}