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);
w1->AddEntry(ksSpectra[7],"K^{0}_{s}","P");
w1->AddEntry(laSpectra[7],"#Lambda/#bar{#Lambda}","P");
w1->AddEntry(xiSpectra[7],"#Xi^{0}","P");
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("MIGRAD", arglist, 1, 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);
gPad->SetLogy(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();
}
}