#ifndef SPLITPOLEHit_H
#define SPLITPOLEHit_H

#include "Isotope.h"

#include <cmath>
#include <random>

// static double randZeroToOne() {
//     return static_cast<double>(rand()) / RAND_MAX;
// }

// // Box-Muller transform to generate random Gaussian numbers
// static double generateGaussian(double mean, double stddev) {
//     static bool hasSpare = false;
//     static double spare;
//     if (hasSpare) {
//         hasSpare = false;
//         return mean + stddev * spare;
//     } else {
//         double u, v, s;
//         do {
//             u = 2.0 * randZeroToOne() - 1.0;
//             v = 2.0 * randZeroToOne() - 1.0;
//             s = u * u + v * v;
//         } while (s >= 1.0 || s == 0.0);

//         s = std::sqrt(-2.0 * std::log(s) / s);
//         spare = v * s;
//         hasSpare = true;
//         return mean + stddev * u * s;
//     }
// }

namespace SPS{
  namespace ChMap{

    const short ScinR = 0;
    const short ScinL = 1;
    const short dFR = 9;
    const short dFL = 8;
    const short dBR = 11;
    const short dBL = 10;
    const short Cathode = 7;
    const short AnodeF = 13;
    const short AnodeB = 15;

  };

  const double c = 299.792458; // mm/ns
  const double pi = M_PI;
  const double deg2rad = pi/180.;

  const double DISPERSION = 1.96; // x-position/rho
  const double MAGNIFICATION = 0.39; // in x-position
  const double X1X2Separation = 36; // cm

}

class SplitPoleHit{

public:
  SplitPoleHit(){
    ClearData();
  }

  unsigned int eSR;   unsigned long long tSR;
  unsigned int eSL;   unsigned long long tSL;
  unsigned int eFR;   unsigned long long tFR;
  unsigned int eFL;   unsigned long long tFL;
  unsigned int eBR;   unsigned long long tBR;
  unsigned int eBL;   unsigned long long tBL;
  unsigned int eCath; unsigned long long tCath;
  unsigned int eAF;   unsigned long long tAF;
  unsigned int eAB;   unsigned long long tAB;

  float eSAvg;
  float x1, x2, theta;
  double xAvg;

  double GetQ0() const {return Q0;}
  double GetRho0() const {return rho0;}
  double GetZoffset() const {return zOffset;}

  void SetMassTablePath(std::string path){
    target.SetMassTablePath(path);
    beam.SetMassTablePath(path);
    recoil.SetMassTablePath(path);
    heavyRecoil.SetMassTablePath(path);
  }

  void CalConstants(std::string targetStr, std::string beamStr, std::string recoilStr, double energyMeV, double angleDeg){
    target.SetIsoByName(targetStr);
    beam.SetIsoByName(beamStr);
    recoil.SetIsoByName(recoilStr);
    heavyRecoil.SetIso(target.A + beam.A - recoil.A, target.Z + beam.Z - recoil.Z);

    angleDegree = angleDeg; // degree
    beamKE = energyMeV; // MeV

    Ei = target.Mass + beamKE + beam.Mass;
    k1 = sqrt( 2*beam.Mass*beamKE + beamKE*beamKE);
    cs = cos(angleDegree * SPS::deg2rad);
    ma = recoil.Mass;
    mb = heavyRecoil.Mass;

    isConstantCal = true;

    printf("============================================\n");
    printf("  Beam Mass : %20.4f MeV/c2\n", beam.Mass);
    printf("    beam KE : %20.4f MeV\n", beamKE);
    printf("Target Mass : %20.4f MeV/c2\n", target.Mass);
    printf("Recoil Mass : %20.4f MeV/c2\n", recoil.Mass);
    printf("H. Rec Mass : %20.4f MeV/c2\n", heavyRecoil.Mass);

    printf("      angle : %20.4f deg\n", angleDegree);
    printf("         k1 : %20.4f MeV/c\n", k1);
    printf("         Ei : %20.4f MeV\n", Ei);

  }

  double CalRecoilMomentum(double Ex){

    if( !isConstantCal ) return 0;

    double p =  Ei*Ei - k1*k1;
    double q = ma*ma - (mb + Ex)*(mb + Ex);

    double x = k1* ( p + q) * cs;
    double y = pow( p, 2) + pow(q, 2)- 2 * Ei * Ei * (ma* ma + (mb + Ex)*(mb + Ex)) + 2 * k1 * k1 * (ma*ma * cos(2* angleDegree * SPS::deg2rad) + (mb + Ex)*(mb + Ex));
    double z = 2 * ( Ei*Ei - k1*k1 * cs * cs) ;

    return (x + Ei * sqrt(y))/z;

  }

  double Momentum2Rho(double ka){
    return ka / (recoil.Z * Bfield * SPS::c);
  }

  double Momentum2Ex(double ka){
    return sqrt( Ei*Ei - k1*k1 + ma*ma + 2 * cs * k1 * ka - 2*Ei*sqrt(ma*ma + ka*ka)) - mb;
  }

  double Rho2Ex(double rhoInM){
    double ka = rhoInM * (recoil.Z * Bfield * SPS::c);
    return Momentum2Ex(ka);
  }

  void CalZoffset(double magFieldinT){

    Bfield = magFieldinT;

    if( !isConstantCal ) return;

    double recoilP = CalRecoilMomentum(0);

    Q0 = target.Mass + beam.Mass - recoil.Mass - heavyRecoil.Mass;

    double recoilKE = sqrt(ma*ma + recoilP* recoilP) - ma;

    printf("Q value : %f \n", Q0);
    printf("recoil enegry for ground state: %f MeV = %f MeV/c\n", recoilKE, recoilP);

    rho0 = recoilP/(recoil.Z * Bfield * SPS::c); // in m

    double haha = sqrt( ma * beam.Mass * beamKE / recoilKE );
    double k = haha * sin(angleDegree * SPS::deg2rad) / ( ma + mb - haha * cs);


    zOffset = -100.0 * rho0 * k * SPS::DISPERSION * SPS::MAGNIFICATION;

    printf("rho: %f m; z-offset: %f cm\n", rho0, zOffset);

  }

  void ClearData(){
    eSR = 0;  tSR = 0;
    eSL = 0;  tSL = 0;
    eFR = 0;  tFR = 0;
    eFL = 0;  tFL = 0;
    eBR = 0;  tBR = 0;
    eBL = 0;  tBL = 0;
    eCath = 0;  tCath = 0;
    eAF = 0;  tAF = 0;
    eAB = 0;  tAB = 0;

    eSAvg = -1;
    x1 = NAN;
    x2 = NAN;
    theta = NAN;
    xAvg = NAN;

    isConstantCal = false;
  }

  void CalData(float scale = 2.){

    if( eSR  > 0 && eSL  > 0 ) eSAvg = (eSR + eSL)/2;
    if( eSR  > 0 && eSL == 0 ) eSAvg = eSR;
    if( eSR == 0 && eSL  > 0 ) eSAvg = eSL;

    if( tFR > 0 && tFL > 0 ) {
      if( tFL > tFR) x1 = (tFL - tFR)/scale/2.1;
      if( tFL < tFR) x1 = (tFR - tFL)/scale/-2.1;
    }
    if( tBR > 0 && tBL > 0 ) {
      if( tBL > tBR) x2 = (tBL - tBR)/scale/1.98;
      if( tBR > tBL) x2 = (tBR - tBL)/scale/-1.98;
    }

    // printf("x1: %f, x2 : %f \n", x1, x2);

    if( !std::isnan(x1)  && !std::isnan(x2)) {

      if( x2 > x1 ) {
        theta = atan((x2-x1)/SPS::X1X2Separation);
      }else if(x2 < x1){
        theta = SPS::pi + atan((x2-x1)/SPS::X1X2Separation);
      }else{
        theta = SPS::pi * 0.5;
      }

      double w1 = 0.5 - zOffset/4.28625;
      xAvg = w1 * x1 + (1-w1)* x2;

    }

  }

private:

  Isotope target;
  Isotope beam;
  Isotope recoil;
  Isotope heavyRecoil;

  double Bfield;
  double angleDegree;
  double beamKE;

  double zOffset;
  double Q0, rho0;

  bool isConstantCal;

  double Ei, k1, cs, ma, mb;

};

#endif