Updates for many plotting macros
authordgangadh <dhevan.raja.gangadharan@cern.ch>
Mon, 17 Nov 2014 09:49:40 +0000 (10:49 +0100)
committerdgangadh <dhevan.raja.gangadharan@cern.ch>
Mon, 17 Nov 2014 09:49:40 +0000 (10:49 +0100)
PWGCF/FEMTOSCOPY/macros/Fit_c3.C
PWGCF/FEMTOSCOPY/macros/Makec3EAfile.C
PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C
PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C
PWGCF/FEMTOSCOPY/macros/Plot_plotsTPR.C

index 04d65a0..fd15e96 100755 (executable)
 using namespace std;
 
 //
-int CollisionType_def=1;// PbPb, pPb, pp
+int CollisionType_def=0;// 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 %
+double G_def = 0.;// coherent %
+double Rcoh_def = 1.;// Radius of Gaussian coherent component in fm
 //
 bool MRC=1;// Momentum Resolution Corrections?
 bool MuonCorrection=1;// correct for Muons?
@@ -107,32 +108,37 @@ double BinCentersQ4[20];
 
 
 
-void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, double G=G_def){
+void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, double G=G_def, double Rcoh=Rcoh_def){
   
   SaveToFile_def=SaveToFile;
   CollisionType_def=CollisionType;
   G /= 100.;
   G_def=G;
+  Rcoh_def=Rcoh;
 
   TFile *_file0;
   if(CollisionType==0){// PbPb
-    _file0 = new TFile("Results/PDC_11h_3Dhistos.root","READ");
+    //_file0 = new TFile("Results/PDC_11h_3Dhistos.root","READ");
+    //_file0 = new TFile("Results/PDC_11h_c3FitBuild.root","READ");
+    _file0 = new TFile("Results/PDC_11h_LowNorm_HighNorm.root","READ");
   }else if(CollisionType==1){// pPb
-    _file0 = new TFile("Results/PDC_13bc_kINT7.root","READ");
+    //_file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm.root","READ");
+    _file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm_HighNorm.root","READ");
   }else{// pp
-    _file0 = new TFile("Results/PDC_10bcde_kMB.root","READ");
+    _file0 = new TFile("Results/PDC_10bcde_kMB_3Dhisto_LowNorm_HighNorm.root","READ");
   }
 
   TList *MyList;
   TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root");
-  MyList=(TList*)tdir->Get("FourPionOutput_1");
+  if(CollisionType==0) MyList=(TList*)tdir->Get("FourPionOutput_1");
+  else MyList=(TList*)tdir->Get("FourPionOutput_2");
   //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
+  if(CollisionType==0) {Q3LimitLow = 0.01; Q3LimitHigh = 0.08;}// 0.01 and 0.08 
+  else {Q3LimitLow = 0.01; Q3LimitHigh = 0.25;}// 0.01 and 0.25
   
   
   //
@@ -487,11 +493,11 @@ void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def,
   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;
+  //          1.0              1.0              10.              0.              0.              0.             1.5  
+  par_c3[0] = 1.0; par_c3[1] = 1.0; 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;
+  minVal_c3[0] = 0.8; minVal_c3[1] = 0.2; 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] = 1000.; 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){ 
@@ -543,7 +549,7 @@ void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def,
     MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]);
   }
   
-  cout<<"Tij Norm = "<<pow(OutputPar_c3[1]/2., 1/3.)<<endl;
+  cout<<"Tij Norm = "<<pow(OutputPar_c3[1], 1/3.)<<endl;
   
   if(FitType!=2){
     c3Fit1D_Expan->FixParameter(0,OutputPar_c3[0]);
@@ -613,12 +619,14 @@ void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def,
        }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;
+       double t12_coh = exp(-pow(Rcoh/FmToGeV * Q12,2)/2.);
+       double t23_coh = exp(-pow(Rcoh/FmToGeV * Q23,2)/2.);
+       double t13_coh = exp(-pow(Rcoh/FmToGeV * Q13,2)/2.);
+       double C = 2*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 += 2*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*t23_coh;
+       C += 2*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*t13_coh;
+       C += 2*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*t12_coh;
        C += 1.0;
-       //if(Q3<0.026) cout<<Q3<<"  "<<C<<endl;
        C *= TERM5;
        C *= OutputPar_c3[0];
        //if(Q3<0.018) continue;
@@ -689,6 +697,8 @@ void Fit_c3(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def,
     *savefilename += CollisionType;
     savefilename->Append("_FT");
     *savefilename += FitType;
+    savefilename->Append("_R");
+    *savefilename += Rcoh;
     savefilename->Append("_G");
     *savefilename += int((G+0.001)/0.02);
     savefilename->Append(".root");
@@ -912,11 +922,13 @@ void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){
          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;
+       double t12_coh = exp(-pow(Rcoh_def/FmToGeV * q12,2)/2.);
+       double t23_coh = exp(-pow(Rcoh_def/FmToGeV * q23,2)/2.);
+       double t13_coh = exp(-pow(Rcoh_def/FmToGeV * q13,2)/2.);
+       C = 2*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 += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg13,par[6]))/2.)*Expan12*Expan13*t23_coh;
+       C += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg12,par[6])+pow(arg23,par[6]))/2.)*Expan12*Expan23*t13_coh;
+       C += 2*pow(par[1],2/3.) * G_def*pow(1-G_def,2)*exp(-(pow(arg13,par[6])+pow(arg23,par[6]))/2.)*Expan13*Expan23*t12_coh;
        C += 1.0;
        C *= par[0];// Norm
        
index 6078a27..29525f0 100755 (executable)
@@ -34,17 +34,18 @@ using namespace std;
 
 void Makec3EAfile(){
 
+  int FT=0;// 0(EW) or 1(LG)
+
   TFile *infile;
   TMinuit *fit;
   //
-  TH3D *PbPbEA = new TH3D("PbPbEA","",2,0.5,2.5, 4,0.5,4.5, 50,-0.5,49.5);
+  TH3D *PbPbEA = new TH3D("PbPbEA","",6,0.5,6.5, 4,0.5,4.5, 50,-0.5,49.5);// Rcoh type, parNum, Gindex
   PbPbEA->SetDirectory(0);
   //
-  TH3D *pPbEA = new TH3D("pPbEA","",2,0.5,2.5, 4,0.5,4.5, 50,-0.5,49.5);
+  TH3D *pPbEA = new TH3D("pPbEA","",6,0.5,6.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);
+  TH3D *ppEA = new TH3D("ppEA","",6,0.5,6.5, 4,0.5,4.5, 50,-0.5,49.5);
   ppEA->SetDirectory(0);
   //
  
@@ -55,12 +56,14 @@ void Makec3EAfile(){
   //////////////////////////////
   double value=0, value_e=0;
   // 
-  for(int Gindex=0; Gindex<46; Gindex++){
+  for(int Gindex=0; Gindex<=25; Gindex++){
     // PbPb
-    for(int FT=0; FT<2; FT++){// EW or LG
+    for(int RT=0; RT<6; RT++){// Rcoh type
       
       TString *name1 = new TString("FitFiles/FitFile_CT0_FT");
       *name1 += FT;
+      name1->Append("_R");
+      *name1 += RT;
       name1->Append("_G");
       *name1 += Gindex;
       name1->Append(".root");
@@ -68,7 +71,7 @@ void Makec3EAfile(){
       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);
+       PbPbEA->SetBinContent(RT+1, parNum+1, Gindex+1, value);
       }
       infile->Close();
    
@@ -76,6 +79,8 @@ void Makec3EAfile(){
       // pPb
       TString *name1 = new TString("FitFiles/FitFile_CT1_FT");
       *name1 += FT;
+      name1->Append("_R");
+      *name1 += RT;
       name1->Append("_G");
       *name1 += Gindex;
       name1->Append(".root");
@@ -83,13 +88,15 @@ void Makec3EAfile(){
       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);
+       pPbEA->SetBinContent(RT+1, parNum+1, Gindex+1, value);
       }
       infile->Close();
       //
       //
       TString *name1 = new TString("FitFiles/FitFile_CT2_FT");
       *name1 += FT;
+      name1->Append("_R");
+      *name1 += RT;
       name1->Append("_G");
       *name1 += Gindex;
       name1->Append(".root");
@@ -97,28 +104,28 @@ void Makec3EAfile(){
       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);
+       ppEA->SetBinContent(RT+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 Gindex=26; Gindex<50; Gindex++){
+    for(int RT=0; RT<6; RT++){// 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));
+       PbPbEA->SetBinContent(RT+1, parNum+1, Gindex+1, PbPbEA->GetBinContent(RT+1, parNum+1, 26));
+       pPbEA->SetBinContent(RT+1, parNum+1, Gindex+1, pPbEA->GetBinContent(RT+1, parNum+1, 26));
+       ppEA->SetBinContent(RT+1, parNum+1, Gindex+1, ppEA->GetBinContent(RT+1, parNum+1, 26));
       }
     }
   }
 
   // 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.));
+    for(int RT=0; RT<6; RT++){// EW or LG
+      PbPbEA->SetBinContent(RT+1, 1, Gindex+1, pow(PbPbEA->GetBinContent(RT+1, 1, Gindex+1), 1/3.));
+      pPbEA->SetBinContent(RT+1, 1, Gindex+1, pow(pPbEA->GetBinContent(RT+1, 1, Gindex+1), 1/3.));
+      ppEA->SetBinContent(RT+1, 1, Gindex+1, pow(ppEA->GetBinContent(RT+1, 1, Gindex+1), 1/3.));
     }
   }
   
index 6679b8c..011eb4f 100755 (executable)
@@ -38,17 +38,17 @@ using namespace std;
 int FileSetting=0;
 //
 bool MCcase_def=0;// MC data?
-int CollisionType=1;// PbPb, pPb, pp
+int CollisionType=0;// PbPb, pPb, pp
 //
 int Mbin_def=0;// 0-9: centrality bin in widths of 5%
 bool SameCharge_def=1;// same-charge?
-int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or -++,  ---+ or -+++
+int CHARGE_def=1;// -1 or +1: + or - pions for same-charge case, --+ or -++,  ---+ or -+++
 int MixedCharge4pionType_def = 2;// 1(---+) or 2(--++)
 //
 int EDbin_def=0;// 0 or 1: Kt3 bin
 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 Gbin = int( (0) /2. ) + 5;// +5 (Rcoh=0), +55 (Rcoh=Rch)
+int c3FitType = 1;// 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?
@@ -251,7 +251,10 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   
   TH2D *TwoParticle_2d[2][2][2];// ch1,ch2,term
   TH1D *TwoParticle[2][2][2];// ch1,ch2,term
+  TH2D *UnitMult_2d[2][2][2];// ch1,ch2,term
+  TH1D *UnitMult[2][2][2];// ch1,ch2,term
   double norm_2[2]={0};
+  double norm_2_UM[2]={0};
   //
   TH1D *ThreeParticle[2][2][2][5];// ch1,ch2,ch3,term
   TProfile *K3avg[2][2][2][4];
@@ -295,7 +298,8 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
       //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_c3FitBuild.root","READ");
+      if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_LowNorm_HighNorm.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");
@@ -305,9 +309,9 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
       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");
+    _file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm_HighNorm.root","READ");
   }else{// pp
-    _file0 = new TFile("Results/PDC_10bcde_kMB.root","READ");
+    _file0 = new TFile("Results/PDC_10bcde_kMB_LowNorm_HighNorm.root","READ");
   }
   
   SetFSIindex(10.);
@@ -332,7 +336,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   TList *MyList;
   if(!MCcase){
     TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root");
-    if(FileSetting==2 || FileSetting==5 || FileSetting==8) MyList=(TList*)tdir->Get("FourPionOutput_2");
+    if(FileSetting==2 || FileSetting==5 || FileSetting==8 || CollisionType!=0) MyList=(TList*)tdir->Get("FourPionOutput_2");
     else MyList=(TList*)tdir->Get("FourPionOutput_1");
     //MyList=(TList*)_file0->Get("MyList");
   }else{
@@ -385,7 +389,22 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
        TwoParticle[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
        TwoParticle[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
        TwoParticle[c1][c2][term]->SetTitle("");
+       //
        
+               TString *nameUM=new TString(name2->Data());
+       nameUM->Append("_UnitMult");
+       UnitMult_2d[c1][c2][term] = (TH2D*)MyList->FindObject(nameUM->Data());
+       UnitMult_2d[c1][c2][term]->Sumw2();
+       TString *pronameUM = new TString(nameUM->Data());
+       pronameUM->Append("_pro");
+       UnitMult[c1][c2][term] = (TH1D*)UnitMult_2d[c1][c2][term]->ProjectionY(pronameUM->Data(),13,13);// 11 means 1000 pions
+       norm_2_UM[term] = UnitMult[c1][c2][term]->Integral(UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.0),UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.2));
+       UnitMult[c1][c2][term]->Scale(norm_2_UM[0]/norm_2_UM[term]);
+       //
+       UnitMult[c1][c2][term]->SetMarkerStyle(20);
+       UnitMult[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
+       UnitMult[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
+       UnitMult[c1][c2][term]->SetTitle("");
       }// term
       
     
@@ -544,8 +563,10 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
            if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
            /////////////////////////////////////////
            norm_4[term] = ((TH1D*)MyList->FindObject(nameNorm4->Data()))->Integral();
+           
            //if( (c1+c2+c3+c4)==4) cout<<"4-pion norms  "<<norm_4[term]<<endl;
            if(norm_4[term] > 0){
+             //if(c1==c2 && c1==c3 && c1==c4) cout<<term<<"  "<<norm_4[0]/norm_4[term]<<endl;
              FourParticle[c1][c2][c3][c4][term] = (TH1D*)MyList->FindObject(name4->Data());
              FourParticle[c1][c2][c3][c4][term]->Sumw2();
              FourParticle[c1][c2][c3][c4][term]->Scale(norm_4[0]/norm_4[term]);
@@ -570,7 +591,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
              }
            }
            if(term==12 && c1==c2 && c1==c3 && c1==c4){
-            
+             
              TPFullWeight_FourParticle_2D[c1] = (TH2D*)MyList->FindObject(nameTPN4->Data());
              TPFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]);
              TPFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin);
@@ -579,7 +600,7 @@ 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){
+             if(c1==0 && !MCcase){
                FullBuildFromFits_3D = (TH3D*)MyList->FindObject(nameFitBuild->Data());
                FullBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]);
                FullBuildFromFits_3D->RebinZ(FourParticleRebin);
@@ -633,13 +654,13 @@ 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){
+             if(c1==0 && !MCcase){
                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);
+               FullBuildFromFits->SetLineColor(4);
                //
                PartialBuildFromFits = (TH1D*)PartialBuildFromFits_3D->ProjectionZ(proNamePartialFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin);
                PartialBuildFromFits->Add(tempDen2);
@@ -732,7 +753,8 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   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);
@@ -748,6 +770,25 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   C2QS->SetMarkerColor(2); C2QS->SetLineColor(2);
   if(!MCcase) C2QS->Draw("same");
 
+
+  // To visualize the Qcut and Norm Regions
+  //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
+  //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);  
+  //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);  
+  //TERM1_2->Scale(1/TERM1_2->Integral());
+  //for(int bin=1; bin<=16; bin++) QcutRegion->SetBinContent(bin,TERM1_2->GetBinContent(bin));
+  //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
+  //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
+  //QcutRegion->SetFillColor(4);
+  //NormRegion1->SetFillColor(2);
+  //NormRegion2->SetFillColor(3);
+  //TERM1_2->GetYaxis()->SetTitle("Pair probability");
+  //TERM1_2->GetYaxis()->SetTitleOffset(1.3);
+  //TERM1_2->GetXaxis()->SetRangeUser(0,1.6);
+  //TERM1_2->Draw();
+  //QcutRegion->Draw("same");
+  //cout<<TERM1_2->Integral(1,16) / TERM1_2->Integral()<<endl;
+
   // Print 2-pion values
   /*for(int bin=1; bin<=20; bin++){
     cout<<C2QS->GetBinContent(bin)<<", ";
@@ -775,7 +816,26 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   //C2ref_0p75->Draw("same");
 
 
-
+  TH1D *UnitMultC2=(TH1D*)UnitMult[0][0][0]->Clone();
+  TH1D *UnitMultC2den=(TH1D*)UnitMult[0][0][1]->Clone();
+  UnitMultC2->Divide(UnitMultC2den);
+  UnitMultC2->GetXaxis()->SetRangeUser(0,0.3);
+  UnitMultC2->SetMinimum(0.97);
+  UnitMultC2->SetMaximum(1.25);
+  UnitMultC2->Draw();
+  
+  TH1D *UnitMultC2QS=(TH1D*)UnitMult[0][0][0]->Clone();
+  UnitMultC2QS->Add(UnitMult[0][0][1], -(1-TwoFrac));
+  UnitMultC2QS->Scale(1/TwoFrac);
+  for(int bin=1; bin<=UnitMultC2QS->GetNbinsX(); bin++) {
+    Float_t K2 = 1.0;
+    if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, UnitMultC2QS->GetXaxis()->GetBinCenter(bin));
+    UnitMultC2QS->SetBinContent(bin, UnitMultC2QS->GetBinContent(bin)/K2);
+    UnitMultC2QS->SetBinError(bin, UnitMultC2QS->GetBinError(bin)/K2);
+  }
+  UnitMultC2QS->Divide(UnitMultC2den);
+  UnitMultC2QS->SetMarkerColor(2); UnitMultC2QS->SetLineColor(2);
+  UnitMultC2QS->Draw("same");
 
   ///////////////////////////////////////////////////////////
   // 3-pion 
@@ -803,14 +863,13 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   TH1D *TERM4_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][3]->Clone();
   TH1D *TERM5_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][4]->Clone();
   //
-  
   if(SameCharge){
     TERM1_3->Multiply(MRC_SC_3[0]);
     TERM2_3->Multiply(MRC_SC_3[1]);
     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]);
@@ -826,7 +885,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
       TERM5_3->Multiply(MRC_MC_3[3]);
     }
   }
+  
   
   TH1D *N3QS = (TH1D*)TERM1_3->Clone();
   N3QS->Add(TERM5_3, -f31);
@@ -842,6 +901,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   if(SameCharge) C3QS->Multiply(C3muonCorrectionSC[0]);
   else C3QS->Multiply(C3muonCorrectionMC[0]);
   
+
   C3QS->SetMarkerStyle(20);
   C3QS->SetMarkerColor(4);
   C3QS->GetXaxis()->SetRangeUser(0,0.1);
@@ -972,7 +1032,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
       //if(NDF>0) cout<<binG<<"  "<<sqrt(Chi2)/NDF<<endl;
     }
     Chi2NDF_PointSize_3->SetMarkerStyle(20); Chi2NDF_FullSize_3->SetMarkerStyle(20);
-    
+   
     C3QSratio->Divide(TPFullWeight_ThreeParticle[ch1_3]);
     //
     C3QSratio->SetMinimum(0.94); C3QSratio->SetMaximum(1.02);
@@ -1037,7 +1097,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   can3->SetFrameBorderMode(0);
   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
   TLegend *legend3 = (TLegend*)legend1->Clone();
-  legend3->SetX1(0.45); legend3->SetX2(0.98);  legend3->SetY1(0.6); legend3->SetY2(0.95);
+  legend3->SetX1(0.47); legend3->SetX2(0.98);  legend3->SetY1(0.7); legend3->SetY2(0.98);
 
   TH1D *TERMS_4[13];
   for(int index=0; index<13; index++){
@@ -1356,11 +1416,13 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   
 
 
-  if(SameCharge) C4QS->SetMaximum(9);
-  else if(MixedCharge4pionType==1) C4QS->SetMaximum(4);
+  if(SameCharge) {
+    if(CollisionType==0) C4QS->SetMaximum(7);
+    else C4QS->SetMaximum(12);
+  }else if(MixedCharge4pionType==1) C4QS->SetMaximum(4);
   else C4QS->SetMaximum(3);
   
-  if(CollisionType!=0) C4QS->GetXaxis()->SetRangeUser(0,0.6);
+  if(CollisionType!=0) C4QS->GetXaxis()->SetRangeUser(0,0.5);
 
   C4QS->Draw();
   c4QS_RemovalStage1->Draw("same");
@@ -1371,7 +1433,7 @@ void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool
   
   //cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(9)<<endl;
 
-  if(SameCharge) {
+  if(SameCharge && !MCcase) {
     //TPFullWeight_FourParticle[ch1_4]->Draw("same");
     FullBuildFromFits->Draw("same");
     PartialBuildFromFits->Draw("same");
index c70b33d..b0dbac0 100755 (executable)
@@ -38,18 +38,18 @@ double kappa4 = 0.5; // 0.5 (default),
 
 using namespace std;
 
-int CollisionType_def=2;// 0=PbPb, 1=pPb, 2=pp
+int CollisionType_def=0;// 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=0;// Gaussian or Exponential?
+bool Gaussian=1;// 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=16;// 0-19 (0=1050-2000 pions, 19=0-5 pions)
+int Mbin_def=2;// 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
index 76b3e4d..6eb2bea 100755 (executable)
@@ -38,7 +38,7 @@ using namespace std;
 
 bool SaveFiles=kFALSE;
 const int KT3Bin=0;// Kt3 bin. 0-1
-int FitType=2;// 0 (Gaussian), 1 (Edgeworth), 2 (Exponential)
+int FitType=1;// 0 (Gaussian), 1 (Edgeworth), 2 (Exponential)
 bool pp_pPb_Comp=0;
 bool AddedCC=kTRUE;// Charge Conjugate already added?
 bool NchOneThirdAxis=0;
@@ -839,7 +839,7 @@ void Plot_plotsTPR(){
   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)<<" TO "<<K0s_C3->GetXaxis()->GetBinUpEdge(i+1)<<"; "<<K0s_C3->GetBinContent(i+1)<<" +- "<<K0s_C3->GetBinError(i+1)<<"  (DSYS="<<Sys_K0s_C3->GetBinError(i+1)<<"); "<<K0s_c3->GetBinContent(i+1)<<" +- "<<K0s_c3->GetBinError(i+1)<<"  (DSYS="<<Sys_K0s_c3->GetBinError(i+1)<<"); "<<endl;
     //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;
   }
 
@@ -1303,36 +1303,50 @@ void Plot_plotsTPR(){
 
   // print radii and lambda
   cout.precision(3);
-  cout<<"Pb--Pb:"<<endl;
-  for(int cb=0; cb<20; cb++){
-    int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
-    if(RadiiPbPb->GetBinContent(binPbPb)==0) continue;
-    //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<<"Radii and Lambda data"<<endl;
+  cout<<"p--p:"<<endl;
+  // new way for HEP
+  for(int cb=19; cb>=0; cb--){
+    int binpp = RadiiC2pp->GetXaxis()->FindBin(meanNchpp[cb]);
+    if(Radiipp->GetBinContent(binpp)==0) continue;
+    cout<<meanNchpp[cb]-(0.05)*meanNchpp[cb]<<" TO "<<meanNchpp[cb]+(0.05)*meanNchpp[cb]<<"; "<<RadiiC2pp->GetBinContent(binpp)<<" +- "<<RadiiC2pp->GetBinError(binpp)<<" (DSYS=+"<<grRadiiC2Sys_pp->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_pp->GetErrorYlow(cb)<<"); ";
+    cout<<LambdaC2pp->GetBinContent(binpp)<<" +- "<<LambdaC2pp->GetBinError(binpp)<<" (DSYS=+"<<grLambdaC2Sys_pp->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_pp->GetErrorYlow(cb)<<"); ";
+    cout<<Radiipp->GetBinContent(binpp)<<" +- "<<Radiipp->GetBinError(binpp)<<" (DSYS="<<grRadiiSys_pp->GetErrorY(cb)<<"); ";
+    cout<<Lambdapp->GetBinContent(binpp)<<" +- "<<Lambdapp->GetBinError(binpp)<<" (DSYS="<<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;
   }
+  cout<<endl;
+  
   cout<<"p--Pb:"<<endl;
-  for(int cb=0; cb<20; cb++){
+  for(int cb=19; cb>=0; cb--){
     int binpPb = RadiiC2pPb->GetXaxis()->FindBin(meanNchpPb[cb]);
     if(RadiipPb->GetBinContent(binpPb)==0) continue;
-    //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;
+    cout<<meanNchpPb[cb]-(0.05)*meanNchpPb[cb]<<" TO "<<meanNchpPb[cb]+(0.05)*meanNchpPb[cb]<<"; "<<RadiiC2pPb->GetBinContent(binpPb)<<" +- "<<RadiiC2pPb->GetBinError(binpPb)<<" (DSYS=+"<<grRadiiC2Sys_pPb->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_pPb->GetErrorYlow(cb)<<"); ";
+    cout<<LambdaC2pPb->GetBinContent(binpPb)<<" +- "<<LambdaC2pPb->GetBinError(binpPb)<<" (DSYS=+"<<grLambdaC2Sys_pPb->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_pPb->GetErrorYlow(cb)<<"); ";
+    cout<<RadiipPb->GetBinContent(binpPb)<<" +- "<<RadiipPb->GetBinError(binpPb)<<" (DSYS="<<grRadiiSys_pPb->GetErrorY(cb)<<"); ";
+    cout<<LambdapPb->GetBinContent(binpPb)<<" +- "<<LambdapPb->GetBinError(binpPb)<<" (DSYS="<<grLambdaSys_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<<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;
+  cout<<endl;
+  cout<<"Pb--Pb:"<<endl;
+  for(int cb=19; cb>=0; cb--){
+    int binPbPb = RadiiC2PbPb->GetXaxis()->FindBin(meanNchPbPb[cb]);
+    if(RadiiPbPb->GetBinContent(binPbPb)==0) continue;
+    cout<<meanNchPbPb[cb]-(0.05)*meanNchPbPb[cb]<<" TO "<<meanNchPbPb[cb]+(0.05)*meanNchPbPb[cb]<<"; "<<RadiiC2PbPb->GetBinContent(binPbPb)<<" +- "<<RadiiC2PbPb->GetBinError(binPbPb)<<" (DSYS=+"<<grRadiiC2Sys_PbPb->GetErrorYhigh(cb)<<", -"<<grRadiiC2Sys_PbPb->GetErrorYlow(cb)<<"); ";
+    cout<<LambdaC2PbPb->GetBinContent(binPbPb)<<" +- "<<LambdaC2PbPb->GetBinError(binPbPb)<<" (DSYS=+"<<grLambdaC2Sys_PbPb->GetErrorYhigh(cb)<<", -"<<grLambdaC2Sys_PbPb->GetErrorYlow(cb)<<"); ";
+    cout<<RadiiPbPb->GetBinContent(binPbPb)<<" +- "<<RadiiPbPb->GetBinError(binPbPb)<<" (DSYS="<<grRadiiSys_PbPb->GetErrorY(cb)<<"); ";
+    cout<<LambdaPbPb->GetBinContent(binPbPb)<<" +- "<<LambdaPbPb->GetBinError(binPbPb)<<" (DSYS="<<grLambdaSys_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<<endl;
   
     
+
   
   can3->cd();
   TPad *pad3_2 = new TPad("pad3_2","pad3_2",0.0,0.0,1.,1.);
@@ -1585,18 +1599,15 @@ void Plot_plotsTPR(){
     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;
+    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)<<" TO "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetXaxis()->GetBinUpEdge(binN)<<"; "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" +- "<<C3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<" (DSYS="<<C3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<"); "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinContent(binN)<<" +- "<<c3[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<" (DSYS="<<c3_Sys[System_proof][0][ChComb_proof][KT3Bin][Mb_proof]->GetBinError(binN)<<"); ";
+      if(System_proof==0){
+       cout<<HIJING_c3_SC->GetBinContent(binN)<<" +- "<<HIJING_c3_SC->GetBinError(binN)<<"; - ;"<<endl;
+      }else{
+       cout<<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.4
     //
@@ -1842,8 +1853,8 @@ void Plot_plotsTPR(){
   
   // 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(pp_pPb_Comp) cout<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" TO "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<"; "<<c3_pp->GetBinContent(binN)<<" +- "<<c3_pp->GetBinError(binN)<<" (DSYS="<<c3_Sys[2][0][0][KT3Bin_CorrComp][Mbin_SysComp_pp]->GetBinError(binN)<<"); "<<c3_pPb->GetBinContent(binN)<<" +- "<<c3_pPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<");"<<endl;
+    else cout<<c3_pPb->GetXaxis()->GetBinLowEdge(binN)<<" TO "<<c3_pPb->GetXaxis()->GetBinUpEdge(binN)<<"; "<<c3_pPb->GetBinContent(binN)<<" +- "<<c3_pPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[1][0][0][KT3Bin_CorrComp][Mbin_SysComp_pPb]->GetBinError(binN)<<"); "<<c3_PbPb->GetBinContent(binN)<<" +- "<<c3_PbPb->GetBinError(binN)<<" (DSYS="<<c3_Sys[0][0][0][KT3Bin][Mbin_SysComp_PbPb]->GetBinError(binN)<<");"<<endl;
   }
 
   //
@@ -2014,26 +2025,26 @@ void Plot_plotsTPR(){
   grShade[2]->SetFillStyle(1000);
   grShade[2]->SetFillColor(kBlue-10);
   //
-  /*grShade[0]->Draw("f same");
+  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");
+  grRadiiC2Sys_PbPb->Draw("E 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");
+  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);
@@ -2084,13 +2095,13 @@ void Plot_plotsTPR(){
   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->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