Macro updates
authordgangadh <dhevan.raja.gangadharan@cern.ch>
Sun, 28 Sep 2014 14:13:24 +0000 (16:13 +0200)
committerdgangadh <dhevan.raja.gangadharan@cern.ch>
Sun, 28 Sep 2014 14:13:24 +0000 (16:13 +0200)
PWGCF/FEMTOSCOPY/macros/Fit_c3.C [new file with mode: 0755]
PWGCF/FEMTOSCOPY/macros/Makec3EAfile.C [new file with mode: 0755]
PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C
PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C
PWGCF/FEMTOSCOPY/macros/Plot_plotsTPR.C

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