//***********************************************************************************************************
// BinToStd_proton_position.C
// Root command file
// Type: root BinToStd_proton_position.C
//
// Read the X-ray output file that is generated by Geant4 tomography
// simulation. It reads gamma information, either at creation, or at exit, and rewrite the events
// in a binary file StimEvent_std.DAT
//
// More information is available in UserGuide
// Created by Z.LI LP2i Bordeaux 2022
//***********************************************************************************************************

#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>

#include <vector>
// using namespace std;

bool IsEqual(double a, double b, double eps, double releps)
{
  if (a == b) {
    return true;
  }

  if (fabs(a - b) <= releps * fabs(b)) {
    return true;
  }

  if (fabs(a - b) < eps) {
    return true;
  }

  return false;
}
double eps = 1e-20;  // absolut difference
double releps = 1e-10;  // relative difference

// Define a structure to read and write each event in the required binary format
struct StimEvent
{
  uint16_t energy_keV;  // different from Pixe Event, it is in keV
  uint16_t pixelIndex;
  uint16_t sliceIndex;
  uint8_t projectionIndex;
};
struct ParticleInfo
{
  float energy_keV;
  float mx;
  float my;
  float mz;
  float x;
  float y;
  float z;
};
struct RunInfo
{
  // uint_16t
  uint8_t projectionIndex;  // 1 byte
  uint16_t sliceIndex;  //
  uint16_t pixelIndex;
  uint32_t nbParticle;  // 4 bytes int
};
struct Point
{
  double m_x;
  double m_y;
  double m_z;
};

bool IsDetected(Point poi1, Point poi2, double theta)
{
  double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
             / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
             / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
  if (a > 1.0) a = 1;
  if (a < -1.0) a = -1;
  double r = acos(a);
  if (r > theta)
    return false;
  else {
    // printf("      acos: %f, radius: %f\n", r, theta);
    return true;
  }
}

bool IsDetected_position(Point poi1, Point poi2, double r)
{
  double a = sqrt((poi1.m_x - poi2.m_x) * (poi1.m_x - poi2.m_x)
                  + (poi1.m_y - poi2.m_y) * (poi1.m_y - poi2.m_y)
                  + (poi1.m_z - poi2.m_z) * (poi1.m_z - poi2.m_z));

  // if(a <= r)    return true;
  if (a > r)
    return false;

  else {
    // printf("      distance of two points: %f, radius: %f\n", a, r);
    return true;
  }
}

void BinToStd_proton_position()
{
  // printf("%f %f %f\n", acos(1), acos(-1), acos(0));
  // return;

  //***********************************************************************
  //**************************Detection parameters (begin)*****************
  //***********************************************************************

  const int nbProjection = 1;
  const int nbSlice = 1;
  const int nbPixel = 1;
  double totalAngleSpan = 180.;  // in degree

  double angleOfDetector = 0.;  // angle of detector relative to the incident

  double distanceObjectDetector = 22000.;  // um

  // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
  // angle of the right circular cone in radian
  double theta = 10.2 * TMath::DegToRad();  // in radian
  double radiusOfDetector = distanceObjectDetector * tan(theta);
  bool usePosition = true;

  //***********************************************************************
  //**************************Detection parameters (end)*******************
  //***********************************************************************

  FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
  FILE* out = fopen("../build/StimEvent_std", "wb");

  if (input == NULL) {
    printf("error for opening the input file\n");
    return;
  }

  RunInfo runInfo;
  StimEvent stimEvent;
  Point centerOfDetector;
  Point protonMomentum;
  Point protonPosition;
  Point intersectionPoint;
  long long count = 0;
  int runID = -1;  // index of simulations, namely runID, starting from 0

  // while(!feof(input)) //if not the end, read
  while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
    runID++;
    int nbParticle = runInfo.nbParticle;

    //(begin)****************************************************************
    // the following codes are used only when in the simulation
    // the index of projection, slice and pixel is not
    // correctly configured
    runInfo.projectionIndex = runID / (nbSlice * nbPixel);
    int remain = runID % (nbSlice * nbPixel);
    runInfo.sliceIndex = remain / nbPixel;
    runInfo.pixelIndex = remain % nbPixel;

    //(end)*******************************************************************

    //***********************************************************************
    //**************************Print information (begin)********************
    //***********************************************************************

    printf(
      "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
      "nbParticle = %d\n",
      runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);

    //***********************************************************************
    //**************************Print information (end)**********************
    //***********************************************************************

    if (!nbParticle) continue;
    std::vector<ParticleInfo> protonAtExit(nbParticle);
    fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input);

    // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
    // the angle between source direction and detector, which should be constant
    // when source is rotating
    double ra = TMath::DegToRad()
                * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
    centerOfDetector.m_x = distanceObjectDetector * cos(ra);
    centerOfDetector.m_y = distanceObjectDetector * sin(ra);
    centerOfDetector.m_z = 0;

    for (int i = 0; i < nbParticle; ++i) {
      // proton selection: energy should be lower than 4095 keV
      if (protonAtExit[i].energy_keV >= 4095) continue;

      protonMomentum.m_x = protonAtExit[i].mx;
      protonMomentum.m_y = protonAtExit[i].my;
      protonMomentum.m_z = protonAtExit[i].mz;

      if (!usePosition) {
        if (!IsDetected(centerOfDetector, protonMomentum, theta)) continue;
      }
      else {
        double c =
          distanceObjectDetector * (protonMomentum.m_x * cos(ra) + protonMomentum.m_y * sin(ra));
        if (IsEqual(0, c, eps, releps)) continue;  // parallel

        protonPosition.m_x = protonAtExit[i].x;
        protonPosition.m_y = protonAtExit[i].y;
        protonPosition.m_z = protonAtExit[i].z;

        double t = (distanceObjectDetector * distanceObjectDetector
                    - protonPosition.m_x * distanceObjectDetector * cos(ra)
                    - protonPosition.m_y * distanceObjectDetector * sin(ra))
                   / c;

        intersectionPoint.m_x = protonPosition.m_x + protonMomentum.m_x * t;
        intersectionPoint.m_y = protonPosition.m_y + protonMomentum.m_y * t;
        intersectionPoint.m_z = protonPosition.m_z + protonMomentum.m_z * t;

        if (!IsDetected_position(centerOfDetector, intersectionPoint, radiusOfDetector)) continue;

        // printf("      t = %f, intersection point: (%f, %f, %f) centor of detector: (%f, %f, %f)
        // 111=%f, 222=%f \n", t, intersectionPoint.m_x,intersectionPoint.m_y,intersectionPoint.m_z,
        // centerOfDetector.m_x,centerOfDetector.m_y,centerOfDetector.m_z,
        // (distanceObjectDetector*distanceObjectDetector-protonPosition.m_x*distanceObjectDetector*cos(ra)
        // -protonPosition.m_y*distanceObjectDetector*sin(ra)), c);

        // printf("      distanceObjectDetector = %f, protonPosition.m_x=%f,
        // distanceObjectDetector*cos(ra)=%f, protonPosition.m_y=%f,
        // distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, protonPosition.m_x,
        // distanceObjectDetector*cos(ra),
        // protonPosition.m_y,
        // distanceObjectDetector*sin(ra));

        double tt = (intersectionPoint.m_x - protonPosition.m_x) * protonMomentum.m_x
                    + (intersectionPoint.m_y - protonPosition.m_y) * protonMomentum.m_y
                    + (intersectionPoint.m_z - protonPosition.m_z) * protonMomentum.m_z;
        if (tt < 0) continue;
      }

      stimEvent.energy_10eV = floor(100 * protonAtExit[i].energy_keV + 0.5);
      stimEvent.projectionIndex = runInfo.projectionIndex;
      stimEvent.sliceIndex = runInfo.sliceIndex;
      stimEvent.pixelIndex = runInfo.pixelIndex;
      fwrite(&stimEvent, 7, 1, out);
      count++;

      //***********************************************************************
      //**************************Print information
      //(begin)********************
      //***********************************************************************

      if (!usePosition) {
        printf(
          "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
          "(%f, %f, %f), energy: %f keV\n",
          i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
          protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
      }
      else {
        // printf("---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
        // momentum: (%f, %f, %f), energy: %f keV, position: (%f, %f, %f)\n", i, runID,
        // runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, protonAtExit[i].mx,
        // protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV, protonAtExit[i].x,
        // protonAtExit[i].y, protonAtExit[i].z);

        printf(
          "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
          "(%f, %f, %f), energy: %f keV\n",
          i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
          protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
      }

      //***********************************************************************
      //**************************Print information (end)**********************
      //***********************************************************************
    }
  }
  printf(
    "\n---------------Number of StimEvent in total: "
    "%lld------------------------\n",
    count);
  fclose(input);
  fclose(out);

  // FILE* input2;
  // input2 = fopen("StimEvent_std.DAT","rb");
  // StimEvent p;
  // double eventId = -1;
  // while(fread(&p, 7, 1, input2))
  // {

  // if(p.projectionIndex == 8 &&p.sliceIndex ==64 && p.pixelIndex==64)
  // {
  // eventId++;
  // printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
  // eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);

  // }

  // }
  // fclose(input2);
}