///////////////////////////////////////////////////////////////////////////////
// File:  MyCherenkov.cc
// Description: Find number of Cherenkov photons in a medium
///////////////////////////////////////////////////////////////////////////////

#include "MyCherenkov.hh"
#include "G4DynamicParticle.hh"
#include "G4SystemOfUnits.hh"
#include "TMath.h"
#include <iostream>

int MyCherenkov::computeNPE(const G4Step* aStep, const double& ref_index, const bool& verbose) {
  const G4Track *aTrack = aStep->GetTrack();
  const G4ParticleDefinition *particleDef = aTrack->GetDefinition();
  if (!isApplicable(pDef, verbose)) {
    return 0;
  }

  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
  double energy = aParticle->GetTotalEnergy();
  double momentum = aParticle->GetTotalMomentum();
  double pBeta = momentum / energy;
  double step_length = aStep->GetStepLength();
  if (pBeta < (1 / ref_index) || step_length < 0.0001) {
    if (verbose)
      std::cout << "MyCherenkov::computeNPE: pBeta " << pBeta << " 1/mu " << (1 / ref_index) << " step_length " << step_length << std::endl;
    return 0;
  }

  int nbOfPhotons = computeNbOfPhotons(pBeta, step_length, ref_index, verbose);
  if (verbose)
    std::cout << "MyCherenkov::computeNPE: pBeta " << pBeta << " step_length " << step_length << " nbOfPhotons " << nbOfPhotons << std::endl;
  return std::max(nbOfPhotons, 0);
}

int MyCherenkov::computeNbOfPhotons(const double& pBeta, const double& step_length, const double& ref_index, const bool& verbose) {
  constexpr double lambda1 = (280.0 / 1.e+7) * CLHEP::cm;
  constexpr double lambda2 = (700.0 / 1.e+7) * CLHEP::cm;
  double alpha = 0.0073;
  double theta_C = acos(1. / (pBeta * ref_index));
  double lambdaDiff = (1. / lambda1 - 1. / lambda2);
  double cherenPhPerLength = 2 * M_PI * alpha * lambdaDiff * CLHEP::cm;
  double d_NOfPhotons = cherenPhPerLength * sin(theta_C) * sin(theta_C) * (step_length / CLHEP::cm);
  int nbOfPhotons = int(d_NOfPhotons);
  if (verbose)
    std::cout << "MyCherenkov::computeNbOfPhotons: StepLength " << step_length << " theta_C " << theta_C << " lambdaDiff " << lambdaDiff << " cherenPhPerLength " << cherenPhPerLength << " Photons " << d_NOfPhotons << " " << nbOfPhotons << std::endl;
  return nbOfPhotons;
}

bool MyCherenkov::isApplicable(const G4ParticleDefinition* aParticleType, const bool& verbose) {
  bool tmp = (aParticleType->GetPDGCharge() != 0);
  if (verbose)
    std::cout << "MyCherenkov::isApplicable: aParticleType " << aParticleType->GetParticleName() << " PDGCharge " << aParticleType->GetPDGCharge() << " Result " << tmp << std::endl;
  return tmp;
}