// GunPMTGamma.cc
// Contact person: Sean Grullon <grullons@hep.upenn.edu>
// See GunPMTGamma.hh for more details
//———————————————————————//

#include <cstdio>
#include <globals.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleTable.hh>
#include <G4IonTable.hh>
#include <G4Ions.hh>
#include <G4OpticalPhoton.hh>
#include <G4Event.hh>
#include <G4Track.hh>
#include <G4HEPEvtParticle.hh>
#include <G4Gamma.hh>

#include <Randomize.hh>
#include <RAT/DB.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/GunPMTGamma.hh>
#include <RAT/Log.hh>

#include <string.h>  // for strcmp
#include <sstream>
#include <TF1.h>
#include <TH2F.h>
#include <TH1F.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <TVector3.h>

using namespace CLHEP;
using namespace std;

GLG4VertexGen_GunPMTGamma::GLG4VertexGen_GunPMTGamma(const char *arg_dbname)
  :GLG4VertexGen(arg_dbname),fKe(0.0)
{
  G4ParticleDefinition *gammaDefinition=G4Gamma::GammaDefinition();
  fPDef= gammaDefinition;
}

GLG4VertexGen_GunPMTGamma::~GLG4VertexGen_GunPMTGamma()
{
}

double PMTGDistFn(double cosbeta, double lambda, double radiusGen, double radiusPSUP);

// Following necessary so I can make the function into a ROOT TF1 function
double PMTGDistFnRoot(double* cosbetaptr, double* par) {
  // &cosbetaptr =  value of cos(beta) at which function is to be calculated
  // par[0] = scale factor
  // par[1] = attenuation length
  // par[2] = inner radius
  // par[3] = PSUP radius

  double lambda = par[1];
  double radiusGen = par[2];
  double radiusPSUP = par[3];
  double cosbeta = *cosbetaptr;
  double fn = PMTGDistFn(cosbeta, lambda, radiusGen, radiusPSUP);
  return  par[0]*fn;
}


// Following is the function that calculates the value of the gamma direction
// distribution given cos(beta).  It does not normalize so the integral is one.
double PMTGDistFn(double cosbeta, double lambda, double radiusGen, double radiusPSUP) {
  // cosbeta =  value of cos(beta) at which function is to be calculated
  // lambda = attenuation length
  // radiusGen = inner radius
  // radiusPSUP = PSUP radius

  double Rsq =radiusPSUP*radiusPSUP;
  double rsq =radiusGen*radiusGen;
  double costheta,difratio,s,decayfactor,f;
  double ssq=0.;
  if(cosbeta==0.) {
    difratio = radiusPSUP/sqrt(Rsq-rsq);
  }else{
    double sign = 1.;
    if(cosbeta<0.) sign = -1.;
    // Rsq, bq, and cq are the a, b, and c in the quadratic formula
    double bq = -2.*radiusGen*radiusPSUP*(1.-cosbeta*cosbeta);
    double cq = rsq-cosbeta*cosbeta*(rsq+Rsq);
    double root = sqrt( bq*bq-4.*Rsq*cq );
    costheta = (-bq+ sign*root)/(2*Rsq);
    ssq = Rsq+rsq-2.*radiusGen*radiusPSUP*costheta;
    difratio=(costheta - (1.-costheta*costheta)*radiusGen*radiusPSUP/ssq)*Rsq/(cosbeta*ssq);
  }
  s = sqrt (ssq);
  decayfactor = exp(-s/lambda);
  f = decayfactor/difratio;

  return  f;
}


void GLG4VertexGen_GunPMTGamma::GeneratePrimaryVertex(G4Event *argEvent, G4ThreeVector &dx, G4double dt)
{

  G4PrimaryVertex* vertex= new G4PrimaryVertex( dx, dt );

  G4double mass= fPDef->GetPDGMass();

  G4PrimaryParticle* particle;

  G4ThreeVector dir;

  //first, find u... vector pointing to center of AV

  TVector3 mcpos(dx.x(),dx.y(),dx.z());
  TVector3 rhat=-mcpos.Unit();

  // used to normalize later for orthogonal vector
  TVector3 rhat2=-mcpos.Unit();
  // next, find random orthogonal vector
  TVector3 randvec(G4UniformRand(),G4UniformRand(),G4UniformRand());
  rhat2*=rhat2.Dot(randvec)/(rhat2.Dot(rhat2));
  TVector3 orthvec=randvec-rhat2;

  // next, sample randomly from Bill's angular distribution to get alpha

  // compton scattering length; best fit from Bill's Model (in mm)
  Double_t lambdapmt = 233.3;

  // The PMT beta/gamma angular distribution specifically depends
  // on the position of the generated particle.  Grab MC radial position
  // and calculate the angular distribution from Bill's function.

  // in mm
  Double_t rpmt=mcpos.Mag();
  // Radius of PSUP in Bill's model fit, in mm
  Double_t Rpmt = 8500.;
  Double_t bvalpmt = 1.0;
  Double_t pvalpmt[4] = {bvalpmt,lambdapmt, rpmt, Rpmt};
  Double_t  maxValuepmt = PMTGDistFnRoot(&bvalpmt, pvalpmt);

  // Angular distribution from PMT beta gamma events as a function
  // of radius. See docdb 1712.
  TF1* fdistalpha = new TF1("fdistalpha",PMTGDistFnRoot,0.0,+1.,4);
  fdistalpha->SetParameters( 1./maxValuepmt, lambdapmt, rpmt, Rpmt);
  fdistalpha->SetMinimum(0.);
  double alph=acos(fdistalpha->GetRandom());

  //finally, rotate rhat by alpha about random orthogonal vector
  rhat.Rotate(alph,orthvec);

  dir = G4ThreeVector(rhat.X(),rhat.Y(),rhat.Z());
  G4ThreeVector rmom( dir * sqrt(fKe*(fKe+2.*mass)) );
  particle= new G4PrimaryParticle(fPDef,        // particle code
				  rmom.x(),     // x component of momentum
				  rmom.y(),     // y component of momentum
				  rmom.z() );   // z component of momentum



  G4ThreeVector fPol(0.0,0.0,0.0);
  // adjust polarization for gammas
  // spin 1 mass 0 particles should have transverse polarization
  G4ThreeVector rpol= fPol - dir * (dir*fPol);
  G4double rpolmag= rpol.mag();
  if (rpolmag > 0.0) {
    // use projected orthogonal pol, normalized
    rpol *= (1.0/rpolmag);
  }
  else {
    // choose random pol
    G4double phi= (G4UniformRand()*2.0-1.0)*pi;
    G4ThreeVector e1= dir.orthogonal().unit();
    G4ThreeVector e2= dir.cross(e1);
    rpol= e1*cos(phi)+e2*sin(phi);
  }
  particle->SetPolarization(rpol.x(), rpol.y(), rpol.z());


  particle->SetMass(mass); // Geant4 is silly.

  vertex->SetPrimary( particle );

  argEvent->AddPrimaryVertex(vertex);

}

///////////////// End Generate Primary Vertex

void GLG4VertexGen_GunPMTGamma::
SetState(G4String newValues)
{
  newValues = util_strip_default(newValues);
  if (newValues.length() == 0) {
    // print help and current state
    RAT::info << "Current state of this GLG4VertexGen_GunPMTGamma:\n"
	      << " \"" << GetState() << "\"\n" << newline;
    RAT::info << "Format of argument to GLG4VertexGen_GunPMTGamma::SetState: \n"
      " \"KE_MeV \"\n"
      " where KE_MeV is the kinetic energy\n"
	      << newline;
    return;
  }

  std::istringstream is(newValues.c_str());

  // set particle
  //  std::string pname;
  //  pname="gamma";

  if (is.fail())
    return;

  //  G4ParticleDefinition * newTestGunG4Code =   G4ParticleTable::GetParticleTable()->FindParticle(pname);

  // set energy
  G4double x;
  is >> x;
  if (is.fail()) {
    return;
  }
  if (x <= 0.0) {
    fKe= 0.0;
  }
  else {
    fKe= x*MeV;
  }

}


G4String GLG4VertexGen_GunPMTGamma::
GetState()
{
  std::ostringstream os;

  os << fPDef->GetParticleName() << "    "
     << fKe << "    " << std::ends;
  G4String rv(os.str());

  return rv;
}


///////////////////////// end GunPMTGamma