--- /dev/null
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+#include <assert.h>
+
+#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<<f33 + 3*f32 + f31<<endl;
+
+
+ cout<<"Mbin = "<<Mbin<<" KT3 = "<<EDbin<<endl;
+
+ if(CollisionType==0){
+ if(Mbin==0) {RbinMRC=10;}
+ else if(Mbin==1) {RbinMRC=9;}
+ else if(Mbin<=3) {RbinMRC=8;}
+ else if(Mbin<=5) {RbinMRC=7;}
+ else {RbinMRC=6;}
+ }else{
+ RbinMRC=2;
+ }
+
+
+ if(CollisionType==0) BINLIMIT_3=20;
+ else BINLIMIT_3=30;
+
+ // bin centers from QS+FSI
+ double BinCenterPbPbCentral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509};
+ double BinCenterpPbAndpp[40]={0.00359275, 0.016357, 0.0257109, 0.035445, 0.045297, 0.0552251, 0.0651888, 0.0751397, 0.0851088, 0.0951108, 0.105084, 0.115079, 0.12507, 0.135059, 0.145053, 0.155049, 0.16505, 0.175038, 0.185039, 0.195034, 0.205023, 0.215027, 0.225024, 0.235023, 0.245011, 0.255017, 0.265017, 0.275021, 0.285021, 0.295017, 0.305018, 0.315018, 0.325013, 0.335011, 0.345016, 0.355019, 0.365012, 0.375016, 0.385017, 0.395016};
+ if(CollisionType==0){
+ for(int i=0; i<40; i++) BinCenters[i] = BinCenterPbPbCentral[i];
+ }else{
+ for(int i=0; i<40; i++) BinCenters[i] = BinCenterpPbAndpp[i];
+ }
+
+ // extend BinCenters for high q
+ for(int index=40; index<400; index++){
+ if(CollisionType==0) BinCenters[index] = (index+0.5)*(0.005);
+ else BinCenters[index] = (index+0.5)*(0.010);
+ }
+ // Set 0's to 3-particle fit arrays
+ for(int i=1; i<=BINLIMIT_3; i++){// bin number
+ for(int j=1; j<=BINLIMIT_3; j++){// bin number
+ for(int k=1; k<=BINLIMIT_3; k++){// bin number
+ A_3[i-1][j-1][k-1]=0;
+ A_3_e[i-1][j-1][k-1]=0;
+ B_3[i-1][j-1][k-1]=0;
+ }
+ }
+ }
+
+
+ //
+ TH3D *ThreeParticle[2][2][2][5];// ch1,ch2,ch3,term
+ TProfile3D *K3avg[2][2][2][4];
+ double norm_3[5]={0};
+ //
+
+ gStyle->SetOptStat(0);
+ gStyle->SetOptDate(0);
+ gStyle->SetOptFit(0111);
+
+
+
+ SetFSIindex(10.);
+ SetFSICorrelations();
+ SetMomResCorrections();
+ SetMuonCorrections();
+ //
+ /////////////////////////////////////////////////////////
+
+
+
+ TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
+ //
+
+ cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
+
+
+
+ ///////////////////////////////////
+ // Get Histograms
+
+ for(int term=0; term<5; term++){
+
+ TString *name3 = new TString("ThreeParticle_Charge1_");
+ *name3 += 0;
+ name3->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 "<<norm_3[term]<<endl;
+ ThreeParticle[0][0][0][term]->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"<<endl;
+
+ TF1 *Unity=new TF1("Unity","1",0,100);
+ Unity->SetLineStyle(2);
+
+
+ int ch1=0,ch2=0,ch3=0;
+
+
+
+ ///////////////////////////////////////////////////////////
+ // 3-pion
+ cout<<"3-pion section"<<endl;
+
+ TCanvas *can2 = new TCanvas("can2", "can2",600,53,700,500);
+ can2->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<<A_3[i-1][j-1][k-1]<<" "<<B_3[i-1][j-1][k-1]<<" "<<A_3_e[i-1][j-1][k-1]<<endl;
+ ///////////////////////////////////////////////////////////
+
+ }
+ }
+ }
+
+
+
+ ////////////////////////////
+
+ // Intermediate steps
+ for(int i=0; i<Q3BINS; i++) {c3_hist->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<npar_c3; i++){
+ MyMinuit_c3.DefineParameter(i, parName_c3[i].c_str(), par_c3[i], stepSize_c3[i], minVal_c3[i], maxVal_c3[i]);
+ }
+ if(FitType==0 || FitType==1) { MyMinuit_c3.FixParameter(6);}
+ if(FitType==2){
+ MyMinuit_c3.FixParameter(3);
+ MyMinuit_c3.FixParameter(4);
+ MyMinuit_c3.FixParameter(5);
+ }
+ MyMinuit_c3.FixParameter(0);
+ //MyMinuit_c3.FixParameter(1);
+ //MyMinuit_c3.FixParameter(2);
+ //MyMinuit_c3.FixParameter(3);
+ //MyMinuit_c3.FixParameter(4);
+ MyMinuit_c3.FixParameter(5);
+
+ /////////////////////////////////////////////////////////////
+ // Do the minimization!
+ cout<<"Start Three-d fit"<<endl;
+ MyMinuit_c3.Migrad();// Minuit's best minimization algorithm
+ cout<<"End Three-d fit"<<endl;
+ /////////////////////////////////////////////////////////////
+ MyMinuit_c3.mnexcm("SHOw PARameters", &arglist_c3, 1, ierflg_c3);
+ cout<<"c3 Fit: Chi2 = "<<Chi2_c3<<" NDF = "<<(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
+ cout<<" Chi2/NDF = "<<Chi2_c3 / (NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
+
+ for(int i=0; i<npar_c3; i++){
+ MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]);
+ }
+
+ cout<<"Tij Norm = "<<pow(OutputPar_c3[1]/2., 1/3.)<<endl;
+
+ if(FitType!=2){
+ c3Fit1D_Expan->FixParameter(0,OutputPar_c3[0]);
+ c3Fit1D_Expan->FixParameter(1,OutputPar_c3[1]);
+ c3Fit1D_Expan->FixParameter(2,OutputPar_c3[2]);
+ c3Fit1D_Expan->FixParameter(3,OutputPar_c3[3]);
+ c3Fit1D_Expan->FixParameter(4,OutputPar_c3[4]);
+ c3Fit1D_Expan->FixParameter(5,OutputPar_c3[5]);
+ }else{// Levy
+ c3Fit1D_Expan->FixParameter(0,OutputPar_c3[0]);
+ c3Fit1D_Expan->FixParameter(1,OutputPar_c3[1]);
+ c3Fit1D_Expan->FixParameter(2,OutputPar_c3[2]);
+ c3Fit1D_Expan->FixParameter(3,OutputPar_c3[6]);
+ }
+ c3Fit1D_Expan->SetLineStyle(1);
+ //c3Fit1D_Expan->Draw("same");
+
+
+ // project 3D EW fit onto 1D histogram
+ for(int i=2; i<=BINLIMIT_3; i++){// bin number
+ double Q12 = BinCenters[i-1];// true center
+ for(int j=2; j<=BINLIMIT_3; j++){// bin number
+ double Q13 = BinCenters[j-1];// true center
+ for(int k=2; k<=BINLIMIT_3; k++){// bin number
+ double Q23 = BinCenters[k-1];// true center
+ //
+ double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
+
+ 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 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<<Q3<<" "<<C<<endl;
+ C *= TERM5;
+ C *= OutputPar_c3[0];
+ //if(Q3<0.018) continue;
+ GenSignalExpected_num->Fill(Q3, C);
+ GenSignalExpected_den->Fill(Q3, TERM5);
+ //if(Q3<0.02) cout<<Q3<<" "<<TERM5<<endl;
+ ///////////////////////////////////////////////////////////
+
+
+
+ }
+ }
+ }
+
+
+ GenSignalExpected_num->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<<c3_hist->GetBinContent(bin)<<" "<<ypoints[grIndex]<<" "<<c3_hist->GetBinError(bin)<<endl;
+ cout<<pow((c3_hist->GetBinContent(bin)-ypoints[grIndex]) / c3_hist->GetBinError(bin),2)<<endl;
+ SumNDF_1D++;
+ }
+ cout<<"1D Chi2/NDF = "<<ChiSqSum_1D / (SumNDF_1D-5.)<<endl;
+ */
+
+
+
+
+ if(SaveToFile){
+ TString *savefilename = new TString("FitFiles/FitFile_CT");
+ *savefilename += CollisionType;
+ savefilename->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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
+ }
+
+ TH1D *temphistoSS[12];
+ TH1D *temphistoOS[12];
+ for(Int_t MB=0; MB<12; MB++) {
+ TString *nameK2SS = new TString("K2ss_");
+ *nameK2SS += MB;
+ temphistoSS[MB] = (TH1D*)fsifile->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"<<endl;}
+ */
+
+ NFitPoints_c3++;
+
+ }
+ }
+ }
+ //f = -2.0*lnL;// log-liklihood minimization
+ f = SumChi2;// Chi2 minimization
+ Chi2_c3 = f;
+ //Chi2_c3 = SumChi2_test;
+
+}
--- /dev/null
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+
+#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();
+
+
+}
// 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)
//
//
//
-bool PbPbcase=kTRUE;// always kTRUE
int fFSIindex=0;
int Ktlowbin;
int Kthighbin;
double NormQcutLow_PbPb = .15;
double NormQcutHigh_PbPb = .2;// was .175
-int ThreeParticleRebin=2;
-int FourParticleRebin=3;
+int ThreeParticleRebin;
+int FourParticleRebin;
int RbinMRC;
ThreeFrac = pow(TwoFrac, 1.5);
FourFrac = pow(TwoFrac, 2.);
- cout<<"Mbin = "<<Mbin<<" KT3 = "<<EDbin<<" Kt = "<<Ktbin<<" SameCharge = "<<SameCharge<<endl;
+ //// 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<<f33 + 3*f32 + f31<<endl;
+
+ //// Core/Halo, 40fm, 70fm, 100fm
+ float ThermShift_f44[4]={0., 0.08741, 0.005284, -0.01064};
+ float ThermShift_f43[4]={0., -0.01126, -0.001486, 0.001686};
+ float ThermShift_f42[4]={0., -0.006466, -7.683e-05, 0.0004572};
+ float ThermShift_f41[4]={0., -0.003556, 0.00112, 0.00115};
+ float f44prime = FourFrac;
+ float f43prime = ThreeFrac*(1-OneFrac);
+ float f42prime = TwoFrac*pow(1-OneFrac,2);
+ float f41prime = pow(1-OneFrac,4) + 4*OneFrac*pow(1-OneFrac,3);
+ f44prime += ThermShift_f44[f_choice];
+ f43prime += ThermShift_f43[f_choice];
+ f42prime += ThermShift_f42[f_choice];
+ f41prime += ThermShift_f41[f_choice];
+ float f44 = f44prime;
+ float f43 = f43prime/f33prime;
+ float f42 = f42prime/TwoFrac;
+ f42 -= 2*f43prime*f32prime/f33prime/TwoFrac;
+ float f41 = f41prime;
+ f41 -= 4*f43prime*f31prime/f33prime;
+ f41 -= 6*f42prime*(1-TwoFrac)/TwoFrac;
+ f41 += 12*f43prime/f33prime*f32prime/TwoFrac*(1-TwoFrac);
+ //cout<<f44 + 4*f43 + 6*f42 + f41<<endl;
+ //cout<<f44<<" "<<f43<<" "<<f42<<" "<<f41<<endl;
+ //
+ // Core/Halo, 40fm, 70fm, 100fm
+ //float TherminatorMod_f33[4]={1., 1.123, 1.022, 1.008};// 1.008
+ //float TherminatorMod_f32[4]={1., 0.844, 0.933, 0.967};// .967
+ //float TherminatorMod_f31[4]={1., 0.828, 0.982, 1.028};// 1.028
+ /*float TherminatorMod_f44[4]={1., 1.188, 1.008, 0.985};// .985
+ float TherminatorMod_f43[4]={1., 0.885, 0.979, 1.026};// 1.026
+ float TherminatorMod_f42[4]={1., 0.687, 0.99, 1.08};// 1.08
+ float TherminatorMod_f41[4]={1., 0.806, 1.33, 1.548};//1.548
+ 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);
+ f44 *= TherminatorMod_f44[f_choice];
+ f43 *= TherminatorMod_f43[f_choice];
+ f42 *= TherminatorMod_f42[f_choice];
+ f41 *= TherminatorMod_f41[f_choice];*/
+ //f44prime *= TherminatorMod_f44[f_choice];
+ //f43prime *= TherminatorMod_f43[f_choice];
+ //f42prime *= TherminatorMod_f42[f_choice];
+ //f41prime *= TherminatorMod_f41[f_choice];
+ //float f44 = f44prime;
+ //float f43 = f43prime/ThreeFrac;
+ //float f42 = f42prime/TwoFrac - 2*pow(1-OneFrac,2);
+ //float f41 = f41prime - 4*pow(1-OneFrac,4) - 12*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2);
+
+ cout<<"Mbin = "<<Mbin<<" KT3 = "<<EDbin<<" SameCharge = "<<SameCharge<<endl;
if(!SameCharge) cout<<"4-pion MixedCharge type = "<<MixedCharge4pionType<<endl;
- if(Mbin==0) {RbinMRC=10;}
- else if(Mbin==1) {RbinMRC=9;}
- else if(Mbin<=3) {RbinMRC=8;}
- else if(Mbin<=5) {RbinMRC=7;}
- else {RbinMRC=6;}
-
-
+ if(CollisionType==0){
+ if(Mbin==0) {RbinMRC=10;}
+ else if(Mbin==1) {RbinMRC=9;}
+ else if(Mbin<=3) {RbinMRC=8;}
+ else if(Mbin<=5) {RbinMRC=7;}
+ else {RbinMRC=6;}
+ }else{
+ RbinMRC=2;
+ }
+
+ if(MCcase) MuonCorrection=kFALSE;
+ if(CollisionType==0){
+ ThreeParticleRebin=2;
+ FourParticleRebin=3;
+ }else{
+ ThreeParticleRebin=2;
+ FourParticleRebin=6;
+ }
float Cutoff_FullWeight_Q3[10]={0.065, 0.065, 0.075, 0.075, 0.075, 0.075, 0.085, 0.085, 0.085, 0.085};
float Cutoff_FullWeight_Q4[10]={0.115, 0.115, 0.115, 0.130, 0.130, 0.130, 0.145, 0.145, 0.145, 0.145};
TH1D *TPFullWeight_FourParticle[2];// charge
TH1D *TPNegFullWeight_FourParticle[2];// charge
TH1D *TPN_FourParticle[2];// charge
+ TH3D *FullBuildFromFits_3D;// charge
+ TH1D *FullBuildFromFits;// charge
+ TH3D *PartialBuildFromFits_3D;// charge
+ TH1D *PartialBuildFromFits;// charge
double norm_4[13]={0};
gStyle->SetOptFit(0111);
TFile *_file0;
- if(PbPbcase){// PbPb
+ if(CollisionType==0){// PbPb
if(MCcase) {
if(Mbin<=1){
- _file0 = new TFile("Results/PDC_12a17a.root","READ");
+ //_file0 = new TFile("Results/PDC_12a17a_GeneratedSignal.root","READ");
+ //_file0 = new TFile("Results/PDC_12a17a_pTSpectrumWeight.root","READ");
+ _file0 = new TFile("Results/PDC_12a17a_Qweights.root","READ");
}
}else{
- if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_extendedGweights.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0_FullRunWrongWeightsNoPadRowTTCandMCTTC_RawWeightFile.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p16to0p25_0p25to1p0.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT0p16to0p25_KT0p35.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_1percentCentEM.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_0p02eta0p045phi_0p03eta0p067phi.root","READ");
+ if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_c3FitBuild.root","READ");
+ //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_extendedGweights.root","READ");// Preliminary results
if(FileSetting==5) _file0 = new TFile("Results/PDC_11h_Cubic_Linear.root","READ");
if(FileSetting==1 || FileSetting==2) _file0 = new TFile("Results/PDC_11h_Lam0p65_Lam0p75.root","READ");
if(FileSetting==3) _file0 = new TFile("Results/PDC_11h_Cubic_Linear_Bminus.root","READ");
if(FileSetting==6) _file0 = new TFile("Results/PDC_10h_Cubic_Linear.root","READ");
if(FileSetting==7 || FileSetting==8) _file0 = new TFile("Results/PDC_11h_MRC10percIncrease_Muon92percent.root","READ");
}
+ }else if(CollisionType==1){// pPb
+ _file0 = new TFile("Results/PDC_13bc_kINT7.root","READ");
}else{// pp
- cout<<"pp not currently supported"<<endl;
- return;
+ _file0 = new TFile("Results/PDC_10bcde_kMB.root","READ");
}
SetFSIindex(10.);
double NormQcutLow;
double NormQcutHigh;
- if(PbPbcase) {
+ if(CollisionType==0) {
NormQcutLow = NormQcutLow_PbPb;
NormQcutHigh = NormQcutHigh_PbPb;
}else {
TwoParticle[c1][c2][term] = (TH1D*)TwoParticle_2d[c1][c2][term]->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 "<<norm_2[term]<<endl;
+ //cout<<"2-pion norms "<<norm_2[term]<<endl;
TwoParticle[c1][c2][term]->Scale(norm_2[0]/norm_2[term]);
TwoParticle[c1][c2][term]->SetMarkerStyle(20);
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){
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<<binx<<" "<<biny<<" "<<q3<<" "<<InterCorr<<endl;
+ TPFullWeight_ThreeParticle_2D[c1]->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);
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);
//
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");
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
FourParticle[c1][c2][c3][c4][term]->SetTitle("");
//
FourParticle[c1][c2][c3][c4][term]->Rebin(FourParticleRebin);
-
+
}else{
- cout<<"normalization = 0."<<endl;
+ cout<<"4-particle normalization = 0."<<endl;
}
if(term<12){
K4avg[c1][c2][c3][c4][term] = (TProfile*)MyList->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){
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");
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());
cout<<temperr->GetBinContent(3)<<endl;
cout<<(temperr->GetBinContent(5) / tempDen->GetBinContent(5))<<" "<<TPFullWeight_FourParticle[c1]->GetBinContent(5)<<endl;
*/
+ if(c1==0){
+ FullBuildFromFits = (TH1D*)FullBuildFromFits_3D->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
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"<<endl;
TF1 *Unity=new TF1("Unity","1",0,100);
TERM1_2->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<<C2QS->GetBinContent(bin)<<", ";
+ }
+ cout<<endl;
+ for(int bin=1; bin<=20; bin++){
+ cout<<C2QS->GetBinError(bin)<<", ";
+ }
+ cout<<endl;*/
+ double C2refPoints_0p65[20]={-0.218603, 1.00087, 1.00215, 1.00209, 1.00202, 1.00215, 1.00162, 1.00129, 1.00118, 1.00076, 1.00049, 1.00021, 0.999996, 0.999733, 0.999622, 0.999554, 0.999764, 0.999679, 0.999677, 0.999641};
+ double C2refPoints_e_0p65[20]={0.309152, 0.000473757, 0.000305599, 0.00022603, 0.000179596, 0.000149295, 0.000128009, 0.000112371, 0.00010047, 9.11363e-05, 8.37061e-05, 7.76871e-05, 7.27445e-05, 6.86217e-05, 6.5155e-05, 6.22148e-05, 5.9718e-05, 5.75613e-05, 5.57066e-05, 5.41028e-05};
+ double C2refPoints_0p75[20]={-0.135326, 0.968257, 0.984681, 0.991465, 0.99503, 0.997249, 0.998093, 0.998645, 0.999135, 0.999172, 0.999237, 0.999213, 0.999197, 0.999095, 0.999105, 0.999123, 0.999385, 0.999366, 0.999408, 0.999411};
+ double C2refPoints_e_0p75[20]={0.19138, 0.000422359, 0.000272193, 0.000201183, 0.000159788, 0.000132797, 0.000113845, 9.99271e-05, 8.93377e-05, 8.10334e-05, 7.44239e-05, 6.90702e-05, 6.46743e-05, 6.10077e-05, 5.79248e-05, 5.53103e-05, 5.30903e-05, 5.11725e-05, 4.95234e-05, 4.80974e-05};
+ TH1D *C2ref_0p65=(TH1D*)C2QS->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");
+
+
///////////////////////////////////////////////////////////
TERM3_3->Multiply(MRC_SC_3[1]);
TERM4_3->Multiply(MRC_SC_3[1]);
TERM5_3->Multiply(MRC_SC_3[2]);
+ //cout<<MRC_SC_3[2]->GetBinContent(2)<<" "<<MRC_SC_3[2]->GetBinContent(3)<<" "<<MRC_SC_3[2]->GetBinContent(4)<<endl;
}else{
if(CHARGE==-1){// --+ but MRC is stored as -++
TERM1_3->Multiply(MRC_MC_3[0]);
}
}
- //
- 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);
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));
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<<TPFullWeight_ThreeParticle[ch1_3]->GetBinContent(2)<<endl;
//TH1D *C3raw = (TH1D*)TERM1_3->Clone();
//C3raw->Divide(TERM5_3);
//C3raw->SetMarkerColor(4);
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);
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);
// 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;
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");
// 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]);
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");
// 4-pion
cout<<"4-pion section"<<endl;
- TCanvas *can3 = new TCanvas("can3", "can3",11,600,700,500);
+ TCanvas *can3 = new TCanvas("can3", "can3",11,600,700,500);// 60 was 600
can3->SetHighLightColor(2);
can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
gStyle->SetOptFit(0111);
}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();
}
//
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<<f44 + 4*f43 + 6*f42 + f41<<endl;
- //cout<<f44<<" "<<f43<<" "<<f42<<" "<<f41<<endl;
+
TH1D *N4QS=(TH1D*)TERMS_4[0]->Clone();
//
N4QS->Add(TERMS_4[1], -f43);
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++){
//
}
// 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;
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);
//
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);
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);
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");
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<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(9)<<endl;
+
+ if(SameCharge) {
+ //TPFullWeight_FourParticle[ch1_4]->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");
//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);
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);
// 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;
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;
// Print 4-pion points
for(int bin=1; bin<=12; bin++){
- cout<<C4QS->GetBinContent(bin)<<", ";
+ //cout<<C4QS->GetBinContent(bin)<<", ";
//cout<<c4QS->GetBinContent(bin)<<", ";
//cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(bin)<<", ";
//cout<<C4raw->GetBinContent(bin)<<", ";
}
cout<<endl;
for(int bin=1; bin<=12; bin++){
+ //cout<<c4QS->GetBinContent(bin)<<", ";
//cout<<C4QS->GetBinError(bin)<<", ";
//cout<<c4QS->GetBinError(bin)<<", ";
//cout<<C4raw->GetBinError(bin)<<", ";
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);
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");
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);
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];
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);
//________________________________________________________________________
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%
#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?
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<<Q2Limit<<endl;
+
// PbPb then pPb then pp
//double Nch_means[3][20]={{1828,1407,1071,896,757,619,496,402,324,252,171,117,83,63,50,36, 36,36,36,36}, {71,71,71,71,71,71,71,71,71,71,71,71, 71,57,45,32,23,16,9.7,5}, {47,47,47,47,47,47,47,47,47,47,47,47,47, 47,37,27,20,15,8.6,4.6}};
- kappa3 = 0.1;
- kappa4 = 0.5;
+ //kappa3 = 0.1;
+ //kappa4 = 0.5;
//kappa3 = 0.05 + 0.25*(Mbin/19.);// linear dependence of kappa3
//double MainFitParams[8]={0};
//
- //double RefMainFitParams[8]={3.58014, 2.84635, 0.652247, 0.831958, 3.13505, 2.33717, 0.761147, 0.909211};
+ //double RefMainFitParams[8]={6.71377, 5.39974, 1.3822, 1.82065, 6.71377, 5.39974, 1.3822, 1.82065};
if(Gaussian) ExpPower=2.;
else ExpPower=1.;
- // old way v5 and before
- /*if(Mbin==0) FSIindex = 0;// 0-5% Therm
- else if(Mbin==1) FSIindex = 1;// 0-10% Therm
- else if(Mbin<=3) FSIindex = 2;// 10-20% Therm
- else if(Mbin<=5) FSIindex = 3;// 20-30% Therm
- else if(Mbin<=7) FSIindex = 4;// 30-40% Therm
- else if(Mbin<=9) FSIindex = 5;// 40-50% Therm
- else if(Mbin<=11) FSIindex = 6;// Gauss R=3.5 fm
- else if(Mbin<=14) FSIindex = 7;// Gauss R=3.0 fm
- else if(Mbin<=16) FSIindex = 8;// Gauss R=2.5 fm
- else FSIindex = 9;// Gauss R=2.0 fm*/
-
if(CollisionType==0){
if(Mbin==0) FSIindex = 0;// 0-5% Therm
else if(Mbin==1) FSIindex = 1;// 0-10% Therm
double BinCenterPbPbPeripheral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509};
double BinCenterpPbAndpp[40]={0.00359275, 0.016357, 0.0257109, 0.035445, 0.045297, 0.0552251, 0.0651888, 0.0751397, 0.0851088, 0.0951108, 0.105084, 0.115079, 0.12507, 0.135059, 0.145053, 0.155049, 0.16505, 0.175038, 0.185039, 0.195034, 0.205023, 0.215027, 0.225024, 0.235023, 0.245011, 0.255017, 0.265017, 0.275021, 0.285021, 0.295017, 0.305018, 0.315018, 0.325013, 0.335011, 0.345016, 0.355019, 0.365012, 0.375016, 0.385017, 0.395016};
if(CollisionType==0){
- Q3Limit = 0.1 + 0.1*Mbin/16.;
for(int i=0; i<40; i++){
if(Mbin<10) BinCenters[i] = BinCenterPbPbCentral[i];
else BinCenters[i] = BinCenterPbPbPeripheral[i];
}
BinWidthQ2 = 0.005;
}else{
- Q3Limit = 0.3 + 0.2*fabs(Mbin-10)/9.;
for(int i=0; i<40; i++){
BinCenters[i] = BinCenterpPbAndpp[i];
}
BinWidthQ2 = 0.01;
}
- if(FitRangeShift) Q3Limit -= 0.3*Q3Limit;
- Q2Limit = Q3Limit/sqrt(2.);
- //Q2Limit = 1.2;
- //cout<<Q2Limit<<" "<<Q3Limit<<endl;
+
// extend BinCenters for high q
for(int index=40; index<400; index++){
//_file0 = new TFile("Results/PDC_10h_dEta0p025dPhi0p055.root","READ");
//_file0 = new TFile("Results/PDC_11h_3Kt3bins_FB5and7overlap.root","READ");
//_file0 = new TFile("Results/PDC_10h_11h_V0binning.root","READ");
+ //_file0 = new TFile("Results/PDC_10h_11h_lowerMultBins.root","READ");
_file0 = new TFile("Results/PDC_10h_11h_NclsFix.root","READ");// standard
}
}else if(CollisionType==1){// pPb
//_file0 = new TFile("Results/PDC_13e_kHighMult.root","READ");
//_file0 = new TFile("Results/PDC_13bcd_kINT7_plus_kHighMult.root","READ");// v5 and before
//_file0 = new TFile("Results/PDC_13bcde_kINT7_NclsFix.root","READ");
+ //_file0 = new TFile("Results/PDC_13e_kHighMult_NclsFix.root","READ");
_file0 = new TFile("Results/PDC_13bcde_kINT7_kHighMult_NclsFix.root","READ");// standard
}else{
_file0 = new TFile("Results/PDC_13b2_efix_p1234_R2_2Ktbins.root","READ");// standard (gen level)
_file0 = new TFile("Results/PDC_10cde_kMB_kHighMult_NclsFix.root","READ");// standard
}else{
_file0 = new TFile("Results/PDC_10f6a_R2_2Ktbins.root","READ");// standard
+ //_file0 = new TFile("Results/PDC_10f6a_R2_AOD161.root","READ");
//_file0 = new TFile("Results/PDC_10f6a_R10_genLevel_2Ktbins.root","READ");
}
}
gPad->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);
legend1->Draw("same");
}
- //for(int i=0; i<Two_ex_clone_mm->GetNbinsX(); i++){
- //cout<<Two_ex_clone_mm->GetBinError(i+1)<<", ";
+ //for(int i=0; i<100; i++){
+ //cout<<Two_ex_clone_mm->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};
//
// 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]);
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";
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];
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; i<npar; i++) {
fitC2ss_Expan->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_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;
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);
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<<C2_ss->GetBinContent(bin)<<", ";
+ //cout<<C2_ss->GetBinError(bin)<<", ";
+ //cout<<fitC2ss_h->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<<endl;
+ temp_hist_ss->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");
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);
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 (-+), <Nch>=36","p");
+ legend1->Draw("same");
}
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<<A_3[i-1][j-1][k-1]<<" "<<B_3[i-1][j-1][k-1]<<" "<<A_3_e[i-1][j-1][k-1]<<endl;
///////////////////////////////////////////////////////////
}
maxVal_c3[0] = 1.1; maxVal_c3[1] = 2.5; maxVal_c3[2] = 15.; maxVal_c3[3] = +1; maxVal_c3[4] = +1;// was minVal_c3[3] = +1; minVal_c3[4] = +1;
parName_c3[0] = "N"; parName_c3[1] = "#lambda_{3}"; parName_c3[2] = "R_{inv}"; parName_c3[3] = "coeff_{3}"; parName_c3[4] = "coeff_{4}";
c3Fit1D_Base->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; i<npar_c3; i++){
MyMinuit_c3.DefineParameter(i, parName_c3[i].c_str(), par_c3[i], stepSize_c3[i], minVal_c3[i], maxVal_c3[i]);
}
//
if(!FixExpansionAvg) {
maxVal_c3[1] = 2.0;
- double RadiusFix;
+ /*double RadiusFix;
if(CollisionType==0) RadiusFix = 10 - 7*Mbin/15.;
else RadiusFix = 2.5 - 1.25*(Mbin-12)/7.;
MyMinuit_c3.DefineParameter(2, parName_c3[2].c_str(), RadiusFix, stepSize_c3[2], minVal_c3[2], maxVal_c3[2]); MyMinuit_c3.FixParameter(2);
MyMinuit_c3.DefineParameter(1, parName_c3[1].c_str(), 2.0, stepSize_c3[1], minVal_c3[1], maxVal_c3[1]); MyMinuit_c3.FixParameter(1);
+ */
}
// kappa_3 and kappa_4
if(ft==0){// Gaussian
MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), 0., stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
}else{// Edgeworth
if(FixExpansionAvg){
- MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), kappa3, stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3);
- MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), kappa4, stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
- }else{
+ if(Gaussian){
+ MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), kappa3, stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3);
+ MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), kappa4, stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
+ }else{
+ MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), 0, stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3);
+ MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), 0, stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
+ }
+ }else{
Int_t err=0;
Double_t tmp[1];
tmp[0] = 3+1;
cout<<"End Three-d fit"<<endl;
/////////////////////////////////////////////////////////////
MyMinuit_c3.mnexcm("SHOw PARameters", &arglist_c3, 1, ierflg);
- cout<<"c3 Fit: Chi2/NDF = "<<Chi2_c3/(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
+ cout<<"c3 Fit: Chi2 = "<<Chi2_c3<<" NDF = "<<(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
for(int i=0; i<npar_c3; i++){
MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]);
c3Fit1D_Expan->SetParameter(3,OutputPar_c3[3]); c3Fit1D_Expan->SetParError(3,OutputPar_c3_e[3]);
c3Fit1D_Expan->SetParameter(4,OutputPar_c3[4]); c3Fit1D_Expan->SetParError(4,OutputPar_c3_e[4]);
c3Fit1D_Expan->SetLineStyle(1);
- c3Fit1D_Expan->Draw("same");
+ //c3Fit1D_Expan->Draw("same");
double lambdaStar=0, lambdaStar_e=0;
if(Gaussian){
lambdaStar = OutputPar_c3[1]*pow(1 + OutputPar_c3[4]/8.,3);
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);
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;
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<<c3_hist->GetBinContent(bin)<<" "<<c3Fit1D_Base->Eval(binCenter)<<" "<<c3_hist->GetBinError(bin)<<endl;
+ //cout<<c3_hist->GetBinContent(bin)<<" "<<ypoints[grIndex]<<" "<<c3_hist->GetBinError(bin)<<endl;
+
+ SumNDF_1D++;
+ }
+ cout<<"1D Chi2/NDF = "<<ChiSqSum_1D / (SumNDF_1D-3.)<<endl;
// uncomment to display fit box
//c3_hist->GetListOfFunctions()->Add(c3Fit1D_Base);
//cout<<endl;
- /*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,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,"<N_{rec,pions}>=103");
//fitBox3->SetNDC(kTRUE);
//fitBox3->Draw();
//for(int ii=0; ii<10; ii++) cout<<c3_hist->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
C3_Extended->SetMarkerColor(4);
c3_Extended->GetXaxis()->SetRangeUser(0,2);
c3_Extended->SetMarkerStyle(20);
- c3_Extended->Draw();*/
+ c3_Extended->Draw();
+ //cout<<c3_Extended->GetBinError(40)<<" "<<Three_1d[CB1][CB2][CB3][SCBin][0]->GetBinError(40)/Three_1d[CB1][CB2][CB3][SCBin][4]->GetBinContent(40)<<" "<<Three_1d[CB1][CB2][CB3][SCBin][1]->GetBinError(40)/Three_1d[CB1][CB2][CB3][SCBin][4]->GetBinContent(40)<<endl;
//legend4->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");
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");
+ */
- */
//////////////////////////////////////////////////////////////////////////////////////
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");
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
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;
lnL += B_3[i][j][k]*log(term2);
}else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
*/
- float DegenerateCount=1;
- if(i==j && i==k) DegenerateCount=1;
- else if(i==j || i==k || j==k) DegenerateCount=3;
- else DegenerateCount=6;
+ float DegenerateCount=1.;
+ //if(i==j && i==k) DegenerateCount=1;
+ //else if(i==j || i==k || j==k) DegenerateCount=3;//3
+ //else DegenerateCount=6;//6
NFitPoints_c3 += 1/DegenerateCount;
}
//f = -2.0*lnL;// log-liklihood minimization
f = SumChi2;// Chi2 minimization
Chi2_c3 = f;
+ //Chi2_c3 = SumChi2_test;
}
//________________________________________________________________________
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TSpline.h"
+#include "TMultiGraph.h"
#define BohrR 1963.6885
#define FmToGeV 0.19733 // conversion to fm
bool SaveFiles=kFALSE;
const int KT3Bin=0;// Kt3 bin. 0-1
-int FitType=1;// 0 (Gaussian), 1 (Edgeworth)
-bool p_pPb_Comp=0;
+int FitType=2;// 0 (Gaussian), 1 (Edgeworth), 2 (Exponential)
+bool pp_pPb_Comp=0;
bool AddedCC=kTRUE;// Charge Conjugate already added?
bool NchOneThirdAxis=0;
//
float SizeTitle=0.06;//
float SizeSpecif=0.05;//
-
+bool RadiusOnly=0;// only show radii, not lambdas
double RightMargin=0.002;// 0.002
//
TF1 *MixedChargeSysFit=new TF1("MixedChargeSysFit","[0]+[1]*exp(-pow([2]*x/0.19733,2))",0,.5);
//
- TF1 *Fit_C2[3][2][3][20];// CollType, Gauss/EW, EDbin, cb
- TH1D *Parameters_C2[3][2][3][6];// CollType, Gauss/EW, EDbin, Parameter#
+ TF1 *Fit_C2[3][3][3][20];// CollType, Gauss/EW/Exp, EDbin, cb
+ TH1D *Parameters_C2[3][3][3][6];// CollType, Gauss/EW/Exp, EDbin, Parameter#
TH1D *C2[3][3][20];// CollType, EDbin, cb
TH1D *C2_Sys[3][3][20];// CollType, EDbin, cb
TH1D *C3[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
TH1D *c3[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
TH1D *C3_Sys[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
TH1D *c3_Sys[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
- TF1 *c3_fit[3][2][3][20];// CollType, Gauss/EW, EDbin, cb
+ TF1 *c3_fit[3][3][3][20];// CollType, Gauss/EW/Exp, EDbin, cb
TGraph *gr_c3Spline[3][3][20];// CollType, EDbin, cb
+ TGraph *gr_c3SplineExpFit[3][3][20];// CollType, EDbin, cb
TF1 *c3_mixedChargeSysFit[3][2][20];// CollType, EDbin, cb
- TH1D *Parameters_c3[3][2][3][6];// CollType, Gaussian/EW, EDbin, Parameter#
-
+ TH1D *Parameters_c3[3][3][3][6];// CollType, Gaussian/EW, EDbin, Parameter#
+ TH1D *Parameters_Bjoern[2][3];// Bjoern's points: Hydro case, CollType
+
for(int ct=0; ct<3; ct++){
- 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++){
TString *name_C2=new TString("Parameters_C2_");
if(NchOneThirdAxis) {
Parameters_C2[ct][ft][kt3][par] = new TH1D(name_C2->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
+ }
}
}
}
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");
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);
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{
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);
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
//
}
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
}// 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){
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}");
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
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; i<Sys_K0s_C3->GetNbinsX(); 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<<K0s_C3->GetXaxis()->GetBinLowEdge(i+1)<<" "<<K0s_C3->GetXaxis()->GetBinUpEdge(i+1)<<" "<<K0s_C3->GetBinContent(i+1)<<" "<<K0s_C3->GetBinError(i+1)<<" "<<Sys_K0s_C3->GetBinError(i+1)<<endl;
+ //cout<<K0s_C3->GetXaxis()->GetBinLowEdge(i+1)<<" "<<K0s_C3->GetXaxis()->GetBinUpEdge(i+1)<<" "<<K0s_c3->GetBinContent(i+1)<<" "<<K0s_c3->GetBinError(i+1)<<" "<<Sys_K0s_c3->GetBinError(i+1)<<endl;
}
+
Sys_K0s_C3->Draw("E2 same");
Sys_K0s_c3->Draw("E2 same");
K0s_C3->Draw("same");
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 <dNch/deta> * 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
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
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();
//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};
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;}
}
- 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);
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");
//
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");
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();
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);
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
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]);
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);
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;}
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);
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");
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);
for(int cb=0; cb<20; cb++){
int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
if(RadiiPbPb->GetBinContent(binPbPb)==0) continue;
- cout<<"Nch="<<meanNchPbPb[cb]<<": R_3 = "<<RadiiPbPb->GetBinContent(binPbPb)<<" +- "<<RadiiPbPb->GetBinError(binPbPb)<<" +- "<<grRadiiSys_PbPb->GetErrorY(cb)<<" lambda_3 = "<<LambdaPbPb->GetBinContent(binPbPb)<<" +- "<<LambdaPbPb->GetBinError(binPbPb)<<" +- "<<grLambdaSys_PbPb->GetErrorY(cb)<<endl;
- cout<<" R_2 = "<<RadiiC2PbPb->GetBinContent(binPbPb)<<" +- "<<RadiiC2PbPb->GetBinError(binPbPb)<<" +- "<<grRadiiC2Sys_PbPb->GetErrorY(cb)<<" lambda_2 = "<<LambdaC2PbPb->GetBinContent(binPbPb)<<" +- "<<LambdaC2PbPb->GetBinError(binPbPb)<<" +- "<<grLambdaC2Sys_PbPb->GetErrorY(cb)<<endl;
+ //cout<<meanNchPbPb[cb]-(0.05)*meanNchPbPb[cb]<<" "<<meanNchPbPb[cb]+(0.05)*meanNchPbPb[cb]<<" "<<RadiiC2PbPb->GetBinContent(binPbPb)<<" "<<RadiiC2PbPb->GetBinError(binPbPb)<<" "<<grRadiiC2Sys_PbPb->GetErrorY(cb)<<" "<<RadiiPbPb->GetBinContent(binPbPb)<<" "<<RadiiPbPb->GetBinError(binPbPb)<<" "<<grRadiiSys_PbPb->GetErrorY(cb)<<endl;
+ //
+ //cout<<meanNchPbPb[cb]-(0.05)*meanNchPbPb[cb]<<" "<<meanNchPbPb[cb]+(0.05)*meanNchPbPb[cb]<<" "<<LambdaC2PbPb->GetBinContent(binPbPb)<<" "<<LambdaC2PbPb->GetBinError(binPbPb)<<" "<<grLambdaC2Sys_PbPb->GetErrorY(cb)<<" "<<LambdaPbPb->GetBinContent(binPbPb)<<" "<<LambdaPbPb->GetBinError(binPbPb)<<" "<<grLambdaSys_PbPb->GetErrorY(cb)<<endl;
}
cout<<"p--Pb:"<<endl;
for(int cb=0; cb<20; cb++){
int binpPb = RadiiC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]);
if(RadiipPb->GetBinContent(binpPb)==0) continue;
- cout<<"Nch="<<meanNchpPb[cb]<<": R_3 = "<<RadiipPb->GetBinContent(binpPb)<<" +- "<<RadiipPb->GetBinError(binpPb)<<" +- "<<grRadiiSys_pPb->GetErrorY(cb)<<" lambda_3 = "<<LambdapPb->GetBinContent(binpPb)<<" +- "<<LambdapPb->GetBinError(binpPb)<<" +- "<<grLambdaSys_pPb->GetErrorY(cb)<<endl;
- cout<<" R_2 = "<<RadiiC2pPb->GetBinContent(binpPb)<<" +- "<<RadiiC2pPb->GetBinError(binpPb)<<" +- "<<grRadiiC2Sys_pPb->GetErrorY(cb)<<" lambda_2 = "<<LambdaC2pPb->GetBinContent(binpPb)<<" +- "<<LambdaC2pPb->GetBinError(binpPb)<<" +- "<<grLambdaC2Sys_pPb->GetErrorY(cb)<<endl;
+ //cout<<meanNchpPb[cb]-(0.05)*meanNchpPb[cb]<<" "<<meanNchpPb[cb]+(0.05)*meanNchpPb[cb]<<" "<<RadiiC2pPb->GetBinContent(binpPb)<<" "<<RadiiC2pPb->GetBinError(binpPb)<<" "<<grRadiiC2Sys_pPb->GetErrorYhigh(cb)<<" "<<grRadiiC2Sys_pPb->GetErrorYlow(cb)<<" "<<RadiipPb->GetBinContent(binpPb)<<" "<<RadiipPb->GetBinError(binpPb)<<" "<<grRadiiSys_pPb->GetErrorY(cb)<<endl;
+ // Gaussian lambdas
+ //cout<<meanNchpPb[cb]-(0.05)*meanNchpPb[cb]<<" "<<meanNchpPb[cb]+(0.05)*meanNchpPb[cb]<<" "<<LambdaC2pPb->GetBinContent(binpPb)<<" "<<LambdaC2pPb->GetBinError(binpPb)<<" "<<grLambdaC2Sys_pPb->GetErrorY(cb)<<" "<<LambdapPb->GetBinContent(binpPb)<<" "<<LambdapPb->GetBinError(binpPb)<<" "<<grLambdaSys_pPb->GetErrorY(cb)<<endl;
+ // Edgeworth lambdas
+ //cout<<meanNchpPb[cb]-(0.05)*meanNchpPb[cb]<<" "<<meanNchpPb[cb]+(0.05)*meanNchpPb[cb]<<" "<<LambdaC2pPb->GetBinContent(binpPb)<<" "<<LambdaC2pPb->GetBinError(binpPb)<<" "<<grLambdaC2Sys_pPb->GetErrorYhigh(cb)<<" "<<grLambdaC2Sys_pPb->GetErrorYlow(cb)<<" "<<LambdapPb->GetBinContent(binpPb)<<" "<<LambdapPb->GetBinError(binpPb)<<" "<<grLambdaSys_pPb->GetErrorY(cb)<<endl;
}
cout<<"p--p:"<<endl;
for(int cb=0; cb<20; cb++){
int binpp = RadiiC2pp->GetXaxis()->FindBin(meanNchpp[cb]);
if(Radiipp->GetBinContent(binpp)==0) continue;
- cout<<"Nch="<<meanNchpp[cb]<<": R_3 = "<<Radiipp->GetBinContent(binpp)<<" +- "<<Radiipp->GetBinError(binpp)<<" +- "<<grRadiiSys_pp->GetErrorY(cb)<<" lambda_3 = "<<Lambdapp->GetBinContent(binpp)<<" +- "<<Lambdapp->GetBinError(binpp)<<" +- "<<grLambdaSys_pp->GetErrorY(cb)<<endl;
- cout<<" R_2 = "<<RadiiC2pp->GetBinContent(binpp)<<" +- "<<RadiiC2pp->GetBinError(binpp)<<" +- "<<grRadiiC2Sys_pp->GetErrorY(cb)<<" lambda_2 = "<<LambdaC2pp->GetBinContent(binpp)<<" +- "<<LambdaC2pp->GetBinError(binpp)<<" +- "<<grLambdaC2Sys_pp->GetErrorY(cb)<<endl;
- //
-
+ //cout<<meanNchpp[cb]-(0.05)*meanNchpp[cb]<<" "<<meanNchpp[cb]+(0.05)*meanNchpp[cb]<<" "<<RadiiC2pp->GetBinContent(binpp)<<" "<<RadiiC2pp->GetBinError(binpp)<<" "<<grRadiiC2Sys_pp->GetErrorYhigh(cb)<<" "<<grRadiiC2Sys_pp->GetErrorYlow(cb)<<" "<<Radiipp->GetBinContent(binpp)<<" "<<Radiipp->GetBinError(binpp)<<" "<<grRadiiSys_pp->GetErrorY(cb)<<endl;
+ // Gaussian lambdas
+ //cout<<meanNchpp[cb]-(0.05)*meanNchpp[cb]<<" "<<meanNchpp[cb]+(0.05)*meanNchpp[cb]<<" "<<LambdaC2pp->GetBinContent(binpp)<<" "<<LambdaC2pp->GetBinError(binpp)<<" "<<grLambdaC2Sys_pp->GetErrorY(cb)<<" "<<Lambdapp->GetBinContent(binpp)<<" "<<Lambdapp->GetBinError(binpp)<<" "<<grLambdaSys_pp->GetErrorY(cb)<<endl;
+ // Edgeworth lambdas
+ //cout<<meanNchpp[cb]-(0.05)*meanNchpp[cb]<<" "<<meanNchpp[cb]+(0.05)*meanNchpp[cb]<<" "<<LambdaC2pp->GetBinContent(binpp)<<" "<<LambdaC2pp->GetBinError(binpp)<<" "<<grLambdaC2Sys_pp->GetErrorYhigh(cb)<<" "<<grLambdaC2Sys_pp->GetErrorYlow(cb)<<" "<<Lambdapp->GetBinContent(binpp)<<" "<<Lambdapp->GetBinError(binpp)<<" "<<grLambdaSys_pp->GetErrorY(cb)<<endl;
}
pad3_2->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);
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);
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// 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
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<<Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinContent(bin)<<", ";
- cout<<Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin)<<", ";
- }else cout<<0<<", ";
+
//cout<<Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin)<<", ";
+ //}else cout<<0<<", ";
- }
- cout<<endl;
- }*/
+ //}
+ //cout<<endl;
+ //}
TF1 *Fit_kappa4_PbPb=new TF1("Fit_kappa4_PbPb","pol0",2,3000);
Combined_kappaPlot_1->Fit(Fit_kappa4_PbPb,"IMEN","",2,2000);
Fit_kappa4_PbPb->Draw("same");
-
+ */
////////////////////////////////////////////////////
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);
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++){
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<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinLowEdge(binN)<<" "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinUpEdge(binN)<<" "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<" "<<C3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<endl;
+ //
+ //cout<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinLowEdge(binN)<<" "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinUpEdge(binN)<<" "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<" "<<c3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<endl;
+ //
+ //if(System_proof==0){
+ //cout<<HIJING_c3_SC->GetXaxis()->GetBinLowEdge(binN)<<" "<<HIJING_c3_SC->GetXaxis()->GetBinUpEdge(binN)<<" "<<HIJING_c3_SC->GetBinContent(binN)<<" "<<HIJING_c3_SC->GetBinError(binN)<<endl;
+ //}else{
+ //cout<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinLowEdge(binN)<<" "<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinUpEdge(binN)<<" "<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" "<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<endl;
+ //}
+ //}
+ //cout<<endl;
C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->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){
}
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");
}
}
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);
//
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();
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");
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<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<" "<<c3_pPb->GetBinContent(binN)<<" "<<c3_pPb->GetBinError(binN)<<" "<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<" "<<c3_pp->GetBinContent(binN)<<" "<<c3_pp->GetBinError(binN)<<" "<<c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->GetBinError(binN)<<endl;
+ else cout<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<" "<<c3_pPb->GetBinContent(binN)<<" "<<c3_pPb->GetBinError(binN)<<" "<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<" "<<c3_PbPb->GetBinContent(binN)<<" "<<c3_PbPb->GetBinError(binN)<<" "<<c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->GetBinError(binN)<<endl;
+ }
//
if(SaveFiles) {
}
+
+
+
+
+
+
+
+
+
+
////////////////////////////////////////////////////
////////////////////////////////////////////////////
+ TCanvas *can6 = new TCanvas("can6", "can6",1700,700,600,600);// 11,53,700,500
+ can6->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();
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff