Kong's Simultaneous Blast Wave fit code

Kong_Sim_Fit.cpp

#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();

  }



}