PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <assert.h>
   7
   8 #include "TVector2.h"
   9 #include "TFile.h"
  10 #include "TString.h"
  11 #include "TF1.h"
  12 #include "TF3.h"
  13 #include "TH1.h"
  14 #include "TH2.h"
  15 #include "TH3.h"
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TMath.h"
  19 #include "TText.h"
  20 #include "TRandom3.h"
  21 #include "TArray.h"
  22 #include "TLegend.h"
  23 #include "TStyle.h"
  24 #include "TMinuit.h"
  25 #include "TCanvas.h"
  26 #include "TLatex.h"
  27 #include "TPaveStats.h"
  28 #include "TASImage.h"
  29
  30 #define BohrR 1963.6885 // Bohr Radius for pions
  31 #define FmToGeV 0.19733 // conversion to fm
  32 #define PI 3.1415926
  33 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  34
  35 using namespace std;
  36
  37 int FileSetting=0;// 0(standard), 1(r3 Lambda), 2(TTC), 3(PID), 4(Pos B), 5(Neg B), 6(GaussR8to5), 7(GaussR11to6), 8(4 kt bins old TTC cuts), 9 (FB5and7overlap), 10 (muon Variation: -0.5<NsigmaPion<2.0)
  38 bool MCcase_def=kFALSE;// MC data?
  39 int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case
  40 bool SameCharge_def=1;// 3-pion same-charge?
  41 int EDbin_def=0;// 0 or 1: Kt3 bin
  42 //
  43 bool MuonCorrection=1;// correct for Muons?
  44 bool IncludeG_def=kFALSE;// Include coherence?
  45 bool IncludeEW_def=kTRUE;// Include EdgeWorth coefficients?
  46 bool IncludeEWfromTherm_def=kFALSE;// Include EdgeWorth coefficients from Therminator?
  47 bool GRS_def=kTRUE;// Generalized RiverSide 3-body FSI correction?
  48 int Mbin_def=0;// 0-9: centrality bin in widths of 5%
  49 int Ktbin_def=10;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
  50 double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase
  51 double KShift=1.0;// K factor shift for testing (1.0 for default).  TURN OFF STRONGSC FOR THIS TESTING
  52 //
  53 //
  54 bool IncludeMJcorrection=kTRUE;// linear Mini-Jet correction for denominator of r3?
  55 bool IncludeStrongSC=kTRUE;// same-charge strong FSI
  56 bool SaveToFile_def=kFALSE;// Save outputs to file?
  57 int SourceType=0;// 0=Therminator, 1=Therminator with 50fm cut (keep at 0 for default), 2=Gaussian
  58 bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin?
  59 bool GofP=kFALSE;// Include momentum dependence of coherent fraction?
  60 bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states?
  61 bool LinkN=kFALSE;// Make N(++)=N(+-)?
  62 bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
  63 bool Tricubic=kFALSE;// Tricubic or Trilinear interpolation?  Tricubic was found to be unstable.
  64 bool IncludeSS=kTRUE;// Include same-charge two-pion correlations in the fit?
  65 bool IncludeOS=kTRUE;// Include mixed-charge two-pion correlations in the fit?
  66 //
  67 //
  68 //
  69 //
  70 const int Sbins_2=1;// 2-particle Species bins. 1=Pion-Pion only. max=6 (pi-pi, pi-k, pi-p, k-p, k-k, p-p)
  71 const int Sbins_3=1;// 3-particle Species bins. 1=Pion-Pion-Pion only. max=10
  72 const int fitlimitLowBin = 3;
  73 const int fitlimitHighBin = 14;// max = 20 (100MeV)
  74 bool LHC11h=kTRUE;// always kTRUE
  75 bool ExplicitLoop=kFALSE;// always kFALSE
  76 bool PairCut=kTRUE;// always kTRUE
  77 bool PbPbcase=kTRUE;// always kTRUE
  78 int KtbinFSI;
  79 int Ktlowbin;
  80 int Kthighbin;
  81 int bBin;// impact parameter bin
  82 int MomResCentBin;// momentum resolution centrality bin
  83 int RValue;// Radius setting for Gaussian source profile
  84 int bValue;// impact parameter for Therminator source profile (default)
  85 float TwoFrac;// Lambda parameter
  86 float OneFrac;// Lambda^{1/2}
  87 float ThreeFrac;// Lambda^{3/2}
  88 float TwoFracMomRes;// Lambda parameter for momentum resolution corrections
  89 float RValueMomRes;// Gaussian radius value for momentum resolution corrections
  90 float TwoFracr3;// Lambda parameter for r3
  91 int TPN_bin;// Two Particle Normalization bin (always 0)
  92 double Qcut_pp = 0.6;// 0.6
  93 double Qcut_PbPb = 0.1;
  94 double NormQcutLow_pp = 1.0;
  95 double NormQcutHigh_pp = 1.5;
  96 double NormQcutLow_PbPb = .15;//1.06
  97 double NormQcutHigh_PbPb = .175;//1.1
  98
  99 // single-pion emitter size (for fits with momentum dependence of coherence)
 100 double RT = 0.72/FmToGeV;
 101 double Temp = 0.135;// 0.135 GeV
 102
 103
 104 const int Nlines = 50;
 105 TH3D *CoulCorrOmega0SS;
 106 TH3D *CoulCorrOmega0OS;
 107 TH1D *CoulCorr2SS;
 108 TH1D *CoulCorr2OS;
 109 TF1 *StrongSC;// same-charge pion strong FSI
 110
 111
 112 void ReadCoulCorrections(int, int, int);
 113 void ReadCoulCorrections_Omega0();
 114 void ReadMuonCorrections(int);
 115 void ReadMomResFile(int, double);
 116 double CoulCorr2(int, double);
 117 double CoulCorrOmega0(bool, double, double, double);
 118 double CoulCorrGRS(bool, double, double, double);
 119 double OSLfitfunction(Double_t *x, Double_t *par);
 120 double D3fitfunction_3(Double_t *x, Double_t *par);
 121 double r3_fitfunction(Double_t *x, Double_t *par);
 122 double Gamov(int, double);
 123 //double LednickyCorr(double);
 124 double C2ssFitFunction(Double_t *x, Double_t *par);
 125 double C2osFitFunction(Double_t *x, Double_t *par);
 126 double cubicInterpolate(double[4], double);
 127 double bicubicInterpolate(double[4][4], double, double);
 128 double tricubicInterpolate(double[4][4][4], double, double, double);
 129 void DrawALICELogo(Bool_t, Float_t, Float_t, Float_t, Float_t);
 130 //
 131 void fcnC2_Global(int&, double*, double&, double[], int);
 132 void fcnOSL(int&, double*, double&, double[], int);
 133 float fact(float);
 134
 135 const int BINRANGE_C2global=20;
 136 double C2ssFitting[BINRANGE_C2global];
 137 double C2osFitting[BINRANGE_C2global];
 138 double C2ssFitting_e[BINRANGE_C2global];
 139 double C2osFitting_e[BINRANGE_C2global];
 140 double K2SS[BINRANGE_C2global];
 141 double K2OS[BINRANGE_C2global];
 142 double K2SS_e[BINRANGE_C2global];
 143 double K2OS_e[BINRANGE_C2global];
 144 double K3SS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
 145 double K3OS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
 146 double Chi2_C2global;
 147 double NFitPoints_C2global;
 148 double Chi2_C3global;
 149 double NFitPoints_C3global;
 150
 151 const int BINS_OSL=40;// out,side,long bins
 152 double A[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 153 double B[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 154 double A_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 155 double B_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 156 double avg_q[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 157 double K_OSL[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 158 double K_OSL_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 159 double Chi2_OSL;
 160 double NFitPoints_OSL;
 161
 162 const int BINRANGE_3=50;// q12,q13,q23
 163 double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 164 double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 165
 166 double C3_base[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 167 double N_term1[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 168 double N_term2[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 169 double N_term3[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 170 double N_term4[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 171 double N_term5[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 172
 173 double MomRes_C2_pp[40];// Qinv bins
 174 double MomRes_C2_mp[40];//
 175 double MomRes_term1_pp[40];
 176 double MomRes_term2_pp[40];
 177 double MomRes_term1_mp[40];
 178 double MomRes_term2_mp[40];
 179
 180 TH3D *MomRes3d[2][5];// SC/MC, term#
 181 TH1D *MomRes1d[2][5];// SC/MC, term#
 182 TH3D *MomRes3d_FVP[2][5];// sum#, term#
 183 TH3D *AvgK3[2];// SC/MC
 184 double AvgP[6][20];// kt bin, qinv bin
 185 int Ktbin_GofP;
 186 //
 187 double MomResDev_C2_pp[40];// Qinv bins
 188 double MomResDev_C2_mp[40];
 189 double MomResDev_term1_pp[40];
 190 double MomResDev_term2_pp[40];
 191
 192
 193 double r3_3d_arr[20][20][20];
 194 double r3_3d_arr_e[20][20][20];
 195 double AvgQinvSS[30];
 196 double AvgQinvOS[30];
 197 //
 198 TH1D *C2muonCorrection_SC;
 199 TH1D *C2muonCorrection_MC;
 200 TH1D *WeightmuonCorrection;
 201 TH1D *C3muonCorrection;
 202
 203
 204
 205
 206 void Plot_PDCumulants(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool IncludeEWfromTherm=IncludeEWfromTherm_def, bool SameCharge=SameCharge_def, bool IncludeG=IncludeG_def, bool IncludeEW=IncludeEW_def, bool GRS=GRS_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
 207
 208   EDbin_def=EDbin;
 209   SaveToFile_def=SaveToFile;
 210   MCcase_def=MCcase;
 211   IncludeEWfromTherm_def=IncludeEWfromTherm;
 212   CHARGE_def=CHARGE;
 213   IncludeG_def=IncludeG;
 214   IncludeEW_def=IncludeEW;
 215   SameCharge_def=SameCharge;// 3-pion same-charge
 216   Mbin_def=Mbin;
 217   Ktbin_def=Ktbin;
 218   if(Ktbin_def==10) Ktbin_GofP=2;
 219   else Ktbin_GofP=Ktbin_def;
 220   //
 221   Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
 222   Kthighbin=(Ktbin)*2+4;
 223
 224   //
 225   if(ConstantFSI) KtbinFSI=10;
 226   else KtbinFSI=Ktbin;
 227
 228
 229   TwoFracr3 = 0.7;// standard lambda is 0.7 for r3
 230   if(FileSetting==1) TwoFracr3 = 0.6;
 231
 232   if(Mbin == 0) {// 0-5%
 233     RValue=8;
 234     TwoFrac=0.68;
 235     TwoFracMomRes=0.68;
 236     RValueMomRes=10;// for C2
 237     bValue=2;
 238     bBin=1;
 239     MomResCentBin=1;// for C3
 240   }else if(Mbin == 1) {// 5-10%
 241     RValue=8;
 242     TwoFrac=0.68;
 243     TwoFracMomRes=0.68;
 244     RValueMomRes=9;
 245     bValue=3;
 246     bBin=2;
 247     MomResCentBin=2;
 248   }else if(Mbin <= 3) {// 10-20%
 249     RValue=7;
 250     TwoFrac=0.68;
 251     TwoFracMomRes=0.68;
 252     RValueMomRes=8;
 253     bValue=5;
 254     bBin=3;
 255     MomResCentBin=3;
 256   }else if(Mbin <= 5) {// 20-30%
 257     RValue=6;
 258     TwoFrac=0.68;
 259     TwoFracMomRes=0.68;
 260     RValueMomRes=7;
 261     bValue=7;
 262     bBin=4;
 263     MomResCentBin=4;
 264   }else if(Mbin <= 7){// 30-40%
 265     RValue=6;
 266     TwoFrac=0.68;
 267     TwoFracMomRes=0.68;
 268     RValueMomRes=6;
 269     bValue=8;
 270     bBin=5;
 271     MomResCentBin=5;
 272   }else {// 40-50%
 273     RValue=5;
 274     TwoFrac=0.68;
 275     TwoFracMomRes=0.68;
 276     RValueMomRes=5;
 277     bValue=9;
 278     bBin=6;
 279     MomResCentBin=6;
 280   }
 281
 282   OneFrac = sqrt(TwoFrac);
 283   ThreeFrac = pow(TwoFrac,3/2.);
 284
 285   // same-charge pion strong FSI (obtained from ratio of CoulStrong to Coul at R=7 Gaussian from Lednicky's code)
 286   StrongSC=new TF1("StrongSC","[0]+[1]*exp(-pow([2]*x,2))",0,50);
 287   StrongSC->FixParameter(0,9.99192e-01);
 288   StrongSC->FixParameter(1,-8.01113e-03);
 289   StrongSC->FixParameter(2,5.35912e-02);
 290
 291
 292   if(SourceType==2 && RValue > 8) {cout<<"Radius value too large!!!"<<endl; return;}
 293
 294   cout<<"Mbin = "<<Mbin<<"   Kt = "<<Ktbin<<"   SameCharge = "<<SameCharge<<endl;
 295
 296   // Core-Halo modeling of single-particle and triple-particle core fraction
 297   //float AvgN[10]={472.56, 390.824, 319.597, 265.933, 218.308, 178.865, 141.2, 115.88, 87.5556, 69.3235};// 10h (avg total pion mult)/2. 2.0sigma
 298   //float AvgN[10]={571.2, 472.7, 400.2, 325.2, 268.721, 220.3, 178.9, 143.4, 113.412, 88.0};// 11h (avg total pion mult)/2.  2.0sigma
 299
 300   //
 301   // avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values.  From unsmeared HIJING 0-5% with input signal
 302   double AvgQinvSS_temp[30]={0.00421164, 0.00796583, 0.0128473, 0.0178262, 0.0228416, 0.0276507, 0.0325368, 0.0376114, 0.0424707, 0.0476274, 0.0526015, 0.0575645, 0.0625478, 0.0675416, 0.0725359, 0.0775219, 0.0825521, 0.0875211, 0.0925303, 0.0975319, 0.102544, 0.10753, 0.112499, 0.117545, 0.122522, 0.127522, 0.132499, 0.137514, 0.142495, 0.147521};
 303   double AvgQinvOS_temp[30]={0.00352789, 0.00797596, 0.0128895, 0.0177198, 0.0226397, 0.0276331, 0.0326309, 0.0376471, 0.0426217, 0.047567, 0.0525659, 0.0575472, 0.0625886, 0.0675261, 0.0725543, 0.077564, 0.0825617, 0.0875465, 0.092539, 0.0975348, 0.102529, 0.107527, 0.112506, 0.117531, 0.122536, 0.1275, 0.132508, 0.137508, 0.14251, 0.147511};
 304   for(int i=0; i<30; i++){
 305     AvgQinvSS[i] = AvgQinvSS_temp[i];
 306     AvgQinvOS[i] = AvgQinvOS_temp[i];
 307   }
 308
 309
 310   // average total momentum in each kt and qinv bin.  Used for a test of momentum dependence of fits including coherence
 311   double AvgP_kt1[20]={0.25, 0.254421, 0.270286, 0.27346, 0.274017, 0.274252, 0.274438, 0.274675, 0.27497, 0.275354, 0.275844, 0.276433, 0.27708, 0.277768, 0.278469, 0.279179, 0.279871, 0.280536, 0.281167, 0.281749};
 312   double AvgP_kt2[20]={0.34, 0.34282, 0.370731, 0.38009, 0.381977, 0.382246, 0.382231, 0.38212, 0.381979, 0.381813, 0.381653, 0.381502, 0.381358, 0.381227, 0.381111, 0.381008, 0.380923, 0.380864, 0.380835, 0.380824};
 313   double AvgP_kt3[20]={0.47, 0.47, 0.469358, 0.486368, 0.492106, 0.49335, 0.493497, 0.49341, 0.493263, 0.493048, 0.492821, 0.492587, 0.49234, 0.492093, 0.491842, 0.491607, 0.491363, 0.491126, 0.490887, 0.490645};
 314   double AvgP_kt4[20]={0.59, 0.59, 0.59, 0.588732, 0.598683, 0.601864, 0.602459, 0.602302, 0.602005, 0.601617, 0.601192, 0.600794, 0.600365, 0.599961, 0.599571, 0.599208, 0.598841, 0.598511, 0.59816, 0.597836};
 315   double AvgP_kt5[20]={0.67, 0.67, 0.67, 0.676098, 0.694986, 0.701248, 0.703524, 0.704307, 0.704402, 0.704432, 0.704255, 0.703942, 0.703521, 0.703031, 0.702484, 0.701869, 0.701206, 0.700521, 0.699822, 0.699084};
 316   double AvgP_kt6[20]={0.79, 0.79, 0.79, 0.79, 0.790909, 0.799252, 0.801984, 0.801945, 0.800688, 0.799177, 0.797658, 0.796156, 0.794652, 0.793256, 0.791942, 0.790697, 0.789588, 0.788608, 0.787698, 0.786941};
 317   for(int ktbin=1; ktbin<=6; ktbin++){
 318     for(int qbin=1; qbin<=20; qbin++){
 319       if(ktbin==1) AvgP[ktbin-1][qbin-1] = AvgP_kt1[qbin-1];
 320       else if(ktbin==2) AvgP[ktbin-1][qbin-1] = AvgP_kt2[qbin-1];
 321       else if(ktbin==3) AvgP[ktbin-1][qbin-1] = AvgP_kt3[qbin-1];
 322       else if(ktbin==4) AvgP[ktbin-1][qbin-1] = AvgP_kt4[qbin-1];
 323       else if(ktbin==5) AvgP[ktbin-1][qbin-1] = AvgP_kt5[qbin-1];
 324       else AvgP[ktbin-1][qbin-1] = AvgP_kt6[qbin-1];
 325     }
 326   }
 327
 328
 329
 330   gStyle->SetOptStat(0);
 331   gStyle->SetOptDate(0);
 332   gStyle->SetOptFit(0111);
 333
 334   TFile *myfile;
 335   if(PbPbcase){// PbPb
 336     if(MCcase) {
 337       if(Mbin<=1){
 338         //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11.root","READ");
 339         //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11_Smeared.root","READ");
 340         //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_KtgenSignal_Rinv11.root","READ");
 341         myfile = new TFile("Results/PDC_pureHIJING_12a17a_fix_KtgenSignal_Rinv11.root","READ");
 342         //myfile = new TFile("Results/PDC_HIJING_10h8.root","READ");
 343       }else{
 344         myfile = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p68R11_C2plots.root","READ");
 345       }
 346     }else{
 347       //myfile = new TFile("Results/PDC_11h_v1paper.root","READ");// v1 paper
 348       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// Old standard (bad r3 Interp, 0.035TTC)
 349       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// v3 Standard
 350       //
 351       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins_FullTPC.root","READ");// FullTPC runlist
 352       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_HighNorm.root","READ");// High Norm variation 1.06<qinv<1.1
 353       if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// Standard
 354       if(FileSetting==1) myfile = new TFile("Results/PDC_11h_Lam0p7_and_Lam0p6.root","READ");// Lam=0.7 and Lam=0.6 (3ktbins, 0.035TTC)
 355       if(FileSetting==2) myfile = new TFile("Results/PDC_11h_3sigmaTTC.root","READ");// TTC variation
 356       if(FileSetting==3) myfile = new TFile("Results/PDC_11h_PID1p5.root","READ");// PID variation (0.035TTC)
 357       if(FileSetting==4) myfile = new TFile("Results/PDC_11h_PosB.root","READ");// Positive B field for r3num (0.035TTC)
 358       if(FileSetting==5) myfile = new TFile("Results/PDC_11h_NegB.root","READ");// Negative B field for r3num (0.035TTC)
 359       if(FileSetting==6) myfile = new TFile("Results/PDC_11h_GaussR8to5.root","READ");// Gaussian
 360       if(FileSetting==7) myfile = new TFile("Results/PDC_11h_GaussR11to6.root","READ");// Gaussian
 361       if(FileSetting==8) myfile = new TFile("Results/PDC_11h_4ktbins.root","READ");// 4 kt bins (0.035TTC)
 362       if(FileSetting==9) myfile = new TFile("Results/PDC_11h_FB5and7overlap_l0p8.root","READ");// FB5and7overlap
 363       if(FileSetting==10) myfile = new TFile("Results/PDC_11h_muonVar.root","READ");
 364     }
 365   }else{// pp
 366     cout<<"pp not currently supported"<<endl;
 367     return;
 368   }
 369
 370   ReadCoulCorrections(RValue, bValue, KtbinFSI);
 371   //ReadLednickyFile(RValue);
 372   ReadMomResFile(RValueMomRes, TwoFracMomRes);
 373   ReadCoulCorrections_Omega0();
 374   ReadMuonCorrections(Mbin);
 375   //
 376   /////////////////////////////////////////////////////////
 377
 378
 379   double NormQcutLow;
 380   double NormQcutHigh;
 381   if(PbPbcase) {
 382     NormQcutLow = NormQcutLow_PbPb;
 383     NormQcutHigh = NormQcutHigh_PbPb;
 384   }else {
 385     NormQcutLow = NormQcutLow_pp;
 386     NormQcutHigh = NormQcutHigh_pp;
 387   }
 388
 389
 390   TList *MyList;
 391   if(!MCcase){
 392     TDirectoryFile *tdir = (TDirectoryFile*)myfile->Get("PWGCF.outputChaoticityAnalysis.root");
 393     if(FileSetting!=1) MyList=(TList*)tdir->Get("ChaoticityOutput_1");
 394     else MyList=(TList*)tdir->Get("ChaoticityOutput_2");
 395   }else{
 396     MyList=(TList*)myfile->Get("MyList");
 397   }
 398   myfile->Close();
 399
 400
 401   TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
 402   //
 403
 404   cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
 405
 406
 407   TH1D *ChiSquaredNDF = new TH1D("ChiSquaredNDF","",2,0.5,2.5);// Chi^2/NDF records
 408
 409   // Explicit Loop Histos
 410   TH2D *Two_ex_2d[2][2][Sbins_2][2];
 411   TH1D *Two_ex[2][2][Sbins_2][2];
 412
 413   // Normalizations
 414   double NormH_pc[5]={0};
 415   TH3D *PionNorm[2];
 416   //TH3D *PionNorm_e[2];
 417
 418   double norm_ex_2[6][2]={{0}};
 419
 420   // 3d method histograms
 421   TH3D *Three_3d[2][2][2][Sbins_3][5];
 422   // 4-vect Product Sum
 423   TH3D *TPNRejects;
 424   ///////////////////////////////////
 425   // Create Histograms
 426
 427   // 2-particle
 428   for(int c1=0; c1<2; c1++){
 429     for(int c2=0; c2<2; c2++){
 430       for(int sc=0; sc<Sbins_2; sc++){
 431         for(int term=0; term<2; term++){
 432
 433           TString *name2_ex = new TString("Explicit2_Charge1_");
 434           *name2_ex += c1;
 435           name2_ex->Append("_Charge2_");
 436           *name2_ex += c2;
 437           name2_ex->Append("_SC_");
 438           *name2_ex += sc;
 439           name2_ex->Append("_M_");
 440           *name2_ex += Mbin;
 441           name2_ex->Append("_ED_");
 442           *name2_ex += 0;// EDbin
 443           name2_ex->Append("_Term_");
 444           *name2_ex += term+1;
 445
 446           if(sc==0 || sc==3 || sc==5){
 447             if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
 448           }
 449
 450
 451
 452           Two_ex_2d[c1][c2][sc][term] = (TH2D*)MyList->FindObject(name2_ex->Data());
 453           Two_ex_2d[c1][c2][sc][term]->Sumw2();
 454           TString *proname = new TString(name2_ex->Data());
 455           proname->Append("_pro");
 456
 457           if(Ktbin==10) {Ktlowbin=1; Kthighbin=Two_ex_2d[c1][c2][sc][term]->GetNbinsX();}
 458           Two_ex[c1][c2][sc][term] = (TH1D*)Two_ex_2d[c1][c2][sc][term]->ProjectionY(proname->Data(),Ktlowbin,Kthighbin);// bins:5-6 (kt:0.2-0.3)
 459
 460           norm_ex_2[sc][term] = Two_ex[c1][c2][sc][term]->Integral(Two_ex[c1][c2][sc][term]->GetXaxis()->FindBin(NormQcutLow),Two_ex[c1][c2][sc][term]->GetXaxis()->FindBin(NormQcutHigh));
 461           Two_ex[c1][c2][sc][term]->Scale(norm_ex_2[sc][0]/norm_ex_2[sc][term]);
 462
 463           Two_ex[c1][c2][sc][term]->SetMarkerStyle(20);
 464           Two_ex[c1][c2][sc][term]->GetXaxis()->SetTitle("q_{inv}");
 465           Two_ex[c1][c2][sc][term]->GetYaxis()->SetTitle("C_{2}");
 466           Two_ex[c1][c2][sc][term]->SetTitle("");
 467
 468         }// term
 469       }// sc
 470
 471       // 3-particle
 472       for(int c3=0; c3<2; c3++){
 473         for(int sc=0; sc<Sbins_3; sc++){
 474           for(int term=0; term<5; term++){
 475
 476             TString *name3_ex = new TString("Explicit3_Charge1_");
 477             *name3_ex += c1;
 478             name3_ex->Append("_Charge2_");
 479             *name3_ex += c2;
 480             name3_ex->Append("_Charge3_");
 481             *name3_ex += c3;
 482             name3_ex->Append("_SC_");
 483             *name3_ex += sc;
 484             name3_ex->Append("_M_");
 485             *name3_ex += Mbin;
 486             name3_ex->Append("_ED_");
 487             *name3_ex += EDbin;
 488             name3_ex->Append("_Term_");
 489             *name3_ex += term+1;
 490
 491
 492             TString *name3_pc = new TString("PairCut3_Charge1_");
 493             *name3_pc += c1;
 494             name3_pc->Append("_Charge2_");
 495             *name3_pc += c2;
 496             name3_pc->Append("_Charge3_");
 497             *name3_pc += c3;
 498             name3_pc->Append("_SC_");
 499             *name3_pc += sc;
 500             name3_pc->Append("_M_");
 501             *name3_pc += Mbin;
 502             name3_pc->Append("_ED_");
 503             *name3_pc += EDbin;
 504             name3_pc->Append("_Term_");
 505             *name3_pc += term+1;
 506
 507             ///////////////////////////////////////
 508             // skip degenerate histograms
 509             if(sc==0 || sc==6 || sc==9){// Identical species
 510               if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
 511               if( (c1+c2+c3)==2) {if(c1!=0) continue;}
 512             }else if(sc!=5){
 513               if( (c1+c2)==1) {if(c1!=0) continue;}
 514             }else {}// do nothing for pi-k-p case
 515
 516             /////////////////////////////////////////
 517
 518
 519             if(ExplicitLoop) {
 520               /*
 521                 Three_ex[c1][c2][c3][sc][term] = (TH1D*)MyList->FindObject(name3_ex->Data());
 522
 523                 Three_ex[c1][c2][c3][sc][term]->SetMarkerStyle(20);
 524                 Three_ex[c1][c2][c3][sc][term]->SetLineColor(2);
 525                 Three_ex[c1][c2][c3][sc][term]->SetMarkerColor(2);
 526                 Three_ex[c1][c2][c3][sc][term]->GetXaxis()->SetTitle("Q_{3}");
 527                 Three_ex[c1][c2][c3][sc][term]->GetYaxis()->SetTitle("C_{3}");
 528                 Three_ex[c1][c2][c3][sc][term]->SetTitle("");
 529
 530                 TString *name = new TString(name3_ex->Data());
 531                 name->Append("_Norm");
 532                 NormH_ex[term] = ((TH1D*)MyList->FindObject(name->Data()))->Integral();
 533
 534                 if(NormH_ex[term] > 0){
 535                   Three_ex[c1][c2][c3][sc][term]->Scale(NormH_ex[0]/NormH_ex[term]);
 536                 }else{ cout<<"normalization = 0.  Skipping this SC."<<endl;}
 537               */
 538             }
 539             if(PairCut){
 540
 541               TString *nameNorm=new TString("PairCut3_Charge1_");
 542               *nameNorm += c1;
 543               nameNorm->Append("_Charge2_");
 544               *nameNorm += c2;
 545               nameNorm->Append("_Charge3_");
 546               *nameNorm += c3;
 547               nameNorm->Append("_SC_");
 548               *nameNorm += sc;
 549               nameNorm->Append("_M_");
 550               *nameNorm += Mbin;
 551               nameNorm->Append("_ED_");
 552               *nameNorm += 0;// EDbin
 553               nameNorm->Append("_Term_");
 554               *nameNorm += term+1;
 555               nameNorm->Append("_Norm");
 556               NormH_pc[term] = ((TH1D*)MyList->FindObject(nameNorm->Data()))->Integral();
 557
 558               if(NormH_pc[term] > 0){
 559
 560                 if(sc<=2) {
 561
 562                   TString *name_3d = new TString(name3_pc->Data());
 563                   name_3d->Append("_3D");
 564                   Three_3d[c1][c2][c3][sc][term] = (TH3D*)MyList->FindObject(name_3d->Data());
 565                   Three_3d[c1][c2][c3][sc][term]->Sumw2();
 566                   Three_3d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]);
 567                   //cout<<"Scale factor = "<<NormH_pc[0]/NormH_pc[term]<<endl;
 568
 569                   // 4-vect Product Sum
 570                   if(c1==c2 && c1==c3 && sc==0 && term==4){
 571
 572                     TString *nameDenType=new TString("PairCut3_Charge1_");
 573                     *nameDenType += c1;
 574                     nameDenType->Append("_Charge2_");
 575                     *nameDenType += c2;
 576                     nameDenType->Append("_Charge3_");
 577                     *nameDenType += c3;
 578                     nameDenType->Append("_SC_");
 579                     *nameDenType += sc;
 580                     nameDenType->Append("_M_");
 581                     *nameDenType += Mbin;
 582                     nameDenType->Append("_ED_");
 583                     *nameDenType += EDbin;
 584                     nameDenType->Append("_TPN_");
 585                     *nameDenType += 0;
 586
 587                     PionNorm[c1] = (TH3D*)MyList->FindObject(nameDenType->Data());
 588                     PionNorm[c1]->Sumw2();
 589                     PionNorm[c1]->Scale(NormH_pc[0]/NormH_pc[term]);
 590                     if(c1==0) {
 591                       TPNRejects = (TH3D*)MyList->FindObject("fTPNRejects4");
 592                       TPNRejects->Scale(NormH_pc[0]/NormH_pc[term]);
 593                     }
 594
 595                   }
 596                 }
 597
 598               }else{
 599                 cout<<"normalization = 0.  Skipping this SC.  SC="<<sc<<"  c1="<<c1<<"  c2="<<c2<<"  c3="<<c3<<endl;
 600               }
 601
 602
 603
 604             }// pair cut
 605           }// term
 606
 607         }// SC
 608
 609
 610       }// c3
 611     }// c2
 612   }// c1
 613
 614
 615   cout<<"Done getting Histograms"<<endl;
 616
 617
 618   TCanvas *can1 = new TCanvas("can1", "can1",11,53,700,500);
 619   can1->SetHighLightColor(2);
 620   can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 621   gStyle->SetOptFit(0111);
 622   can1->SetFillColor(10);
 623   can1->SetBorderMode(0);
 624   can1->SetBorderSize(2);
 625   can1->SetGridx();
 626   can1->SetGridy();
 627   can1->SetFrameFillColor(0);
 628   can1->SetFrameBorderMode(0);
 629   can1->SetFrameBorderMode(0);
 630
 631   TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
 632   legend1->SetBorderSize(1);
 633   legend1->SetTextSize(.04);// small .03; large .036
 634   legend1->SetFillColor(0);
 635
 636   /////////////////////////////////////////////////////////
 637   /////////////////////////////////////////////////////////
 638   // This part for plotting track splitting/merging effects in MC data only
 639   /*
 640   TH3F *Merge3d_num=(TH3F*)MyList->FindObject("fPairsDetaDPhiNum");
 641   TH3F *Merge3d_den=(TH3F*)MyList->FindObject("fPairsDetaDPhiDen");
 642   //TH3F *Merge3d_num=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_num");
 643   //TH3F *Merge3d_den=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_den");
 644
 645   TH1D *Merge1d_num[10];
 646   TH1D *Merge1d_den[10];
 647   TString *newnamenum[10];
 648   TString *newnameden[10];
 649   TF1 *MergedGaus=new TF1("MergedGaus","1-[0]*exp(-pow(x/[1],2))",-0.1, 0.1);
 650   MergedGaus->SetParName(0,"Amplitude");
 651   MergedGaus->SetParName(1,"width");
 652   MergedGaus->SetParameter(0,0.06);
 653   MergedGaus->SetParameter(1,0.01);
 654   MergedGaus->SetParLimits(0,0.001,0.5);
 655   MergedGaus->SetParLimits(1,0.001,0.1);
 656
 657   for(int i=2; i<10; i++){
 658     if(i!=5 && i!=8) continue;// 5 and 8
 659     newnamenum[i]=new TString("namenum_");
 660     *newnamenum[i] += i;
 661     newnameden[i]=new TString("nameden_");
 662     *newnameden[i] += i;
 663
 664     Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionZ(newnamenum[i]->Data(),i+1,i+1,90,110);//90,110 (phi)
 665     Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionZ(newnameden[i]->Data(),i+1,i+1,90,110);// (phi)
 666     //Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionY(newnamenum[i]->Data(),i+1,i+1,190,410);// 190,410 (eta)
 667     //Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionY(newnameden[i]->Data(),i+1,i+1,190,410);// (eta)
 668     //Merge1d_num[i]->Rebin(2);
 669     //Merge1d_den[i]->Rebin(2);
 670     Merge1d_num[i]->Sumw2();
 671     Merge1d_den[i]->Sumw2();
 672     Merge1d_num[i]->SetMarkerStyle(20);
 673
 674     if(Merge1d_den[i]->Integral(1,100)<=0) continue;
 675     double SF_merge = Merge1d_num[i]->Integral(1,100)/Merge1d_den[i]->Integral(1,100);// Z projection (phi)
 676     //double SF_merge = Merge1d_num[i]->Integral(1,50)/Merge1d_den[i]->Integral(1,50);// Y projection (eta)
 677     Merge1d_den[i]->Scale(SF_merge);
 678     Merge1d_num[i]->Divide(Merge1d_den[i]);
 679
 680
 681     if(i<9){
 682       Merge1d_num[i]->SetLineColor(i+1);
 683       Merge1d_num[i]->SetMarkerColor(i+1);
 684     }else{
 685       Merge1d_num[i]->SetLineColor(11);
 686       Merge1d_num[i]->SetMarkerColor(11);
 687     }
 688     if(i==4) {
 689       Merge1d_num[i]->SetLineColor(2);
 690       Merge1d_num[i]->SetMarkerColor(2);
 691     }
 692     if(i==5) {
 693       Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#phi*");
 694       //Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#eta");
 695       Merge1d_num[i]->GetYaxis()->SetTitle("C_{2}^{HIJING}");
 696       Merge1d_num[i]->GetXaxis()->SetRangeUser(-.1,.1);
 697       Merge1d_num[i]->SetMinimum(.91);
 698       Merge1d_num[i]->SetMaximum(1.1);
 699       Merge1d_num[i]->DrawCopy();
 700
 701       //Merge1d_num[i]->Fit(MergedGaus,"IME","",-0.1,0.1);
 702     }else{
 703       Merge1d_num[i]->DrawCopy("same");
 704     }
 705
 706     TString *Dname=new TString("D=0.2*");
 707     *Dname +=i;
 708     Dname->Append(" m");
 709     legend1->AddEntry(newnamenum[i]->Data(),Dname->Data(),"lpf");
 710   }
 711   legend1->Draw("same");
 712   gStyle->SetOptFit(111);
 713   Merge1d_num[8]->Fit(MergedGaus,"IME","",-0.1,0.1);
 714   MergedGaus->Draw("same");
 715   */
 716   /*TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsPadRowNum");// kt, shfrac, qinv
 717   TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsPadRowDen");// kt, shfrac, qinv
 718   PadRowDen3->Scale(PadRowNum3->Integral(1,20,1,159, 35,40)/PadRowDen3->Integral(1,20,1,159, 35,40));
 719   PadRowNum3->GetYaxis()->SetRangeUser(0,0.01);
 720   PadRowDen3->GetYaxis()->SetRangeUser(0,0.01);
 721   TH1D *PadRowNum=(TH1D*)PadRowNum3->Project3D("z");
 722   TH1D *PadRowDen=(TH1D*)PadRowDen3->Project3D("z");
 723   PadRowNum->Divide(PadRowDen);
 724   PadRowNum->Draw();*/
 725   /*
 726   TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiNum");// r, shfrac, deltaphi
 727   TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiDen");// r, shfrac, deltaphi
 728   PadRowDen3->Scale(PadRowNum3->Integral(1,10,1,159, 90,100)/PadRowDen3->Integral(1,10,1,159, 90,100));
 729   PadRowNum3->GetXaxis()->SetRange(5,5);
 730   PadRowDen3->GetXaxis()->SetRange(5,5);
 731   TH2D *PadRowNum=(TH2D*)PadRowNum3->Project3D("zy");
 732   TH2D *PadRowDen=(TH2D*)PadRowDen3->Project3D("zy");
 733   PadRowNum->Divide(PadRowDen);
 734   PadRowNum->Draw("lego");
 735   */
 736   /////////////////////////
 737   // 2-particle legend
 738   // 0 = pi-pi
 739   // 1 = pi-k
 740   // 2 = pi-p
 741   // 3 = k-k
 742   // 4 = k-p
 743   // 5 = p-p
 744   /////////////////////////
 745   TH1D *Two_pipi_mp = (TH1D*)Two_ex[0][1][0][0]->Clone();
 746   Two_pipi_mp->Divide(Two_ex[0][1][0][1]);
 747
 748   // Normalization range counting.  Just to evaluate required normalization width in qinv.
 749   //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.06), Two_ex[0][0][0][0]->GetXaxis()->FindBin(1.1))<<endl;
 750   //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.15), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.175))<<endl;
 751   //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.15), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.175))<<endl;
 752
 753
 754   const int SCOI_2=0;
 755
 756   TH1D *Two_ex_clone_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 757   Two_ex_clone_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 758   TH1D *Two_ex_clone_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 759   Two_ex_clone_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 760   TH1D *Two_ex_clone_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 761   Two_ex_clone_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 762
 763   // Mini-jet ++ background linear estimation.
 764   TF1 *MJ_bkg_ss=new TF1("MJ_bkg_ss","pol1",0,1);
 765   Two_ex_clone_mm->Fit(MJ_bkg_ss,"IMENQ","",0.2,0.4);
 766   cout<<"Non-femto bkg C2(q=0.01) = "<<MJ_bkg_ss->Eval(0.01)<<endl;
 767
 768   Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}");
 769   Two_ex_clone_mm->SetTitle("");
 770   Two_ex_clone_mp->GetYaxis()->SetTitle("C_{2}");
 771   Two_ex_clone_mm->SetMarkerColor(2);
 772   Two_ex_clone_mm->SetLineColor(2);
 773   Two_ex_clone_mp->SetMarkerColor(1);
 774   Two_ex_clone_pp->SetMarkerColor(4);
 775   Two_ex_clone_pp->SetLineColor(4);
 776   Two_ex_clone_mm->GetXaxis()->SetRangeUser(0,0.1);
 777   Two_ex_clone_mm->SetMinimum(0.95);
 778   Two_ex_clone_mm->SetMaximum(1.4);
 779
 780   if(MCcase){
 781     Two_ex_clone_mp->SetMarkerColor(4);
 782     Two_ex_clone_mp->SetLineColor(4);
 783     Two_ex_clone_mm->Add(Two_ex_clone_pp);
 784     Two_ex_clone_mm->Scale(0.5);
 785     Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}^{HIJING}");
 786     Two_ex_clone_mm->GetYaxis()->SetTitleOffset(1.3);
 787     Two_ex_clone_mm->SetMinimum(0.97);
 788     Two_ex_clone_mm->SetMaximum(1.02);
 789     Two_ex_clone_mm->DrawCopy();
 790     //Two_ex_clone_pp->DrawCopy("same");
 791     Two_ex_clone_mp->DrawCopy("same");
 792     legend1->AddEntry(Two_ex_clone_mm,"same-charge","p");
 793     //legend1->AddEntry(Two_ex_clone_pp,"++","p");
 794     legend1->AddEntry(Two_ex_clone_mp,"mixed-charge","p");
 795     legend1->Draw("same");
 796   }
 797
 798
 799
 800   ////////////////////////////////////////////
 801   // Get Therminator EdgeWorth coefficients
 802   TString *ThermName = new TString("../ThermData/Therm_FSI_b");
 803   *ThermName += bValue;
 804   ThermName->Append(".root");
 805   TFile *Thermfile = new TFile(ThermName->Data(),"READ");
 806   TH2D *thermNum2d_ss=(TH2D*)Thermfile->Get("Num_CosFSI_ss");
 807   TH2D *thermNumSq2d_ss=(TH2D*)Thermfile->Get("NumSq_CosFSI_ss");
 808   TH2D *thermDen2d_ss=(TH2D*)Thermfile->Get("Den_ss");
 809   TH2D *thermLargeR2D_ss=(TH2D*)Thermfile->Get("LargeRpairs_ss");
 810   TH1D *thermNum_ss=(TH1D*)thermNum2d_ss->ProjectionY("thermNum_ss",Ktlowbin,Kthighbin);
 811   TH1D *thermNumSq_ss=(TH1D*)thermNumSq2d_ss->ProjectionY("thermNumSq_ss",Ktlowbin,Kthighbin);
 812   TH1D *thermDen_ss=(TH1D*)thermDen2d_ss->ProjectionY("thermDen_ss",Ktlowbin,Kthighbin);
 813   TH1D *thermLargeR_ss=(TH1D*)thermLargeR2D_ss->ProjectionY("thermLargeR_ss",Ktlowbin,Kthighbin);
 814   //
 815   TH2D *thermNum2d_os=(TH2D*)Thermfile->Get("Num_CosFSI_os");
 816   TH2D *thermNumSq2d_os=(TH2D*)Thermfile->Get("NumSq_CosFSI_os");
 817   TH2D *thermDen2d_os=(TH2D*)Thermfile->Get("Den_os");
 818   TH2D *thermLargeR2D_os=(TH2D*)Thermfile->Get("LargeRpairs_os");
 819   TH1D *thermNum_os=(TH1D*)thermNum2d_os->ProjectionY("thermNum_os",Ktlowbin,Kthighbin);
 820   TH1D *thermNumSq_os=(TH1D*)thermNumSq2d_os->ProjectionY("thermNumSq_os",Ktlowbin,Kthighbin);
 821   TH1D *thermDen_os=(TH1D*)thermDen2d_os->ProjectionY("thermDen_os",Ktlowbin,Kthighbin);
 822   TH1D *thermLargeR_os=(TH1D*)thermLargeR2D_os->ProjectionY("thermLargeR_os",Ktlowbin,Kthighbin);
 823   TH1D *C2Therm_ss = (TH1D*)thermNum_ss->Clone();
 824   TH1D *C2Therm_os = (TH1D*)thermNum_os->Clone();
 825   TH1D *C2Den_ss = (TH1D*)thermDen_ss->Clone();
 826   TH1D *C2Den_os = (TH1D*)thermDen_os->Clone();
 827   C2Therm_ss->Add(thermLargeR_ss);
 828   C2Den_ss->Add(thermLargeR_ss);
 829   C2Therm_os->Add(thermLargeR_os);
 830   C2Den_os->Add(thermLargeR_os);
 831   C2Therm_ss->Sumw2();
 832   C2Therm_os->Sumw2();
 833
 834   C2Therm_ss->Divide(C2Den_ss);
 835   C2Therm_os->Divide(C2Den_os);
 836
 837
 838   for(int i=1; i<=thermDen_ss->GetNbinsX(); i++){
 839     if(thermDen_ss->GetBinContent(i) > 0){
 840       double err = thermNumSq_ss->GetBinContent(i)/thermDen_ss->GetBinContent(i);
 841       err -= pow(thermNum_ss->GetBinContent(i)/thermDen_ss->GetBinContent(i),2);
 842       err /= thermDen_ss->GetBinContent(i);
 843       err = sqrt(fabs(err));
 844       C2Therm_ss->SetBinError(i, err);
 845     }
 846     if(thermDen_os->GetBinContent(i) > 0){
 847       double err = thermNumSq_os->GetBinContent(i)/thermDen_os->GetBinContent(i);
 848       err -= pow(thermNum_os->GetBinContent(i)/thermDen_os->GetBinContent(i),2);
 849       err /= thermDen_os->GetBinContent(i);
 850       err = sqrt(fabs(err));
 851       C2Therm_os->SetBinError(i, err);
 852     }
 853   }
 854
 855   C2Therm_ss->SetMarkerStyle(20);
 856   C2Therm_os->SetMarkerStyle(25);
 857   C2Therm_ss->SetMarkerColor(2);
 858   C2Therm_os->SetMarkerColor(4);
 859   C2Therm_ss->SetLineColor(2);
 860   C2Therm_os->SetLineColor(4);
 861   C2Therm_ss->SetMarkerSize(1.5);
 862   C2Therm_os->SetMarkerSize(1.5);
 863   C2Therm_ss->SetDirectory(0);
 864   C2Therm_os->SetDirectory(0);
 865   C2Therm_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 866   C2Therm_ss->GetYaxis()->SetTitle("1+<cos(qx)>");
 867   C2Therm_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 868   C2Therm_os->GetYaxis()->SetTitle("C_{2}^{+-}");
 869   C2Therm_os->SetMinimum(0.98);
 870   C2Therm_os->SetMaximum(1.25);
 871   //C2Therm_ss->DrawCopy();
 872   //C2Therm_os->DrawCopy("same");
 873   C2Therm_ss->SetDirectory(0);
 874   C2Therm_os->SetDirectory(0);
 875   //C2Therm->Fit(GaussEW,"IME","",0.005,0.1);
 876   Thermfile->Close();
 877
 878
 879
 880   /////////////////////////////////////////////////////
 881   // Global fitting C2os and C2ss
 882   double ThermEW[4]={0};
 883   double C2ssSys_e[BINRANGE_C2global]={0};
 884   double C2osSys_e[BINRANGE_C2global]={0};
 885   TF1 *fitC2ssTherm;
 886   TF1 *fitC2osTherm;
 887   TF1 *fitC2ssThermGaus;
 888   TF1 *fitC2osThermGaus;
 889   //
 890   const int npar=11;// 10
 891   TMinuit MyMinuit(npar);
 892   MyMinuit.SetFCN(fcnC2_Global);
 893   double OutputPar[npar]={0};
 894   double OutputPar_e[npar]={0};
 895
 896   double par[npar];               // the start values
 897   double stepSize[npar];          // step sizes
 898   double minVal[npar];            // minimum bound on parameter
 899   double maxVal[npar];            // maximum bound on parameter
 900   string parName[npar];
 901
 902   TH1D *temp_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 903   temp_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 904   TH1D *temp_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 905   temp_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 906   TH1D *temp_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 907   temp_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 908   TH1D *C2ssRaw=(TH1D*)temp_mm->Clone();// MRC uncorrected
 909   TH1D *C2osRaw=(TH1D*)temp_mp->Clone();// MRC uncorrected
 910   C2ssRaw->SetMarkerStyle(24);
 911   C2osRaw->SetMarkerStyle(21);//21
 912   C2ssRaw->SetMarkerColor(2);
 913   C2osRaw->SetMarkerColor(4);
 914
 915
 916   int iteration=0;
 917   while(iteration<3){// 0=pure Gaussian Therminator fits, 1=EW Therminator fits, 2=real data fits with Therm EW
 918     if(!IncludeEWfromTherm) iteration=2;// skip Therminator
 919
 920     for(int i=0; i<BINRANGE_C2global; i++){
 921
 922       C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.;
 923       if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= (MomRes_C2_pp[i]-1)*MRCShift+1;
 924       C2ssFitting_e[i] = pow(MomRes_C2_pp[i]*sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2);
 925       C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.);
 926       C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2));
 927       C2ssFitting_e[i] += pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
 928       C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]);
 929       C2osFitting[i] = temp_mp->GetBinContent(i+1);
 930       if(!GeneratedSignal && !MCcase) C2osFitting[i] *= (MomRes_C2_mp[i]-1)*MRCShift+1;
 931       C2osFitting_e[i] = pow(MomRes_C2_mp[i]*temp_mp->GetBinError(i+1),2);
 932       C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1));
 933       C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1));
 934       C2osFitting_e[i] += pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 935       C2osFitting_e[i] = sqrt(C2osFitting_e[i]);
 936       //
 937       C2ssSys_e[i] = pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
 938       C2ssSys_e[i] = sqrt(C2ssSys_e[i]);
 939       C2osSys_e[i] = pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 940       C2osSys_e[i] = sqrt(C2osSys_e[i]);
 941       //
 942
 943       K2SS[i] = CoulCorr2(+1, AvgQinvSS[i]);
 944       K2OS[i] = CoulCorr2(-1, AvgQinvOS[i]);
 945       //K2SS[i] = 1;
 946       //K2OS[i] = 1;
 947       //
 948       K2SS_e[i] = sqrt(pow((K2SS[i]-1)*0.02,2) + pow((K2SS[i]-Gamov(+1, AvgQinvSS[i]))*0.02,2));
 949       K2OS_e[i] = sqrt(pow((K2OS[i]-1)*0.02,2) + pow((K2OS[i]-Gamov(-1, AvgQinvOS[i]))*0.02,2));
 950       //K2SS_e[i] = 0.0001;
 951       //K2OS_e[i] = 0.0001;
 952
 953       //
 954       if(iteration<2){
 955         C2ssFitting[i] = C2Therm_ss->GetBinContent(i+1);// Therminator fit
 956         C2ssFitting_e[i] = C2Therm_ss->GetBinError(i+1);// Therminator fit
 957         C2osFitting[i] = C2Therm_os->GetBinContent(i+1);// Therminator fit
 958         C2osFitting_e[i] = C2Therm_os->GetBinError(i+1);// Therminator fit
 959         K2SS_e[i] = 0.;
 960         K2OS_e[i] = 0.;
 961       }
 962
 963     }
 964
 965
 966
 967     C2ssFitting[0]=-1000; C2osFitting[0]=-1000;
 968     C2ssFitting_e[0]=1000; C2osFitting_e[0]=1000;
 969     K2SS_e[0]=1000; K2OS_e[0]=1000;
 970
 971
 972
 973     par[0] = 1.0; par[1]=0.5; par[2]=0.5; par[3]=9.2; par[4] = .1; par[5] = .2; par[6] = .0; par[7] = 0.; par[8] = 0.; par[9] = 1.; par[10] = 0.01;
 974     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; stepSize[10]=0.01;
 975     minVal[0] = 0.995; minVal[1] = 0.2; 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[10] = 0;
 976     maxVal[0] = 1.03; maxVal[1] = 1.0; 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.03; maxVal[10]=1.;
 977     parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh";
 978     parName[5] = "kappa_3"; parName[6] = "kappa_4"; parName[7] = "kappa_5"; parName[8]="kappa_6"; parName[9]="Norm_2"; parName[10]="P_{coh}";
 979
 980     MyMinuit.DefineParameter(10, parName[10].c_str(), 0., stepSize[10], minVal[10], maxVal[10]); MyMinuit.FixParameter(10);// extra parameter
 981
 982     for (int i=0; i<npar; i++){
 983       MyMinuit.DefineParameter(i, parName[i].c_str(), par[i], stepSize[i], minVal[i], maxVal[i]);
 984     }
 985     //MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
 986     //MyMinuit.DefineParameter(1, parName[1].c_str(), .68, stepSize[1], minVal[1], maxVal[1]); MyMinuit.FixParameter(1);// lambda
 987     //MyMinuit.DefineParameter(2, parName[2].c_str(), 0.9999, stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
 988
 989     if(IncludeG==kFALSE || iteration<2) {
 990       MyMinuit.DefineParameter(2, parName[2].c_str(), 0., stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
 991       MyMinuit.DefineParameter(4, parName[4].c_str(), 0., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
 992     }else{
 993       Int_t err=0;
 994       Double_t tmp[1];
 995       tmp[0] = 2+1;
 996       MyMinuit.mnexcm( "RELease", tmp,  1,  err );// G
 997       tmp[0] = 4+1;
 998       MyMinuit.mnexcm( "RELease", tmp,  1,  err );// Rcoh
 999     }
1000
1001     if(!IncludeEW){
1002       if(!IncludeEWfromTherm){
1003         // kappa3=0.24, kappa4=0.16 for testing
1004         MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1005         MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1006         MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
1007         MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1008       }else{
1009         if(iteration==0){
1010           MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1011           MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1012         }else if(iteration==1){
1013           Int_t err=0;
1014           Double_t tmp[1];
1015           tmp[0] = 5+1;
1016           MyMinuit.mnexcm( "RELease", tmp,  1,  err );// k3
1017           tmp[0] = 6+1;
1018           MyMinuit.mnexcm( "RELease", tmp,  1,  err );// k4
1019         }else{// iteration=2
1020           MyMinuit.DefineParameter(5, parName[5].c_str(), ThermEW[0], stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1021           MyMinuit.DefineParameter(6, parName[6].c_str(), ThermEW[1], stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1022           MyMinuit.DefineParameter(7, parName[7].c_str(), ThermEW[2], stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
1023           MyMinuit.DefineParameter(8, parName[8].c_str(), ThermEW[3], stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1024         }
1025       }// IncludeEWfromTherm
1026     }// IncludeEW
1027
1028     if(IncludeSS==kFALSE){
1029       MyMinuit.DefineParameter(3, parName[3].c_str(), 9.1, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
1030       MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
1031       MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1032       MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1033     }
1034     if(IncludeOS==kFALSE){
1035       MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2
1036     }
1037
1038     //MyMinuit.DefineParameter(3, parName[3].c_str(), 10.5, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
1039     //MyMinuit.DefineParameter(4, parName[4].c_str(), 5., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
1040     //MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1041     //MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1042     MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
1043     MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1044     //
1045     int ierflg=0;
1046     double arglist[10];
1047     //arglist[0]=2;// improve Minimization Strategy (1 is default)
1048     //MyMinuit.mnexcm("SET STR",arglist,1,ierflg);
1049     //arglist[0] = 0;
1050     //MyMinuit.mnexcm("SCAN", arglist,1,ierflg);
1051     arglist[0] = 5000;
1052     MyMinuit.mnexcm("MIGRAD", arglist ,1,ierflg);
1053     // Do the minimization!
1054     cout<<"Start C2 Global fit"<<endl;
1055     MyMinuit.Migrad();// generally the best minimization algorithm
1056     cout<<"End C2 Global fit"<<endl;
1057
1058     for (int i=0; i<npar; i++){
1059       MyMinuit.GetParameter(i,OutputPar[i],OutputPar_e[i]);
1060     }
1061
1062     cout<<"Chi2/NDF = "<<Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars())<<"   Chi^2="<<Chi2_C2global<<endl;
1063     if(iteration==2) {
1064       ChiSquaredNDF->Fill(1, Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars()));
1065     }if(iteration==1) {
1066       ThermEW[0]=OutputPar[5]; ThermEW[1]=OutputPar[6]; ThermEW[2]=OutputPar[7]; ThermEW[3]=OutputPar[8];
1067       fitC2ssTherm = new TF1("fitC2ssTherm",C2ssFitFunction, 0.005,0.2, npar);
1068       for(int i=0; i<npar; i++) {
1069         fitC2ssTherm->FixParameter(i,OutputPar[i]);
1070         fitC2ssTherm->SetParError(i,OutputPar_e[i]);
1071       }
1072       fitC2osTherm = new TF1("fitC2osTherm",C2osFitFunction, 0.005,0.2, npar);
1073       for(int i=0; i<npar; i++) {
1074         fitC2osTherm->FixParameter(i,OutputPar[i]);
1075         fitC2osTherm->SetParError(i,OutputPar_e[i]);
1076       }
1077       fitC2osTherm->SetLineColor(4);
1078     }
1079     if(iteration==0){
1080       fitC2ssThermGaus = new TF1("fitC2ssThermGaus",C2ssFitFunction, 0.005,0.2, npar);
1081       for(int i=0; i<npar; i++) {
1082         fitC2ssThermGaus->FixParameter(i,OutputPar[i]);
1083         fitC2ssThermGaus->SetParError(i,OutputPar_e[i]);
1084       }
1085       fitC2osThermGaus = new TF1("fitC2osThermGaus",C2osFitFunction, 0.005,0.2, npar);
1086       for(int i=0; i<npar; i++) {
1087         fitC2osThermGaus->FixParameter(i,OutputPar[i]);
1088         fitC2osThermGaus->SetParError(i,OutputPar_e[i]);
1089       }
1090       fitC2osThermGaus->SetLineColor(4);
1091     }
1092
1093     iteration++;
1094
1095   }// iteration loop
1096
1097
1098   TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone();
1099   TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone();
1100   TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone();
1101   TH1D *C2_os_Momsys=(TH1D*)Two_ex_clone_mp->Clone();
1102   TH1D *C2_ss_Ksys=(TH1D*)Two_ex_clone_mm->Clone();
1103   TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone();
1104   TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone();
1105   TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone();
1106
1107   for(int i=0; i<BINRANGE_C2global; i++){
1108     C2_ss->SetBinContent(i+1, C2ssFitting[i]);
1109     C2_os->SetBinContent(i+1, C2osFitting[i]);
1110     C2_ss_Momsys->SetBinContent(i+1, C2ssFitting[i]);
1111     C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]);
1112     C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]);
1113     C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]);
1114     double Toterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
1115     double Toterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
1116     C2_ss->SetBinError(i+1, Toterror_ss);
1117     C2_os->SetBinError(i+1, Toterror_os);
1118     C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]);
1119     C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]);
1120     C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]);
1121     C2_os_Ksys->SetBinError(i+1, OutputPar[1]*K2OS_e[i]);
1122     //
1123     K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0);
1124     K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0);
1125   }
1126
1127   C2_ss_Momsys->SetMarkerSize(0);
1128   C2_ss_Momsys->SetFillStyle(1000);
1129   C2_ss_Momsys->SetFillColor(kRed-10);
1130   C2_os_Momsys->SetMarkerSize(0);
1131   C2_os_Momsys->SetFillStyle(1000);
1132   C2_os_Momsys->SetFillColor(kBlue-10);
1133   C2_ss_Ksys->SetMarkerSize(0);
1134   C2_ss_Ksys->SetFillStyle(3004);
1135   C2_ss_Ksys->SetFillColor(kRed);
1136   C2_os_Ksys->SetMarkerSize(0);
1137   C2_os_Ksys->SetFillStyle(3004);
1138   C2_os_Ksys->SetFillColor(kBlue);
1139
1140
1141
1142   C2_ss->GetXaxis()->SetRangeUser(0,0.098);
1143   C2_ss->GetYaxis()->SetRangeUser(0.98,1.3);
1144   C2_ss->GetXaxis()->SetTitleOffset(.8);
1145   C2_ss->GetYaxis()->SetTitleOffset(.8);
1146   C2_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1147   C2_ss->GetXaxis()->SetTitleSize(0.05);
1148   C2_ss->GetYaxis()->SetTitleSize(0.05);
1149   C2_os->GetXaxis()->SetRangeUser(0,0.098);
1150   C2_os->GetYaxis()->SetRangeUser(0.98,1.3);
1151   C2_os->GetXaxis()->SetTitleOffset(.8);
1152   C2_os->GetYaxis()->SetTitleOffset(.8);
1153   C2_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1154   C2_os->GetXaxis()->SetTitleSize(0.05);
1155   C2_os->GetYaxis()->SetTitleSize(0.05);
1156
1157   C2_ss->SetMarkerSize(1.5);
1158   C2_os->SetMarkerSize(1.5);
1159   C2_os->SetMarkerStyle(25);
1160   C2_os->SetMarkerColor(4);
1161
1162   TF1 *fitC2ss = new TF1("fitC2ss",C2ssFitFunction, 0.001,0.2, npar);
1163   TF1 *fitC2os = new TF1("fitC2os",C2osFitFunction, 0.001,0.2, npar);
1164   for(int i=0; i<npar; i++) {
1165     fitC2ss->FixParameter(i,OutputPar[i]);
1166     fitC2ss->SetParError(i,OutputPar_e[i]);
1167   }
1168   for(int i=0; i<npar; i++) {
1169     fitC2os->FixParameter(i,OutputPar[i]);
1170     fitC2os->SetParError(i,OutputPar_e[i]);
1171   }
1172   TH1D *fitC2ss_h=new TH1D("fitC2ss_h","",C2_ss->GetNbinsX(),0,C2_ss->GetXaxis()->GetBinUpEdge(C2_ss->GetNbinsX()));
1173   TH1D *fitC2os_h=new TH1D("fitC2os_h","",C2_os->GetNbinsX(),0,C2_os->GetXaxis()->GetBinUpEdge(C2_os->GetNbinsX()));
1174   fitC2ss_h->SetLineWidth(2); fitC2os_h->SetLineWidth(2);
1175   fitC2ss_h->SetLineColor(2); fitC2os_h->SetLineColor(4);
1176
1177   for(int bin=1; bin<=fitC2ss_h->GetNbinsX(); bin++){
1178     double qinv=(bin-0.5)*0.005;
1179     fitC2ss_h->SetBinContent(bin, fitC2ss->Eval(qinv));
1180     fitC2os_h->SetBinContent(bin, fitC2os->Eval(qinv));
1181   }
1182
1183   if(!MCcase){
1184     C2_ss->DrawCopy();
1185     legend1->AddEntry(C2_ss,"same-charge","p");
1186     C2_os->DrawCopy("same");
1187     legend1->AddEntry(C2_os,"opp-charge","p");
1188     C2_ss_Momsys->DrawCopy("E2 same");
1189     C2_os_Momsys->DrawCopy("E2 same");
1190     C2_ss_Ksys->DrawCopy("E2 same");
1191     C2_os_Ksys->DrawCopy("E2 same");
1192     C2_ss->DrawCopy("same");
1193     C2_os->DrawCopy("same");
1194     fitC2os->SetLineColor(4);
1195     fitC2ss_h->DrawCopy("C same");
1196     fitC2os_h->DrawCopy("C same");
1197     MJ_bkg_ss->SetLineColor(1);
1198     MJ_bkg_ss->Draw("same");
1199     legend1->AddEntry(MJ_bkg_ss,"Non-femto bkg","l");
1200     legend1->Draw("same");
1201   }
1202   //C2Therm_ss->DrawCopy();
1203   //C2Therm_os->DrawCopy("same");
1204
1205   /*
1206   ////////// C2 systematics
1207   // C2 -- base, M0, Pos B (0.035 TTC for all)
1208   //double C2ss_base[20]={0, 1.17121, 1.17272, 1.13947, 1.09819, 1.06512, 1.04115, 1.02563, 1.01558, 1.00937, 1.00551, 1.00309, 1.00157, 1.00065, 1.00023, 0.999911, 0.99976, 0.999733, 0.999593, 0.999641};
1209   //double C2ss_base_e[20]={0, 0.00135764, 0.000430769, 0.000242626, 0.000167895, 0.00012913, 0.000105724, 9.01743e-05, 7.90587e-05, 7.07036e-05, 6.41813e-05, 5.89418e-05, 5.46406e-05, 5.10458e-05, 4.80009e-05, 4.53857e-05, 4.31177e-05, 4.11375e-05, 3.93869e-05, 3.78413e-05};
1210   // --, M0,
1211   double C2ss_base[20]={0, 1.16926, 1.17306, 1.13932, 1.09834, 1.06475, 1.04096, 1.02546, 1.01547, 1.00931, 1.00542, 1.00311, 1.00155, 1.00071, 1.00018, 0.999786, 0.999742, 0.999612, 0.999548, 0.999631};
1212   double C2ss_base_e[20]={0, 0.00128558, 0.000410643, 0.000231512, 0.000160337, 0.000123314, 0.000100987, 8.61441e-05, 7.55245e-05, 6.75439e-05, 6.13117e-05, 5.63119e-05, 5.21976e-05, 4.87656e-05, 4.58515e-05, 4.33492e-05, 4.11861e-05, 3.92891e-05, 3.76189e-05, 3.61428e-05};
1213   // -+, M0,
1214   //double C2ss_base[20]={1.97333, 1.17447, 1.09412, 1.05561, 1.03613, 1.02507, 1.01809, 1.01363, 1.01048, 1.00816, 1.00645, 1.00526, 1.00427, 1.0035, 1.00289, 1.00241, 1.00209, 1.0018, 1.00148, 1.00127};
1215   //double C2ss_base_e[20]={1.61122, 0.00030054, 0.000174753, 0.000123051, 9.53599e-05, 7.81537e-05, 6.64323e-05, 5.79583e-05, 5.15436e-05, 4.65277e-05, 4.25025e-05, 3.92115e-05, 3.64672e-05, 3.4149e-05, 3.21673e-05, 3.04575e-05, 2.89698e-05, 2.76642e-05, 2.65087e-05, 2.54847e-05};
1216   // --, M0, kt4
1217   //double C2ss_base[20]={0, 0, 0, 1.15387, 1.14361, 1.10213, 1.06552, 1.04357, 1.02804, 1.01826, 1.01182, 1.00767, 1.00509, 1.00294, 1.00169, 1.0008, 1.00063, 1.00022, 1.00001, 0.999986};
1218   //double C2ss_base_e[20]={0, 0, 0, 0.00229526, 0.000813871, 0.000518539, 0.00039803, 0.000328676, 0.000281681, 0.00024771, 0.000221807, 0.000201394, 0.000184824, 0.000171048, 0.000159468, 0.000149554, 0.000141058, 0.000133583, 0.00012702, 0.000121218};
1219   // --, M0, kt5
1220   //double C2ss_base[20]={0, 0, 0, 1.11133, 1.13824, 1.11806, 1.08288, 1.05648, 1.03774, 1.02516, 1.01728, 1.01166, 1.00841, 1.00553, 1.00357, 1.00194, 1.00132, 1.00104, 1.00058, 1.00039};
1221   //double C2ss_base_e[20]={0, 0, 0, 0.0147639, 0.00188068, 0.000910766, 0.000637117, 0.000510753, 0.000432397, 0.000377855, 0.000337624, 0.000306333, 0.000281495, 0.00026099, 0.000243881, 0.000229331, 0.000216961, 0.000206243, 0.000196785, 0.000188487};
1222   // --, M0, kt6
1223   //double C2ss_base[20]={0, 0, 0, 0, 1.08106, 1.12109, 1.09838, 1.07193, 1.05087, 1.03534, 1.02433, 1.01785, 1.01339, 1.00928, 1.00662, 1.00441, 1.00326, 1.00192, 1.00175, 1.00143};
1224   //double C2ss_base_e[20]={0, 0, 0, 0, 0.00722887, 0.00211091, 0.00123591, 0.000924715, 0.000765679, 0.000662277, 0.000588573, 0.000533669, 0.000490831, 0.000456049, 0.000427965, 0.000404612, 0.000385285, 0.000368671, 0.00035474, 0.000342708};
1225
1226   cout<<"C2 values:"<<endl;
1227   for(int ii=0; ii<20; ii++){
1228     cout<<Two_ex_clone_mm->GetBinContent(ii+1)<<", ";
1229   }
1230   cout<<endl;
1231   cout<<"C2 errors:"<<endl;
1232   for(int ii=0; ii<20; ii++){
1233     cout<<Two_ex_clone_mm->GetBinError(ii+1)<<", ";
1234   }
1235   cout<<endl;
1236   TH1D *C2ss_Sys=(TH1D*)Two_ex_clone_mm->Clone();
1237   for(int ii=1; ii<20; ii++){
1238     if(C2ss_base[ii] ==0) {
1239       C2ss_Sys->SetBinContent(ii+1, 100.);
1240       C2ss_Sys->SetBinError(ii+1, 100.);
1241       continue;
1242     }
1243     C2ss_Sys->SetBinContent(ii+1, (C2ss_base[ii]-C2ss_Sys->GetBinContent(ii+1))/C2ss_base[ii]);
1244     C2ss_Sys->SetBinError(ii+1, sqrt(fabs(pow(C2ss_Sys->GetBinError(ii+1),2) - pow(C2ss_base_e[ii],2))));
1245   }
1246   C2ss_Sys->GetXaxis()->SetRangeUser(0,0.099);
1247   C2ss_Sys->GetYaxis()->SetTitle("(C_{2}^{s}-C_{2}^{v})/C_{2}^{s}");
1248   C2ss_Sys->GetYaxis()->SetTitleOffset(1.3);
1249   C2ss_Sys->SetMinimum(-0.01);
1250   C2ss_Sys->SetMaximum(0.01);
1251   C2ss_Sys->Draw();
1252   TF1 *C2lineSys=new TF1("C2lineSys","pol0",0,0.099);
1253   //C2ss_Sys->Fit(C2lineSys,"MEQ","",0,0.099);
1254   */
1255
1256   // Momentum resolution uncorrected C2
1257   /*C2ssRaw->Draw("same");
1258   C2osRaw->Draw("same");
1259   legend1->AddEntry(C2ssRaw,"same-charge, Raw","p");
1260   legend1->AddEntry(C2osRaw,"opp-charge, Raw","p");
1261   legend1->Draw("same");
1262   */
1263
1264   // FSI corrected C2+-
1265   /*TH1D *C2_os_FSIcorrected=(TH1D*)C2_os->Clone();
1266   for(int ii=2; ii<20; ii++){
1267     double value = (C2_os_FSIcorrected->GetBinContent(ii) - (1-OutputPar[1]))/(OutputPar[1]*K2OS[ii-1]);
1268     C2_os_FSIcorrected->SetBinContent(ii,value);
1269   }
1270   C2_os_FSIcorrected->SetMinimum(0.95);
1271   C2_os_FSIcorrected->SetMarkerStyle(24);
1272   C2_os_FSIcorrected->Draw();
1273   */
1274
1275
1276   /*Two_ex_clone_mm->SetMarkerStyle(25);
1277   Two_ex_clone_mp->SetMarkerStyle(25);
1278   Two_ex_clone_mm->SetMarkerColor(2);
1279   Two_ex_clone_mp->SetMarkerColor(4);
1280   Two_ex_clone_mm->Draw("same");
1281   Two_ex_clone_mp->Draw("same");
1282   legend1->AddEntry(C2_ss,"same-charge; Momentum Resolution Corrected","p");
1283   legend1->AddEntry(C2_os,"opp-charge; Momentum Resolution Corrected","p");
1284   legend1->AddEntry(Two_ex_clone_mm,"same-charge; Raw","p");
1285   legend1->AddEntry(Two_ex_clone_mp,"opp-charge; Raw","p");
1286   legend1->Draw("same");
1287   */
1288   //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1289   //TGraph *gr2 = (TGraph*)MyMinuit.Contour(15,1,2);
1290   //gr2->GetXaxis()->SetTitle("lambda");
1291   //gr2->GetYaxis()->SetTitle("G");
1292   //gr2->SetTitle("");
1293   //gr2->SetFillColor(kYellow);
1294   //gr2->Draw("alf");
1295   //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1296   //MyMinuit.SetErrorDef(1);
1297   //TGraph *gr1 = (TGraph*)MyMinuit.Contour(15,1,2);
1298   //gr1->SetFillColor(kGreen);//38
1299   //gr1->Draw("lf");
1300
1301
1302
1303   ///////////////////////////
1304   // STAR comparison
1305   /* double C2_star_mm[50]={0, 0, 0.566654, 1.11195, 1.09964, 1.12566, 1.1479, 1.15596, 1.15658, 1.14583, 1.13182, 1.11544, 1.09851, 1.08326, 1.0701, 1.05686, 1.048, 1.03924, 1.03122, 1.02544, 1.02019, 1.01617, 1.01239, 1.00942, 1.00736, 1.00485, 1.00296, 1.00152, 1.00039, 0.999461, 0.998703, 0.997563, 0.997294, 0.996726, 0.996739, 0.996355, 0.996426, 0.996581, 0.996336, 0.996157, 0.996312, 0.996413, 0.996501, 0.996204, 0.996816, 0.996654, 0.99695, 0.997494, 0.997135, 0.997188};
1306   double C2_star_mm_e[50]={0, 0, 0.454334, 0.0216601, 0.00590084, 0.00297232, 0.001939, 0.00142608, 0.00112851, 0.000929702, 0.000792097, 0.000690265, 0.000613054, 0.000552755, 0.000504494, 0.000464031, 0.000431371, 0.000403169, 0.000378672, 0.000357685, 0.000339133, 0.00032286, 0.000308251, 0.000295149, 0.000283433, 0.000272615, 0.000262784, 0.000253889, 0.000245693, 0.00023813, 0.000231135, 0.00022455, 0.000218593, 0.000212961, 0.000207822, 0.000202906, 0.000198339, 0.000194087, 0.000189971, 0.000186181, 0.000182575, 0.000179202, 0.000175989, 0.000172894, 0.000170072, 0.000167293, 0.000164706, 0.000162286, 0.000159867, 0.000157606};
1307   TH1D *Two_star_mm=(TH1D*)Two_ex_clone_mm->Clone();
1308   Two_star_mm->Add(Two_star_mm,-1);
1309   for(int i=0; i<50; i++) {Two_star_mm->SetBinContent(i+1, C2_star_mm[i]); Two_star_mm->SetBinError(i+1, C2_star_mm_e[i]);}
1310   Two_star_mm->SetMarkerColor(4);
1311   Two_star_mm->SetLineColor(4);
1312   //Two_star_mm->Draw("same");
1313   //legend1->AddEntry(Two_star_mm,"--, 200 GeV","p");
1314
1315   //Two_ex_clone_mp->Multiply(Two_ex_clone_pp);
1316   //Two_ex_clone_mp->Draw();
1317   //Two_ex_clone_mm->Multiply(Two_ex_clone_mp);
1318   //Two_ex_clone_mm->Draw();
1319   //Two_ex_clone_pp->Draw("same");
1320
1321   //legend1->Draw("same");
1322   */
1323
1324
1325
1326
1327
1328   /////////////////////////////////////////////////////////////////////////
1329   // 3-d fitting (out,side,long).  Not used for paper.
1330   /*
1331   TString *name1 = new TString("Explicit2_Charge1_0_Charge2_0_SC_0_M_");
1332   *name1 += Mbin;
1333   TString *name2 = new TString(name1->Data());
1334   TString *name3 = new TString(name1->Data());
1335   name1->Append("_ED_0_Term_1_osl_b2");// b1 (0.2<kt<0.3). b2 (0.6<kt<0.7)
1336   name2->Append("_ED_0_Term_2_osl_b2");
1337   name3->Append("_ED_0_Term_1_osl_b2_QW");
1338
1339
1340   TH3D *num_osl = (TH3D*)MyList->FindObject(name1->Data());
1341   TH3D *den_osl = (TH3D*)MyList->FindObject(name2->Data());
1342   den_osl->Scale(num_osl->Integral(28,40, 28,40, 28,40)/den_osl->Integral(28,40, 28,40, 28,40));
1343   TH3D *num_osl_QW = (TH3D*)MyList->FindObject(name3->Data());
1344
1345   for(int i=0; i<BINS_OSL; i++){
1346     for(int j=0; j<BINS_OSL; j++){
1347       for(int k=0; k<BINS_OSL; k++){
1348         if(num_osl->GetBinContent(i+1,j+1,k+1) < 1 || den_osl->GetBinContent(i+1,j+1,k+1) < 1) continue;
1349
1350         avg_q[i][j][k] = num_osl_QW->GetBinContent(i+1,j+1,k+1)/num_osl->GetBinContent(i+1,j+1,k+1);
1351         int QinvBin = int((avg_q[i][j][k])/0.005);
1352         if(QinvBin >=20) QinvBin=19;
1353         if(MomRes_term1_pp[QinvBin] ==0) continue;
1354         if(MomRes_term2_pp[QinvBin] ==0) continue;
1355         //
1356         num_osl->SetBinContent(i+1,j+1,k+1, num_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term1_pp[QinvBin]);
1357         den_osl->SetBinContent(i+1,j+1,k+1, den_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term2_pp[QinvBin]);
1358         //
1359         A[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1360         B[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1361         //
1362         A_e[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1363         A_e[i][j][k] += pow(num_osl->GetBinContent(i+1,j+1,k+1)/MomRes_term1_pp[QinvBin] *fabs(MomRes_term1_pp[QinvBin]-1)*0.1,2);
1364         A_e[i][j][k] = sqrt(A_e[i][j][k]);
1365         B_e[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1366         B_e[i][j][k] += pow(den_osl->GetBinContent(i+1,j+1,k+1)/MomRes_term2_pp[QinvBin] *fabs(MomRes_term2_pp[QinvBin]-1)*0.1,2);
1367         B_e[i][j][k] = sqrt(B_e[i][j][k]);
1368         //
1369         K_OSL[i][j][k] = CoulCorr2(+1, avg_q[i][j][k]);
1370         K_OSL_e[i][j][k] = sqrt(pow((K_OSL[i][j][k]-1)*0.02,2) + pow((K_OSL[i][j][k]-Gamov(+1,avg_q[i][j][k]))*0.01,2));
1371         //K_OSL_e[i][j][k] = 0;
1372       }
1373     }
1374   }
1375
1376
1377   const int npar_osl=6;//5
1378   TMinuit MyMinuit_osl(npar_osl);
1379   MyMinuit_osl.SetFCN(fcnOSL);
1380   double OutputPar_osl[npar_osl]={0};
1381   double OutputPar_osl_e[npar_osl]={0};
1382
1383   double par_osl[npar_osl];               // the start values
1384   double stepSize_osl[npar_osl];          // step sizes
1385   double minVal_osl[npar_osl];            // minimum bound on parameter
1386   double maxVal_osl[npar_osl];            // maximum bound on parameter
1387   string parName_osl[npar_osl];
1388
1389   par_osl[0] = 1.0; par_osl[1]=0.5; par_osl[2]=0.; par_osl[3]=5.; par_osl[4] = 5.; par_osl[5]=5.;
1390   stepSize_osl[0] = 0.01; stepSize_osl[1] = 0.01; stepSize_osl[2] = 0.02; stepSize_osl[3] = 0.2; stepSize_osl[4] = 0.2; stepSize_osl[5] = 0.2;
1391   minVal_osl[0] = 0.8; minVal_osl[1] = 0.01; minVal_osl[2] = 0.0; minVal_osl[3] = 1.; minVal_osl[4] = 1.; minVal_osl[5] = 1.;
1392   maxVal_osl[0] = 1.2; maxVal_osl[1] = 1.0; maxVal_osl[2] = 0.99; maxVal_osl[3] = 20.; maxVal_osl[4] = 20.; maxVal_osl[5] = 20.;
1393   parName_osl[0] = "Norm"; parName_osl[1] = "Lambda"; parName_osl[2] = "G"; parName_osl[3] = "R_out"; parName_osl[4] = "R_side"; parName_osl[5] = "R_long";
1394
1395   for (int i=0; i<npar_osl; i++){
1396     MyMinuit_osl.DefineParameter(i, parName_osl[i].c_str(), par_osl[i], stepSize_osl[i], minVal_osl[i], maxVal_osl[i]);
1397   }
1398   MyMinuit_osl.FixParameter(2);// G
1399   //MyMinuit.FixParameter(1);// lambda
1400   cout<<"here!!"<<endl;
1401
1402   // Do the minimization!
1403   cout<<"Start Three-d fit"<<endl;
1404   MyMinuit_osl.Migrad();// Minuit's best minimization algorithm
1405   cout<<"End Three-d fit"<<endl;
1406   cout<<"Chi2/NDF = "<<Chi2_OSL/(NFitPoints_OSL-MyMinuit_osl.GetNumFreePars())<<endl;
1407
1408   for (int i=0; i<npar_osl; i++){
1409     MyMinuit_osl.GetParameter(i,OutputPar_osl[i],OutputPar_osl_e[i]);
1410   }
1411
1412
1413   TF3 *fit3d = new TF3("fit3d",OSLfitfunction, 0,0.2, 0,0.2, 0,0.2, npar_osl);
1414   for(int i=0; i<npar_osl; i++) {fit3d->FixParameter(i,OutputPar_osl[i]);}
1415   Int_t BinsOSL = num_osl->GetNbinsX();
1416   double LimitOSL = num_osl->GetXaxis()->GetBinUpEdge(BinsOSL);
1417   TH3D *den_timesFit = new TH3D("den_timesFit","",BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL);
1418   for(int i=0; i<BinsOSL; i++){
1419     for(int j=0; j<BinsOSL; j++){
1420       for(int k=0; k<BinsOSL; k++){
1421
1422         double qout=(i+.5)*0.005;
1423         double qside=(j+.5)*0.005;
1424         double qlong=(k+.5)*0.005;
1425         den_timesFit->SetBinContent(i+1,j+1,k+1, fit3d->Eval(qout,qside,qlong)*den_osl->GetBinContent(i+1,j+1,k+1));
1426         den_timesFit->SetBinError(i+1,j+1,k+1, 0);
1427       }
1428     }
1429   }
1430
1431   int binL=1, binH=4;
1432   TH1D *num_pro=(TH1D*)num_osl->ProjectionX("num_pro",binL,binH, binL,binH);
1433   TH1D *den_pro=(TH1D*)den_osl->ProjectionX("den_pro",binL,binH, binL,binH);
1434   TH1D *dentimesFit_pro=(TH1D*)den_timesFit->ProjectionX("dentimesFit_pro",binL,binH, binL,binH);
1435   num_pro->Sumw2();
1436   den_pro->Sumw2();
1437   num_pro->Divide(den_pro);
1438   num_pro->SetMarkerStyle(20);
1439   num_pro->SetTitle("");
1440   num_pro->GetXaxis()->SetTitle("q_{out}");
1441   num_pro->GetYaxis()->SetTitle("C_{2}");
1442   num_pro->SetMinimum(0.97);
1443   num_pro->SetMaximum(1.48);
1444   num_pro->DrawCopy();
1445   //
1446   dentimesFit_pro->Divide(den_pro);
1447   dentimesFit_pro->SetLineColor(2);
1448   dentimesFit_pro->DrawCopy("same");
1449   //
1450   */
1451   //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1452   //TGraph *gr2 = (TGraph*)MyMinuit.Contour(10,1,2);
1453   //gr2->SetFillColor(kYellow);
1454   //gr2->Draw("alf");
1455   //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1456   //MyMinuit.SetErrorDef(1);
1457   //TGraph *gr1 = (TGraph*)MyMinuit.Contour(10,1,2);
1458   //gr1->SetFillColor(kGreen);//38
1459   //gr1->Draw("lf");
1460
1461
1462
1463   //////////////////////////////////////////////////////////////////////////////////////
1464
1465
1466   // To visualize the Qcut and Norm Regions
1467   //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
1468   //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);
1469   //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);
1470   //for(int bin=1; bin<=20; bin++) QcutRegion->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1471   //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1472   //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1473   //QcutRegion->SetFillColor(4);
1474   //NormRegion1->SetFillColor(2);
1475   //NormRegion2->SetFillColor(3);
1476   //Two_ex[0][0][0][0]->Draw();
1477   //QcutRegion->Draw("same");
1478   //NormRegion1->Draw("same");
1479   //NormRegion2->Draw("same");
1480
1481
1482
1483
1484   TCanvas *can2 = new TCanvas("can2", "can2",800,0,900,900);//800,0,600,900
1485   can2->SetHighLightColor(2);
1486   gStyle->SetOptFit(0111);
1487   can2->SetFillColor(10);
1488   can2->SetBorderMode(0);
1489   can2->SetBorderSize(2);
1490   can2->SetGridx();
1491   can2->SetGridy();
1492   can2->SetFrameFillColor(0);
1493   can2->SetFrameBorderMode(0);
1494   can2->SetFrameBorderMode(0);
1495   can2->cd();
1496   gPad->SetLeftMargin(0.14);
1497   gPad->SetRightMargin(0.02);
1498
1499   TLegend *legend2 = new TLegend(.4,.67,1.,.87,NULL,"brNDC");
1500   legend2->SetBorderSize(1);
1501   legend2->SetTextSize(.06);// small .03; large .06
1502   legend2->SetFillColor(0);
1503
1504   /////////////////////////////////////////////
1505   TH3D *C3QS_3d = new TH3D("C3QS_3d","",20,0,.98, 20,0,.1, 20,0,.1);
1506   TH3D *Combinatorics_3d = new TH3D("Combinatorics_3d","",20,0,.1, 20,0,.1, 20,0,.1);
1507   //
1508   const int Q3BINS = 20;
1509   TH1D *num_withRS = new TH1D("num_withRS","",Q3BINS,0,0.2);
1510   TH1D *num_withGRS = new TH1D("num_withGRS","",Q3BINS,0,0.2);
1511   TH1D *num_withTM = new TH1D("num_withTM","",Q3BINS,0,0.2);
1512   TH1D *Cterm1 = new TH1D("Cterm1","",Q3BINS,0,0.2);
1513   TH1D *Cterm1_noMRC = new TH1D("Cterm1_noMRC","",Q3BINS,0,0.2);
1514   TH1D *Cterm2 = new TH1D("Cterm2","",Q3BINS,0,0.2);
1515   TH1D *Cterm3 = new TH1D("Cterm3","",Q3BINS,0,0.2);
1516   TH1D *Cterm4 = new TH1D("Cterm4","",Q3BINS,0,0.2);
1517   TH1D *num_QS = new TH1D("num_QS","",Q3BINS,0,0.2);
1518   TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,0.2);
1519   TH1D *Combinatorics_1d_noMRC = new TH1D("Combinatorics_1d_noMRC","",Q3BINS,0,0.2);
1520   TH1D *Coul_Riverside = new TH1D("Coul_Riverside","",Q3BINS,0,0.2);
1521   TH1D *Coul_GRiverside = new TH1D("Coul_GRiverside","",Q3BINS,0,0.2);
1522   TH1D *Coul_Omega0 = new TH1D("Coul_Omega0","",Q3BINS,0,0.2);
1523   TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,0.2);
1524   TH1D *c3_hist_STAR = new TH1D("c3_hist_STAR","",Q3BINS,0,0.2);
1525   //TProfile *MomRes_2d = new TProfile("MomRes_2d","",Q3BINS,0,0.2, 0,20.,"");
1526   TProfile *MomRes_3d_term1 = new TProfile("MomRes_3d_term1","",Q3BINS,0,0.2, 0,20.,"");
1527   TProfile *MomRes_3d_term2 = new TProfile("MomRes_3d_term2","",Q3BINS,0,0.2, 0,20.,"");
1528   TProfile *MomRes_3d_term5 = new TProfile("MomRes_3d_term5","",Q3BINS,0,0.2, 0,20.,"");
1529   TH1D *r3_num_Q3 = new TH1D("r3_num_Q3","",Q3BINS,0,0.2);
1530   TH1D *r3_den_Q3 = new TH1D("r3_den_Q3","",Q3BINS,0,0.2);
1531   r3_num_Q3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1532   r3_num_Q3->GetYaxis()->SetTitle("r_{3}");
1533   r3_num_Q3->GetYaxis()->SetTitleOffset(1.3);
1534   r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
1535   r3_num_Q3->SetMinimum(0);
1536   r3_num_Q3->SetMaximum(2.4);
1537   r3_den_Q3->SetLineColor(2);
1538   r3_num_Q3->SetMarkerStyle(20);
1539   Coul_Omega0->GetXaxis()->SetRangeUser(0,0.099);
1540   Coul_Omega0->GetXaxis()->SetLabelSize(0.04);
1541   Coul_Omega0->GetYaxis()->SetLabelSize(0.04);
1542   c3_hist_STAR->GetXaxis()->SetRangeUser(0,0.099);
1543   c3_hist_STAR->SetMinimum(0.8); c3_hist_STAR->SetMaximum(1.02);
1544   c3_hist_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1545   c3_hist_STAR->GetYaxis()->SetTitle("c_{3}*^{++-}");
1546   c3_hist_STAR->GetYaxis()->SetTitleOffset(1.6);
1547   if(SameCharge) {Cterm1_noMRC->Sumw2(); Combinatorics_1d_noMRC->Sumw2();}
1548   //
1549   double num_QS_e[Q3BINS]={0};
1550   double c3_e[Q3BINS]={0};
1551   double r3_e[Q3BINS]={0};
1552
1553   // CB=Charge Bin. 0 means pi-, 1 means pi+
1554   int CB1=0, CB2=0, CB3=0;
1555   int CP12=1, CP13=1, CP23=1;
1556   if(CHARGE==-1) {
1557     if(SameCharge) {CB1=0; CB2=0; CB3=0;}
1558     else {CB1=0; CB2=0; CB3=1; CP12=+1, CP13=-1, CP23=-1;}
1559   }else {
1560     if(SameCharge) {CB1=1; CB2=1; CB3=1;}
1561     else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;}
1562   }
1563
1564   cout<<"here"<<endl;
1565   // SC = species combination.  Always 0 meaning pi-pi-pi.
1566   int SCBin=0;
1567   //
1568   ReadCoulCorrections(RValue, bValue, 10);// switch to full kt range, 10.
1569   //ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation (STAR method testing)
1570   TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",20,0,0.2);
1571   TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",20,0,0.2);
1572   //
1573   double value_num[2]={0};
1574   double value_num_e[2]={0};
1575   double value_den[2]={0};
1576   double N3_QS[2]={0};
1577   double N3_QS_e[2]={0};
1578   double OutTriplets=0, InTriplets=0;
1579   for(int i=2; i<=20; i++){// bin number
1580     //double Q12_m = (i-0.5)*(0.005);// geometric center
1581     double Q12 = AvgQinvSS[i-1]; if(!SameCharge && CHARGE==+1) {Q12 = AvgQinvOS[i-1];}// true center
1582     int Q12bin = int(Q12/0.005)+1;
1583     //
1584     for(int j=2; j<=20; j++){// bin number
1585       //double Q13_m = (j-0.5)*(0.005);// geometric center
1586       double Q13 = AvgQinvSS[j-1]; if(!SameCharge) {Q13 = AvgQinvOS[j-1];}// true center
1587       int Q13bin = int(Q13/0.005)+1;
1588       //
1589       for(int k=2; k<=20; k++){// bin number
1590         //double Q23_m = (k-0.5)*(0.005);// geometric center
1591         double Q23 = AvgQinvSS[k-1]; if(!SameCharge && CHARGE==-1) {Q23 = AvgQinvOS[k-1];}// true center
1592         int Q23bin = int(Q23/0.005)+1;
1593         //
1594         //if(fabs(i-j)>3 || fabs(i-k)>3 || fabs(j-k)>3) continue;// testing
1595
1596         double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
1597         int Q3bin = Cterm1->GetXaxis()->FindBin(Q3);
1598
1599         if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
1600         if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
1601
1602
1603         // point-source Coulomb correlation
1604         double G_12 = Gamov(CP12, Q12);
1605         double G_13 = Gamov(CP13, Q13);
1606         double G_23 = Gamov(CP23, Q23);
1607         // finite-source Coulomb+Strong correlation from Therminator.
1608         double K2_12 = CoulCorr2(CP12, Q12);
1609         double K2_13 = CoulCorr2(CP13, Q13);
1610         double K2_23 = CoulCorr2(CP23, Q23);
1611         double K3 = 1.;// 3-body Coulomb+Strong correlation
1612         if(GRS) {// GRS approach
1613           if(SameCharge || CHARGE==-1) K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
1614           else K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
1615         }else {
1616           if(SameCharge || CHARGE==-1) K3 = CoulCorrOmega0(SameCharge, Q12, Q13, Q23);
1617           else K3 = CoulCorrOmega0(SameCharge, Q23, Q12, Q13);
1618         }
1619         if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;}
1620         if(K3==0) continue;
1621         if(GeneratedSignal) K3 = K2_12*K2_13*K2_23;// no interpolation for Generated signal
1622
1623         double TERM1=Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k);// N3 (3-pion yield per q12,q13,q23 cell). 3-pions from same-event
1624         double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
1625         double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
1626         double TERM4=Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event
1627         double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
1628         if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
1629         if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
1630
1631         double muonCorr_C3=1.0, muonCorr_C2_12=1.0, muonCorr_C2_13=1.0, muonCorr_C2_23=1.0;
1632         double muonCorr_W12=1.0, muonCorr_W13=1.0, muonCorr_W23=1.0;
1633         if(MuonCorrection){
1634           if(SameCharge){
1635             muonCorr_C3 = C3muonCorrection->GetBinContent(Q3bin);
1636             muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(Q12bin);
1637             muonCorr_C2_13 = C2muonCorrection_SC->GetBinContent(Q13bin);
1638             muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(Q23bin);
1639             muonCorr_W12 = WeightmuonCorrection->GetBinContent(Q12bin);
1640             muonCorr_W13 = WeightmuonCorrection->GetBinContent(Q13bin);
1641             muonCorr_W23 = WeightmuonCorrection->GetBinContent(Q23bin);
1642           }else{
1643             muonCorr_C3 = C3muonCorrection->GetBinContent(Q3bin);
1644             if(CHARGE==-1){// pi-pi-pi+
1645               muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(Q12bin);
1646               muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(Q13bin);
1647               muonCorr_C2_23 = C2muonCorrection_MC->GetBinContent(Q23bin);
1648             }else{// pi-pi+pi+
1649               muonCorr_C2_12 = C2muonCorrection_MC->GetBinContent(Q12bin);
1650               muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(Q13bin);
1651               muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(Q23bin);
1652             }
1653           }
1654         }
1655
1656         if(Q3>0.08 && Q3<0.09){// just for testing
1657           if(Q12>0.08 || Q13>0.08 || Q23>0.08) OutTriplets++;
1658           else InTriplets++;
1659         }
1660
1661         // apply momentum resolution correction
1662         if(!MCcase && !GeneratedSignal){
1663           if(SameCharge){
1664             // 3d momentum resolution corrections
1665             TERM1 *= (MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1666             TERM2 *= (MomRes3d[0][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1667             TERM3 *= (MomRes3d[0][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1668             TERM4 *= (MomRes3d[0][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1669             TERM5 *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1670             MomRes_3d_term1->Fill(Q3, MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM1);
1671             MomRes_3d_term2->Fill(Q3, MomRes3d[0][1]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM2);
1672             MomRes_3d_term5->Fill(Q3, MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM5);
1673           }else{
1674             if(CHARGE==-1){// pi-pi-pi+
1675               TERM1 *= (MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1676               TERM2 *= (MomRes3d[1][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1677               TERM3 *= (MomRes3d[1][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1678               TERM4 *= (MomRes3d[1][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1679               TERM5 *= (MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1680             }else {// pi-pi+pi+
1681               TERM1 *= (MomRes3d[1][0]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1682               TERM2 *= (MomRes3d[1][3]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1683               TERM3 *= (MomRes3d[1][2]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1684               TERM4 *= (MomRes3d[1][1]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1685               TERM5 *= (MomRes3d[1][4]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1686             }
1687             MomRes_3d_term1->Fill(Q3, MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin), TERM1);
1688             MomRes_3d_term2->Fill(Q3, MomRes3d[1][1]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM2);
1689             MomRes_3d_term5->Fill(Q3, MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin),TERM5);
1690           }
1691         }
1692         //
1693
1694         for(int LamType=0; LamType<2; LamType++){
1695
1696           if(LamType==1) TwoFrac=TwoFracr3;// r3
1697           else TwoFrac = 0.7;//OutputPar[1];// c3 (newly fixed to 0.7)
1698
1699           if(FileSetting==9) TwoFrac=0.8;// for FB5and7overlap choose a fixed higher lambda
1700
1701           OneFrac=pow(TwoFrac,0.5); // 0.495 to bring baseline up
1702           ThreeFrac=pow(TwoFrac,1.5);
1703
1704           // finite-multiplicity method
1705           //OneFrac = (1+sqrt(1 + 4*AvgN[Mbin]*TwoFrac*(AvgN[Mbin]-1)))/(2*AvgN[Mbin]);
1706           //ThreeFrac = (OneFrac*AvgN[Mbin]*(OneFrac*AvgN[Mbin]-1)*(OneFrac*AvgN[Mbin]-2))/(AvgN[Mbin]*(AvgN[Mbin]-1)*(AvgN[Mbin]-2));
1707
1708
1709           // Purify.  Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
1710           N3_QS[LamType] = TERM1;
1711           N3_QS[LamType] -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1712           N3_QS[LamType] -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
1713           N3_QS[LamType] /= ThreeFrac;
1714           N3_QS[LamType] /= K3;
1715           N3_QS[LamType] *= muonCorr_C3;
1716
1717           if(LamType==0) num_QS->Fill(Q3, N3_QS[LamType]);
1718
1719           // Isolate 3-pion cumulant
1720           value_num[LamType] = N3_QS[LamType];
1721           value_num[LamType] -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac * muonCorr_C2_12;
1722           value_num[LamType] -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac * muonCorr_C2_13;
1723           value_num[LamType] -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac * muonCorr_C2_23;
1724           value_num[LamType] += 2*TERM5;
1725
1726           //value_num[LamType] += 0.004*TERM5;
1727
1728           // r3 denominator
1729           if(LamType==1 && SameCharge) {
1730             value_den[LamType] = PionNorm[CB1]->GetBinContent(i,j,k);// Raw r3 denominator
1731             if(!MCcase){
1732               // Momentum Resolution Correction Systematic variation. Only important when MRCShift != 1.0.
1733               double denMRC1 = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1734               denMRC1 *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1735               denMRC1 *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1736               double denMRC2 = (C2ssRaw->GetBinContent(i)*((MomRes_C2_pp[i-1]-1)*MRCShift+1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1737               denMRC2 *= (C2ssRaw->GetBinContent(j)*((MomRes_C2_pp[j-1]-1)*MRCShift+1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1738               denMRC2 *= (C2ssRaw->GetBinContent(k)*((MomRes_C2_pp[k-1]-1)*MRCShift+1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1739               // Non-femto background correction (Mini-Jet).
1740               double denMJ = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - (MJ_bkg_ss->Eval(Q12)-1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1741               denMJ *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - (MJ_bkg_ss->Eval(Q13)-1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1742               denMJ *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - (MJ_bkg_ss->Eval(Q23)-1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1743               // Strong same-charge correction
1744               double Coul_12 = K2_12/StrongSC->Eval(Q12*1000/2.);// Coulomb used online
1745               double Coul_13 = K2_13/StrongSC->Eval(Q13*1000/2.);// Coulomb used online
1746               double Coul_23 = K2_23/StrongSC->Eval(Q23*1000/2.);// Coulomb used online
1747               double denCoulOnly = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*Coul_12 - (1-TwoFrac))/(TwoFrac*Coul_12);
1748               denCoulOnly *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*Coul_13 - (1-TwoFrac))/(TwoFrac*Coul_13);
1749               denCoulOnly *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*Coul_23 - (1-TwoFrac))/(TwoFrac*Coul_23);
1750               // K shift testing
1751               double K2back_12 = (K2_12-1)/KShift+1;
1752               double K2back_13 = (K2_13-1)/KShift+1;
1753               double K2back_23 = (K2_23-1)/KShift+1;
1754               double denKshiftBack = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*K2back_12 - (1-TwoFrac))/(TwoFrac*K2back_12);
1755               denKshiftBack *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*K2back_13 - (1-TwoFrac))/(TwoFrac*K2back_13);
1756               denKshiftBack *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*K2back_23 - (1-TwoFrac))/(TwoFrac*K2back_23);
1757
1758               double den_ratio = sqrt(fabs(denMRC2))/sqrt(fabs(denMRC1));
1759               value_den[LamType] *= den_ratio;// apply Momentum Resolution correction systematic variation if any.
1760               double den_ratioMJ = sqrt(fabs(denMJ))/sqrt(fabs(denMRC1));
1761               if(IncludeMJcorrection) value_den[LamType] *= den_ratioMJ;// Non-femto bkg correction
1762               if(IncludeStrongSC){
1763                 double den_ratioStrong = sqrt(fabs(denMRC1))/sqrt(fabs(denCoulOnly));
1764                 value_den[LamType] *= den_ratioStrong;
1765               }
1766               double den_ratioKShift = sqrt(fabs(denMRC1))/sqrt(fabs(denKshiftBack));
1767               value_den[LamType] *= den_ratioKShift;
1768               //
1769               value_den[LamType] *= sqrt(muonCorr_W12*muonCorr_W13*muonCorr_W23);// muon correction
1770               //
1771               //value_den[LamType] -= TPNRejects->GetBinContent(i,j,k);// Testing for 0-5% only to estimate the effect when C2^QS < 1.0
1772               value_den[LamType] *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;// momentum resolution correction to combinatorics
1773             }// !MCcase
1774           }
1775
1776
1777           // errors
1778           N3_QS_e[LamType] = fabs(TERM1);
1779           N3_QS_e[LamType] += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(fabs(TERM5)),2);
1780           N3_QS_e[LamType] += (pow((1-OneFrac),2)*fabs(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(fabs(TERM5)),2));
1781           N3_QS_e[LamType] /= pow(ThreeFrac,2);
1782           N3_QS_e[LamType] /= pow(K3,2);
1783           if(LamType==0) num_QS_e[Q3bin-1] += N3_QS_e[LamType];
1784           value_num_e[LamType] = N3_QS_e[LamType];
1785           value_num_e[LamType] += (pow(1/K2_12/TwoFrac*sqrt(fabs(TERM2)),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(fabs(TERM5)),2));
1786           value_num_e[LamType] += (pow(1/K2_13/TwoFrac*sqrt(fabs(TERM3)),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(fabs(TERM5)),2));
1787           value_num_e[LamType] += (pow(1/K2_23/TwoFrac*sqrt(fabs(TERM4)),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(fabs(TERM5)),2));
1788           value_num_e[LamType] += pow(2*sqrt(fabs(TERM5)),2);
1789           if(LamType==0) c3_e[Q3bin-1] += value_num_e[LamType] + TERM5;// add baseline stat error
1790           else r3_e[Q3bin-1] += value_num_e[LamType];
1791
1792         }
1793
1794         // Fill histograms
1795         c3_hist->Fill(Q3, value_num[0] + TERM5);// for cumulant correlation function
1796         C3QS_3d->SetBinContent(i,j,k, N3_QS[0]);
1797         C3QS_3d->SetBinError(i,j,k, N3_QS_e[0]);
1798         //
1799         Coul_Riverside->Fill(Q3, G_12*G_13*G_23*TERM5);
1800         Coul_GRiverside->Fill(Q3, TERM5*CoulCorrGRS(SameCharge, Q12, Q13, Q23));
1801         Coul_Omega0->Fill(Q3, K3*TERM5);
1802         num_withGRS->Fill(Q3, N3_QS[0]);
1803         Cterm1->Fill(Q3, TERM1);
1804         Cterm2->Fill(Q3, TERM2);
1805         Cterm3->Fill(Q3, TERM3);
1806         Cterm4->Fill(Q3, TERM4);
1807         Combinatorics_1d->Fill(Q3, TERM5);
1808         Combinatorics_3d->Fill(Q12,Q13,Q23, TERM5);
1809         r3_num_Q3->Fill(Q3, value_num[1]);
1810         r3_den_Q3->Fill(Q3, value_den[1]);
1811         double Q3_m = sqrt(pow((i-0.5)*(0.005),2) + pow((j-0.5)*(0.005),2) + pow((k-0.5)*(0.005),2));
1812         Cterm1_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k));
1813         Combinatorics_1d_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1814         double cumulant_STAR = Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k)/(K2_12*K2_13*K2_23);
1815         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k)/(K2_12);
1816         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k)/(K2_13);
1817         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k)/(K2_23);
1818         c3_hist_STAR->Fill(Q3_m, cumulant_STAR + 3*Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1819
1820         ///////////////////////////////////////////////////////////
1821         // Edgeworth 3-pion expection.  Not important for r3.
1822         //double radius_exp = (11-(MomResCentBin-1))/FmToGeV;
1823         //TwoFrac=0.68; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1824         double radius_exp = (OutputPar[3])/FmToGeV;
1825         TwoFrac=OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1826         double arg12 = Q12*radius_exp;
1827         double arg13 = Q13*radius_exp;
1828         double arg23 = Q23*radius_exp;
1829         Float_t EW12 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1830         EW12 += 0.45/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1831         Float_t EW13 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1832         EW13 += 0.45/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1833         Float_t EW23 = 1 + 0.2/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1834         EW23 += 0.45/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
1835
1836         //Float_t EW12 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1837         //EW12 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1838         //Float_t EW13 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1839         //EW13 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1840         //Float_t EW23 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1841         //EW23 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
1842         //
1843         double cumulant_exp=0, term1_exp=0;
1844         if(SameCharge) {
1845           cumulant_exp = (exp(-pow(radius_exp*Q12,2))*pow(EW12,2)+exp(-pow(radius_exp*Q13,2))*pow(EW13,2)+exp(-pow(radius_exp*Q23,2))*pow(EW23,2))*TERM5;
1846           cumulant_exp += 2*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5;
1847           term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1848           term1_exp += TwoFrac*(1-OneFrac)*(K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2))+K2_13*(1+exp(-pow(radius_exp*Q13,2))*pow(EW13,2))+K2_23*(1+exp(-pow(radius_exp*Q23,2))*pow(EW23,2)))*TERM5;
1849           term1_exp += ThreeFrac*K3*cumulant_exp;
1850           //term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1851         }else {
1852           cumulant_exp = (exp(-pow(radius_exp*Q12,2))*pow(EW12,2))*TERM5 + TERM5;
1853           term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1854           term1_exp += TwoFrac*(1-OneFrac)*(K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)) + K2_13 + K2_23)*TERM5;
1855           term1_exp += ThreeFrac*K3*cumulant_exp;
1856           term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1857           //term1_exp = ((1-TwoFrac) + TwoFrac*K2_13)*TERM5;
1858         }
1859
1860         GenSignalExpected_num->Fill(Q3, term1_exp);
1861         GenSignalExpected_den->Fill(Q3, TERM5);
1862         ///////////////////////////////////////////////////////////
1863
1864       }
1865     }
1866   }
1867
1868   cout<<"OutTriplets: "<<OutTriplets<<"   InTriplets: "<<InTriplets<<endl;
1869   ////////////////////////////
1870
1871   // Intermediate steps
1872   num_withRS->Sumw2();
1873   num_withGRS->Sumw2();
1874   num_withTM->Sumw2();
1875   Cterm1->Sumw2();
1876   Cterm2->Sumw2();
1877   Cterm3->Sumw2();
1878   Cterm4->Sumw2();
1879   Combinatorics_1d->Sumw2();
1880   Combinatorics_3d->Sumw2();
1881   r3_num_Q3->Sumw2();
1882   r3_den_Q3->Sumw2();
1883   for(int i=0; i<Q3BINS; i++) {c3_hist->SetBinError(i+1, sqrt(c3_e[i])); num_QS->SetBinError(i+1, sqrt(num_QS_e[i]));}
1884   for(int i=0; i<Q3BINS; i++) {c3_hist_STAR->SetBinError(i+1, sqrt(c3_e[i]));}// approximate error
1885   num_withRS->Divide(Combinatorics_1d);
1886   num_withGRS->Divide(Combinatorics_1d);
1887   num_withTM->Divide(Combinatorics_1d);
1888   for(int q3bin=1; q3bin<=r3_num_Q3->GetNbinsX(); q3bin++){
1889     r3_num_Q3->SetBinError(q3bin, sqrt(r3_e[q3bin-1]));
1890   }
1891   Cterm1->Divide(Combinatorics_1d);
1892   Cterm1_noMRC->Divide(Combinatorics_1d_noMRC);
1893   Cterm2->Divide(Combinatorics_1d);
1894   Cterm3->Divide(Combinatorics_1d);
1895   Cterm4->Divide(Combinatorics_1d);
1896   c3_hist->Divide(Combinatorics_1d);
1897   c3_hist_STAR->Divide(Combinatorics_1d_noMRC);
1898   num_QS->Divide(Combinatorics_1d);
1899   r3_num_Q3->Divide(r3_den_Q3);
1900   GenSignalExpected_num->Sumw2();
1901   GenSignalExpected_num->Divide(GenSignalExpected_den);
1902
1903   //
1904   //
1905   Coul_Riverside->Divide(Combinatorics_1d);
1906   Coul_GRiverside->Divide(Combinatorics_1d);
1907   Coul_Omega0->Divide(Combinatorics_1d);
1908   for(int ii=1; ii<=Coul_Omega0->GetNbinsX(); ii++){
1909     Coul_Riverside->SetBinError(ii,0.000001);
1910     Coul_GRiverside->SetBinError(ii,0.000001);
1911     Coul_Omega0->SetBinError(ii,0.000001);
1912   }
1913   ////////////////////////////
1914
1915
1916
1917   ///////////////////////////////////////////////////////////////////////////////////////////////////
1918
1919   Coul_Riverside->SetLineColor(1);
1920   Coul_GRiverside->SetLineColor(2);
1921   Coul_Omega0->SetLineColor(4);
1922
1923   //Coul_Riverside->Draw();
1924   //legend1->AddEntry(Coul_Riverside,"Riverside","l");
1925   //Coul_GRiverside->Draw("same");
1926   //legend1->AddEntry(Coul_GRiverside,"Generalized Riverside","l");
1927   //Coul_Omega0->Draw("same");
1928   //legend1->AddEntry(Coul_Omega0,"Omega0","l");
1929
1930
1931   Cterm1->SetMarkerStyle(20);
1932   Cterm1->SetMarkerColor(4);
1933   Cterm1->SetLineColor(4);
1934   Cterm1->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1935   Cterm1->GetYaxis()->SetTitle("C_{3}");
1936   //Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}");
1937   Cterm1->SetMinimum(0.99);
1938   Cterm1->SetMaximum(1.08);// 6.1 for same-charge
1939   Cterm1->GetXaxis()->SetRangeUser(0,.095);
1940   Cterm1->GetYaxis()->SetTitleOffset(1.4);
1941   Cterm1->Draw();
1942   legend2->AddEntry(Cterm1,"C_{3} raw","p");
1943   //
1944   Cterm2->SetMarkerStyle(20);
1945   Cterm2->SetMarkerColor(7);
1946
1947   Cterm3->SetMarkerStyle(25);
1948   Cterm3->SetMarkerColor(8);
1949   Cterm4->SetMarkerStyle(26);
1950   Cterm4->SetMarkerColor(9);
1951   //Cterm2->Draw("same");
1952   //Cterm3->Draw("same");
1953   //Cterm4->Draw("same");
1954   //legend2->AddEntry(Cterm1,"++-","p");
1955
1956
1957   if(MCcase){
1958     double C3points_HIJING_mmm[10]={0, 0.961834, 1.01827, 0.999387, 1.00202, 1.00081, 1.00082, 1.00037, 0.999074, 0.999099};
1959     double C3points_HIJING_mmm_e[10]={0, 0.0833895, 0.015158, 0.0047978, 0.00235067, 0.00138155, 0.00087485, 0.000612203, 0.000450162, 0.000346943};
1960     double C3points_HIJING_ppp[10]={0, 1.13015, 1.00623, 1.00536, 1.00293, 0.999964, 1.00116, 1.0007, 1.00037, 1.00105};
1961     double C3points_HIJING_ppp_e[10]={0, 0.0977058, 0.0150694, 0.0048196, 0.00235518, 0.00138376, 0.000877562, 0.000614069, 0.000452051, 0.00034856};
1962     TH1D *C3_HIJING_mmm=(TH1D*)Cterm1->Clone();
1963     TH1D *C3_HIJING_ppp=(TH1D*)Cterm1->Clone();
1964     for(int i=0; i<10; i++){
1965       C3_HIJING_mmm->SetBinContent(i+1, C3points_HIJING_mmm[i]);
1966       C3_HIJING_mmm->SetBinError(i+1, C3points_HIJING_mmm_e[i]);
1967       C3_HIJING_ppp->SetBinContent(i+1, C3points_HIJING_ppp[i]);
1968       C3_HIJING_ppp->SetBinError(i+1, C3points_HIJING_ppp_e[i]);
1969     }
1970     C3_HIJING_mmm->SetMarkerColor(2);
1971     C3_HIJING_mmm->SetLineColor(2);
1972     C3_HIJING_ppp->SetMarkerColor(4);
1973     C3_HIJING_ppp->SetLineColor(4);
1974     //C3_HIJING_mmm->Draw("same");
1975     //C3_HIJING_ppp->Draw("same");
1976     //legend2->AddEntry(C3_HIJING_mmm,"---","p");
1977     //legend2->AddEntry(C3_HIJING_ppp,"+++","p");
1978   }
1979
1980   num_QS->SetMarkerStyle(24);
1981   num_QS->SetMarkerColor(4);
1982   num_QS->SetLineColor(4);
1983   num_QS->GetXaxis()->SetTitle("Q_{3}");
1984   num_QS->GetYaxis()->SetTitle("C_{3}^{QS}");
1985   num_QS->GetXaxis()->SetRangeUser(0,.12);
1986   num_QS->SetMaximum(6);
1987   //num_QS->Draw("same");
1988   //legend2->AddEntry(num_QS,"C_{3}^{QS}","p");
1989
1990   c3_hist->GetXaxis()->SetTitle("Q_{3}");
1991   c3_hist->GetYaxis()->SetTitle("c_{3}^{++-}");
1992   c3_hist->GetXaxis()->SetRangeUser(0,.15);
1993   c3_hist->SetMarkerStyle(25);
1994   c3_hist->SetMarkerColor(2);
1995   c3_hist->SetLineColor(2);
1996   c3_hist->SetMaximum(3);
1997   c3_hist->SetMinimum(.9);
1998   if(!MCcase && !GeneratedSignal) c3_hist->Draw("same");
1999   //legend2->AddEntry(c3_hist,"#font[12]{c}_{3} (cumulant correlation)","p");
2000
2001   if(SameCharge){
2002     //for(int ii=0; ii<10; ii++) cout<<c3_hist->GetBinContent(ii+1)<<", ";
2003     TF1 *c3Fit=new TF1("c3Fit","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,0.2);
2004     //TF1 *c3Fit=new TF1("c3Fit","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + (-0.12/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + (0.43/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1);
2005     c3Fit->SetParameter(0,1);
2006     c3Fit->SetParameter(1,2);
2007     c3Fit->SetParameter(2,6);
2008     //c3Fit->FixParameter(1,0.8545);
2009     //c3Fit->FixParameter(2,8.982);
2010     //c3_hist->Fit(c3Fit,"IME","",0,0.11);
2011     //cout<<"Chi2/NDF = "<<c3Fit->GetChisquare()/c3Fit->GetNDF()<<endl;
2012   }
2013
2014   GenSignalExpected_num->SetMarkerStyle(20);
2015   //GenSignalExpected_num->Draw("same");
2016   //legend2->AddEntry(GenSignalExpected_num,"#kappa_{3}=0.24, #kappa_{4}=0.16, #lambda=0.68, R=6 fm","p");
2017   //legend2->AddEntry(GenSignalExpected_num,"Edeworth expectation (fully chaotic)","p");
2018
2019   //MomRes_2d->SetMarkerStyle(20);
2020   //MomRes_3d->SetMarkerStyle(20);
2021   //MomRes_3d->SetMarkerColor(4);
2022   //MomRes_2d->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
2023   //MomRes_2d->GetYaxis()->SetTitle("< MRC >");
2024   //MomRes_2d->SetTitle("");
2025   //MomRes_2d->Draw();
2026   //legend2->AddEntry(MomRes_2d,"2D: Triple pair product","p");
2027   MomRes_3d_term2->SetLineColor(2);
2028   MomRes_3d_term5->SetLineColor(4);
2029   //MomRes_3d_term1->Draw();
2030   //MomRes_3d_term2->Draw("same");
2031   //MomRes_3d_term5->Draw("same");
2032   //legend2->AddEntry(MomRes_3d,"3D","p");
2033
2034   //legend2->Draw("same");
2035
2036   //cout<<c3_hist->Integral(8,10)/3.<<endl;
2037
2038   /*
2039   ////////// C3 systematics
2040   // C3 --+ base, M0, Kt3 0
2041   double C3_base[10]={0, 1.64715, 1.40709, 1.24344, 1.15809, 1.1071, 1.07544, 1.0534, 1.03881, 1.02974};
2042   double C3_base_e[10]={0, 0.00348782, 0.000841611, 0.00031699, 0.000163779, 9.95652e-05, 6.43811e-05, 4.60347e-05, 3.45978e-05, 2.68396e-05};
2043   //
2044   // HIJING C3 --- base, M0
2045   //double C3_base[10]={0, 0.970005, 1.00655, 1.00352, 1.00291, 1.00071, 1.0002, 0.999524, 0.999404, 0.999397};
2046   //double C3_base_e[10]={0, 0.050845, 0.0099602, 0.00334862, 0.00138008, 0.000841743, 0.000531776, 0.000371712, 0.000274716, 0.00021};
2047   // HIJING C3 --- base, M1
2048   //double C3_base[10]={0, 1.03657, 1.00199, 0.997984, 1.00067, 1.0006, 0.999901, 0.999967, 0.999792, 0.999777};
2049   //double C3_base_e[10]={0, 0.0634232, 0.0117204, 0.0039446, 0.00163131, 0.000996638, 0.000629662, 0.000440266, 0.00032534, 0.000249};
2050   // HIJING C3 --- base, M2
2051   //double C3_base[10]={0, 1.34345, 1.04226, 1.0278, 0.99582, 1.00554, 1.00296, 1.00057, 1.00271, 1.00152};
2052   //double C3_base_e[10]={0, 0.363559, 0.0503531, 0.0170117, 0.00679185, 0.00419035, 0.00264603, 0.00184587, 0.00136663, 0.00104772};
2053   // HIJING C3 --- base, M3
2054   //double C3_base[10]={0, 0.890897, 1.05222, 1.00461, 1.01364, 0.998981, 1.00225, 1.00305, 1.00235, 1.00043};
2055   //double C3_base_e[10]={0, 0.297124, 0.0604397, 0.0195066, 0.00812906, 0.00490835, 0.00310751, 0.00217408, 0.00160575, 0.00123065};
2056   TH1D *C3baseHisto=(TH1D*)Cterm1->Clone();
2057   for(int ii=0; ii<10; ii++){
2058     C3baseHisto->SetBinContent(ii+1, C3_base[ii]);
2059     C3baseHisto->SetBinError(ii+1, C3_base_e[ii]);
2060   }
2061
2062   cout<<"C3 values:"<<endl;
2063   for(int ii=0; ii<10; ii++){
2064     cout<<Cterm1->GetBinContent(ii+1)<<", ";
2065   }
2066   cout<<endl;
2067   cout<<"C3 errors:"<<endl;
2068   for(int ii=0; ii<10; ii++){
2069     cout<<Cterm1->GetBinError(ii+1)<<", ";
2070   }
2071   cout<<endl;
2072   TH1D *C3_Sys=(TH1D*)Cterm1->Clone();
2073   for(int ii=1; ii<10; ii++){
2074     if(C3_base[ii] ==0) {
2075       C3_Sys->SetBinContent(ii+1, 100.);
2076       C3_Sys->SetBinError(ii+1, 100.);
2077       continue;
2078     }
2079     C3_Sys->SetBinContent(ii+1, (C3_base[ii]-C3_Sys->GetBinContent(ii+1))/C3_base[ii]);
2080     C3_Sys->SetBinError(ii+1, sqrt(fabs(pow(C3_Sys->GetBinError(ii+1),2) - pow(C3_base_e[ii],2))));
2081   }
2082   gStyle->SetOptFit(111);
2083   C3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2084   C3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2085   C3_Sys->GetYaxis()->SetTitleOffset(2);
2086   C3_Sys->SetMinimum(-0.01);
2087   C3_Sys->SetMaximum(0.01);
2088   //C3_Sys->Draw();
2089   TF1 *C3lineSys=new TF1("C3lineSys","pol3",0,0.099);
2090   C3_Sys->Fit(C3lineSys,"MEQ","",0,0.099);
2091
2092   if(MCcase){
2093     // Plotting +++ added with --- for HIJING
2094     TLegend *legendC3Hijing = new TLegend(.5,.15,.9,.3,NULL,"brNDC");
2095     legendC3Hijing->SetBorderSize(1);
2096     legendC3Hijing->SetTextSize(.03);// small .03; large .06
2097     legendC3Hijing->SetFillColor(0);
2098
2099     Cterm1->Add(C3baseHisto);
2100     Cterm1->Scale(0.5);
2101     Cterm1->SetMinimum(0.9);
2102     Cterm1->SetMaximum(1.1);
2103     Cterm1->Draw();
2104     legendC3Hijing->AddEntry(Cterm1,"same-charge, HIJING","p");
2105     legendC3Hijing->Draw("same");
2106   }
2107   */
2108   /*
2109   ////////// c3 systematics
2110   // c3 --- base, M0, (0.04 TTC )
2111   //double c3_base[10]={0, 1.86014, 1.47533, 1.23733, 1.09944, 1.04145, 1.01693, 1.00715, 1.00253, 1.00111};
2112   //double c3_base_e[10]={0, 0.104645, 0.0120917, 0.00333303, 0.00118126, 0.0006692, 0.000405246, 0.000274163, 0.000198507, 0.000150258};
2113   // c3 --- base, M0, Kt3 0 (New Standard)
2114   double c3_base[10]={0, 1.00906, 1.00013, 1.00203, 1.00001, 1.00017, 1.00001, 0.999721, 0.999778, 0.999844};
2115   double c3_base_e[10]={0, 0.00726286, 0.00196376, 0.00081047, 0.00044086, 0.000276571, 0.000182574, 0.000132467, 0.000100557, 7.85009e-05};
2116
2117   cout<<"c3 values:"<<endl;
2118   for(int ii=0; ii<10; ii++){
2119     cout<<c3_hist->GetBinContent(ii+1)<<", ";
2120   }
2121   cout<<endl;
2122   cout<<"c3 errors:"<<endl;
2123   for(int ii=0; ii<10; ii++){
2124     cout<<c3_hist->GetBinError(ii+1)<<", ";
2125   }
2126   cout<<endl;
2127   TH1D *c3_Sys=(TH1D*)c3_hist->Clone();
2128   for(int ii=1; ii<10; ii++){
2129     if(c3_base[ii] ==0) {
2130       c3_Sys->SetBinContent(ii+1, 100.);
2131       c3_Sys->SetBinError(ii+1, 100.);
2132       continue;
2133     }
2134     c3_Sys->SetBinContent(ii+1, (c3_base[ii]-c3_Sys->GetBinContent(ii+1))/c3_base[ii]);
2135     c3_Sys->SetBinError(ii+1, sqrt(fabs(pow(c3_Sys->GetBinError(ii+1),2) - pow(c3_base_e[ii],2))));
2136   }
2137   gStyle->SetOptFit(111);
2138   c3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2139   c3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2140   c3_Sys->GetYaxis()->SetTitleOffset(2);
2141   c3_Sys->SetMinimum(-0.01);
2142   c3_Sys->SetMaximum(0.01);
2143   c3_Sys->Draw();
2144   TF1 *c3lineSys=new TF1("c3lineSys","pol3",0,0.099);
2145   c3_Sys->Fit(c3lineSys,"MEQ","",0,0.099);
2146   */
2147
2148
2149   /*TPad *pad1 = new TPad("pad1","pad1",0.0,0.0,1.,1.);
2150   gPad->SetGridx(0);
2151   gPad->SetGridy(0);
2152   gPad->SetTickx();
2153   gPad->SetTicky();
2154   pad1->SetTopMargin(0.02);//0.05
2155   //pad1->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2156   pad1->SetRightMargin(0.01);//1e-2
2157   pad1->SetBottomMargin(0.06);//0.12
2158   pad1->Divide(1,2,0,0);
2159   pad1->Draw();
2160   pad1->cd(1);
2161   gPad->SetLeftMargin(0.14);
2162   gPad->SetRightMargin(0.02);
2163   if(SameCharge){
2164     Coul_Omega0->SetMinimum(0.48);
2165     Coul_Omega0->SetMaximum(1.01);
2166   }else{
2167     Coul_Omega0->SetMinimum(0.99);
2168     Coul_Omega0->SetMaximum(1.32);
2169   }
2170   Coul_Omega0->SetMarkerStyle(20);
2171   Coul_Omega0->SetMarkerColor(4);
2172   Coul_GRiverside->SetMarkerStyle(25);
2173   Coul_GRiverside->SetMarkerColor(2);
2174   Coul_Riverside->SetMarkerStyle(23);
2175   Coul_Omega0->Draw();
2176   //Coul_Riverside->Draw("same");
2177   Coul_GRiverside->Draw("same");
2178   legend2->AddEntry(Coul_Omega0,"#Omega_{0}","p");
2179   legend2->AddEntry(Coul_GRiverside,"Generalized Riverside","p");
2180   //legend2->AddEntry(Coul_Riverside,"Riverside","p");
2181   legend2->Draw("same");
2182   TLatex *K3Label = new TLatex(-0.011,0.92,"K_{3}");// -0.011,0.92 (ss), 1.26 (os)
2183   K3Label->SetTextSize(0.08);
2184   K3Label->SetTextAngle(90);
2185   K3Label->Draw();
2186   //
2187   pad1->cd(2);
2188   gPad->SetLeftMargin(0.14);
2189   gPad->SetRightMargin(0.02);
2190   gPad->SetBottomMargin(0.13);
2191   TH1D *K3_compOmega0 = (TH1D*)Coul_Omega0->Clone();
2192   TH1D *K3_compGRS = (TH1D*)Coul_GRiverside->Clone();
2193   for(int ii=0; ii<K3_compOmega0->GetNbinsX(); ii++){
2194     K3_compOmega0->SetBinContent(ii+1, K3_compOmega0->GetBinContent(ii+1)-1.0);
2195     K3_compGRS->SetBinContent(ii+1, K3_compGRS->GetBinContent(ii+1)-1.0);
2196   }
2197   K3_compOmega0->Add(K3_compGRS,-1);
2198   K3_compOmega0->Divide(K3_compGRS);
2199   K3_compOmega0->SetMarkerStyle(20);
2200   K3_compOmega0->SetMarkerColor(1);
2201   K3_compOmega0->SetLineColor(1);
2202   K3_compOmega0->SetMinimum(-0.12);// -0.021
2203   K3_compOmega0->SetMaximum(0.12);// 0.021
2204   K3_compOmega0->SetBinContent(1,-100);
2205   K3_compOmega0->Draw();
2206   TLatex *RatioLabel = new TLatex(-.011,-0.05,"#Delta (K_{3}-1)/(K_{3}(#Omega_{0})-1)");// -0.011,0.04
2207   RatioLabel->SetTextSize(0.08);
2208   RatioLabel->SetTextAngle(90);
2209   RatioLabel->Draw();
2210   TLatex *Q3Label = new TLatex(.065,-0.147,"Q_{3} (GeV/c)");//.065,-0.147
2211   Q3Label->SetTextSize(0.08);
2212   Q3Label->Draw();
2213   */
2214
2215
2216   /////////////////////////////
2217   // 3d C3QS fitting
2218   /*
2219   for(int i=0; i<BINRANGE_C2global; i++){
2220     for(int j=0; j<BINRANGE_C2global; j++){
2221       for(int k=0; k<BINRANGE_C2global; k++){
2222         C3ssFitting[i][j][k] = C3QS_3d->GetBinContent(i+1,j+1,k+1);
2223         C3ssFitting_e[i][j][k] = C3QS_3d->GetBinError(i+1,j+1,k+1);
2224       }
2225     }
2226   }
2227
2228   TF3 *C3ssFitResult=new TF3("C3ssFitResult",C3ssFitFunction,0,0.2,0,0.2,0,0.2,npar);
2229   for(int i=0; i<npar; i++) {
2230     if(i<=5) C3ssFitResult->FixParameter(i, OutputPar[i]);
2231     else C3ssFitResult->FixParameter(i, OutputPar[i]-OutputPar_e[i]);
2232   }
2233
2234
2235   const int NbinsC3_1d=Q3BINS;
2236   TH1D *C3fit_1d=new TH1D("C3fit_1d","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2237   TH1D *C3fit_1d_den=new TH1D("C3fit_1d_den","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2238   double C3_err[NbinsC3_1d]={0};
2239
2240
2241   for(int i=0; i<C3QS_3d->GetNbinsX(); i++){// bin number
2242     for(int j=0; j<C3QS_3d->GetNbinsY(); j++){// bin number
2243       for(int k=0; k<C3QS_3d->GetNbinsZ(); k++){// bin number
2244
2245         if(i==0 || j==0 || k==0) continue;
2246         //if(i==1 || j==1 || k==1) continue;
2247         double Q3 = sqrt(pow(AvgQinvSS[i],2) + pow(AvgQinvSS[j],2) + pow(AvgQinvSS[k],2));
2248         C3fit_1d->Fill(Q3, C3ssFitResult->Eval(AvgQinvSS[i],AvgQinvSS[j],AvgQinvSS[k])*Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2249         C3fit_1d_den->Fill(Q3, Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2250       }
2251     }
2252   }
2253
2254   C3fit_1d->Divide(C3fit_1d_den);
2255   for(int q3bin=1; q3bin<=num_QS->GetNbinsX(); q3bin++) C3fit_1d->SetBinError(q3bin, 0.001);// non-zero to display marker
2256   C3fit_1d->SetMarkerStyle(22);
2257   C3fit_1d->SetLineColor(4);
2258   //C3fit_1d->Draw("same");
2259
2260   TF1 *lineC3 = new TF1("lineC3","6",0,10);
2261   lineC3->SetLineColor(4);
2262   lineC3->SetLineStyle(2);
2263   //lineC3->Draw("same");
2264   //legend2->AddEntry(lineC3,"Chaotic Limit C_{3}^{QS}","l");
2265
2266   TF1 *linec3 = new TF1("linec3","3",0,10);
2267   linec3->SetLineColor(2);
2268   linec3->SetLineStyle(2);
2269   //linec3->Draw("same");
2270   //legend2->AddEntry(linec3,"Chaotic Limit #font[12]{c}_{3}","l");
2271   */
2272
2273   //legend2->Draw("same");
2274
2275
2276   /////////////////////////////////////////////////////////////////
2277   /////////////////////////////////////////////////////////////////
2278
2279   TCanvas *can3 = new TCanvas("can3", "can3",1600,0,700,700);
2280   can3->SetHighLightColor(2);
2281   can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
2282   //gStyle->SetOptFit(0111);
2283   can3->SetFillColor(10);
2284   can3->SetBorderMode(0);
2285   can3->SetBorderSize(2);
2286   can3->SetGridx();
2287   can3->SetGridy();
2288   can3->SetFrameFillColor(0);
2289   can3->SetFrameBorderMode(0);
2290   can3->SetFrameBorderMode(0);
2291   can3->cd();
2292   TPad *pad3 = new TPad("pad3","pad3",0.0,0.0,1.,1.);
2293   gPad->SetGridx(1);
2294   gPad->SetGridy(1);
2295   gPad->SetTickx();
2296   gPad->SetTicky();
2297   pad3->SetTopMargin(0.02);//0.05
2298   //pad3->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2299   pad3->SetRightMargin(0.01);//1e-2
2300   pad3->SetBottomMargin(0.06);//0.12
2301   pad3->Draw();
2302   pad3->cd(1);
2303   gPad->SetLeftMargin(0.14);
2304   gPad->SetRightMargin(0.02);
2305
2306   TF1 *QuarticFit = new TF1("QuarticFit","[0]*(1-[1]*pow(x,4))",0,.1);
2307   QuarticFit->SetParameter(0,2); QuarticFit->SetParameter(1,0);
2308   QuarticFit->SetLineColor(4);
2309   QuarticFit->SetParName(0,"I^{Quartic}"); QuarticFit->SetParName(1,"a^{Quartic}");
2310   TF1 *QuadraticFit = new TF1("QuadraticFit","[0]*(1-[1]*pow(x,2))",0,.1);
2311   QuadraticFit->SetParameter(0,2); QuadraticFit->SetParameter(1,0);
2312   QuadraticFit->SetParName(0,"I^{Quadratic}"); QuadraticFit->SetParName(1,"a^{Quadratic}");
2313   TLegend *legend3 = new TLegend(.2,.85,.5,.95,NULL,"brNDC");
2314   legend3->SetBorderSize(1);
2315   legend3->SetTextSize(.04);// small .03; large .06
2316   legend3->SetFillColor(0);
2317
2318
2319   if(SameCharge){
2320
2321     r3_num_Q3->SetMinimum(0); r3_num_Q3->SetMaximum(2.5);// 0 to 2.5
2322     r3_num_Q3->GetXaxis()->SetRangeUser(0.0,0.1);
2323     r3_num_Q3->Draw();
2324
2325
2326     // HIJING standalone
2327     // Kt3 1
2328     /*double Gen_r3_m_M0[10]={0, 1.97963, 2.10248, 2.04465, 1.96697, 2.02295, 1.92281, 2.10031, 2.22338, 2.49729};
2329     double Gen_r3_m_M0_e[10]={0, 0.132726, 0.0668652, 0.0451347, 0.042838, 0.0546967, 0.0754362, 0.133624, 0.286307, 0.628381};
2330     double Gen_r3_p_M0[10]={0, 1.91771, 1.9653, 2.00742, 2.02393, 1.90624, 1.93554, 1.66, 1.79584, 0.301761};
2331     double Gen_r3_p_M0_e[10]={0, 0.133654, 0.0667213, 0.0451512, 0.0428925, 0.0547591, 0.0754764, 0.13365, 0.28638, 0.628441};*/
2332     // Kt3 2
2333     /*double Gen_r3_m_M0[10]={0, 1.12993, 2.09715, 1.91886, 2.08493, 2.10931, 2.00286, 1.99898, 1.78549, 1.91861};
2334     double Gen_r3_m_M0_e[10]={0, 0.237903, 0.115454, 0.0725178, 0.0630867, 0.0730326, 0.0886163, 0.12885, 0.213654, 0.379273};
2335     double Gen_r3_p_M0[10]={0, 2.05766, 1.97408, 2.05182, 2.02431, 2.11783, 1.93294, 1.97525, 2.21833, 2.31318};
2336     double Gen_r3_p_M0_e[10]={0, 0.24238, 0.11515, 0.0725077, 0.0629938, 0.0729842, 0.0886386, 0.12884, 0.213737, 0.379196};*/
2337
2338
2339     // HIJING + ALICE
2340     // Kt3 1
2341     /*double Gen_r3_m_M0[10]={0, 2.03468, 2.05783, 1.97757, 2.03809, 1.95703, 2.02915, 1.80055, 1.97664, 1.49573};
2342     double Gen_r3_m_M0_e[10]={0, 0.284164, 0.0923288, 0.0543837, 0.045952, 0.0565222, 0.0748221, 0.128994, 0.271954, 0.599077};
2343     double Gen_r3_p_M0[10]={0, 1.74611, 1.92369, 2.11024, 2.02823, 2.06235, 2.00127, 1.91551, 1.90576, 2.521};
2344     double Gen_r3_p_M0_e[10]={0, 0.279962, 0.0928866, 0.0548168, 0.0462253, 0.0569129, 0.0752683, 0.129803, 0.2736, 0.602535};
2345     double Gen_r3_m_M1[10]={0, 1.71341, 1.8735, 1.9352, 1.96466, 1.93852, 1.69509, 1.51402, 1.6014, -0.802941};
2346     double Gen_r3_m_M1_e[10]={0, 0.334562, 0.108368, 0.0640739, 0.0540579, 0.0663174, 0.0876532, 0.149302, 0.307982, 0.640548};
2347     double Gen_r3_p_M1[10]={0, 1.50861, 2.09508, 2.05338, 1.96178, 1.97999, 2.07162, 1.8607, 1.91805, 1.10468};
2348     double Gen_r3_p_M1_e[10]={0, 0.3399, 0.109449, 0.0646241, 0.0544419, 0.0667711, 0.0882695, 0.150195, 0.309844, 0.643923};*/
2349     // Kt3 2
2350     /*double Gen_r3_m_M0[10]={0, 1.04272, 2.05961, 1.93025, 2.08203, 2.07565, 2.29192, 1.93009, 2.68715, 1.71175};
2351     double Gen_r3_m_M0_e[10]={0, 3.09343, 0.30016, 0.117634, 0.0757208, 0.0794904, 0.090823, 0.128235, 0.21308, 0.40023};
2352     double Gen_r3_p_M0[10]={0, 1.83276, 1.88969, 2.03778, 2.0907, 2.11919, 2.04992, 2.02593, 1.95209, 1.68264};
2353     double Gen_r3_p_M0_e[10]={0, 3.59, 0.298056, 0.118354, 0.0762849, 0.079909, 0.0912179, 0.128866, 0.213999, 0.402173};
2354     double Gen_r3_m_M1[10]={0, 5.40628, 1.52822, 1.93258, 2.13338, 2.05811, 2.02963, 2.1204, 2.04906, 1.9021};
2355     double Gen_r3_m_M1_e[10]={0, 4.3025, 0.350163, 0.13871, 0.0897272, 0.0938973, 0.107045, 0.151557, 0.25029, 0.446325};
2356     double Gen_r3_p_M1[10]={0, 3.15883, 1.86772, 2.24914, 2.12118, 2.09175, 2.01447, 2.36802, 2.57239, 2.68729};
2357     double Gen_r3_p_M1_e[10]={0, 4.6622, 0.348237, 0.140294, 0.0903003, 0.09446, 0.107692, 0.152575, 0.251939, 0.448919};*/
2358
2359     //double r3Sys_M1[10]={0, 0.097, 0.056, 0.063, 0.097, 0.17, 0.32, 0.66, 1.4, 3.4};
2360     // muon variation
2361     /*double r3_muonVar_M1[10]={0, 1.73244, 1.80921, 1.77852, 1.7192, 1.62059, 1.50122, 1.37656, 1.01344, 0.755781};
2362     double r3_muonVar_M1_e[10]={0, 0.160786, 0.051075, 0.031982, 0.0293467, 0.0390241, 0.0514494, 0.0760959, 0.112944, 0.154592};
2363     TH1D *r3_muonVar=(TH1D*)r3_num_Q3->Clone();
2364     for(int i=0; i<10; i++){
2365       r3_muonVar->SetBinContent(i+1,r3_muonVar_M1[i]);
2366       r3_muonVar->SetBinError(i+1,sqrt(pow(r3_muonVar_M1_e[i],2)+pow(r3Sys_M1[i],2)));
2367     }
2368     r3_muonVar->SetMarkerStyle(21);
2369     r3_muonVar->SetMarkerColor(2); r3_muonVar->SetLineColor(2);
2370     r3_muonVar->SetFillStyle(1000); r3_muonVar->SetFillColor(kRed-10);
2371     r3_muonVar->Draw("E2 same");
2372     r3_num_Q3->Draw("same");
2373     legend3->AddEntry(r3_num_Q3,"-2.0<N#sigma_{Pion}<2.0 (default)","p");
2374     legend3->AddEntry(r3_muonVar,"-0.5<N#sigma_{Pion}<2.0","p");
2375     legend3->Draw("same");*/
2376
2377     // muon correction
2378     /*double r3_muonCorr_M1[10]={0, 1.8462, 1.84371, 1.79934, 1.77217, 1.76725, 1.79545, 1.7986, 2.11717, 2.86177};
2379     double r3_muonCorr_M1_e[10]={0, 0.0838277, 0.0266269, 0.0168108, 0.0156166, 0.0211333, 0.0286836, 0.0450595, 0.075618, 0.141419};
2380     TH1D *r3_muonCorr=(TH1D*)r3_num_Q3->Clone();
2381     for(int i=0; i<10; i++){
2382       r3_muonCorr->SetBinContent(i+1,r3_muonCorr_M1[i]);
2383       r3_muonCorr->SetBinError(i+1,sqrt(pow(r3_muonCorr_M1_e[i],2) + pow(r3Sys_M1[i],2)));
2384     }
2385     r3_muonCorr->SetMarkerStyle(21);
2386     r3_muonCorr->SetMarkerColor(2); r3_muonCorr->SetLineColor(2);
2387     r3_muonCorr->SetFillStyle(1000); r3_muonCorr->SetFillColor(kRed-10);
2388     r3_muonCorr->Draw("E2 same");
2389     r3_num_Q3->Draw("same");
2390     legend3->AddEntry(r3_num_Q3,"No Muon Correction","p");
2391     legend3->AddEntry(r3_muonCorr,"Correction Applied","p");
2392     legend3->Draw("same");*/
2393
2394     // muon correction for muon variation data
2395     /*double r3_muonCorr_M1[10]={0, 1.71115, 1.8128, 1.79761, 1.76112, 1.70699, 1.63908, 1.63417, 1.49861, 1.75565};
2396     double r3_muonCorr_M1_e[10]={0, 0.155228, 0.0491558, 0.0308125, 0.0283357, 0.0378035, 0.0501422, 0.0756057, 0.119169, 0.202093};
2397     TH1D *r3_muonCorr=(TH1D*)r3_num_Q3->Clone();
2398     for(int i=0; i<10; i++){
2399       r3_muonCorr->SetBinContent(i+1,r3_muonCorr_M1[i]);
2400       r3_muonCorr->SetBinError(i+1,sqrt(pow(r3_muonCorr_M1_e[i],2) + pow(r3Sys_M1[i],2)));
2401     }
2402     r3_muonCorr->SetMarkerStyle(21);
2403     r3_muonCorr->SetMarkerColor(2); r3_muonCorr->SetLineColor(2);
2404     r3_muonCorr->SetFillStyle(1000); r3_muonCorr->SetFillColor(kRed-10);
2405     r3_muonCorr->Draw("E2 same");
2406     r3_num_Q3->Draw("same");
2407     legend3->AddEntry(r3_num_Q3,"Muon Corrected.  Default PID","p");
2408     legend3->AddEntry(r3_muonCorr,"Muon Corrected.  Varied PID","p");
2409     legend3->Draw("same");*/
2410
2411     /*for(int i=1; i<=10; i++){
2412       cout<<r3_num_Q3->GetBinContent(i)<<", ";
2413     }
2414     cout<<endl;
2415     for(int i=1; i<=10; i++){
2416       cout<<r3_num_Q3->GetBinError(i)<<", ";
2417       }*/
2418
2419     /*
2420     TH1D *Merged_SanityCheck=(TH1D*)r3_num_Q3->Clone();
2421     for(int i=1; i<=10; i++){
2422       double value = (Gen_r3_m_M0[i-1]+Gen_r3_p_M0[i-1])/2.;
2423       double value_e = sqrt(pow(Gen_r3_m_M0_e[i-1],2)+pow(Gen_r3_p_M0_e[i-1],2))/2.;
2424       //double value = (Gen_r3_m_M0[i-1]+Gen_r3_p_M0[i-1]+Gen_r3_m_M1[i-1]+Gen_r3_p_M1[i-1])/4.;
2425       //double value_e = sqrt(pow(Gen_r3_m_M0_e[i-1],2)+pow(Gen_r3_p_M0_e[i-1],2)+pow(Gen_r3_m_M1_e[i-1],2)+pow(Gen_r3_p_M1_e[i-1],2))/4.;
2426       Merged_SanityCheck->SetBinContent(i,value);
2427       Merged_SanityCheck->SetBinError(i,value_e);
2428     }
2429     gPad->SetTopMargin(0.02); gPad->SetLeftMargin(0.1);
2430     gPad->SetGridx(0); gPad->SetGridy(0);
2431     Merged_SanityCheck->GetXaxis()->SetTitleOffset(1.2);
2432     Merged_SanityCheck->GetYaxis()->SetTitleOffset(1.3);
2433     Merged_SanityCheck->SetMinimum(0.3); Merged_SanityCheck->SetMaximum(2.68);
2434     Merged_SanityCheck->Draw();
2435
2436
2437     //r3_num_Q3->Fit(QuarticFit,"IME","",0,0.1);
2438     Merged_SanityCheck->Fit(QuarticFit,"IME","",0,0.1);
2439     gPad->Update();
2440     //TPaveStats *Quartic_stats=(TPaveStats*)r3_num_Q3->FindObject("stats");
2441     TPaveStats *Quartic_stats=(TPaveStats*)Merged_SanityCheck->FindObject("stats");
2442     Quartic_stats->SetX1NDC(0.15);
2443     Quartic_stats->SetX2NDC(0.52);
2444     Quartic_stats->SetY1NDC(.2);
2445     Quartic_stats->SetY2NDC(.3);
2446
2447     //TH1D *r3_clone=(TH1D*)r3_num_Q3->Clone();
2448     TH1D *r3_clone=(TH1D*)Merged_SanityCheck->Clone();
2449     r3_clone->Fit(QuadraticFit,"IME","",0,0.1);
2450     gPad->Update();
2451     TPaveStats *Quadratic_stats=(TPaveStats*)r3_clone->FindObject("stats");
2452     Quadratic_stats->SetX1NDC(0.54);
2453     Quadratic_stats->SetX2NDC(0.91);
2454     Quadratic_stats->SetY1NDC(.2);
2455     Quadratic_stats->SetY2NDC(.3);
2456
2457     QuarticFit->Draw("same");
2458     QuadraticFit->Draw("same");
2459     Quartic_stats->Draw("same");
2460     Quadratic_stats->Draw("same");
2461
2462     TF1 *ChaoticLimit = new TF1("ChaoticLimit","2.0",0,1);
2463     ChaoticLimit->SetLineStyle(3);
2464     ChaoticLimit->Draw("same");
2465
2466     TLatex *Specif_SanityCheck=new TLatex(0.2,0.35,"HIJING With Simulated QS+FSI Weights");
2467     Specif_SanityCheck->SetNDC();
2468     Specif_SanityCheck->SetTextFont(42);
2469     Specif_SanityCheck->SetTextSize(0.04);
2470     Specif_SanityCheck->Draw();
2471     //TLatex *Specif_Kt3=new TLatex(0.2,0.45,"0.16<K_{t,3}<0.3 GeV/c");
2472     TLatex *Specif_Kt3=new TLatex(0.2,0.45,"0.3<K_{t,3}<1.0 GeV/c");
2473     Specif_Kt3->SetNDC();
2474     Specif_Kt3->SetTextFont(42);
2475     Specif_Kt3->SetTextSize(0.04);
2476     Specif_Kt3->Draw();
2477
2478     legend3->AddEntry(QuarticFit,"Quartic fit","l");
2479     legend3->AddEntry(QuadraticFit,"Quadratic fit","l");
2480     legend3->Draw("same");
2481     //DrawALICELogo(kFALSE, .72, .4, .87, .55);
2482     */
2483     /*
2484     //////////
2485     // r3 Sys
2486     cout<<"r3 values:"<<endl;
2487     for(int ii=0; ii<10; ii++){
2488       cout<<r3_num_Q3->GetBinContent(ii+1)<<", ";
2489     }
2490     cout<<endl;
2491     cout<<"r3 errors:"<<endl;
2492     for(int ii=0; ii<10; ii++){
2493       cout<<r3_num_Q3->GetBinError(ii+1)<<", ";
2494     }
2495     cout<<endl;
2496
2497     // M0,pi- base values with TTC=0.035 (old standard)
2498     //double r3base[10]={0, 2.02584, 1.82659, 1.85384, 1.84294, 1.82431, 1.72205, 1.5982, 1.17379, 0.881535};
2499     //double r3base_e[10]={0, 0.15115, 0.038311, 0.0231936, 0.0209217, 0.0291552, 0.0406093, 0.0624645, 0.0933104, 0.122683};
2500     // M0,pi+ base values with TTC=0.035 (old standard)
2501     //double r3base[10]={0, 1.84738, 1.81805, 1.80057, 1.82563, 1.76585, 1.6749, 1.37454, 1.37558, 0.74562};
2502     //double r3base_e[10]={0, 0.147285, 0.0381688, 0.0231643, 0.0209571, 0.0292255, 0.0407515, 0.0627228, 0.0937758, 0.123392};
2503     // M1,pi- base values with TTC=0.035 (old standard)
2504     //double r3base[10]={0, 1.72641, 1.82547, 1.79747, 1.84141, 1.83557, 1.7994, 1.67511, 1.44238, 1.03004};
2505     //double r3base_e[10]={0, 0.192385, 0.0471262, 0.0270385, 0.0229849, 0.0302041, 0.0397806, 0.0595701, 0.0902909, 0.123976};
2506     // M2,pi- base values with TTC=0.035 (old standard)
2507     //double r3base[10]={0, 1.61726, 1.58973, 1.7834, 1.8914, 1.75445, 1.76511, 1.81367, 1.55617, 1.17797};
2508     //double r3base_e[10]={0, 0.417924, 0.0973277, 0.0527296, 0.0421998, 0.0522792, 0.065144, 0.0923891, 0.135867, 0.190967};
2509     // M9,pi- base values with TTC=0.035 (old standard)
2510     //double r3base[10]={0, 2.51235, 2.68042, 1.85333, 1.66794, 1.39832, 1.717, 1.5337, 1.80686, 1.42286};
2511     //double r3base_e[10]={0, 3.09403, 0.650924, 0.272827, 0.160562, 0.145494, 0.137767, 0.146212, 0.164592, 0.189126};
2512
2513     // M0 PosB, pi- base values with TTC=0.035 (old standard)
2514     //double r3base[10]={0, 1.81361, 1.84779, 1.81693, 1.8252, 1.78366, 1.71822, 1.41243, 1.31656, 0.719842};
2515     //double r3base_e[10]={0, 0.158656, 0.0405, 0.0244075, 0.0219925, 0.0306193, 0.0426494, 0.0655837, 0.0980259, 0.129011};
2516     // M1 PosB, pi- base values with TTC=0.035 (old standard)
2517     //double r3base[10]={0, 1.7208, 1.79016, 1.76858, 1.85008, 1.80416, 1.67274, 1.6501, 1.45098, 1.1052};
2518     //double r3base_e[10]={0, 0.202391, 0.0496228, 0.0282103, 0.0239023, 0.0313499, 0.041234, 0.061679, 0.0934039, 0.128249};
2519
2520     // M1 PID1p5 special case (M0 MomRes used),pi- base values with TTC=0.035 (old standard)
2521     //double r3base[10]={0, 1.77035, 1.82471, 1.77934, 1.80615, 1.77819, 1.71105, 1.60505, 1.44563, 1.09473};
2522     //double r3base_e[10]={0, 0.193326, 0.0471058, 0.0269928, 0.0229581, 0.0301777, 0.0397553, 0.0595472, 0.0902768, 0.123968};
2523
2524     // M0,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2525     //double r3base[10]={0, 1.93236, 1.77265, 1.71324, 1.63426, 1.56167, 1.38348, 1.25064, 0.897615, 0.691928};
2526     //double r3base_e[10]={0, 0.221991, 0.0433103, 0.023234, 0.0187678, 0.0243522, 0.0318668, 0.0459759, 0.066614, 0.0873847};
2527     // M1,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2528     //double r3base[10]={0, 1.46815, 1.76361, 1.68282, 1.6334, 1.56039, 1.44628, 1.27343, 1.01443, 0.691937};
2529     //double r3base_e[10]={0, 0.277164, 0.053134, 0.0271048, 0.0207043, 0.0253691, 0.0315642, 0.0443283, 0.0641466, 0.0858159};
2530     // M2,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2531     //double r3base[10]={0, 1.43331, 1.72759, 1.71155, 1.68957, 1.53846, 1.45209, 1.40698, 1.16202, 1.04672};
2532     //double r3base_e[10]={0, 0.594776, 0.109832, 0.0531013, 0.0383514, 0.0444564, 0.0525963, 0.0705969, 0.0998679, 0.135742};
2533     // M9,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2534     //double r3base[10]={0, 6.98181, 2.42845, 1.71962, 1.57641, 1.23102, 1.49543, 1.48063, 1.50938, 1.27714};
2535     //double r3base_e[10]={0, 7.09254, 0.731841, 0.277336, 0.152199, 0.130775, 0.118981, 0.121864, 0.132875, 0.148409};
2536     //
2537     // M0,pi- base values with 3sigma
2538     //double r3base[10]={0, 1.71114, 1.72142, 1.74749, 1.69154, 1.60411, 1.47805, 1.27561, 1.02355, 0.769853};
2539     //double r3base_e[10]={0, 0.251731, 0.0398122, 0.0196354, 0.0153904, 0.0197609, 0.0258499, 0.0377931, 0.0560878, 0.0774703};
2540     // M1,pi- base values with 3sigma
2541     //double r3base[10]={0, 2.02471, 1.78737, 1.68681, 1.69639, 1.62634, 1.57337, 1.42229, 1.28352, 0.958121};
2542     //double r3base_e[10]={0, 0.314475, 0.048507, 0.022855, 0.0169939, 0.0207322, 0.0260424, 0.0375715, 0.0569927, 0.0830874};
2543     // M2,pi- base values with 3sigma
2544     //double r3base[10]={0, 1.50807, 1.75138, 1.71876, 1.76974, 1.6322, 1.5883, 1.5079, 1.4551, 1.49369};
2545     //double r3base_e[10]={0, 0.658863, 0.0979719, 0.0436736, 0.0308217, 0.035654, 0.0425953, 0.0585648, 0.0865817, 0.128056};
2546     // M9,pi- base values with 3sigma
2547     //double r3base[10]={0, 1.17509, 1.5117, 1.98348, 1.63032, 1.57195, 1.48828, 1.48218, 1.53689, 1.39533};
2548     //double r3base_e[10]={0, 4.89181, 0.60212, 0.230712, 0.11983, 0.102686, 0.0937373, 0.0971281, 0.107949, 0.123923};
2549     //
2550     // M0,pi-,EDbin=0, Femto Plus+Minus
2551     //double r3base[10]={0, 1.77055, 1.75926, 1.7711, 1.67159, 1.57851, 1.41187, 1.11004, 0.870851, 0.442019};
2552     //double r3base_e[10]={0, 0.0688956, 0.0232182, 0.0156323, 0.0155991, 0.0225749, 0.0324584, 0.0516114, 0.0797203, 0.112197};
2553     // M1,pi-,EDbin=0, Femto Plus+Minus
2554     //double r3base[10]={0, 1.83169, 1.79042, 1.7275, 1.67116, 1.6032, 1.5423, 1.30426, 1.13238, 0.656902};
2555     //double r3base_e[10]={0, 0.0891247, 0.0286553, 0.0182713, 0.0171956, 0.0235918, 0.0324304, 0.0510863, 0.0823637, 0.125496};
2556     // M2,pi-,EDbin=0, Femto Plus+Minus
2557     //double r3base[10]={0, 1.60769, 1.81565, 1.75909, 1.74026, 1.62865, 1.58101, 1.50583, 1.3503, 1.27367};
2558     //double r3base_e[10]={0, 0.188422, 0.0581482, 0.0348698, 0.0309363, 0.0401393, 0.0519907, 0.0771148, 0.120583, 0.187679};
2559     // M3,pi-,EDbin=0, Femto Plus+Minus
2560     double r3base[10]={0, 1.52778, 1.81442, 1.78982, 1.69978, 1.70934, 1.59516, 1.56605, 1.45676, 1.17468};
2561     double r3base_e[10]={0, 0.288329, 0.0882351, 0.0510576, 0.0432015, 0.0536294, 0.0667149, 0.0950229, 0.143021, 0.217322};
2562     // M8,pi-,EDbin=0, Femto Plus+Minus
2563     //double r3base[10]={0, 2.27921, 2.40794, 1.72409, 1.66853, 1.70966, 1.52681, 1.64425, 1.18747, 1.35988};
2564     //double r3base_e[10]={0, 1.08315, 0.289559, 0.136951, 0.0919134, 0.0900869, 0.0901754, 0.101887, 0.121289, 0.148063};
2565     // M9,pi-,EDbin=0, Femto Plus+Minus
2566     //double r3base[10]={0, 2.1271, 2.12837, 1.95071, 1.4832, 1.50308, 1.41068, 1.29826, 1.10739, 0.755899};
2567     //double r3base_e[10]={0, 1.52325, 0.385705, 0.181328, 0.116165, 0.10896, 0.105705, 0.115163, 0.131994, 0.156299};
2568
2569     //
2570     // M9,pi+ (No MRC)
2571     double r3base_noMRC[10]={0, 4.16425, 0.429681, 1.56506, 1.64596, 1.73785, 1.57181, 1.51971, 1.59096, 1.54403};
2572     double r3base_noMRC_e[10]={0, 2.00855, 0.437928, 0.187872, 0.1087, 0.097367, 0.0893871, 0.0933358, 0.103603, 0.11839};
2573     // M9,pi+ (MRC normal)
2574     double r3base_MRC[10]={0, 3.84986, 0.558354, 1.60733, 1.67855, 1.75362, 1.59637, 1.54972, 1.62045, 1.57307};
2575     double r3base_MRC_e[10]={0, 1.93428, 0.402333, 0.171664, 0.0996222, 0.0903509, 0.0840615, 0.0889731, 0.100023, 0.115704};
2576     // M9,pi+ (MRC 50% larger)
2577     double r3base_MRC_over[10]={0, 3.70355, 0.614999, 1.62641, 1.69371, 1.76209, 1.60931, 1.56565, 1.63663, 1.58935};
2578     double r3base_MRC_over_e[10]={0, 1.90276, 0.387084, 0.164703, 0.0956857, 0.0872585, 0.0816823, 0.0870033, 0.0984005, 0.114503};
2579
2580
2581
2582     TH1D *r3_noMRC=(TH1D*)r3_num_Q3->Clone();
2583     TH1D *r3_MRC=(TH1D*)r3_num_Q3->Clone();
2584     TH1D *r3_MRC_over=(TH1D*)r3_num_Q3->Clone();
2585
2586     TH1D *r3Sys=(TH1D*)r3_num_Q3->Clone();
2587     for(int ii=1; ii<10; ii++){
2588       double Stat_RelativeError = sqrt(fabs(pow(r3_num_Q3->GetBinError(ii+1),2) - pow(r3base_e[ii],2)));// + for independent stat's
2589       r3Sys->SetBinContent(ii+1, (r3base[ii]-r3_num_Q3->GetBinContent(ii+1))/r3base[ii]);
2590       r3Sys->SetBinError(ii+1, Stat_RelativeError);
2591       //
2592       r3_noMRC->SetBinContent(ii+1, r3base_noMRC[ii]);
2593       r3_MRC->SetBinContent(ii+1, r3base_MRC[ii]);
2594       r3_MRC_over->SetBinContent(ii+1, r3base_MRC_over[ii]);
2595       r3_noMRC->SetBinError(ii+1, r3base_noMRC_e[ii]);
2596       r3_MRC->SetBinError(ii+1, r3base_MRC_e[ii]);
2597       r3_MRC_over->SetBinError(ii+1, r3base_MRC_over_e[ii]);
2598     }
2599     r3_noMRC->SetMarkerStyle(20); r3_noMRC->SetMarkerColor(1); r3_noMRC->SetLineColor(1);
2600     r3_MRC->SetMarkerStyle(23); r3_MRC->SetMarkerColor(4); r3_MRC->SetLineColor(4);
2601     r3_MRC_over->SetMarkerStyle(24); r3_MRC_over->SetMarkerColor(2); r3_MRC_over->SetLineColor(2);
2602
2603     //r3_MRC->Draw();
2604     //r3_noMRC->Draw("same");
2605     //r3_MRC_over->Draw("same");
2606     TLegend *legendMRC = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2607     legendMRC->SetBorderSize(1);
2608     legendMRC->SetTextSize(.03);// small .03; large .06
2609     legendMRC->SetFillColor(0);
2610     //legendMRC->AddEntry(r3_MRC,"Momentum Resolution Corrected","p");
2611     //legendMRC->AddEntry(r3_noMRC,"No Correction","p");
2612     //legendMRC->AddEntry(r3_MRC_over,"50% Increased Correction","p");
2613     //legendMRC->Draw("same");
2614
2615     r3Sys->GetYaxis()->SetTitle("(r_{3}^{t1}-r_{3}^{t2})/r_{3}^{t1}");
2616     r3Sys->GetYaxis()->SetTitleOffset(1.8);
2617     r3Sys->SetMinimum(-0.1);
2618     r3Sys->SetMaximum(0.1);
2619     r3Sys->GetXaxis()->SetRangeUser(0,0.099);
2620     r3Sys->Draw();
2621     TF1 *lineFit=new TF1("lineFit","pol0",0,0.02);
2622     gStyle->SetOptFit(111);
2623     r3Sys->Fit(lineFit,"MEQ","",0,0.1);
2624     */
2625
2626
2627     /*
2628     ////////////////////////////////////////////////////////////////////////
2629     // The STAR method
2630     ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation
2631
2632     TLegend *legendSTAR = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2633     legendSTAR->SetBorderSize(1);
2634     legendSTAR->SetTextSize(.03);// small .03; large .06
2635     legendSTAR->SetFillColor(0);
2636
2637     TH1D *r3_STAR = new TH1D("r3_STAR","",Q3BINS,0,0.2);
2638     TH1D *r3_den_STAR = new TH1D("r3_den_STAR","",Q3BINS,0,0.2);
2639     double r3_STAR_num_e[20]={0};
2640     double r3_STAR_den_e[20]={0};
2641
2642     for(int i=2; i<=20; i++){// bin number
2643       double Q12 = (i-0.5)*(0.005);
2644       int Q12bin = int(Q12/0.005)+1;
2645       //
2646       for(int j=2; j<=20; j++){// bin number
2647         double Q13 = (j-0.5)*(0.005);
2648         int Q13bin = int(Q13/0.005)+1;
2649         //
2650         for(int k=2; k<=20; k++){// bin number
2651           double Q23 = (k-0.5)*(0.005);
2652           int Q23bin = int(Q23/0.005)+1;
2653           //
2654           double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
2655           int Q3bin = Cterm1_noMRC->GetXaxis()->FindBin(Q3);
2656           if(Q3bin>19) continue;
2657
2658           if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;
2659           if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;
2660
2661           double K2_12 = CoulCorr2(CP12, Q12);
2662           double K2_13 = CoulCorr2(CP13, Q13);
2663           double K2_23 = CoulCorr2(CP23, Q23);
2664
2665           double r3num = Cterm1_noMRC->GetBinContent(Q3bin)/(K2_12*K2_13*K2_23) - 1;
2666           r3num -= C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1;
2667           r3num -= C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1;
2668           r3num -= C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1;
2669           double r3den = (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1);
2670           r3den *= (C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1);
2671           r3den *= (C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1);
2672           if(r3den<0) continue;
2673           r3den = sqrt(r3den);
2674           //
2675           r3_STAR_num_e[Q3bin-1] += pow(Cterm1_noMRC->GetBinError(Q3bin)/(K2_12*K2_13*K2_23),2);
2676           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q12bin)/(K2_12),2);
2677           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q13bin)/(K2_13),2);
2678           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q23bin)/(K2_23),2);
2679           r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q12bin)/(K2_12) * (C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1)*(C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1) /r3den,2);
2680           r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q13bin)/(K2_13) * (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1)*(C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1) /r3den,2);
2681           r3_STAR_den_e[Q3bin-1] += pow(0.5*C2ssRaw->GetBinError(Q23bin)/(K2_23) * (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1)*(C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1) /r3den,2);
2682           //
2683           r3_STAR->Fill(Q3, r3num);
2684           r3_den_STAR->Fill(Q3, r3den);
2685         }
2686       }
2687     }
2688
2689
2690     for(int i=0; i<20; i++){
2691       r3_STAR->SetBinError(i+1, sqrt( r3_STAR_num_e[i]));
2692       r3_den_STAR->SetBinError(i+1, sqrt( r3_STAR_den_e[i]));
2693     }
2694     r3_STAR->Divide(r3_den_STAR);
2695     r3_STAR->SetMaximum(5);
2696     r3_STAR->SetMinimum(-10);
2697     r3_STAR->GetXaxis()->SetRangeUser(0,0.099);
2698     r3_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
2699     r3_STAR->GetYaxis()->SetTitle("r_{3}*");
2700     r3_STAR->SetLineColor(2); r3_STAR->SetMarkerColor(2);
2701     r3_STAR->SetMarkerStyle(kFullStar);
2702     r3_STAR->SetMarkerSize(1.2);
2703     r3_STAR->Draw();
2704     r3_STAR->Fit(QuarticFit,"IME","",0.01,0.1);
2705     gPad->Update();
2706     TPaveStats *Quartic_stats_STAR=(TPaveStats*)r3_STAR->FindObject("stats");
2707     Quartic_stats_STAR->SetX1NDC(0.2);
2708     Quartic_stats_STAR->SetX2NDC(0.5);
2709     Quartic_stats_STAR->SetY1NDC(.8);
2710     Quartic_stats_STAR->SetY2NDC(.9);
2711
2712     TH1D *r3_clone_STAR=(TH1D*)r3_STAR->Clone();
2713     r3_STAR->Fit(QuadraticFit,"IME","",0,0.1);
2714     gPad->Update();
2715     TPaveStats *Quadratic_stats_STAR=(TPaveStats*)r3_clone_STAR->FindObject("stats");
2716     Quadratic_stats_STAR->SetX1NDC(0.55);
2717     Quadratic_stats_STAR->SetX2NDC(0.85);
2718     Quadratic_stats_STAR->SetY1NDC(.8);
2719     Quadratic_stats_STAR->SetY2NDC(.9);
2720
2721     QuarticFit->Draw("same");
2722     QuadraticFit->Draw("same");
2723     Quartic_stats_STAR->Draw("same");
2724     Quadratic_stats_STAR->Draw("same");
2725     //Cterm1_noMRC->Draw();
2726
2727     /////////////////////////////////////////////////////////////////////////
2728     if(SameCharge) r3_num_Q3->Draw("same");
2729     legendSTAR->AddEntry(r3_STAR,"STAR method","p");
2730     legendSTAR->AddEntry(r3_num_Q3,"ALICE method","p");
2731     legendSTAR->Draw("same");
2732     */
2733   }// SameCharge
2734
2735   //c3_hist_STAR->Draw();
2736   //Cterm1_noMRC->Draw();
2737
2738   //r3_num_Q3->SetMarkerStyle(20);
2739   //r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
2740   //r3_num_Q3->SetMinimum(-7);
2741   //r3_num_Q3->SetMaximum(2.2);
2742
2743
2744   /*
2745     double C3_star_noCoul[20]={0, 1.59478, 1.49099, 1.4239, 1.31166, 1.21337, 1.14344, 1.09671, 1.06525, 1.04391, 1.0306, 1.02045, 1.01359, 1.0035, 0.992264, 0.989587, 0.989968, 0.996795, 0, 0};
2746   double C3_star_noCoul_e[20]={0, 0.260834, 0.0191839, 0.00524928, 0.0021307, 0.00110448, 0.000661211, 0.000438894, 0.000311247, 0.00023452, 0.000182819, 0.000147496, 0.000126547, 0.000130428, 0.000159494, 0.000240658, 0.000489113, 0.00387954, 0, 0};
2747   TH1D *C3_star=(TH1D*)Cterm1->Clone();
2748   C3_star->Add(C3_star,-1);
2749   for(int i=0; i<20; i++) {C3_star->SetBinContent(i+1, C3_star_noCoul[i]); C3_star->SetBinError(i+1, C3_star_noCoul_e[i]);}
2750   C3_star->SetMarkerColor(4);
2751   C3_star->SetLineColor(4);
2752   //C3_star->Draw("same");
2753   //legend1->AddEntry(C3_star,"---; 200 GeV","p");
2754
2755   //legend1->Draw("same");
2756
2757   //TGraph *gr_MomRes_term1=new TGraph(10,
2758   */
2759
2760
2761
2762
2763
2764   //////////////////////////////////////////////////////////////////////////////////////
2765
2766
2767
2768
2769
2770
2771   if(SaveToFile){
2772     TString *savefilename = new TString("OutFiles/OutFile");
2773     if(SameCharge) savefilename->Append("SC");
2774     else savefilename->Append("MC");
2775     //
2776     if(CHARGE==1) savefilename->Append("Pos");
2777     else savefilename->Append("Neg");
2778     //
2779     if(IncludeG) savefilename->Append("YesG");
2780     else savefilename->Append("NoG");
2781     //
2782     if(IncludeEWfromTherm) savefilename->Append("EWfromTherm");
2783     if(IncludeEW) savefilename->Append("EW");
2784     if(GRS) savefilename->Append("GRS");
2785     else savefilename->Append("Omega0");
2786
2787     if(EDbin==0) savefilename->Append("Kt3_1");
2788     else savefilename->Append("Kt3_2");
2789
2790     savefilename->Append("_Kt");
2791     *savefilename += Ktbin;
2792     savefilename->Append("_M");
2793     *savefilename += Mbin;
2794     savefilename->Append(".root");
2795     TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
2796     can1->Write("can1");
2797     C2_ss->Write("C2_ss");
2798     C2_os->Write("C2_os");
2799     C2_ss_Momsys->Write("C2_ss_Momsys");
2800     C2_os_Momsys->Write("C2_os_Momsys");
2801     C2_ss_Ksys->Write("C2_ss_Ksys");
2802     C2_os_Ksys->Write("C2_os_Ksys");
2803     fitC2ss->Write("fitC2ss");
2804     fitC2os->Write("fitC2os");
2805     fitC2ss_h->Write("fitC2ss_h");
2806     fitC2os_h->Write("fitC2os_h");
2807     if(IncludeEWfromTherm) {
2808       fitC2ssTherm->Write("fitC2ssTherm");
2809       fitC2osTherm->Write("fitC2osTherm");
2810       fitC2ssThermGaus->Write("fitC2ssThermGaus");
2811       fitC2osThermGaus->Write("fitC2osThermGaus");
2812     }
2813     C2Therm_ss->Write("C2Therm_ss");
2814     C2Therm_os->Write("C2Therm_os");
2815     K2_ss->Write("K2_ss");
2816     K2_os->Write("K2_os");
2817     Cterm1->Write("C3");
2818     Coul_Omega0->Write("Coul_Omega0");
2819     Coul_GRiverside->Write("Coul_GRiverside");
2820     GenSignalExpected_num->Write("C3_EWexpectation");
2821     num_QS->Write("num_QS");
2822     c3_hist->Write("c3");
2823     Combinatorics_1d->Write("Combinatorics_1d");
2824     ChiSquaredNDF->Write("ChiSquaredNDF");
2825     r3_num_Q3->Write("r3_Q3");
2826
2827     //
2828     savefile->Close();
2829   }
2830
2831
2832 }
2833
2834 void ReadCoulCorrections(int RVal, int bVal, int kt){
2835   ///////////////////////
2836
2837   TString *fname2;
2838   if(FileSetting==6) fname2 = new TString("KFile_GaussR8to5.root");
2839   else if(FileSetting==7) fname2 = new TString("KFile_GaussR11to6.root");
2840   else{
2841     if(SourceType==0) fname2 = new TString("KFile.root");
2842     if(SourceType==1) fname2 = new TString("KFile_50fmCut.root");
2843     if(SourceType==2) fname2 = new TString("KFile_GaussR11to6.root");
2844   }
2845
2846   TFile *File=new TFile(fname2->Data(),"READ");
2847   if(bVal!=2 && bVal!=3 && bVal!=5 && bVal!=7 && bVal!=8 && bVal!=9) cout<<"Therminator bVal not acceptable in 2-particle Coulomb read"<<endl;
2848
2849   if(kt==10){// kt integrated
2850     TH2D *tempT_ss = (TH2D*)File->Get("K2ssT");
2851     TH2D *tempT_os = (TH2D*)File->Get("K2osT");
2852     CoulCorr2SS = (TH1D*)tempT_ss->ProjectionY("CoulCorr2SS",bBin, bBin);
2853     CoulCorr2OS = (TH1D*)tempT_os->ProjectionY("CoulCorr2OS",bBin, bBin);
2854   }else{
2855     if(kt < 1 || kt > 6) cout<<"kt bin out of range in 2-particle Coulomb read"<<endl;
2856     TH3D *tempT3_ss = (TH3D*)File->Get("K2ssT_kt");
2857     TH3D *tempT3_os = (TH3D*)File->Get("K2osT_kt");
2858     CoulCorr2SS = (TH1D*)tempT3_ss->ProjectionZ("CoulCorr2SS",bBin, bBin, kt,kt);
2859     CoulCorr2OS = (TH1D*)tempT3_os->ProjectionZ("CoulCorr2OS",bBin, bBin, kt,kt);
2860   }
2861
2862   if(IncludeStrongSC){// include strong FSI for pi+pi+ as per Lednicky code with R=7 (ratio of Coul+Strong to Coul)
2863     for(int ii=1; ii<=CoulCorr2SS->GetNbinsX(); ii++){
2864       double newValue = CoulCorr2SS->GetBinContent(ii) * StrongSC->Eval(CoulCorr2SS->GetBinCenter(ii)*1000/2.);// k* not qinv
2865       CoulCorr2SS->SetBinContent(ii, newValue);
2866     }
2867   }
2868   // K factor shift for testing
2869   for(int ii=1; ii<=CoulCorr2SS->GetNbinsX(); ii++){
2870     double newValueSS = (CoulCorr2SS->GetBinContent(ii)-1)*KShift+1;// k* not qinv
2871     CoulCorr2SS->SetBinContent(ii, newValueSS);
2872     double newValueOS = (CoulCorr2OS->GetBinContent(ii)-1)*KShift+1;// k* not qinv
2873     CoulCorr2OS->SetBinContent(ii, newValueOS);
2874   }
2875
2876   CoulCorr2SS->SetDirectory(0);
2877   CoulCorr2OS->SetDirectory(0);
2878   File->Close();
2879 }
2880 double CoulCorr2(int chargeproduct, double Q2){// returns K2
2881   int indexL=0;
2882   int indexH=0;
2883   double slope=0;
2884
2885
2886   indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2887   indexH = indexL+1;
2888
2889   if(indexH >= CoulCorr2SS->GetNbinsX()) return 1.0;
2890   if(chargeproduct==1){
2891     slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
2892     slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
2893     return slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1);
2894   }else {
2895     slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1));
2896     slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1));
2897     return slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
2898   }
2899
2900 }
2901
2902 //----------------------------------------------------------------------
2903
2904
2905 //________________________________________________________________________
2906 void fcnC2_Global(int& npar, double* deriv, double& f, double par[], int flag){
2907
2908   double qinvSS=0;
2909
2910   double Rch=par[3]/FmToGeV;
2911   double Rcoh=par[4]/FmToGeV;
2912   double Dp=0;
2913   double CSS=0, COS=0;
2914   double SumChi2=0;
2915   double EW=0;
2916   double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2917   NFitPoints_C2global=0;
2918   if(LinkN) par[9]=par[0];// Link N factors
2919
2920   for(int i=1; i<BINRANGE_C2global; i++){
2921
2922     qinvSS = AvgQinvSS[i];
2923     //if(qinvSS>0.08) continue;
2924
2925     if(!GofP) Dp=fabs(par[2])/(1-fabs(par[2]));// p independent
2926     else Dp = fabs(par[2])/(1-fabs(par[2]));
2927
2928     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2) - pow(RT,2))*pow(AvgP[Ktbin_GofP-1][i]-qinvSS/2.,2));
2929     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2) - pow(RT,2))*pow(AvgP[Ktbin_GofP-1][i]+qinvSS/2.,2));
2930     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][i]-qinvSS/2.,2) - 2*MT/Temp);
2931     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][i]+qinvSS/2.,2) - 2*MT/Temp);
2932     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]-qinvSS/2.,2) - 2*MT/Temp);
2933     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]+qinvSS/2.,2) - 2*MT/Temp);
2934     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2935     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2936     double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2937     double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2938
2939     if(!GofP) {Dp1=Dp; Dp2=Dp;}
2940
2941     //
2942     double arg=qinvSS*Rch;
2943     EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2944     EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2945     EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2946     EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2947     //
2948     double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
2949     double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
2950     CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
2951     if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2952     else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
2953     //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
2954
2955     CSS *= par[1]*K2SS[i];
2956     if(MuonCorrection) CSS /= C2muonCorrection_SC->GetBinContent(i+1);
2957     CSS += 1-par[1];
2958     CSS *= par[0];
2959     //
2960     COS = 1;
2961     if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2962     COS *= par[1]*K2OS[i];
2963     if(MuonCorrection) COS /= C2muonCorrection_MC->GetBinContent(i+1);
2964     COS += 1-par[1];
2965     COS *= par[9];// different Norm
2966     //
2967     if(C2ssFitting[i] > 0){
2968       if(IncludeSS) {
2969         double errorSS = sqrt(pow( (CSS/par[0] - (1-par[1]))/K2SS[i] * K2SS_e[i] * par[0],2) + pow(C2ssFitting_e[i],2));
2970         //double errorSS = C2ssFitting_e[i];
2971         SumChi2 += pow((CSS - C2ssFitting[i])/errorSS,2);
2972         NFitPoints_C2global++;
2973       }
2974     }
2975     if(IncludeOS) {
2976       double errorOS = sqrt(pow( (COS/par[9] - (1-par[1]))/K2OS[i] * K2OS_e[i] * par[9],2) + pow(C2osFitting_e[i],2));
2977       //double errorOS = C2osFitting_e[i];
2978       SumChi2 += pow((COS - C2osFitting[i])/errorOS,2);
2979       NFitPoints_C2global++;
2980     }
2981   }
2982
2983
2984   f = SumChi2;
2985   Chi2_C2global = f;
2986
2987 }
2988 //________________________________________________________________________
2989 double C2ssFitFunction(Double_t *x, Double_t *par)
2990 {
2991   double Rch=par[3]/FmToGeV;
2992   double Rcoh=par[4]/FmToGeV;
2993   double Dp=0;
2994   int qbin = int(fabs(x[0]/0.005));
2995   if(qbin >= BINRANGE_C2global) return 1.0;
2996   double qinvSS = AvgQinvSS[qbin];
2997   double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2998
2999   if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
3000   else Dp = fabs(par[2])/(1-fabs(par[2]));
3001
3002   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]-qinvSS/2.,2));
3003   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]+qinvSS/2.,2));
3004   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][qbin]-qinvSS/2.,2) - 2*MT/Temp);
3005   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(AvgP[Ktbin_GofP-1][qbin]+qinvSS/2.,2) - 2*MT/Temp);
3006   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
3007   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
3008   double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
3009   double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
3010
3011   if(!GofP) {Dp1=Dp; Dp2=Dp;}
3012   double arg=qinvSS*Rch;
3013   double EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
3014   EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3015   EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
3016   EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
3017   double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
3018   double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
3019   double CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
3020   if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
3021   else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
3022   //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
3023   CSS *= par[1]*K2SS[qbin];
3024   if(MuonCorrection) CSS /= C2muonCorrection_SC->GetBinContent(qbin+1);
3025   CSS += 1-par[1];
3026   CSS *= par[0];
3027
3028   //cout<<qinvSS<<"  "<<Dp1/(1+Dp1) * Dp2/(1+Dp2)<<endl;
3029
3030   return CSS;
3031 }
3032 //________________________________________________________________________
3033 double C2osFitFunction(Double_t *x, Double_t *par)
3034 {
3035   if(LinkN) par[9]=par[0];// Link N factors
3036   double Rcoh=par[4]/FmToGeV;
3037   double Dp=0;
3038   int qbin = int(fabs(x[0]/0.005));
3039   if(qbin >= BINRANGE_C2global) return 1.0;
3040   double qinvOS = AvgQinvOS[qbin];
3041
3042   if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
3043   else Dp = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin],2));
3044   //Dp = fabs(par[2])/(1-fabs(par[2]));
3045   double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]-qinvOS/2.,2));
3046   double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin]+qinvOS/2.,2));
3047   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
3048   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
3049
3050   if(!GofP) {Dp1=Dp; Dp2=Dp;}
3051   double COS = 1;
3052   if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
3053   COS *= par[1]*K2OS[qbin];
3054   if(MuonCorrection) COS /= C2muonCorrection_MC->GetBinContent(qbin+1);
3055   COS += 1-par[1];
3056   COS *= par[9];
3057   return COS;
3058 }
3059 //________________________________________________________________________
3060 void fcnC3_Global(int& npar, double* deriv, double& f, double par[], int flag){
3061
3062   double qinvSS_12=0;
3063   double qinvSS_13=0;
3064   double qinvSS_23=0;
3065
3066   //double D=fabs(par[2])/(1-fabs(par[2]));
3067   double Rch=par[3]/FmToGeV;
3068   //double Rcoh=par[4]/FmToGeV;
3069   double CSS=0;
3070   double SumChi2=0;
3071   double EW_12=0, EW_13=0, EW_23=0;
3072   NFitPoints_C3global=0;
3073
3074   for(int i=0; i<BINRANGE_C2global; i++){
3075     qinvSS_12 = AvgQinvSS[i];
3076     double arg_12=qinvSS_12*Rch;
3077     EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
3078     EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
3079     EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
3080     EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
3081
3082     for(int j=0; j<BINRANGE_C2global; j++){
3083       qinvSS_13 = AvgQinvSS[j];
3084       double arg_13=qinvSS_13*Rch;
3085       EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
3086       EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
3087       EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
3088       EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
3089
3090       for(int k=0; k<BINRANGE_C2global; k++){
3091         qinvSS_23 = AvgQinvSS[k];
3092         double arg_23=qinvSS_23*Rch;
3093         EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
3094         EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
3095         EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
3096         EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
3097         //
3098         if(i==0 || j==0 || k==0) continue;
3099         CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
3100         CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
3101         CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
3102         CSS += 2*EW_12*EW_13*EW_23 * exp(-1/2.*pow(Rch,2)*(pow(qinvSS_12,2)+pow(qinvSS_13,2)+pow(qinvSS_23,2)));
3103         //
3104         if(IncludeSS) {
3105           /*double errorSS = sqrt(pow(C3ssFitting_e[i],2));
3106           if(errorSS > 0){
3107             SumChi2 += pow((CSS - C3ssFitting[i])/errorSS,2);
3108             NFitPoints_C3global++;
3109             }*/
3110         }
3111         if(IncludeOS) {
3112           //double errorOS = sqrt(pow(C3osFitting_e[i],2));
3113           //SumChi2 += pow((COS - C3osFitting[i])/errorOS,2);
3114           //NFitPoints_C3global++;
3115         }
3116
3117
3118       }// k
3119     }// j
3120   }// i
3121
3122   f = SumChi2;
3123   Chi2_C3global = f;
3124
3125 }
3126 //______________________________________________________________________
3127 /*void C3ssFitFunction(Double_t *x, Double_t *par){
3128
3129   double qinvSS_12=0, qinvOS_12=0;
3130   double qinvSS_13=0, qinvOS_13=0;
3131   double qinvSS_23=0, qinvOS_23=0;
3132
3133   double D=fabs(par[2])/(1-fabs(par[2]));
3134   double Rch=par[3]/FmToGeV;
3135   double Rcoh=par[4]/FmToGeV;
3136   double CSS=0, COS=0;
3137   double term1=0, term2=0;
3138   double SumChi2=0;
3139   double EW_12=0, EW_13=0, EW_23=0;
3140   NFitPoints_C3global=0;
3141
3142
3143   int qbin_12 = x[0]*1000/5;
3144   qinvSS_12 = AvgQinvSS[qbin_12];
3145   qinvOS_12 = AvgQinvOS[qbin_12];
3146   double arg_12=qinvSS_12*Rch;
3147   EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
3148   EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
3149   EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
3150   EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
3151   //
3152   int qbin_13 = x[1]*1000/5;
3153   qinvSS_13 = AvgQinvSS[qbin_13];
3154   qinvOS_13 = AvgQinvOS[qbin_13];
3155   double arg_13=qinvSS_13*Rch;
3156   EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
3157   EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
3158   EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
3159   EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
3160   //
3161   int qbin_23 = x[2]*1000/5;
3162   qinvSS_23 = AvgQinvSS[qbin_23];
3163   qinvOS_23 = AvgQinvOS[qbin_23];
3164   double arg_23=qinvSS_23*Rch;
3165   EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
3166   EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
3167   EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
3168   EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
3169   //
3170   if(qbin_12==0 || qbin_13==0 || qbin_23==0) return 0;
3171   CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
3172   CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
3173   CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
3174   CSS += 2*EW_12*EW_13*EW_23 * exp(-1/2.*pow(Rch,2)*(pow(qinvSS_12,2)+pow(qinvSS_13,2)+pow(qinvSS_23,2)));
3175   //
3176
3177   return CSS;
3178   }*/
3179 //----------------------------------------------------------------------
3180 void fcnOSL(int& npar, double* deriv, double& f, double par[], int flag){
3181
3182   double qout=0, qside=0, qlong=0;
3183   double Gaus_ch=0;
3184   double C=0;
3185   //double lnL=0;
3186   double SumChi2=0;
3187   NFitPoints_OSL=0;
3188
3189   for(int i=0; i<BINS_OSL-10; i++){// out, 0 to BINS_OSL-10
3190     for(int j=0; j<BINS_OSL-10; j++){// side, 0 to BINS_OSL-10
3191       for(int k=0; k<BINS_OSL-10; k++){// long, 0 to BINS_OSL-10
3192
3193         if(B[i][j][k] < 1. ) continue;
3194
3195         qout = (i+0.5)*0.005;
3196         qside = (j+0.5)*0.005;
3197         qlong = (k+0.5)*0.005;
3198
3199         Gaus_ch = exp(-pow(par[3]*qout/FmToGeV,2) - pow(par[4]*qside/FmToGeV,2) - pow(par[5]*qlong/FmToGeV,2));
3200
3201         //C = par[0]*(1 + par[1]*(K_OSL[i][j][k]-1) + K_OSL[i][j][k]*par[1]*Gaus_ch);
3202         C = par[0]*( (1-par[1]) + par[1]*K_OSL[i][j][k]*(1+Gaus_ch));
3203
3204         if(C < 0) continue;
3205
3206         double error = pow(A_e[i][j][k]/B[i][j][k],2) + pow(B_e[i][j][k]*A[i][j][k]/pow(B[i][j][k],2),2);
3207         error += pow((C/par[0] - (1-par[1]))/K_OSL[i][j][k] * K_OSL_e[i][j][k]*par[0],2);
3208         error = sqrt(error);
3209         SumChi2 += pow( (C - A[i][j][k]/B[i][j][k])/error,2);
3210
3211         /*
3212         if(A[i][j][k] >= 1) term1 = C*(A[i][j][k]+B[i][j][k])/(A[i][j][k]*(C+1));
3213         else term1 = 0;
3214         term2 = (A[i][j][k]+B[i][j][k])/(B[i][j][k]*(C+1));
3215
3216         if(term1 > 0.0 && term2 > 0.0){
3217           lnL += A[i][j][k]*log(term1) + B[i][j][k]*log(term2);
3218         }else if(term1==0 && term2 > 0.0){
3219           lnL += B[i][j][k]*log(term2);
3220         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
3221         */
3222         NFitPoints_OSL++;
3223       }
3224     }
3225   }
3226   f = SumChi2;
3227   //f = -2.0*lnL;
3228   Chi2_OSL = f;
3229
3230 }
3231 //________________________________________________________________________
3232 double OSLfitfunction(Double_t *x, Double_t *par)
3233 {
3234   int bin_out = int(1000*(x[0])/5.);
3235   int bin_side = int(1000*(x[1])/5.);
3236   int bin_long = int(1000*(x[2])/5.);
3237
3238   double K = CoulCorr2(+1, avg_q[bin_out][bin_side][bin_long]);
3239   double Gaus_ch = exp(-pow(par[3]*x[0]/FmToGeV,2) - pow(par[4]*x[1]/FmToGeV,2) - pow(par[5]*x[2]/FmToGeV,2));
3240   return par[0]*(1 + par[1]*(K-1) + K*par[1]*Gaus_ch);
3241
3242 }
3243 //__________________________________________________________________________
3244
3245 void fcn_3(int& npar, double* deriv, double& f, double par[], int flag){
3246
3247   double q12=0, q13=0, q23=0;
3248   double K12=0, K13=0, K23=0, K123=0;
3249   double G1=0, G2=0, G3=0, G4=0;
3250   double EW12=0, EW13=0, EW23=0;
3251   double C=0;
3252   double term1=0, term2=0;
3253   double lnL=0;
3254
3255   // start from 2-4 MeV bins
3256   for(int i=1; i<BINRANGE_3-10; i++){// q12
3257     for(int j=1; j<BINRANGE_3-10; j++){// q13
3258       for(int k=1; k<BINRANGE_3-10; k++){// q23
3259
3260         if(B_3[i][j][k] < 1. ) continue;
3261
3262         q12 = (i+0.5)*0.002;
3263         q13 = (j+0.5)*0.002;
3264         q23 = (k+0.5)*0.002;
3265
3266         //K12 = 1/CoulCorr_CsorgoMate(2, q12);
3267         //K13 = 1/CoulCorr_CsorgoMate(2, q13);
3268         //K23 = 1/CoulCorr_CsorgoMate(2, q23);
3269         //K12 = 1/Coul_pipi_data(q12);
3270         //K13 = 1/Coul_pipi_data(q13);
3271         //K23 = 1/Coul_pipi_data(q23);
3272         K12=1.0; K13=1.0; K23=1.0;
3273         K123 = K12*K13*K23;
3274         //
3275         //K123 = 1 - pow(par[1],3)*(1-K12*K13*K23);
3276         //K12 = 1 - pow(par[1],3)*(1-K12);
3277         //K13 = 1 - pow(par[1],3)*(1-K13);
3278         //K23 = 1 - pow(par[1],3)*(1-K23);
3279         //
3280         EW12 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q12*par[2]/FmToGeV,3) - 12*pow(q12*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q12*par[2]/FmToGeV,4) -48*pow(q12*par[2]/FmToGeV,2) + 12);
3281         EW13 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q13*par[2]/FmToGeV,3) - 12*pow(q13*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q13*par[2]/FmToGeV,4) -48*pow(q13*par[2]/FmToGeV,2) + 12);
3282         EW23 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(q23*par[2]/FmToGeV,3) - 12*pow(q23*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(q23*par[2]/FmToGeV,4) -48*pow(q23*par[2]/FmToGeV,2) + 12);
3283         //EW12=1; EW13=1; EW23=1;
3284         //
3285         G1 = K12*exp(-pow(par[2]*q12/FmToGeV,2)/2.)*EW12;
3286         G1 += K13*exp(-pow(par[2]*q13/FmToGeV,2)/2.)*EW13;
3287         G1 += K23*exp(-pow(par[2]*q23/FmToGeV,2)/2.)*EW23;
3288         //
3289         G2 = K12*exp(-pow(par[2]*q12/FmToGeV,2))*pow(EW12,2);
3290         G2 += K13*exp(-pow(par[2]*q13/FmToGeV,2))*pow(EW13,2);
3291         G2 += K23*exp(-pow(par[2]*q23/FmToGeV,2))*pow(EW23,2);
3292         //
3293         G3 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q23*q23)/2.)*EW12*EW23;
3294         G3 += exp(-pow(par[2]/FmToGeV,2)*(q13*q13 + q23*q23)/2.)*EW13*EW23;
3295         G3 += exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13)/2.)*EW12*EW13;
3296         G3 *= K123;
3297         //
3298         G4 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13 + q23*q23)/2.);
3299         G4 *= EW12*EW13*EW23;
3300         G4 *= K123;
3301         //
3302         C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3303         //C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3304         C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3305         C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3306         C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3307         C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3308         C *= par[0];// Norm
3309
3310
3311
3312         if(C < 0) continue;
3313
3314         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));
3315         else term1 = 0;
3316         term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1));
3317
3318         if(term1 > 0.0 && term2 > 0.0){
3319           lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2);
3320         }else if(term1==0 && term2 > 0.0){
3321           lnL += B_3[i][j][k]*log(term2);
3322         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
3323
3324       }
3325     }
3326   }
3327
3328   f = -2.0*lnL;
3329
3330 }
3331 //________________________________________________________________________
3332 double D3fitfunction_3(Double_t *x, Double_t *par)
3333 {
3334   double K12=CoulCorr2(+1, x[0]);
3335   double K13=CoulCorr2(+1, x[1]);
3336   double K23=CoulCorr2(+1, x[2]);
3337   K12=1.0; K13=1.0; K23=1.0;
3338   double K123=K12*K13*K23;
3339   //
3340   //K123 = 1 - pow(par[1],3)*(1 - K123);
3341   //K12 = 1 - pow(par[1],3)*(1-K12);
3342   //K13 = 1 - pow(par[1],3)*(1-K13);
3343   //K23 = 1 - pow(par[1],3)*(1-K23);
3344   //
3345   double EW12 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[0]*par[2]/FmToGeV,3) - 12*pow(x[0]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[0]*par[2]/FmToGeV,4) -48*pow(x[0]*par[2]/FmToGeV,2) + 12);
3346   double EW13 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[1]*par[2]/FmToGeV,3) - 12*pow(x[1]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[1]*par[2]/FmToGeV,4) -48*pow(x[1]*par[2]/FmToGeV,2) + 12);
3347   double EW23 = 1 + par[4]/(6.*pow(2,1.5))*(8*pow(x[2]*par[2]/FmToGeV,3) - 12*pow(x[2]*par[2]/FmToGeV,1)) + par[5]/(24.*pow(2,2))*(16*pow(x[2]*par[2]/FmToGeV,4) -48*pow(x[2]*par[2]/FmToGeV,2) + 12);
3348   //EW12=1; EW13=1; EW23=1;
3349   //
3350   double G1=0, G2=0, G3=0, G4=0;
3351
3352   G1 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2)/2.)*EW12;
3353   G1 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2)/2.)*EW13;
3354   G1 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2)/2.)*EW23;
3355   //
3356   G2 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2))*pow(EW12,2);
3357   G2 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2))*pow(EW13,2);
3358   G2 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2))*pow(EW23,2);
3359   //
3360   G3 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[2]*x[2])/2.)*EW12*EW23;
3361   G3 += exp(-pow(par[2]/FmToGeV,2)*(x[1]*x[1] + x[2]*x[2])/2.)*EW13*EW23;
3362   G3 += exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1])/2.)*EW12*EW13;
3363   G3 *= K123;
3364   //
3365   G4 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/2.);
3366   G4 *= EW12*EW13*EW23;
3367   G4 *= K123;
3368   //
3369   double C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3370   //double C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3371   C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3372   C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3373   C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3374   C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3375   C *= par[0];// Norm
3376
3377   return C;
3378 }
3379
3380 void ReadCoulCorrections_Omega0(){
3381   // read in 3d 3-particle coulomb correlations = K3
3382
3383   TString *fname;
3384   if(FileSetting==6) fname = new TString("KFile_GaussR8to5.root");
3385   else if(FileSetting==7) fname = new TString("KFile_GaussR11to6.root");
3386   else{
3387     if(SourceType==0) fname = new TString("KFile.root");
3388     if(SourceType==1) fname = new TString("KFile_50fmCut.root");
3389     if(SourceType==2) fname = new TString("KFile_GaussR11to6.root");
3390   }
3391   TFile *coulfile=new TFile(fname->Data(),"READ");
3392
3393   TString *name=new TString("K3ss_");
3394   *name += bBin-1;
3395   CoulCorrOmega0SS = (TH3D*)coulfile->Get(name->Data());
3396   CoulCorrOmega0SS->SetDirectory(0);
3397   name=new TString("K3os_");
3398   *name += bBin-1;
3399   CoulCorrOmega0OS = (TH3D*)coulfile->Get(name->Data());
3400   CoulCorrOmega0OS->SetDirectory(0);
3401   coulfile->Close();
3402
3403 }
3404 double CoulCorrGRS(bool SC, double Q_12, double Q_13, double Q_23){
3405   /*int Q12bin = CoulCorr2SS->GetXaxis()->FindBin(Q_12);
3406   int Q13bin = CoulCorr2SS->GetXaxis()->FindBin(Q_13);
3407   int Q23bin = CoulCorr2SS->GetXaxis()->FindBin(Q_23);*/
3408   int index12L = int(fabs(Q_12 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3409   int index12H = index12L+1;
3410   int index13L = int(fabs(Q_13 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3411   int index13H = index13L+1;
3412   int index23L = int(fabs(Q_23 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3413   int index23H = index23L+1;
3414
3415   if(Tricubic){
3416     // Tricubic Interpolation
3417     double arr[4][4][4]={{{0}}};
3418     for(int x=0; x<4; x++){
3419       for(int y=0; y<4; y++){
3420         for(int z=0; z<4; z++){
3421           if(SC){
3422             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2SS->GetBinContent((index23L)+y)*CoulCorr2SS->GetBinContent((index13L)+z);
3423           }else{
3424             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2OS->GetBinContent((index23L)+y)*CoulCorr2OS->GetBinContent((index13L)+z);
3425           }
3426
3427         }
3428       }
3429     }
3430     return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3431   }else{
3432     // Trilinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3433     //
3434     double xd = (Q_12-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3435     xd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3436     double yd = (Q_13-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3437     yd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index13H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3438     double zd = (Q_23-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3439     zd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index23H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3440     //
3441     if(SC){
3442       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3443       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3444       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3445       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3446       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3447       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3448       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3449       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3450       //
3451       double c0 = c00*(1-yd) + c10*yd;
3452       double c1 = c01*(1-yd) + c11*yd;
3453       return (c0*(1-zd) + c1*zd);
3454     }else{
3455       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3456       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3457       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3458       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3459       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3460       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3461       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3462       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3463       //
3464       double c0 = c00*(1-yd) + c10*yd;
3465       double c1 = c01*(1-yd) + c11*yd;
3466       return (c0*(1-zd) + c1*zd);
3467     }
3468   }
3469
3470 }
3471 double CoulCorrOmega0(bool SC, double Q_12, double Q_13, double Q_23){// 12 is same-charge pair
3472   int Q12bin = CoulCorrOmega0SS->GetXaxis()->FindBin(Q_12);
3473   int Q13bin = CoulCorrOmega0SS->GetZaxis()->FindBin(Q_13);
3474   int Q23bin = CoulCorrOmega0SS->GetYaxis()->FindBin(Q_23);
3475   int index12L = int(fabs(Q_12 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3476   int index12H = index12L+1;
3477   int index13L = int(fabs(Q_13 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3478   int index13H = index13L+1;
3479   int index23L = int(fabs(Q_23 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3480   int index23H = index23L+1;
3481
3482
3483   if(SC) {if(CoulCorrOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic same-charge Omega0 bin"<<endl;}}
3484   else {if(CoulCorrOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic mixed-charge Omega0 bin"<<endl;}}
3485
3486   if(Tricubic){
3487     //////////////////////////
3488     // Tricubic Interpolation
3489     double arr[4][4][4]={{{0}}};
3490     for(int x=0; x<4; x++){
3491       for(int y=0; y<4; y++){
3492         for(int z=0; z<4; z++){
3493           if(SC) arr[x][y][z] = CoulCorrOmega0SS->GetBinContent((index12L)+x, (index23L)+y, (index13L)+z);
3494           else arr[x][y][z] = CoulCorrOmega0OS->GetBinContent((Q12bin)+x, (Q23bin)+y, (Q13bin)+z);
3495         }
3496       }
3497     }
3498     return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3499   }else{
3500     ///////////////////////////
3501     // Trilinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3502     //
3503     double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3504     xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3505     double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3506     yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3507     double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3508     zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3509
3510
3511     double c00=0, c10=0, c01=0, c11=0;
3512     if(SC){
3513       // interpolate along x
3514       if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
3515       if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
3516       if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
3517       if(CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0SS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0SS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
3518     }else {
3519       if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)>0) c00 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13L+1)*xd;
3520       if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)>0) c10 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13L+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13L+1)*xd;
3521       if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)>0) c01 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23L+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23L+1, index13H+1)*xd;
3522       if(CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)>0 && CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)>0) c11 = CoulCorrOmega0OS->GetBinContent(index12L+1, index23H+1, index13H+1)*(1-xd) + CoulCorrOmega0OS->GetBinContent(index12H+1, index23H+1, index13H+1)*xd;
3523     }
3524     //
3525     double c0=0, c1=0;
3526     if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
3527     if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
3528
3529     if(c0>0 && c1>0) {
3530       double valueOm = (c0*(1-zd) + c1*zd);
3531       if(SC) valueOm *= StrongSC->Eval(Q_12*1000/2.) * StrongSC->Eval(Q_23*1000/2.) * StrongSC->Eval(Q_13*1000/2.);
3532       else valueOm *= StrongSC->Eval(Q_12*1000/2.);
3533       return valueOm;
3534     }else {
3535       cout<<"Not all Omega0 bins well defined.  Use GRS instead"<<endl;
3536       return CoulCorrGRS(kTRUE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
3537     }
3538   }
3539 }
3540
3541 double Gamov(int chargeProduct, double qinv){
3542
3543   double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
3544
3545   return arg/(exp(arg)-1);
3546 }
3547 /*double LednickyCorr(double qinv){
3548
3549   int indexL = int(fabs(1000*qinv-Lednicky_qinv[0])/(2*2.027027));// starts from 2.027 MeV
3550   int indexH = indexL+1;
3551   if( (indexL >= 73)) return 1.0;
3552   double slope = (Lednicky_CoulStrong[indexL] - Lednicky_CoulStrong[indexH])/(Lednicky_qinv[indexL]-Lednicky_qinv[indexH]);
3553   return (slope*(1000*qinv - Lednicky_qinv[indexL]) + Lednicky_CoulStrong[indexL]);
3554
3555 }
3556 void ReadLednickyFile(int Rvalue){
3557
3558   TString *filename=new TString("../Strong_FSI/converted_CoulStrong_");
3559   *filename += Rvalue;
3560   filename->Append("fm.root");
3561   TFile *Ledfile = new TFile(filename->Data(),"READ");
3562   TH1F *LedHisto = (TH1F*)Ledfile->Get("h27");
3563   for(int i=0; i<74; i++){
3564     Lednicky_qinv[i] = 2.*(LedHisto->GetXaxis()->GetBinLowEdge(i+1) + LedHisto->GetXaxis()->GetBinUpEdge(i+1))/2.;
3565     Lednicky_CoulStrong[i] = LedHisto->GetBinContent(i+1);
3566   }
3567   Ledfile->Close();
3568   }*/
3569 void ReadMomResFile(int rvalue, double lambda){
3570
3571   for(int i=0; i<40; i++){
3572     MomRes_C2_pp[i] = 1.;
3573     MomRes_C2_mp[i] = 1.;
3574     MomRes_term1_pp[i] = 1.;
3575     MomRes_term2_pp[i] = 1.;
3576     MomRes_term1_mp[i] = 1.;
3577     MomRes_term2_mp[i] = 1.;
3578     //
3579     MomResDev_C2_pp[i] = 1.;
3580     MomResDev_C2_mp[i] = 1.;
3581     MomResDev_term1_pp[i] = 1.;
3582     MomResDev_term2_pp[i] = 1.;
3583   }
3584
3585   TString *momresfilename = new TString("MomResFile");
3586   momresfilename->Append("_Offline_");
3587   if(FileSetting==3) momresfilename->Append("PID1p5");
3588   else if(FileSetting==9) momresfilename->Append("FB5and7overlap");
3589   else momresfilename->Append("11h");
3590   momresfilename->Append(".root");// no corresponding file for 10h
3591
3592   TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
3593
3594   TH2D *MomResWeights_pp = (TH2D*)MomResFile->Get("MomResHisto_pp");
3595   TH2D *MomResWeights_mp = (TH2D*)MomResFile->Get("MomResHisto_mp");
3596   TH2D *MomResWeights_pp_term1 = (TH2D*)MomResFile->Get("MomResHisto_pp_term1");
3597   TH2D *MomResWeights_pp_term2 = (TH2D*)MomResFile->Get("MomResHisto_pp_term2");
3598   TH2D *MomResWeights_mp_term1 = (TH2D*)MomResFile->Get("MomResHisto_mp_term1");
3599   TH2D *MomResWeights_mp_term2 = (TH2D*)MomResFile->Get("MomResHisto_mp_term2");
3600   //
3601   //
3602   TString *names3[2][5];// SC/MC, term#
3603   TString *names1D[2][5];// SC/MC, term#
3604   TString *names3_AvgK3[2];// SC/MC
3605   for(int ChProd=0; ChProd<2; ChProd++){
3606     if(ChProd==0) names3_AvgK3[ChProd] = new TString("AvgK3_SC_M");
3607     else names3_AvgK3[ChProd] = new TString("AvgK3_MC_M");
3608     *names3_AvgK3[ChProd] += MomResCentBin-1;
3609     //AvgK3[ChProd] = (TH3D*)MomResFile->Get(names3_AvgK3[ChProd]->Data());
3610     //AvgK3[ChProd]->SetDirectory(0);
3611     //
3612     for(int term=0; term<5; term++){
3613       //
3614       if(ChProd==0) {names3[ChProd][term] = new TString("MomResHisto3_SC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_SC_term");}
3615       else {names3[ChProd][term] = new TString("MomResHisto3_MC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_MC_term");}
3616       *names3[ChProd][term] += term+1;
3617       *names1D[ChProd][term] += term+1;
3618       names3[ChProd][term]->Append("_M");
3619       names1D[ChProd][term]->Append("_M");
3620       *names3[ChProd][term] += MomResCentBin-1;
3621       *names1D[ChProd][term] += MomResCentBin-1;
3622       MomRes3d[ChProd][term] = (TH3D*)MomResFile->Get(names3[ChProd][term]->Data());
3623       MomRes1d[ChProd][term] = (TH1D*)MomResFile->Get(names1D[ChProd][term]->Data());
3624       MomRes3d[ChProd][term]->SetDirectory(0);
3625       MomRes1d[ChProd][term]->SetDirectory(0);
3626
3627     }
3628   }
3629
3630
3631   Int_t LambdaRbin;
3632   if(rvalue<=5) LambdaRbin = 1 + int(float(lambda-0.5+0.0001)/0.02);
3633   else LambdaRbin = 1 + 16*(rvalue-5) + int(float(lambda-0.5+0.0001)/0.02);
3634   Int_t LambdaRbinDev;
3635   if(rvalue<=5) LambdaRbinDev = 1 + 16*(1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3636   else LambdaRbinDev = 1 + 16*(rvalue-5 - 1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3637
3638   for(int i=0; i<40; i++){
3639     MomRes_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbin, i+1);
3640     MomRes_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbin, i+1);
3641     MomRes_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbin, i+1);
3642     MomRes_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbin, i+1);
3643     MomRes_term1_mp[i] = MomResWeights_mp_term1->GetBinContent(LambdaRbin, i+1);
3644     MomRes_term2_mp[i] = MomResWeights_mp_term2->GetBinContent(LambdaRbin, i+1);
3645     //
3646     MomResDev_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbinDev, i+1);
3647     MomResDev_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbinDev, i+1);
3648     MomResDev_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbinDev, i+1);
3649     MomResDev_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbinDev, i+1);
3650   }
3651
3652   MomResFile->Close();
3653
3654 }
3655
3656 void fcn_r3_3d(int& npar, double* deriv, double& f, double par[], int flag){
3657
3658   double r3=0;
3659   double Chi2sum=0;
3660
3661   for(int i=fitlimitLowBin-1; i<fitlimitHighBin; i++){
3662     double q12=(i+0.5)*0.005;
3663     for(int j=fitlimitLowBin-1; j<fitlimitHighBin; j++){
3664       double q13=(j+0.5)*0.005;
3665       for(int k=fitlimitLowBin-1; k<fitlimitHighBin; k++){
3666         double q23=(k+0.5)*0.005;
3667
3668         //r3 = 1 - par[2]*pow(q12*q13 + q12*q23 + q13*q23,2) + 1/5.*pow(par[0],2)*(3-par[1]*(q12*q12+q13*q13+q23*q23)) - 112/70.*pow(par[0],3);
3669         //r3 /= pow( (1-par[1]*q12*q12 - 4/5.*pow(par[0],2))*(1-par[1]*q13*q13 - 4/5.*pow(par[0],2))*(1-par[1]*q23*q23 - 4/5.*pow(par[0],2)),0.5);
3670         r3 = par[0] - par[1]*(q12*q13 + q12*q23 + q13*q23) - par[2]*pow(q12*q13 + q12*q23 + q13*q23,2);
3671         if(r3_3d_arr_e[i][j][k] > 0) Chi2sum += pow((r3_3d_arr[i][j][k] - r3)/r3_3d_arr_e[i][j][k],2);
3672
3673       }
3674     }
3675   }
3676
3677   f = Chi2sum;
3678 }
3679
3680 double r3_fitfunction(Double_t *x, Double_t *par){
3681
3682   //double r3 = 1 - par[2]*pow(x[0]*x[1] + x[0]*x[2] + x[1]*x[2],2) + 1/5.*pow(par[0],2)*(3-par[1]*(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])) - 112/70.*pow(par[0],3);
3683   //r3 /= pow( (1-par[1]*x[0]*x[0] - 4/5.*pow(par[0],2))*(1-par[1]*x[1]*x[1] - 4/5.*pow(par[0],2))*(1-par[1]*x[2]*x[2] - 4/5.*pow(par[0],2)),0.5);
3684   double r3 = par[0] - par[1]*(x[0]*x[1] + x[0]*x[2] + x[1]*x[2]) - par[2]*pow(x[0]*x[1] + x[0]*x[2] + x[1]*x[2],2);
3685
3686   return r3;
3687
3688 }
3689
3690 double cubicInterpolate (double p[4], double x) {
3691   return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));// Paulinternet
3692 }
3693 double bicubicInterpolate (double p[4][4], double x, double y) {
3694         double arr[4];
3695         arr[0] = cubicInterpolate(p[0], y);
3696         arr[1] = cubicInterpolate(p[1], y);
3697         arr[2] = cubicInterpolate(p[2], y);
3698         arr[3] = cubicInterpolate(p[3], y);
3699         return cubicInterpolate(arr, x);
3700 }
3701 double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
3702         double arr[4];
3703         arr[0] = bicubicInterpolate(p[0], y, z);
3704         arr[1] = bicubicInterpolate(p[1], y, z);
3705         arr[2] = bicubicInterpolate(p[2], y, z);
3706         arr[3] = bicubicInterpolate(p[3], y, z);
3707         return cubicInterpolate(arr, x);
3708 }
3709 float fact(float n){
3710   return (n < 1.00001) ? 1 : fact(n - 1) * n;
3711 }
3712 void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2)
3713 {
3714   // revision on July 24th or 25th 2012
3715   // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
3716   x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
3717   // Printf("%f %f %f %f", x1, x2, y1, y2);
3718
3719   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
3720   myPadLogo->SetLeftMargin(0);
3721   myPadLogo->SetTopMargin(0);
3722   myPadLogo->SetRightMargin(0);
3723   myPadLogo->SetBottomMargin(0);
3724   myPadLogo->Draw();
3725   myPadLogo->cd();
3726   TASImage *myAliceLogo = new TASImage((prel) ? "~/Pictures/2011-Nov-24-ALICE_PRELIMARY_logo_BLACK_small_usage_design.eps" : "~/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.eps");
3727   myAliceLogo->Draw();
3728 }
3729 //________________________________________________________________________________________
3730 void ReadMuonCorrections(int mbin){
3731   TString *name = new TString("MuonCorrection_");
3732   if(mbin<=1) *name += 0;
3733   else if(mbin<=3) *name += 1;
3734   else if(mbin<=5) *name += 2;
3735   else if(mbin<=7) *name += 3;
3736   else if(mbin<=9) *name += 3;
3737   else *name += 4;
3738   name->Append(".root");
3739   TFile *muonfile=new TFile(name->Data(),"READ");
3740   C2muonCorrection_SC = (TH1D*)muonfile->Get("C2muonCorrection_SC");
3741   C2muonCorrection_MC = (TH1D*)muonfile->Get("C2muonCorrection_MC");
3742   WeightmuonCorrection = (TH1D*)muonfile->Get("WeightmuonCorrection");
3743   if(SameCharge_def) C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_SC");
3744   else C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_MC");
3745   //
3746   C2muonCorrection_SC->SetDirectory(0);
3747   C2muonCorrection_MC->SetDirectory(0);
3748   WeightmuonCorrection->SetDirectory(0);
3749   C3muonCorrection->SetDirectory(0);
3750   //
3751   muonfile->Close();
3752 }