//
// computes probability of no hadronic interaction as a function of b
// uses as input T_AA computed by program glauber.C
//


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

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


// my headers
#include "Global.h"

Double_t SigmaPP_INEL(Double_t sqrt_s)
{
  // from https://inspirehep.net/record/1631672?ln=en
  Double_t A = 25.0; // +/-0.9
  Double_t B = 0.146; // +/-0.004
  Double_t s = sqrt_s*sqrt_s;
  return (A+B*TMath::Log(s)*TMath::Log(s));
}


Double_t SigmaPP_SL(Double_t sqrt_s)
{
  // from PDG 2016
  Double_t H = 0.272;
  Double_t P = 34.41;
  Double_t R1 = 13.07;
  Double_t R2 = 7.39;
  Double_t eta1 = 0.4473;
  Double_t eta2 = 0.5486;
  Double_t M = 2.12;
  Double_t s_abM = TMath::Power(2.0*ProtonMass+M,2.0);
  Double_t sR = sqrt_s*sqrt_s/s_abM;
  Double_t F1 = H*TMath::Power(TMath::Log(sR),2.0)+P;
  Double_t F2 = R1*TMath::Power(sR,-eta1);
  Double_t F3 = R2*TMath::Power(sR,-eta2);  
  return (F1+F2+F3);
}

void ProbNoInt(Int_t opt = 0, Double_t sqrt_s = 2760)
// if opt == 0, Pb; opt == 1, Au; opt == 2, Xe
// sqrt_s of interaction

{
  // get profile from file
  TFile *fIn = NULL;
  if (opt == 0) fIn = new TFile("TPbPb.root");
  else if (opt == 1) fIn = new TFile("TAuAu.root");
  else if (opt == 2) fIn = new TFile("TXeXe.root");
  TH1 *h = (TH1*) fIn->Get("OverlapFunction");
  TH1 *hclon = (TH1*) h->Clone("hclon");
  TH1 *hProbNoInt = (TH1*) h->Clone("hProbNoInt");  


  // compute prob of no interaction as a function of b
  Int_t nb = hclon->GetXaxis()->GetNbins();   // get number of bins
  Double_t sigmaNNmb = SigmaPP_INEL(sqrt_s);
  Double_t sigmaNN = 0.1*sigmaNNmb; // now in fm^2
  for(Int_t i=0;i<nb;i++) {
    Double_t b = hclon->GetBinCenter(i+1);
    Double_t TAAb = hclon->GetBinContent(i+1);
    Double_t ProbNoInt = TMath::Exp(-TAAb*sigmaNN);
    hProbNoInt->SetBinContent(i+1,ProbNoInt);
  }


  // save
  char tmp[120];
  if (opt == 0) sprintf(tmp,"ProbNoInt_PbPb_%i.root",((Int_t) sqrt_s));
  else if (opt == 1) sprintf(tmp,"ProbNoInt_AuAu_%i.root",((Int_t) sqrt_s));
  else if (opt == 2) sprintf(tmp,"ProbNoInt_XeXe_%i.root",((Int_t) sqrt_s));  
  TFile *fOut = new TFile(tmp,"recreate");
  
  fOut->cd();
  hProbNoInt->Write();
  fOut->Close();
  
  // plot
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  TCanvas *c = new TCanvas("c","c",600,400);
  gPad->SetLeftMargin(0.09);
  gPad->SetRightMargin(0.05);
  gPad->SetTopMargin(0.025);//0.025
  gPad->SetBottomMargin(0.13);

  hProbNoInt->SetLineColor(TColor::GetColor("#000099"));

  hProbNoInt->Draw("l");
  hProbNoInt->SetLineWidth(2);
  hProbNoInt->SetXTitle("b (GeV)");
  hProbNoInt->SetYTitle("P_{NH}(b)");

  if (opt == 0) sprintf(tmp,"ProbNoInt_PbPb_%i.eps",((Int_t) sqrt_s));
  else if (opt == 1) sprintf(tmp,"ProbNoInt_AuAu_%i.eps",((Int_t) sqrt_s));
  else if (opt == 2) sprintf(tmp,"ProbNoInt_XeXe_%i.eps",((Int_t) sqrt_s));  
  c->Print(tmp);

}