//
// computes overlap and thickness functions for Au or Pb
// Stores them in histograms which later are used
// to compute no-interaction probabilities
//

// 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"

static TF1* f_wdsx_z;
static TF1* f_ta;
static TF2* f_ta_rfi;

const Double_t bmax = 20.0;

Double_t wdsx_z(Double_t* x, Double_t* par)
{
//
//  Wood Saxon Parameterisation
//  as a function of z for fixed b
//
    Double_t bb  = par[5];
    Double_t zz  = x[0];
    Double_t r0  = par[0];
    Double_t d   = par[1];
    Double_t w   = par[2];
    Double_t n   = par[3];
    Double_t b2  = par[4];

    Double_t r = TMath::Sqrt(bb*bb+zz*zz);
    Double_t costh = zz/r;
    Double_t Y20 = TMath::Sqrt(5.0/(16.0*TMath::Pi()))*(3.0*costh*costh-1.0);
    Double_t rho  = n * (1.+w*(r/r0)*(r/r0))/(1.+TMath::Exp((r-r0*(1+b2*Y20))/d));
    return rho;
}

Double_t ta(Double_t* x, Double_t* par)
{
//
//  Thickness function 
//
    Double_t b  = x[0];
    f_wdsx_z->SetParameter(5,b);
    Double_t y  = 2.*f_wdsx_z->Integral(0.,bmax);
    return y;
}

Double_t ta_rfi(Double_t* x, Double_t* par)
{
//
//  Kernel for overlap function
//
    Double_t b    = par[0];
    Double_t r1   = x[0];
    Double_t phi  = x[1];
    Double_t r2   = TMath::Sqrt(r1*r1+b*b-2.*r1*b*TMath::Cos(phi));
    Double_t y    = r1*f_ta->Eval(r1)*f_ta->Eval(r2);
//    Double_t y    = r1*f_ta->Eval(r1);
    return y;
}


Double_t taa(Double_t* x, Double_t* par)
{
//
//  Overlap function
//

    
    Double_t b    = x[0];
    f_ta_rfi->SetParameter(0,b);
    f_ta_rfi->Update();
    
    Double_t y    = 2.*
	f_ta_rfi->Integral(0.,bmax, 0., TMath::Pi(), 0.000001);
    //    printf("taa %f %f\n", x[0], y);
    return y;
}


void glauber(Int_t opt = 0)
// if opt == 0, Pb; opt == 1,  Au; opt == 2, Xe
{
  Int_t nBin = 200;
  
//
//  Calculates some geometrical properties of PbPb collisions 
//  in the Glauber Model

//
//  Wood Saxon-nuclear density (z, for fixed b)
//
    TCanvas *c3 = new TCanvas("c3","Wood Saxon",400,10,600,700);
    f_wdsx_z = new TF1("f_wdsx_z", wdsx_z, 0, bmax, 6);
    f_wdsx_z->SetNpx(nBin);
    if (opt == 0) {
      f_wdsx_z->SetParameter(0,WSc_Pb);
      f_wdsx_z->SetParameter(1,WSz_Pb);
      f_wdsx_z->SetParameter(3,WSrho0_Pb);
      f_wdsx_z->SetParameter(4,0.);
    } else if (opt == 1){
      f_wdsx_z->SetParameter(0,WSc_Au);
      f_wdsx_z->SetParameter(1,WSz_Au);
      f_wdsx_z->SetParameter(3,WSrho0_Au);
      f_wdsx_z->SetParameter(4,0.);
    }  else if (opt == 2){
      f_wdsx_z->SetParameter(0,Xe_R);
      f_wdsx_z->SetParameter(1,Xe_a);
      f_wdsx_z->SetParameter(3,Xe_rho0);
      f_wdsx_z->SetParameter(4,Xe_b2);
    }
    f_wdsx_z->SetParameter(2,0.000);
    f_wdsx_z->SetParameter(5,0.); // b=0
    f_wdsx_z->Draw();
    
//
//  Thickness function 
//
    TCanvas *c4 = new TCanvas("c4","T_A",400,10,600,700);
    f_ta = new TF1("f_ta", ta, 0, bmax, 0);
    f_ta->SetNpx(nBin);
    f_ta->Draw();
//
//  Kernel of overlap function
//    

    TCanvas *c5 = new TCanvas("c5","T_A_rfti",400,10,600,700);
    f_ta_rfi = new TF2("f_ta_rfi", ta_rfi, 0, bmax, 0., TMath::Pi(), 1);
    f_ta_rfi->SetNpx(nBin);
    f_ta_rfi->SetNpy(nBin);    
    f_ta_rfi->SetParameter(0,0.);
    f_ta_rfi->Draw();
  
//
//  Overlap Function
//
    TCanvas *c6 = new TCanvas("c6","T_AA",400,10,600,700);
    TF1* f_taa = new TF1("f_taa", taa,0.,bmax, 0);
    f_taa->SetNpx(nBin);
    f_taa->Draw();


    TH1F *h = (TH1F*) f_taa->GetHistogram();
    h->SetName("OverlapFunction");
    h->SetTitle("OverlapFunction");
    TH1F *h1 = (TH1F*) f_ta->GetHistogram();
    h1->SetName("ThicknessFunction");
    h1->SetTitle("ThicknessFunction");

    // save
    TFile *f = NULL;
    if (opt == 0) f = new TFile("TPbPb.root","recreate");
    else if (opt == 1) f = new TFile("TAuAu.root","recreate");
    else if (opt == 2) f = new TFile("TXeXe.root","recreate");    
    f->cd();
    h->Write();
    h1->Write();    
    f->Close();

}