//
// computes photon flux according to SL for different systems and
// impact parameter ranges
// Needs results from ProbNoInt.C

// c++ headers
#include <iostream>
#include <fstream>

// root headers
#include <TMath.h>
#include <TH2D.h>
#include <TFile.h>
#include <TF2.h>
#include <TCanvas.h>

// my headers
#include "Global.h"

// see FitFluxRatio
// corrects for using an approximate analytical results
// instead of a numerical calculation
Double_t p1 = -1.17; // middle between -1.20 and -1.14
Double_t p2 = 5.6; // middle between 5.8 and 5.4



//////////////////////////////////////////////
// Functions
//////////////////////////////////////////////


  
Double_t PhotonEnergy(Double_t Mass, Double_t Rapidity)
// 
// computes photon energy from the mass and rapidity of vector meson
//
{
  return 0.5*Mass*TMath::Exp(-Rapidity);
}

Double_t HardSphereFlux(Double_t PhotonEnergy, Double_t bmin=2.0*IonRadius)

//
// computes the flux from a nuclei using a hard sphere approximation
// return is (dN/dk) where k = Photon Energy
// see eq. (7) in arXiv 9902259

{
  Double_t F = 2.0*IonCharge*IonCharge*AlphaEM/(TMath::Pi()*PhotonEnergy);
  Double_t X = bmin*PhotonEnergy/(HbarC*LorentzBoost);
  Double_t K0 = TMath::BesselK0(X);
  Double_t K1 = TMath::BesselK1(X);
  
  return F*(X*K0*K1-0.5*X*X*(K1*K1-K0*K0));
}

Double_t GetN(Double_t k, Double_t r)
// get dN(k,r)/dkd^2r according to eq. 1 of arXiv 9902259
// k is in GeV and r in fermi
{
  Double_t x = k*r/(HbarC*LorentzBoost);
  Double_t F = x*x*IonCharge*IonCharge*AlphaEM/(TMath::Pi()*TMath::Pi()*k*r*r);
  Double_t K1 = TMath::BesselK1(x);
  Double_t K2 = TMath::BesselK0(x);  
  return F*(K1*K1+K2*K2/LorentzBoost);
}


//////////////////////////////////////////////
// Main
//////////////////////////////////////////////

Double_t Flux(Int_t opt = 0, Double_t y = 0.0, Bool_t kUPC = kTRUE)
// compute flux at rapidity y, and in the impact parameter range (b_min, b_max) in fm

{
  TFile *fIn = NULL;
  // Pb-Pb centrality class for Run1
  Double_t b_min = 13.05; // fm
  Double_t b_max = 14.96; // fm  
  if (opt == 0) { // lead, run1
  // load probability of no interaction
    fIn = new TFile("ProbNoInt_PbPb_2760.root");
    pBeamEnergy = 3500;
    BeamEnergy = pBeamEnergy*IonChargePb/IonAtomicNumberPb;
    IonRadius =  WSc_Pb;
    IonCharge = IonChargePb;
    b_min = 13.05; // fm
    b_max = 14.96; // fm
    mVM = JPsiMass; // GeV
  } else if (opt == 1) { // lead, run2
  // load probability of no interaction
    fIn = new TFile("ProbNoInt_PbPb_5020.root");
    pBeamEnergy = 6370;
    BeamEnergy = pBeamEnergy*IonChargePb/IonAtomicNumberPb;
    IonRadius =  WSc_Pb;
    IonCharge = IonChargePb;
    b_min = 13.1; // fm
    b_max = 15.0; // fm
    mVM = JPsiMass; // GeV
  } else if (opt == 2) { // gold
    fIn = new TFile("ProbNoInt_AuAu_200.root");
    pBeamEnergy = 250;
    BeamEnergy = pBeamEnergy*IonChargeAu/IonAtomicNumberAu;
    IonRadius =  WSc_Au;
    IonCharge = IonChargeAu;
    b_min = -1; // fm 
    b_max = -2; // fm
    mVM = JPsiMass; // GeV
  } else if (opt == 3) { // gold
    fIn = new TFile("ProbNoInt_AuAu_200.root");
    pBeamEnergy = 250;
    BeamEnergy = pBeamEnergy*IonChargeAu/IonAtomicNumberAu;
    IonRadius =  WSc_Au;
    IonCharge = IonChargeAu;
    b_min = -1; // fm 
    b_max = -2; // fm
    mVM = RhoMass; // GeV
  } else if (opt == 4) { // gold
    fIn = new TFile("ProbNoInt_AuAu_130.root");
    pBeamEnergy = 162;
    BeamEnergy = pBeamEnergy*IonChargeAu/IonAtomicNumberAu;
    IonRadius =  WSc_Au;
    IonCharge = IonChargeAu;
    b_min = -1; // fm 
    b_max = -2; // fm
    mVM = RhoMass; // GeV
  } else if (opt == 5) { // gold
    fIn = new TFile("ProbNoInt_AuAu_62.root");
    pBeamEnergy = 77.8;
    BeamEnergy = pBeamEnergy*IonChargeAu/IonAtomicNumberAu;
    IonRadius =  WSc_Au;
    IonCharge = IonChargeAu;
    b_min = -1; // fm 
    b_max = -2; // fm
    mVM = RhoMass; // GeV
  } else if (opt == 6) { // Xe
    fIn = new TFile("ProbNoInt_XeXe_5440.root");
    pBeamEnergy = 6500;
    BeamEnergy = pBeamEnergy*IonChargeAu/IonAtomicNumberAu;
    IonRadius =  Xe_R;
    IonCharge = IonChargeXe;
    b_min = -1; // fm 
    b_max = -2; // fm
    mVM = RhoMass; // GeV
  } else if (opt == 7) { // lead, run2
  // load probability of no interaction
    fIn = new TFile("ProbNoInt_PbPb_5020.root");
    pBeamEnergy = 6370;
    BeamEnergy = pBeamEnergy*IonChargePb/IonAtomicNumberPb;
    IonRadius =  WSc_Pb;
    IonCharge = IonChargePb;
    b_min = 13.1; // fm
    b_max = 15.0; // fm
    mVM = RhoMass; // GeV
  } else if (opt == 8) { // lead, run1
  // load probability of no interaction
    fIn = new TFile("ProbNoInt_PbPb_2760.root");
    pBeamEnergy = 3500;
    BeamEnergy = pBeamEnergy*IonChargePb/IonAtomicNumberPb;
    IonRadius =  WSc_Pb;
    IonCharge = IonChargePb;
    b_min = 13.05; // fm
    b_max = 14.96; // fm
    mVM = RhoMass; // GeV
  } else {
    cout << " option not known. bye. " << endl;
    return -1;
  }
  LorentzBoost = BeamEnergy/ProtonMass;
  TH1F *hProbNoInt = (TH1F*) fIn->Get("hProbNoInt");
  // read histogram info
  Double_t nb_h =  hProbNoInt->GetNbinsX();
  Double_t bmin_h = hProbNoInt->GetBinLowEdge(1);
  Double_t bmax_h = hProbNoInt->GetBinLowEdge(1+nb_h);
  Double_t bwidth_h = hProbNoInt->GetBinWidth(1);

  // define integral over r
  Int_t nr = 100;
  Double_t stepr = IonRadius/nr;
  // define integral over phi
  Int_t np = 50;
  Double_t stepp = TMath::Pi()/np;

  
  // to skip low probs
  Double_t prob_min = 0.001;

  // convert to photon energy y
  Double_t k = PhotonEnergy(mVM,y);
  Double_t b_integral = 0.0;
  for(Int_t ib=1;ib<=nb_h;ib++) { // loop over b
    Double_t prob = 0;
    if (kUPC) prob = hProbNoInt->GetBinContent(ib);
    if (!kUPC) prob = 1.0-hProbNoInt->GetBinContent(ib);
    Double_t b = hProbNoInt->GetBinCenter(ib);
    if (prob<prob_min) continue;  // skip low probabilities
    if (!kUPC && (b_min > b)) continue;
    if (!kUPC && (b_max < b)) break;
    Double_t r_integral=0;
    for(Int_t ir=0;ir<nr;ir++){
      Double_t r = ir*stepr+0.5*stepr; // center of r-bin
      Double_t phi_integral=0;
      for(Int_t ip=0;ip<np;ip++){
	Double_t phi = ip*stepp+0.5*stepp;
	Double_t d = b+r*TMath::Cos(phi);
	// correct for using analytical instead of numerical flux
	Double_t s = 1.0/(1+TMath::Exp(p1*(d-p2))); 
	Double_t Nkd = s*GetN(k,d);
	phi_integral+=Nkd;
      } // loop over phi
      phi_integral*=(2.0*stepp); // factor of two because the integral is from 0 to 2pi
      r_integral+=(r*phi_integral);
    } // loop over r
    r_integral*=(stepr/(TMath::Pi()*IonRadius*IonRadius));
    b_integral+=(b*prob*r_integral);
  } // loop over b
  b_integral*=(2*TMath::Pi()*bwidth_h);
  
  // use analytic result for the rest  
  Double_t hsf  = 0;
  if (kUPC) hsf = HardSphereFlux(k,bmax_h);
  return k*(b_integral+hsf);
} 

void Make_Table()
{

  // run1 LHC
  cout << " LHC Run1: " << endl;
  for(Double_t y=-4.1; y<4.11;y+=0.1) {
    Double_t nUPC = Flux(0,y,kTRUE);
    Double_t nPER = Flux(0,y,kFALSE);    
    cout <<  y << " " << nUPC << " " << nPER << endl;
  }


  
  // run2 LHC
  cout << " LHC Run2: " << endl;
  for(Double_t y=-4.1; y<4.11;y+=0.1) {
    Double_t nUPC = Flux(1,y,kTRUE);
    Double_t nPER = Flux(1,y,kFALSE);    
    cout <<  y << " " << nUPC << " " << nPER << endl;
  }
  
}