PWGCF/FEMTOSCOPY/macros/Plot_plotsTPR.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6
   7 #include "TVector2.h"
   8 #include "TFile.h"
   9 #include "TString.h"
  10 #include "TLatex.h"
  11 #include "TCanvas.h"
  12 #include "TPad.h"
  13 #include "TF1.h"
  14 #include "TH1.h"
  15 #include "TH2.h"
  16 #include "TH3.h"
  17 #include "TProfile.h"
  18 #include "TProfile2D.h"
  19 #include "TMath.h"
  20 #include "TText.h"
  21 #include "TRandom3.h"
  22 #include "TArray.h"
  23 #include "TLegend.h"
  24 #include "TStyle.h"
  25 #include "TMinuit.h"
  26 #include "TASImage.h"
  27 #include "TGraphErrors.h"
  28 #include "TGraphAsymmErrors.h"
  29 #include "TSpline.h"
  30 #include "TMultiGraph.h"
  31
  32 #define BohrR 1963.6885
  33 #define FmToGeV 0.19733 // conversion to fm
  34 #define PI 3.1415926
  35 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  36
  37 using namespace std;
  38
  39 bool SaveFiles=kFALSE;
  40 const int KT3Bin=0;// Kt3 bin. 0-1
  41 int FitType=1;// 0 (Gaussian), 1 (Edgeworth), 2 (Exponential)
  42 bool pp_pPb_Comp=0;
  43 bool AddedCC=kTRUE;// Charge Conjugate already added?
  44 bool NchOneThirdAxis=0;
  45 //
  46 int MbinMaxPbPb=15;// 15
  47 int MbinMinpPb=12;// 13
  48 int MbinMinpp=13;// 14
  49 int MbinMinPbPb=0;// 0
  50 int MbinMaxpPb=18;// 18
  51 //
  52 //
  53 //
  54 const int MaxKT3Bins=2;
  55 int TextFont=42;// 63, or 42
  56 float SizeLabel=0.06;// 20(63 font), 0.08(42 font)
  57 float SizeLegend=0.05;//
  58 float SizeTitle=0.06;//
  59 float SizeSpecif=0.05;//
  60
  61 bool RadiusOnly=0;// only show radii, not lambdas
  62
  63
  64 double RightMargin=0.002;// 0.002
  65 //
  66 double Chi2_C2global;
  67 double NFitPoints_C2global;
  68 TH1D *C2_ss[20];
  69 TH1D *C2_os[20];
  70
  71
  72
  73 void DrawALICELogo(Bool_t, Float_t, Float_t, Float_t, Float_t);
  74 TCanvas *make_canvas(const Char_t*,const Char_t*,Int_t,Int_t,Int_t,Int_t,Int_t);
  75
  76 void Plot_plotsTPR(){
  77
  78   gStyle->SetOptStat(0);
  79   gStyle->SetOptDate(0);
  80   //gStyle->SetOptFit(0111);
  81
  82   ////////////////////////////////////
  83   // Get Nrec to Nch mapping
  84   double meanNchPbPb[20]={0};
  85   //double meanNchPbPb_e[20]={0};
  86   double meanNchpPb[20]={0};
  87   //double meanNchpPb_e[20]={0};
  88   double meanNchpp[20]={0};
  89   //double meanNchpp_e[20]={0};
  90   TFile *NrecMapFile;
  91   TH2D *NrecMap;
  92   TList *MyList;
  93   //
  94   NrecMapFile = new TFile("Results/NrecMapping_12a17a_NclsFix.root","READ");// standard
  95   //NrecMapFile = new TFile("Results/NrecMapping_12a17a.root","READ");// v5 and before (with P < 1.0 cut)
  96   //NrecMapFile = new TFile("Results/NrecMapping_12a17a_TuneOnData.root","READ");
  97   //NrecMapFile = new TFile("Results/Old_NrecMappingFiles/NrecMapping_12a17a.root","READ");// paper v1 file (without P < 1.0 cut)
  98   //NrecMapFile = new TFile("Results/NrecMapping_12a11a.root","READ");// MC variation
  99   MyList=(TList*)NrecMapFile->Get("MyList");
 100   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");// Nch
 101   //NrecMap = (TH2D*)MyList->FindObject("fNpionTrueDist");// Npions
 102   for(int bin=1; bin<=2; bin++){// 1 to 2 (FB7),  1 to 1 (AMPT),  1 to 4 (FB5and7overlap)
 103     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 104     if(NrecMap->GetMean(2)>0) {
 105       meanNchPbPb[bin-1] = NrecMap->GetMean(2);
 106       //cout<<NrecMap->GetMean(2)<<"  "<<fabs(NrecMap->GetRMS(2))/NrecMap->GetMean(2)<<endl;
 107     }
 108     if(NchOneThirdAxis) if(NrecMap->GetMean(2)>0) meanNchPbPb[bin-1] = pow(NrecMap->GetMean(2),1/3.);
 109   }
 110   NrecMapFile->Close();
 111   //
 112   NrecMapFile = new TFile("Results/NrecMapping_12a17e_NclsFix.root","READ");// standard
 113   //NrecMapFile = new TFile("Results/NrecMapping_12a17e.root","READ");// v5 and before (with P < 1.0 cut)
 114   //NrecMapFile = new TFile("Results/NrecMapping_12a17e_TuneOnData.root","READ");
 115   //NrecMapFile = new TFile("Results/Old_NrecMappingFiles/NrecMapping_12a17e.root","READ");// paper v1 file (without P < 1.0 cut)
 116   //NrecMapFile = new TFile("Results/NrecMapping_12a11b.root","READ");// MC variation
 117   MyList=(TList*)NrecMapFile->Get("MyList");
 118   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");
 119   //NrecMap = (TH2D*)MyList->FindObject("fNpionTrueDist");
 120   for(int bin=3; bin<=10; bin++){// 3 to 10 (FB7),  2 to 3 (AMPT), 5 to 12 (FB5and7overlap)
 121     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 122     if(NrecMap->GetMean(2)>0) {
 123       meanNchPbPb[bin-1] = NrecMap->GetMean(2);
 124       //cout<<NrecMap->GetMean(2)<<"  "<<fabs(NrecMap->GetRMS(2))/NrecMap->GetMean(2)<<endl;
 125     }
 126     if(NchOneThirdAxis) if(NrecMap->GetMean(2)>0) meanNchPbPb[bin-1] = pow(NrecMap->GetMean(2),1/3.);
 127   }
 128   NrecMapFile->Close();
 129   //
 130   ////////////////////// extra for AMPT
 131   /*NrecMapFile = new TFile("Results/NrecMapping_12a11d.root","READ");
 132   MyList=(TList*)NrecMapFile->Get("MyList");
 133   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");
 134   for(int bin=4; bin<=7; bin++){// 4 to 7 (AMPT)
 135     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 136     if(NrecMap->GetMean(2)>0) meanNchPbPb[bin-1] = log10(NrecMap->GetMean(2));
 137   }
 138   NrecMapFile->Close();*/
 139   //////////////////////////////////////////
 140   //
 141   NrecMapFile = new TFile("Results/NrecMapping_12a17c_NclsFix.root","READ");// standard
 142   //NrecMapFile = new TFile("Results/NrecMapping_12a17c.root","READ");// v5 and before (with P < 1.0 cut)
 143   //NrecMapFile = new TFile("Results/NrecMapping_12a17c_TuneOnData.root","READ");
 144   //NrecMapFile = new TFile("Results/Old_NrecMappingFiles/NrecMapping_12a17c.root","READ");// paper v1 file (without P < 1.0 cut)
 145   //NrecMapFile = new TFile("Results/NrecMapping_12a11g.root","READ");// MC variation
 146   MyList=(TList*)NrecMapFile->Get("MyList");
 147   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");
 148   //NrecMap = (TH2D*)MyList->FindObject("fNpionTrueDist");
 149   for(int bin=11; bin<=19; bin++){// 11 to 19 (FB7),  1 to 1 (AMPT), 13 to 19 (FB5and7overlap)
 150     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 151     if(NrecMap->GetMean(2)>0) {
 152       meanNchPbPb[bin-1] = NrecMap->GetMean(2);
 153       //cout<<NrecMap->GetMean(2)<<"  "<<fabs(NrecMap->GetRMS(2))/NrecMap->GetMean(2)<<endl;
 154     }
 155     if(NchOneThirdAxis) if(NrecMap->GetMean(2)>0) meanNchPbPb[bin-1] = pow(NrecMap->GetMean(2),1/3.);
 156   }
 157   NrecMapFile->Close();
 158   ///cout<<endl;
 159   //
 160   NrecMapFile = new TFile("Results/NrecMapping_13b2_efix_NclsFix.root","READ");// standard
 161   //NrecMapFile = new TFile("Results/NrecMapping_13b2_efix_p1.root","READ");// v5 and before (with P < 1.0 cut)
 162   //NrecMapFile = new TFile("Results/NrecMapping_13b2_TuneOnData.root","READ");
 163   //NrecMapFile = new TFile("Results/PDC_13b2_efix_p1_R2.root","READ");// paper v1 file (without P < 1.0 cut)
 164   //NrecMapFile = new TFile("Results/NrecMapping_13b3.root","READ");// MC variation
 165   MyList=(TList*)NrecMapFile->Get("MyList");
 166   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");
 167   //NrecMap = (TH2D*)MyList->FindObject("fNpionTrueDist");
 168   for(int bin=10; bin<=20; bin++){
 169     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 170     if(NrecMap->GetMean(2)>0) {
 171       meanNchpPb[bin-1] = NrecMap->GetMean(2);
 172       //cout<<NrecMap->GetMean(2)<<"  "<<fabs(NrecMap->GetRMS(2))/NrecMap->GetMean(2)<<endl;
 173     }
 174     if(NchOneThirdAxis) if(NrecMap->GetMean(2)>0) meanNchpPb[bin-1] = pow(NrecMap->GetMean(2),1/3.);
 175   }
 176   NrecMapFile->Close();
 177   //cout<<endl;
 178   //
 179   NrecMapFile = new TFile("Results/NrecMapping_10f6a_NclsFix.root","READ");// standard
 180   //NrecMapFile = new TFile("Results/NrecMapping_10f6a.root","READ");// v5 (with P < 1.0 cut)
 181   //NrecMapFile = new TFile("Results/NrecMapping_10f6a_TuneOnData.root","READ");
 182   //NrecMapFile = new TFile("Results/PDC_10f6a_R2.root","READ");// paper v1 file (without P < 1.0 cut)
 183   //NrecMapFile = new TFile("Results/NrecMapping_10f6.root","READ");// MC variation
 184   MyList=(TList*)NrecMapFile->Get("MyList");
 185   NrecMap = (TH2D*)MyList->FindObject("fNchTrueDist");
 186   //NrecMap = (TH2D*)MyList->FindObject("fNpionTrueDist");
 187   for(int bin=10; bin<=20; bin++){
 188     NrecMap->GetXaxis()->SetRangeUser(bin,bin);
 189     if(NrecMap->GetMean(2)>0) {
 190       meanNchpp[bin-1] = NrecMap->GetMean(2);
 191       //cout<<NrecMap->GetMean(2)<<"  "<<fabs(NrecMap->GetRMS(2))/NrecMap->GetMean(2)<<endl;
 192     }
 193     if(NchOneThirdAxis) if(NrecMap->GetMean(2)>0) meanNchpp[bin-1] = pow(NrecMap->GetMean(2),1/3.);
 194   }
 195   NrecMapFile->Close();
 196   //cout<<endl;
 197
 198   //for(int i=0; i<20; i++) cout<<pow(10,meanNchPbPb[i])<<endl;
 199   //cout<<"+++++++++++++++"<<endl;
 200   //for(int i=0; i<20; i++) cout<<pow(10,meanNchpPb[i])<<endl;
 201   //cout<<"+++++++++++++++"<<endl;
 202   //for(int i=0; i<20; i++) cout<<meanNchpp[i]<<endl;
 203
 204   TFile *ExRangeFile=new TFile("Results/ExtendedQ3rangeM0.root","READ");
 205   TList *ExList=(TList*)ExRangeFile->Get("MyList");
 206   TH1D *ExRangeTerm1=(TH1D*)ExList->FindObject("fExtendedQ3Histo_term1");
 207   TH1D *ExRangeTerm2=(TH1D*)ExList->FindObject("fExtendedQ3Histo_term2");
 208   TH1D *ExRangeTerm5=(TH1D*)ExList->FindObject("fExtendedQ3Histo_term5");
 209   ExRangeTerm1->SetDirectory(0); ExRangeTerm2->SetDirectory(0); ExRangeTerm5->SetDirectory(0);
 210   ExRangeFile->Close();
 211   ExRangeTerm1->Sumw2(); ExRangeTerm2->Sumw2(); ExRangeTerm5->Sumw2();
 212   ExRangeTerm2->Scale(ExRangeTerm1->Integral(25,29)/ExRangeTerm2->Integral(25,29));
 213   ExRangeTerm5->Scale(ExRangeTerm1->Integral(25,29)/ExRangeTerm5->Integral(25,29));
 214   float TwoFrac=0.7;
 215   float OneFrac=pow(TwoFrac,.5);
 216   float ThreeFrac=pow(TwoFrac,1.5);
 217   // Purify.  Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
 218   ExRangeTerm1->Add(ExRangeTerm5, -(pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2)));
 219   ExRangeTerm1->Add(ExRangeTerm2, -3*(1-OneFrac));
 220   ExRangeTerm1->Add(ExRangeTerm5, (1-OneFrac)*3*(1-TwoFrac));
 221   ExRangeTerm1->Scale(1/ThreeFrac);
 222   // Isolate 3-pion cumulant
 223   ExRangeTerm1->Add(ExRangeTerm2, -3/TwoFrac);
 224   ExRangeTerm1->Add(ExRangeTerm5, 3*(1-TwoFrac)/TwoFrac);
 225   ExRangeTerm1->Add(ExRangeTerm5, 3);
 226   ExRangeTerm1->Divide(ExRangeTerm5);
 227   ExRangeTerm1->SetMarkerStyle(20);
 228   ExRangeTerm1->SetMarkerColor(1);
 229
 230   //
 231   TF1 *MixedChargeSysFit=new TF1("MixedChargeSysFit","[0]+[1]*exp(-pow([2]*x/0.19733,2))",0,.5);
 232   //
 233   TF1 *Fit_C2[3][3][3][20];// CollType, Gauss/EW/Exp, EDbin, cb
 234   TH1D *Parameters_C2[3][3][3][6];// CollType, Gauss/EW/Exp, EDbin, Parameter#
 235   TH1D *C2[3][3][20];// CollType, EDbin, cb
 236   TH1D *C2_Sys[3][3][20];// CollType, EDbin, cb
 237   TH1D *C3[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
 238   TH1D *c3[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
 239   TH1D *C3_Sys[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
 240   TH1D *c3_Sys[3][2][2][3][20];// CollType, Real/MonteCarlo, SC/MC, EDbin, cb
 241   TF1 *c3_fit[3][3][3][20];// CollType, Gauss/EW/Exp, EDbin, cb
 242   TGraph *gr_c3Spline[3][3][20];// CollType, EDbin, cb
 243   TGraph *gr_c3SplineExpFit[3][3][20];// CollType, EDbin, cb
 244   TF1 *c3_mixedChargeSysFit[3][2][20];// CollType, EDbin, cb
 245   TH1D *Parameters_c3[3][3][3][6];// CollType, Gaussian/EW, EDbin, Parameter#
 246   TH1D *Parameters_Bjoern[2][3];// Bjoern's points: Hydro case, CollType
 247
 248   for(int ct=0; ct<3; ct++){
 249     for(int ft=0; ft<3; ft++){// Gaussian or EW or Exp
 250       for(int kt3=0; kt3<3; kt3++){
 251         for(int par=0; par<6; par++){
 252           TString *name_C2=new TString("Parameters_C2_");
 253           *name_C2 += ct;
 254           *name_C2 += ft;
 255           *name_C2 += kt3;
 256           *name_C2 += par;
 257           TString *name_c3=new TString("Parameters_c3_");
 258           *name_c3 += ct;
 259           *name_c3 += ft;
 260           *name_c3 += kt3;
 261           *name_c3 += par;
 262
 263           if(NchOneThirdAxis) {
 264             Parameters_C2[ct][ft][kt3][par] = new TH1D(name_C2->Data(),"",3000,1,13.5);// Nch^(1/3)
 265             Parameters_c3[ct][ft][kt3][par] = new TH1D(name_c3->Data(),"",3000,1,13.5);// Nch^(1/3)
 266             if(ft==0 && kt3==0 && par==0){
 267               TString *name_Bjoern = new TString("Bjoern_"); *name_Bjoern += ct;
 268               Parameters_Bjoern[0][ct] = new TH1D(name_Bjoern->Data(),"",3000,1,13.5);// Nch^(1/3)
 269               name_Bjoern->Append("_hydro");
 270               Parameters_Bjoern[1][ct] = new TH1D(name_Bjoern->Data(),"",3000,1,13.5);// Nch^(1/3)
 271             }
 272           }else{
 273             Parameters_C2[ct][ft][kt3][par] = new TH1D(name_C2->Data(),"",30000,1,3001);// Nch
 274             Parameters_c3[ct][ft][kt3][par] = new TH1D(name_c3->Data(),"",30000,1,3001);// Nch
 275             if(ft==0 && kt3==0 && par==0){
 276               TString *name_Bjoern = new TString("Bjoern_"); *name_Bjoern += ct;
 277               Parameters_Bjoern[0][ct] = new TH1D(name_Bjoern->Data(),"",30000,1,3001);// Nch
 278               name_Bjoern->Append("_hydro");
 279               Parameters_Bjoern[1][ct] = new TH1D(name_Bjoern->Data(),"",30000,1,3001);// Nch
 280             }
 281           }
 282         }
 283       }
 284     }
 285   }
 286
 287   TH1D *RadiiC2pp_Published = new TH1D("RadiiC2pp_Published","",3000,1.0,3001);
 288   //
 289   double N_1 = 0, N_1_e=0;
 290   double lambda_1 = 0, lambda_1_e=0;
 291   double radius_1 = 0, radius_1_e=0;
 292   double EW1_1 = 0, EW1_1_e=0;
 293   double EW2_1 = 0, EW2_1_e=0;
 294
 295
 296   // Start File access
 297   for(int ct=0; ct<3; ct++){
 298     for(int dt=0; dt<2; dt++){// data type (Real or Monte Carlo)
 299       if(ct==0 && dt==1) continue; // no MC for PbPb
 300       for(int cb=0; cb<20; cb++){
 301         if(ct==0 && cb<MbinMinPbPb) continue;
 302         if(ct==0 && cb>MbinMaxPbPb) continue;
 303         if(ct==1 && cb>MbinMaxpPb) continue;
 304         if(ct==1 && cb<MbinMinpPb) continue;
 305         if(ct==2 && cb<MbinMinpp) continue;
 306         if(ct==2 && dt==1 && cb<=13) continue;// no Pythia data for cb=13
 307         for(int ChComb=0; ChComb<2; ChComb++) {// SC or MC
 308           for(int ch=0; ch<1; ch++) {// - or +
 309             for(int KT3=0; KT3<MaxKT3Bins; KT3++) {// Kt3 bin
 310               if(dt==1 && KT3>0) continue;// no MC data yet for higher kt3
 311               TString *name3 = new TString("OutFiles/OutFile");
 312               if(ct==0) name3->Append("PbPb");
 313               if(ct==1) name3->Append("pPb");
 314               if(ct==2) name3->Append("pp");
 315               if(dt==1) name3->Append("MonteCarlo");
 316               if(ChComb==0) name3->Append("SC");
 317               else name3->Append("MC");
 318               if(ch==0) name3->Append("Neg");
 319               else name3->Append("Pos");
 320
 321               name3->Append("Kt3_");
 322               *name3 += KT3+1;
 323               name3->Append("_M");
 324               *name3 += cb;
 325               name3->Append(".root");
 326
 327               TFile *file=new TFile(name3->Data(),"READ");
 328
 329               TString *name3_2 = new TString("OutFiles/OutFileExpFit");
 330               if(ct==0) name3_2->Append("PbPb");
 331               if(ct==1) name3_2->Append("pPb");
 332               if(ct==2) name3_2->Append("pp");
 333               name3_2->Append("SC");
 334               if(ch==0) name3_2->Append("Neg");
 335               else name3_2->Append("Pos");
 336
 337               name3_2->Append("Kt3_");
 338               *name3_2 += KT3+1;
 339               name3_2->Append("_M");
 340               *name3_2 += cb;
 341               name3_2->Append(".root");
 342
 343               TFile *fileExpFit=new TFile(name3_2->Data(),"READ");
 344               //
 345
 346
 347               if(ch==0) {
 348                 C3[ct][dt][ChComb][KT3][cb]=(TH1D*)file->Get("C3");
 349                 C3[ct][dt][ChComb][KT3][cb]->SetDirectory(0);
 350                 C3[ct][dt][ChComb][KT3][cb]->SetMarkerStyle(24+ct);
 351                 c3[ct][dt][ChComb][KT3][cb]=(TH1D*)file->Get("c3");
 352                 c3[ct][dt][ChComb][KT3][cb]->SetDirectory(0);
 353                 c3[ct][dt][ChComb][KT3][cb]->SetMarkerStyle(20+ct);
 354                 if(dt==1) {c3[ct][dt][ChComb][KT3][cb]->SetMarkerStyle(28);}
 355                 C3[ct][dt][ChComb][KT3][cb]->GetXaxis()->SetRangeUser(0,0.5);
 356
 357                 C3[ct][dt][ChComb][KT3][cb]->GetXaxis()->SetRangeUser(0,0.5);
 358                 C3[ct][dt][ChComb][KT3][cb]->GetXaxis()->SetLabelFont(TextFont); C3[ct][dt][ChComb][KT3][cb]->GetYaxis()->SetLabelFont(TextFont);
 359                 C3[ct][dt][ChComb][KT3][cb]->GetYaxis()->SetTitleFont(TextFont);
 360                 c3[ct][dt][ChComb][KT3][cb]->GetXaxis()->SetRangeUser(0,0.5);
 361                 c3[ct][dt][ChComb][KT3][cb]->GetXaxis()->SetLabelFont(TextFont); c3[ct][dt][ChComb][KT3][cb]->GetYaxis()->SetLabelFont(TextFont);
 362                 c3[ct][dt][ChComb][KT3][cb]->GetYaxis()->SetTitleFont(TextFont);
 363                 if(ct==0) {
 364                   C3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(1); C3[ct][dt][ChComb][KT3][cb]->SetLineColor(1);
 365                   c3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(1); c3[ct][dt][ChComb][KT3][cb]->SetLineColor(1);
 366                 }else if(ct==1){
 367                   C3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(2); C3[ct][dt][ChComb][KT3][cb]->SetLineColor(2);
 368                   c3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(2); c3[ct][dt][ChComb][KT3][cb]->SetLineColor(2);
 369                   if(dt==1) {c3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(1); c3[ct][dt][ChComb][KT3][cb]->SetLineColor(1);}
 370                 }else {
 371                   C3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(4); C3[ct][dt][ChComb][KT3][cb]->SetLineColor(4);
 372                   c3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(4); c3[ct][dt][ChComb][KT3][cb]->SetLineColor(4);
 373                   if(dt==1) {c3[ct][dt][ChComb][KT3][cb]->SetMarkerColor(1); c3[ct][dt][ChComb][KT3][cb]->SetLineColor(1);}
 374                 }
 375
 376                 if(dt==0) {
 377                   if(ChComb==0) {
 378                     C2[ct][KT3][cb]=(TH1D*)file->Get("C2_ss");
 379                     Fit_C2[ct][0][KT3][cb]=(TF1*)file->Get("fitC2ss_Base");// was fitC2ss_Gauss
 380                     Fit_C2[ct][1][KT3][cb]=(TF1*)file->Get("fitC2ss_Expan");// fitC2ss_EW
 381                     Fit_C2[ct][2][KT3][cb]=(TF1*)fileExpFit->Get("fitC2ss_Base");// Exp fit
 382
 383                     C2_Sys[ct][KT3][cb]=(TH1D*)C2[ct][KT3][cb]->Clone();
 384                     if(ct==0){
 385                       C2[ct][KT3][cb]->SetMarkerColor(1); C2[ct][KT3][cb]->SetLineColor(1);
 386                       C2_Sys[ct][KT3][cb]->SetFillColor(kGray);
 387                     }else if(ct==1){
 388                       C2[ct][KT3][cb]->SetMarkerColor(2); C2[ct][KT3][cb]->SetLineColor(2);
 389                       C2_Sys[ct][KT3][cb]->SetFillColor(kRed-10);
 390                     }else{
 391                       C2[ct][KT3][cb]->SetMarkerColor(4); C2[ct][KT3][cb]->SetLineColor(4);
 392                       C2_Sys[ct][KT3][cb]->SetFillColor(kBlue-10);
 393                     }
 394                     C2_Sys[ct][KT3][cb]->SetMarkerSize(0); C2_Sys[ct][KT3][cb]->SetFillStyle(1000);
 395                     // C2 Systematics
 396                     for(int bin=1; bin<=C2_Sys[ct][KT3][cb]->GetNbinsX(); bin++){
 397                       C2_Sys[ct][KT3][cb]->SetBinError(bin, 0.01*C2_Sys[ct][KT3][cb]->GetBinContent(bin));
 398                     }
 399                     C2[ct][KT3][cb]->SetDirectory(0); //Fit_h_C2[ct][KT3][cb]->SetDirectory(0);
 400                     C2_Sys[ct][KT3][cb]->SetDirectory(0);
 401                     C2[ct][KT3][cb]->SetMarkerStyle(24+ct);
 402                     C2[ct][KT3][cb]->GetXaxis()->SetRangeUser(0,0.33);
 403                     C2[ct][KT3][cb]->GetXaxis()->SetLabelFont(TextFont); C2[ct][KT3][cb]->GetYaxis()->SetLabelFont(TextFont);
 404                     C2[ct][KT3][cb]->GetXaxis()->SetTitleOffset(1.2); C2[ct][KT3][cb]->GetYaxis()->SetTitleOffset(1.2);
 405                     C2[ct][KT3][cb]->GetXaxis()->SetTitleFont(TextFont); C2[ct][KT3][cb]->GetYaxis()->SetTitleFont(TextFont);
 406                     C2[ct][KT3][cb]->GetXaxis()->SetTitleSize(SizeTitle); //C2[ct][KT3][cb]->GetYaxis()->SetTitleFont(SizeTitle*SF2);
 407                     //
 408
 409                     c3_fit[ct][0][KT3][cb]=(TF1*)file->Get("c3Fit1D_Base");// was c3Fit1D_Gauss
 410                     c3_fit[ct][1][KT3][cb]=(TF1*)file->Get("c3Fit1D_Expan");// was c3Fit1D_EW
 411                     c3_fit[ct][2][KT3][cb]=(TF1*)fileExpFit->Get("c3Fit1D_Base");
 412                     gr_c3Spline[ct][KT3][cb] = (TGraph*)file->Get("gr_c3Spline");// Spline of a spline + TF1
 413                     gr_c3SplineExpFit[ct][KT3][cb] = (TGraph*)fileExpFit->Get("gr_c3Spline");// Exp fit
 414                     c3_fit[ct][0][KT3][cb]->SetLineStyle(6);
 415                     gr_c3Spline[ct][KT3][cb]->SetLineStyle(7); c3_fit[ct][1][KT3][cb]->SetLineStyle(7);
 416                     gr_c3SplineExpFit[ct][KT3][cb]->SetLineStyle(1); c3_fit[ct][2][KT3][cb]->SetLineStyle(1);
 417                     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);}
 418                     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);}
 419                     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);}
 420                   }// ChComb==0
 421
 422                 }
 423               }else{
 424
 425                 if(ChComb==0){
 426                   N_1 = 0; N_1_e=0;
 427                   lambda_1 = 0; lambda_1_e=0;
 428                   radius_1 = 0; radius_1_e=0;
 429                   EW1_1 = 0; EW1_1_e=0;
 430                   EW2_1 = 0; EW2_1_e=0;
 431                   //
 432                   if(!AddedCC){
 433                     cout<<"Not Supported!!"<<endl;
 434                   }
 435                 }
 436                 //
 437                 if(!AddedCC){
 438                   cout<<"Not Supported!!"<<endl;
 439                 }
 440               }// ch==1
 441
 442
 443               for(int bin=1; bin<10; bin++){// Remove large error bins
 444                 if(C3[ct][dt][ChComb][KT3][cb]->GetBinError(bin) > 0.33*C3[ct][dt][ChComb][KT3][cb]->GetBinContent(bin)){
 445                   C3[ct][dt][ChComb][KT3][cb]->SetBinContent(bin,10); C3[ct][dt][ChComb][KT3][cb]->SetBinError(bin,10);
 446                 }
 447                 if(c3[ct][dt][ChComb][KT3][cb]->GetBinError(bin) > 0.33*c3[ct][dt][ChComb][KT3][cb]->GetBinContent(bin)){
 448                   c3[ct][dt][ChComb][KT3][cb]->SetBinContent(bin,10); c3[ct][dt][ChComb][KT3][cb]->SetBinError(bin,10);
 449                 }
 450
 451               }
 452
 453               if(AddedCC && dt==0){
 454                 if(ct==0 || ct==1) c3[ct][dt][ChComb][KT3][cb]->SetMarkerSize(1.12*C3[ct][dt][ChComb][KT3][cb]->GetMarkerSize());
 455                 else c3[ct][dt][ChComb][KT3][cb]->SetMarkerSize(1.2*C3[ct][dt][ChComb][KT3][cb]->GetMarkerSize());
 456
 457                 //
 458                 if(ChComb==0){
 459                   double logNch=0;
 460                   if(ct==0) logNch=meanNchPbPb[cb];
 461                   else if(ct==1) logNch=meanNchpPb[cb];
 462                   else logNch=meanNchpp[cb];
 463                   int logNchBin = Parameters_c3[ct][0][KT3][0]->GetXaxis()->FindBin(logNch);
 464                   for(int ft=0; ft<3; ft++){// Gaussian or EW or Exp
 465                     N_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(0); N_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(0);
 466                     lambda_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(1); lambda_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(1);
 467                     radius_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(2); radius_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(2);
 468                     if(ft==2) {radius_1 /= sqrt(TMath::Pi()); radius_1_e /= sqrt(TMath::Pi());}
 469                     EW1_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(3); EW1_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(3);
 470                     EW2_1 = c3_fit[ct][ft][KT3][cb]->GetParameter(4); EW2_1_e = c3_fit[ct][ft][KT3][cb]->GetParError(4);
 471
 472                     if(ft!=1) {EW1_1=0; EW2_1=0; EW1_1_e=0; EW2_1_e=0;}// make sure they are zero
 473                     //
 474                     Parameters_c3[ct][ft][KT3][0]->SetBinContent(logNchBin, N_1); Parameters_c3[ct][ft][KT3][0]->SetBinError(logNchBin, N_1_e);
 475                     Parameters_c3[ct][ft][KT3][1]->SetBinContent(logNchBin, lambda_1); Parameters_c3[ct][ft][KT3][1]->SetBinError(logNchBin, lambda_1_e);
 476                     Parameters_c3[ct][ft][KT3][2]->SetBinContent(logNchBin, radius_1); Parameters_c3[ct][ft][KT3][2]->SetBinError(logNchBin, radius_1_e);
 477                     Parameters_c3[ct][ft][KT3][3]->SetBinContent(logNchBin, EW1_1); Parameters_c3[ct][ft][KT3][3]->SetBinError(logNchBin, EW1_1_e);
 478                     Parameters_c3[ct][ft][KT3][4]->SetBinContent(logNchBin, EW2_1); Parameters_c3[ct][ft][KT3][4]->SetBinError(logNchBin, EW2_1_e);
 479                     // lambda_3* parameter
 480                     Parameters_c3[ct][ft][KT3][5]->SetBinContent(logNchBin, lambda_1*pow(1 + EW2_1/8.,3));
 481                     Parameters_c3[ct][ft][KT3][5]->SetBinError(logNchBin, lambda_1_e*pow(1 + EW2_1/8.,3));
 482                     // remove unstable c3 Fit points
 483                     bool badbin=kFALSE;
 484                     if(ct==0 && cb>12) badbin=kTRUE;
 485                     if(ct==1 && cb<12) badbin=kTRUE;
 486                     if(ct==2 && cb<14) badbin=kTRUE;
 487                     if(badbin){
 488                       Parameters_c3[ct][ft][KT3][0]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][0]->SetBinError(logNchBin, 100);
 489                       Parameters_c3[ct][ft][KT3][1]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][1]->SetBinError(logNchBin, 100);
 490                       Parameters_c3[ct][ft][KT3][2]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][2]->SetBinError(logNchBin, 100);
 491                       Parameters_c3[ct][ft][KT3][3]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][3]->SetBinError(logNchBin, 100);
 492                       Parameters_c3[ct][ft][KT3][4]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][4]->SetBinError(logNchBin, 100);
 493                       Parameters_c3[ct][ft][KT3][5]->SetBinContent(logNchBin, 100); Parameters_c3[ct][ft][KT3][5]->SetBinError(logNchBin, 100);
 494                     }
 495                     //
 496                     Parameters_C2[ct][ft][KT3][0]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(0));// N
 497                     Parameters_C2[ct][ft][KT3][0]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(0));
 498                     Parameters_C2[ct][ft][KT3][1]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(1));// lambda
 499                     Parameters_C2[ct][ft][KT3][1]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(1));
 500                     radius_1 = Fit_C2[ct][ft][KT3][cb]->GetParameter(3); radius_1_e = Fit_C2[ct][ft][KT3][cb]->GetParError(3);
 501                     if(ft==2) {radius_1 /= sqrt(TMath::Pi()); radius_1_e /= sqrt(TMath::Pi());}
 502                     Parameters_C2[ct][ft][KT3][2]->SetBinContent(logNchBin, radius_1);// R
 503                     Parameters_C2[ct][ft][KT3][2]->SetBinError(logNchBin, radius_1_e);
 504                     Parameters_C2[ct][ft][KT3][3]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(5));// kappa3
 505                     Parameters_C2[ct][ft][KT3][3]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(5));
 506                     Parameters_C2[ct][ft][KT3][4]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(6));// kappa4
 507                     Parameters_C2[ct][ft][KT3][4]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(6));
 508                     // lambda_* parameter
 509                     Parameters_C2[ct][ft][KT3][5]->SetBinContent(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParameter(1)*pow(1 + EW2_1/8.,2));
 510                     Parameters_C2[ct][ft][KT3][5]->SetBinError(logNchBin, Fit_C2[ct][ft][KT3][cb]->GetParError(1)*pow(1 + EW2_1/8.,2));
 511                   }// ft
 512                 }// ChComb==0
 513                 //
 514                 // Sys errors
 515                 C3_Sys[ct][0][ChComb][KT3][cb] = (TH1D*)C3[ct][0][ChComb][KT3][cb]->Clone();
 516                 c3_Sys[ct][0][ChComb][KT3][cb] = (TH1D*)c3[ct][0][ChComb][KT3][cb]->Clone();
 517                 if(ct==0){
 518                   C3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kGray);
 519                   c3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kGray);
 520                 }else if(ct==1){
 521                   C3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kRed-10);
 522                   c3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kRed-10);
 523                 }else {
 524                   C3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); C3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kBlue-10);
 525                   c3_Sys[ct][0][ChComb][KT3][cb]->SetMarkerSize(0); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillStyle(1000); c3_Sys[ct][0][ChComb][KT3][cb]->SetFillColor(kBlue-10);
 526                 }
 527
 528                 if(ChComb==1){
 529                   MixedChargeSysFit->SetParameter(0,1);
 530                   MixedChargeSysFit->SetParameter(1,.1);
 531                   MixedChargeSysFit->SetParameter(2,1);
 532                   c3[ct][0][ChComb][KT3][cb]->Fit(MixedChargeSysFit,"IMNQ","",0.01,0.5);
 533                   c3_mixedChargeSysFit[ct][KT3][cb] = (TF1*)MixedChargeSysFit->Clone();
 534                   for(int i=1; i<=c3[ct][0][ChComb][KT3][cb]->GetNbinsX(); i++) {
 535                     float Q3=(i-0.5)*0.01;
 536                     // SameCharge
 537                     C3_Sys[ct][0][0][KT3][cb]->SetBinError(i, 0.01 * C3_Sys[ct][0][0][KT3][cb]->GetBinContent(i));
 538                     c3_Sys[ct][0][0][KT3][cb]->SetBinError(i, sqrt(pow(MixedChargeSysFit->Eval(Q3)-1.0,2) + pow(0.1*(c3_Sys[ct][0][0][KT3][cb]->GetBinContent(i)-1.0),2)));// residue + lambda undilution variation (0.7 to 0.65)
 539                     // MixedCharge
 540                     C3_Sys[ct][0][1][KT3][cb]->SetBinError(i, 0.01 * C3_Sys[ct][0][1][KT3][cb]->GetBinContent(i));// correlation function uncertainty
 541                     c3_Sys[ct][0][1][KT3][cb]->SetBinError(i, sqrt(pow(0.01 * c3_Sys[ct][0][1][KT3][cb]->GetBinContent(i),2) + pow(0.1*(c3_Sys[ct][0][1][KT3][cb]->GetBinContent(i)-1.0),2)));// correlation function uncertainty + lambda undilution variation (0.7 to 0.65)
 542                   }
 543                 }
 544                 C3_Sys[ct][0][ChComb][KT3][cb]->SetDirectory(0); c3_Sys[ct][0][ChComb][KT3][cb]->SetDirectory(0);
 545               }// AddedCC and dt==0
 546
 547               file->Close();
 548               fileExpFit->Close();
 549
 550
 551             }// Kt3
 552           }// ch
 553         }// ChComb
 554
 555       }// cb
 556
 557       if(dt==0){
 558         for(int ft=0; ft<3; ft++){// Gaussian or EW or Exp
 559           for(int KT3=0; KT3<3; KT3++){
 560             for(int par=0; par<6; par++){
 561               if(ct<2){
 562                 Parameters_C2[ct][ft][KT3][par]->SetMarkerStyle(24+ct);
 563                 Parameters_c3[ct][ft][KT3][par]->SetMarkerStyle(20+ct);
 564               }else{
 565                 Parameters_C2[ct][ft][KT3][par]->SetMarkerStyle(28);
 566                 Parameters_c3[ct][ft][KT3][par]->SetMarkerStyle(34);
 567                 RadiiC2pp_Published->SetMarkerStyle(28);
 568               }
 569
 570               if(ct==0){
 571                 Parameters_C2[ct][ft][KT3][par]->SetMarkerColor(1);
 572                 Parameters_C2[ct][ft][KT3][par]->SetLineColor(1);
 573                 Parameters_c3[ct][ft][KT3][par]->SetMarkerColor(1);
 574                 Parameters_c3[ct][ft][KT3][par]->SetLineColor(1);
 575               }else if(ct==1){
 576                 Parameters_C2[ct][ft][KT3][par]->SetMarkerColor(2);
 577                 Parameters_C2[ct][ft][KT3][par]->SetLineColor(2);
 578                 Parameters_c3[ct][ft][KT3][par]->SetMarkerColor(2);
 579                 Parameters_c3[ct][ft][KT3][par]->SetLineColor(2);
 580               }else {
 581                 Parameters_C2[ct][ft][KT3][par]->SetMarkerColor(4);
 582                 Parameters_C2[ct][ft][KT3][par]->SetLineColor(4);
 583                 Parameters_c3[ct][ft][KT3][par]->SetMarkerColor(4);
 584                 Parameters_c3[ct][ft][KT3][par]->SetLineColor(4);
 585                 RadiiC2pp_Published->SetMarkerColor(4);
 586                 RadiiC2pp_Published->SetLineColor(4);
 587               }
 588
 589               if(par==0) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("N");
 590               if(par==1) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#lambda_{e} or #lambda_{e,3}");
 591               if(par==2) {
 592                 if(FitType==0) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R}_{inv,2} or #font[12]{R}_{inv,3} (fm)");
 593                 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)");
 594                 else Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#font[12]{R}^{Exp}_{inv,2} or #font[12]{R}^{Exp}_{inv,3} (fm)");
 595               }
 596               if(par==3) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#kappa_{3}");
 597               if(par==4) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#kappa_{4}");
 598               if(par==5) Parameters_c3[ct][ft][KT3][par]->GetYaxis()->SetTitle("#lambda^{#font[12]{G}}_{e} or #lambda^{#font[12]{G}}_{e,3}");
 599               if(NchOneThirdAxis) Parameters_c3[ct][ft][KT3][par]->GetXaxis()->SetTitle("N_{ch}^{1/3}");
 600               else Parameters_c3[ct][ft][KT3][par]->GetXaxis()->SetTitle("#LT#font[12]{N}_{ch}#GT");
 601             }// par
 602           }// KT3
 603         }// ft
 604       }// dt==0
 605
 606     }// dt
 607   }// ct
 608
 609   cout<<"Done Getting Histograms"<<endl;
 610
 611   TF1 *Unity = new TF1("Unity","1",0,100);
 612   Unity->SetLineStyle(2);
 613   Unity->SetLineColor(1);
 614
 615
 616   /////////////////////////////////////
 617   // Bjoern's points
 618   double Bjoern_xaxis_pp[4] = {4.1, 9.2, 16, 26};// Nch
 619   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
 620   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
 621   //
 622   double Bjoern_xaxis_pPb[6] = {5, 11, 20, 32, 45, 69};// Nch
 623   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}};
 624   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}};
 625   //
 626   double Bjoern_xaxis_PbPb[10] = {22, 36, 50, 77, 98, 137, 172, 326, 498, 760};// Nch
 627   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}};
 628   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}};
 629   //
 630   int mchoice = 1;// 1 or 2 indicating Bjoern's cutoff
 631   double SF_Bjoern = 1.15;// 1.15 (m=0.1)
 632   //
 633
 634
 635   for(int index=0; index<4; index++){// pp
 636     for(int m=0; m<2; m++) {Bjoern_Ri_pp[m][index] *= SF_Bjoern; Bjoern_Rhydro_pp[m][index] *= SF_Bjoern;}
 637     double xpoint = Bjoern_xaxis_pp[index];
 638     if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_pp[index], 1/3.);
 639     int bin = Parameters_Bjoern[0][2]->GetXaxis()->FindBin(xpoint);
 640     Parameters_Bjoern[0][2]->SetBinContent(bin, Bjoern_Ri_pp[mchoice-1][index]);
 641     Parameters_Bjoern[0][2]->SetBinError(bin, 0.001);
 642     Parameters_Bjoern[1][2]->SetBinContent(bin, Bjoern_Rhydro_pp[mchoice-1][index]);
 643     Parameters_Bjoern[1][2]->SetBinError(bin, 0.001);
 644   }
 645   for(int index=0; index<6; index++){// pPb
 646     for(int m=0; m<2; m++) {Bjoern_Ri_pPb[m][index] *= SF_Bjoern; Bjoern_Rhydro_pPb[m][index] *= SF_Bjoern;}
 647     double xpoint = Bjoern_xaxis_pPb[index];
 648     if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_pPb[index], 1/3.);
 649     int bin = Parameters_Bjoern[0][1]->GetXaxis()->FindBin(xpoint);
 650     Parameters_Bjoern[0][1]->SetBinContent(bin, Bjoern_Ri_pPb[mchoice-1][index]);
 651     Parameters_Bjoern[0][1]->SetBinError(bin, 0.001);
 652     Parameters_Bjoern[1][1]->SetBinContent(bin, Bjoern_Rhydro_pPb[mchoice-1][index]);
 653     Parameters_Bjoern[1][1]->SetBinError(bin, 0.001);
 654   }
 655   for(int index=0; index<10; index++){// PbPb
 656     for(int m=0; m<2; m++) {Bjoern_Ri_PbPb[m][index] *= SF_Bjoern; Bjoern_Rhydro_PbPb[m][index] *= SF_Bjoern;}
 657     double xpoint = Bjoern_xaxis_PbPb[index];
 658     if(NchOneThirdAxis) xpoint = pow(Bjoern_xaxis_PbPb[index], 1/3.);
 659     int bin = Parameters_Bjoern[0][0]->GetXaxis()->FindBin(xpoint);
 660     Parameters_Bjoern[0][0]->SetBinContent(bin, Bjoern_Ri_PbPb[mchoice-1][index]);
 661     Parameters_Bjoern[0][0]->SetBinError(bin, 0.001);
 662     Parameters_Bjoern[1][0]->SetBinContent(bin, Bjoern_Rhydro_PbPb[mchoice-1][index]);
 663     Parameters_Bjoern[1][0]->SetBinError(bin, 0.001);
 664   }
 665   Parameters_Bjoern[0][0]->SetMarkerColor(1); Parameters_Bjoern[0][0]->SetMarkerStyle(20);
 666   Parameters_Bjoern[0][1]->SetMarkerColor(2); Parameters_Bjoern[0][1]->SetMarkerStyle(21);
 667   Parameters_Bjoern[0][2]->SetMarkerColor(4); Parameters_Bjoern[0][2]->SetMarkerStyle(34);
 668   //
 669   Parameters_Bjoern[1][0]->SetMarkerColor(1); Parameters_Bjoern[1][0]->SetMarkerStyle(20);
 670   Parameters_Bjoern[1][1]->SetMarkerColor(2); Parameters_Bjoern[1][1]->SetMarkerStyle(21);
 671   Parameters_Bjoern[1][2]->SetMarkerColor(4); Parameters_Bjoern[1][2]->SetMarkerStyle(34);
 672
 673   ////////////////////////////////////////////////////
 674   ////////////////////////////////////////////////////
 675   // Progaganda Plot
 676   /*
 677   TCanvas *can2 = new TCanvas("can2", "can2",10,0,600,600);// 11,53,700,500
 678   can2->SetHighLightColor(2);
 679   gStyle->SetOptFit(0111);
 680   can2->SetFillColor(0);//10
 681   can2->SetBorderMode(0);
 682   can2->SetBorderSize(2);
 683   can2->SetFrameFillColor(0);
 684   can2->SetFrameBorderMode(0);
 685   can2->SetFrameBorderMode(0);
 686   can2->cd();
 687   TPad *pad2 = new TPad("pad2","pad2",0.0,0.0,1.,1.);
 688   gPad->SetGridx(0);
 689   gPad->SetGridy(0);
 690   gPad->SetTickx(1);
 691   gPad->SetTicky(1);
 692   pad2->SetTopMargin(0.02);//0.05
 693   pad2->SetRightMargin(0.01);//3e-2
 694   pad2->SetBottomMargin(0.07);//0.12
 695   pad2->Draw();
 696   pad2->cd();
 697   TLegend *legend3 = new TLegend(.35,.75, .94,.95,NULL,"brNDC");//.45 or .4 for x1
 698   legend3->SetBorderSize(0);
 699   legend3->SetFillColor(0);
 700   legend3->SetTextFont(TextFont);
 701   legend3->SetTextSize(SizeLegend);
 702
 703   gPad->SetRightMargin(0.03); gPad->SetLeftMargin(0.12);
 704   gPad->SetBottomMargin(0.1); gPad->SetTopMargin(0.02);
 705   gPad->SetTickx(1);
 706   gPad->SetTicky(1);
 707   int Mb_pp=18;
 708   int Mb_pPb=13;
 709   int Mb_PbPb=0;
 710   c3[2][0][0][0][Mb_pp]->GetXaxis()->SetLabelSize(SizeLabel); c3[2][0][0][0][Mb_pp]->GetYaxis()->SetLabelSize(SizeLabel);
 711   c3[2][0][0][0][Mb_pp]->GetXaxis()->SetNdivisions(808);
 712   //
 713   c3[2][0][0][0][Mb_pp]->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}");
 714   c3[2][0][0][0][Mb_pp]->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})");
 715   c3[2][0][0][0][Mb_pp]->GetYaxis()->SetTitleOffset(1.1);
 716   c3[2][0][0][0][Mb_pp]->GetYaxis()->SetTitleSize(0.05);
 717   c3[2][0][0][0][Mb_pp]->SetMinimum(0.9);
 718   c3[2][0][0][0][Mb_pp]->Draw();
 719   //legend3->AddEntry(c3[2][0][0][Mb_pp],"N_{ch} = 9#pm0.2, pp #sqrt{s}=7 TeV","p");
 720   //legend3->AddEntry(c3[1][0][0][Mb_pPb],"N_{ch} = 59#pm2, p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV","p");
 721   //legend3->AddEntry(c3[0][0][0][Mb_PbPb],"N_{ch} = 1922#pm135, Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV","p");
 722   legend3->AddEntry(c3[2][0][0][0][Mb_pp],"Low N_{ch}, pp #sqrt{#font[12]{s}}=7 TeV","p");
 723   legend3->AddEntry(c3[1][0][0][0][Mb_pPb],"Mid N_{ch}, p-Pb #sqrt{#font[12]{#font[12]{s}_{NN}}}=5.02 TeV","p");
 724   legend3->AddEntry(c3[0][0][0][0][Mb_PbPb],"High N_{ch}, Pb-Pb #sqrt{#font[12]{#font[12]{s}_{NN}}}=2.76 TeV","p");
 725   //
 726   //TH1D *FillerHist = new TH1D("FillerHist","",50,0,0.5);
 727   for(int i=1; i<=50; i++) {
 728     if(i<13) {
 729       ExRangeTerm1->SetBinContent(i,10.0);
 730       //FillerHist->SetBinContent(i,10.0); FillerHist->SetBinError(i,.0001);
 731     }else {
 732       //FillerHist->SetBinContent(i,1.0); FillerHist->SetBinError(i,.0001);
 733       c3[0][0][0][0][Mb_PbPb]->SetBinContent(i,10.0);
 734       c3_Sys[0][0][0][0][Mb_PbPb]->SetBinContent(i,10.0);
 735     }
 736   }
 737   //FillerHist->SetMarkerStyle(24);
 738   //
 739   c3_Sys[2][0][0][0][Mb_pp]->Draw("E2 same");
 740   c3_Sys[1][0][0][0][Mb_pPb]->Draw("E2 same");
 741   c3_Sys[0][0][0][0][Mb_PbPb]->Draw("E2 same");
 742   c3[2][0][0][0][Mb_pp]->Draw("same");
 743   c3[1][0][0][0][Mb_pPb]->Draw("same");
 744   c3[0][0][0][0][Mb_PbPb]->Draw("same");
 745   //FillerHist->Draw("same");
 746   ExRangeTerm1->Draw("same");
 747
 748   TLatex *Specif_Kt3_1 = new TLatex(0.5,0.59,"0.16<#font[12]{K}_{T,3}<1.0 GeV/#font[12]{c}");
 749   Specif_Kt3_1->SetNDC();
 750   Specif_Kt3_1->SetTextFont(TextFont);
 751   Specif_Kt3_1->SetTextSize(SizeSpecif);
 752   Specif_Kt3_1->Draw("same");
 753   TLatex *Specif_Kinematics_1 = new TLatex(0.5,0.52,"|#eta|<0.8, 0.16<#font[12]{p}_{T}<1.0 GeV/#font[12]{c}");
 754   Specif_Kinematics_1->SetNDC();
 755   Specif_Kinematics_1->SetTextFont(TextFont);
 756   Specif_Kinematics_1->SetTextSize(SizeSpecif);
 757   Specif_Kinematics_1->Draw("same");
 758   TLatex *Species_1 = new TLatex(0.5,0.45,"#pi^{+}#pi^{+}#pi^{+} & #pi^{-}#pi^{-}#pi^{-} combined");
 759   Species_1->SetNDC();
 760   Species_1->SetTextFont(TextFont);
 761   Species_1->SetTextSize(SizeSpecif);
 762   Species_1->Draw("same");
 763
 764   legend3->Draw("same");
 765   DrawALICELogo(kTRUE, .75, .25, .9, .4);
 766
 767   if(SaveFiles) can2->SaveAs("ThreePionCorrelation_Evolution.eps");
 768   */
 769
 770
 771
 772   ////////////////////////////////////////////////////
 773   ////////////////////////////////////////////////////
 774   // K0s removal plot
 775   TCanvas *can1 = new TCanvas("can1", "can1",10,700,600,600);// 11,53,700,500
 776   can1->SetHighLightColor(2);
 777   gStyle->SetOptFit(0111);
 778   can1->SetFillColor(0);//10
 779   can1->SetBorderMode(0);
 780   can1->SetBorderSize(2);
 781   can1->SetFrameFillColor(0);
 782   can1->SetFrameBorderMode(0);
 783   can1->SetFrameBorderMode(0);
 784   can1->cd();
 785   TPad *pad1 = new TPad("pad1","pad1",0.,0.,1.,1.);
 786   gPad->SetTickx();
 787   gPad->SetTicky();
 788   pad1->SetTopMargin(0.0);//0.05
 789   pad1->SetRightMargin(0.0);//3e-2
 790   pad1->SetBottomMargin(0.0);//0.12
 791   pad1->SetLeftMargin(0.0);
 792   pad1->Draw();
 793   pad1->cd();
 794   TLegend *legend1 = new TLegend(.2,.65, .55,.85,NULL,"brNDC");//.45 or .4 for x1
 795   legend1->SetBorderSize(0);
 796   legend1->SetFillColor(0);
 797   legend1->SetTextFont(TextFont);
 798   legend1->SetTextSize(SizeLegend);
 799
 800   gPad->SetRightMargin(0.01); gPad->SetLeftMargin(0.14);
 801   gPad->SetBottomMargin(0.14); gPad->SetTopMargin(0.02);
 802   gPad->SetTickx(); gPad->SetTicky();
 803   // pp, M18
 804   double points_K0s_C3[50]={1.44542, 1.33366, 1.25518, 1.21051, 1.20484, 1.18742, 1.18246, 1.17024, 1.16383, 1.1555, 1.15536, 1.15272, 1.14571, 1.14822, 1.14708, 1.14066, 1.14133, 1.13866, 1.13755, 1.13308, 1.1386, 1.13821, 1.13429, 1.13707, 1.13294, 1.13545, 1.13929, 1.13374, 1.13546, 1.13205, 1.12846, 1.1281, 1.12947, 1.13162, 1.13822, 1.13626, 1.13912, 1.14379, 1.15594, 1.19387, 1.29685, 1.21203, 1.14452, 1.14288, 1.14161, 1.1408, 1.13796, 1.13661, 1.13565, 1.13312};
 805   double points_K0s_C3_e[50]={0.0270571, 0.00940436, 0.0058043, 0.00430829, 0.00355541, 0.00306969, 0.00275566, 0.00253068, 0.00236966, 0.00224498, 0.00215682, 0.00208189, 0.00201555, 0.00196896, 0.00192381, 0.0018782, 0.00184374, 0.0018119, 0.00178103, 0.00175109, 0.00173313, 0.00171119, 0.00168873, 0.00167507, 0.00165525, 0.00164399, 0.00163488, 0.00161686, 0.00160668, 0.00159163, 0.00157479, 0.00155996, 0.00154251, 0.00152405, 0.00150469, 0.00147687, 0.00144939, 0.00142223, 0.0013988, 0.00139214, 0.00142684, 0.0013476, 0.00128368, 0.00126468, 0.0012497, 0.00123827, 0.001228, 0.00122253, 0.00121921, 0.00121782};
 806   double points_K0s_c3[50]={0.965148, 1.05077, 0.983133, 0.972848, 0.993922, 0.994064, 1.00618, 0.993434, 1.00322, 0.999747, 0.982221, 1.00029, 0.992184, 0.993566, 1.00007, 0.995263, 0.998559, 0.987037, 0.987919, 0.991035, 0.991488, 0.99592, 0.99677, 0.990709, 0.990352, 0.994706, 0.99606, 0.998525, 0.994121, 0.986511, 0.991009, 0.990309, 1.00161, 0.997478, 1.00117, 0.989375, 0.996592, 0.990832, 0.998989, 0.995124, 0.996301, 1.00189, 0.995304, 0.992511, 0.994084, 0.994777, 1.00074, 0.996959, 0.998713, 0.996205};
 807   double points_K0s_c3_e[50]={0.0630244, 0.0223583, 0.0140589, 0.0105446, 0.00871274, 0.00755392, 0.00678715, 0.00625387, 0.00586333, 0.00556614, 0.00535136, 0.00516553, 0.00501148, 0.00489185, 0.00478054, 0.00467542, 0.00458851, 0.00451377, 0.00443821, 0.00436787, 0.00431762, 0.00426282, 0.0042107, 0.00417461, 0.00412901, 0.00409783, 0.00407137, 0.0040317, 0.00400527, 0.00397212, 0.00393349, 0.00389671, 0.00385012, 0.00380258, 0.00374821, 0.0036816, 0.0036096, 0.00353852, 0.00346903, 0.00342075, 0.00342477, 0.00329675, 0.0031932, 0.00314748, 0.00311129, 0.00308337, 0.00305968, 0.00304743, 0.00303989, 0.00303859};
 808   TH1D *K0s_C3=new TH1D("K0s_C3","",50,0,0.5);
 809   TH1D *K0s_c3=new TH1D("K0s_c3","",50,0,0.5);
 810   for(int i=0; i<50; i++){
 811     K0s_C3->SetBinContent(i+1, points_K0s_C3[i]);
 812     K0s_C3->SetBinError(i+1, points_K0s_C3_e[i]);
 813     K0s_c3->SetBinContent(i+1, points_K0s_c3[i]);
 814     K0s_c3->SetBinError(i+1, points_K0s_c3_e[i]);
 815   }
 816   K0s_c3->SetMarkerStyle(20);
 817   K0s_c3->SetMarkerColor(4);
 818   K0s_C3->SetMarkerStyle(24);
 819   K0s_C3->SetMarkerColor(4);
 820   K0s_c3->SetMarkerSize(K0s_C3->GetMarkerSize() * 1.12);
 821   K0s_C3->SetMinimum(0.92);
 822   K0s_C3->GetXaxis()->SetTitle("#font[12]{q_{31}^{#pm#mp}} (GeV/#font[12]{c})");
 823   K0s_C3->GetYaxis()->SetTitle("#font[12]{C_{3}} or #font[12]{#bf{c}_{3}}");
 824   K0s_C3->GetXaxis()->SetTitleOffset(1.05); K0s_C3->GetYaxis()->SetTitleOffset(1.12);
 825   K0s_C3->GetXaxis()->SetTitleSize(SizeTitle); K0s_C3->GetYaxis()->SetTitleSize(SizeTitle);
 826   K0s_C3->GetXaxis()->SetLabelSize(SizeTitle); K0s_C3->GetYaxis()->SetLabelSize(SizeTitle);
 827   K0s_C3->GetXaxis()->SetNdivisions(404);
 828   K0s_C3->GetYaxis()->SetNdivisions(404);
 829   //K0s_C3->GetXaxis()->SetRangeUser(-0.001,0.5);
 830
 831   K0s_C3->Draw();
 832   // Sys errors
 833   TH1D *Sys_K0s_C3=(TH1D*)K0s_C3->Clone();
 834   TH1D *Sys_K0s_c3=(TH1D*)K0s_c3->Clone();
 835   Sys_K0s_C3->SetMarkerSize(0); Sys_K0s_c3->SetMarkerSize(0);
 836   Sys_K0s_C3->SetFillColor(kBlue-10); Sys_K0s_c3->SetFillColor(kBlue-10);
 837   Sys_K0s_C3->SetFillStyle(1000); Sys_K0s_c3->SetFillStyle(1000);
 838   cout.precision(3);
 839   for(int i=0; i<Sys_K0s_C3->GetNbinsX(); i++) {
 840     Sys_K0s_C3->SetBinError(i+1, 0.01 * Sys_K0s_C3->GetBinContent(i+1));
 841     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)));
 842     cout<<K0s_C3->GetXaxis()->GetBinLowEdge(i+1)<<" TO "<<K0s_C3->GetXaxis()->GetBinUpEdge(i+1)<<"; "<<K0s_C3->GetBinContent(i+1)<<" +- "<<K0s_C3->GetBinError(i+1)<<"  (DSYS="<<Sys_K0s_C3->GetBinError(i+1)<<"); "<<K0s_c3->GetBinContent(i+1)<<" +- "<<K0s_c3->GetBinError(i+1)<<"  (DSYS="<<Sys_K0s_c3->GetBinError(i+1)<<"); "<<endl;
 843     //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;
 844   }
 845
 846   Sys_K0s_C3->Draw("E2 same");
 847   Sys_K0s_c3->Draw("E2 same");
 848   K0s_C3->Draw("same");
 849   K0s_c3->Draw("same");
 850
 851   legend1->AddEntry(K0s_C3,"#font[12]{C_{3}^{#pm#pm#mp}} projection","p");
 852   legend1->AddEntry(K0s_c3,"#font[12]{#bf{c}_{3}^{#pm#pm#mp}} projection","p");
 853   legend1->Draw("same");
 854
 855   TLatex *Specif_qRange = new TLatex(0.25,0.6,"0.2 < #font[12]{q_{12}^{#pm#pm},q_{23}^{#pm#mp}} < 0.5 GeV/#font[12]{c}");
 856   Specif_qRange->SetNDC();
 857   Specif_qRange->SetTextFont(42);
 858   Specif_qRange->SetTextSize(SizeSpecif);
 859   Specif_qRange->Draw("same");
 860   TLatex *Specif_System = new TLatex(0.25,0.9,"pp #sqrt{#font[12]{s}}=7 TeV, #LT#font[12]{N}_{ch}#GT = 8.6 #pm 0.4");
 861   Specif_System->SetNDC();
 862   Specif_System->SetTextFont(42);
 863   Specif_System->SetTextSize(SizeSpecif);
 864   Specif_System->Draw("same");
 865
 866
 867
 868
 869   ////////////////////////////////////////////////////
 870   ////////////////////////////////////////////////////
 871   // Radii and lambda plot
 872   float SF_2panel=1.2;
 873   float RadiiC2ppPubPoints[2][8]={{0.88,1.09,1.23,1.28,1.34,1.37,1.42,1.48},{0.86,1.06,1.16,1.23,1.23,1.28,1.32,1.38}};
 874   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}};
 875   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
 876
 877   int yExtent = 900;
 878   if(RadiusOnly) yExtent = 600;
 879   TCanvas *can3 = new TCanvas("can3", "can3",1000,0,600,yExtent);// 1000,0,600,900
 880   can3->SetHighLightColor(2);
 881   gStyle->SetOptFit(0111);
 882   can3->SetFillColor(0);//10
 883   can3->SetBorderMode(0);
 884   can3->SetBorderSize(2);
 885   can3->SetFrameFillColor(0);
 886   can3->SetFrameBorderMode(0);
 887   can3->SetFrameBorderMode(0);
 888   can3->cd();
 889   float PadLeftMargin=0.0, PadBottomMargin=0.;
 890   TPad *pad3 = new TPad("pad3","pad3",PadLeftMargin,PadBottomMargin,1.,1.);
 891   gPad->SetTickx();
 892   gPad->SetTicky();
 893   pad3->SetTopMargin(0.0);//0.05
 894   pad3->SetRightMargin(0.0);//3e-2
 895   pad3->SetBottomMargin(0.0);//0.12
 896   pad3->SetLeftMargin(0.0);
 897   if(!RadiusOnly) pad3->Divide(1,2,0,0);
 898   pad3->Draw();
 899   pad3->cd();
 900   TLegend *legend4 = new TLegend(.15,.65, .55,.95,NULL,"brNDC");//.45 or .4 for x1
 901   legend4->SetBorderSize(0);
 902   legend4->SetFillColor(0);
 903   legend4->SetTextFont(TextFont);
 904   legend4->SetTextSize(SizeLegend);
 905   //
 906   pad3->cd(1);
 907   gPad->SetLeftMargin(0.14);
 908   gPad->SetRightMargin(0.01);
 909   gPad->SetTopMargin(0.01);
 910   gPad->SetBottomMargin(0.001);// 0.16
 911   if(RadiusOnly){
 912     gPad->SetLeftMargin(0.16);
 913     gPad->SetBottomMargin(0.16);
 914   }
 915
 916   gPad->SetTickx(); gPad->SetTicky();
 917   if(!NchOneThirdAxis)gPad->SetLogx();
 918
 919   TH1D *RadiiPbPb=(TH1D*)Parameters_c3[0][FitType][KT3Bin][2]->Clone();
 920   TH1D *RadiipPb=(TH1D*)Parameters_c3[1][FitType][KT3Bin][2]->Clone();
 921   TH1D *Radiipp=(TH1D*)Parameters_c3[2][FitType][KT3Bin][2]->Clone();
 922   RadiiPbPb->GetXaxis()->SetLabelFont(TextFont); RadiiPbPb->GetYaxis()->SetLabelFont(TextFont);
 923   RadiiPbPb->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); RadiiPbPb->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel);
 924   RadiiPbPb->GetXaxis()->SetNdivisions(808);
 925   RadiiPbPb->GetXaxis()->SetTitleOffset(0.95);//1.3
 926   RadiiPbPb->GetYaxis()->SetTitleOffset(100);//1.1
 927   RadiiPbPb->GetXaxis()->SetTitleFont(TextFont); RadiiPbPb->GetXaxis()->SetTitleSize(0);// SizeTitle*SF_2panel
 928   RadiiPbPb->GetYaxis()->SetTitleFont(TextFont); RadiiPbPb->GetYaxis()->SetTitleSize(0);// SizeTitle*SF_2panel
 929   RadiiPbPb->SetMinimum(0.01); RadiiPbPb->SetMaximum(11.9);// 0 and 11.9
 930   //gStyle->SetErrorX(0);
 931   if(NchOneThirdAxis) RadiiPbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000
 932   else RadiiPbPb->GetXaxis()->SetRangeUser(3,3000);// 3,3000
 933   if(RadiusOnly){
 934     RadiiPbPb->GetXaxis()->SetTitleSize(SizeTitle);
 935     RadiiPbPb->GetXaxis()->SetTitleOffset(1.25);
 936     RadiiPbPb->GetYaxis()->SetTitleSize(SizeTitle);
 937     RadiiPbPb->GetYaxis()->SetTitleOffset(1.2);
 938     if(FitType==0) RadiiPbPb->GetYaxis()->SetTitle("#font[12]{R}^{#font[12]{G}}_{inv} or #font[12]{R}^{#font[12]{G}}_{inv,3} (fm)");
 939     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)");
 940     else RadiiPbPb->GetYaxis()->SetTitle("#font[12]{R}^{#font[12]{Exp}_{inv} or #font[12]{R}^{#font[12]{Exp}_{inv,3} (fm)");
 941   }
 942
 943   RadiiPbPb->Draw();
 944
 945   //
 946   double xAxisC2_e[20]={0};
 947   double xAxis_e[20]={0};
 948   double yAxisPbPb[20]={0};
 949   double yAxisPbPb_eL[20]={0};
 950   double yAxisPbPb_eH[20]={0};
 951   double yAxispPb[20]={0};
 952   double yAxispPb_eL[20]={0};
 953   double yAxispPb_eH[20]={0};
 954   double yAxispp[20]={0};
 955   double yAxispp_eL[20]={0};
 956   double yAxispp_eH[20]={0};
 957
 958   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)
 959   // Narrow Fit Range for C2 fits
 960   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
 961   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
 962   // Wide Fit Range for C2 fits
 963   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
 964   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
 965   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
 966
 967   for(int cb=0; cb<20; cb++){// 3-particle
 968     int binPbPb = RadiiPbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
 969     int binpPb = RadiipPb->GetXaxis()->FindBin(meanNchpPb[cb]);
 970     int binpp = Radiipp->GetXaxis()->FindBin(meanNchpp[cb]);
 971     double RsysPbPb = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][1],2) + pow(1,2)) *  RadiiPbPb->GetBinContent(binPbPb);// fit Variation + MRC
 972     double RsyspPb = 0.01*sqrt(pow(SysPercent_pPb_NFR[FitType][1],2) + pow(1,2)) *  RadiipPb->GetBinContent(binpPb);// fit Variation + MRC
 973     double Rsyspp = 0.01*sqrt(pow(SysPercent_pp_NFR[FitType][1],2) + pow(1,2)) *  Radiipp->GetBinContent(binpp);// fit Variation + MRC
 974     //
 975     if(cb==19 && FitType==2) Rsyspp =  0.01*sqrt(pow(16,2) + pow(1,2)) *  Radiipp->GetBinContent(binpp);// larger fit var in this bin
 976
 977     if(RadiiPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100
 978     else {yAxisPbPb[cb]=RadiiPbPb->GetBinContent(binPbPb); yAxisPbPb_eL[cb]=RsysPbPb; yAxisPbPb_eH[cb]=RsysPbPb;}
 979     //
 980     if(RadiipPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;}
 981     else {yAxispPb[cb]=RadiipPb->GetBinContent(binpPb); yAxispPb_eL[cb]=RsyspPb; yAxispPb_eH[cb]=RsyspPb;}
 982     //
 983     if(Radiipp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=100; yAxispp_eH[cb]=100;}
 984     else {yAxispp[cb]=Radiipp->GetBinContent(binpp); yAxispp_eL[cb]=Rsyspp; yAxispp_eH[cb]=Rsyspp;}
 985     //
 986
 987     if(NchOneThirdAxis) {
 988       if(cb<13) xAxis_e[cb]=fabs(meanNchPbPb[cb] - pow(0.95,1/3.)*meanNchPbPb[cb]);
 989       else xAxis_e[cb]=fabs(meanNchpPb[cb] - pow(0.95,1/3.)*meanNchpPb[cb]);
 990     }else {
 991       if(cb<13) xAxis_e[cb]=fabs(meanNchPbPb[cb] - (0.95)*meanNchPbPb[cb]);
 992       else xAxis_e[cb]=fabs(meanNchpPb[cb] - (0.95)*meanNchpPb[cb]);
 993     }
 994   }
 995
 996
 997   TGraphErrors *grRadiiSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_eL);
 998   TGraphAsymmErrors *grRadiiSys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH);
 999   TGraphAsymmErrors *grRadiiSys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH);
1000   grRadiiSys_pp->SetMarkerSize(0); grRadiiSys_pp->SetFillStyle(1000); grRadiiSys_pp->SetFillColor(kBlue-10);
1001   grRadiiSys_pPb->SetMarkerSize(0); grRadiiSys_pPb->SetFillStyle(1000); grRadiiSys_pPb->SetFillColor(kRed-10);
1002   grRadiiSys_PbPb->SetMarkerSize(0); grRadiiSys_PbPb->SetFillStyle(1000); grRadiiSys_PbPb->SetFillColor(kGray);
1003   grRadiiSys_pp->SetMarkerColor(kBlue-10); grRadiiSys_pp->SetMarkerColor(kRed-10); grRadiiSys_pp->SetMarkerColor(kGray);
1004   // C2
1005   TH1D *RadiiC2PbPb=(TH1D*)Parameters_C2[0][FitType][KT3Bin][2]->Clone();
1006   TH1D *RadiiC2pPb=(TH1D*)Parameters_C2[1][FitType][KT3Bin][2]->Clone();
1007   TH1D *RadiiC2pp=(TH1D*)Parameters_C2[2][FitType][KT3Bin][2]->Clone();
1008   RadiiC2pp_Published->SetMarkerStyle(30);
1009   //if(FitType==0) RadiiC2pp->SetMarkerStyle(30);// for legend marker
1010
1011   for(int mbin=0; mbin<8; mbin++){
1012     int bin = RadiiC2pp_Published->GetXaxis()->FindBin(MeanPubNch[mbin]);
1013     RadiiC2pp_Published->SetBinContent(bin, RadiiC2ppPubPoints[KT3Bin][mbin]);
1014     RadiiC2pp_Published->SetBinError(bin, RadiiC2ppPubPoints_e[KT3Bin][mbin]);
1015   }
1016
1017   for(int cb=0; cb<20; cb++){// 2-particle
1018     int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
1019     int binpPb = RadiiC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]);
1020     int binpp = RadiiC2pp->GetXaxis()->FindBin(meanNchpp[cb]);
1021     double RsysPbPb_L = 0.01*sqrt(pow(SysPercent_PbPb_WFR[FitType][0],2) + pow(1,2)) *  RadiiC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC
1022     double RsysPbPb_H = 0.01*sqrt(pow(SysPercent_PbPb_NFR[FitType][0],2) + pow(1,2)) *  RadiiC2PbPb->GetBinContent(binPbPb);// fit Variation + MRC
1023     double RsyspPb_L = 0.01*sqrt(pow(SysPercent_pPb_WFR[FitType][0],2) + pow(1,2)) *  RadiiC2pPb->GetBinContent(binpPb);// fit Variation + MRC
1024     double RsyspPb_H = 0.01*sqrt(pow(SysPercent_pPb_NFR[FitType][0],2) + pow(1,2)) *  RadiiC2pPb->GetBinContent(binpPb);// fit Variation + MRC
1025     double Rsyspp_L = 0.01*sqrt(pow(SysPercent_pp_WFR[FitType][0],2) + pow(1,2)) *  RadiiC2pp->GetBinContent(binpp);// fit Variation + MRC
1026     double Rsyspp_H = 0.01*sqrt(pow(SysPercent_pp_NFR[FitType][0],2) + pow(1,2)) *  RadiiC2pp->GetBinContent(binpp);// fit Variation + MRC
1027     //
1028     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
1029     //
1030
1031     if(RadiiC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=1000; yAxisPbPb_eH[cb]=1000;}// errors were 1000
1032     else {
1033       yAxisPbPb[cb]=RadiiC2PbPb->GetBinContent(binPbPb);
1034       if(cb>=13) yAxisPbPb_eL[cb]=RsysPbPb_L;
1035       else yAxisPbPb_eL[cb]=RsysPbPb_H;
1036       yAxisPbPb_eH[cb]=RsysPbPb_H;
1037     }
1038     if(RadiiC2pPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=1000; yAxispPb_eH[cb]=1000;}
1039     else {yAxispPb[cb]=RadiiC2pPb->GetBinContent(binpPb); yAxispPb_eL[cb]=RsyspPb_L; yAxispPb_eH[cb]=RsyspPb_H;}
1040     if(RadiiC2pp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=1000; yAxispp_eH[cb]=1000;}
1041     else {yAxispp[cb]=RadiiC2pp->GetBinContent(binpp); yAxispp_eL[cb]=Rsyspp_L; yAxispp_eH[cb]=Rsyspp_H;}
1042     //
1043   }
1044
1045   TGraphAsymmErrors *grRadiiC2Sys_PbPb=new TGraphAsymmErrors(20,meanNchPbPb,yAxisPbPb, xAxisC2_e,xAxisC2_e, yAxisPbPb_eL,yAxisPbPb_eH);
1046   TGraphAsymmErrors *grRadiiC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxisC2_e,xAxisC2_e, yAxispPb_eL,yAxispPb_eH);
1047   TGraphAsymmErrors *grRadiiC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxisC2_e,xAxisC2_e, yAxispp_eL,yAxispp_eH);
1048   grRadiiC2Sys_pp->SetMarkerSize(0);
1049   grRadiiC2Sys_pPb->SetMarkerSize(0);
1050   grRadiiC2Sys_PbPb->SetMarkerSize(0);
1051   grRadiiC2Sys_pp->SetLineColor(4); grRadiiC2Sys_pPb->SetLineColor(2); grRadiiC2Sys_PbPb->SetLineColor(1);
1052   grRadiiC2Sys_pp->SetLineWidth(2.*grRadiiC2Sys_pp->GetLineWidth());
1053   grRadiiC2Sys_pPb->SetLineWidth(2.*grRadiiC2Sys_pPb->GetLineWidth());
1054   grRadiiC2Sys_PbPb->SetLineWidth(1.*grRadiiC2Sys_PbPb->GetLineWidth()); grRadiiC2Sys_PbPb->SetLineStyle(2);
1055   //
1056   grRadiiC2Sys_pPb->Draw("|| p");
1057   grRadiiC2Sys_pp->Draw("|| p");
1058
1059   grRadiiSys_pp->Draw("E2 p");
1060   grRadiiSys_pPb->Draw("E2 p");
1061   grRadiiSys_PbPb->Draw("E2 p");
1062   RadiiPbPb->Draw("same");
1063   RadiipPb->Draw("same");
1064   Radiipp->Draw("same");
1065   grRadiiC2Sys_PbPb->Draw("E");// E2 or "|| p" to visualize pol2 fit below
1066   //if(FitType==0) RadiiC2pp_Published->Draw("same");
1067   //
1068
1069   RadiiC2PbPb->Draw("same");
1070   RadiiC2pPb->Draw("same");
1071   RadiiC2pp->Draw("same");
1072
1073
1074   legend4->AddEntry(Radiipp,"pp #sqrt{s}=7 TeV","p");
1075   legend4->AddEntry(RadiipPb,"p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV","p");
1076   legend4->AddEntry(RadiiPbPb,"Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV","p");
1077
1078   TF1 *ppLine = new TF1("ppLine","pol1",0,13);
1079   ppLine->SetLineColor(4);
1080   TF1 *pPbLine = new TF1("pPbLine","pol1",0,13);
1081   TF1 *PbPbLine = new TF1("PbPbLine","pol1",0,13);
1082   PbPbLine->SetLineColor(1);
1083   if(NchOneThirdAxis){
1084     Radiipp->Fit(ppLine,"IMEN","",1,4.);
1085     ppLine->Draw("same");
1086     RadiipPb->Fit(pPbLine,"IMEN","",2,4.5);
1087     pPbLine->Draw("same");
1088     RadiiC2PbPb->Fit(PbPbLine,"IMEN","",4,13);
1089     PbPbLine->Draw("same");
1090   }
1091
1092   TLatex *Specif_Kt3;
1093   TLatex *Specif_kt;
1094   if(KT3Bin==0) {
1095     Specif_Kt3 = new TLatex(0.17, 0.57,"0.16<#font[12]{K}_{T,3}<0.3 GeV/#font[12]{c}");
1096     Specif_kt = new TLatex(0.17, 0.47,"0.2<#font[12]{k}_{T}<0.3 GeV/#font[12]{c}");
1097     // KT3:  #LT#font[12]{k}_{T}#GT=0.24 GeV/#font[12]{c}
1098     // kT:  #LT#font[12]{k}_{T}#GT=0.25 GeV/#font[12]{c}
1099   }
1100   if(KT3Bin==1) {
1101     Specif_Kt3 = new TLatex(0.17, 0.57,"0.3<#font[12]{K}_{T,3}<1.0 GeV/#font[12]{c}");
1102     Specif_kt = new TLatex(0.17, 0.47,"0.3<#font[12]{k}_{T}<1.0 GeV/#font[12]{c}");
1103     // KT3:  #LT#font[12]{k}_{T}#GT=0.39 GeV/#font[12]{c}
1104     // kT:  #LT#font[12]{k}_{T}#GT=0.43 GeV/#font[12]{c}
1105   }
1106   //if(KT3Bin==0) {Specif_Kt3 = new TLatex(0.57, 0.83,"0.16<K_{T,3}<0.25 GeV/#font[12]{c}"); Specif_kt = new TLatex(0.57, 0.76,"0.2<k_{T}<0.3 GeV/#font[12]{c}");}
1107   //if(KT3Bin==1) {Specif_Kt3 = new TLatex(0.57, 0.83,"0.25<K_{T,3}<0.35 GeV/#font[12]{c}"); Specif_kt = new TLatex(0.57, 0.76,"0.3<k_{T}<0.4 GeV/#font[12]{c}");}
1108   //if(KT3Bin==2) {Specif_Kt3 = new TLatex(0.57, 0.83,"0.35<K_{T,3}<1.0 GeV/#font[12]{c}"); Specif_kt = new TLatex(0.57, 0.76,"0.4<k_{T}<0.5 GeV/#font[12]{c}");}
1109   Specif_Kt3->SetTextFont(TextFont); Specif_kt->SetTextFont(TextFont);
1110   Specif_Kt3->SetTextSize(SizeSpecif*SF_2panel); Specif_kt->SetTextSize(SizeSpecif*SF_2panel);
1111   Specif_Kt3->SetNDC(); Specif_kt->SetNDC();
1112   if(!RadiusOnly){
1113     Specif_Kt3->Draw("same");
1114     Specif_kt->Draw("same");
1115   }
1116   legend4->SetTextFont(TextFont);
1117   legend4->SetTextSize(SizeLegend*SF_2panel);
1118   if(RadiusOnly) legend4->SetTextSize(SizeLegend*SF_2panel*0.8);
1119   legend4->Draw("same");
1120
1121   TH1D *MarkerPbPb_3=(TH1D*)RadiiPbPb->Clone();
1122   TH1D *MarkerpPb_3=(TH1D*)RadiipPb->Clone();
1123   TH1D *Markerpp_3=(TH1D*)Radiipp->Clone();
1124   TH1D *MarkerPbPb_2=(TH1D*)RadiiC2PbPb->Clone();
1125   TH1D *MarkerpPb_2=(TH1D*)RadiiC2pPb->Clone();
1126   TH1D *Markerpp_2=(TH1D*)RadiiC2pp->Clone();
1127   for(int i=1; i<=MarkerPbPb_3->GetNbinsX(); i++){
1128     MarkerPbPb_3->SetBinContent(i,1000); MarkerpPb_3->SetBinContent(i,1000); Markerpp_3->SetBinContent(i,1000);
1129     MarkerPbPb_2->SetBinContent(i,1000); MarkerpPb_2->SetBinContent(i,1000); Markerpp_2->SetBinContent(i,1000);
1130     MarkerPbPb_3->SetBinError(i,0.001); MarkerpPb_3->SetBinError(i,0.001); Markerpp_3->SetBinError(i,0.001);
1131     MarkerPbPb_2->SetBinError(i,0.001); MarkerpPb_2->SetBinError(i,0.001); Markerpp_2->SetBinError(i,0.001);
1132   }
1133   if(!NchOneThirdAxis){
1134     MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(450), 1.25);// 1.25, 1.45 for RadiusOnly
1135     MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(600), 1.25);// 1.25, 1.45 for RadiusOnly
1136     Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(800), 1.25);// 1.25, 1.45 for RadiusOnly
1137     MarkerPbPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(450), 3.1);// 3.1
1138     MarkerpPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(600), 3.1);// 3.1
1139     Markerpp_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(800), 3.1);// 3.1
1140   }else{
1141     MarkerPbPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 1.25);// 1.25
1142     MarkerpPb_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 1.25);// 1.25
1143     Markerpp_3->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 1.25);// 1.25
1144     MarkerPbPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10), 3.1);// 3.1
1145     MarkerpPb_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(10.5), 3.1);// 3.1
1146     Markerpp_2->SetBinContent(MarkerPbPb_3->GetXaxis()->FindBin(11), 3.1);// 3.1
1147   }
1148
1149   MarkerPbPb_3->Draw("same"); MarkerpPb_3->Draw("same"); Markerpp_3->Draw("same");
1150   MarkerPbPb_2->Draw("same"); MarkerpPb_2->Draw("same"); Markerpp_2->Draw("same");
1151
1152   TLatex *TwoPionText;
1153   if(!RadiusOnly) TwoPionText = new TLatex(0.74,0.3,"Two-Pions");// 0.74,0.3
1154   else TwoPionText = new TLatex(0.67,0.31,"Two-Pions");//
1155   TLatex *ThreePionText;
1156   if(!RadiusOnly) ThreePionText = new TLatex(0.74,0.15,"Three-Pions");// 0.74,0.15
1157   else ThreePionText = new TLatex(0.65,0.2,"Three-Pions");//
1158   TwoPionText->SetNDC(); ThreePionText->SetNDC();
1159   TwoPionText->SetTextFont(TextFont); ThreePionText->SetTextFont(TextFont);
1160   TwoPionText->SetTextSize(SizeSpecif*SF_2panel); ThreePionText->SetTextSize(SizeSpecif*SF_2panel);
1161   TwoPionText->Draw("same");
1162   ThreePionText->Draw("same");
1163
1164   TLatex *Specif_kappas = new TLatex(0.42,0.05,"#kappa_{3}=0.1, #kappa_{4}=0.5");// 0.42,0.2
1165   //TLatex *Specif_kappas = new TLatex(0.42,0.05,"#kappa_{3}(N_{ch}), #kappa_{4}=0.5");// 0.42,0.2
1166   Specif_kappas->SetNDC();
1167   Specif_kappas->SetTextFont(TextFont);
1168   Specif_kappas->SetTextSize(SizeSpecif*SF_2panel);
1169   if(FitType==1 && !RadiusOnly) Specif_kappas->Draw("same");
1170
1171   if(RadiusOnly) return;
1172
1173   ///////////////////////////////////////////////////////////////////
1174   pad3->cd(2);
1175   gPad->SetLeftMargin(0.14);
1176   gPad->SetRightMargin(0.01);
1177   gPad->SetTopMargin(0.0);// 0.01
1178   gPad->SetBottomMargin(0.16);
1179   gPad->SetTickx(); gPad->SetTicky();
1180   if(!NchOneThirdAxis) gPad->SetLogx();
1181   //gPad->SetGridx(); gPad->SetGridy();
1182   TH1D *LambdaPbPb=(TH1D*)Parameters_c3[0][FitType][KT3Bin][5]->Clone();
1183   TH1D *LambdapPb=(TH1D*)Parameters_c3[1][FitType][KT3Bin][5]->Clone();
1184   TH1D *Lambdapp=(TH1D*)Parameters_c3[2][FitType][KT3Bin][5]->Clone();
1185
1186   LambdaPbPb->GetXaxis()->SetLabelFont(TextFont); LambdaPbPb->GetYaxis()->SetLabelFont(TextFont);
1187   LambdaPbPb->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); LambdaPbPb->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel);
1188   LambdaPbPb->GetXaxis()->SetNdivisions(808);
1189   LambdaPbPb->GetYaxis()->SetNdivisions(604);
1190   LambdaPbPb->GetXaxis()->SetTitleFont(TextFont); LambdaPbPb->GetXaxis()->SetTitleSize(SizeTitle*SF_2panel);
1191   LambdaPbPb->GetYaxis()->SetTitleFont(TextFont); LambdaPbPb->GetYaxis()->SetTitleSize(SizeTitle*SF_2panel);
1192   LambdaPbPb->SetMaximum(2.8);// 2.8
1193   if(FitType==2) LambdaPbPb->SetMaximum(5.8);
1194   LambdaPbPb->GetXaxis()->SetTitleOffset(0.95);
1195   LambdaPbPb->GetYaxis()->SetTitleOffset(100);//1.1
1196   if(NchOneThirdAxis) LambdaPbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000
1197   else LambdaPbPb->GetXaxis()->SetRangeUser(3,3000);// 3,3000
1198   LambdaPbPb->Draw();
1199
1200   TF1 *ChaoticLimit_C2 = new TF1("ChaoticLimit_C2","1.0",0,5000);
1201   TF1 *ChaoticLimit_c3 = new TF1("ChaoticLimit_c3","2.0",0,5000);
1202   ChaoticLimit_C2->SetLineColor(1); ChaoticLimit_c3->SetLineColor(1);
1203   ChaoticLimit_C2->SetLineStyle(7); ChaoticLimit_c3->SetLineStyle(6);
1204   ChaoticLimit_C2->Draw("same");
1205   ChaoticLimit_c3->Draw("same");
1206
1207   for(int cb=0; cb<20; cb++){// 3-particle
1208     int binPbPb = LambdaPbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
1209     int binpPb = LambdapPb->GetXaxis()->FindBin(meanNchpPb[cb]);
1210     int binpp = Lambdapp->GetXaxis()->FindBin(meanNchpp[cb]);
1211     double f_syst_PbPb = 0, f_syst_pPb=0, f_syst_pp=0;
1212     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
1213     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
1214     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
1215     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
1216     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
1217     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
1218     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
1219     if(LambdaPbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100
1220     else {yAxisPbPb[cb]=LambdaPbPb->GetBinContent(binPbPb); yAxisPbPb_eL[cb]=LambdasysPbPb; yAxisPbPb_eH[cb]=LambdasysPbPb;}
1221     //
1222     if(LambdapPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;}
1223     else {yAxispPb[cb]=LambdapPb->GetBinContent(binpPb); yAxispPb_eL[cb]=LambdasyspPb; yAxispPb_eH[cb]=LambdasyspPb;}
1224     //
1225     if(Lambdapp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=100; yAxispp_eH[cb]=100;}
1226     else {yAxispp[cb]=Lambdapp->GetBinContent(binpp); yAxispp_eL[cb]=Lambdasyspp; yAxispp_eH[cb]=Lambdasyspp;}
1227
1228
1229     if(NchOneThirdAxis) {
1230       if(cb<13) xAxis_e[cb]=fabs(meanNchPbPb[cb] - pow(0.95,1/3.)*meanNchPbPb[cb]);
1231       else xAxis_e[cb]=fabs(meanNchpPb[cb] - pow(0.95,1/3.)*meanNchpPb[cb]);
1232     }else {
1233       if(cb<13) xAxis_e[cb]=fabs(meanNchPbPb[cb] - (0.95)*meanNchPbPb[cb]);
1234       else xAxis_e[cb]=fabs(meanNchpPb[cb] - (0.95)*meanNchpPb[cb]);
1235     }
1236   }
1237
1238   TGraphErrors *grLambdaSys_PbPb=new TGraphErrors(20,meanNchPbPb,yAxisPbPb,xAxis_e,yAxisPbPb_eL);
1239   TGraphAsymmErrors *grLambdaSys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxis_e,xAxis_e, yAxispPb_eL,yAxispPb_eH);
1240   TGraphAsymmErrors *grLambdaSys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxis_e,xAxis_e, yAxispp_eL,yAxispp_eH);
1241   grLambdaSys_pp->SetMarkerSize(0); grLambdaSys_pp->SetFillStyle(1000); grLambdaSys_pp->SetFillColor(kBlue-10);
1242   grLambdaSys_pPb->SetMarkerSize(0); grLambdaSys_pPb->SetFillStyle(1000); grLambdaSys_pPb->SetFillColor(kRed-10);
1243   grLambdaSys_PbPb->SetMarkerSize(0); grLambdaSys_PbPb->SetFillStyle(1000); grLambdaSys_PbPb->SetFillColor(kGray);
1244   grLambdaSys_pp->SetMarkerColor(kBlue-10); grLambdaSys_pPb->SetMarkerColor(kRed-10); grLambdaSys_PbPb->SetMarkerColor(kGray);
1245   // C2 sys
1246   TH1D *LambdaC2PbPb=(TH1D*)Parameters_C2[0][FitType][KT3Bin][5]->Clone();
1247   TH1D *LambdaC2pPb=(TH1D*)Parameters_C2[1][FitType][KT3Bin][5]->Clone();
1248   TH1D *LambdaC2pp=(TH1D*)Parameters_C2[2][FitType][KT3Bin][5]->Clone();
1249   for(int cb=0; cb<20; cb++){// 2-particle
1250     int binPbPb = LambdaC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
1251     int binpPb = LambdaC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]);
1252     int binpp = LambdaC2pp->GetXaxis()->FindBin(meanNchpp[cb]);
1253     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
1254     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
1255     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
1256     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
1257     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
1258     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
1259     //
1260     if(LambdaC2PbPb->GetBinContent(binPbPb)==0) {yAxisPbPb[cb]=100; yAxisPbPb_eL[cb]=100; yAxisPbPb_eH[cb]=100;}// errors were 100
1261     else {
1262       yAxisPbPb[cb]=LambdaC2PbPb->GetBinContent(binPbPb);
1263       if(cb>=13) yAxisPbPb_eL[cb]=LambdasysPbPb_L;
1264       else yAxisPbPb_eL[cb]=LambdasysPbPb_H;
1265       yAxisPbPb_eH[cb]=LambdasysPbPb_H;
1266     }
1267     //
1268     if(LambdaC2pPb->GetBinContent(binpPb)==0) {yAxispPb[cb]=100; yAxispPb_eL[cb]=100; yAxispPb_eH[cb]=100;}
1269     else {yAxispPb[cb]=LambdaC2pPb->GetBinContent(binpPb); yAxispPb_eL[cb]=LambdasyspPb_L; yAxispPb_eH[cb]=LambdasyspPb_H;}
1270     //
1271     if(LambdaC2pp->GetBinContent(binpp)==0) {yAxispp[cb]=100; yAxispp_eL[cb]=100; yAxispp_eH[cb]=100;}
1272     else {yAxispp[cb]=LambdaC2pp->GetBinContent(binpp); yAxispp_eL[cb]=Lambdasyspp_L; yAxispp_eH[cb]=Lambdasyspp_H;}
1273     //xAxis_e[cb]=10000;
1274   }
1275   TGraphAsymmErrors *grLambdaC2Sys_PbPb=new TGraphAsymmErrors(20,meanNchPbPb,yAxisPbPb,xAxisC2_e,xAxisC2_e, yAxisPbPb_eL,yAxisPbPb_eH);
1276   TGraphAsymmErrors *grLambdaC2Sys_pPb=new TGraphAsymmErrors(20,meanNchpPb,yAxispPb, xAxisC2_e,xAxisC2_e, yAxispPb_eL,yAxispPb_eH);
1277   TGraphAsymmErrors *grLambdaC2Sys_pp=new TGraphAsymmErrors(20,meanNchpp,yAxispp, xAxisC2_e,xAxisC2_e, yAxispp_eL,yAxispp_eH);
1278   grLambdaC2Sys_pp->SetMarkerSize(0); grLambdaC2Sys_pp->SetFillStyle(3001); grLambdaC2Sys_pp->SetFillColor(0);
1279   grLambdaC2Sys_pPb->SetMarkerSize(0); grLambdaC2Sys_pPb->SetFillStyle(3001); grLambdaC2Sys_pPb->SetFillColor(0);
1280   grLambdaC2Sys_PbPb->SetMarkerSize(0); grLambdaC2Sys_PbPb->SetFillStyle(3001); grLambdaC2Sys_PbPb->SetFillColor(0);
1281   grLambdaC2Sys_pp->SetLineColor(4); grLambdaC2Sys_pPb->SetLineColor(2); grLambdaC2Sys_PbPb->SetLineColor(1);
1282   grLambdaC2Sys_pp->SetLineWidth(2.*grLambdaC2Sys_pp->GetLineWidth());
1283   grLambdaC2Sys_pPb->SetLineWidth(2.*grLambdaC2Sys_pPb->GetLineWidth());
1284   grLambdaC2Sys_PbPb->SetLineWidth(2.*grLambdaC2Sys_PbPb->GetLineWidth());
1285   //
1286
1287
1288   grLambdaSys_pp->Draw("E2 p");
1289   grLambdaSys_pPb->Draw("E2 p");
1290   grLambdaSys_PbPb->Draw("E2 p");
1291   //
1292   LambdaPbPb->Draw("same");
1293   LambdapPb->Draw("same");
1294   Lambdapp->Draw("same");
1295   //
1296   LambdaC2PbPb->Draw("same");
1297   LambdaC2pPb->Draw("same");
1298   LambdaC2pp->Draw("same");
1299
1300   grLambdaC2Sys_pp->Draw("|| p");
1301   grLambdaC2Sys_pPb->Draw("|| p");
1302   grLambdaC2Sys_PbPb->Draw("|| p");
1303
1304   // print radii and lambda
1305   cout.precision(3);
1306   cout<<"Radii and Lambda data"<<endl;
1307   cout<<"p--p:"<<endl;
1308   // new way for HEP
1309   for(int cb=19; cb>=0; cb--){
1310     int binpp = RadiiC2pp->GetXaxis()->FindBin(meanNchpp[cb]);
1311     if(Radiipp->GetBinContent(binpp)==0) continue;
1312     cout<<meanNchpp[cb]-(0.05)*meanNchpp[cb]<<" TO "<<meanNchpp[cb]+(0.05)*meanNchpp[cb]<<"; "<<RadiiC2pp->GetBinContent(binpp)<<" +- "<<RadiiC2pp->GetBinError(binpp)<<" (DSYS=+"<<grRadiiC2Sys_pp->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_pp->GetErrorYlow(cb)<<"); ";
1313     cout<<LambdaC2pp->GetBinContent(binpp)<<" +- "<<LambdaC2pp->GetBinError(binpp)<<" (DSYS=+"<<grLambdaC2Sys_pp->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_pp->GetErrorYlow(cb)<<"); ";
1314     cout<<Radiipp->GetBinContent(binpp)<<" +- "<<Radiipp->GetBinError(binpp)<<" (DSYS="<<grRadiiSys_pp->GetErrorY(cb)<<"); ";
1315     cout<<Lambdapp->GetBinContent(binpp)<<" +- "<<Lambdapp->GetBinError(binpp)<<" (DSYS="<<grLambdaSys_pp->GetErrorY(cb)<<");"<<endl;
1316     // Edgeworth lambdas
1317     //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;
1318   }
1319   cout<<endl;
1320
1321   cout<<"p--Pb:"<<endl;
1322   for(int cb=19; cb>=0; cb--){
1323     int binpPb = RadiiC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]);
1324     if(RadiipPb->GetBinContent(binpPb)==0) continue;
1325     cout<<meanNchpPb[cb]-(0.05)*meanNchpPb[cb]<<" TO "<<meanNchpPb[cb]+(0.05)*meanNchpPb[cb]<<"; "<<RadiiC2pPb->GetBinContent(binpPb)<<" +- "<<RadiiC2pPb->GetBinError(binpPb)<<" (DSYS=+"<<grRadiiC2Sys_pPb->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_pPb->GetErrorYlow(cb)<<"); ";
1326     cout<<LambdaC2pPb->GetBinContent(binpPb)<<" +- "<<LambdaC2pPb->GetBinError(binpPb)<<" (DSYS=+"<<grLambdaC2Sys_pPb->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_pPb->GetErrorYlow(cb)<<"); ";
1327     cout<<RadiipPb->GetBinContent(binpPb)<<" +- "<<RadiipPb->GetBinError(binpPb)<<" (DSYS="<<grRadiiSys_pPb->GetErrorY(cb)<<"); ";
1328     cout<<LambdapPb->GetBinContent(binpPb)<<" +- "<<LambdapPb->GetBinError(binpPb)<<" (DSYS="<<grLambdaSys_pPb->GetErrorY(cb)<<");"<<endl;
1329     // Gaussian lambdas
1330     //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;
1331     // Edgeworth lambdas
1332     //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;
1333   }
1334   cout<<endl;
1335   cout<<"Pb--Pb:"<<endl;
1336   for(int cb=19; cb>=0; cb--){
1337     int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
1338     if(RadiiPbPb->GetBinContent(binPbPb)==0) continue;
1339     cout<<meanNchPbPb[cb]-(0.05)*meanNchPbPb[cb]<<" TO "<<meanNchPbPb[cb]+(0.05)*meanNchPbPb[cb]<<"; "<<RadiiC2PbPb->GetBinContent(binPbPb)<<" +- "<<RadiiC2PbPb->GetBinError(binPbPb)<<" (DSYS=+"<<grRadiiC2Sys_PbPb->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_PbPb->GetErrorYlow(cb)<<"); ";
1340     cout<<LambdaC2PbPb->GetBinContent(binPbPb)<<" +- "<<LambdaC2PbPb->GetBinError(binPbPb)<<" (DSYS=+"<<grLambdaC2Sys_PbPb->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_PbPb->GetErrorYlow(cb)<<"); ";
1341     cout<<RadiiPbPb->GetBinContent(binPbPb)<<" +- "<<RadiiPbPb->GetBinError(binPbPb)<<" (DSYS="<<grRadiiSys_PbPb->GetErrorY(cb)<<"); ";
1342     cout<<LambdaPbPb->GetBinContent(binPbPb)<<" +- "<<LambdaPbPb->GetBinError(binPbPb)<<" (DSYS="<<grLambdaSys_PbPb->GetErrorY(cb)<<");"<<endl;
1343     //
1344     //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;
1345   }
1346   cout<<endl;
1347
1348
1349
1350
1351   can3->cd();
1352   TPad *pad3_2 = new TPad("pad3_2","pad3_2",0.0,0.0,1.,1.);
1353   pad3_2->SetFillStyle(0);
1354   pad3_2->Draw();
1355   pad3_2->cd();
1356   TLatex *RinvTitle;
1357   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)");
1358   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)");
1359   else RinvTitle=new TLatex(0.062,0.61,"(#font[12]{R}^{Exp}_{inv} or #font[12]{R}^{Exp}_{inv,3}) / #sqrt{#pi} (fm)");
1360   RinvTitle->SetNDC();
1361   RinvTitle->SetTextFont(TextFont);
1362   RinvTitle->SetTextSize(SizeTitle);
1363   RinvTitle->SetTextAngle(90);
1364   RinvTitle->Draw("same");
1365   TLatex *LambdaTitle;
1366   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
1367   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
1368   else LambdaTitle=new TLatex(0.064,0.31,"#lambda^{Exp}_{e} or #lambda^{Exp}_{e,3}");// 0.064,0.33
1369   LambdaTitle->SetNDC();
1370   LambdaTitle->SetTextFont(TextFont);
1371   LambdaTitle->SetTextSize(SizeTitle);
1372   LambdaTitle->SetTextAngle(90);
1373   LambdaTitle->Draw("same");
1374
1375
1376   if(SaveFiles && FitType==0) can3->SaveAs("ThreePionFitParametersGauss.eps");
1377   if(SaveFiles && FitType==1) can3->SaveAs("ThreePionFitParametersEW.eps");
1378
1379
1380
1381   ///////////////////////////////////////////////////////////////////////////////
1382   ///////////////////////////////////////////////////////////////////////////////
1383   // kappa plots
1384   /*
1385   TCanvas *can7 = new TCanvas("can7", "can7",1700,700,600,600);// 11,53,700,500
1386   can7->SetHighLightColor(2);
1387   gStyle->SetOptFit(0);// 0111 to show fit stat box
1388   can7->SetFillColor(0);//10
1389   can7->SetBorderMode(0);
1390   can7->SetBorderSize(2);
1391   can7->SetFrameFillColor(0);
1392   can7->SetFrameBorderMode(0);
1393   can7->SetFrameBorderMode(0);
1394   can7->cd();
1395   TPad *pad7 = new TPad("pad7","pad7",0.0,0.0,1.,1.);
1396   gPad->SetGridx(0);
1397   gPad->SetGridy(1);
1398   pad7->SetTopMargin(0.02);//0.05
1399   pad7->SetRightMargin(0.02);//3e-2
1400   pad7->SetBottomMargin(0.1);//0.12
1401   pad7->SetLeftMargin(0.1);//0.12
1402   pad7->Draw();
1403   pad7->cd();
1404   gPad->SetLogx();
1405   gPad->SetGridy(1);
1406   TLegend *legend8 = new TLegend(.2,.70, .4,.95,NULL,"brNDC");//.45 or .4 for x1
1407   legend8->SetBorderSize(0);
1408   legend8->SetFillColor(0);
1409   legend8->SetTextFont(TextFont);
1410   legend8->SetTextSize(SizeLegend);
1411   // CollType, Gaussian/EW, EDbin, Parameter#
1412   int paramNum=4;
1413   Parameters_c3[0][1][KT3Bin][paramNum]->GetXaxis()->SetTitleOffset(1.2); Parameters_c3[0][1][KT3Bin][paramNum]->GetYaxis()->SetTitleOffset(1.4);
1414   if(paramNum==3) {Parameters_c3[0][1][KT3Bin][paramNum]->SetMinimum(-0.1); Parameters_c3[0][1][KT3Bin][paramNum]->SetMaximum(.4);}
1415   if(paramNum==4) {Parameters_c3[0][1][KT3Bin][paramNum]->SetMinimum(-0.1); Parameters_c3[0][1][KT3Bin][paramNum]->SetMaximum(1.0);}
1416   Parameters_c3[0][1][KT3Bin][paramNum]->Draw();
1417   Parameters_c3[1][1][KT3Bin][paramNum]->Draw("same");
1418   Parameters_c3[2][1][KT3Bin][paramNum]->Draw("same");
1419   legend8->AddEntry(Parameters_c3[2][1][KT3Bin][paramNum],"pp #sqrt{s}=7 TeV","p");
1420   legend8->AddEntry(Parameters_c3[1][1][KT3Bin][paramNum],"p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV","p");
1421   legend8->AddEntry(Parameters_c3[0][1][KT3Bin][paramNum],"Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV","p");
1422   //legend8->Draw("same");
1423   //for(int ct=0; ct<3; ct++){
1424     //for(int cb=0; cb<20; cb++){
1425       //int bin = 0;
1426       //if(ct==0) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->GetXaxis()->FindBin(meanNchPbPb[cb]);
1427       //else if(ct==1) bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpPb[cb]);
1428       //else bin = Parameters_c3[ct][1][KT3Bin][paramNum]->FindBin(meanNchpp[cb]);
1429       //
1430       //if(Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin) >0) {
1431         //cout<<Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinContent(bin)<<", ";
1432         //cout<<Parameters_c3[ct][1][KT3Bin][paramNum]->GetBinError(bin)<<", ";
1433       //}else cout<<0<<", ";
1434
1435     //}
1436     //cout<<endl;
1437     //}
1438
1439
1440
1441   TH1D *Combined_kappaPlot_1=(TH1D*)Parameters_c3[0][1][KT3Bin][paramNum]->Clone();
1442   Combined_kappaPlot_1->Add(Parameters_c3[1][1][KT3Bin][paramNum]);
1443   Combined_kappaPlot_1->Add(Parameters_c3[2][1][KT3Bin][paramNum]);
1444
1445   //TF1 *Fit_kappa3_PbPb=new TF1("Fit_kappa3_PbPb","[0]+[1]*log(x)",2,3000);
1446   //Fit_kappa3_PbPb->SetParameter(0, 0.05);
1447   //Fit_kappa3_PbPb->SetParameter(1, -0.01);
1448   //Fit_kappa3_PbPb->SetLineColor(1);
1449
1450   //
1451   //TF1 *Fit_kappa3_pp=new TF1("Fit_kappa3_pp","[0]+[1]*log(x)",1,3000);
1452   //Fit_kappa3_pp->SetParameter(0, 0.05);
1453   //Fit_kappa3_pp->SetParameter(1, -0.01);
1454   //Fit_kappa3_pp->SetLineColor(3);
1455   //Combined_kappaPlot_1->Fit(Fit_kappa3_pp,"IMEN","",2,80);
1456   //Fit_kappa3_pp->Draw("same");
1457   //
1458   TF1 *Fit_kappa4_PbPb=new TF1("Fit_kappa4_PbPb","pol0",2,3000);
1459   Combined_kappaPlot_1->Fit(Fit_kappa4_PbPb,"IMEN","",2,2000);
1460   Fit_kappa4_PbPb->Draw("same");
1461   */
1462
1463
1464   ////////////////////////////////////////////////////
1465   ////////////////////////////////////////////////////
1466   // Radii ratios
1467   if(NchOneThirdAxis){
1468     TCanvas *can6 = new TCanvas("can6", "can6",1700,700,600,600);// 11,53,700,500
1469     can6->SetHighLightColor(2);
1470     gStyle->SetOptFit(0);// 0111 to show fit stat box
1471     can6->SetFillColor(0);//10
1472     can6->SetBorderMode(0);
1473     can6->SetBorderSize(2);
1474     can6->SetFrameFillColor(0);
1475     can6->SetFrameBorderMode(0);
1476     can6->SetFrameBorderMode(0);
1477     can6->cd();
1478     TPad *pad6 = new TPad("pad6","pad6",0.0,0.0,1.,1.);
1479     gPad->SetGridx(0);
1480     gPad->SetGridy(0);
1481     pad6->SetTopMargin(0.0);//0.05
1482     pad6->SetRightMargin(0.0);//3e-2
1483     pad6->SetBottomMargin(0.0);//0.12
1484     pad6->SetLeftMargin(0.0);//0.12
1485     pad6->Draw();
1486     pad6->cd();
1487     TLegend *legend8 = new TLegend(.52,.3, .9,.5,NULL,"brNDC");//.45 or .4 for x1
1488     legend8->SetBorderSize(0);
1489     legend8->SetFillColor(0);
1490     legend8->SetTextFont(TextFont);
1491     legend8->SetTextSize(SizeLegend);
1492     gPad->SetRightMargin(0.01); gPad->SetLeftMargin(0.14);
1493     gPad->SetBottomMargin(0.14); gPad->SetTopMargin(0.02);
1494     gPad->SetTickx(); gPad->SetTicky();
1495
1496     TH1D *Ratio_pPb_to_pp=(TH1D*)RadiipPb->Clone();
1497     TH1D *Ratio_PbPb_to_pPb=(TH1D*)RadiiC2PbPb->Clone();
1498     double avgRatio_pPb_to_pp=0, avgRatioSq_pPb_to_pp=0, avgRatioStat_pPb_to_pp=0,  avgRatioEn_pPb_to_pp=0;
1499     double avgRatio_PbPb_to_pPb=0, avgRatioSq_PbPb_to_pPb=0, avgRatioStat_PbPb_to_pPb=0, avgRatioEn_PbPb_to_pPb=0;
1500     for(int cb=0; cb<20; cb++){// 3-particle
1501       int binPbPb = RadiiPbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
1502       int binpPb = RadiipPb->GetXaxis()->FindBin(meanNchpPb[cb]);
1503       //
1504       Ratio_pPb_to_pp->SetBinContent(binpPb, Ratio_pPb_to_pp->GetBinContent(binpPb) / ppLine->Eval(meanNchpPb[cb]));
1505       Ratio_PbPb_to_pPb->SetBinContent(binPbPb, Ratio_PbPb_to_pPb->GetBinContent(binPbPb) / pPbLine->Eval(meanNchPbPb[cb]));
1506       //
1507       if(cb<=18 && cb>=14){
1508         avgRatio_pPb_to_pp += Ratio_pPb_to_pp->GetBinContent(binpPb);
1509         avgRatioSq_pPb_to_pp += pow(Ratio_pPb_to_pp->GetBinContent(binpPb),2);
1510         avgRatioStat_pPb_to_pp += pow(Ratio_pPb_to_pp->GetBinError(binpPb),2);
1511         avgRatioEn_pPb_to_pp++;
1512       }
1513       if(cb<=15 && cb>=13){
1514         avgRatio_PbPb_to_pPb += Ratio_PbPb_to_pPb->GetBinContent(binPbPb);
1515         avgRatioSq_PbPb_to_pPb += pow(Ratio_PbPb_to_pPb->GetBinContent(binPbPb),2);
1516         avgRatioStat_PbPb_to_pPb += pow(Ratio_PbPb_to_pPb->GetBinError(binPbPb),2);
1517         avgRatioEn_PbPb_to_pPb++;
1518       }
1519     }
1520     Ratio_pPb_to_pp->SetMinimum(0.9); Ratio_pPb_to_pp->SetMaximum(1.65);
1521     Ratio_pPb_to_pp->GetYaxis()->SetTitle("Radius ratio");
1522     Ratio_pPb_to_pp->GetXaxis()->SetLabelSize(SizeLabel); Ratio_pPb_to_pp->GetYaxis()->SetLabelSize(SizeLabel);
1523     Ratio_pPb_to_pp->GetXaxis()->SetTitleSize(SizeTitle); Ratio_pPb_to_pp->GetYaxis()->SetTitleSize(SizeTitle);
1524     Ratio_pPb_to_pp->GetXaxis()->SetTitleOffset(1.0);
1525     Ratio_pPb_to_pp->GetYaxis()->SetTitleOffset(1.2);
1526     Ratio_pPb_to_pp->GetXaxis()->SetNdivisions(808);
1527     Ratio_pPb_to_pp->GetYaxis()->SetNdivisions(505);
1528     Ratio_pPb_to_pp->Draw();
1529     Ratio_PbPb_to_pPb->Draw("same");
1530     legend8->AddEntry(Ratio_pPb_to_pp,"p-Pb over pp","p");
1531     legend8->AddEntry(Ratio_PbPb_to_pPb,"Pb-Pb over p-Pb","p");
1532     legend8->Draw("same");
1533     Unity->Draw("same");
1534
1535     double avgStat_pPb_to_pp = sqrt(avgRatioStat_pPb_to_pp/avgRatioEn_pPb_to_pp);
1536     double RMS_pPb_to_pp = sqrt( (avgRatioSq_pPb_to_pp/avgRatioEn_pPb_to_pp - pow(avgRatio_pPb_to_pp/avgRatioEn_pPb_to_pp,2)) / avgRatioEn_pPb_to_pp );
1537     double avgStat_PbPb_to_pPb = sqrt(avgRatioStat_PbPb_to_pPb/avgRatioEn_PbPb_to_pPb);
1538     double RMS_PbPb_to_pPb = sqrt( (avgRatioSq_PbPb_to_pPb/avgRatioEn_PbPb_to_pPb - pow(avgRatio_PbPb_to_pPb/avgRatioEn_PbPb_to_pPb,2)) / avgRatioEn_PbPb_to_pPb );
1539     cout.precision(4);
1540     cout<<"avg Ratio of pPb to pp = "<<avgRatio_pPb_to_pp/avgRatioEn_pPb_to_pp<<" +- "<<sqrt(pow(avgStat_pPb_to_pp,2) + pow(RMS_pPb_to_pp,2))<<endl;
1541     cout<<"avg Ratio of PbPb to pPb = "<<avgRatio_PbPb_to_pPb/avgRatioEn_PbPb_to_pPb<<" +- "<<sqrt(pow(avgStat_PbPb_to_pPb,2) + pow(RMS_PbPb_to_pPb,2))<<endl;
1542   }
1543
1544
1545   if(KT3Bin>0) {cout<<"Skip the rest for this setting"<<endl; return;}
1546
1547   ////////////////////////////////////////////////////
1548   ////////////////////////////////////////////////////
1549   // Correlation functions and Monte Carlo
1550   TCanvas *can4 = (TCanvas*)make_canvas("can4","can4",3,2,0,900,600);
1551   can4->Draw();
1552
1553   TLegend *legend5[6];
1554   legend5[0] = new TLegend(.3,.52, .97,.99,NULL,"brNDC");//.45 or .4 for x1
1555   legend5[0]->SetBorderSize(0);
1556   legend5[0]->SetFillColor(0);
1557   legend5[0]->SetTextFont(TextFont);
1558   legend5[1]=(TLegend*)legend5[0]->Clone();
1559   legend5[2]=(TLegend*)legend5[0]->Clone();
1560   legend5[3]=(TLegend*)legend5[0]->Clone();
1561   legend5[4]=(TLegend*)legend5[0]->Clone();
1562   legend5[5]=(TLegend*)legend5[0]->Clone();
1563   TLegend *legendFitTypes = (TLegend*)legend5[0]->Clone();
1564
1565
1566   double HIJING_c3_SC_K1_M3[15]={0, 0.851168, 0.979088, 1.0209, 0.976797, 1.01555, 0.992262, 1.00773, 0.991634, 0.991504, 0.997317, 0.993006, 0.99694, 0.999369, 0.998404};
1567   double HIJING_c3_SC_K1_M3_e[15]={0, 0.546937, 0.118551, 0.0436675, 0.0226652, 0.0139659, 0.00906562, 0.00649369, 0.00488794, 0.00380819, 0.00306916, 0.00255166, 0.00219781, 0.00235171, 0.00292962};
1568   double HIJING_c3_MC_K1_M3[15]={0, 0.886712, 1.02583, 0.985831, 1.00453, 1.01572, 1.00153, 0.991872, 0.997636, 0.997151, 0.996838, 0.999562, 0.998487, 0.996162, 1.001};
1569   double HIJING_c3_MC_K1_M3_e[15]={0, 0.190786, 0.0527107, 0.0220311, 0.0119954, 0.00756783, 0.00499146, 0.00360927, 0.0027377, 0.00214376, 0.00173608, 0.00144723, 0.00124891, 0.00133898, 0.0016771};
1570   TH1D *HIJING_c3_SC=(TH1D*)c3[1][1][0][0][17]->Clone();// choose any one that exists to copy histo attributes
1571   TH1D *HIJING_c3_MC=(TH1D*)c3[1][1][0][0][17]->Clone();// choose any one that exists to copy histo attributes
1572   for(int i=1; i<=15; i++){
1573     HIJING_c3_SC->SetBinContent(i, HIJING_c3_SC_K1_M3[i-1]);
1574     HIJING_c3_SC->SetBinError(i, HIJING_c3_SC_K1_M3_e[i-1]);
1575     HIJING_c3_MC->SetBinContent(i, HIJING_c3_MC_K1_M3[i-1]);
1576     HIJING_c3_MC->SetBinError(i, HIJING_c3_MC_K1_M3_e[i-1]);
1577   }
1578   cout.precision(4);
1579   //
1580   for(int padNum=1; padNum<=6; padNum++){
1581
1582     can4->cd(padNum);
1583     if(padNum==3 || padNum==6) gPad->SetRightMargin(0.005);
1584     float SF_6pannel=2;
1585     double SF_correction=1.0;
1586     if(padNum==3) SF_correction=1.1;
1587     if(padNum==4) SF_correction=0.8;
1588     if(padNum==5 || padNum==6) SF_correction=0.92;
1589
1590     //
1591     int System_proof=0;
1592     int ChComb_proof=0;
1593     int Mb_proof=0;
1594     if(padNum==1) {System_proof=2; ChComb_proof=0; Mb_proof=18;}
1595     if(padNum==2) {System_proof=1; ChComb_proof=0; Mb_proof=14;}
1596     if(padNum==3) {System_proof=0; ChComb_proof=0; Mb_proof=3;}
1597     if(padNum==4) {System_proof=2; ChComb_proof=1; Mb_proof=18;}
1598     if(padNum==5) {System_proof=1; ChComb_proof=1; Mb_proof=14;}
1599     if(padNum==6) {System_proof=0; ChComb_proof=1; Mb_proof=3;}
1600
1601     // print out data points
1602     for(int binN=1; binN<=C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetNbinsX(); binN++){
1603       cout<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinLowEdge(binN)<<" TO "<<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)<<" (DSYS="<<C3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<"); "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" +- "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<" (DSYS="<<c3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<"); ";
1604       if(System_proof==0){
1605         cout<<HIJING_c3_SC->GetBinContent(binN)<<" +- "<<HIJING_c3_SC->GetBinError(binN)<<"; - ;"<<endl;
1606       }else{
1607         cout<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" +- "<<c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<"; - ;"<<endl;
1608       }
1609     }
1610     cout<<endl;
1611     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->SetMinimum(0.9);
1612     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->SetMaximum(3.4);// 3.4
1613     //
1614     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetTitle("#font[12]{C}_{3} or #font[12]{#bf{c}}_{3} ");
1615     if(padNum<=5){
1616       C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetTitleOffset(10);
1617     }else C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetTitleOffset(0.88);
1618     if(padNum==1) C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetTitleOffset(0.55);
1619     else C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetTitleOffset(10.);
1620     if(padNum>=5) C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetLabelOffset(-.0);
1621
1622     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetNdivisions(504);
1623     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetNdivisions(504);
1624     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetTitleSize(SizeTitle*SF_6pannel*SF_correction);
1625     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetTitleSize(SizeTitle*SF_6pannel*SF_correction);
1626     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetLabelSize(SizeTitle*SF_6pannel*SF_correction);
1627     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetYaxis()->SetLabelSize(SizeTitle*SF_6pannel*SF_correction);
1628     double Q3Limit;
1629     if(System_proof==1 || System_proof==2) Q3Limit = 0.49;
1630     else Q3Limit = 0.1099;
1631     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->SetRangeUser(0,Q3Limit);
1632     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy();
1633
1634
1635     C3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy("E2 same");
1636     c3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy("E2 same");
1637     C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy("same");
1638     c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy("same");
1639     if(System_proof!=0) c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof]->DrawCopy("same");// MC cumulants
1640     else{
1641       if(ChComb_proof==0) HIJING_c3_SC->DrawCopy("same");
1642       else HIJING_c3_MC->DrawCopy("same");
1643     }
1644     if(padNum<=3){
1645       legend5[padNum-1]->AddEntry(C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof],"#font[12]{C}_{3}^{#pm#pm#pm}","p");
1646       legend5[padNum-1]->AddEntry(c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof],"#font[12]{#bf{c}}_{3}^{#pm#pm#pm}","p");
1647     }else{
1648       legend5[padNum-1]->AddEntry(C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof],"#font[12]{C}_{3}^{#pm#pm#mp}","p");
1649       legend5[padNum-1]->AddEntry(c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof],"#font[12]{#bf{c}}_{3}^{#pm#pm#mp}","p");
1650       if(System_proof==0) legend5[padNum-1]->AddEntry(HIJING_c3_MC,"HIJING #font[12]{#bf{c}}_{3}","p");
1651       if(System_proof==1) legend5[padNum-1]->AddEntry(c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof],"DPMJET #font[12]{#bf{c}}_{3}","p");
1652       if(System_proof==2) legend5[padNum-1]->AddEntry(c3[System_proof][1][ChComb_proof][KT3Bin][Mb_proof],"PYTHIA #font[12]{#bf{c}}_{3}","p");
1653     }
1654
1655     if(ChComb_proof==0) {
1656       c3_fit[System_proof][0][KT3Bin][Mb_proof]->Draw("c same");// Gauss
1657       gr_c3Spline[System_proof][KT3Bin][Mb_proof]->Draw("c same");// EW with spline for mid-q and high q
1658       gr_c3SplineExpFit[System_proof][KT3Bin][Mb_proof]->Draw("c same");// Exp
1659       //c3_fit[System_proof][1][KT3Bin][Mb_proof]->Draw("c same");// old approximation
1660       if(padNum==3){
1661         legendFitTypes->AddEntry(c3_fit[System_proof][0][KT3Bin][Mb_proof],"Gaussian","l");
1662         legendFitTypes->AddEntry(c3_fit[System_proof][1][KT3Bin][Mb_proof],"Edgeworth","l");
1663         legendFitTypes->AddEntry(c3_fit[System_proof][2][KT3Bin][Mb_proof],"Exponential","l");
1664       }
1665     }
1666
1667     TLatex *CTLabel;
1668     if(System_proof==0) CTLabel = new TLatex(0.12,0.9,"Pb-Pb #sqrt{#font[12]{s}_{NN}}=2.76 TeV");// 0.003,.988
1669     if(System_proof==1) CTLabel = new TLatex(0.15,0.9,"p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV");// 0.003,.988
1670     if(System_proof==2) CTLabel = new TLatex(0.4,0.9,"pp #sqrt{s}=7 TeV");// 0.003,.988
1671     CTLabel->SetNDC();
1672     CTLabel->SetTextFont(TextFont);
1673     CTLabel->SetTextSize(SizeSpecif*SF_6pannel*SF_correction);
1674     if(padNum>=4) CTLabel->Draw("same");
1675
1676     TString *nameCh=new TString("#LT#font[12]{N}_{ch}#GT = ");
1677     float Nch=1;
1678     if(System_proof==0) Nch = meanNchPbPb[Mb_proof];
1679     else if(System_proof==1) Nch = meanNchpPb[Mb_proof];
1680     else Nch = meanNchpp[Mb_proof];
1681     *nameCh += int(Nch + 0.5);
1682     nameCh->Append(" #pm ");
1683     float SysPercent = 0.05;
1684     int SigFig=0;
1685     if(SysPercent*Nch < 1) {nameCh->Append("0."); SigFig=int(10*SysPercent*Nch + 0.5);}
1686     else {SigFig=int(SysPercent*Nch + 0.5);}
1687
1688     *nameCh += SigFig;
1689     TLatex *MLabel;
1690     if(padNum==1) MLabel = new TLatex(0.45,0.6,"#LT#font[12]{N}_{ch}#GT = 8.6 #pm 0.4");// was nameCh->Data()
1691     if(padNum==2) MLabel = new TLatex(0.4,0.6,nameCh->Data());
1692     if(padNum==3) MLabel = new TLatex(0.4,0.6,nameCh->Data());
1693
1694     MLabel->SetNDC();
1695     MLabel->SetTextFont(TextFont);
1696     MLabel->SetTextSize(SizeSpecif*SF_6pannel*SF_correction);
1697     if(padNum<=3) MLabel->Draw("same");
1698
1699
1700     legend5[padNum-1]->SetTextSize(SizeLegend*SF_6pannel*SF_correction);
1701     if(padNum==1) {legend5[padNum-1]->SetX1(0.45); legend5[padNum-1]->SetY1(0.69);}
1702     if(padNum==2) {legend5[padNum-1]->SetX1(0.32); legend5[padNum-1]->SetY1(0.69);}
1703     if(padNum==3) {
1704       legend5[padNum-1]->SetX1(0.32); legend5[padNum-1]->SetY1(0.69);
1705       legendFitTypes->SetX1(0.44);
1706       legendFitTypes->SetY1(0.22); legendFitTypes->SetY2(0.56);
1707       legendFitTypes->Draw("same");
1708     }
1709     if(padNum==4) {legend5[padNum-1]->SetX1(0.45); legend5[padNum-1]->SetY1(0.45); legend5[padNum-1]->SetY2(0.85);}
1710     if(padNum==5) {legend5[padNum-1]->SetX1(0.32); legend5[padNum-1]->SetY1(0.45); legend5[padNum-1]->SetY2(0.85);}
1711     if(padNum==6) {legend5[padNum-1]->SetX1(0.32); legend5[padNum-1]->SetY1(0.45); legend5[padNum-1]->SetY2(0.85);}
1712     legend5[padNum-1]->Draw("same");
1713   }
1714
1715
1716
1717   can4->cd();
1718
1719   TPad *pad4_2 = new TPad("pad4_2","pad4_2",0.0,0.0,1.,1.);
1720   pad4_2->SetFillStyle(0);
1721   pad4_2->Draw();
1722   pad4_2->cd();
1723   TBox *CoverUp1 = new TBox(0.35,0.05,0.42,.115);
1724   CoverUp1->SetFillColor(0);
1725   CoverUp1->Draw();
1726   TBox *CoverUp2 = new TBox(0.66,0.05,0.73,.115);
1727   CoverUp2->SetFillColor(0);
1728   CoverUp2->Draw();
1729
1730
1731
1732
1733
1734   ////////////////////////////////////////////////////
1735   ////////////////////////////////////////////////////
1736   // 2 system correlation function comparison
1737   TCanvas *can5 = new TCanvas("can5", "can5",1700,700,600,600);// 11,53,700,500
1738   can5->SetHighLightColor(2);
1739   gStyle->SetOptFit(0);// 0111 to show fit stat box
1740   can5->SetFillColor(0);//10
1741   can5->SetBorderMode(0);
1742   can5->SetBorderSize(2);
1743   can5->SetFrameFillColor(0);
1744   can5->SetFrameBorderMode(0);
1745   can5->SetFrameBorderMode(0);
1746   can5->cd();
1747   TPad *pad5 = new TPad("pad5","pad5",0.0,0.0,1.,1.);
1748   gPad->SetGridx(0);
1749   gPad->SetGridy(0);
1750   pad5->SetTopMargin(0.0);//0.05
1751   pad5->SetRightMargin(0.0);//3e-2
1752   pad5->SetBottomMargin(0.0);//0.12
1753   pad5->SetLeftMargin(0.0);//0.12
1754   pad5->Draw();
1755   pad5->cd();
1756   TLegend *legend6 = new TLegend(.42,.6, .9,.95,NULL,"brNDC");//.45 or .4 for x1
1757   legend6->SetBorderSize(0);
1758   legend6->SetFillColor(0);
1759   legend6->SetTextFont(TextFont);
1760   legend6->SetTextSize(SizeLegend);
1761   TLegend *legend7 = new TLegend(.66,.36, .97,.53,NULL,"brNDC");//.45 or .4 for x1
1762   legend7->SetBorderSize(0);
1763   legend7->SetFillColor(0);
1764   legend7->SetTextFont(TextFont);
1765   legend7->SetTextSize(SizeLegend);
1766   //
1767   gPad->SetRightMargin(0.01); gPad->SetLeftMargin(0.10);
1768   gPad->SetBottomMargin(0.14); gPad->SetTopMargin(0.02);
1769   gPad->SetTickx(); gPad->SetTicky();
1770   //
1771   int KT3Bin_CorrComp=0;
1772   int Mbin_SysComp_PbPb=12;// 12
1773   int Mbin_SysComp_pPb;
1774   int Mbin_SysComp_pp=15;// 15
1775   if(pp_pPb_Comp) Mbin_SysComp_pPb=16;// 16
1776   else Mbin_SysComp_pPb=12;// 12
1777
1778   TH1D *c3_PbPb=(TH1D*)c3[0][0][0][KT3Bin_CorrComp][Mbin_SysComp_PbPb]->Clone();
1779   TH1D *c3_pPb=(TH1D*)c3[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->Clone();
1780   TH1D *c3_pp=(TH1D*)c3[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->Clone();
1781
1782   c3_pPb->GetYaxis()->SetTitle("#bf{c}_{3}^{#pm#pm#pm}");
1783   c3_pPb->GetXaxis()->SetRangeUser(0,0.27);
1784
1785   c3_pPb->GetXaxis()->SetLabelSize(SizeLabel); c3_pPb->GetYaxis()->SetLabelSize(SizeLabel);
1786   c3_pPb->GetXaxis()->SetTitleSize(SizeTitle); c3_pPb->GetYaxis()->SetTitleSize(SizeTitle);
1787
1788   c3_pPb->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}^{#pm#pm#pm}");
1789   c3_pPb->GetXaxis()->SetTitleOffset(1.0);
1790   c3_pPb->GetYaxis()->SetTitleOffset(0.7);
1791   c3_pPb->SetMinimum(0.9);
1792   c3_pPb->SetMaximum(3.3);// 3.3
1793   c3_pPb->GetXaxis()->SetNdivisions(503);
1794   c3_pPb->GetYaxis()->SetNdivisions(503);
1795   c3_pPb->SetMarkerStyle(25);
1796
1797   c3_pPb->Draw();
1798   //
1799   if(pp_pPb_Comp) c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->Draw("E2 same");
1800   c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->Draw("E2 same");
1801   if(!pp_pPb_Comp) c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("E2 same");
1802   if(pp_pPb_Comp) c3_pp->Draw("same");
1803   c3_pPb->Draw("same");
1804   if(!pp_pPb_Comp) c3_PbPb->Draw("same");
1805
1806   if(pp_pPb_Comp) {
1807     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
1808     legend6->AddEntry(c3_pp,"#splitline{pp #sqrt{s}=7 TeV}{#LT#font[12]{N}_{ch}#GT = 27 #pm 1}","p");// Mpp=15, Nch=27
1809   }else{
1810     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
1811     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
1812   }
1813
1814   //
1815
1816   if(!pp_pPb_Comp) c3_fit[0][0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same");
1817   c3_fit[1][0][KT3Bin][Mbin_SysComp_pPb]->Draw("c same");
1818   if(pp_pPb_Comp) c3_fit[2][0][KT3Bin][Mbin_SysComp_pp]->Draw("c same");
1819   //
1820   TF1 *GaussFit_forLegend=(TF1*)c3_fit[1][0][KT3Bin][Mbin_SysComp_pPb]->Clone();
1821   TF1 *EwFit_forLegend=(TF1*)c3_fit[1][1][KT3Bin][Mbin_SysComp_pPb]->Clone();
1822   TF1 *ExpFit_forLegend=(TF1*)c3_fit[1][2][KT3Bin][Mbin_SysComp_pPb]->Clone();
1823   GaussFit_forLegend->SetLineColor(1);
1824   EwFit_forLegend->SetLineColor(1);
1825   ExpFit_forLegend->SetLineColor(1);
1826   legend7->AddEntry(GaussFit_forLegend,"Gaussian","l");
1827   legend7->AddEntry(EwFit_forLegend,"Edgeworth","l");
1828   legend7->AddEntry(ExpFit_forLegend,"Exponential","l");
1829   //
1830
1831
1832   // spline draw
1833   if(!pp_pPb_Comp) gr_c3Spline[0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same");
1834   gr_c3Spline[1][KT3Bin][Mbin_SysComp_pPb]->Draw("c same");
1835   if(pp_pPb_Comp) gr_c3Spline[2][KT3Bin][Mbin_SysComp_pp]->Draw("c same");
1836   //
1837
1838   //exp draw
1839   if(!pp_pPb_Comp) gr_c3SplineExpFit[0][KT3Bin][Mbin_SysComp_PbPb]->Draw("c same");
1840   gr_c3SplineExpFit[1][KT3Bin][Mbin_SysComp_pPb]->Draw("c same");
1841   if(pp_pPb_Comp) gr_c3SplineExpFit[2][KT3Bin][Mbin_SysComp_pp]->Draw("c same");
1842   //
1843
1844
1845   legend6->Draw("same");
1846   legend7->Draw("same");
1847
1848   TLatex *Specif_Kt3_4 = new TLatex(0.52,0.55,"0.16<#font[12]{K}_{T,3}<0.3 GeV/#font[12]{c}");
1849   Specif_Kt3_4->SetNDC();
1850   Specif_Kt3_4->SetTextFont(TextFont);
1851   Specif_Kt3_4->SetTextSize(SizeSpecif);
1852   Specif_Kt3_4->Draw("same");
1853
1854   // print out data points
1855   for(int binN=1; binN<=c3_pPb->GetNbinsX(); binN++){
1856     if(pp_pPb_Comp) cout<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" TO "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<"; "<<c3_pp->GetBinContent(binN)<<" +- "<<c3_pp->GetBinError(binN)<<" (DSYS="<<c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->GetBinError(binN)<<"); "<<c3_pPb->GetBinContent(binN)<<" +- "<<c3_pPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<");"<<endl;
1857     else cout<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" TO "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<"; "<<c3_pPb->GetBinContent(binN)<<" +- "<<c3_pPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<"); "<<c3_PbPb->GetBinContent(binN)<<" +- "<<c3_PbPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->GetBinError(binN)<<");"<<endl;
1858   }
1859
1860   //
1861   if(SaveFiles) {
1862     TString *name = new TString("c3SystemComp");
1863     if(FitType==1) name->Append("EW");
1864     name->Append("_MpPb");
1865     *name += Mbin_SysComp_pPb;
1866     name->Append(".eps");
1867     //
1868     can5->SaveAs(name->Data());
1869   }
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884   ////////////////////////////////////////////////////
1885   ////////////////////////////////////////////////////
1886
1887   TCanvas *can6 = new TCanvas("can6", "can6",1700,700,600,600);// 11,53,700,500
1888   can6->SetHighLightColor(2);
1889   gStyle->SetOptFit(0);// 0111 to show fit stat box
1890   can6->SetFillColor(0);//10
1891   can6->SetBorderMode(0);
1892   can6->SetBorderSize(2);
1893   can6->SetFrameFillColor(0);
1894   can6->SetFrameBorderMode(0);
1895   can6->SetFrameBorderMode(0);
1896   can6->cd();
1897   TPad *pad6 = new TPad("pad6","pad6",0.0,0.0,1.,1.);
1898   gPad->SetGridx(0);
1899   gPad->SetGridy(0);
1900   pad6->SetTopMargin(0.0);//0.05
1901   pad6->SetRightMargin(0.0);//3e-2
1902   pad6->SetBottomMargin(0.0);//0.12
1903   pad6->SetLeftMargin(0.0);//0.12
1904   pad6->Draw();
1905   pad6->cd();
1906   TLegend *legend8 = new TLegend(.17,.4, .5,.6,NULL,"brNDC");//.45 or .4 for x1
1907   legend8->SetBorderSize(0);
1908   legend8->SetFillColor(0);
1909   legend8->SetTextFont(TextFont);
1910   legend8->SetTextSize(0.8*SizeLegend);
1911   TLegend *legend9 = new TLegend(.17,.6, .6,.98,NULL,"brNDC");//.45 or .4 for x1
1912   legend9->SetBorderSize(0);
1913   legend9->SetFillColor(0);
1914   legend9->SetTextFont(TextFont);
1915   legend9->SetTextSize(0.8*SizeLegend);
1916   //
1917   gPad->SetLeftMargin(0.14);
1918   gPad->SetRightMargin(0.01);
1919   gPad->SetTopMargin(0.01);
1920   gPad->SetBottomMargin(0.14);
1921
1922   gPad->SetTickx(); gPad->SetTicky();
1923   if(!NchOneThirdAxis) gPad->SetLogx();
1924   RadiiC2PbPb->GetXaxis()->SetLabelFont(TextFont); RadiiC2PbPb->GetYaxis()->SetLabelFont(TextFont);
1925   RadiiC2PbPb->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); RadiiC2PbPb->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel);
1926   RadiiC2PbPb->GetXaxis()->SetLabelOffset(-0.01);
1927   RadiiC2PbPb->GetXaxis()->SetNdivisions(808);
1928   RadiiC2PbPb->GetXaxis()->SetTitleOffset(1.05);//1.3
1929   RadiiC2PbPb->GetYaxis()->SetTitleOffset(1.1);//1.1
1930   RadiiC2PbPb->GetXaxis()->SetTitleFont(TextFont); RadiiC2PbPb->GetXaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel
1931   RadiiC2PbPb->GetYaxis()->SetTitleFont(TextFont); RadiiC2PbPb->GetYaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel
1932   RadiiC2PbPb->SetMinimum(0.01); RadiiC2PbPb->SetMaximum(11.9);// 0 and 11.9
1933   //gStyle->SetErrorX(0);
1934   if(NchOneThirdAxis) RadiiC2PbPb->GetXaxis()->SetRangeUser(0,3000);// 0,3000
1935   else RadiiC2PbPb->GetXaxis()->SetRangeUser(3,3000);// 3,3000
1936   //
1937   //
1938   Parameters_Bjoern[0][0]->GetXaxis()->SetLabelFont(TextFont); Parameters_Bjoern[0][0]->GetYaxis()->SetLabelFont(TextFont);
1939   Parameters_Bjoern[0][0]->GetXaxis()->SetLabelSize(SizeLabel*SF_2panel); Parameters_Bjoern[0][0]->GetYaxis()->SetLabelSize(SizeLabel*SF_2panel);
1940   Parameters_Bjoern[0][0]->GetXaxis()->SetLabelOffset(-0.01);
1941   Parameters_Bjoern[0][0]->GetXaxis()->SetNdivisions(808);
1942   Parameters_Bjoern[0][0]->GetXaxis()->SetTitleOffset(1.05);//1.3
1943   Parameters_Bjoern[0][0]->GetYaxis()->SetTitleOffset(1.1);//1.1
1944   Parameters_Bjoern[0][0]->GetXaxis()->SetTitleFont(TextFont); Parameters_Bjoern[0][0]->GetXaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel
1945   Parameters_Bjoern[0][0]->GetYaxis()->SetTitleFont(TextFont); Parameters_Bjoern[0][0]->GetYaxis()->SetTitleSize(SizeTitle);// SizeTitle*SF_2panel
1946   Parameters_Bjoern[0][0]->SetMinimum(0.01); Parameters_Bjoern[0][0]->SetMaximum(8.5);// 0 and 11.9
1947   //gStyle->SetErrorX(0);
1948   if(NchOneThirdAxis) Parameters_Bjoern[0][0]->GetXaxis()->SetRangeUser(0,3000);// 0,3000
1949   else Parameters_Bjoern[0][0]->GetXaxis()->SetRangeUser(3,3000);// 3,3000
1950
1951
1952   RadiiC2PbPb->GetXaxis()->SetTitle("#LT#font[12]{N}_{ch}#GT");
1953   RadiiC2PbPb->GetYaxis()->SetTitle("Radius (fm)");
1954   Parameters_Bjoern[0][0]->GetXaxis()->SetTitle("#LT#font[12]{N}_{ch}#GT");
1955   Parameters_Bjoern[0][0]->GetYaxis()->SetTitle("Radius (fm)");
1956
1957   RadiiC2PbPb->Draw();
1958
1959
1960
1961   //int HydroCase=0;// 0 or 1
1962   //Parameters_Bjoern[HydroCase][0]->Draw("same");
1963   //Parameters_Bjoern[HydroCase][1]->Draw("same");
1964   //Parameters_Bjoern[HydroCase][2]->Draw("same");
1965
1966   legend8->AddEntry(RadiiC2pp,"ALICE pp","p");
1967   legend8->AddEntry(RadiiC2pPb,"ALICE p-Pb","p");
1968   legend8->AddEntry(RadiiC2PbPb,"ALICE Pb-Pb","p");
1969   //
1970   //legend8->AddEntry(Parameters_Bjoern[HydroCase][2],"IP-GLASMA pp (w/o hydro)","p");
1971   //legend8->AddEntry(Parameters_Bjoern[HydroCase][1],"IP-GLASMA p-Pb (w/o hydro)","p");
1972   //legend8->AddEntry(Parameters_Bjoern[HydroCase][0],"IP-GLASMA Pb-Pb (w/o hydro)","p");
1973
1974
1975   Parameters_Bjoern[0][0]->SetMarkerStyle(24); Parameters_Bjoern[0][1]->SetMarkerStyle(25); Parameters_Bjoern[0][2]->SetMarkerStyle(28);
1976
1977   //TMultiGraph *mg = new TMultiGraph();
1978   //mg->SetTitle("");
1979   TGraph *Radii_Bjoern[2][3];// Hydro case, CollType
1980   TGraph *grShade[3];// CollType
1981
1982   Radii_Bjoern[0][0] = new TGraph(10, Bjoern_xaxis_PbPb, Bjoern_Ri_PbPb[mchoice-1]);
1983   Radii_Bjoern[0][0]->SetLineWidth(2);
1984
1985   Radii_Bjoern[1][0] = new TGraph(10, Bjoern_xaxis_PbPb, Bjoern_Rhydro_PbPb[mchoice-1]);
1986   Radii_Bjoern[1][0]->SetLineWidth(2);
1987   Radii_Bjoern[1][0]->SetLineStyle(2);
1988   Radii_Bjoern[1][0]->SetFillStyle(1000);
1989   Radii_Bjoern[1][0]->SetFillColor(kGray);
1990   grShade[0] = new TGraph(10);
1991   for(int i=0; i<10; i++){
1992     grShade[0]->SetPoint(i,Bjoern_xaxis_PbPb[i],Bjoern_Rhydro_PbPb[mchoice-1][i]);
1993     grShade[0]->SetPoint(10+i,Bjoern_xaxis_PbPb[10-i-1],Bjoern_Ri_PbPb[mchoice-1][10-i-1]);
1994   }
1995   grShade[0]->SetFillStyle(1000);
1996   grShade[0]->SetFillColor(kGray);
1997   //
1998   Radii_Bjoern[0][1] = new TGraph(6, Bjoern_xaxis_pPb, Bjoern_Ri_pPb[mchoice-1]);
1999   Radii_Bjoern[0][1]->SetLineWidth(2);
2000   Radii_Bjoern[0][1]->SetLineColor(2);
2001   Radii_Bjoern[1][1] = new TGraph(6, Bjoern_xaxis_pPb, Bjoern_Rhydro_pPb[mchoice-1]);
2002   Radii_Bjoern[1][1]->SetLineWidth(2);
2003   Radii_Bjoern[1][1]->SetLineColor(2);
2004   Radii_Bjoern[1][1]->SetLineStyle(2);
2005   grShade[1] = new TGraph(6);
2006   for(int i=0; i<6; i++){
2007     grShade[1]->SetPoint(i,Bjoern_xaxis_pPb[i],Bjoern_Rhydro_pPb[mchoice-1][i]);
2008     grShade[1]->SetPoint(6+i,Bjoern_xaxis_pPb[6-i-1],Bjoern_Ri_pPb[mchoice-1][6-i-1]);
2009   }
2010   grShade[1]->SetFillStyle(1000);
2011   grShade[1]->SetFillColor(kRed-10);
2012   //
2013   Radii_Bjoern[0][2] = new TGraph(4, Bjoern_xaxis_pp, Bjoern_Ri_pp[mchoice-1]);
2014   Radii_Bjoern[0][2]->SetLineWidth(2);
2015   Radii_Bjoern[0][2]->SetLineColor(4);
2016   Radii_Bjoern[1][2] = new TGraph(4, Bjoern_xaxis_pp, Bjoern_Rhydro_pp[mchoice-1]);
2017   Radii_Bjoern[1][2]->SetLineWidth(2);
2018   Radii_Bjoern[1][2]->SetLineColor(4);
2019   Radii_Bjoern[1][2]->SetLineStyle(2);
2020   grShade[2] = new TGraph(4);
2021   for(int i=0; i<4; i++){
2022     grShade[2]->SetPoint(i,Bjoern_xaxis_pp[i],Bjoern_Rhydro_pp[mchoice-1][i]);
2023     grShade[2]->SetPoint(4+i,Bjoern_xaxis_pp[4-i-1],Bjoern_Ri_pp[mchoice-1][4-i-1]);
2024   }
2025   grShade[2]->SetFillStyle(1000);
2026   grShade[2]->SetFillColor(kBlue-10);
2027   //
2028   grShade[0]->Draw("f same");
2029   grShade[1]->Draw("f same");
2030   grShade[2]->Draw("f same");
2031   Radii_Bjoern[0][0]->Draw("l same");
2032   Radii_Bjoern[0][1]->Draw("l same");
2033   Radii_Bjoern[0][2]->Draw("l same");
2034
2035   grRadiiC2Sys_pp->Draw("|| p");
2036   grRadiiC2Sys_pPb->Draw("|| p");
2037   grRadiiC2Sys_PbPb->Draw("E p");
2038   RadiiC2PbPb->Draw("same");
2039   RadiiC2pPb->Draw("same");
2040   RadiiC2pp->Draw("same");
2041
2042   legend9->AddEntry(Radii_Bjoern[0][2],"GLASMA pp R_{initial}","l");
2043   legend9->AddEntry(Radii_Bjoern[0][1],"GLASMA p-Pb R_{initial}","l");
2044   legend9->AddEntry(Radii_Bjoern[0][0],"GLASMA Pb-Pb R_{initial}","l");
2045   legend9->AddEntry(grShade[2],"GLASMA pp R_{hydro}","f");
2046   legend9->AddEntry(grShade[1],"GLASMA p-Pb R_{hydro}","f");
2047   legend9->AddEntry(grShade[0],"GLASMA Pb-Pb R_{hydro}","f");
2048   /*
2049   TF1 *fit1=new TF1("fit1","pol5",0,1000);
2050   fit1->SetLineColor(1);
2051   Radii_Bjoern[0][0]->Fit(fit1,"Q","",30,800);
2052   //fit1->Draw("same");
2053   //
2054   TF1 *fit2=new TF1("fit2","pol4",0,1000);
2055   Radii_Bjoern[0][1]->Fit(fit2,"Q","",4,70);
2056   //fit2->Draw("same");
2057   //
2058   TF1 *fit3=new TF1("fit3","pol4",0,1000);
2059   fit3->SetLineColor(1);
2060   Radii_Bjoern[0][2]->Fit(fit3,"Q","",4,27);
2061   //fit3->Draw("same");
2062
2063   TH1D *BjoernRatio_PbPb = (TH1D*)RadiiC2PbPb->Clone();
2064   TH1D *BjoernRatio_pPb = (TH1D*)RadiiC2pPb->Clone();
2065   TH1D *BjoernRatio_pp = (TH1D*)RadiiC2pp->Clone();
2066   for(int point=0; point<20; point++){
2067     int bin = BjoernRatio_PbPb->GetXaxis()->FindBin(meanNchPbPb[point]);
2068     if(meanNchPbPb[point] > Bjoern_xaxis_PbPb[8]) {BjoernRatio_PbPb->SetBinContent(bin,0); continue;}
2069     if(BjoernRatio_PbPb->GetBinContent(bin) < 0.5 || BjoernRatio_PbPb->GetBinContent(bin) > 12) BjoernRatio_PbPb->SetBinContent(bin,0);
2070     else{
2071       //      double xx =
2072       BjoernRatio_PbPb->SetBinContent(bin, BjoernRatio_PbPb->GetBinContent(bin) / fit1->Eval(meanNchPbPb[point]));
2073     }
2074   }
2075   for(int point=0; point<20; point++){
2076     int bin = BjoernRatio_pPb->GetXaxis()->FindBin(meanNchpPb[point]);
2077     if(meanNchpPb[point] > Bjoern_xaxis_pPb[5]) {BjoernRatio_pPb->SetBinContent(bin,0); continue;}
2078     if(BjoernRatio_pPb->GetBinContent(bin) < 0.5 || BjoernRatio_pPb->GetBinContent(bin) > 12) BjoernRatio_pPb->SetBinContent(bin,0);
2079     else{
2080       BjoernRatio_pPb->SetBinContent(bin, BjoernRatio_pPb->GetBinContent(bin) / fit2->Eval(meanNchpPb[point]));
2081     }
2082   }
2083   for(int point=0; point<20; point++){
2084     int bin = BjoernRatio_pp->GetXaxis()->FindBin(meanNchpp[point]);
2085     if(meanNchpp[point] > Bjoern_xaxis_pp[3]) {BjoernRatio_pp->SetBinContent(bin,0); continue;}
2086     if(BjoernRatio_pp->GetBinContent(bin) < 0.5 || BjoernRatio_pp->GetBinContent(bin) > 12) BjoernRatio_pp->SetBinContent(bin,0);
2087     else{
2088       BjoernRatio_pp->SetBinContent(bin, BjoernRatio_pp->GetBinContent(bin) / fit2->Eval(meanNchpp[point]));
2089     }
2090   }
2091   BjoernRatio_PbPb->SetMinimum(0.55); BjoernRatio_PbPb->SetMaximum(1.45);
2092   BjoernRatio_PbPb->GetYaxis()->SetTitle("Radius Ratio (ALICE/GLASMA)");
2093   BjoernRatio_PbPb->Draw();
2094   BjoernRatio_pPb->Draw("same");
2095   BjoernRatio_pp->Draw("same");
2096   */
2097
2098   /*legend8->AddEntry(Parameters_Bjoern[0][2],"pp R_{initial} (no hydro)","p");
2099   legend8->AddEntry(Parameters_Bjoern[0][1],"p-Pb R_{initial} (no hydro)","p");
2100   legend8->AddEntry(Parameters_Bjoern[0][0],"Pb-Pb R_{initial} (no hydro)","p");
2101   legend9->AddEntry(Parameters_Bjoern[1][2],"pp R_{max} (hydro)","p");
2102   legend9->AddEntry(Parameters_Bjoern[1][1],"p-Pb R_{max} (hydro)","p");
2103   legend9->AddEntry(Parameters_Bjoern[1][0],"Pb-Pb R_{max} (hydro)","p");
2104   */
2105   legend8->Draw("same");
2106   legend9->Draw("same");
2107   // Normalization plots
2108   // ct, ft, KT3, par
2109   //Parameters_C2[2][0][0][0]->Draw();
2110
2111
2112
2113
2114 }
2115 //____________________________________________________________________________________________________
2116 void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2)
2117 {
2118   // revision on July 24th or 25th 2012
2119   // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
2120   x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
2121   // Printf("%f %f %f %f", x1, x2, y1, y2);
2122
2123   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
2124   myPadLogo->SetLeftMargin(0);
2125   myPadLogo->SetTopMargin(0);
2126   myPadLogo->SetRightMargin(0);
2127   myPadLogo->SetBottomMargin(0);
2128   myPadLogo->Draw();
2129   myPadLogo->cd();
2130   TASImage *myAliceLogo = new TASImage((prel) ? "~/Pictures/2011-Nov-24-ALICE_PRELIMARY_logo_BLACK_small_usage_design.eps" : "~/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.eps");
2131   myAliceLogo->Draw();
2132 }
2133 //____________________________________________________________________________________________________
2134 TCanvas *make_canvas(const Char_t *name, const Char_t *title, Int_t n_x, Int_t n_y, Int_t share_axes, Int_t ww, Int_t wh){
2135   const Int_t width=350;
2136   const Int_t height=350;
2137   TCanvas *canvas;
2138
2139   if (share_axes==0) {
2140     if (ww==0)
2141       ww=width*n_x;
2142     if (wh==0)
2143       wh=height*n_y;
2144     canvas=new TCanvas(name,title,10,0,ww,wh);
2145     canvas->SetTopMargin(0.01);
2146     canvas->SetRightMargin(0.01);
2147     canvas->SetLeftMargin(0.14);
2148     canvas->SetBottomMargin(0.18);
2149
2150     canvas->Divide(n_x,n_y,0,0);
2151   }
2152   else {
2153     Float_t pix_width=(1-gStyle->GetPadLeftMargin()-gStyle->GetPadRightMargin())*width;
2154     Float_t pix_height=(1-gStyle->GetPadTopMargin()-gStyle->GetPadBottomMargin())*height;
2155     if (ww==0)
2156       ww=width+(n_x-1)*pix_width;
2157     if (wh==0)
2158       wh=height+(n_y-1)*pix_height;
2159
2160     canvas=new TCanvas(name,title,ww,wh);
2161     canvas->SetTopMargin(0); canvas->SetBottomMargin(0);
2162     canvas->SetLeftMargin(0); canvas->SetRightMargin(0);
2163     Float_t tot_width;
2164     if (n_x>1)
2165       tot_width=(n_x-2)+1./(1-gStyle->GetPadLeftMargin())
2166         +1./(1-gStyle->GetPadRightMargin());
2167     else
2168       tot_width=1./(1-gStyle->GetPadLeftMargin()-gStyle->GetPadRightMargin());
2169     Float_t tot_height;
2170     if (n_y>1)
2171       tot_height=(n_y-2)+1./(1-gStyle->GetPadTopMargin())
2172         +1./(1-gStyle->GetPadBottomMargin());
2173     else
2174       tot_height=1./(1-gStyle->GetPadTopMargin()-gStyle->GetPadBottomMargin());
2175
2176     //Int_t idx=n_x*n_y;
2177
2178     for (Int_t j=0; j<n_y; j++) {
2179       for (Int_t i=0; i<n_x; i++) {
2180         //for (Int_t j=n_y-1; j>=0; j--) {
2181         //for (Int_t i=n_x-1; i>=0; i--) {
2182         Char_t tmp_str[256];
2183         Char_t p_title[256];
2184         Int_t idx=n_x*j+i+1;
2185         sprintf(tmp_str,"%s_%d",canvas->GetName(),idx);
2186         sprintf(p_title,"Pad %d",idx);
2187         Float_t x1=0,y1=0;
2188         Float_t x2=1,y2=1;
2189         if (n_x>1) {
2190           if (i==0)
2191             x2=1./(1-gStyle->GetPadLeftMargin())/tot_width;
2192           else {
2193             x1=(1./(1-gStyle->GetPadLeftMargin())+i-1)/tot_width;
2194             if (i<n_x-1)
2195               x2=(1./(1-gStyle->GetPadLeftMargin())+i)/tot_width;
2196           }
2197         }
2198         if (n_y>1) {
2199           if (j==0)
2200             y1=1-1./(1-gStyle->GetPadTopMargin())/tot_height;
2201           else {
2202             y2=1-(1./(1-gStyle->GetPadTopMargin())+j-1)/tot_height;
2203             if (j<n_y-1)
2204               y1=1-(1./(1-gStyle->GetPadTopMargin())+j)/tot_height;
2205           }
2206         }
2207         //cout << "x1 " << x1 << ", x2 " << x2 << endl;
2208         TPad *pad=new TPad(tmp_str,title,x1,y1,x2,y2);
2209         //pad->SetFillColor(idx+1);
2210         if (i>0)
2211           pad->SetLeftMargin(0.001);
2212         if (i<n_x-1)
2213           pad->SetRightMargin(0.001);
2214         if (j>0)
2215           pad->SetTopMargin(0.001);
2216         if (j<n_y-1)
2217           pad->SetBottomMargin(0.001);
2218
2219
2220         pad->SetNumber(idx);
2221         //pad->SetNumber(n_x*j+i+1);
2222         pad->Draw();
2223         //idx--;
2224         //idx++;
2225         //pad->SetP
2226       }
2227     }
2228   }
2229
2230   return canvas;
2231 }