// // 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 #include // root headers #include #include #include #include #include // 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(); }