From cbf4f1cb440def23bedbcd39ef7492555ae03995 Mon Sep 17 00:00:00 2001 From: dgangadh Date: Sun, 28 Sep 2014 16:13:24 +0200 Subject: [PATCH] Macro updates --- PWGCF/FEMTOSCOPY/macros/Fit_c3.C | 952 ++++++++++++++++++ PWGCF/FEMTOSCOPY/macros/Makec3EAfile.C | 134 +++ PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C | 461 ++++++--- .../FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C | 232 +++-- PWGCF/FEMTOSCOPY/macros/Plot_plotsTPR.C | 688 ++++++++++--- 5 files changed, 2131 insertions(+), 336 deletions(-) create mode 100755 PWGCF/FEMTOSCOPY/macros/Fit_c3.C create mode 100755 PWGCF/FEMTOSCOPY/macros/Makec3EAfile.C diff --git a/PWGCF/FEMTOSCOPY/macros/Fit_c3.C b/PWGCF/FEMTOSCOPY/macros/Fit_c3.C new file mode 100755 index 00000000000..04d65a079ff --- /dev/null +++ b/PWGCF/FEMTOSCOPY/macros/Fit_c3.C @@ -0,0 +1,952 @@ +#include +#include +#include +#include +#include +#include + +#include "TVector2.h" +#include "TFile.h" +#include "TString.h" +#include "TF1.h" +#include "TF3.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TProfile3D.h" +#include "TMath.h" +#include "TText.h" +#include "TRandom3.h" +#include "TArray.h" +#include "TLegend.h" +#include "TStyle.h" +#include "TMinuit.h" +#include "TCanvas.h" +#include "TLatex.h" +#include "TPaveStats.h" +#include "TASImage.h" +#include "TGraph.h" +#include "TSpline.h" +#include "TVirtualFitter.h" + +#define BohrR 1963.6885 // Bohr Radius for pions +#define FmToGeV 0.19733 // conversion to fm +#define PI 3.1415926 +#define masspiC 0.1395702 // pi+ mass (GeV/c^2) + +using namespace std; + +// +int CollisionType_def=1;// PbPb, pPb, pp +int FitType=0;// (0)Edgeworth, (1)Laguerre, (2)Levy +// +int Mbin=0;// 0-9: centrality bin in widths of 5% +int CHARGE=-1;// -1 or +1: + or - pions for same-charge case, --+ or -++, ---+ or -+++ +// +int EDbin=0;// 0 or 1: Kt3 bin +double G_def = 90.;// coherent % +// +bool MRC=1;// Momentum Resolution Corrections? +bool MuonCorrection=1;// correct for Muons? +// +int f_choice=0;// 0(Core/Halo), 1(40fm), 2(70fm), 3(100fm) +// +// +// +// +bool SaveToFile_def=0; +int fFSIindex=0; +float TwoFrac;// Lambda parameter +float OneFrac;// Lambda^{1/2} +float ThreeFrac;// Lambda^{3/2} +double Qcut_pp = 0.6;// 0.6 +double Qcut_PbPb = 0.1; +double NormQcutLow_pp = 1.0; +double NormQcutHigh_pp = 1.5; +double NormQcutLow_PbPb = .15; +double NormQcutHigh_PbPb = .2;// was .175 + + +const int BINRANGE_3=60;// q12,q13,q23 +int BINLIMIT_3; +double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 +double A_3_e[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 +double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23 +double BinCenters[400]; +double Chi2_c3; +double NFitPoints_c3; +double Q3LimitLow; +double Q3LimitHigh; +void fcn_c3(int&, double*, double&, double[], int); + +int RbinMRC; + +TH1D *fFSIss[12]; +TH1D *fFSIos[12]; + +void SetFSICorrelations(); +void SetFSIindex(Float_t); +Float_t FSICorrelation(Int_t, Int_t, Float_t); +void SetMuonCorrections(); +void SetMomResCorrections(); +double Gamov(int, double); + +// +float fact(float); + + + +TH1D *MRC_SC_3[3]; +TH1D *C3muonCorrectionSC[2]; + +double AvgQinvSS[30]; +double AvgQinvOS[30]; +double BinCentersQ4[20]; + + + +void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, double G=G_def){ + + SaveToFile_def=SaveToFile; + CollisionType_def=CollisionType; + G /= 100.; + G_def=G; + + TFile *_file0; + if(CollisionType==0){// PbPb + _file0 = new TFile("Results/PDC_11h_3Dhistos.root","READ"); + }else if(CollisionType==1){// pPb + _file0 = new TFile("Results/PDC_13bc_kINT7.root","READ"); + }else{// pp + _file0 = new TFile("Results/PDC_10bcde_kMB.root","READ"); + } + + TList *MyList; + TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root"); + MyList=(TList*)tdir->Get("FourPionOutput_1"); + //MyList=(TList*)_file0->Get("MyList"); + + _file0->Close(); + + + if(CollisionType==0) {Q3LimitLow = 0.02; Q3LimitHigh = 0.06;}// 0.01 and 0.08 + else {Q3LimitLow = 0.08; Q3LimitHigh = 0.2;}// 0.01 and 0.25 + + + // + TwoFrac=0.7; + OneFrac = sqrt(TwoFrac); + ThreeFrac = pow(TwoFrac, 1.5); + + //// Core/Halo, 40fm, 70fm, 100fm + float ThermShift_f33[4]={0., 0.06933, 0.01637, 0.006326}; + float ThermShift_f32[4]={0., -0.0185, -0.005301, -0.002286}; + float ThermShift_f31[4]={0., -0.01382, -0.0004682, 0.0005337}; + float f33prime = ThreeFrac; + float f32prime = TwoFrac*(1-OneFrac); + float f31prime = pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2); + f33prime += ThermShift_f33[f_choice]; + f32prime += ThermShift_f32[f_choice]; + f31prime += ThermShift_f31[f_choice]; + float f33 = f33prime; + float f32 = f32prime/TwoFrac; + float f31 = f31prime - 3*((1-TwoFrac)*(1-OneFrac) + ThermShift_f32[f_choice]*(1-TwoFrac)/TwoFrac); + //cout<Append("_Charge2_"); + *name3 += 0; + name3->Append("_Charge3_"); + *name3 += 0; + name3->Append("_M_"); + *name3 += Mbin; + name3->Append("_ED_"); + *name3 += EDbin; + name3->Append("_Term_"); + *name3 += term+1; + + TString *nameNorm3=new TString(name3->Data()); + nameNorm3->Append("_Norm"); + // + TString *nameK3=new TString(name3->Data()); + nameK3->Append("_Kfactor3D"); + // + name3->Append("_3D"); + + + + norm_3[term] = ((TH1D*)MyList->FindObject(nameNorm3->Data()))->Integral(); + ThreeParticle[0][0][0][term] = (TH3D*)MyList->FindObject(name3->Data()); + ThreeParticle[0][0][0][term]->Sumw2(); + //if(0==0 && 0==0) cout<<"3-pion norms "<Scale(norm_3[0]/norm_3[term]); + ThreeParticle[0][0][0][term]->SetMarkerStyle(20); + ThreeParticle[0][0][0][term]->SetTitle(""); + // + + // + if(term<4){ + K3avg[0][0][0][term] = (TProfile3D*)MyList->FindObject(nameK3->Data()); + } + + // + }// term + + + + + cout<<"Done getting Histograms"<SetLineStyle(2); + + + int ch1=0,ch2=0,ch3=0; + + + + /////////////////////////////////////////////////////////// + // 3-pion + cout<<"3-pion section"<SetHighLightColor(2); + can2->Range(-0.7838115,-1.033258,0.7838115,1.033258); + gStyle->SetOptFit(0111); + can2->SetFillColor(10); + can2->SetBorderMode(0); + can2->SetBorderSize(2); + can2->SetGridx(); + can2->SetGridy(); + can2->SetFrameFillColor(0); + can2->SetFrameBorderMode(0); + can2->SetFrameBorderMode(0); + gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02); + + int Q3BINS=12; + float Q3HistoLimit=0.12; + if(CollisionType!=0){ Q3BINS=60; Q3HistoLimit=0.6;} + + TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,Q3HistoLimit); + TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,Q3HistoLimit); + TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",Q3BINS,0,Q3HistoLimit); + TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",Q3BINS,0,Q3HistoLimit); + double c3_e[100]={0}; + + + double value_num; + double value_num_e; + double N3_QS; + double N3_QS_e; + + for(int i=2; i<=ThreeParticle[0][0][0][0]->GetNbinsX(); i++){// bin number + double Q12 = BinCenters[i-1];// true center + //int Q12bin = int(Q12/0.01)+1; + // + for(int j=2; j<=ThreeParticle[0][0][0][0]->GetNbinsY(); j++){// bin number + double Q13 = BinCenters[j-1];// true center + //int Q13bin = int(Q13/0.01)+1; + // + for(int k=2; k<=ThreeParticle[0][0][0][0]->GetNbinsZ(); k++){// bin number + double Q23 = BinCenters[k-1];// true center + //int Q23bin = int(Q23/0.01)+1; + // + + if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible + if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible + + double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2)); + int Q3bin = c3_hist->GetXaxis()->FindBin(Q3); + + // + double K3 = K3avg[0][0][0][0]->GetBinContent(i,j,k); + double K2_12 = K3avg[0][0][0][1]->GetBinContent(i,j,k); + double K2_13 = K3avg[0][0][0][2]->GetBinContent(i,j,k); + double K2_23 = K3avg[0][0][0][3]->GetBinContent(i,j,k); + + + if(K3==0) continue; + + double TERM1=ThreeParticle[ch1][ch2][ch3][0]->GetBinContent(i,j,k);// N3 (3-pion yield per q12,q13,q23 cell). 3-pions from same-event + double TERM2=ThreeParticle[ch1][ch2][ch3][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event + double TERM3=ThreeParticle[ch1][ch2][ch3][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event + double TERM4=ThreeParticle[ch1][ch2][ch3][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event + double TERM5=ThreeParticle[ch1][ch2][ch3][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics) + + + if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue; + if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue; + + // + if(MRC){ + TERM1 *= MRC_SC_3[0]->GetBinContent(MRC_SC_3[0]->GetXaxis()->FindBin(Q3)); + TERM2 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3)); + TERM3 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3)); + TERM4 *= MRC_SC_3[1]->GetBinContent(MRC_SC_3[1]->GetXaxis()->FindBin(Q3)); + TERM5 *= MRC_SC_3[2]->GetBinContent(MRC_SC_3[2]->GetXaxis()->FindBin(Q3)); + } + double MuonCorr1=1, MuonCorr2=1, MuonCorr3=1, MuonCorr4=1; + if(MuonCorrection){ + MuonCorr1 = C3muonCorrectionSC[0]->GetBinContent(C3muonCorrectionSC[0]->GetXaxis()->FindBin(Q3)); + MuonCorr2 = C3muonCorrectionSC[1]->GetBinContent(C3muonCorrectionSC[0]->GetXaxis()->FindBin(Q3)); + MuonCorr3 = MuonCorr2; + MuonCorr4 = MuonCorr2; + } + + + + // Purify. Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations) + N3_QS = TERM1; + N3_QS -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5; + N3_QS -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5); + N3_QS /= ThreeFrac; + N3_QS *= K3; + N3_QS *= MuonCorr1; + + + // Isolate 3-pion cumulant + value_num = N3_QS; + value_num -= (TERM2 - (1-TwoFrac)*TERM5)*K2_12/TwoFrac * MuonCorr2; + value_num -= (TERM3 - (1-TwoFrac)*TERM5)*K2_13/TwoFrac * MuonCorr3; + value_num -= (TERM4 - (1-TwoFrac)*TERM5)*K2_23/TwoFrac * MuonCorr4; + value_num += 2*TERM5; + + + + + // errors + N3_QS_e = TERM1; + N3_QS_e += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2); + N3_QS_e += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2)); + N3_QS_e /= pow(ThreeFrac,2); + N3_QS_e *= pow(K3,2); + // + value_num_e = N3_QS_e; + value_num_e += (pow(K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)*K2_12/TwoFrac*sqrt(TERM5),2)); + value_num_e += (pow(K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)*K2_13/TwoFrac*sqrt(TERM5),2)); + value_num_e += (pow(K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)*K2_23/TwoFrac*sqrt(TERM5),2)); + value_num_e += pow(2*sqrt(TERM5),2); + + c3_e[Q3bin-1] += value_num_e + TERM5;// add baseline stat error + + + // Fill histograms + c3_hist->Fill(Q3, value_num + TERM5);// for cumulant correlation function + Combinatorics_1d->Fill(Q3, TERM5); + + // + A_3[i-1][j-1][k-1] = value_num + TERM5; + B_3[i-1][j-1][k-1] = TERM5; + A_3_e[i-1][j-1][k-1] = sqrt(value_num_e + TERM5); + //if(i==j && i==k && i==4) cout<SetBinError(i+1, sqrt(c3_e[i]));} + + c3_hist->Divide(Combinatorics_1d); + + /////////////////////////////////////////////////////////////////////////////////////////////////// + + + + c3_hist->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})"); + c3_hist->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}"); + c3_hist->GetYaxis()->SetTitleSize(0.045); c3_hist->GetXaxis()->SetTitleSize(0.045); + c3_hist->GetYaxis()->SetTitleOffset(1.1); + c3_hist->GetXaxis()->SetRangeUser(0,Q3LimitHigh); + c3_hist->GetYaxis()->SetRangeUser(0.9,4); + c3_hist->SetMarkerStyle(25); + c3_hist->SetMarkerColor(2); + c3_hist->SetLineColor(2); + c3_hist->SetMaximum(3); + c3_hist->SetMinimum(.7); + c3_hist->Draw(); + //legend2->AddEntry(c3_hist,"#font[12]{#bf{c}}_{3} (cumulant correlation)","p"); + + + + const int npar_c3=7; + TMinuit MyMinuit_c3(npar_c3); + MyMinuit_c3.SetFCN(fcn_c3); + int ierflg_c3=0; + double arglist_c3 = 0; + MyMinuit_c3.mnexcm("SET NOWarnings",&arglist_c3,1, ierflg_c3); + arglist_c3 = -1; + MyMinuit_c3.mnexcm("SET PRint",&arglist_c3,1, ierflg_c3); + //arglist_c3=2;// improve Minimization Strategy (1 is default) + //MyMinuit_c3.mnexcm("SET STR",&arglist_c3,1,ierflg_c3); + //arglist_c3 = 0; + //MyMinuit_c3.mnexcm("SCAN", &arglist_c3,1,ierflg_c3); + arglist_c3 = 5000; + MyMinuit_c3.mnexcm("MIGRAD", &arglist_c3 ,1,ierflg_c3); + + + TF1 *c3Fit1D_Expan; + if(FitType==0) { + c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2.,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2.,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)) + [5]/(120.*pow(2.,2.5))*(32.*pow(x/sqrt(3.)*[2]/0.19733,5) - 160.*pow(x/sqrt(3.)*[2]/0.19733,3) + 120*x/sqrt(3.)*[2]/0.19733) ,3))",0,1); + }else if(FitType==1){ + c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-[2]*x/0.19733 * sqrt(3.)/2.) * pow(1 + [3]*([2]*x/sqrt(3.)/0.19733 - 1) + [4]/2*(pow([2]*x/sqrt(3.)/0.19733,2) - 4*[2]*x/sqrt(3.)/0.19733 + 2) + [5]/6.*(-pow(x/sqrt(3.)*[1]/0.19733,3) + 9*pow(x/sqrt(3.)*[1]/0.19733,2) - 18*x/sqrt(3.)*[1]/0.19733 + 6),3))",0,1); + }else{ + c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733, [3])))",0,1); + } + + + double OutputPar_c3[npar_c3]={0}; + double OutputPar_c3_e[npar_c3]={0}; + + double par_c3[npar_c3]; // the start values + double stepSize_c3[npar_c3]; // step sizes + double minVal_c3[npar_c3]; // minimum bound on parameter + double maxVal_c3[npar_c3]; // maximum bound on parameter + string parName_c3[npar_c3]; + // 1.0 1.5 10. 0. 0. 0. 1.5 + par_c3[0] = 1.0; par_c3[1] = 1.5; par_c3[2] = 10.; par_c3[3] = 0.; par_c3[4] = 0.; par_c3[5] = 0; par_c3[6] = 1.5; + stepSize_c3[0] = 0.01; stepSize_c3[1] = 0.1; stepSize_c3[2] = 0.1; stepSize_c3[3] = 0.01; stepSize_c3[4] = 0.01; stepSize_c3[5] = 0.01; stepSize_c3[6] = 0.1; + minVal_c3[0] = 0.8; minVal_c3[1] = 0.4; minVal_c3[2] = 4.; minVal_c3[3] = -2; minVal_c3[4] = -1; minVal_c3[5] = -1; minVal_c3[6] = 0.5; + maxVal_c3[0] = 1.1; maxVal_c3[1] = 2000.; maxVal_c3[2] = 100.; maxVal_c3[3] = +2; maxVal_c3[4] = +1; maxVal_c3[5] = +1; maxVal_c3[6] = 2.5; + parName_c3[0] = "N"; parName_c3[1] = "#lambda_{3}"; parName_c3[2] = "R_{inv}"; parName_c3[3] = "coeff_{3}"; parName_c3[4] = "coeff_{4}"; parName_c3[5] = "coeff_{5}"; parName_c3[6] = "#alpha"; + + if(CollisionType==0){ + if(FitType!=0) { + par_c3[2]=15.; minVal_c3[2] = 8.; + } + }else{ + if(FitType==0) {par_c3[2] = 2.0; minVal_c3[2] = 1.0;} + else { + par_c3[1] = 4.0; minVal_c3[1] = 1.0; + // + par_c3[2] = 1.3; minVal_c3[2] = .9; maxVal_c3[2] = 10.; + } + } + + if(FitType==0) {par_c3[6] = 2.;} + if(FitType==1) {par_c3[6] = 1.;} + if(FitType==2) {par_c3[3] = 0; par_c3[4] = 0; par_c3[5] = 0;} + + if(FitType==2) {maxVal_c3[1] = 10.;} + + for (int i=0; i sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible + + double TERM5=ThreeParticle[ch1][ch2][ch3][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics) + if(TERM5==0) continue; + + + if(MRC) TERM5 *= MRC_SC_3[2]->GetBinContent(MRC_SC_3[2]->GetXaxis()->FindBin(Q3)); + // + double radius = OutputPar_c3[2]/FmToGeV; + double arg12 = Q12*radius; + double arg13 = Q13*radius; + double arg23 = Q23*radius; + double Expan12=1, Expan13=1, Expan23=1; + if(FitType==0){ + Expan12 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg12,3) - 12*pow(arg12,1)); + Expan12 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg12,4) -48*pow(arg12,2) + 12); + Expan12 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg12,5) - 160.*pow(arg12,3) + 120*arg12); + // + Expan13 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg13,3) - 12*pow(arg13,1)); + Expan13 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg13,4) -48*pow(arg13,2) + 12); + Expan13 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg13,5) - 160.*pow(arg13,3) + 120*arg13); + // + Expan23 += OutputPar_c3[3]/pow(2.,3/2.)/(6.)*(8*pow(arg23,3) - 12*pow(arg23,1)); + Expan23 += OutputPar_c3[4]/pow(2.,4/2.)/(24.)*(16*pow(arg23,4) -48*pow(arg23,2) + 12); + Expan23 += OutputPar_c3[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg23,5) - 160.*pow(arg23,3) + 120*arg23); + }else if(FitType==1){ + Expan12 += OutputPar_c3[3]*(arg12 - 1); + Expan12 += OutputPar_c3[4]/2.*(pow(arg12,2) - 4*arg12 + 2); + Expan12 += OutputPar_c3[5]/6.*(-pow(arg12,3) + 9*pow(arg12,2) - 18*arg12 + 6); + // + Expan13 += OutputPar_c3[3]*(arg13 - 1); + Expan13 += OutputPar_c3[4]/2.*(pow(arg13,2) - 4*arg13 + 2); + Expan13 += OutputPar_c3[5]/6.*(-pow(arg13,3) + 9*pow(arg13,2) - 18*arg13 + 6); + // + Expan23 += OutputPar_c3[3]*(arg23 - 1); + Expan23 += OutputPar_c3[4]/2.*(pow(arg23,2) - 4*arg23 + 2); + Expan23 += OutputPar_c3[5]/6.*(-pow(arg23,3) + 9*pow(arg23,2) - 18*arg23 + 6); + }else{} + + // + double C = OutputPar_c3[1] * pow(1-G,3)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg13,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan12*Expan13*Expan23; + C += pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg13,OutputPar_c3[6]))/2.)*Expan12*Expan13; + C += pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg12,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan12*Expan23; + C += pow(OutputPar_c3[1], 2/3.) * G*pow(1-G,2)*exp(-(pow(arg13,OutputPar_c3[6])+pow(arg23,OutputPar_c3[6]))/2.)*Expan13*Expan23; + C += 1.0; + //if(Q3<0.026) cout<Fill(Q3, C); + GenSignalExpected_den->Fill(Q3, TERM5); + //if(Q3<0.02) cout<Sumw2(); + GenSignalExpected_num->Divide(GenSignalExpected_den); + + TSpline3 *c3Fit1D_ExpanSpline = new TSpline3(GenSignalExpected_num); + c3Fit1D_ExpanSpline->SetLineWidth(2); + double xpoints[1000], ypoints[1000]; + bool splineOnly=kFALSE; + for(int iii=0; iii<1000; iii++){ + xpoints[iii] = 0 + (iii+0.5)*0.001; + //ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]);// to skip spline + splineOnly=kTRUE;// to skip 1D approximation + if(CollisionType==0) splineOnly=kTRUE; + if(CollisionType!=0 && xpoints[iii] > 0.06) splineOnly=kTRUE; + if(!splineOnly){ + ypoints[iii] = c3Fit1D_Expan->Eval(xpoints[iii]); + if(c3Fit1D_Expan->Eval(xpoints[iii])<2. && fabs(c3Fit1D_Expan->Eval(xpoints[iii])-c3Fit1D_ExpanSpline->Eval(xpoints[iii])) < 0.01) splineOnly=kTRUE; + } + else {ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]); splineOnly=kTRUE;} + } + TGraph *gr_c3Spline = new TGraph(1000,xpoints,ypoints); + gr_c3Spline->SetLineWidth(2); + if(CollisionType==0) gr_c3Spline->SetLineColor(1); + if(CollisionType==1) gr_c3Spline->SetLineColor(2); + if(CollisionType==2) gr_c3Spline->SetLineColor(4); + + gr_c3Spline->Draw("c same"); + + /* + double ChiSqSum_1D=0, SumNDF_1D=0; + for(int bin=1; bin<=300; bin++){ + double binCenter = c3_hist->GetXaxis()->GetBinCenter(bin); + if(binCenter > Q3Limit) continue; + if(c3_hist->GetBinError(bin)==0) continue; + if(binCenter < 0.01) continue; + int grIndex=1; + for(int gr=0; gr<999; gr++){ + if(binCenter > xpoints[gr] && (binCenter < xpoints[gr+1])) {grIndex=gr; break;} + } + + ChiSqSum_1D += pow((c3_hist->GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2); + //cout<GetBinContent(bin)<<" "<GetBinError(bin)<GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2)<Append("_FT"); + *savefilename += FitType; + savefilename->Append("_G"); + *savefilename += int((G+0.001)/0.02); + savefilename->Append(".root"); + TFile *savefile = new TFile(savefilename->Data(),"RECREATE"); + MyMinuit_c3.Write("MyMinuit_c3"); + // + savefile->Close(); + } + + +} + +//________________________________________________________________________ +void SetFSICorrelations(){ + // read in 2-particle and 3-particle FSI correlations = K2 & K3 + // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair + TFile *fsifile = new TFile("KFile.root","READ"); + if(!fsifile->IsOpen()) { + cout<<"No FSI file found!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<Get(nameK2SS->Data()); + // + TString *nameK2OS = new TString("K2os_"); + *nameK2OS += MB; + temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data()); + // + fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone(); + fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone(); + fFSIss[MB]->SetDirectory(0); + fFSIos[MB]->SetDirectory(0); + } + // + + fsifile->Close(); + +} + +double Gamov(int chargeProduct, double qinv){ + + double arg = chargeProduct*2.*PI/(BohrR*qinv/2.); + + return arg/(exp(arg)-1); +} + +void SetMomResCorrections(){ + + TString *momresfilename = new TString("MomResFile"); + if(CollisionType_def!=0) momresfilename->Append("_ppAndpPb"); + momresfilename->Append(".root"); + + TFile *MomResFile = new TFile(momresfilename->Data(),"READ"); + + TString *proName[28]; + for(int ii=0; ii<28; ii++){ + proName[ii] = new TString("MRC_pro_"); + *proName[ii] += ii; + } + + + // + MRC_SC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term1"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC)); + MRC_SC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term2"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC)); + MRC_SC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term3"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC)); + MRC_SC_3[0]->SetDirectory(0); MRC_SC_3[1]->SetDirectory(0); MRC_SC_3[2]->SetDirectory(0); + // + + if(!MRC){ + for(int bin=1; bin<=MRC_SC_3[0]->GetNbinsX(); bin++){ + MRC_SC_3[0]->SetBinContent(bin, 1.0); MRC_SC_3[1]->SetBinContent(bin, 1.0); MRC_SC_3[2]->SetBinContent(bin, 1.0); + } + + } + MomResFile->Close(); + +} + + +float fact(float n){ + return (n < 1.00001) ? 1 : fact(n - 1) * n; +} +//________________________________________________________________________________________ +void SetMuonCorrections(){ + TString *name = new TString("MuonCorrection"); + if(CollisionType_def!=0) name->Append("_ppAndpPb"); + + name->Append(".root"); + TFile *MuonFile=new TFile(name->Data(),"READ"); + TString *proName[22]; + for(int ii=0; ii<22; ii++){ + proName[ii] = new TString("MuonCorr_pro_"); + *proName[ii] += ii; + } + + // + C3muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term1"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC)); + C3muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term2"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC)); + + C3muonCorrectionSC[0]->SetDirectory(0); C3muonCorrectionSC[1]->SetDirectory(0); + + // + if(!MuonCorrection){ + for(int bin=1; bin<=C3muonCorrectionSC[0]->GetNbinsX(); bin++){ + C3muonCorrectionSC[0]->SetBinContent(bin, 1.0); C3muonCorrectionSC[1]->SetBinContent(bin, 1.0); + } + } + + MuonFile->Close(); +} +//________________________________________________________________________ +void SetFSIindex(Float_t R){ + if(CollisionType_def==0){ + if(Mbin==0) fFSIindex = 0;//0-5% + else if(Mbin==1) fFSIindex = 1;//5-10% + else if(Mbin<=3) fFSIindex = 2;//10-20% + else if(Mbin<=5) fFSIindex = 3;//20-30% + else if(Mbin<=7) fFSIindex = 4;//30-40% + else if(Mbin<=9) fFSIindex = 5;//40-50% + else if(Mbin<=12) fFSIindex = 6;//40-50% + else if(Mbin<=15) fFSIindex = 7;//40-50% + else if(Mbin<=18) fFSIindex = 8;//40-50% + else fFSIindex = 8;//90-100% + }else fFSIindex = 9;// pp and pPb + +} +//________________________________________________________________________ +Float_t FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){ + // returns 2-particle Coulomb correlations = K2 + Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.); + Int_t qbinH = qbinL+1; + if(qbinL <= 0) return 1.0; + if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) {// Scaled Gamov approximation + int chargeproduct = 1; + if(charge1!=charge2) { + chargeproduct = -1; + Float_t ScaleFac = (fFSIos[fFSIindex]->GetBinContent(fFSIos[fFSIindex]->GetNbinsX()-1) - 1); + ScaleFac /= (Gamov(chargeproduct, fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(fFSIos[fFSIindex]->GetNbinsX()-1)) - 1); + return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1); + }else{ + Float_t ScaleFac = (fFSIss[fFSIindex]->GetBinContent(fFSIss[fFSIindex]->GetNbinsX()-1) - 1); + ScaleFac /= (Gamov(chargeproduct, fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(fFSIss[fFSIindex]->GetNbinsX()-1)) - 1); + return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1); + } + } + + Float_t slope=0; + if(charge1==charge2){ + slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH); + slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH); + return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL)); + }else { + slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH); + slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH); + return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL)); + } +} +//__________________________________________________________________________ +void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){ + + double q12=0, q13=0, q23=0; + double Expan12=0, Expan13=0, Expan23=0; + double C=0; + double Rch=par[2]/FmToGeV; + double SumChi2=0; + //double lnL=0, term1=0, term2=0; + NFitPoints_c3=0; + //double SumChi2_test=0; + + for(int i=0; i<=BINLIMIT_3; i++){// q12 + for(int j=0; j<=BINLIMIT_3; j++){// q13 + for(int k=0; k<=BINLIMIT_3; k++){// q23 + + if(B_3[i][j][k] == 0) continue; + if(A_3[i][j][k] == 0) continue; + if(A_3_e[i][j][k] == 0) continue; + + q12 = BinCenters[i]; + q13 = BinCenters[j]; + q23 = BinCenters[k]; + double q3 = sqrt(pow(q12,2)+pow(q13,2)+pow(q23,2)); + if(q3 > Q3LimitHigh) continue; + if(q3 < Q3LimitLow) continue; + // + double arg12 = q12*Rch; + double arg13 = q13*Rch; + double arg23 = q23*Rch; + if(FitType==0){// Edgeworth expansion + Expan12 = 1; + Expan12 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg12,3) - 12*pow(arg12,1)); + Expan12 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg12,4) -48*pow(arg12,2) + 12); + Expan12 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg12,5) - 160.*pow(arg12,3) + 120*arg12); + // + Expan13 = 1; + Expan13 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg13,3) - 12*pow(arg13,1)); + Expan13 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg13,4) -48*pow(arg13,2) + 12); + Expan13 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg13,5) - 160.*pow(arg13,3) + 120*arg13); + // + Expan23 = 1; + Expan23 += par[3]/pow(2.,3/2.)/(6.)*(8*pow(arg23,3) - 12*pow(arg23,1)); + Expan23 += par[4]/pow(2.,4/2.)/(24.)*(16*pow(arg23,4) -48*pow(arg23,2) + 12); + Expan23 += par[5]/pow(2.,5/2.)/(120.)*(32.*pow(arg23,5) - 160.*pow(arg23,3) + 120*arg23); + }else if(FitType==1){// Laguerre expansion + Expan12 = 1; + Expan12 += par[3]*(arg12 - 1); + Expan12 += par[4]/2.*(pow(arg12,2) - 4*arg12 + 2); + Expan12 += par[5]/6.*(-pow(arg12,3) + 9*pow(arg12,2) - 18*arg12 + 6); + // + Expan13 = 1; + Expan13 += par[3]*(arg13 - 1); + Expan13 += par[4]/2.*(pow(arg13,2) - 4*arg13 + 2); + Expan13 += par[5]/6.*(-pow(arg13,3) + 9*pow(arg13,2) - 18*arg13 + 6); + // + Expan23 = 1; + Expan23 += par[3]*(arg23 - 1); + Expan23 += par[4]/2.*(pow(arg23,2) - 4*arg23 + 2); + Expan23 += par[5]/6.*(-pow(arg23,3) + 9*pow(arg23,2) - 18*arg23 + 6); + }else{Expan12=1.0; Expan13=1.0; Expan23=1.0;} + // + + C = par[1] * pow(1-G_def,3)*exp(-(pow(arg12,par[6])+pow(arg13,par[6])+pow(arg23,par[6]))/2.)*Expan12*Expan13*Expan23; + C += pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg13,par[6]))/2.)*Expan12*Expan13; + C += pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg23,par[6]))/2.)*Expan12*Expan23; + C += pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg13,par[6])+pow(arg23,par[6]))/2.)*Expan13*Expan23; + C += 1.0; + C *= par[0];// Norm + + + double error = pow(A_3_e[i][j][k]/B_3[i][j][k],2); + error += pow(sqrt(B_3[i][j][k])*A_3[i][j][k]/pow(B_3[i][j][k],2),2); + error = sqrt(error); + SumChi2 += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2); + + //if(q3<0.05) SumChi2_test += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2); + // + /*if(A_3[i][j][k] >= 1) term1 = C*(A_3[i][j][k]+B_3[i][j][k])/(A_3[i][j][k]*(C+1)); + else term1 = 0; + term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1)); + + if(term1 > 0.0 && term2 > 0.0){ + lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2); + }else if(term1==0 && term2 > 0.0){ + lnL += B_3[i][j][k]*log(term2); + }else {cout<<"WARNING -- term1 and term2 are negative"< +#include +#include +#include +#include + +#include "TVector2.h" +#include "TFile.h" +#include "TString.h" +#include "TF1.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TMath.h" +#include "TText.h" +#include "TRandom3.h" +#include "TArray.h" +#include "TLegend.h" +#include "TStyle.h" +#include "TMinuit.h" +#include "TCanvas.h" +#include "TPad.h" + +#define BohrR 1963.6885 +#define FmToGeV 0.19733 // conversion to fm +#define PI 3.1415926 +#define masspiC 0.1395702 // pi+ mass (GeV/c^2) + +using namespace std; + + + +void Makec3EAfile(){ + + + TFile *infile; + TMinuit *fit; + // + TH3D *PbPbEA = new TH3D("PbPbEA","",2,0.5,2.5, 4,0.5,4.5, 50,-0.5,49.5); + PbPbEA->SetDirectory(0); + // + TH3D *pPbEA = new TH3D("pPbEA","",2,0.5,2.5, 4,0.5,4.5, 50,-0.5,49.5); + pPbEA->SetDirectory(0); + // + TH3D *ppEA = new TH3D("ppEA","",2,0.5,2.5, 4,0.5,4.5, 50,-0.5,49.5); + ppEA->SetDirectory(0); + // + + + + + // + ////////////////////////////// + double value=0, value_e=0; + // + for(int Gindex=0; Gindex<46; Gindex++){ + // PbPb + for(int FT=0; FT<2; FT++){// EW or LG + + TString *name1 = new TString("FitFiles/FitFile_CT0_FT"); + *name1 += FT; + name1->Append("_G"); + *name1 += Gindex; + name1->Append(".root"); + infile = new TFile(name1->Data(),"READ"); + fit = (TMinuit*)infile->Get("MyMinuit_c3"); + for(int parNum=0; parNum<4; parNum++){ + fit->GetParameter(parNum+1, value,value_e); + PbPbEA->SetBinContent(FT+1, parNum+1, Gindex+1, value); + } + infile->Close(); + + // + // pPb + TString *name1 = new TString("FitFiles/FitFile_CT1_FT"); + *name1 += FT; + name1->Append("_G"); + *name1 += Gindex; + name1->Append(".root"); + infile = new TFile(name1->Data(),"READ"); + fit = (TMinuit*)infile->Get("MyMinuit_c3"); + for(int parNum=0; parNum<4; parNum++){ + fit->GetParameter(parNum+1, value,value_e); + pPbEA->SetBinContent(FT+1, parNum+1, Gindex+1, value); + } + infile->Close(); + // + // + TString *name1 = new TString("FitFiles/FitFile_CT2_FT"); + *name1 += FT; + name1->Append("_G"); + *name1 += Gindex; + name1->Append(".root"); + infile = new TFile(name1->Data(),"READ"); + fit = (TMinuit*)infile->Get("MyMinuit_c3"); + for(int parNum=0; parNum<4; parNum++){ + fit->GetParameter(parNum+1, value,value_e); + ppEA->SetBinContent(FT+1, parNum+1, Gindex+1, value); + } + infile->Close(); + } + } + // blank for the rest + for(int Gindex=46; Gindex<50; Gindex++){ + for(int FT=0; FT<2; FT++){// EW or LG + for(int parNum=0; parNum<4; parNum++){ + PbPbEA->SetBinContent(FT+1, parNum+1, Gindex+1, PbPbEA->GetBinContent(FT+1, parNum+1, 46)); + pPbEA->SetBinContent(FT+1, parNum+1, Gindex+1, pPbEA->GetBinContent(FT+1, parNum+1, 46)); + ppEA->SetBinContent(FT+1, parNum+1, Gindex+1, ppEA->GetBinContent(FT+1, parNum+1, 46)); + } + } + } + + // Convert Lam_3 to proper EA normalization + for(int Gindex=0; Gindex<50; Gindex++){ + for(int FT=0; FT<2; FT++){// EW or LG + PbPbEA->SetBinContent(FT+1, 1, Gindex+1, pow(PbPbEA->GetBinContent(FT+1, 1, Gindex+1)/ 2., 1/3.)); + pPbEA->SetBinContent(FT+1, 1, Gindex+1, pow(pPbEA->GetBinContent(FT+1, 1, Gindex+1)/ 2., 1/3.)); + ppEA->SetBinContent(FT+1, 1, Gindex+1, pow(ppEA->GetBinContent(FT+1, 1, Gindex+1)/ 2., 1/3.)); + } + } + + + TFile *outfile=new TFile("c3EAfile_temp.root","RECREATE"); + PbPbEA->Write(); + pPbEA->Write(); + ppEA->Write(); + // + outfile->Close(); + + +} diff --git a/PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C b/PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C index cb4024de5de..6679b8cc8e7 100755 --- a/PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C +++ b/PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C @@ -37,21 +37,24 @@ using namespace std; // 5(Linear instead of Cubic,also fc^2=0.75), 6(10h), 7(MRC 10% increase), 8(MuonCorr 92% pion-pair purity) int FileSetting=0; // -bool MCcase_def=kFALSE;// MC data? +bool MCcase_def=0;// MC data? +int CollisionType=1;// PbPb, pPb, pp // int Mbin_def=0;// 0-9: centrality bin in widths of 5% -bool SameCharge_def=0;// same-charge? +bool SameCharge_def=1;// same-charge? int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or -++, ---+ or -+++ -int MixedCharge4pionType_def = 1;// 1(---+) or 2(--++) +int MixedCharge4pionType_def = 2;// 1(---+) or 2(--++) // int EDbin_def=0;// 0 or 1: Kt3 bin -int TPNbin=3;// TPN bin for r3 and r4 +int TPNbin=0;// TPN bin for r3 and r4 int Gbin = int( (0) /2. ) + 55;// +5 (Rcoh=0), +25 (Rcoh=Rch) or +55 for extended G range +int c3FitType = 2;// EW(1), LG(2) int Ktbin_def=1;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range // bool MRC=1;// Momentum Resolution Corrections? bool MuonCorrection=1;// correct for Muons? -bool IncludeEW_def=kTRUE;// Include EdgeWorth coefficients? +bool FSICorrection=1;// correct For Final-State-Interactions +bool InterpCorrection=1;// correct for finite bin interpolation // int f_choice=0;// 0(Core/Halo), 1(40fm), 2(70fm), 3(100fm) // @@ -63,7 +66,6 @@ bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC? // // -bool PbPbcase=kTRUE;// always kTRUE int fFSIindex=0; int Ktlowbin; int Kthighbin; @@ -78,8 +80,8 @@ double NormQcutHigh_pp = 1.5; double NormQcutLow_PbPb = .15; double NormQcutHigh_PbPb = .2;// was .175 -int ThreeParticleRebin=2; -int FourParticleRebin=3; +int ThreeParticleRebin; +int FourParticleRebin; int RbinMRC; @@ -146,17 +148,92 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool ThreeFrac = pow(TwoFrac, 1.5); FourFrac = pow(TwoFrac, 2.); - cout<<"Mbin = "<ProjectionY(proname->Data(),Ktlowbin,Kthighbin);// bins:5-6 (kt:0.2-0.3) norm_2[term] = TwoParticle[c1][c2][term]->Integral(TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutLow),TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutHigh)); - //if(c1==c2) cout<<"2-pion norms "<Scale(norm_2[0]/norm_2[term]); TwoParticle[c1][c2][term]->SetMarkerStyle(20); @@ -351,6 +442,12 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool if(term<4){ K3avg[c1][c2][c3][term] = (TProfile*)MyList->FindObject(nameK3->Data()); K3avg[c1][c2][c3][term]->Rebin(ThreeParticleRebin); + if(MCcase || !FSICorrection) { + K3avg[c1][c2][c3][term]->Reset(); + for(int ii=1; ii<=K3avg[c1][c2][c3][term]->GetNbinsX(); ii++) { + K3avg[c1][c2][c3][term]->Fill(K3avg[c1][c2][c3][term]->GetXaxis()->GetBinCenter(ii), 1); + } + } } // if(term==4 && c1==c2 && c1==c3){ @@ -365,8 +462,20 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool for(int binx=1; binx<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsX(); binx++) { for(int biny=1; biny<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsY(); biny++) { TPFullWeight_ThreeParticle_2D[c1]->SetBinError(binx, biny, 0); + if(binx!=4){ + TPFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny)); + TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny)); + if(InterpCorrection){ + double q3 = TPFullWeight_ThreeParticle_2D[c1]->GetYaxis()->GetBinCenter(biny); + double InterCorr = 1.024 - 0.2765*q3; + //if(q3<0.1) cout<SetBinContent(binx, biny, InterCorr * TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny)); + TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny)); + } + } } } + TString *proName=new TString(nameTPN3->Data()); TString *proNameNeg=new TString(nameNegTPN3->Data()); proName->Append("_pro"); proNameNeg->Append("_pro"); TPN_ThreeParticle[c1] = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin); @@ -381,7 +490,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool tempDen->Add(tempDenNeg); TPFullWeight_ThreeParticle[c1]->Add(TPNegFullWeight_ThreeParticle[c1]); // - TPFullWeight_ThreeParticle[c1]->Add(tempDen); + //TPFullWeight_ThreeParticle[c1]->Add(tempDen);// now added above in Interp section TPFullWeight_ThreeParticle[c1]->Divide(tempDen); TPFullWeight_ThreeParticle[c1]->SetLineColor(1); // @@ -413,7 +522,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool nameNorm4->Append("_Norm"); // TString *nameK4=new TString(name4->Data()); - nameK4->Append("_Kfactor"); + nameK4->Append("_Kfactor");// or Weighted // TString *nameTPN4=new TString(name4->Data()); nameTPN4->Append("_TwoPartNorm"); @@ -421,6 +530,12 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TString *nameNegTPN4=new TString(name4->Data()); nameNegTPN4->Append("_TwoPartNegNorm"); // + TString *nameFitBuild=new TString(name4->Data()); + nameFitBuild->Append("_FullBuildFromFits"); + // + TString *namePartialFitBuild=new TString(name4->Data()); + namePartialFitBuild->Append("_PartialBuildFromFits"); + // name4->Append("_1D"); /////////////////////////////////////// // skip degenerate histograms @@ -440,13 +555,19 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool FourParticle[c1][c2][c3][c4][term]->SetTitle(""); // FourParticle[c1][c2][c3][c4][term]->Rebin(FourParticleRebin); - + }else{ - cout<<"normalization = 0."<FindObject(nameK4->Data()); K4avg[c1][c2][c3][c4][term]->Rebin(FourParticleRebin); + if(MCcase || !FSICorrection) { + K4avg[c1][c2][c3][c4][term]->Reset(); + for(int ii=1; ii<=K4avg[c1][c2][c3][c4][term]->GetNbinsX(); ii++) { + K4avg[c1][c2][c3][c4][term]->Fill(K4avg[c1][c2][c3][c4][term]->GetXaxis()->GetBinCenter(ii), 1); + } + } } if(term==12 && c1==c2 && c1==c3 && c1==c4){ @@ -458,13 +579,35 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TPNegFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]); TPNegFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin); // + if(c1==0){ + FullBuildFromFits_3D = (TH3D*)MyList->FindObject(nameFitBuild->Data()); + FullBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]); + FullBuildFromFits_3D->RebinZ(FourParticleRebin); + // + PartialBuildFromFits_3D = (TH3D*)MyList->FindObject(namePartialFitBuild->Data()); + PartialBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]); + PartialBuildFromFits_3D->RebinZ(FourParticleRebin); + } + // for(int binx=1; binx<=TPFullWeight_FourParticle_2D[c1]->GetNbinsX(); binx++) { for(int biny=1; biny<=TPFullWeight_FourParticle_2D[c1]->GetNbinsY(); biny++) { TPFullWeight_FourParticle_2D[c1]->SetBinError(binx, biny, 0); + if(binx!=4){ + TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny)); + TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny)); + + if(InterpCorrection){// Interpolator correction + double q4 = TPFullWeight_FourParticle_2D[c1]->GetYaxis()->GetBinCenter(biny); + double InterCorr = 1.032 - 0.2452*q4; + TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny)); + TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny)); + } + } } } TString *proName=new TString(nameTPN4->Data()); TString *proNegName=new TString(nameNegTPN4->Data()); - proName->Append("_pro"); proNegName->Append("_pro"); + TString *proNameFitBuild=new TString(nameFitBuild->Data()); TString *proNamePartialFitBuild=new TString(namePartialFitBuild->Data()); + proName->Append("_pro"); proNegName->Append("_pro"); proNameFitBuild->Append("_pro"); proNamePartialFitBuild->Append("_pro"); TPN_FourParticle[c1] = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin); // proName->Append("_FullWeight"); proNegName->Append("_FullWeight"); @@ -478,7 +621,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool tempDen->Add(tempDenNeg); TPFullWeight_FourParticle[c1]->Add(TPNegFullWeight_FourParticle[c1]); // - TPFullWeight_FourParticle[c1]->Add(tempDen); + //TPFullWeight_FourParticle[c1]->Add(tempDen);// now added above in Interp section TPFullWeight_FourParticle[c1]->Divide(tempDen); TPFullWeight_FourParticle[c1]->SetLineColor(1); /*TString *ErrName=new TString(nameTPN4->Data()); @@ -490,6 +633,19 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool cout<GetBinContent(3)<GetBinContent(5) / tempDen->GetBinContent(5))<<" "<GetBinContent(5)<ProjectionZ(proNameFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin); + TH1D *tempDen2 = (TH1D*)FullBuildFromFits_3D->ProjectionZ("tempDen2", c3FitType, c3FitType, 4, 4); + tempDen2->Scale(1/100.);// It was filled 100 times with the same value + FullBuildFromFits->Add(tempDen2); + FullBuildFromFits->Divide(tempDen2); + FullBuildFromFits->SetLineColor(1); + // + PartialBuildFromFits = (TH1D*)PartialBuildFromFits_3D->ProjectionZ(proNamePartialFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin); + PartialBuildFromFits->Add(tempDen2); + PartialBuildFromFits->Divide(tempDen2); + PartialBuildFromFits->SetLineColor(2); + } } }// term 4-particle @@ -500,6 +656,12 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TH1D *fTPNRejects3pion = (TH1D*)MyList->FindObject("fTPNRejects3pion2"); TH1D *fTPNRejects4pion1 = (TH1D*)MyList->FindObject("fTPNRejects4pion1"); + + //TH2D *c4QSFitNum2D = (TH2D*)MyList->FindObject("fc4QSFitNum"); + //TH2D *c4QSFitDen2D = (TH2D*)MyList->FindObject("fc4QSFitDen"); + //c4QSFitNum2D->RebinY(3); + //c4QSFitDen2D->RebinY(3); + cout<<"Done getting Histograms"<Multiply(MRC_MC_2[0]); TERM2_2->Multiply(MRC_MC_2[1]); } + TH1D *C2=(TH1D*)TERM1_2->Clone(); C2->Divide(TERM2_2); C2->GetXaxis()->SetRangeUser(0,0.1); - C2->SetMinimum(0.98); C2->SetMaximum(1.8); + C2->SetMinimum(0.98); C2->SetMaximum(2.5); + /*C2->GetXaxis()->SetTitleSize(0.06); + C2->GetXaxis()->SetLabelSize(0.06); + C2->GetYaxis()->SetTitleSize(0.06); + C2->GetYaxis()->SetLabelSize(0.06); + C2->GetXaxis()->SetTitleOffset(1.05); + C2->GetYaxis()->SetTitleOffset(1.1); + C2->GetXaxis()->SetNdivisions(606); + C2->GetYaxis()->SetNdivisions(505);*/ C2->Draw(); TH1D *C2QS=(TH1D*)TERM1_2->Clone(); C2QS->Add(TERM2_2, -(1-TwoFrac)); C2QS->Scale(1/TwoFrac); for(int bin=1; bin<=C2QS->GetNbinsX(); bin++) { - Float_t K2 = FSICorrelation(ch1_2, ch2_2, C2QS->GetXaxis()->GetBinCenter(bin)); + Float_t K2 = 1.0; + if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, C2QS->GetXaxis()->GetBinCenter(bin)); C2QS->SetBinContent(bin, C2QS->GetBinContent(bin)/K2); C2QS->SetBinError(bin, C2QS->GetBinError(bin)/K2); } C2QS->Divide(TERM2_2); if(SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionSC); if(!SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionMC); - C2QS->SetMarkerColor(2); - C2QS->Draw("same"); + C2QS->SetMarkerColor(2); C2QS->SetLineColor(2); + if(!MCcase) C2QS->Draw("same"); + + // Print 2-pion values + /*for(int bin=1; bin<=20; bin++){ + cout<GetBinContent(bin)<<", "; + } + cout<GetBinError(bin)<<", "; + } + cout<Clone(); + TH1D *C2ref_0p75=(TH1D*)C2QS->Clone(); + for(int bin=1; bin<=20; bin++){ + C2ref_0p65->SetBinContent(bin, C2refPoints_0p65[bin-1]); + C2ref_0p65->SetBinError(bin, C2refPoints_e_0p65[bin-1]); + C2ref_0p75->SetBinContent(bin, C2refPoints_0p75[bin-1]); + C2ref_0p75->SetBinError(bin, C2refPoints_e_0p75[bin-1]); + } + C2ref_0p65->SetMarkerColor(4); C2ref_0p65->SetLineColor(4); + C2ref_0p75->SetMarkerColor(6); C2ref_0p75->SetLineColor(6); + //C2ref_0p65->Draw("same"); + //C2ref_0p75->Draw("same"); + + /////////////////////////////////////////////////////////// @@ -610,6 +810,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TERM3_3->Multiply(MRC_SC_3[1]); TERM4_3->Multiply(MRC_SC_3[1]); TERM5_3->Multiply(MRC_SC_3[2]); + //cout<GetBinContent(2)<<" "<GetBinContent(3)<<" "<GetBinContent(4)<Multiply(MRC_MC_3[0]); @@ -626,16 +827,6 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool } } - // - float f33 = ThreeFrac; - float f32 = (1-OneFrac); - float f31 = pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) - 3*(1-TwoFrac)*(1-OneFrac); - float TherminatorMod_f33[4]={1., 1.123, 1.022, 1.008};// Core/Halo, 40fm, 70fm, 100fm - float TherminatorMod_f32[4]={1., 0.844, 0.933, 0.967}; - float TherminatorMod_f31[4]={1., 0.828, 0.982, 1.028}; - f33 *= TherminatorMod_f33[f_choice]; - f32 *= TherminatorMod_f32[f_choice]; - f31 *= TherminatorMod_f31[f_choice]; TH1D *N3QS = (TH1D*)TERM1_3->Clone(); N3QS->Add(TERM5_3, -f31); @@ -681,7 +872,8 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool value -= ( TERM2_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][1]->GetBinContent(bin) / TwoFrac * MuonCorr2; value -= ( TERM3_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][2]->GetBinContent(bin) / TwoFrac * MuonCorr3; value -= ( TERM4_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][3]->GetBinContent(bin) / TwoFrac * MuonCorr4; - value -= TERM5_3->GetBinContent(bin)*(MuonCorr1 - MuonCorr2 - MuonCorr3 - MuonCorr4); + //value -= TERM5_3->GetBinContent(bin)*(MuonCorr1 - MuonCorr2 - MuonCorr3 - MuonCorr4); + value += 2*TERM5_3->GetBinContent(bin); // n3QS->SetBinContent(bin, value); c3QS->SetBinContent(bin, value + TERM5_3->GetBinContent(bin)); @@ -690,12 +882,14 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool C3QS->Draw(); c3QS->Draw("same"); - int lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); - int highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); - double SF=C3QS->Integral(lowBin, highBin); SF /= TPFullWeight_ThreeParticle[ch1_3]->Integral(lowBin, highBin); + //int lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); + //int highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); + //double SF=C3QS->Integral(lowBin, highBin); SF /= TPFullWeight_ThreeParticle[ch1_3]->Integral(lowBin, highBin); //TPFullWeight_ThreeParticle[ch1_3]->Scale(SF); + // if(SameCharge) TPFullWeight_ThreeParticle[ch1_3]->Draw("same"); + //cout<GetBinContent(2)<Clone(); //C3raw->Divide(TERM5_3); //C3raw->SetMarkerColor(4); @@ -706,7 +900,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool if(SameCharge) legend2->AddEntry(TPFullWeight_ThreeParticle[ch1_3],"C_{3}^{QS} built","l"); legend2->Draw("same"); - + /////////////////////////////////////// // C3QS built ratio TCanvas *can2_2 = new TCanvas("can2_2", "can2_2",1200,53,700,500); @@ -733,9 +927,9 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TH1D *Chi2NDF_FullSize_3 = new TH1D("Chi2NDF_FullSize_3","",100,-0.5,99.5); Chi2NDF_PointSize_3->SetLineColor(4); Chi2NDF_FullSize_3->SetLineColor(2); Chi2NDF_PointSize_3->SetMarkerColor(4); Chi2NDF_FullSize_3->SetMarkerColor(2); - Chi2NDF_PointSize_3->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_3->GetYaxis()->SetTitle("3-pion #chi^{2} / NDF"); + Chi2NDF_PointSize_3->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_3->GetYaxis()->SetTitle("#sqrt{#chi^{2}}"); - if(SameCharge){ + if(SameCharge && CollisionType==0){ TH1D *tempDen = (TH1D*)TPFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPFullWeight3_Den", 4, 4); TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPNegFullWeight3_Den", 4, 4); @@ -751,12 +945,12 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // Add Pos and Neg Num tempNum->Add(tempNumNeg); // - tempNum->Add(tempDen); + //tempNum->Add(tempDen);// now added in histogram read section tempNum->Divide(tempDen); - lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); - highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); - SF=C3QS->Integral(lowBin, highBin); - SF /= tempNum->Integral(lowBin, highBin); + //lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); + //highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]); + //SF=C3QS->Integral(lowBin, highBin); + //SF /= tempNum->Integral(lowBin, highBin); //tempNum->Scale(SF); double Chi2=0; @@ -786,7 +980,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool C3QSratio->GetYaxis()->SetTitle("C_{3}^{QS} / C_{3}(built)"); C3QSratio->Draw(); Unity->Draw("same"); - + Chi2NDF_PointSize_3->Draw(); Chi2NDF_FullSize_3->Draw("same"); legend2_2->AddEntry(Chi2NDF_PointSize_3,"R_{coh}=0 (Point Source)","p"); @@ -796,7 +990,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // r3 TH1D *r3; - if(SameCharge){ + if(SameCharge && CollisionType==0){ r3 = (TH1D*)n3QS->Clone(); TPN_ThreeParticle[ch1_3]->Multiply(MRC_SC_3[2]); r3->Divide(TPN_ThreeParticle[ch1_3]); @@ -804,7 +998,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool r3->SetMinimum(0.5); r3->SetMaximum(2.5); r3->GetYaxis()->SetTitle("r_{3}"); // - //r3->Draw(); + r3->Draw(); //TPN_ThreeParticle[ch1_3]->Draw(); //fTPNRejects3pion->SetLineColor(2); //fTPNRejects3pion->Draw("same"); @@ -829,7 +1023,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // 4-pion cout<<"4-pion section"<SetHighLightColor(2); can3->Range(-0.7838115,-1.033258,0.7838115,1.033258); gStyle->SetOptFit(0111); @@ -873,9 +1067,9 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool }else{// --++ if(index==0) MRC_temp = (TH1D*)MRC_MC2_4[0]->Clone(); else if(index<=4) MRC_temp = (TH1D*)MRC_MC2_4[1]->Clone(); - else if(index==5 || index==10) MRC_temp = (TH1D*)MRC_MC2_4[2]->Clone(); - else if(index<=9) MRC_temp = (TH1D*)MRC_MC2_4[3]->Clone(); - else if(index==11) MRC_temp = (TH1D*)MRC_MC2_4[4]->Clone(); + else if(index==5 || index==10) MRC_temp = (TH1D*)MRC_MC2_4[2]->Clone();//2 + else if(index<=9) MRC_temp = (TH1D*)MRC_MC2_4[3]->Clone();// 3 + else if(index==11) MRC_temp = (TH1D*)MRC_MC2_4[4]->Clone();// 4 else MRC_temp = (TH1D*)MRC_MC2_4[5]->Clone(); } // @@ -885,22 +1079,8 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TH1D *C4raw=(TH1D*)TERMS_4[2]->Clone(); C4raw->Divide(TERMS_4[12]); C4raw->GetXaxis()->SetRangeUser(0,0.2); - - float f44 = FourFrac; - float f43 = (1-OneFrac); - float f42 = -pow(1-OneFrac,2); - float f41 = -3*pow(1-OneFrac,4) - 8*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2); - float TherminatorMod_f44[4]={1., 1.188, 1.008, 0.985};// Core/Halo, 40fm, 70fm, 100fm - float TherminatorMod_f43[4]={1., 0.885, 0.979, 1.026}; - float TherminatorMod_f42[4]={1., 0.687, 0.99, 1.08}; - float TherminatorMod_f41[4]={1., 0.806, 1.33, 1.548}; - f44 *= TherminatorMod_f44[f_choice]; - f43 *= TherminatorMod_f43[f_choice]; - f42 *= TherminatorMod_f42[f_choice]; - f41 *= TherminatorMod_f41[f_choice]; - //cout<Clone(); // N4QS->Add(TERMS_4[1], -f43); @@ -918,16 +1098,17 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool N4QS->Add(TERMS_4[12], -f41); // N4QS->Scale(1/f44); - //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Scale(1.005);// K factorization + //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Scale(1.005);// K scale variation N4QS->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]); - //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]); - //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]); + //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]);// K factorization test (MC1) + //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]);// K factorization test (MC2) TH1D *C4QS=(TH1D*)N4QS->Clone(); C4QS->GetXaxis()->SetRangeUser(0,0.179); C4QS->SetMinimum(0.5); C4QS->SetMarkerColor(4); + // /*TH1D *C4QS_basic=(TH1D*)TERMS_4[0]->Clone(); for(int ii=1; ii<=C4QS_basic->GetNbinsX(); ii++){ @@ -1028,7 +1209,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // } // 2-pion terms - SubtractionTerm[4] = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac; + SubtractionTerm[4] = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin); ErrorTerm[5] = sqrt( TERMS_4[5]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac; // SubtractionTerm[5] = ( TERMS_4[6]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][6]->GetBinContent(bin) / TwoFrac; @@ -1052,10 +1233,13 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool SubtractionTerm[10] = TERMS_4[11]->GetBinContent(bin);// N22 SubtractionTerm[10] -= 2*(1-TwoFrac) * TERMS_4[5]->GetBinContent(bin); SubtractionTerm[10] += pow(1-TwoFrac,2) * TERMS_4[12]->GetBinContent(bin); - // SubtractionTerm[10] /= pow(TwoFrac,2); SubtractionTerm[10] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]->GetBinContent(bin); - SubtractionTerm[10] -= 2*( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac; + SubtractionTerm[10] *= C4muonCorrectionSC[3]->GetBinContent(bin); + double intermed = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac; + intermed *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin); + intermed = intermed * C4muonCorrectionSC[2]->GetBinContent(bin); + SubtractionTerm[10] -= 2*intermed; SubtractionTerm[10] += 2*TERMS_4[12]->GetBinContent(bin); // ErrorTerm[11] = TERMS_4[11]->GetBinContent(bin); @@ -1072,18 +1256,31 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // double FinalValue[4]={0}; double FinalValue_e[4]={0}; - if(SameCharge) {N4QSvalue *= C4muonCorrectionSC[0]->GetBinContent(bin); ErrorTerm[0] *= C4muonCorrectionSC[0]->GetBinContent(bin);} - else if(MixedCharge4pionType==1) {N4QSvalue *= C4muonCorrectionMC1[0]->GetBinContent(bin); ErrorTerm[0] *= C4muonCorrectionMC1[0]->GetBinContent(bin);} - else {N4QSvalue *= C4muonCorrectionMC2[0]->GetBinContent(bin); ErrorTerm[0] *= C4muonCorrectionMC2[0]->GetBinContent(bin);} + if(SameCharge) { + //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin); + N4QSvalue *= C4muonCorrectionSC[0]->GetBinContent(bin); + ErrorTerm[0] *= C4muonCorrectionSC[0]->GetBinContent(bin); + }else if(MixedCharge4pionType==1) { + //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin); + N4QSvalue *= C4muonCorrectionMC1[0]->GetBinContent(bin); + ErrorTerm[0] *= C4muonCorrectionMC1[0]->GetBinContent(bin); + }else { + //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin); + N4QSvalue *= C4muonCorrectionMC2[0]->GetBinContent(bin); + ErrorTerm[0] *= C4muonCorrectionMC2[0]->GetBinContent(bin); + } for(int ii=0; ii<4; ii++) { FinalValue[ii] = N4QSvalue; FinalValue_e[ii] = pow(ErrorTerm[0],2); } if(SameCharge) { - FinalValue[0] -= 4*(SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[1]->GetBinContent(bin); - FinalValue[0] += 6*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[2]->GetBinContent(bin); - FinalValue[0] -= 3*(SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[3]->GetBinContent(bin); + //FinalValue[0] -= 4*(SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[1]->GetBinContent(bin); + //FinalValue[0] += 6*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[2]->GetBinContent(bin); + FinalValue[0] -= 4*(SubtractionTerm[0] * C4muonCorrectionSC[1]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); + FinalValue[0] += 6*(SubtractionTerm[4] * C4muonCorrectionSC[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); + FinalValue[0] -= 3*(SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin)); + FinalValue_e[0] += 4*pow(ErrorTerm[1]*C4muonCorrectionSC[1]->GetBinContent(bin),2); FinalValue_e[0] += 6*pow(ErrorTerm[5]*C4muonCorrectionSC[2]->GetBinContent(bin),2); FinalValue_e[0] += 3*pow(ErrorTerm[11]*C4muonCorrectionSC[3]->GetBinContent(bin),2); @@ -1096,18 +1293,26 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool FinalValue[3] = SubtractionTerm[10]; }else{ if(MixedCharge4pionType==1 && CHARGE==-1){ - FinalValue[0] -= (SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin); + //FinalValue[0] -= (SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin);// [2] is +++ MuonCorr + FinalValue[0] -= SubtractionTerm[0] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin);// [2] is +++ MuonCorr + FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); + // FinalValue_e[0] += pow(ErrorTerm[1]*C4muonCorrectionMC1[2]->GetBinContent(bin),2); FinalValue[1] -= 3*SubtractionTerm[4] - 3*TERMS_4[12]->GetBinContent(bin); FinalValue_e[1] += 3*pow(ErrorTerm[5],2) + 3*TERMS_4[12]->GetBinContent(bin); }else if(MixedCharge4pionType==1 && CHARGE==+1){ - FinalValue[0] -= (SubtractionTerm[8] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin); + //FinalValue[0] -= (SubtractionTerm[3] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin); + FinalValue[0] -= SubtractionTerm[3] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin); + FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); + // FinalValue_e[0] += pow(ErrorTerm[4]*C4muonCorrectionMC1[2]->GetBinContent(bin),2); FinalValue[1] -= 3*SubtractionTerm[9] - 3*TERMS_4[12]->GetBinContent(bin); FinalValue_e[1] += 3*pow(ErrorTerm[10],2) + 3*TERMS_4[12]->GetBinContent(bin); }else{// --++ case - FinalValue[0] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin); - FinalValue[0] -= (SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[4]->GetBinContent(bin); + //FinalValue[0] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin);// old way + FinalValue[0] -= 2*(SubtractionTerm[4] * C4muonCorrectionMC2[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); + FinalValue[0] -= (SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin)); + FinalValue[0] -= 4*(SubtractionTerm[6] * C4muonCorrectionMC2[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin)); FinalValue_e[0] += 2*pow(ErrorTerm[5],2) + pow(ErrorTerm[11],2) + 2*TERMS_4[12]->GetBinContent(bin); // FinalValue[1] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin); @@ -1148,13 +1353,15 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool c4QS_RemovalStage3->SetMarkerColor(7); c4QS_RemovalStage3->SetLineColor(7); // // - + if(SameCharge) C4QS->SetMaximum(9); else if(MixedCharge4pionType==1) C4QS->SetMaximum(4); else C4QS->SetMaximum(3); + if(CollisionType!=0) C4QS->GetXaxis()->SetRangeUser(0,0.6); + C4QS->Draw(); c4QS_RemovalStage1->Draw("same"); if(SameCharge) c4QS_RemovalStage2->Draw("same"); @@ -1162,12 +1369,13 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool c4QS->Draw("same"); - lowBin = C4QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q4[Mbin]); - highBin = C4QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q4[Mbin]); - SF=C4QS->Integral(lowBin, highBin); SF /= TPFullWeight_FourParticle[ch1_4]->Integral(lowBin, highBin); - //TPFullWeight_FourParticle[ch1_4]->Scale(SF); - if(SameCharge) TPFullWeight_FourParticle[ch1_4]->Draw("same"); - + //cout<GetBinContent(9)<Draw("same"); + FullBuildFromFits->Draw("same"); + PartialBuildFromFits->Draw("same"); + } legend3->AddEntry(C4QS,"C_{4}^{QS}","p"); legend3->AddEntry(c4QS_RemovalStage1,"c_{4}^{QS}{2-pion removal}","p"); if(SameCharge) legend3->AddEntry(c4QS_RemovalStage2,"c_{4}^{QS}{2-pion+2-pair removal}","p"); @@ -1181,20 +1389,27 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Draw(); - TH1D *C4QS_Syst = (TH1D*)C4QS->Clone(); - C4QS_Syst->SetMarkerSize(0); - C4QS_Syst->SetFillStyle(1000); - C4QS_Syst->SetFillColor(kBlue-10); - for(int bin=1; bin<=C4QS_Syst->GetNbinsX(); bin++){ - double SysPercent = 0.06*C4QS_Syst->GetBinCenter(bin);// f coefficients - C4QS_Syst->SetBinError(bin, SysPercent * C4QS_Syst->GetBinContent(bin)); - } + + + //TH1D *c4QSFitNum = (TH1D*)c4QSFitNum2D->ProjectionY("c4QSFitNum",5,5); + //TH1D *c4QSFitDen = (TH1D*)c4QSFitDen2D->ProjectionY("c4QSFitDen",5,5); + //c4QSFitNum->Divide(c4QSFitDen); + //for(int bin=1; bin<20; bin++){ + //double fc4 = + //c4QSFitNum->SetBinContent(bin, (c4QSFitNum->GetBinContent(bin) - 1.0)* + //} + //c4QSFitNum->Draw("same"); + //C4QS_Syst->Draw("E2 same"); //C4QS->Draw("same"); //C4raw->Draw(); - + //TH1D *testTerm = (TH1D*)TERMS_4[11]->Clone(); + //testTerm->Divide(TERMS_4[12]); + //testTerm->SetMinimum(1); testTerm->SetMaximum(1.4); testTerm->GetXaxis()->SetRangeUser(0,0.14); + //testTerm->Draw(); + //////////////////////////////////////////////////////////////// - if(SameCharge){ + if(SameCharge && CollisionType==0){ TCanvas *can4 = new TCanvas("can4", "can4",600,600,700,500); can4->SetHighLightColor(2); can4->Range(-0.7838115,-1.033258,0.7838115,1.033258); @@ -1217,7 +1432,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TH1D *Chi2NDF_FullSize_4 = new TH1D("Chi2NDF_FullSize_4","",100,-0.5,99.5); Chi2NDF_PointSize_4->SetLineColor(4); Chi2NDF_FullSize_4->SetLineColor(2); Chi2NDF_PointSize_4->SetMarkerColor(4); Chi2NDF_FullSize_4->SetMarkerColor(2); - Chi2NDF_PointSize_4->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_4->GetYaxis()->SetTitle("4-pion #chi^{2} / NDF"); + Chi2NDF_PointSize_4->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_4->GetYaxis()->SetTitle("#sqrt{#chi^{2}}"); TH1D *tempDen = (TH1D*)TPFullWeight_FourParticle_2D[ch1_4]->ProjectionY("TPFullWeight4_Den", 4, 4); @@ -1234,13 +1449,9 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // Add Pos and Neg Weights tempNum->Add(tempNumNeg); // - tempNum->Add(tempDen); + //tempNum->Add(tempDen);// now added in InterpCorr section tempNum->Divide(tempDen); - lowBin = C4QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q4[Mbin]); - highBin = C4QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q4[Mbin]); - SF=C4QS->Integral(lowBin, highBin); - SF /= tempNum->Integral(lowBin, highBin); - //tempNum->Scale(SF); + double Chi2=0; double NDF=0; @@ -1248,7 +1459,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool if(tempNum->GetBinContent(binQ4) <=0) continue; double value = C4QS->GetBinContent(binQ4) - tempNum->GetBinContent(binQ4); double err = pow(C4QS->GetBinError(binQ4),2); - //err += pow(C4QS_Syst->GetBinError(binQ4),2); + err = sqrt(err); if(err<=0) continue; @@ -1276,7 +1487,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool // Print 4-pion points for(int bin=1; bin<=12; bin++){ - cout<GetBinContent(bin)<<", "; + //cout<GetBinContent(bin)<<", "; //cout<GetBinContent(bin)<<", "; //cout<GetBinContent(bin)<<", "; //cout<GetBinContent(bin)<<", "; @@ -1284,6 +1495,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool } cout<GetBinContent(bin)<<", "; //cout<GetBinError(bin)<<", "; //cout<GetBinError(bin)<<", "; //cout<GetBinError(bin)<<", "; @@ -1295,7 +1507,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool TF1 *ChaoticLimit_r42 = new TF1("ChaoticLimit_r42","6",0,10); ChaoticLimit_r42->SetLineStyle(2); TH1D *r42; - if(SameCharge){ + if(SameCharge && CollisionType==0){ /*TCanvas *can5 = new TCanvas("can5", "can5",1200,600,700,500); can5->SetHighLightColor(2); can5->Range(-0.7838115,-1.033258,0.7838115,1.033258); @@ -1430,6 +1642,7 @@ void SetMomResCorrections(){ TString *momresfilename = new TString("MomResFile"); if(FileSetting==7) momresfilename->Append("_10percentIncrease"); + if(CollisionType!=0) momresfilename->Append("_ppAndpPb"); momresfilename->Append(".root"); TFile *MomResFile = new TFile(momresfilename->Data(),"READ"); @@ -1484,7 +1697,7 @@ void SetMomResCorrections(){ MRC_MC2_4[0]->SetDirectory(0); MRC_MC2_4[1]->SetDirectory(0); MRC_MC2_4[2]->SetDirectory(0); MRC_MC2_4[3]->SetDirectory(0); MRC_MC2_4[4]->SetDirectory(0); MRC_MC2_4[5]->SetDirectory(0); - if(!MRC){ + if(!MRC || MCcase_def){ for(int bin=1; bin<=MRC_SC_2[0]->GetNbinsX(); bin++){ MRC_SC_2[0]->SetBinContent(bin, 1.0); MRC_SC_2[1]->SetBinContent(bin, 1.0); MRC_MC_2[0]->SetBinContent(bin, 1.0); MRC_MC_2[1]->SetBinContent(bin, 1.0); @@ -1535,7 +1748,7 @@ float fact(float n){ void SetMuonCorrections(){ TString *name = new TString("MuonCorrection"); if(FileSetting==8) name->Append("_92percent"); - + if(CollisionType!=0) name->Append("_ppAndpPb"); name->Append(".root"); TFile *MuonFile=new TFile(name->Data(),"READ"); TString *proName[22]; @@ -1580,7 +1793,7 @@ void SetMuonCorrections(){ C4muonCorrectionMC2[2]->SetDirectory(0); C4muonCorrectionMC2[3]->SetDirectory(0); C4muonCorrectionMC2[4]->SetDirectory(0); // // - if(!MuonCorrection){ + if(!MuonCorrection || MCcase_def){ for(int bin=1; bin<=C2muonCorrectionSC->GetNbinsX(); bin++){ C2muonCorrectionSC->SetBinContent(bin, 1.0); C2muonCorrectionMC->SetBinContent(bin, 1.0); @@ -1606,7 +1819,7 @@ void SetMuonCorrections(){ //________________________________________________________________________ void SetFSIindex(Float_t R){ if(!MCcase_def){ - if(PbPbcase){ + if(CollisionType==0){ if(Mbin_def==0) fFSIindex = 0;//0-5% else if(Mbin_def==1) fFSIindex = 1;//5-10% else if(Mbin_def<=3) fFSIindex = 2;//10-20% diff --git a/PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C b/PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C index 67a01964cdc..c70b33d9d8b 100755 --- a/PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C +++ b/PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C @@ -33,27 +33,28 @@ #define FmToGeV 0.19733 // conversion to fm #define PI 3.1415926 #define masspiC 0.1395702 // pi+ mass (GeV/c^2) -double kappa3 = 0.1; // 0.16 (default), 0.05 or 0.3 as variations -double kappa4 = 0.5; // 0.4 (default), +double kappa3 = 0.1; // 0.1 (default), +double kappa4 = 0.5; // 0.5 (default), using namespace std; -int CollisionType_def=1;// 0=PbPb, 1=pPb, 2=pp +int CollisionType_def=2;// 0=PbPb, 1=pPb, 2=pp bool MCcase_def=0;// MC data? int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case bool SameCharge_def=kTRUE;// 3-pion same-charge? bool AddCC=kTRUE; -bool Gaussian=1;// Gaussian or Exponential? +bool Gaussian=0;// Gaussian or Exponential? bool UseC2Bkg=1;// use Pythia/DPMJET bkg? // bool MuonCorrection=1;// apply Muon correction? bool IncludeExpansion=kTRUE;// Include EdgeWorth coefficients? bool FixExpansionAvg=1; -int Mbin_def=18;// 0-19 (0=1050-2000 pions, 19=0-5 pions) +int Mbin_def=16;// 0-19 (0=1050-2000 pions, 19=0-5 pions) int EDbin_def=0;// 0-2: Kt3 bin int Ktbin_def=1;// 1-6, 10=full range double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase bool FitRangeShift=0;// 30% reduction in Fit Range +bool ExtremeFitRangeC2=0;// 1.2 GeV/c fit range for C2? // // bool SaveToFile_def=kFALSE;// Save outputs to file? @@ -186,33 +187,32 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll SameCharge_def=SameCharge;// 3-pion same-charge Mbin_def=Mbin; // + if(CollisionType==0){ + Q3Limit = 0.1 + 0.1*Mbin/16.; + }else{ + Q3Limit = 0.3 + 0.2*fabs(Mbin-10)/9.; + } + if(FitRangeShift) Q3Limit -= 0.3*Q3Limit; + Q2Limit = Q3Limit/sqrt(2.); + if(ExtremeFitRangeC2) Q2Limit = 1.2; + //Q2Limit = (0.3 + 0.2*fabs(Mbin-10)/9.) / sqrt(2.); + //cout<SetRightMargin(0.025); gPad->SetLeftMargin(0.1); gPad->SetTopMargin(0.02); gPad->SetGridx(0); gPad->SetGridy(0); - TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC"); + TLegend *legend1 = new TLegend(.4,.67,.87,.87,NULL,"brNDC"); legend1->SetBorderSize(1); legend1->SetTextSize(.04);// small .03; large .036 legend1->SetFillColor(0); @@ -740,10 +738,15 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll legend1->Draw("same"); } - //for(int i=0; iGetNbinsX(); i++){ - //cout<GetBinError(i+1)<<", "; + //for(int i=0; i<100; i++){ + //cout<GetBinContent(i+1)<<", "; //} + TF1 *C2osFit_CMS = new TF1("C2osFit_CMS","[0] + [1]*exp(-pow([2]*x,2))",0,2); + C2osFit_CMS->FixParameter(0, 9.93558e-01); + C2osFit_CMS->FixParameter(1, 5.19578e-02); + C2osFit_CMS->FixParameter(2, 1.89425e+00); + ///////////////////////////////////////////////////// // Global fitting C2os and C2ss double C2ssSys_e[BINRANGE_C2global]={0}; @@ -808,6 +811,8 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll // // method with undilution for plotting purposes //C2ssFitting[i] = (C2ssFitting[i] - (1-lambda_PM)) / (CoulCorr2(+1, BinCenters[i]) * lambda_PM); + // divide CMS background + //C2ssFitting[i] /= C2osFit_CMS->Eval(BinCenters[i]); // K2SS[i] = CoulCorr2(+1, BinCenters[i]); K2OS[i] = CoulCorr2(-1, BinCenters[i]); @@ -845,7 +850,7 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll stepSize[0] = 0.01; stepSize[1] = 0.01; stepSize[2] = 0.02; stepSize[3] = 0.2; stepSize[4] = 0.01; stepSize[5] = 0.001; stepSize[6] = 0.001; stepSize[7] = 0.001; stepSize[8]=0.001; stepSize[9]=0.01; minVal[0] = 0.955; minVal[1] = 0.5; minVal[2] = 0.; minVal[3] = 0.1; minVal[4] = 0.001; minVal[5] = -10.; minVal[6] = -10.; minVal[7] = -10.; minVal[8]=-10; minVal[9] = 0.995;// minVal[1] was 0.2 maxVal[0] = 1.1; maxVal[1] = 1.3; maxVal[2] = 0.99; maxVal[3] = 15.; maxVal[4] = 2.; maxVal[5] = 10.; maxVal[6] = 10.; maxVal[7] = 10.; maxVal[8]=10.; maxVal[9]=1.1;// maxVal[1] was 1.0 - if(Gaussian==kFALSE) maxVal[1]=2.0; + if(Gaussian==kFALSE) {maxVal[1]=2.0; maxVal[3]=20.0;} parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh"; parName[5] = "coeff_3"; parName[6] = "coeff_4"; parName[7] = "coeff_5"; parName[8]="coeff_6"; parName[9]="Norm_2"; @@ -867,9 +872,13 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 }else{// IncludeExpansion if(FixExpansionAvg){ - MyMinuit.DefineParameter(5, parName[5].c_str(), kappa3, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 - MyMinuit.DefineParameter(6, parName[6].c_str(), kappa4, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 - + if(Gaussian){ + MyMinuit.DefineParameter(5, parName[5].c_str(), kappa3, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 + MyMinuit.DefineParameter(6, parName[6].c_str(), kappa4, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 + }else{ + MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3 + MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4 + } }else{ Int_t err=0; Double_t tmp[1]; @@ -909,6 +918,10 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll fitC2ss_Base->FixParameter(i,OutputPar[i]); fitC2ss_Base->SetParError(i,OutputPar_e[i]); } + for(int bin=1; bin<=Two_ex_clone_mm->GetNbinsX(); bin++){ + double qinv = Two_ex_clone_mm->GetXaxis()->GetBinCenter(bin); + if(!MCcase) fitC2ss_h->SetBinContent(bin, fitC2ss_Base->Eval(qinv));// Base fit + } }else{// Edgeworth or Laguerre for(int i=0; iSetParError(i,OutputPar_e[i]); } for(int bin=1; bin<=Two_ex_clone_mm->GetNbinsX(); bin++){ - double qinv = Two_ex_clone_mm->GetXaxis()->GetBinCenter(bin); - if(!MCcase) fitC2ss_h->SetBinContent(bin, fitC2ss_Expan->Eval(qinv));// uncomment to show expansion + //double qinv = Two_ex_clone_mm->GetXaxis()->GetBinCenter(bin); + //if(!MCcase) fitC2ss_h->SetBinContent(bin, fitC2ss_Expan->Eval(qinv));// uncomment to show expansion } double lambdaStar=0, lambdaStar_e=0; @@ -985,7 +998,7 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll C2_ss->GetXaxis()->SetTitleOffset(.8); C2_ss->GetYaxis()->SetTitleOffset(.85); C2_ss->GetXaxis()->SetTitle("#font[12]{q} (GeV/c)"); - C2_ss->GetYaxis()->SetTitle("#font[12]{C}_{2}^{#pm#pm}"); + C2_ss->GetYaxis()->SetTitle("#font[12]{C}_{2}"); C2_ss->GetXaxis()->SetTitleSize(0.05); C2_ss->GetYaxis()->SetTitleSize(0.05); C2_os->GetXaxis()->SetRangeUser(0,0.6); @@ -1004,15 +1017,43 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll fitC2ss_h->SetLineWidth(2); fitC2ss_h->SetLineColor(2); - + TH1D *temp_hist_ss=new TH1D("temp_hist_ss","",100,0,0.5); + TH1D *temp_hist_os=new TH1D("temp_hist_os","",100,0,0.5); + TH1D *temp_fit=new TH1D("temp_fit","",100,0,0.5); + + // PbPb, M16 + double temp_points_ss[100]={-1000, 1.31587, 1.40841, 1.4426, 1.36602, 1.37548, 1.34059, 1.30318, 1.27143, 1.28487, 1.26105, 1.21905, 1.2025, 1.19743, 1.15297, 1.16947, 1.14259, 1.13562, 1.13135, 1.10088, 1.10227, 1.08037, 1.0824, 1.08069, 1.0657, 1.06808, 1.05341, 1.0493, 1.05081, 1.04546, 1.04644, 1.03703, 1.03975, 1.04479, 1.03218, 1.02657, 1.02588, 1.02756, 1.02345, 1.02306, 1.02024, 1.01551, 1.01133, 1.01025, 1.00697, 1.01819, 1.00651, 1.01293, 1.00529, 0.997971, 1.00141, 1.00665, 1.01479, 1.01668, 1.00657, 1.01188, 1.0042, 0.995877, 1.01094, 1.0036, 1.00312, 0.994031, 1.0084, 1.00604, 1.00329, 1.00677, 0.997553, 0.992874, 0.993799, 0.995368, 0.999215, 1.00508, 0.993053, 0.996744, 0.990297, 0.990715, 1.00114, 0.996376, 0.984527, 0.99474, 0.993457, 1.00723, 0.999042, 0.996589, 0.988385, 0.9916, 0.995184, 0.988197, 1.00019, 0.998052, 0.985506, 0.99024, 0.996102, 0.995797, 0.991887, 0.99276, 0.992561, 0.991261, 0.996243}; + double temp_points_ss_e[100]={1000, 0.206703, 0.0920405, 0.0637848, 0.0465513, 0.0380257, 0.0318382, 0.0270753, 0.0236593, 0.0216127, 0.0195971, 0.0176654, 0.0163346, 0.0153927, 0.0141902, 0.0136715, 0.0128887, 0.0124019, 0.0119795, 0.0113642, 0.0111615, 0.0106718, 0.0104775, 0.0102842, 0.0100296, 0.00990822, 0.00962737, 0.0094588, 0.00935572, 0.00923423, 0.0091266, 0.00897368, 0.00889453, 0.00883246, 0.00864815, 0.00853381, 0.00846384, 0.00839173, 0.00829825, 0.00822044, 0.00817065, 0.00805352, 0.00798793, 0.00791266, 0.00784862, 0.00788769, 0.00775376, 0.00772075, 0.00766719, 0.00757909, 0.00754696, 0.00755056, 0.00756756, 0.00754031, 0.00746571, 0.00745114, 0.00739226, 0.00729679, 0.00737659, 0.00730319, 0.00730321, 0.0072228, 0.00728131, 0.00726075, 0.0072364, 0.00723868, 0.00716992, 0.00712198, 0.00713157, 0.00713724, 0.00714391, 0.007175, 0.00708838, 0.00711573, 0.00709965, 0.00707786, 0.00713124, 0.00712164, 0.00703573, 0.00712383, 0.00709407, 0.00720098, 0.00713966, 0.00716466, 0.00710145, 0.00711087, 0.00714815, 0.00710751, 0.00717636, 0.00718936, 0.00712875, 0.00715855, 0.00718908, 0.00719846, 0.00719133, 0.00720147, 0.00721384, 0.00720763, 0.0072536}; + + + for(int bin=1; bin<100; bin++){ + //cout<GetBinContent(bin)<<", "; + //cout<GetBinError(bin)<<", "; + //cout<GetBinContent(bin)<<", "; + temp_hist_ss->SetBinContent(bin, temp_points_ss[bin-1]); + temp_hist_ss->SetBinError(bin, temp_points_ss_e[bin-1]); + //temp_hist_os->SetBinContent(bin, temp_points_os[bin-1]); + //temp_hist_os->SetBinError(bin, temp_points_os_e[bin-1]); + //temp_fit->SetBinContent(bin, temp_points_fit[bin-1]); + } + //cout<SetMarkerStyle(20); + temp_hist_ss->SetMarkerColor(1); + temp_hist_os->SetMarkerStyle(24); + temp_hist_os->SetMarkerColor(1); + temp_fit->SetLineColor(1); + C2_os->SetMarkerStyle(24); + C2_os->SetMarkerColor(2); + if(!MCcase){ + //C2_ss->SetMaximum(1.6); C2_ss->DrawCopy(); - legend1->AddEntry(C2_ss,"same-charge","p"); - //C2_os->DrawCopy("same"); - legend1->AddEntry(C2_os,"mixed-charge","p"); - + //legend1->AddEntry(C2_ss,"same-charge","p"); + C2_os->DrawCopy("same"); + //legend1->AddEntry(C2_os,"mixed-charge","p"); + /* Specif_System->Draw("same"); Specif_kT->Draw("same"); Specif_FSI->Draw("same"); @@ -1025,6 +1066,7 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll TLatex *Stats2 = new TLatex(.5,.45, "R_{inv} = 1.59 #pm 0.01 fm"); Stats2->SetNDC(); Stats2->SetTextSize(0.06); + */ //Stats2->Draw("same"); /*TF1 *BkgFitC2 = new TF1("BkgFitC2","1+[0]*exp(-pow([1]*x/0.19733,2))",0,1); BkgFitC2->SetParameter(0,0.08); @@ -1033,13 +1075,19 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll C2_ss->Fit(BkgFitC2,"IME","",0.2,0.8); BkgFitC2->Draw("same");*/ + //fitC2ss_h->SetLineColor(1); Unity->Draw("same"); - //fitC2ss_h->Draw("C same"); + fitC2ss_h->Draw("C same"); + //temp_hist_ss->Draw("same"); + //temp_hist_os->Draw("same"); + //temp_fit->Draw("C same"); //fitC2ss_Base->Draw("C same"); //fitC2ss_Expan->Draw("C same"); - //legend1->AddEntry(fitC2ss_Base,"Gaussian","l"); - //legend1->AddEntry(fitC2ss_Expan,"Edgeworth","l"); - //legend1->Draw("same"); + //legend1->AddEntry(C2_ss,"p-Pb (++)","p"); + //legend1->AddEntry(C2_os,"p-Pb (-+)","p"); + //legend1->AddEntry(temp_hist_ss,"Pb-Pb (++), Mbin 16","p"); + //legend1->AddEntry(temp_hist_os,"Pb-Pb (-+), =36","p"); + legend1->Draw("same"); } @@ -1295,7 +1343,8 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll A_3[i-1][j-1][k-1] = value_num + TERM5; B_3[i-1][j-1][k-1] = TERM5; A_3_e[i-1][j-1][k-1] = sqrt(value_num_e + TERM5); - + + //if(i==j && i==k && i==4) cout<SetParName(0,"N"); c3Fit1D_Base->SetParName(1,"#lambda_{3}"); c3Fit1D_Base->SetParName(2,"R_{inv}"); + if(!Gaussian) {maxVal_c3[1] = 5.0; maxVal_c3[2] = 20.;} for (int i=0; iDraw("same"); double lambdaStar=0, lambdaStar_e=0; if(Gaussian){ lambdaStar = OutputPar_c3[1]*pow(1 + OutputPar_c3[4]/8.,3); @@ -1570,7 +1626,7 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll EW23 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12); // double cumulant_fitvalue=0; - cumulant_fitvalue = c3Fit1D_Expan->GetParameter(0)*(c3Fit1D_Expan->GetParameter(1)*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5); + cumulant_fitvalue = c3Fit1D_Expan->GetParameter(0)*(c3Fit1D_Expan->GetParameter(1)*EW12*EW13*EW23*exp(-pow(radius_exp,ExpPower)/2. * (pow(Q12,ExpPower)+pow(Q13,ExpPower)+pow(Q23,ExpPower)))*TERM5 + TERM5); GenSignalExpected_num->Fill(Q3, cumulant_fitvalue); GenSignalExpected_den->Fill(Q3, TERM5); @@ -1594,6 +1650,8 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll bool splineOnly=kFALSE; for(int iii=0; iii<1000; iii++){ xpoints[iii] = 0 + (iii+0.5)*0.001; + if(!Gaussian && CollisionType!=0 && c3Fit1D_Expan->Eval(xpoints[iii]) < 2.2) splineOnly=kTRUE; + if(!Gaussian && CollisionType==0 && c3Fit1D_Expan->Eval(xpoints[iii]) < 2.0) splineOnly=kTRUE;// 2.0 or 1.4 for Mbin=3 in Fig 2 if(!splineOnly){ ypoints[iii] = c3Fit1D_Expan->Eval(xpoints[iii]); if(c3Fit1D_Expan->Eval(xpoints[iii])<2 && fabs(c3Fit1D_Expan->Eval(xpoints[iii])-c3Fit1D_ExpanSpline->Eval(xpoints[iii])) < 0.01) splineOnly=kTRUE; @@ -1605,9 +1663,27 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll if(CollisionType==0) gr_c3Spline->SetLineColor(1); if(CollisionType==1) gr_c3Spline->SetLineColor(2); if(CollisionType==2) gr_c3Spline->SetLineColor(4); + gr_c3Spline->Draw("c same"); - //c3Fit1D_ExpanSpline->Draw("same"); + double ChiSqSum_1D=0, SumNDF_1D=0; + for(int bin=1; bin<=300; bin++){ + double binCenter = c3_hist->GetXaxis()->GetBinCenter(bin); + if(binCenter > Q3Limit) continue; + if(c3_hist->GetBinError(bin)==0) continue; + int grIndex=1; + for(int gr=0; gr<999; gr++){ + if(binCenter > xpoints[gr] && (binCenter < xpoints[gr+1])) {grIndex=gr; break;} + } + + //ChiSqSum_1D += pow((c3_hist->GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2); + ChiSqSum_1D += pow((c3_hist->GetBinContent(bin)-c3Fit1D_Base->Eval(binCenter)) / c3_hist->GetBinError(bin),2); + //cout<GetBinContent(bin)<<" "<Eval(binCenter)<<" "<GetBinError(bin)<GetBinContent(bin)<<" "<GetBinError(bin)<GetParError(1))/100.; - *RName += round(100*c3Fit1D_Base->GetParameter(2))/100; RName->Append(" #pm "); *RName += int(100*c3Fit1D_Base->GetParError(2))/100.; RName->Append(" fm"); - TLatex *fitBox1=new TLatex(0.5,0.9,lamName->Data()); - fitBox1->SetNDC(kTRUE); - fitBox1->Draw(); - TLatex *fitBox2=new TLatex(0.5,0.8,RName->Data()); - fitBox2->SetNDC(kTRUE); - fitBox2->Draw();*/ + //TString *lamName=new TString("#lambda_{3}="); + //TString *RName=new TString("#R_{inv,3}="); + //*lamName += int(100*c3Fit1D_Base->GetParameter(1))/100.; lamName->Append(" #pm "); *lamName += int(100*c3Fit1D_Base->GetParError(1))/100.; + //*RName += round(100*c3Fit1D_Base->GetParameter(2))/100; RName->Append(" #pm "); *RName += int(100*c3Fit1D_Base->GetParError(2))/100.; RName->Append(" fm"); + //TLatex *fitBox1=new TLatex(0.5,0.9,lamName->Data()); + //fitBox1->SetNDC(kTRUE); + //fitBox1->Draw(); + //TLatex *fitBox2=new TLatex(0.5,0.8,"R_{inv,3} = 2.62 #pm 0.12"); + //fitBox2->SetNDC(kTRUE); + //fitBox2->Draw(); //TLatex *fitBox3=new TLatex(0.5,0.5,"=103"); //fitBox3->SetNDC(kTRUE); //fitBox3->Draw(); @@ -1705,7 +1781,8 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll //for(int ii=0; ii<10; ii++) cout<GetBinContent(ii+1)<<", "; //Coul_GRiverside->Draw(); //Coul_Riverside->Draw("same"); - /*TLegend *legend4 = new TLegend(.3,.8, .5,.95,NULL,"brNDC"); + /* + TLegend *legend4 = new TLegend(.3,.8, .5,.95,NULL,"brNDC"); legend4->SetBorderSize(0); legend4->SetTextFont(42);//42 legend4->SetTextSize(.04);// small .03; large .06 @@ -1751,12 +1828,13 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll C3_Extended->SetMarkerColor(4); c3_Extended->GetXaxis()->SetRangeUser(0,2); c3_Extended->SetMarkerStyle(20); - c3_Extended->Draw();*/ + c3_Extended->Draw(); + //cout<GetBinError(40)<<" "<GetBinError(40)/Three_1d[CB1][CB2][CB3][SCBin][4]->GetBinContent(40)<<" "<GetBinError(40)/Three_1d[CB1][CB2][CB3][SCBin][4]->GetBinContent(40)<AddEntry(c3_Extended,"c_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p"); //legend4->AddEntry(C3_Extended,"C_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p"); //legend4->Draw("same"); //Unity->Draw("same"); - + */ //Specif_System->Draw("same"); //Specif_KT3->Draw("same"); //Specif_FSI->Draw("same"); @@ -1824,12 +1902,12 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll K0sProjection->Draw("same"); K0sProjection_term1->Draw("same"); - legend4->AddEntry(K0sProjection,"#font[12]{C_{3}^{#pm#pm#mp}} projection","p"); - legend4->AddEntry(K0sProjection_term1,"#font[12]{#bf{c}_{3}^{#pm#pm#mp}} projection","p"); - legend4->Draw("same"); + //legend4->AddEntry(K0sProjection,"#font[12]{C_{3}^{#pm#pm#mp}} projection","p"); + //legend4->AddEntry(K0sProjection_term1,"#font[12]{#bf{c}_{3}^{#pm#pm#mp}} projection","p"); + //legend4->Draw("same"); + */ - */ ////////////////////////////////////////////////////////////////////////////////////// @@ -1841,6 +1919,7 @@ void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=Coll if(SaveToFile){ TString *savefilename = new TString("OutFiles/OutFile"); + if(!Gaussian) savefilename->Append("ExpFit"); if(CollisionType==0) savefilename->Append("PbPb"); else if(CollisionType==1) savefilename->Append("pPb"); else savefilename->Append("pp"); @@ -2140,6 +2219,7 @@ void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){ double SumChi2=0; //double lnL=0, term1=0, term2=0; NFitPoints_c3=0; + //double SumChi2_test=0; for(int i=0; i<=BINLIMIT_3; i++){// q12 for(int j=0; j<=BINLIMIT_3; j++){// q13 @@ -2175,6 +2255,7 @@ void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){ error += pow(sqrt(B_3[i][j][k])*A_3[i][j][k]/pow(B_3[i][j][k],2),2); error = sqrt(error); SumChi2 += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2); + //if(q3<0.05) SumChi2_test += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2); // /*if(A_3[i][j][k] >= 1) term1 = C*(A_3[i][j][k]+B_3[i][j][k])/(A_3[i][j][k]*(C+1)); else term1 = 0; @@ -2186,10 +2267,10 @@ void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){ lnL += B_3[i][j][k]*log(term2); }else {cout<<"WARNING -- term1 and term2 are negative"<Data(),"",3000,1,13.5);// Nch^(1/3) Parameters_c3[ct][ft][kt3][par] = new TH1D(name_c3->Data(),"",3000,1,13.5);// Nch^(1/3) + if(ft==0 && kt3==0 && par==0){ + TString *name_Bjoern = new TString("Bjoern_"); *name_Bjoern += ct; + Parameters_Bjoern[0][ct] = new TH1D(name_Bjoern->Data(),"",3000,1,13.5);// Nch^(1/3) + name_Bjoern->Append("_hydro"); + Parameters_Bjoern[1][ct] = new TH1D(name_Bjoern->Data(),"",3000,1,13.5);// Nch^(1/3) + } }else{ Parameters_C2[ct][ft][kt3][par] = new TH1D(name_C2->Data(),"",30000,1,3001);// Nch Parameters_c3[ct][ft][kt3][par] = new TH1D(name_c3->Data(),"",30000,1,3001);// Nch + if(ft==0 && kt3==0 && par==0){ + TString *name_Bjoern = new TString("Bjoern_"); *name_Bjoern += ct; + Parameters_Bjoern[0][ct] = new TH1D(name_Bjoern->Data(),"",30000,1,3001);// Nch + name_Bjoern->Append("_hydro"); + Parameters_Bjoern[1][ct] = new TH1D(name_Bjoern->Data(),"",30000,1,3001);// Nch + } } } } @@ -310,8 +325,24 @@ void Plot_plotsTPR(){ name3->Append(".root"); TFile *file=new TFile(name3->Data(),"READ"); + + TString *name3_2 = new TString("OutFiles/OutFileExpFit"); + if(ct==0) name3_2->Append("PbPb"); + if(ct==1) name3_2->Append("pPb"); + if(ct==2) name3_2->Append("pp"); + name3_2->Append("SC"); + if(ch==0) name3_2->Append("Neg"); + else name3_2->Append("Pos"); + + name3_2->Append("Kt3_"); + *name3_2 += KT3+1; + name3_2->Append("_M"); + *name3_2 += cb; + name3_2->Append(".root"); + + TFile *fileExpFit=new TFile(name3_2->Data(),"READ"); // - + if(ch==0) { C3[ct][dt][ChComb][KT3][cb]=(TH1D*)file->Get("C3"); @@ -347,6 +378,8 @@ void Plot_plotsTPR(){ C2[ct][KT3][cb]=(TH1D*)file->Get("C2_ss"); Fit_C2[ct][0][KT3][cb]=(TF1*)file->Get("fitC2ss_Base");// was fitC2ss_Gauss Fit_C2[ct][1][KT3][cb]=(TF1*)file->Get("fitC2ss_Expan");// fitC2ss_EW + Fit_C2[ct][2][KT3][cb]=(TF1*)fileExpFit->Get("fitC2ss_Base");// Exp fit + C2_Sys[ct][KT3][cb]=(TH1D*)C2[ct][KT3][cb]->Clone(); if(ct==0){ C2[ct][KT3][cb]->SetMarkerColor(1); C2[ct][KT3][cb]->SetLineColor(1); @@ -375,14 +408,17 @@ void Plot_plotsTPR(){ c3_fit[ct][0][KT3][cb]=(TF1*)file->Get("c3Fit1D_Base");// was c3Fit1D_Gauss c3_fit[ct][1][KT3][cb]=(TF1*)file->Get("c3Fit1D_Expan");// was c3Fit1D_EW - c3_fit[ct][0][KT3][cb]->SetLineStyle(2); - c3_fit[ct][1][KT3][cb]->SetLineStyle(1); - if(ct==0) {c3_fit[ct][0][KT3][cb]->SetLineColor(1); c3_fit[ct][1][KT3][cb]->SetLineColor(1);} - if(ct==1) {c3_fit[ct][0][KT3][cb]->SetLineColor(2); c3_fit[ct][1][KT3][cb]->SetLineColor(2);} - if(ct==2) {c3_fit[ct][0][KT3][cb]->SetLineColor(4); c3_fit[ct][1][KT3][cb]->SetLineColor(4);} + c3_fit[ct][2][KT3][cb]=(TF1*)fileExpFit->Get("c3Fit1D_Base"); gr_c3Spline[ct][KT3][cb] = (TGraph*)file->Get("gr_c3Spline");// Spline of a spline + TF1 + gr_c3SplineExpFit[ct][KT3][cb] = (TGraph*)fileExpFit->Get("gr_c3Spline");// Exp fit + c3_fit[ct][0][KT3][cb]->SetLineStyle(6); + gr_c3Spline[ct][KT3][cb]->SetLineStyle(7); c3_fit[ct][1][KT3][cb]->SetLineStyle(7); + gr_c3SplineExpFit[ct][KT3][cb]->SetLineStyle(1); c3_fit[ct][2][KT3][cb]->SetLineStyle(1); + if(ct==0) {c3_fit[ct][0][KT3][cb]->SetLineColor(1); c3_fit[ct][1][KT3][cb]->SetLineColor(1); c3_fit[ct][2][KT3][cb]->SetLineColor(1);} + if(ct==1) {c3_fit[ct][0][KT3][cb]->SetLineColor(2); c3_fit[ct][1][KT3][cb]->SetLineColor(2); c3_fit[ct][2][KT3][cb]->SetLineColor(2);} + if(ct==2) {c3_fit[ct][0][KT3][cb]->SetLineColor(4); c3_fit[ct][1][KT3][cb]->SetLineColor(4); c3_fit[ct][2][KT3][cb]->SetLineColor(4);} }// ChComb==0 - + } }else{ @@ -425,13 +461,15 @@ void Plot_plotsTPR(){ else if(ct==1) logNch=meanNchpPb[cb]; else logNch=meanNchpp[cb]; int logNchBin = Parameters_c3[ct][0][KT3][0]->GetXaxis()->FindBin(logNch); - for(int ft=0; ft<2; ft++){// Gaussian or EW + for(int ft=0; ft<3; ft++){// Gaussian or EW or Exp N_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(0); N_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(0); lambda_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(1); lambda_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(1); radius_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(2); radius_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(2); + if(ft==2) {radius_1 /= sqrt(TMath::Pi()); radius_1_e /= sqrt(TMath::Pi());} EW1_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(3); EW1_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(3); EW2_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(4); EW2_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(4); - if(ft==0) {EW1_1=0; EW2_1=0; EW1_1_e=0; EW2_1_e=0;}// make sure they are zero + + if(ft!=1) {EW1_1=0; EW2_1=0; EW1_1_e=0; EW2_1_e=0;}// make sure they are zero // Parameters_c3[ct][ft][KT3][0]->SetBinContent(logNchBin, N_1); Parameters_c3[ct][ft][KT3][0]->SetBinError(logNchBin, N_1_e); Parameters_c3[ct][ft][KT3][1]->SetBinContent(logNchBin, lambda_1); Parameters_c3[ct][ft][KT3][1]->SetBinError(logNchBin, lambda_1_e); @@ -459,15 +497,17 @@ void Plot_plotsTPR(){ Parameters_C2[ct][ft][KT3][0]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(0)); Parameters_C2[ct][ft][KT3][1]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(1));// lambda Parameters_C2[ct][ft][KT3][1]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(1)); - Parameters_C2[ct][ft][KT3][2]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(3));// R - Parameters_C2[ct][ft][KT3][2]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(3)); + radius_1 = Fit_C2[ct][ft][KT3][cb]->GetParameter(3); radius_1_e = Fit_C2[ct][ft][KT3][cb]->GetParError(3); + if(ft==2) {radius_1 /= sqrt(TMath::Pi()); radius_1_e /= sqrt(TMath::Pi());} + Parameters_C2[ct][ft][KT3][2]->SetBinContent(logNchBin, radius_1);// R + Parameters_C2[ct][ft][KT3][2]->SetBinError(logNchBin, radius_1_e); Parameters_C2[ct][ft][KT3][3]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(5));// kappa3 Parameters_C2[ct][ft][KT3][3]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(5)); Parameters_C2[ct][ft][KT3][4]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(6));// kappa4 Parameters_C2[ct][ft][KT3][4]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(6)); // lambda_* parameter - Parameters_C2[ct][ft][KT3][5]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(1)*pow(1 + Fit_C2[ct][ft][KT3][cb]->GetParameter(6)/8.,2)); - Parameters_C2[ct][ft][KT3][5]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(1)*pow(1 + Fit_C2[ct][ft][KT3][cb]->GetParameter(6)/8.,2)); + Parameters_C2[ct][ft][KT3][5]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(1)*pow(1 + EW2_1/8.,2)); + Parameters_C2[ct][ft][KT3][5]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(1)*pow(1 + EW2_1/8.,2)); }// ft }// ChComb==0 // @@ -503,8 +543,9 @@ void Plot_plotsTPR(){ } C3_Sys[ct][0][ChComb][KT3][cb]->SetDirectory(0); c3_Sys[ct][0][ChComb][KT3][cb]->SetDirectory(0); }// AddedCC and dt==0 + file->Close(); - + fileExpFit->Close(); }// Kt3 @@ -514,7 +555,7 @@ void Plot_plotsTPR(){ }// cb if(dt==0){ - for(int ft=0; ft<2; ft++){// Gaussian or EW + for(int ft=0; ft<3; ft++){// Gaussian or EW or Exp for(int KT3=0; KT3<3; KT3++){ for(int par=0; par<6; par++){ if(ct<2){ @@ -549,7 +590,8 @@ void Plot_plotsTPR(){ if(par==1) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#lambda_{e} or #lambda_{e,3}"); if(par==2) { if(FitType==0) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R}_{inv,2} or #font[12]{R}_{inv,3} (fm)"); - else Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R^{E_{w}}}_{inv,2} or #font[12]{R^{E_{w}}}_{inv,3} (fm)"); + else if(FitType==1) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R^{E_{w}}}_{inv,2} or #font[12]{R^{E_{w}}}_{inv,3} (fm)"); + else Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R}^{Exp}_{inv,2} or #font[12]{R}^{Exp}_{inv,3} (fm)"); } if(par==3) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#kappa_{3}"); if(par==4) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#kappa_{4}"); @@ -571,8 +613,63 @@ void Plot_plotsTPR(){ Unity->SetLineColor(1); - //return; - + ///////////////////////////////////// + // Bjoern's points + double Bjoern_xaxis_pp[4] = {4.1, 9.2, 16, 26};// Nch + double Bjoern_Ri_pp[2][4] = {{0.87, 1.09, 1.22, 1.33},{0.837, 0.97, 1.03, 1.18}};// m01 or m02 (infrared cutoff of calc), Nch + double Bjoern_Rhydro_pp[2][4] = {{0.87, 1.11, 1.35, 1.57},{0.845, 1.04, 1.22, 1.56}};// m01 or m02 (infrared cutoff of calc), Nch + // + double Bjoern_xaxis_pPb[6] = {5, 11, 20, 32, 45, 69};// Nch + double Bjoern_Ri_pPb[2][6] = {{0.98, 1.27, 1.46, 1.59, 1.75, 2},{0.96, 1.18, 1.28, 1.36, 1.44, 1.56}}; + double Bjoern_Rhydro_pPb[2][6] = {{0.98, 1.28, 1.62, 1.96, 2.32, 2.79},{0.96, 1.22, 1.5, 1.76, 2.09, 2.58}}; + // + double Bjoern_xaxis_PbPb[10] = {22, 36, 50, 77, 98, 137, 172, 326, 498, 760};// Nch + double Bjoern_Ri_PbPb[2][10] = {{1.6, 2.05, 2.47, 2.86, 3.17, 3.52, 3.88, 4.54, 5.19, 5.85},{1.57, 1.87, 2.27, 2.65, 3, 3.3, 3.52, 4.17, 4.82, 5.49}}; + double Bjoern_Rhydro_PbPb[2][10] = {{1.6, 2.06, 2.49, 2.88, 3.27, 3.63, 4.03, 4.91, 5.88, 6.88},{1.57, 1.88, 2.32, 2.75, 3.14, 3.53, 3.84, 4.65, 5.52, 6.4}}; + // + int mchoice = 1;// 1 or 2 indicating Bjoern's cutoff + double SF_Bjoern = 1.15;// 1.15 (m=0.1) + // + + + for(int index=0; index<4; index++){// pp + for(int m=0; m<2; m++) {Bjoern_Ri_pp[m][index] *= SF_Bjoern; Bjoern_Rhydro_pp[m][index] *= SF_Bjoern;} + double xpoint = Bjoern_xaxis_pp[index]; + if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_pp[index], 1/3.); + int bin = Parameters_Bjoern[0][2]->GetXaxis()->FindBin(xpoint); + Parameters_Bjoern[0][2]->SetBinContent(bin, Bjoern_Ri_pp[mchoice-1][index]); + Parameters_Bjoern[0][2]->SetBinError(bin, 0.001); + Parameters_Bjoern[1][2]->SetBinContent(bin, Bjoern_Rhydro_pp[mchoice-1][index]); + Parameters_Bjoern[1][2]->SetBinError(bin, 0.001); + } + for(int index=0; index<6; index++){// pPb + for(int m=0; m<2; m++) {Bjoern_Ri_pPb[m][index] *= SF_Bjoern; Bjoern_Rhydro_pPb[m][index] *= SF_Bjoern;} + double xpoint = Bjoern_xaxis_pPb[index]; + if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_pPb[index], 1/3.); + int bin = Parameters_Bjoern[0][1]->GetXaxis()->FindBin(xpoint); + Parameters_Bjoern[0][1]->SetBinContent(bin, Bjoern_Ri_pPb[mchoice-1][index]); + Parameters_Bjoern[0][1]->SetBinError(bin, 0.001); + Parameters_Bjoern[1][1]->SetBinContent(bin, Bjoern_Rhydro_pPb[mchoice-1][index]); + Parameters_Bjoern[1][1]->SetBinError(bin, 0.001); + } + for(int index=0; index<10; index++){// PbPb + for(int m=0; m<2; m++) {Bjoern_Ri_PbPb[m][index] *= SF_Bjoern; Bjoern_Rhydro_PbPb[m][index] *= SF_Bjoern;} + double xpoint = Bjoern_xaxis_PbPb[index]; + if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_PbPb[index], 1/3.); + int bin = Parameters_Bjoern[0][0]->GetXaxis()->FindBin(xpoint); + Parameters_Bjoern[0][0]->SetBinContent(bin, Bjoern_Ri_PbPb[mchoice-1][index]); + Parameters_Bjoern[0][0]->SetBinError(bin, 0.001); + Parameters_Bjoern[1][0]->SetBinContent(bin, Bjoern_Rhydro_PbPb[mchoice-1][index]); + Parameters_Bjoern[1][0]->SetBinError(bin, 0.001); + } + Parameters_Bjoern[0][0]->SetMarkerColor(1); Parameters_Bjoern[0][0]->SetMarkerStyle(20); + Parameters_Bjoern[0][1]->SetMarkerColor(2); Parameters_Bjoern[0][1]->SetMarkerStyle(21); + Parameters_Bjoern[0][2]->SetMarkerColor(4); Parameters_Bjoern[0][2]->SetMarkerStyle(34); + // + Parameters_Bjoern[1][0]->SetMarkerColor(1); Parameters_Bjoern[1][0]->SetMarkerStyle(20); + Parameters_Bjoern[1][1]->SetMarkerColor(2); Parameters_Bjoern[1][1]->SetMarkerStyle(21); + Parameters_Bjoern[1][2]->SetMarkerColor(4); Parameters_Bjoern[1][2]->SetMarkerStyle(34); + //////////////////////////////////////////////////// //////////////////////////////////////////////////// // Progaganda Plot @@ -738,10 +835,14 @@ void Plot_plotsTPR(){ Sys_K0s_C3->SetMarkerSize(0); Sys_K0s_c3->SetMarkerSize(0); Sys_K0s_C3->SetFillColor(kBlue-10); Sys_K0s_c3->SetFillColor(kBlue-10); Sys_K0s_C3->SetFillStyle(1000); Sys_K0s_c3->SetFillStyle(1000); + cout.precision(3); for(int i=0; iGetNbinsX(); i++) { Sys_K0s_C3->SetBinError(i+1, 0.01 * Sys_K0s_C3->GetBinContent(i+1)); Sys_K0s_c3->SetBinError(i+1, sqrt(pow(0.01 * Sys_K0s_c3->GetBinContent(i+1),2) + pow(0.1*(Sys_K0s_c3->GetBinContent(i+1)-1),2))); + //cout<GetXaxis()->GetBinLowEdge(i+1)<<" "<GetXaxis()->GetBinUpEdge(i+1)<<" "<GetBinContent(i+1)<<" "<GetBinError(i+1)<<" "<GetBinError(i+1)<GetXaxis()->GetBinLowEdge(i+1)<<" "<GetXaxis()->GetBinUpEdge(i+1)<<" "<GetBinContent(i+1)<<" "<GetBinError(i+1)<<" "<GetBinError(i+1)<Draw("E2 same"); Sys_K0s_c3->Draw("E2 same"); K0s_C3->Draw("same"); @@ -773,7 +874,9 @@ void Plot_plotsTPR(){ float RadiiC2ppPubPoints_e[2][8]={{0.058,0.064,0.071,0.071,0.078,0.078,0.086,0.098},{0.12,0.12,0.13,0.13,0.13,0.13,0.13,0.13}}; float MeanPubNch[8]={3.36, 7.92, 11.04, 14.4, 17.88, 21.48, 25.68, 33.12};// Adam's * 1.6 * 0.75. Factor of 0.75 for low,high pt extrap - TCanvas *can3 = new TCanvas("can3", "can3",1000,0,600,900);// 11,53,700,500 + int yExtent = 900; + if(RadiusOnly) yExtent = 600; + TCanvas *can3 = new TCanvas("can3", "can3",1000,0,600,yExtent);// 1000,0,600,900 can3->SetHighLightColor(2); gStyle->SetOptFit(0111); can3->SetFillColor(0);//10 @@ -791,7 +894,7 @@ void Plot_plotsTPR(){ pad3->SetRightMargin(0.0);//3e-2 pad3->SetBottomMargin(0.0);//0.12 pad3->SetLeftMargin(0.0); - pad3->Divide(1,2,0,0); + if(!RadiusOnly) pad3->Divide(1,2,0,0); pad3->Draw(); pad3->cd(); TLegend *legend4 = new TLegend(.15,.65, .55,.95,NULL,"brNDC");//.45 or .4 for x1 @@ -805,7 +908,11 @@ void Plot_plotsTPR(){ gPad->SetRightMargin(0.01); gPad->SetTopMargin(0.01); gPad->SetBottomMargin(0.001);// 0.16 - + if(RadiusOnly){ + gPad->SetLeftMargin(0.16); + gPad->SetBottomMargin(0.16); + } + gPad->SetTickx(); gPad->SetTicky(); if(!NchOneThirdAxis)gPad->SetLogx(); @@ -823,11 +930,24 @@ void Plot_plotsTPR(){ //gStyle->SetErrorX(0); if(NchOneThirdAxis) RadiiPbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000 else RadiiPbPb->GetXaxis()->SetRangeUser(3,3000);// 3,3000 + if(RadiusOnly){ + RadiiPbPb->GetXaxis()->SetTitleSize(SizeTitle); + RadiiPbPb->GetXaxis()->SetTitleOffset(1.25); + RadiiPbPb->GetYaxis()->SetTitleSize(SizeTitle); + RadiiPbPb->GetYaxis()->SetTitleOffset(1.2); + if(FitType==0) RadiiPbPb->GetYaxis()->SetTitle("#font[12]{R}^{#font[12]{G}}_{inv} or #font[12]{R}^{#font[12]{G}}_{inv,3} (fm)"); + else if(FitType==1) RadiiPbPb->GetYaxis()->SetTitle("#font[12]{R}^{#font[12]{E}_{w}}_{inv} or #font[12]{R}^{#font[12]{E}_{w}}_{inv,3} (fm)"); + else RadiiPbPb->GetYaxis()->SetTitle("#font[12]{R}^{#font[12]{Exp}_{inv} or #font[12]{R}^{#font[12]{Exp}_{inv,3} (fm)"); + } + RadiiPbPb->Draw(); + // + double xAxisC2_e[20]={0}; double xAxis_e[20]={0}; double yAxisPbPb[20]={0}; - double yAxisPbPb_e[20]={0}; + double yAxisPbPb_eL[20]={0}; + double yAxisPbPb_eH[20]={0}; double yAxispPb[20]={0}; double yAxispPb_eL[20]={0}; double yAxispPb_eH[20]={0}; @@ -835,24 +955,27 @@ void Plot_plotsTPR(){ double yAxispp_eL[20]={0}; double yAxispp_eH[20]={0}; - double SysPercent_PbPb[2][4]={{12,10,11,16},{5,7,5,10}};// Gauss/EW, parameter#(Rinv2, Rinv3, lambda, lambda_3) + double SysPercent_PbPb_NFR[3][4]={{12,10,11,16}, {5,7,5,10}, {1,3,1,5}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) // Narrow Fit Range for C2 fits - double SysPercent_pPb_NFR[2][4]={{15,10,6,10},{6,2,1,4}};// Gauss/EW, parameter#(Rinv2, Rinv3, lambda, lambda_3) Old values with 30% fit range variation - double SysPercent_pp_NFR[2][4]={{11,9,2,5},{12,5,2,5}};// Gauss/EW, parameter#(Rinv2, Rinv3, lambda, lambda_3) Old values with 30% fit range variation + double SysPercent_pPb_NFR[3][4]={{15,10,6,10}, {6,2,1,4}, {5,3,1,4}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) Old values with 30% fit range variation + double SysPercent_pp_NFR[3][4]={{11,9,2,5}, {12,5,2,5}, {4,5,1,7}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) Old values with 30% fit range variation // Wide Fit Range for C2 fits - double SysPercent_pPb_WFR[2][4]={{24,10,6,10},{16,2,7,4}};// Gauss/EW, parameter#(Rinv2, Rinv3, lambda, lambda_3) New values with 1.2 GeV/c fit range variation - double SysPercent_pp_WFR[2][4]={{21,9,2,5},{6,5,1,5}};// Gauss/EW, parameter#(Rinv2, Rinv3, lambda, lambda_3) New values with 1.2 GeV/c fit range variation + double SysPercent_PbPb_WFR[3][4]={{25,10,12,16}, {8,7,3,10}, {10,3,6,5}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) with pPb fit range + double SysPercent_pPb_WFR[3][4]={{24,10,6,10}, {16,2,7,4}, {15,3,9,4}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) New values with 1.2 GeV/c fit range variation + double SysPercent_pp_WFR[3][4]={{21,9,2,5}, {6,5,1,5}, {2,5,1,7}};// Gauss/EW/Exp, parameter#(Rinv2, Rinv3, lambda, lambda_3) New values with 1.2 GeV/c fit range variation for(int cb=0; cb<20; cb++){// 3-particle int binPbPb = RadiiPbPb->GetXaxis()->FindBin(meanNchPbPb[cb]); int binpPb = RadiipPb->GetXaxis()->FindBin(meanNchpPb[cb]); int binpp = Radiipp->GetXaxis()->FindBin(meanNchpp[cb]); - - double RsysPbPb = 0.01*sqrt(pow(SysPercent_PbPb[FitType][1],2) + pow(1,2)) * RadiiPbPb->GetBinContent(binPbPb);// fit Variation + MRC + double RsysPbPb = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][1],2) + pow(1,2)) * RadiiPbPb->GetBinContent(binPbPb);// fit Variation + MRC double RsyspPb = 0.01*sqrt(pow(SysPercent_pPb_NFR[FitType][1],2) + pow(1,2)) * RadiipPb->GetBinContent(binpPb);// fit Variation + MRC double Rsyspp = 0.01*sqrt(pow(SysPercent_pp_NFR[FitType][1],2) + pow(1,2)) * Radiipp->GetBinContent(binpp);// fit Variation + MRC - if(RadiiPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_e[cb]=100;} - else {yAxisPbPb[cb]=RadiiPbPb->GetBinContent(binPbPb); yAxisPbPb_e[cb]=RsysPbPb;} + // + if(cb==19 && FitType==2) Rsyspp = 0.01*sqrt(pow(16,2) + pow(1,2)) * Radiipp->GetBinContent(binpp);// larger fit var in this bin + + if(RadiiPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100 + else {yAxisPbPb[cb]=RadiiPbPb->GetBinContent(binPbPb); yAxisPbPb_eL[cb]=RsysPbPb; yAxisPbPb_eH[cb]=RsysPbPb;} // if(RadiipPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;} else {yAxispPb[cb]=RadiipPb->GetBinContent(binpPb); yAxispPb_eL[cb]=RsyspPb; yAxispPb_eH[cb]=RsyspPb;} @@ -871,7 +994,7 @@ void Plot_plotsTPR(){ } - TGraphErrors *grRadiiSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_e); + TGraphErrors *grRadiiSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_eL); TGraphAsymmErrors *grRadiiSys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH); TGraphAsymmErrors *grRadiiSys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH); grRadiiSys_pp->SetMarkerSize(0); grRadiiSys_pp->SetFillStyle(1000); grRadiiSys_pp->SetFillColor(kBlue-10); @@ -895,40 +1018,51 @@ void Plot_plotsTPR(){ int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]); int binpPb = RadiiC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]); int binpp = RadiiC2pp->GetXaxis()->FindBin(meanNchpp[cb]); - double RsysPbPb = 0.01*sqrt(pow(SysPercent_PbPb[FitType][0],2) + pow(1,2)) * RadiiC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC + double RsysPbPb_L = 0.01*sqrt(pow(SysPercent_PbPb_WFR[FitType][0],2) + pow(1,2)) * RadiiC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC + double RsysPbPb_H = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][0],2) + pow(1,2)) * RadiiC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC double RsyspPb_L = 0.01*sqrt(pow(SysPercent_pPb_WFR[FitType][0],2) + pow(1,2)) * RadiiC2pPb->GetBinContent(binpPb);// fit Variation + MRC double RsyspPb_H = 0.01*sqrt(pow(SysPercent_pPb_NFR[FitType][0],2) + pow(1,2)) * RadiiC2pPb->GetBinContent(binpPb);// fit Variation + MRC double Rsyspp_L = 0.01*sqrt(pow(SysPercent_pp_WFR[FitType][0],2) + pow(1,2)) * RadiiC2pp->GetBinContent(binpp);// fit Variation + MRC double Rsyspp_H = 0.01*sqrt(pow(SysPercent_pp_NFR[FitType][0],2) + pow(1,2)) * RadiiC2pp->GetBinContent(binpp);// fit Variation + MRC - if(RadiiC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_e[cb]=1000;} - else {yAxisPbPb[cb]=RadiiC2PbPb->GetBinContent(binPbPb); yAxisPbPb_e[cb]=RsysPbPb;} + // + if(cb==15) RsysPbPb_L = 0.01*sqrt(pow(1.5*SysPercent_PbPb_WFR[FitType][0],2) + pow(1,2)) * RadiiC2PbPb->GetBinContent(binPbPb);// larger error + // + + if(RadiiC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=1000; yAxisPbPb_eH[cb]=1000;}// errors were 1000 + else { + yAxisPbPb[cb]=RadiiC2PbPb->GetBinContent(binPbPb); + if(cb>=13) yAxisPbPb_eL[cb]=RsysPbPb_L; + else yAxisPbPb_eL[cb]=RsysPbPb_H; + yAxisPbPb_eH[cb]=RsysPbPb_H; + } if(RadiiC2pPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=1000; yAxispPb_eH[cb]=1000;} else {yAxispPb[cb]=RadiiC2pPb->GetBinContent(binpPb); yAxispPb_eL[cb]=RsyspPb_L; yAxispPb_eH[cb]=RsyspPb_H;} if(RadiiC2pp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=1000; yAxispp_eH[cb]=1000;} else {yAxispp[cb]=RadiiC2pp->GetBinContent(binpp); yAxispp_eL[cb]=Rsyspp_L; yAxispp_eH[cb]=Rsyspp_H;} - //xAxis_e[cb]=10000; + // } - TGraphErrors *grRadiiC2Sys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_e); - TGraphAsymmErrors *grRadiiC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH); - TGraphAsymmErrors *grRadiiC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH); + + TGraphAsymmErrors *grRadiiC2Sys_PbPb=new TGraphAsymmErrors(20,meanNchPbPb,yAxisPbPb, xAxisC2_e,xAxisC2_e, yAxisPbPb_eL,yAxisPbPb_eH); + TGraphAsymmErrors *grRadiiC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxisC2_e,xAxisC2_e, yAxispPb_eL,yAxispPb_eH); + TGraphAsymmErrors *grRadiiC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxisC2_e,xAxisC2_e, yAxispp_eL,yAxispp_eH); grRadiiC2Sys_pp->SetMarkerSize(0); grRadiiC2Sys_pPb->SetMarkerSize(0); grRadiiC2Sys_PbPb->SetMarkerSize(0); grRadiiC2Sys_pp->SetLineColor(4); grRadiiC2Sys_pPb->SetLineColor(2); grRadiiC2Sys_PbPb->SetLineColor(1); grRadiiC2Sys_pp->SetLineWidth(2.*grRadiiC2Sys_pp->GetLineWidth()); grRadiiC2Sys_pPb->SetLineWidth(2.*grRadiiC2Sys_pPb->GetLineWidth()); - grRadiiC2Sys_PbPb->SetLineWidth(2.*grRadiiC2Sys_PbPb->GetLineWidth()); + grRadiiC2Sys_PbPb->SetLineWidth(1.*grRadiiC2Sys_PbPb->GetLineWidth()); grRadiiC2Sys_PbPb->SetLineStyle(2); // grRadiiC2Sys_pPb->Draw("|| p"); grRadiiC2Sys_pp->Draw("|| p"); - + grRadiiSys_pp->Draw("E2 p"); grRadiiSys_pPb->Draw("E2 p"); grRadiiSys_PbPb->Draw("E2 p"); RadiiPbPb->Draw("same"); RadiipPb->Draw("same"); Radiipp->Draw("same"); - grRadiiC2Sys_PbPb->Draw("|| p");// E2 or || to visualize pol2 fit below + grRadiiC2Sys_PbPb->Draw("E");// E2 or "|| p" to visualize pol2 fit below //if(FitType==0) RadiiC2pp_Published->Draw("same"); // @@ -936,7 +1070,7 @@ void Plot_plotsTPR(){ RadiiC2pPb->Draw("same"); RadiiC2pp->Draw("same"); - + legend4->AddEntry(Radiipp,"pp #sqrt{s}=7 TeV","p"); legend4->AddEntry(RadiipPb,"p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV","p"); legend4->AddEntry(RadiiPbPb,"Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV","p"); @@ -975,11 +1109,13 @@ void Plot_plotsTPR(){ Specif_Kt3->SetTextFont(TextFont); Specif_kt->SetTextFont(TextFont); Specif_Kt3->SetTextSize(SizeSpecif*SF_2panel); Specif_kt->SetTextSize(SizeSpecif*SF_2panel); Specif_Kt3->SetNDC(); Specif_kt->SetNDC(); - Specif_Kt3->Draw("same"); - Specif_kt->Draw("same"); - + if(!RadiusOnly){ + Specif_Kt3->Draw("same"); + Specif_kt->Draw("same"); + } legend4->SetTextFont(TextFont); legend4->SetTextSize(SizeLegend*SF_2panel); + if(RadiusOnly) legend4->SetTextSize(SizeLegend*SF_2panel*0.8); legend4->Draw("same"); TH1D *MarkerPbPb_3=(TH1D*)RadiiPbPb->Clone(); @@ -995,39 +1131,45 @@ void Plot_plotsTPR(){ MarkerPbPb_2->SetBinError(i,0.001); MarkerpPb_2->SetBinError(i,0.001); Markerpp_2->SetBinError(i,0.001); } if(!NchOneThirdAxis){ - MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(450), 1.25);// 1 - MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(600), 1.25);// 1 - Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(800), 1.25);// 1 + MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(450), 1.25);// 1.25, 1.45 for RadiusOnly + MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(600), 1.25);// 1.25, 1.45 for RadiusOnly + Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(800), 1.25);// 1.25, 1.45 for RadiusOnly MarkerPbPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(450), 3.1);// 3.1 MarkerpPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(600), 3.1);// 3.1 Markerpp_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(800), 3.1);// 3.1 }else{ - MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 1.25);// - MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 1.25);// - Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 1.25);// - MarkerPbPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 3.1);// - MarkerpPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 3.1);// - Markerpp_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 3.1);// + MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 1.25);// 1.25 + MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 1.25);// 1.25 + Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 1.25);// 1.25 + MarkerPbPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 3.1);// 3.1 + MarkerpPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 3.1);// 3.1 + Markerpp_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 3.1);// 3.1 } MarkerPbPb_3->Draw("same"); MarkerpPb_3->Draw("same"); Markerpp_3->Draw("same"); MarkerPbPb_2->Draw("same"); MarkerpPb_2->Draw("same"); Markerpp_2->Draw("same"); - TLatex *TwoPionText = new TLatex(0.74,0.3,"Two-Pions");// 0.74,0.41 - TLatex *ThreePionText = new TLatex(0.74,0.15,"Three-Pions");// 0.74,0.265 + TLatex *TwoPionText; + if(!RadiusOnly) TwoPionText = new TLatex(0.74,0.3,"Two-Pions");// 0.74,0.3 + else TwoPionText = new TLatex(0.67,0.31,"Two-Pions");// + TLatex *ThreePionText; + if(!RadiusOnly) ThreePionText = new TLatex(0.74,0.15,"Three-Pions");// 0.74,0.15 + else ThreePionText = new TLatex(0.65,0.2,"Three-Pions");// TwoPionText->SetNDC(); ThreePionText->SetNDC(); TwoPionText->SetTextFont(TextFont); ThreePionText->SetTextFont(TextFont); TwoPionText->SetTextSize(SizeSpecif*SF_2panel); ThreePionText->SetTextSize(SizeSpecif*SF_2panel); TwoPionText->Draw("same"); ThreePionText->Draw("same"); - + TLatex *Specif_kappas = new TLatex(0.42,0.05,"#kappa_{3}=0.1, #kappa_{4}=0.5");// 0.42,0.2 //TLatex *Specif_kappas = new TLatex(0.42,0.05,"#kappa_{3}(N_{ch}), #kappa_{4}=0.5");// 0.42,0.2 Specif_kappas->SetNDC(); Specif_kappas->SetTextFont(TextFont); Specif_kappas->SetTextSize(SizeSpecif*SF_2panel); - if(FitType==1) Specif_kappas->Draw("same"); + if(FitType==1 && !RadiusOnly) Specif_kappas->Draw("same"); + if(RadiusOnly) return; + /////////////////////////////////////////////////////////////////// pad3->cd(2); gPad->SetLeftMargin(0.14); @@ -1048,6 +1190,7 @@ void Plot_plotsTPR(){ LambdaPbPb->GetXaxis()->SetTitleFont(TextFont); LambdaPbPb->GetXaxis()->SetTitleSize(SizeTitle*SF_2panel); LambdaPbPb->GetYaxis()->SetTitleFont(TextFont); LambdaPbPb->GetYaxis()->SetTitleSize(SizeTitle*SF_2panel); LambdaPbPb->SetMaximum(2.8);// 2.8 + if(FitType==2) LambdaPbPb->SetMaximum(5.8); LambdaPbPb->GetXaxis()->SetTitleOffset(0.95); LambdaPbPb->GetYaxis()->SetTitleOffset(100);//1.1 if(NchOneThirdAxis) LambdaPbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000 @@ -1069,18 +1212,20 @@ void Plot_plotsTPR(){ if(cb<=12) f_syst_PbPb = 100 * (c3_mixedChargeSysFit[0][KT3Bin][cb]->Eval(0.025)-1.0) / (c3_fit[0][1][KT3Bin][cb]->Eval(0.025)-1.0);// residue / EW fit at Q3=0.025 if(cb>=12 && cb<19) f_syst_pPb = 100 * (c3_mixedChargeSysFit[1][KT3Bin][cb]->Eval(0.075)-1.0) / (c3_fit[1][1][KT3Bin][cb]->Eval(0.075)-1.0);// residue / EW fit at Q3=0.075 if(cb>=14) f_syst_pp = 100 * (c3_mixedChargeSysFit[2][KT3Bin][cb]->Eval(0.075)-1.0) / (c3_fit[2][1][KT3Bin][cb]->Eval(0.075)-1.0);// residue / EW fit at Q3=0.075 - double LambdasysPbPb = 0.01*sqrt(pow(SysPercent_PbPb[FitType][3],2) + pow(1,2) + pow(5,2) + pow(10,2) + pow(f_syst_PbPb,2)) * LambdaPbPb->GetBinContent(binPbPb);// fit Variation + MRC + TTC + undilution + f1,f2,f3 uncertainties + double LambdasysPbPb = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][3],2) + pow(1,2) + pow(5,2) + pow(10,2) + pow(f_syst_PbPb,2)) * LambdaPbPb->GetBinContent(binPbPb);// fit Variation + MRC + TTC + undilution + f1,f2,f3 uncertainties double LambdasyspPb = 0.01*sqrt(pow(SysPercent_pPb_WFR[FitType][3],2) + pow(1,2) + pow(10,2) + pow(f_syst_pPb,2)) * LambdapPb->GetBinContent(binpPb);// fit Variation + MRC + undilution + f1,f2,f3 uncertainties double Lambdasyspp = 0.01*sqrt(pow(SysPercent_pp_WFR[FitType][3],2) + pow(1,2) + pow(10,2) + pow(f_syst_pp,2)) * Lambdapp->GetBinContent(binpp);// fit Variation + MRC + undilution + f1,f2,f3 uncertainties - if(LambdaPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_e[cb]=100;} - else {yAxisPbPb[cb]=LambdaPbPb->GetBinContent(binPbPb); yAxisPbPb_e[cb]=LambdasysPbPb;} + if(cb==19 && FitType==2) Lambdasyspp = 0.01*sqrt(pow(25,2) + pow(1,2) + pow(10,2) + pow(f_syst_pp,2)) * Lambdapp->GetBinContent(binpp);// larger fit var in this bin + if(LambdaPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100 + else {yAxisPbPb[cb]=LambdaPbPb->GetBinContent(binPbPb); yAxisPbPb_eL[cb]=LambdasysPbPb; yAxisPbPb_eH[cb]=LambdasysPbPb;} // if(LambdapPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;} else {yAxispPb[cb]=LambdapPb->GetBinContent(binpPb); yAxispPb_eL[cb]=LambdasyspPb; yAxispPb_eH[cb]=LambdasyspPb;} // if(Lambdapp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=100; yAxispp_eH[cb]=100;} else {yAxispp[cb]=Lambdapp->GetBinContent(binpp); yAxispp_eL[cb]=Lambdasyspp; yAxispp_eH[cb]=Lambdasyspp;} - + + if(NchOneThirdAxis) { if(cb<13) xAxis_e[cb]=fabs(meanNchPbPb[cb] - pow(0.95,1/3.)*meanNchPbPb[cb]); else xAxis_e[cb]=fabs(meanNchpPb[cb] - pow(0.95,1/3.)*meanNchpPb[cb]); @@ -1089,7 +1234,8 @@ void Plot_plotsTPR(){ else xAxis_e[cb]=fabs(meanNchpPb[cb] - (0.95)*meanNchpPb[cb]); } } - TGraphErrors *grLambdaSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_e); + + TGraphErrors *grLambdaSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_eL); TGraphAsymmErrors *grLambdaSys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH); TGraphAsymmErrors *grLambdaSys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH); grLambdaSys_pp->SetMarkerSize(0); grLambdaSys_pp->SetFillStyle(1000); grLambdaSys_pp->SetFillColor(kBlue-10); @@ -1104,13 +1250,20 @@ void Plot_plotsTPR(){ int binPbPb = LambdaC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]); int binpPb = LambdaC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]); int binpp = LambdaC2pp->GetXaxis()->FindBin(meanNchpp[cb]); - double LambdasysPbPb = 0.01*sqrt(pow(SysPercent_PbPb[FitType][2],2) + pow(1,2) + pow(5,2) + pow(7,2)) * LambdaC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC + TTC + undilution + double LambdasysPbPb_H = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][2],2) + pow(1,2) + pow(5,2) + pow(7,2)) * LambdaC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC + TTC + undilution + double LambdasysPbPb_L = 0.01*sqrt(pow(SysPercent_PbPb_WFR[FitType][2],2) + pow(1,2) + pow(5,2) + pow(7,2)) * LambdaC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC + TTC + undilution double LambdasyspPb_H = 0.01*sqrt(pow(SysPercent_pPb_NFR[FitType][2],2) + pow(1,2) + pow(7,2)) * LambdaC2pPb->GetBinContent(binpPb);// fit Variation + MRC + undilution double LambdasyspPb_L = 0.01*sqrt(pow(SysPercent_pPb_WFR[FitType][2],2) + pow(1,2) + pow(7,2)) * LambdaC2pPb->GetBinContent(binpPb);// fit Variation + MRC + undilution double Lambdasyspp_H = 0.01*sqrt(pow(SysPercent_pp_NFR[FitType][2],2) + pow(1,2) + pow(7,2)) * LambdaC2pp->GetBinContent(binpp);// fit Variation + MRC + undilution double Lambdasyspp_L = 0.01*sqrt(pow(SysPercent_pp_WFR[FitType][2],2) + pow(1,2) + pow(7,2)) * LambdaC2pp->GetBinContent(binpp);// fit Variation + MRC + undilution - if(LambdaC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_e[cb]=100;} - else {yAxisPbPb[cb]=LambdaC2PbPb->GetBinContent(binPbPb); yAxisPbPb_e[cb]=LambdasysPbPb;} + // + if(LambdaC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100 + else { + yAxisPbPb[cb]=LambdaC2PbPb->GetBinContent(binPbPb); + if(cb>=13) yAxisPbPb_eL[cb]=LambdasysPbPb_L; + else yAxisPbPb_eL[cb]=LambdasysPbPb_H; + yAxisPbPb_eH[cb]=LambdasysPbPb_H; + } // if(LambdaC2pPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;} else {yAxispPb[cb]=LambdaC2pPb->GetBinContent(binpPb); yAxispPb_eL[cb]=LambdasyspPb_L; yAxispPb_eH[cb]=LambdasyspPb_H;} @@ -1119,9 +1272,9 @@ void Plot_plotsTPR(){ else {yAxispp[cb]=LambdaC2pp->GetBinContent(binpp); yAxispp_eL[cb]=Lambdasyspp_L; yAxispp_eH[cb]=Lambdasyspp_H;} //xAxis_e[cb]=10000; } - TGraphErrors *grLambdaC2Sys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_e); - TGraphAsymmErrors *grLambdaC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH); - TGraphAsymmErrors *grLambdaC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH); + TGraphAsymmErrors *grLambdaC2Sys_PbPb=new TGraphAsymmErrors(20,meanNchPbPb,yAxisPbPb,xAxisC2_e,xAxisC2_e, yAxisPbPb_eL,yAxisPbPb_eH); + TGraphAsymmErrors *grLambdaC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxisC2_e,xAxisC2_e, yAxispPb_eL,yAxispPb_eH); + TGraphAsymmErrors *grLambdaC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxisC2_e,xAxisC2_e, yAxispp_eL,yAxispp_eH); grLambdaC2Sys_pp->SetMarkerSize(0); grLambdaC2Sys_pp->SetFillStyle(3001); grLambdaC2Sys_pp->SetFillColor(0); grLambdaC2Sys_pPb->SetMarkerSize(0); grLambdaC2Sys_pPb->SetFillStyle(3001); grLambdaC2Sys_pPb->SetFillColor(0); grLambdaC2Sys_PbPb->SetMarkerSize(0); grLambdaC2Sys_PbPb->SetFillStyle(3001); grLambdaC2Sys_PbPb->SetFillColor(0); @@ -1131,10 +1284,6 @@ void Plot_plotsTPR(){ grLambdaC2Sys_PbPb->SetLineWidth(2.*grLambdaC2Sys_PbPb->GetLineWidth()); // - grLambdaC2Sys_pp->Draw("|| p"); - grLambdaC2Sys_pPb->Draw("|| p"); - - grLambdaC2Sys_PbPb->Draw("|| p"); grLambdaSys_pp->Draw("E2 p"); grLambdaSys_pPb->Draw("E2 p"); @@ -1148,7 +1297,9 @@ void Plot_plotsTPR(){ LambdaC2pPb->Draw("same"); LambdaC2pp->Draw("same"); - + grLambdaC2Sys_pp->Draw("|| p"); + grLambdaC2Sys_pPb->Draw("|| p"); + grLambdaC2Sys_PbPb->Draw("|| p"); // print radii and lambda cout.precision(3); @@ -1156,24 +1307,29 @@ void Plot_plotsTPR(){ for(int cb=0; cb<20; cb++){ int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]); if(RadiiPbPb->GetBinContent(binPbPb)==0) continue; - cout<<"Nch="<GetBinError(binPbPb)<<" +- "<GetErrorY(cb)<<" lambda_3 = "<GetBinContent(binPbPb)<<" +- "<GetBinError(binPbPb)<<" +- "<GetErrorY(cb)<GetBinError(binPbPb)<<" +- "<GetErrorY(cb)<<" lambda_2 = "<GetBinContent(binPbPb)<<" +- "<GetBinError(binPbPb)<<" +- "<GetErrorY(cb)<GetBinContent(binPbPb)<<" "<GetBinError(binPbPb)<<" "<GetErrorY(cb)<<" "<GetBinContent(binPbPb)<<" "<GetBinError(binPbPb)<<" "<GetErrorY(cb)<GetBinContent(binPbPb)<<" "<GetBinError(binPbPb)<<" "<GetErrorY(cb)<<" "<GetBinContent(binPbPb)<<" "<GetBinError(binPbPb)<<" "<GetErrorY(cb)<GetXaxis()->FindBin(meanNchpPb[cb]); if(RadiipPb->GetBinContent(binpPb)==0) continue; - cout<<"Nch="<GetBinError(binpPb)<<" +- "<GetErrorY(cb)<<" lambda_3 = "<GetBinContent(binpPb)<<" +- "<GetBinError(binpPb)<<" +- "<GetErrorY(cb)<GetBinError(binpPb)<<" +- "<GetErrorY(cb)<<" lambda_2 = "<GetBinContent(binpPb)<<" +- "<GetBinError(binpPb)<<" +- "<GetErrorY(cb)<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorYhigh(cb)<<" "<GetErrorYlow(cb)<<" "<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorY(cb)<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorY(cb)<<" "<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorY(cb)<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorYhigh(cb)<<" "<GetErrorYlow(cb)<<" "<GetBinContent(binpPb)<<" "<GetBinError(binpPb)<<" "<GetErrorY(cb)<GetXaxis()->FindBin(meanNchpp[cb]); if(Radiipp->GetBinContent(binpp)==0) continue; - cout<<"Nch="<GetBinError(binpp)<<" +- "<GetErrorY(cb)<<" lambda_3 = "<GetBinContent(binpp)<<" +- "<GetBinError(binpp)<<" +- "<GetErrorY(cb)<GetBinError(binpp)<<" +- "<GetErrorY(cb)<<" lambda_2 = "<GetBinContent(binpp)<<" +- "<GetBinError(binpp)<<" +- "<GetErrorY(cb)<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorYhigh(cb)<<" "<GetErrorYlow(cb)<<" "<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorY(cb)<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorY(cb)<<" "<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorY(cb)<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorYhigh(cb)<<" "<GetErrorYlow(cb)<<" "<GetBinContent(binpp)<<" "<GetBinError(binpp)<<" "<GetErrorY(cb)<cd(); TLatex *RinvTitle; if(FitType==0) RinvTitle=new TLatex(0.062,0.72,"#font[12]{R}^{#font[12]{G}}_{inv} or #font[12]{R}^{#font[12]{G}}_{inv,3} (fm)"); - else RinvTitle=new TLatex(0.062,0.72,"#font[12]{R}^{#font[12]{E}_{w}}_{inv} or #font[12]{R}^{#font[12]{E}_{w}}_{inv,3} (fm)"); + else if(FitType==1) RinvTitle=new TLatex(0.062,0.72,"#font[12]{R}^{#font[12]{E}_{w}}_{inv} or #font[12]{R}^{#font[12]{E}_{w}}_{inv,3} (fm)"); + else RinvTitle=new TLatex(0.062,0.61,"(#font[12]{R}^{Exp}_{inv} or #font[12]{R}^{Exp}_{inv,3}) / #sqrt{#pi} (fm)"); RinvTitle->SetNDC(); RinvTitle->SetTextFont(TextFont); RinvTitle->SetTextSize(SizeTitle); @@ -1193,7 +1350,8 @@ void Plot_plotsTPR(){ RinvTitle->Draw("same"); TLatex *LambdaTitle; if(FitType==0) LambdaTitle=new TLatex(0.064,0.31,"#lambda^{#font[12]{G}}_{e} or #lambda^{#font[12]{G}}_{e,3}");// 0.064,0.33 - else LambdaTitle=new TLatex(0.064,0.31,"#lambda^{#font[12]{E}_{w}}_{e} or #lambda^{#font[12]{E}_{w}}_{e,3}");// 0.064,0.33 + else if(FitType==1) LambdaTitle=new TLatex(0.064,0.31,"#lambda^{#font[12]{E}_{w}}_{e} or #lambda^{#font[12]{E}_{w}}_{e,3}");// 0.064,0.33 + else LambdaTitle=new TLatex(0.064,0.31,"#lambda^{Exp}_{e} or #lambda^{Exp}_{e,3}");// 0.064,0.33 LambdaTitle->SetNDC(); LambdaTitle->SetTextFont(TextFont); LambdaTitle->SetTextSize(SizeTitle); @@ -1209,7 +1367,7 @@ void Plot_plotsTPR(){ /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // kappa plots - + /* TCanvas *can7 = new TCanvas("can7", "can7",1700,700,600,600);// 11,53,700,500 can7->SetHighLightColor(2); gStyle->SetOptFit(0);// 0111 to show fit stat box @@ -1248,21 +1406,21 @@ void Plot_plotsTPR(){ legend8->AddEntry(Parameters_c3[1][1][KT3Bin][paramNum],"p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV","p"); legend8->AddEntry(Parameters_c3[0][1][KT3Bin][paramNum],"Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV","p"); //legend8->Draw("same"); - /*for(int ct=0; ct<3; ct++){ - for(int cb=0; cb<20; cb++){ - int bin = 0; - if(ct==0) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->GetXaxis()->FindBin(meanNchPbPb[cb]); - else if(ct==1) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpPb[cb]); - else bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpp[cb]); + //for(int ct=0; ct<3; ct++){ + //for(int cb=0; cb<20; cb++){ + //int bin = 0; + //if(ct==0) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->GetXaxis()->FindBin(meanNchPbPb[cb]); + //else if(ct==1) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpPb[cb]); + //else bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpp[cb]); // - if(Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin) >0) { + //if(Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin) >0) { //cout<GetBinContent(bin)<<", "; - cout<GetBinError(bin)<<", "; - }else cout<<0<<", "; + //cout<GetBinError(bin)<<", "; + //}else cout<<0<<", "; - } - cout<Fit(Fit_kappa4_PbPb,"IMEN","",2,2000); Fit_kappa4_PbPb->Draw("same"); - + */ //////////////////////////////////////////////////// @@ -1358,7 +1516,6 @@ void Plot_plotsTPR(){ legend8->AddEntry(Ratio_pPb_to_pp,"p-Pb over pp","p"); legend8->AddEntry(Ratio_PbPb_to_pPb,"Pb-Pb over p-Pb","p"); legend8->Draw("same"); - Unity->Draw("same"); double avgStat_pPb_to_pp = sqrt(avgRatioStat_pPb_to_pp/avgRatioEn_pPb_to_pp); @@ -1404,6 +1561,7 @@ void Plot_plotsTPR(){ HIJING_c3_MC->SetBinContent(i, HIJING_c3_MC_K1_M3[i-1]); HIJING_c3_MC->SetBinError(i, HIJING_c3_MC_K1_M3_e[i-1]); } + cout.precision(4); // for(int padNum=1; padNum<=6; padNum++){ @@ -1426,8 +1584,21 @@ void Plot_plotsTPR(){ if(padNum==5) {System_proof=1; ChComb_proof=1; Mb_proof=14;} if(padNum==6) {System_proof=0; ChComb_proof=1; Mb_proof=3;} + // print out data points + //for(int binN=1; binN<=C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetNbinsX(); binN++){ + //cout<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<SetMinimum(0.9); - C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->SetMaximum(3.4);// 3.3 + C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->SetMaximum(3.4);// 3.4 // C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetTitle("#font[12]{C}_{3} or #font[12]{#bf{c}}_{3} "); if(padNum<=5){ @@ -1471,12 +1642,14 @@ void Plot_plotsTPR(){ } if(ChComb_proof==0) { - c3_fit[System_proof][0][KT3Bin][Mb_proof]->Draw("c same"); + c3_fit[System_proof][0][KT3Bin][Mb_proof]->Draw("c same");// Gauss gr_c3Spline[System_proof][KT3Bin][Mb_proof]->Draw("c same");// EW with spline for mid-q and high q + gr_c3SplineExpFit[System_proof][KT3Bin][Mb_proof]->Draw("c same");// Exp //c3_fit[System_proof][1][KT3Bin][Mb_proof]->Draw("c same");// old approximation if(padNum==3){ legendFitTypes->AddEntry(c3_fit[System_proof][0][KT3Bin][Mb_proof],"Gaussian","l"); legendFitTypes->AddEntry(c3_fit[System_proof][1][KT3Bin][Mb_proof],"Edgeworth","l"); + legendFitTypes->AddEntry(c3_fit[System_proof][2][KT3Bin][Mb_proof],"Exponential","l"); } } @@ -1574,7 +1747,7 @@ void Plot_plotsTPR(){ legend6->SetFillColor(0); legend6->SetTextFont(TextFont); legend6->SetTextSize(SizeLegend); - TLegend *legend7 = new TLegend(.67,.35, .98,.52,NULL,"brNDC");//.45 or .4 for x1 + TLegend *legend7 = new TLegend(.66,.36, .97,.53,NULL,"brNDC");//.45 or .4 for x1 legend7->SetBorderSize(0); legend7->SetFillColor(0); legend7->SetTextFont(TextFont); @@ -1586,10 +1759,10 @@ void Plot_plotsTPR(){ // int KT3Bin_CorrComp=0; int Mbin_SysComp_PbPb=12;// 12 - int Mbin_SysComp_pPb;// 12, 16, - int Mbin_SysComp_pp=15;// , 15, - if(p_pPb_Comp) Mbin_SysComp_pPb=16; - else Mbin_SysComp_pPb=12; + int Mbin_SysComp_pPb; + int Mbin_SysComp_pp=15;// 15 + if(pp_pPb_Comp) Mbin_SysComp_pPb=16;// 16 + else Mbin_SysComp_pPb=12;// 12 TH1D *c3_PbPb=(TH1D*)c3[0][0][0][KT3Bin_CorrComp][Mbin_SysComp_PbPb]->Clone(); TH1D *c3_pPb=(TH1D*)c3[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->Clone(); @@ -1604,50 +1777,59 @@ void Plot_plotsTPR(){ c3_pPb->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}^{#pm#pm#pm}"); c3_pPb->GetXaxis()->SetTitleOffset(1.0); c3_pPb->GetYaxis()->SetTitleOffset(0.7); - c3_pPb->SetMinimum(0.9); c3_pPb->SetMaximum(3.3); + c3_pPb->SetMinimum(0.9); + c3_pPb->SetMaximum(3.3);// 3.3 c3_pPb->GetXaxis()->SetNdivisions(503); c3_pPb->GetYaxis()->SetNdivisions(503); - + c3_pPb->SetMarkerStyle(25); + c3_pPb->Draw(); // - if(p_pPb_Comp) c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->Draw("E2 same"); + if(pp_pPb_Comp) c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->Draw("E2 same"); c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->Draw("E2 same"); - if(!p_pPb_Comp) c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("E2 same"); - if(p_pPb_Comp) c3_pp->Draw("same"); + if(!pp_pPb_Comp) c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("E2 same"); + if(pp_pPb_Comp) c3_pp->Draw("same"); c3_pPb->Draw("same"); - if(!p_pPb_Comp) c3_PbPb->Draw("same"); + if(!pp_pPb_Comp) c3_PbPb->Draw("same"); - if(p_pPb_Comp) { - legend6->AddEntry(c3_pPb,"#splitline{p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV}{#LT#font[12]{N}_{ch}#GT = 23 #pm 1}","p");// MpPb=16 - legend6->AddEntry(c3_pp,"#splitline{pp #sqrt{s}=7 TeV}{#LT#font[12]{N}_{ch}#GT = 27 #pm 1}","p");// Mpp=15 + if(pp_pPb_Comp) { + legend6->AddEntry(c3_pPb,"#splitline{p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV}{#LT#font[12]{N}_{ch}#GT = 23 #pm 1}","p");// MpPb=16, Nch=23 + legend6->AddEntry(c3_pp,"#splitline{pp #sqrt{s}=7 TeV}{#LT#font[12]{N}_{ch}#GT = 27 #pm 1}","p");// Mpp=15, Nch=27 }else{ - legend6->AddEntry(c3_pPb,"#splitline{p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV}{#LT#font[12]{N}_{ch}#GT = 71 #pm 4}","p");// MpPb=12 - legend6->AddEntry(c3_PbPb,"#splitline{Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV}{#LT#font[12]{N}_{ch}#GT = 84 #pm 4}","p");// MPbPb=12 + legend6->AddEntry(c3_pPb,"#splitline{p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV}{#LT#font[12]{N}_{ch}#GT = 71 #pm 4}","p");// MpPb=12, Nch=71 + legend6->AddEntry(c3_PbPb,"#splitline{Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV}{#LT#font[12]{N}_{ch}#GT = 84 #pm 4}","p");// MPbPb=12, Nch=84 } // - if(!p_pPb_Comp) c3_fit[0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); + if(!pp_pPb_Comp) c3_fit[0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); c3_fit[1][0][KT3Bin][Mbin_SysComp_pPb]->Draw("c same"); - if(p_pPb_Comp) c3_fit[2][0][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); + if(pp_pPb_Comp) c3_fit[2][0][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); // TF1 *GaussFit_forLegend=(TF1*)c3_fit[1][0][KT3Bin][Mbin_SysComp_pPb]->Clone(); TF1 *EwFit_forLegend=(TF1*)c3_fit[1][1][KT3Bin][Mbin_SysComp_pPb]->Clone(); + TF1 *ExpFit_forLegend=(TF1*)c3_fit[1][2][KT3Bin][Mbin_SysComp_pPb]->Clone(); GaussFit_forLegend->SetLineColor(1); EwFit_forLegend->SetLineColor(1); + ExpFit_forLegend->SetLineColor(1); legend7->AddEntry(GaussFit_forLegend,"Gaussian","l"); legend7->AddEntry(EwFit_forLegend,"Edgeworth","l"); + legend7->AddEntry(ExpFit_forLegend,"Exponential","l"); // - // old way with EW plotting approximation - //if(!p_pPb_Comp) c3_fit[0][1][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); - //c3_fit[1][1][KT3Bin][Mbin_SysComp_pPb]->Draw("c same"); - //if(p_pPb_Comp) c3_fit[2][1][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); - // new way with spline - if(!p_pPb_Comp) gr_c3Spline[0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); + + // spline draw + if(!pp_pPb_Comp) gr_c3Spline[0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); gr_c3Spline[1][KT3Bin][Mbin_SysComp_pPb]->Draw("c same"); - if(p_pPb_Comp) gr_c3Spline[2][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); + if(pp_pPb_Comp) gr_c3Spline[2][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); + // + + //exp draw + if(!pp_pPb_Comp) gr_c3SplineExpFit[0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same"); + gr_c3SplineExpFit[1][KT3Bin][Mbin_SysComp_pPb]->Draw("c same"); + if(pp_pPb_Comp) gr_c3SplineExpFit[2][KT3Bin][Mbin_SysComp_pp]->Draw("c same"); // + legend6->Draw("same"); legend7->Draw("same"); @@ -1658,6 +1840,11 @@ void Plot_plotsTPR(){ Specif_Kt3_4->SetTextSize(SizeSpecif); Specif_Kt3_4->Draw("same"); + // print out data points + for(int binN=1; binN<=c3_pPb->GetNbinsX(); binN++){ + if(pp_pPb_Comp) cout<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<GetXaxis()->GetBinLowEdge(binN)<<" "<GetXaxis()->GetBinUpEdge(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<<" "<GetBinContent(binN)<<" "<GetBinError(binN)<<" "<GetBinError(binN)<SetHighLightColor(2); + gStyle->SetOptFit(0);// 0111 to show fit stat box + can6->SetFillColor(0);//10 + can6->SetBorderMode(0); + can6->SetBorderSize(2); + can6->SetFrameFillColor(0); + can6->SetFrameBorderMode(0); + can6->SetFrameBorderMode(0); + can6->cd(); + TPad *pad6 = new TPad("pad6","pad6",0.0,0.0,1.,1.); + gPad->SetGridx(0); + gPad->SetGridy(0); + pad6->SetTopMargin(0.0);//0.05 + pad6->SetRightMargin(0.0);//3e-2 + pad6->SetBottomMargin(0.0);//0.12 + pad6->SetLeftMargin(0.0);//0.12 + pad6->Draw(); + pad6->cd(); + TLegend *legend8 = new TLegend(.17,.4, .5,.6,NULL,"brNDC");//.45 or .4 for x1 + legend8->SetBorderSize(0); + legend8->SetFillColor(0); + legend8->SetTextFont(TextFont); + legend8->SetTextSize(0.8*SizeLegend); + TLegend *legend9 = new TLegend(.17,.6, .6,.98,NULL,"brNDC");//.45 or .4 for x1 + legend9->SetBorderSize(0); + legend9->SetFillColor(0); + legend9->SetTextFont(TextFont); + legend9->SetTextSize(0.8*SizeLegend); + // + gPad->SetLeftMargin(0.14); + gPad->SetRightMargin(0.01); + gPad->SetTopMargin(0.01); + gPad->SetBottomMargin(0.14); - + gPad->SetTickx(); gPad->SetTicky(); + if(!NchOneThirdAxis) gPad->SetLogx(); + RadiiC2PbPb->GetXaxis()->SetLabelFont(TextFont); RadiiC2PbPb->GetYaxis()->SetLabelFont(TextFont); + RadiiC2PbPb->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); RadiiC2PbPb->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel); + RadiiC2PbPb->GetXaxis()->SetLabelOffset(-0.01); + RadiiC2PbPb->GetXaxis()->SetNdivisions(808); + RadiiC2PbPb->GetXaxis()->SetTitleOffset(1.05);//1.3 + RadiiC2PbPb->GetYaxis()->SetTitleOffset(1.1);//1.1 + RadiiC2PbPb->GetXaxis()->SetTitleFont(TextFont); RadiiC2PbPb->GetXaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel + RadiiC2PbPb->GetYaxis()->SetTitleFont(TextFont); RadiiC2PbPb->GetYaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel + RadiiC2PbPb->SetMinimum(0.01); RadiiC2PbPb->SetMaximum(11.9);// 0 and 11.9 + //gStyle->SetErrorX(0); + if(NchOneThirdAxis) RadiiC2PbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000 + else RadiiC2PbPb->GetXaxis()->SetRangeUser(3,3000);// 3,3000 + // + // + Parameters_Bjoern[0][0]->GetXaxis()->SetLabelFont(TextFont); Parameters_Bjoern[0][0]->GetYaxis()->SetLabelFont(TextFont); + Parameters_Bjoern[0][0]->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); Parameters_Bjoern[0][0]->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel); + Parameters_Bjoern[0][0]->GetXaxis()->SetLabelOffset(-0.01); + Parameters_Bjoern[0][0]->GetXaxis()->SetNdivisions(808); + Parameters_Bjoern[0][0]->GetXaxis()->SetTitleOffset(1.05);//1.3 + Parameters_Bjoern[0][0]->GetYaxis()->SetTitleOffset(1.1);//1.1 + Parameters_Bjoern[0][0]->GetXaxis()->SetTitleFont(TextFont); Parameters_Bjoern[0][0]->GetXaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel + Parameters_Bjoern[0][0]->GetYaxis()->SetTitleFont(TextFont); Parameters_Bjoern[0][0]->GetYaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel + Parameters_Bjoern[0][0]->SetMinimum(0.01); Parameters_Bjoern[0][0]->SetMaximum(8.5);// 0 and 11.9 + //gStyle->SetErrorX(0); + if(NchOneThirdAxis) Parameters_Bjoern[0][0]->GetXaxis()->SetRangeUser(0,3000);// 0,3000 + else Parameters_Bjoern[0][0]->GetXaxis()->SetRangeUser(3,3000);// 3,3000 + + + RadiiC2PbPb->GetXaxis()->SetTitle("#LT#font[12]{N}_{ch}#GT"); + RadiiC2PbPb->GetYaxis()->SetTitle("Radius (fm)"); + Parameters_Bjoern[0][0]->GetXaxis()->SetTitle("#LT#font[12]{N}_{ch}#GT"); + Parameters_Bjoern[0][0]->GetYaxis()->SetTitle("Radius (fm)"); + + RadiiC2PbPb->Draw(); + + + //int HydroCase=0;// 0 or 1 + //Parameters_Bjoern[HydroCase][0]->Draw("same"); + //Parameters_Bjoern[HydroCase][1]->Draw("same"); + //Parameters_Bjoern[HydroCase][2]->Draw("same"); + + legend8->AddEntry(RadiiC2pp,"ALICE pp","p"); + legend8->AddEntry(RadiiC2pPb,"ALICE p-Pb","p"); + legend8->AddEntry(RadiiC2PbPb,"ALICE Pb-Pb","p"); + // + //legend8->AddEntry(Parameters_Bjoern[HydroCase][2],"IP-GLASMA pp (w/o hydro)","p"); + //legend8->AddEntry(Parameters_Bjoern[HydroCase][1],"IP-GLASMA p-Pb (w/o hydro)","p"); + //legend8->AddEntry(Parameters_Bjoern[HydroCase][0],"IP-GLASMA Pb-Pb (w/o hydro)","p"); + + + Parameters_Bjoern[0][0]->SetMarkerStyle(24); Parameters_Bjoern[0][1]->SetMarkerStyle(25); Parameters_Bjoern[0][2]->SetMarkerStyle(28); + + //TMultiGraph *mg = new TMultiGraph(); + //mg->SetTitle(""); + TGraph *Radii_Bjoern[2][3];// Hydro case, CollType + TGraph *grShade[3];// CollType + + Radii_Bjoern[0][0] = new TGraph(10, Bjoern_xaxis_PbPb, Bjoern_Ri_PbPb[mchoice-1]); + Radii_Bjoern[0][0]->SetLineWidth(2); + + Radii_Bjoern[1][0] = new TGraph(10, Bjoern_xaxis_PbPb, Bjoern_Rhydro_PbPb[mchoice-1]); + Radii_Bjoern[1][0]->SetLineWidth(2); + Radii_Bjoern[1][0]->SetLineStyle(2); + Radii_Bjoern[1][0]->SetFillStyle(1000); + Radii_Bjoern[1][0]->SetFillColor(kGray); + grShade[0] = new TGraph(10); + for(int i=0; i<10; i++){ + grShade[0]->SetPoint(i,Bjoern_xaxis_PbPb[i],Bjoern_Rhydro_PbPb[mchoice-1][i]); + grShade[0]->SetPoint(10+i,Bjoern_xaxis_PbPb[10-i-1],Bjoern_Ri_PbPb[mchoice-1][10-i-1]); + } + grShade[0]->SetFillStyle(1000); + grShade[0]->SetFillColor(kGray); + // + Radii_Bjoern[0][1] = new TGraph(6, Bjoern_xaxis_pPb, Bjoern_Ri_pPb[mchoice-1]); + Radii_Bjoern[0][1]->SetLineWidth(2); + Radii_Bjoern[0][1]->SetLineColor(2); + Radii_Bjoern[1][1] = new TGraph(6, Bjoern_xaxis_pPb, Bjoern_Rhydro_pPb[mchoice-1]); + Radii_Bjoern[1][1]->SetLineWidth(2); + Radii_Bjoern[1][1]->SetLineColor(2); + Radii_Bjoern[1][1]->SetLineStyle(2); + grShade[1] = new TGraph(6); + for(int i=0; i<6; i++){ + grShade[1]->SetPoint(i,Bjoern_xaxis_pPb[i],Bjoern_Rhydro_pPb[mchoice-1][i]); + grShade[1]->SetPoint(6+i,Bjoern_xaxis_pPb[6-i-1],Bjoern_Ri_pPb[mchoice-1][6-i-1]); + } + grShade[1]->SetFillStyle(1000); + grShade[1]->SetFillColor(kRed-10); + // + Radii_Bjoern[0][2] = new TGraph(4, Bjoern_xaxis_pp, Bjoern_Ri_pp[mchoice-1]); + Radii_Bjoern[0][2]->SetLineWidth(2); + Radii_Bjoern[0][2]->SetLineColor(4); + Radii_Bjoern[1][2] = new TGraph(4, Bjoern_xaxis_pp, Bjoern_Rhydro_pp[mchoice-1]); + Radii_Bjoern[1][2]->SetLineWidth(2); + Radii_Bjoern[1][2]->SetLineColor(4); + Radii_Bjoern[1][2]->SetLineStyle(2); + grShade[2] = new TGraph(4); + for(int i=0; i<4; i++){ + grShade[2]->SetPoint(i,Bjoern_xaxis_pp[i],Bjoern_Rhydro_pp[mchoice-1][i]); + grShade[2]->SetPoint(4+i,Bjoern_xaxis_pp[4-i-1],Bjoern_Ri_pp[mchoice-1][4-i-1]); + } + grShade[2]->SetFillStyle(1000); + grShade[2]->SetFillColor(kBlue-10); + // + /*grShade[0]->Draw("f same"); + grShade[1]->Draw("f same"); + grShade[2]->Draw("f same"); + Radii_Bjoern[0][0]->Draw("l same"); + Radii_Bjoern[0][1]->Draw("l same"); + Radii_Bjoern[0][2]->Draw("l same"); + */ + grRadiiC2Sys_pp->Draw("|| p"); + grRadiiC2Sys_pPb->Draw("|| p"); + grRadiiC2Sys_PbPb->Draw("|| p"); + RadiiC2PbPb->Draw("same"); + RadiiC2pPb->Draw("same"); + RadiiC2pp->Draw("same"); + + //legend9->AddEntry(Radii_Bjoern[0][2],"GLASMA pp R_{initial}","l"); + //legend9->AddEntry(Radii_Bjoern[0][1],"GLASMA p-Pb R_{initial}","l"); + //legend9->AddEntry(Radii_Bjoern[0][0],"GLASMA Pb-Pb R_{initial}","l"); + //legend9->AddEntry(grShade[2],"GLASMA pp R_{hydro}","f"); + //legend9->AddEntry(grShade[1],"GLASMA p-Pb R_{hydro}","f"); + //legend9->AddEntry(grShade[0],"GLASMA Pb-Pb R_{hydro}","f"); + /* + TF1 *fit1=new TF1("fit1","pol5",0,1000); + fit1->SetLineColor(1); + Radii_Bjoern[0][0]->Fit(fit1,"Q","",30,800); + //fit1->Draw("same"); + // + TF1 *fit2=new TF1("fit2","pol4",0,1000); + Radii_Bjoern[0][1]->Fit(fit2,"Q","",4,70); + //fit2->Draw("same"); + // + TF1 *fit3=new TF1("fit3","pol4",0,1000); + fit3->SetLineColor(1); + Radii_Bjoern[0][2]->Fit(fit3,"Q","",4,27); + //fit3->Draw("same"); + + TH1D *BjoernRatio_PbPb = (TH1D*)RadiiC2PbPb->Clone(); + TH1D *BjoernRatio_pPb = (TH1D*)RadiiC2pPb->Clone(); + TH1D *BjoernRatio_pp = (TH1D*)RadiiC2pp->Clone(); + for(int point=0; point<20; point++){ + int bin = BjoernRatio_PbPb->GetXaxis()->FindBin(meanNchPbPb[point]); + if(meanNchPbPb[point] > Bjoern_xaxis_PbPb[8]) {BjoernRatio_PbPb->SetBinContent(bin,0); continue;} + if(BjoernRatio_PbPb->GetBinContent(bin) < 0.5 || BjoernRatio_PbPb->GetBinContent(bin) > 12) BjoernRatio_PbPb->SetBinContent(bin,0); + else{ + // double xx = + BjoernRatio_PbPb->SetBinContent(bin, BjoernRatio_PbPb->GetBinContent(bin) / fit1->Eval(meanNchPbPb[point])); + } + } + for(int point=0; point<20; point++){ + int bin = BjoernRatio_pPb->GetXaxis()->FindBin(meanNchpPb[point]); + if(meanNchpPb[point] > Bjoern_xaxis_pPb[5]) {BjoernRatio_pPb->SetBinContent(bin,0); continue;} + if(BjoernRatio_pPb->GetBinContent(bin) < 0.5 || BjoernRatio_pPb->GetBinContent(bin) > 12) BjoernRatio_pPb->SetBinContent(bin,0); + else{ + BjoernRatio_pPb->SetBinContent(bin, BjoernRatio_pPb->GetBinContent(bin) / fit2->Eval(meanNchpPb[point])); + } + } + for(int point=0; point<20; point++){ + int bin = BjoernRatio_pp->GetXaxis()->FindBin(meanNchpp[point]); + if(meanNchpp[point] > Bjoern_xaxis_pp[3]) {BjoernRatio_pp->SetBinContent(bin,0); continue;} + if(BjoernRatio_pp->GetBinContent(bin) < 0.5 || BjoernRatio_pp->GetBinContent(bin) > 12) BjoernRatio_pp->SetBinContent(bin,0); + else{ + BjoernRatio_pp->SetBinContent(bin, BjoernRatio_pp->GetBinContent(bin) / fit2->Eval(meanNchpp[point])); + } + } + BjoernRatio_PbPb->SetMinimum(0.55); BjoernRatio_PbPb->SetMaximum(1.45); + BjoernRatio_PbPb->GetYaxis()->SetTitle("Radius Ratio (ALICE/GLASMA)"); + BjoernRatio_PbPb->Draw(); + BjoernRatio_pPb->Draw("same"); + BjoernRatio_pp->Draw("same"); + */ + + //legend8->AddEntry(Parameters_Bjoern[0][2],"pp R_{initial} (no hydro)","p"); + //legend8->AddEntry(Parameters_Bjoern[0][1],"p-Pb R_{initial} (no hydro)","p"); + //legend8->AddEntry(Parameters_Bjoern[0][0],"Pb-Pb R_{initial} (no hydro)","p"); + //legend9->AddEntry(Parameters_Bjoern[1][2],"pp R_{max} (hydro)","p"); + //legend9->AddEntry(Parameters_Bjoern[1][1],"p-Pb R_{max} (hydro)","p"); + //legend9->AddEntry(Parameters_Bjoern[1][0],"Pb-Pb R_{max} (hydro)","p"); + + legend8->Draw("same"); + legend9->Draw("same"); // Normalization plots // ct, ft, KT3, par //Parameters_C2[2][0][0][0]->Draw(); -- 2.43.0