jdbrice
4/22/2015 - 6:57 PM

## Kong's Simultaneous Blast Wave fit code

Kong's Simultaneous Blast Wave fit code

#include <iostream>
#include <fstream>
#include <math.h>
#include "TMath.h"
#include <stdio.h>
#include <iomanip>
#include <vector>
#include "TGraph.h"
#include "TF1.h"
#include "TH1D.h"
#include "TFile.h"
#include "makeMultiPanelCanvas.C"

#include "fitting.h"

using namespace RooFit;

using namespace std;

vector<double> x,ex, y,ey,x1,ex1,y1,ey1,x2,ex2,y2,ey2;

double getChi2;
int getNdf;

const double R = 1.0;
const double dr = 1e-2; // FIXME

double integral(double beta_T, double T, double n, double pt, double mt)
{

double s = 0.;
for(double r = dr/2; r < R; r += dr)
{
double beta_r = (beta_T*(n+2)/2) * TMath::Power((r/R),n);
double rho = TMath::ATanH(beta_r);

s += r * dr * TMath::BesselK1((mt*cosh(rho))/T) * TMath::BesselI0((pt*sinh(rho))/T);
}

return s;
}

void function(int &npar, double *gin, double &f, double *par, int flag)
{
double aka    = par[0];
double aka1   = par[4];
double aka2   = par[5];
double T      = par[2];
double beta_T = par[1];
double n = par[3];

double chi2 = 0.;

for(int i = 0; i < int(x.size()); i++)
{
double m = 0.497;

double pt = x[i];
double mt = sqrt(m*m + pt*pt);

double a = 0.;
a = aka;

double theo = a * mt * integral(beta_T, T, n, pt, mt);
double q = (y[i] - theo)/ey[i];

chi2 += q*q;
}

for(i = 0; i < int(x1.size()); i++){

double m1 = 1.115;

double pt1 = x1[i];
double mt1 = sqrt(m1*m1 + pt1*pt1);

double a1 = 0.;
a1 = aka1;

double theo1 = a1 * mt1 * integral(beta_T,T,n,pt1,mt1);
double q1 = (y1[i]-theo1)/ey1[i];

chi2 += q1*q1;
}

for(i = 0; i < int(x2.size()); i++){

double m2 = 1.314;

double pt2 = x2[i];
double mt2 = sqrt(m2*m2 + pt2*pt2);

double a2 = 0.;
a2 = aka2;

double theo2 = a2 * mt2 * integral(beta_T,T,n,pt2,mt2);
double q2 = (y2[i]-theo2)/ey2[i];

chi2 += q2*q2;
}

f = chi2;

getChi2 = f;
getNdf  = x.size() + x1.size() + x2.size() - 6 - 1;
}

double MyFunc( double *pt, double *p){

double mass = 0.497;

double temp = 0.;
double mt = sqrt(pt[0]*pt[0]+mass*mass);

temp = p[2] * mt * integral(p[0],p[1],p[3],pt[0],mt);

return temp;

}

double MyFunc1( double *pt, double *p){

double mass = 1.115;

double temp = 0.;
double mt = sqrt(pt[0]*pt[0]+mass*mass);

temp = p[2] * mt * integral(p[0],p[1],p[3],pt[0],mt);

return temp;

}

double MyFunc2( double *pt, double *p){

double mass = 1.314;

double temp = 0.;
double mt = sqrt(pt[0]*pt[0]+mass*mass);

temp = p[2] * mt * integral(p[0],p[1],p[3],pt[0],mt);

return temp;

}

double beta_T(double beta_s,double n){

return (beta_s*2)/(n+2);

}

void MeanPtFromMidRapidityCombinedBlastWaveFit(){

gStyle->SetErrorX(0);

TFile* file = new TFile("/Users/kongkong/2015Research/Spectra/pPb2013/rpyDependence/new_files/MidRpy0_0.8_v1.root");

TH1D* ksSpectra[8];
TH1D* laSpectra[8];

stringstream ksHistName;
stringstream laHistName;

for (int mult = 0; mult < 8; mult++){

ksHistName.str("");
laHistName.str("");

ksHistName << "ksSpectra_vtx_";
ksHistName << mult+1;

laHistName << "laSpectra_vtx_";
laHistName << mult+1;

ksSpectra[mult] = (TH1D*)file->Get(ksHistName.str().c_str());
laSpectra[mult] = (TH1D*)file->Get(laHistName.str().c_str());

}

TFile* file1 = new TFile("/Users/kongkong/2015Research/Spectra/pPb2013/XiSpectra/SpectraWithVtxWeight_Combine_V18.root");

TH1D* xiSpectra[8];
TH1D* xiSpectra[8];

xiSpectra[0] = (TH1D*)file1->Get("Spectra_NtrkOffline0_35");
xiSpectra[1] = (TH1D*)file1->Get("Spectra_NtrkOffline35_60");
xiSpectra[2] = (TH1D*)file1->Get("Spectra_NtrkOffline60_90");
xiSpectra[3] = (TH1D*)file1->Get("Spectra_NtrkOffline90_120");
xiSpectra[4] = (TH1D*)file1->Get("Spectra_NtrkOffline120_150");
xiSpectra[5] = (TH1D*)file1->Get("Spectra_NtrkOffline150_185");
xiSpectra[6] = (TH1D*)file1->Get("Spectra_NtrkOffline185_220");
xiSpectra[7] = (TH1D*)file1->Get("Spectra_NtrkOffline220_inf");

double ks_pTbinsBound[34] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22,24,26,28,30,34,38,42,46,50,56,66,90};
double ks_ptbins[34] = {0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
double ks_binwidth[33] = {0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.2,0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.6,1.0,2.4};
double ks_ptbincenter[33] = {0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1.05,1.15,1.25,1.35,1.45,1.55,1.65,1.75,1.85,1.95,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};

double la_pTbinsBound[22] = {4,6,8,10,12,14,16,18,20,22,24,26,28,30,34,38,42,46,50,56,66,90};
double la_ptbins[22] = {0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
double la_ptbincenter[21] = {0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};
double la_binwidth[21] = {0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.6,1.0,2.4};

double xi_pTbinsBound[21] = {6,8,10,12,14,16,18,20,22,24,26,28,30,34,38,42,46,50,56,66,90};
double xi_ptbins[21] = {0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
double xi_ptbincenter[20] = {0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};
double xi_binwidth[20] = {0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.6,1.0,2.4};

TCanvas* c1 = new TCanvas();
c1->Divide(4,2,0,0);

double beta_T[8],Tkin[8],aka[8],aka1[8],aka2[8],n[8];
double ebeta_T[8],eTkin[8],eaka[8],eaka1[8],eaka2[8],en[8];

stringstream f1Name;
stringstream f2Name;
stringstream f3Name;

TF1* f1[8];
TF1* f2[8];
TF1* f3[8];

TH1D* f1_hist[8];
TH1D* f2_hist[8];
TH1D* f3_hist[8];

stringstream f1HistName;
stringstream f2HistName;
stringstream f3HistName;

double ks_ptbins_histo[34] = {0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
double la_ptbins_histo[24] = {0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
double xi_ptbins_histo[24] = {0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};

for(int mult = 0; mult < 8; mult++){

f1HistName.str("");
f2HistName.str("");

f1HistName << "ks_fitFuncHist_";
f1HistName << mult+1;

f2HistName << "la_fitFuncHist_";
f2HistName << mult+1;

f1_hist[mult] = new TH1D(f1HistName.str().c_str(),f1HistName.str().c_str(),33,ks_ptbins_histo);
f2_hist[mult] = new TH1D(f2HistName.str().c_str(),f2HistName.str().c_str(),23,la_ptbins_histo);

}

for(mult = 0; mult < 8; mult++){

f3HistName.str("");

f3HistName << "xi_fitFuncHist_";
f3HistName << mult+1;

f3_hist[mult] = new TH1D(f3HistName.str().c_str(),f3HistName.str().c_str(),23,xi_ptbins_histo);
}

TLatex* ratio[8];
ratio[0] = new TLatex(2.5,45,"2 < N^{offline}_{trk} < 35");
ratio[1] = new TLatex(2.5,45,"35 < N^{offline}_{trk} < 60");
ratio[2] = new TLatex(2.5,45,"60 < N^{offline}_{trk} < 90");
ratio[3] = new TLatex(2.5,45,"90 < N^{offline}_{trk} < 120");
ratio[4] = new TLatex(2.5,45,"120 < N^{offline}_{trk} < 150");
ratio[5] = new TLatex(2.5,45,"150 < N^{offline}_{trk} < 185");
ratio[6] = new TLatex(2.5,45,"185 < N^{offline}_{trk} < 220");
ratio[7] = new TLatex(2.5,45,"220 < N^{offline}_{trk} < #infty");

TLine* l1 = new TLine(0,1,9.0,1.0);
l1->SetLineWidth(2);
l1->SetLineColor(kRed);
l1->SetLineStyle(2);

TLegend *w1 = new TLegend(0.25,0.4,0.5,0.5);
w1->SetLineColor(kWhite);
w1->SetFillColor(0);

TGraph* test[8];

for(int mult = 0; mult < 8; mult++){

x.clear();
ex.clear();
y.clear();
ey.clear();

x1.clear();
ex1.clear();
y1.clear();
ey1.clear();

x2.clear();
ex2.clear();
y2.clear();
ey2.clear();

for (int pt = 3; pt < 18; pt++){

x.push_back( ks_ptbincenter[pt] );
ex.push_back(0.0);
y.push_back( ksSpectra[mult]->GetBinContent(pt+1) );
ey.push_back( (ksSpectra[mult]->GetBinError(pt+1)) );

}

for(pt = 1; pt < 12; pt++){

x1.push_back( la_ptbincenter[pt] );
ex1.push_back(0.0);
y1.push_back( laSpectra[mult]->GetBinContent(pt+2) );
ey1.push_back( (laSpectra[mult]->GetBinError(pt+2)) );

}

for(pt = 0; pt < 11; pt++){

x2.push_back( xi_ptbincenter[pt] );
ex2.push_back(0.0);
y2.push_back( xiSpectra[mult]->GetBinContent(pt+1) );
ey2.push_back( (xiSpectra[mult]->GetBinError(pt+1)) );

}

/**
* simultaneous fit for K0s and Lambda:
*/

TMinuit * gMinuit[8];
gMinuit[mult] = new TMinuit(6);

//set the function to minimize with minuit;
gMinuit[mult]->SetFCN(function);

double arglist[10];
arglist[0] = -1;

gMinuit[mult]->mnexcm("SET ERR", arglist, 1, 0);

gMinuit[mult]->mnparm(0, "aka",    10,    0.1,   1,    10000, 0);
gMinuit[mult]->mnparm(1, "beta_T", 0.32,   0.0001,  0.005, 1.0,   0);
gMinuit[mult]->mnparm(2, "Tkin",   0.17,  0.0001,  0.005, 1.0,   0);
gMinuit[mult]->mnparm(3, "n",      1.0,   0.01,  0.1,  10.0,  0);
gMinuit[mult]->mnparm(4, "aka1",   10,    0.1,   1,    10000, 0);
gMinuit[mult]->mnparm(5, "aka2",   10,    0.1,   1,    10000, 0);

gMinuit[mult]->mnexcm("MINOS",     arglist,  1,   0);
gMinuit[mult]->mnexcm("CALL FCN",  arglist,  1,   0);
//gMinuit[mult]->mnexcm("HESSE",     arglist,  1,   0);

gMinuit[mult]->GetParameter(0, aka[mult],    eaka[mult]);
gMinuit[mult]->GetParameter(1, beta_T[mult], ebeta_T[mult]);
gMinuit[mult]->GetParameter(2, Tkin[mult],   eTkin[mult]);
gMinuit[mult]->GetParameter(3, n[mult],      en[mult]);
gMinuit[mult]->GetParameter(4, aka1[mult],   eaka1[mult]);
gMinuit[mult]->GetParameter(5, aka2[mult],   eaka2[mult]);

/**
* define the fit function by using the parameters from fit:
*/

double xmin = 0.3;
double xmax = 1.8;

double xmin1 = 0.6;
double xmax1 = 3.0;

double xmin2 = 0.6;
double xmax2 = 3.0;

f1Name.str("");
f2Name.str("");
f3Name.str("");

f1Name << "f1_";
f1Name << mult+1;

f2Name << "f2_";
f2Name << mult+1;

f3Name << "f3_";
f3Name << mult+1;

f1[mult] = new TF1(f1Name.str().c_str(),MyFunc,xmin,xmax,4);
f1[mult]->SetParameter(0,beta_T[mult]);
f1[mult]->SetParameter(1,Tkin[mult]);
f1[mult]->SetParameter(2,aka[mult]);
f1[mult]->SetParameter(3,n[mult]);

f2[mult] = new TF1(f2Name.str().c_str(),MyFunc1,xmin1,xmax1,4);
f2[mult]->SetParameter(0,beta_T[mult]);
f2[mult]->SetParameter(1,Tkin[mult]);
f2[mult]->SetParameter(2,aka1[mult]);
f2[mult]->SetParameter(3,n[mult]);

f3[mult] = new TF1(f3Name.str().c_str(),MyFunc2,xmin2,xmax2,4);
f3[mult]->SetParameter(0,beta_T[mult]);
f3[mult]->SetParameter(1,Tkin[mult]);
f3[mult]->SetParameter(2,aka2[mult]);
f3[mult]->SetParameter(3,n[mult]);

c1->cd(mult+1);

ksSpectra[mult]->SetMarkerSize(1.3);
ksSpectra[mult]->SetMarkerStyle(20);
ksSpectra[mult]->SetMarkerColor(kBlue);
ksSpectra[mult]->SetStats(kFALSE);
ksSpectra[mult]->SetTitle("");
ksSpectra[mult]->SetXTitle("P^{}_{T,V0}(GeV/c)");
ksSpectra[mult]->SetYTitle("1/N^{}_{ev}1/(2#PiP^{}_{T}d^{2}N/(dP^{}_{T}dy) [(GeV/c)^{-2}]");

ksSpectra[mult]->GetYaxis()->SetRangeUser(0.0000001,1000);

laSpectra[mult]->SetMarkerSize(1.3);
laSpectra[mult]->SetMarkerStyle(20);
laSpectra[mult]->SetMarkerColor(kBlack);

xiSpectra[mult]->SetMarkerSize(1.3);
xiSpectra[mult]->SetMarkerStyle(20);
xiSpectra[mult]->SetMarkerColor(kYellow-2);

ksSpectra[mult]->Draw("P");
laSpectra[mult]->Draw("Psame");
xiSpectra[mult]->Draw("Psame");
ratio[mult]->Draw("same");

f2[mult]->Draw("same");
f1[mult]->Draw("same");
f3[mult]->Draw("same");

double la_ptbincenter_extrapolate[23] = {0.1,0.3,0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};
double ks_ptbincenter_extrapolate[33] = {0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1.05,1.15,1.25,1.35,1.45,1.55,1.65,1.75,1.85,1.95,
2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};

double xi_ptbincenter_extrapolate[23] = {0.1,0.3,0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};

for(int pt = 0; pt < 2; pt++){

double mass = 1.115;

double temp = 0.;
double mt = sqrt(la_ptbincenter_extrapolate[pt]*la_ptbincenter_extrapolate[pt]+mass*mass);

temp = aka1[mult] * mt * integral(beta_T[mult],Tkin[mult],n[mult],la_ptbincenter_extrapolate[pt],mt);

f2_hist[mult]->SetBinContent(pt+1, temp);
f2_hist[mult]->SetBinError(pt+1, 0.01*temp);

}

for(int pt = 0; pt < 3; pt++){

double mass = 1.314;
double temp = 0.;
double mt = sqrt(xi_ptbincenter_extrapolate[pt]*xi_ptbincenter_extrapolate[pt]+mass*mass);

temp = aka2[mult] * mt * integral(beta_T[mult],Tkin[mult],n[mult],xi_ptbincenter_extrapolate[pt],mt);

f3_hist[mult]->SetBinContent(pt+1, temp);
f3_hist[mult]->SetBinError(pt+1, 0.01*temp);

}

}

for(mult = 0; mult < 8; mult++){

for(pt = 0; pt < 33; pt++){

double temp = ksSpectra[mult]->GetBinContent(pt+1);
double temp_err = ksSpectra[mult]->GetBinError(pt+1);

f1_hist[mult]->SetBinContent(pt+1, temp);
f1_hist[mult]->SetBinError(pt+1, temp_err);

}

for(pt = 0; pt < 21; pt++){

double temp = laSpectra[mult]->GetBinContent(pt+2);
double temp_err = laSpectra[mult]->GetBinError(pt+2);

f2_hist[mult]->SetBinContent(pt+3, temp);
f2_hist[mult]->SetBinError(pt+3, temp_err);

}

for(pt = 0; pt < 20; pt++){

double temp = xiSpectra[mult]->GetBinContent(pt+1);
double temp_err = xiSpectra[mult]->GetBinError(pt+1);

f3_hist[mult]->SetBinContent(pt+4, temp);
f3_hist[mult]->SetBinError(pt+4, temp_err);

}

}

TFile savefile("../files/meanPtSpectra_3particles.root","new");

for(mult = 0; mult < 8; mult++){

f1_hist[mult]->Write();
f2_hist[mult]->Write();
f3_hist[mult]->Write();

}

}