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
  29 #define BohrR 1963.6885 // Bohr Radius for pions
  30 #define FmToGeV 0.19733 // conversion to fm
  31 #define PI 3.1415926
  32 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  33
  34 using namespace std;
  35
  36 int FileSetting=0;// 0(standard), 1(r3 Lambda), 2(TTC), 3(PID), 4(Pos B), 5(Neg B), 6(Gauss), 7(4 kt bins old TTC cuts)
  37 bool MCcase_def=kFALSE;// MC data?
  38 int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case
  39 bool SameCharge_def=kTRUE;// 3-pion same-charge?
  40 int EDbin_def=0;// 0 or 1: Kt3 bin
  41 //
  42 bool IncludeG_def=kFALSE;// Include coherence?
  43 bool IncludeEW_def=kTRUE;// Include EdgeWorth coefficients?
  44 bool IncludeEWfromTherm_def=kFALSE;// Include EdgeWorth coefficients from Therminator?
  45 bool GRS_def=kTRUE;// Generalized RiverSide 3-body FSI correction?
  46 int Mbin_def=0;// 0-9: centrality bin in widths of 5%
  47 int Ktbin_def=10;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
  48 double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase
  49 //
  50 //
  51 bool IncludeMJcorrection=kTRUE;// linear Mini-Jet correction for denominator of r3?
  52 bool SaveToFile_def=kFALSE;// Save outputs to file?
  53 int SourceType=0;// 0=Therminator, 1=Gaussian (keep at 0 for default)
  54 bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin?
  55 bool GofP=kFALSE;// Include momentum dependence of coherent fraction?
  56 bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states?
  57 bool LinkN=kFALSE;// Make N(++)=N(+-)?
  58 bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
  59 bool Tricubic=kFALSE;// Tricubic or Trilinear interpolation?  Tricubic was found to be unstable.
  60 bool IncludeSS=kTRUE;// Include same-charge two-pion correlations in the fit?
  61 bool IncludeOS=kTRUE;// Include mixed-charge two-pion correlations in the fit?
  62 //
  63 //
  64 //
  65 //
  66 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)
  67 const int Sbins_3=1;// 3-particle Species bins. 1=Pion-Pion-Pion only. max=10
  68 const int fitlimitLowBin = 3;
  69 const int fitlimitHighBin = 14;// max = 20 (100MeV)
  70 bool LHC11h=kTRUE;// always kTRUE
  71 bool ExplicitLoop=kFALSE;// always kFALSE
  72 bool PairCut=kTRUE;// always kTRUE
  73 bool PbPbcase=kTRUE;// always kTRUE
  74 int KtbinFSI;
  75 int Ktlowbin;
  76 int Kthighbin;
  77 int bBin;// impact parameter bin
  78 int MomResCentBin;// momentum resolution centrality bin
  79 int RValue;// Radius setting for Gaussian source profile
  80 int bValue;// impact parameter for Therminator source profile (default)
  81 float TwoFrac;// Lambda parameter
  82 float OneFrac;// Lambda^{1/2}
  83 float ThreeFrac;// Lambda^{3/2}
  84 float TwoFracMomRes;// Lambda parameter for momentum resolution corrections
  85 float RValueMomRes;// Gaussian radius value for momentum resolution corrections
  86 float TwoFracr3;// Lambda parameter for r3
  87 int TPN_bin;// Two Particle Normalization bin (always 0)
  88 double Qcut_pp = 0.6;// 0.6
  89 double Qcut_PbPb = 0.1;
  90 double NormQcutLow_pp = 1.0;
  91 double NormQcutHigh_pp = 1.5;
  92 double NormQcutLow_PbPb = .15;//1.06
  93 double NormQcutHigh_PbPb = .175;//1.1
  94
  95 // single-pion emitter size (for fits with momentum dependence of coherence)
  96 double RT = 0.72/FmToGeV;
  97 double Temp = 0.135;// 0.135 GeV
  98
  99
 100 const int Nlines = 50;
 101 TH3D *CoulCorrOmega0SS;
 102 TH3D *CoulCorrOmega0OS;
 103 TH1D *CoulCorr2SS;
 104 TH1D *CoulCorr2OS;
 105 //
 106 //static double Lednicky_qinv[74];
 107 //static double Lednicky_CoulStrong[74];
 108
 109 void ReadCoulCorrections(int, int, int);
 110 void ReadCoulCorrections_Omega0();
 111 //void ReadLednickyFile(int);
 112 void ReadMomResFile(int, double);
 113 double CoulCorr2(int, double);
 114 double CoulCorrOmega0(bool, double, double, double);
 115 double CoulCorrGRS(bool, double, double, double);
 116 double OSLfitfunction(Double_t *x, Double_t *par);
 117 double D3fitfunction_3(Double_t *x, Double_t *par);
 118 double r3_fitfunction(Double_t *x, Double_t *par);
 119 double Gamov(int, double);
 120 //double LednickyCorr(double);
 121 double C2ssFitFunction(Double_t *x, Double_t *par);
 122 double C2osFitFunction(Double_t *x, Double_t *par);
 123 double cubicInterpolate(double[4], double);
 124 double bicubicInterpolate(double[4][4], double, double);
 125 double tricubicInterpolate(double[4][4][4], double, double, double);
 126 //
 127 void fcnC2_Global(int&, double*, double&, double[], int);
 128 void fcnOSL(int&, double*, double&, double[], int);
 129
 130 const int BINRANGE_C2global=20;
 131 double C2ssFitting[BINRANGE_C2global];
 132 double C2osFitting[BINRANGE_C2global];
 133 double C2ssFitting_e[BINRANGE_C2global];
 134 double C2osFitting_e[BINRANGE_C2global];
 135 double K2SS[BINRANGE_C2global];
 136 double K2OS[BINRANGE_C2global];
 137 double K2SS_e[BINRANGE_C2global];
 138 double K2OS_e[BINRANGE_C2global];
 139 double K3SS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
 140 double K3OS_e[BINRANGE_C2global][BINRANGE_C2global][BINRANGE_C2global];
 141 double Chi2_C2global;
 142 double NFitPoints_C2global;
 143 double Chi2_C3global;
 144 double NFitPoints_C3global;
 145
 146 const int BINS_OSL=40;// out,side,long bins
 147 double A[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 148 double B[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 149 double A_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 150 double B_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 151 double avg_q[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 152 double K_OSL[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 153 double K_OSL_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 154 double Chi2_OSL;
 155 double NFitPoints_OSL;
 156
 157 const int BINRANGE_3=50;// q12,q13,q23
 158 double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 159 double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 160
 161 double C3_base[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 162 double N_term1[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 163 double N_term2[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 164 double N_term3[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 165 double N_term4[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 166 double N_term5[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 167
 168 double MomRes_C2_pp[40];// Qinv bins
 169 double MomRes_C2_mp[40];//
 170 double MomRes_term1_pp[40];
 171 double MomRes_term2_pp[40];
 172 double MomRes_term1_mp[40];
 173 double MomRes_term2_mp[40];
 174
 175 TH3D *MomRes3d[2][5];// SC/MC, term#
 176 TH1D *MomRes1d[2][5];// SC/MC, term#
 177 TH3D *MomRes3d_FVP[2][5];// sum#, term#
 178 TH3D *AvgK3[2];// SC/MC
 179 double AvgP[6][20];// kt bin, qinv bin
 180 int Ktbin_GofP;
 181 //
 182 double MomResDev_C2_pp[40];// Qinv bins
 183 double MomResDev_C2_mp[40];
 184 double MomResDev_term1_pp[40];
 185 double MomResDev_term2_pp[40];
 186
 187
 188 double r3_3d_arr[20][20][20];
 189 double r3_3d_arr_e[20][20][20];
 190 double AvgQinvSS[30];
 191 double AvgQinvOS[30];
 192 //
 193
 194
 195 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){
 196
 197   EDbin_def=EDbin;
 198   SaveToFile_def=SaveToFile;
 199   MCcase_def=MCcase;
 200   IncludeEWfromTherm_def=IncludeEWfromTherm;
 201   CHARGE_def=CHARGE;
 202   IncludeG_def=IncludeG;
 203   IncludeEW_def=IncludeEW;
 204   SameCharge_def=SameCharge;// 3-pion same-charge
 205   Mbin_def=Mbin;
 206   Ktbin_def=Ktbin;
 207   if(Ktbin_def==10) Ktbin_GofP=2;
 208   else Ktbin_GofP=Ktbin_def;
 209   //
 210   Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
 211   Kthighbin=(Ktbin)*2+4;
 212
 213   //
 214   if(ConstantFSI) KtbinFSI=10;
 215   else KtbinFSI=Ktbin;
 216
 217
 218   TwoFracr3 = 0.7;// standard lambda is 0.7 for r3
 219   if(FileSetting==1) TwoFracr3 = 0.6;
 220
 221   if(Mbin == 0) {// 0-5%
 222     RValue=8;
 223     TwoFrac=0.68;
 224     TwoFracMomRes=0.68;
 225     RValueMomRes=10;// for C2
 226     bValue=2;
 227     bBin=1;
 228     MomResCentBin=1;// for C3
 229   }else if(Mbin == 1) {// 5-10%
 230     RValue=8;
 231     TwoFrac=0.68;
 232     TwoFracMomRes=0.68;
 233     RValueMomRes=9;
 234     bValue=3;
 235     bBin=2;
 236     MomResCentBin=2;
 237   }else if(Mbin <= 3) {// 10-20%
 238     RValue=7;
 239     TwoFrac=0.68;
 240     TwoFracMomRes=0.68;
 241     RValueMomRes=8;
 242     bValue=5;
 243     bBin=3;
 244     MomResCentBin=3;
 245   }else if(Mbin <= 5) {// 20-30%
 246     RValue=6;
 247     TwoFrac=0.68;
 248     TwoFracMomRes=0.68;
 249     RValueMomRes=7;
 250     bValue=7;
 251     bBin=4;
 252     MomResCentBin=4;
 253   }else if(Mbin <= 7){// 30-40%
 254     RValue=6;
 255     TwoFrac=0.68;
 256     TwoFracMomRes=0.68;
 257     RValueMomRes=6;
 258     bValue=8;
 259     bBin=5;
 260     MomResCentBin=5;
 261   }else {// 40-50%
 262     RValue=5;
 263     TwoFrac=0.68;
 264     TwoFracMomRes=0.68;
 265     RValueMomRes=5;
 266     bValue=9;
 267     bBin=6;
 268     MomResCentBin=6;
 269   }
 270
 271   OneFrac = sqrt(TwoFrac);
 272   ThreeFrac = pow(TwoFrac,3/2.);
 273
 274
 275   if(SourceType==1 && RValue > 8) {cout<<"Radius value too large!!!"<<endl; return;}
 276
 277   cout<<"Mbin = "<<Mbin<<"   Kt = "<<Ktbin<<"   R input = "<<RValue<<"   lambda input = "<<TwoFrac<<endl;
 278
 279   // Core-Halo modeling of single-particle and triple-particle core fraction
 280   //float AvgN[10]={472.56, 390.824, 319.597, 265.933, 218.308, 178.865, 141.2, 115.88, 87.5556, 69.3235};// (avg total pion mult)/2.  2.0sigma
 281   /*
 282   float fc = (1+sqrt(1 + 4*AvgN[Mbin]*TwoFrac*(AvgN[Mbin]-1)))/(2*AvgN[Mbin]);// Two component model (core + non-interacting halo)
 283   cout<<"Before: "<<OneFrac<<"  "<<ThreeFrac<<endl;
 284   OneFrac = fc;
 285   ThreeFrac = (fc*AvgN[Mbin]*(fc*AvgN[Mbin]-1)*(fc*AvgN[Mbin]-2))/(AvgN[Mbin]*(AvgN[Mbin]-1)*(AvgN[Mbin]-2));
 286   cout<<"After:  "<<OneFrac<<"  "<<ThreeFrac<<"  "<<endl;
 287   */
 288   //
 289   // avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values.  From unsmeared HIJING 0-5% with input signal
 290   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};
 291   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};
 292   for(int i=0; i<30; i++){
 293     AvgQinvSS[i] = AvgQinvSS_temp[i];
 294     AvgQinvOS[i] = AvgQinvOS_temp[i];
 295   }
 296
 297
 298   // average total momentum in each kt and qinv bin.  Used for a test of momentum dependence of fits including coherence
 299   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};
 300   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};
 301   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};
 302   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};
 303   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};
 304   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};
 305   for(int ktbin=1; ktbin<=6; ktbin++){
 306     for(int qbin=1; qbin<=20; qbin++){
 307       if(ktbin==1) AvgP[ktbin-1][qbin-1] = AvgP_kt1[qbin-1];
 308       else if(ktbin==2) AvgP[ktbin-1][qbin-1] = AvgP_kt2[qbin-1];
 309       else if(ktbin==3) AvgP[ktbin-1][qbin-1] = AvgP_kt3[qbin-1];
 310       else if(ktbin==4) AvgP[ktbin-1][qbin-1] = AvgP_kt4[qbin-1];
 311       else if(ktbin==5) AvgP[ktbin-1][qbin-1] = AvgP_kt5[qbin-1];
 312       else AvgP[ktbin-1][qbin-1] = AvgP_kt6[qbin-1];
 313     }
 314   }
 315
 316
 317
 318   gStyle->SetOptStat(0);
 319   gStyle->SetOptDate(0);
 320   gStyle->SetOptFit(0111);
 321
 322   TFile *myfile;
 323   if(PbPbcase){// PbPb
 324     if(MCcase) {
 325       if(Mbin<=1){
 326         myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11.root","READ");
 327         //myfile = new TFile("Results/PDC_HIJING_12a17ad_fix_genSignal_Rinv11_Smeared.root","READ");
 328       }else{
 329         myfile = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p68R11_C2plots.root","READ");
 330       }
 331     }else{
 332       //myfile = new TFile("Results/PDC_11h_v1paper.root","READ");// v1 paper
 333       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// Old standard (bad r3 Interp, 0.035TTC)
 334       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// v3 Standard
 335       //
 336       //if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins_FullTPC.root","READ");// FullTPC runlist
 337       if(FileSetting==0) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// Standard
 338       if(FileSetting==1) myfile = new TFile("Results/PDC_11h_Lam0p7_and_Lam0p6.root","READ");// Lam=0.7 and Lam=0.6 (3ktbins, 0.035TTC)
 339       if(FileSetting==2) myfile = new TFile("Results/PDC_11h_2sigmaTTC_3sigmaTTC.root","READ");// TTC variation
 340       if(FileSetting==3) myfile = new TFile("Results/PDC_11h_PID1p5.root","READ");// PID variation (0.035TTC)
 341       if(FileSetting==4) myfile = new TFile("Results/PDC_11h_PosB.root","READ");// Positive B field for r3num (0.035TTC)
 342       if(FileSetting==5) myfile = new TFile("Results/PDC_11h_NegB.root","READ");// Negative B field for r3num (0.035TTC)
 343       if(FileSetting==6) myfile = new TFile("Results/PDC_11h_3sigmaTTCDenextrap_GaussDenextrap.root","READ");// Gaussian, 3sigmaTTC)
 344       if(FileSetting==7) myfile = new TFile("Results/PDC_11h_4ktbins.root","READ");// 4 kt bins (0.035TTC)
 345       if(FileSetting==8) myfile = new TFile("Results/PDC_11h_2Kt3bins.root","READ");// 2 Kt3 bins
 346     }
 347   }else{// pp
 348     cout<<"pp not currently supported"<<endl;
 349     return;
 350   }
 351
 352   ReadCoulCorrections(RValue, bValue, KtbinFSI);
 353   //ReadLednickyFile(RValue);
 354   ReadMomResFile(RValueMomRes, TwoFracMomRes);
 355   ReadCoulCorrections_Omega0();
 356   //
 357   /////////////////////////////////////////////////////////
 358
 359
 360   double NormQcutLow;
 361   double NormQcutHigh;
 362   if(PbPbcase) {
 363     NormQcutLow = NormQcutLow_PbPb;
 364     NormQcutHigh = NormQcutHigh_PbPb;
 365   }else {
 366     NormQcutLow = NormQcutLow_pp;
 367     NormQcutHigh = NormQcutHigh_pp;
 368   }
 369
 370
 371   TList *MyList;
 372   if(!MCcase){
 373     TDirectoryFile *tdir = (TDirectoryFile*)myfile->Get("PWGCF.outputChaoticityAnalysis.root");
 374     if(FileSetting!=1 && FileSetting !=6) MyList=(TList*)tdir->Get("ChaoticityOutput_1");
 375     else MyList=(TList*)tdir->Get("ChaoticityOutput_2");
 376   }else{
 377     MyList=(TList*)myfile->Get("MyList");
 378   }
 379   myfile->Close();
 380
 381
 382   TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
 383   //
 384
 385   cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
 386
 387
 388   TH1D *ChiSquaredNDF = new TH1D("ChiSquaredNDF","",2,0.5,2.5);// Chi^2/NDF records
 389
 390   // Explicit Loop Histos
 391   TH2D *Two_ex_2d[2][2][Sbins_2][2];
 392   TH1D *Two_ex[2][2][Sbins_2][2];
 393
 394   // Normalizations
 395   double NormH_pc[5]={0};
 396   TH3D *PionNorm[2];
 397   //TH3D *PionNorm_e[2];
 398
 399   double norm_ex_2[6][2]={{0}};
 400
 401   // 3d method histograms
 402   TH3D *Three_3d[2][2][2][Sbins_3][5];
 403   // 4-vect Product Sum
 404   TH3D *TPNRejects;
 405   ///////////////////////////////////
 406   // Create Histograms
 407
 408   // 2-particle
 409   for(int c1=0; c1<2; c1++){
 410     for(int c2=0; c2<2; c2++){
 411       for(int sc=0; sc<Sbins_2; sc++){
 412         for(int term=0; term<2; term++){
 413
 414           TString *name2_ex = new TString("Explicit2_Charge1_");
 415           *name2_ex += c1;
 416           name2_ex->Append("_Charge2_");
 417           *name2_ex += c2;
 418           name2_ex->Append("_SC_");
 419           *name2_ex += sc;
 420           name2_ex->Append("_M_");
 421           *name2_ex += Mbin;
 422           name2_ex->Append("_ED_");
 423           *name2_ex += 0;// EDbin
 424           name2_ex->Append("_Term_");
 425           *name2_ex += term+1;
 426
 427           if(sc==0 || sc==3 || sc==5){
 428             if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
 429           }
 430
 431
 432
 433           Two_ex_2d[c1][c2][sc][term] = (TH2D*)MyList->FindObject(name2_ex->Data());
 434           Two_ex_2d[c1][c2][sc][term]->Sumw2();
 435           TString *proname = new TString(name2_ex->Data());
 436           proname->Append("_pro");
 437
 438           if(Ktbin==10) {Ktlowbin=1; Kthighbin=Two_ex_2d[c1][c2][sc][term]->GetNbinsX();}
 439           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)
 440
 441           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));
 442           Two_ex[c1][c2][sc][term]->Scale(norm_ex_2[sc][0]/norm_ex_2[sc][term]);
 443
 444           Two_ex[c1][c2][sc][term]->SetMarkerStyle(20);
 445           Two_ex[c1][c2][sc][term]->GetXaxis()->SetTitle("q_{inv}");
 446           Two_ex[c1][c2][sc][term]->GetYaxis()->SetTitle("C_{2}");
 447           Two_ex[c1][c2][sc][term]->SetTitle("");
 448
 449         }// term
 450       }// sc
 451
 452       // 3-particle
 453       for(int c3=0; c3<2; c3++){
 454         for(int sc=0; sc<Sbins_3; sc++){
 455           for(int term=0; term<5; term++){
 456
 457             TString *name3_ex = new TString("Explicit3_Charge1_");
 458             *name3_ex += c1;
 459             name3_ex->Append("_Charge2_");
 460             *name3_ex += c2;
 461             name3_ex->Append("_Charge3_");
 462             *name3_ex += c3;
 463             name3_ex->Append("_SC_");
 464             *name3_ex += sc;
 465             name3_ex->Append("_M_");
 466             *name3_ex += Mbin;
 467             name3_ex->Append("_ED_");
 468             *name3_ex += EDbin;
 469             name3_ex->Append("_Term_");
 470             *name3_ex += term+1;
 471
 472
 473             TString *name3_pc = new TString("PairCut3_Charge1_");
 474             *name3_pc += c1;
 475             name3_pc->Append("_Charge2_");
 476             *name3_pc += c2;
 477             name3_pc->Append("_Charge3_");
 478             *name3_pc += c3;
 479             name3_pc->Append("_SC_");
 480             *name3_pc += sc;
 481             name3_pc->Append("_M_");
 482             *name3_pc += Mbin;
 483             name3_pc->Append("_ED_");
 484             *name3_pc += EDbin;
 485             name3_pc->Append("_Term_");
 486             *name3_pc += term+1;
 487
 488             ///////////////////////////////////////
 489             // skip degenerate histograms
 490             if(sc==0 || sc==6 || sc==9){// Identical species
 491               if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
 492               if( (c1+c2+c3)==2) {if(c1!=0) continue;}
 493             }else if(sc!=5){
 494               if( (c1+c2)==1) {if(c1!=0) continue;}
 495             }else {}// do nothing for pi-k-p case
 496
 497             /////////////////////////////////////////
 498
 499
 500             if(ExplicitLoop) {
 501               /*
 502                 Three_ex[c1][c2][c3][sc][term] = (TH1D*)MyList->FindObject(name3_ex->Data());
 503
 504                 Three_ex[c1][c2][c3][sc][term]->SetMarkerStyle(20);
 505                 Three_ex[c1][c2][c3][sc][term]->SetLineColor(2);
 506                 Three_ex[c1][c2][c3][sc][term]->SetMarkerColor(2);
 507                 Three_ex[c1][c2][c3][sc][term]->GetXaxis()->SetTitle("Q_{3}");
 508                 Three_ex[c1][c2][c3][sc][term]->GetYaxis()->SetTitle("C_{3}");
 509                 Three_ex[c1][c2][c3][sc][term]->SetTitle("");
 510
 511                 TString *name = new TString(name3_ex->Data());
 512                 name->Append("_Norm");
 513                 NormH_ex[term] = ((TH1D*)MyList->FindObject(name->Data()))->Integral();
 514
 515                 if(NormH_ex[term] > 0){
 516                   Three_ex[c1][c2][c3][sc][term]->Scale(NormH_ex[0]/NormH_ex[term]);
 517                 }else{ cout<<"normalization = 0.  Skipping this SC."<<endl;}
 518               */
 519             }
 520             if(PairCut){
 521
 522               TString *nameNorm=new TString("PairCut3_Charge1_");
 523               *nameNorm += c1;
 524               nameNorm->Append("_Charge2_");
 525               *nameNorm += c2;
 526               nameNorm->Append("_Charge3_");
 527               *nameNorm += c3;
 528               nameNorm->Append("_SC_");
 529               *nameNorm += sc;
 530               nameNorm->Append("_M_");
 531               *nameNorm += Mbin;
 532               nameNorm->Append("_ED_");
 533               *nameNorm += 0;// EDbin
 534               nameNorm->Append("_Term_");
 535               *nameNorm += term+1;
 536               nameNorm->Append("_Norm");
 537               NormH_pc[term] = ((TH1D*)MyList->FindObject(nameNorm->Data()))->Integral();
 538
 539               if(NormH_pc[term] > 0){
 540
 541                 if(sc<=2) {
 542
 543                   TString *name_3d = new TString(name3_pc->Data());
 544                   name_3d->Append("_3D");
 545                   Three_3d[c1][c2][c3][sc][term] = (TH3D*)MyList->FindObject(name_3d->Data());
 546                   Three_3d[c1][c2][c3][sc][term]->Sumw2();
 547                   Three_3d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]);
 548                   //cout<<"Scale factor = "<<NormH_pc[0]/NormH_pc[term]<<endl;
 549
 550                   // 4-vect Product Sum
 551                   if(c1==c2 && c1==c3 && sc==0 && term==4){
 552
 553                     TString *nameDenType=new TString("PairCut3_Charge1_");
 554                     *nameDenType += c1;
 555                     nameDenType->Append("_Charge2_");
 556                     *nameDenType += c2;
 557                     nameDenType->Append("_Charge3_");
 558                     *nameDenType += c3;
 559                     nameDenType->Append("_SC_");
 560                     *nameDenType += sc;
 561                     nameDenType->Append("_M_");
 562                     *nameDenType += Mbin;
 563                     nameDenType->Append("_ED_");
 564                     *nameDenType += EDbin;
 565                     nameDenType->Append("_TPN_");
 566                     *nameDenType += 0;
 567
 568                     PionNorm[c1] = (TH3D*)MyList->FindObject(nameDenType->Data());
 569                     PionNorm[c1]->Sumw2();
 570                     PionNorm[c1]->Scale(NormH_pc[0]/NormH_pc[term]);
 571                     if(c1==0) {
 572                       TPNRejects = (TH3D*)MyList->FindObject("fTPNRejects4");
 573                       TPNRejects->Scale(NormH_pc[0]/NormH_pc[term]);
 574                     }
 575
 576                   }
 577                 }
 578
 579               }else{
 580                 cout<<"normalization = 0.  Skipping this SC.  SC="<<sc<<"  c1="<<c1<<"  c2="<<c2<<"  c3="<<c3<<endl;
 581               }
 582
 583
 584
 585             }// pair cut
 586           }// term
 587
 588         }// SC
 589
 590
 591       }// c3
 592     }// c2
 593   }// c1
 594
 595
 596   cout<<"Done getting Histograms"<<endl;
 597
 598
 599   TCanvas *can1 = new TCanvas("can1", "can1",11,53,700,500);
 600   can1->SetHighLightColor(2);
 601   can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 602   gStyle->SetOptFit(0111);
 603   can1->SetFillColor(10);
 604   can1->SetBorderMode(0);
 605   can1->SetBorderSize(2);
 606   can1->SetGridx();
 607   can1->SetGridy();
 608   can1->SetFrameFillColor(0);
 609   can1->SetFrameBorderMode(0);
 610   can1->SetFrameBorderMode(0);
 611
 612   TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
 613   legend1->SetBorderSize(1);
 614   legend1->SetTextSize(.04);// small .03; large .036
 615   legend1->SetFillColor(0);
 616
 617   /////////////////////////////////////////////////////////
 618   /////////////////////////////////////////////////////////
 619   // This part for plotting track splitting/merging effects in MC data only
 620   /*
 621   TH3F *Merge3d_num=(TH3F*)MyList->FindObject("fPairsDetaDPhiNum");
 622   TH3F *Merge3d_den=(TH3F*)MyList->FindObject("fPairsDetaDPhiDen");
 623   //TH3F *Merge3d_num=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_num");
 624   //TH3F *Merge3d_den=(TH3F*)MyList->FindObject("Pairs_dEtadPhi_UNI_den");
 625
 626   TH1D *Merge1d_num[10];
 627   TH1D *Merge1d_den[10];
 628   TString *newnamenum[10];
 629   TString *newnameden[10];
 630   TF1 *MergedGaus=new TF1("MergedGaus","1-[0]*exp(-pow(x/[1],2))",-0.1, 0.1);
 631   MergedGaus->SetParName(0,"Amplitude");
 632   MergedGaus->SetParName(1,"width");
 633   MergedGaus->SetParameter(0,0.06);
 634   MergedGaus->SetParameter(1,0.01);
 635   MergedGaus->SetParLimits(0,0.001,0.5);
 636   MergedGaus->SetParLimits(1,0.001,0.1);
 637
 638   for(int i=2; i<10; i++){
 639     if(i!=5 && i!=8) continue;// 5 and 8
 640     newnamenum[i]=new TString("namenum_");
 641     *newnamenum[i] += i;
 642     newnameden[i]=new TString("nameden_");
 643     *newnameden[i] += i;
 644
 645     Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionZ(newnamenum[i]->Data(),i+1,i+1,90,110);//90,110 (phi)
 646     Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionZ(newnameden[i]->Data(),i+1,i+1,90,110);// (phi)
 647     //Merge1d_num[i]=(TH1D*)Merge3d_num->ProjectionY(newnamenum[i]->Data(),i+1,i+1,190,410);// 190,410 (eta)
 648     //Merge1d_den[i]=(TH1D*)Merge3d_den->ProjectionY(newnameden[i]->Data(),i+1,i+1,190,410);// (eta)
 649     //Merge1d_num[i]->Rebin(2);
 650     //Merge1d_den[i]->Rebin(2);
 651     Merge1d_num[i]->Sumw2();
 652     Merge1d_den[i]->Sumw2();
 653     Merge1d_num[i]->SetMarkerStyle(20);
 654
 655     if(Merge1d_den[i]->Integral(1,100)<=0) continue;
 656     double SF_merge = Merge1d_num[i]->Integral(1,100)/Merge1d_den[i]->Integral(1,100);// Z projection (phi)
 657     //double SF_merge = Merge1d_num[i]->Integral(1,50)/Merge1d_den[i]->Integral(1,50);// Y projection (eta)
 658     Merge1d_den[i]->Scale(SF_merge);
 659     Merge1d_num[i]->Divide(Merge1d_den[i]);
 660
 661
 662     if(i<9){
 663       Merge1d_num[i]->SetLineColor(i+1);
 664       Merge1d_num[i]->SetMarkerColor(i+1);
 665     }else{
 666       Merge1d_num[i]->SetLineColor(11);
 667       Merge1d_num[i]->SetMarkerColor(11);
 668     }
 669     if(i==4) {
 670       Merge1d_num[i]->SetLineColor(2);
 671       Merge1d_num[i]->SetMarkerColor(2);
 672     }
 673     if(i==5) {
 674       Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#phi*");
 675       //Merge1d_num[i]->GetXaxis()->SetTitle("#Delta#eta");
 676       Merge1d_num[i]->GetYaxis()->SetTitle("C_{2}^{HIJING}");
 677       Merge1d_num[i]->GetXaxis()->SetRangeUser(-.1,.1);
 678       Merge1d_num[i]->SetMinimum(.91);
 679       Merge1d_num[i]->SetMaximum(1.1);
 680       Merge1d_num[i]->DrawCopy();
 681
 682       //Merge1d_num[i]->Fit(MergedGaus,"IME","",-0.1,0.1);
 683     }else{
 684       Merge1d_num[i]->DrawCopy("same");
 685     }
 686
 687     TString *Dname=new TString("D=0.2*");
 688     *Dname +=i;
 689     Dname->Append(" m");
 690     legend1->AddEntry(newnamenum[i]->Data(),Dname->Data(),"lpf");
 691   }
 692   legend1->Draw("same");
 693   gStyle->SetOptFit(111);
 694   Merge1d_num[8]->Fit(MergedGaus,"IME","",-0.1,0.1);
 695   MergedGaus->Draw("same");
 696   */
 697   /*TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsPadRowNum");// kt, shfrac, qinv
 698   TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsPadRowDen");// kt, shfrac, qinv
 699   PadRowDen3->Scale(PadRowNum3->Integral(1,20,1,159, 35,40)/PadRowDen3->Integral(1,20,1,159, 35,40));
 700   PadRowNum3->GetYaxis()->SetRangeUser(0,0.01);
 701   PadRowDen3->GetYaxis()->SetRangeUser(0,0.01);
 702   TH1D *PadRowNum=(TH1D*)PadRowNum3->Project3D("z");
 703   TH1D *PadRowDen=(TH1D*)PadRowDen3->Project3D("z");
 704   PadRowNum->Divide(PadRowDen);
 705   PadRowNum->Draw();*/
 706   /*
 707   TH3D *PadRowNum3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiNum");// r, shfrac, deltaphi
 708   TH3D *PadRowDen3= (TH3D*)MyList->FindObject("fPairsShareFracDPhiDen");// r, shfrac, deltaphi
 709   PadRowDen3->Scale(PadRowNum3->Integral(1,10,1,159, 90,100)/PadRowDen3->Integral(1,10,1,159, 90,100));
 710   PadRowNum3->GetXaxis()->SetRange(5,5);
 711   PadRowDen3->GetXaxis()->SetRange(5,5);
 712   TH2D *PadRowNum=(TH2D*)PadRowNum3->Project3D("zy");
 713   TH2D *PadRowDen=(TH2D*)PadRowDen3->Project3D("zy");
 714   PadRowNum->Divide(PadRowDen);
 715   PadRowNum->Draw("lego");
 716   */
 717   /////////////////////////
 718   // 2-particle legend
 719   // 0 = pi-pi
 720   // 1 = pi-k
 721   // 2 = pi-p
 722   // 3 = k-k
 723   // 4 = k-p
 724   // 5 = p-p
 725   /////////////////////////
 726   TH1D *Two_pipi_mp = (TH1D*)Two_ex[0][1][0][0]->Clone();
 727   Two_pipi_mp->Divide(Two_ex[0][1][0][1]);
 728
 729   // Normalization range counting.  Just to evaluate required normalization width in qinv.
 730   //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.1))<<endl;
 731   //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;
 732   //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;
 733
 734
 735
 736   const int SCOI_2=0;
 737
 738   TH1D *Two_ex_clone_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 739   Two_ex_clone_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 740   TH1D *Two_ex_clone_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 741   Two_ex_clone_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 742   TH1D *Two_ex_clone_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 743   Two_ex_clone_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 744
 745   // Mini-jet ++ background linear estimation.
 746   TF1 *MJ_bkg_ss=new TF1("MJ_bkg_ss","pol1",0,1);
 747   Two_ex_clone_mm->Fit(MJ_bkg_ss,"IMENQ","",0.2,0.4);
 748   cout<<"Non-femto bkg C2(q=0.01) = "<<MJ_bkg_ss->Eval(0.01)<<endl;
 749
 750   Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}");
 751   Two_ex_clone_mm->SetTitle("");
 752   Two_ex_clone_mp->GetYaxis()->SetTitle("C_{2}");
 753   Two_ex_clone_mm->SetMarkerColor(2);
 754   Two_ex_clone_mm->SetLineColor(2);
 755   Two_ex_clone_mp->SetMarkerColor(1);
 756   Two_ex_clone_pp->SetMarkerColor(4);
 757   Two_ex_clone_pp->SetLineColor(4);
 758   Two_ex_clone_mm->GetXaxis()->SetRangeUser(0,0.1);
 759   Two_ex_clone_mm->SetMinimum(0.95);
 760   Two_ex_clone_mm->SetMaximum(1.4);
 761
 762   if(MCcase){
 763     Two_ex_clone_mp->SetMarkerColor(4);
 764     Two_ex_clone_mp->SetLineColor(4);
 765     Two_ex_clone_mm->Add(Two_ex_clone_pp);
 766     Two_ex_clone_mm->Scale(0.5);
 767     Two_ex_clone_mm->GetYaxis()->SetTitle("C_{2}^{HIJING}");
 768     Two_ex_clone_mm->GetYaxis()->SetTitleOffset(1.3);
 769     Two_ex_clone_mm->SetMinimum(0.97);
 770     Two_ex_clone_mm->SetMaximum(1.02);
 771     Two_ex_clone_mm->DrawCopy();
 772     //Two_ex_clone_pp->DrawCopy("same");
 773     Two_ex_clone_mp->DrawCopy("same");
 774     legend1->AddEntry(Two_ex_clone_mm,"same-charge","p");
 775     //legend1->AddEntry(Two_ex_clone_pp,"++","p");
 776     legend1->AddEntry(Two_ex_clone_mp,"mixed-charge","p");
 777     legend1->Draw("same");
 778   }
 779
 780
 781
 782   ////////////////////////////////////////////
 783   // Get Therminator EdgeWorth coefficients
 784   TString *ThermName = new TString("../ThermData/Therm_FSI_b");
 785   *ThermName += bValue;
 786   ThermName->Append(".root");
 787   TFile *Thermfile = new TFile(ThermName->Data(),"READ");
 788   TH2D *thermNum2d_ss=(TH2D*)Thermfile->Get("Num_CosFSI_ss");
 789   TH2D *thermNumSq2d_ss=(TH2D*)Thermfile->Get("NumSq_CosFSI_ss");
 790   TH2D *thermDen2d_ss=(TH2D*)Thermfile->Get("Den_ss");
 791   TH2D *thermLargeR2D_ss=(TH2D*)Thermfile->Get("LargeRpairs_ss");
 792   TH1D *thermNum_ss=(TH1D*)thermNum2d_ss->ProjectionY("thermNum_ss",Ktlowbin,Kthighbin);
 793   TH1D *thermNumSq_ss=(TH1D*)thermNumSq2d_ss->ProjectionY("thermNumSq_ss",Ktlowbin,Kthighbin);
 794   TH1D *thermDen_ss=(TH1D*)thermDen2d_ss->ProjectionY("thermDen_ss",Ktlowbin,Kthighbin);
 795   TH1D *thermLargeR_ss=(TH1D*)thermLargeR2D_ss->ProjectionY("thermLargeR_ss",Ktlowbin,Kthighbin);
 796   //
 797   TH2D *thermNum2d_os=(TH2D*)Thermfile->Get("Num_CosFSI_os");
 798   TH2D *thermNumSq2d_os=(TH2D*)Thermfile->Get("NumSq_CosFSI_os");
 799   TH2D *thermDen2d_os=(TH2D*)Thermfile->Get("Den_os");
 800   TH2D *thermLargeR2D_os=(TH2D*)Thermfile->Get("LargeRpairs_os");
 801   TH1D *thermNum_os=(TH1D*)thermNum2d_os->ProjectionY("thermNum_os",Ktlowbin,Kthighbin);
 802   TH1D *thermNumSq_os=(TH1D*)thermNumSq2d_os->ProjectionY("thermNumSq_os",Ktlowbin,Kthighbin);
 803   TH1D *thermDen_os=(TH1D*)thermDen2d_os->ProjectionY("thermDen_os",Ktlowbin,Kthighbin);
 804   TH1D *thermLargeR_os=(TH1D*)thermLargeR2D_os->ProjectionY("thermLargeR_os",Ktlowbin,Kthighbin);
 805   TH1D *C2Therm_ss = (TH1D*)thermNum_ss->Clone();
 806   TH1D *C2Therm_os = (TH1D*)thermNum_os->Clone();
 807   TH1D *C2Den_ss = (TH1D*)thermDen_ss->Clone();
 808   TH1D *C2Den_os = (TH1D*)thermDen_os->Clone();
 809   C2Therm_ss->Add(thermLargeR_ss);
 810   C2Den_ss->Add(thermLargeR_ss);
 811   C2Therm_os->Add(thermLargeR_os);
 812   C2Den_os->Add(thermLargeR_os);
 813   C2Therm_ss->Sumw2();
 814   C2Therm_os->Sumw2();
 815
 816   C2Therm_ss->Divide(C2Den_ss);
 817   C2Therm_os->Divide(C2Den_os);
 818
 819
 820   for(int i=1; i<=thermDen_ss->GetNbinsX(); i++){
 821     if(thermDen_ss->GetBinContent(i) > 0){
 822       double err = thermNumSq_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i);
 823       err -= pow(thermNum_ss->GetBinContent(i)/C2Den_ss->GetBinContent(i),2);
 824       err /= C2Den_ss->GetBinContent(i);
 825       err = err + pow(sqrt(thermLargeR_ss->GetBinContent(i))/C2Den_ss->GetBinContent(i),2);
 826       err += pow(sqrt(C2Den_ss->GetBinContent(i))*thermLargeR_ss->GetBinContent(i)/pow(C2Den_ss->GetBinContent(i),2),2);
 827       err = sqrt(err);
 828       C2Therm_ss->SetBinError(i, err);
 829     }
 830     if(thermDen_os->GetBinContent(i) > 0){
 831       double err = thermNumSq_os->GetBinContent(i)/C2Den_os->GetBinContent(i);
 832       err -= pow(thermNum_os->GetBinContent(i)/C2Den_os->GetBinContent(i),2);
 833       err /= C2Den_os->GetBinContent(i);
 834       err = err + pow(sqrt(thermLargeR_os->GetBinContent(i))/C2Den_os->GetBinContent(i),2);
 835       err += pow(sqrt(C2Den_os->GetBinContent(i))*thermLargeR_os->GetBinContent(i)/pow(C2Den_os->GetBinContent(i),2),2);
 836       err = sqrt(err);
 837       C2Therm_os->SetBinError(i, err);
 838     }
 839   }
 840
 841   C2Therm_ss->SetMarkerStyle(20);
 842   C2Therm_os->SetMarkerStyle(25);
 843   C2Therm_ss->SetMarkerColor(2);
 844   C2Therm_os->SetMarkerColor(4);
 845   C2Therm_ss->SetLineColor(2);
 846   C2Therm_os->SetLineColor(4);
 847   C2Therm_ss->SetMarkerSize(1.5);
 848   C2Therm_os->SetMarkerSize(1.5);
 849   C2Therm_ss->SetDirectory(0);
 850   C2Therm_os->SetDirectory(0);
 851   C2Therm_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 852   C2Therm_ss->GetYaxis()->SetTitle("1+<cos(qx)>");
 853   C2Therm_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 854   C2Therm_os->GetYaxis()->SetTitle("C_{2}^{+-}");
 855   C2Therm_os->SetMinimum(0.98);
 856   C2Therm_os->SetMaximum(1.25);
 857   //C2Therm_ss->DrawCopy();
 858   //C2Therm_os->DrawCopy("same");
 859   C2Therm_ss->SetDirectory(0);
 860   C2Therm_os->SetDirectory(0);
 861   //C2Therm->Fit(GaussEW,"IME","",0.005,0.1);
 862   Thermfile->Close();
 863
 864
 865
 866   /////////////////////////////////////////////////////
 867   // Global fitting C2os and C2ss
 868   double ThermEW[4]={0};
 869   double C2ssSys_e[BINRANGE_C2global]={0};
 870   double C2osSys_e[BINRANGE_C2global]={0};
 871   TF1 *fitC2ssTherm;
 872   TF1 *fitC2osTherm;
 873   TF1 *fitC2ssThermGaus;
 874   TF1 *fitC2osThermGaus;
 875   //
 876   const int npar=11;// 10
 877   TMinuit MyMinuit(npar);
 878   MyMinuit.SetFCN(fcnC2_Global);
 879   double OutputPar[npar]={0};
 880   double OutputPar_e[npar]={0};
 881
 882   double par[npar];               // the start values
 883   double stepSize[npar];          // step sizes
 884   double minVal[npar];            // minimum bound on parameter
 885   double maxVal[npar];            // maximum bound on parameter
 886   string parName[npar];
 887
 888   TH1D *temp_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 889   temp_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 890   TH1D *temp_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 891   temp_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 892   TH1D *temp_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 893   temp_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 894   TH1D *C2ssRaw=(TH1D*)temp_mm->Clone();// MRC uncorrected
 895   TH1D *C2osRaw=(TH1D*)temp_mp->Clone();// MRC uncorrected
 896   C2ssRaw->SetMarkerStyle(24);
 897   C2osRaw->SetMarkerStyle(21);//21
 898   C2ssRaw->SetMarkerColor(2);
 899   C2osRaw->SetMarkerColor(4);
 900
 901
 902   int iteration=0;
 903   while(iteration<3){// 0=pure Gaussian Therminator fits, 1=EW Therminator fits, 2=real data fits with Therm EW
 904     if(!IncludeEWfromTherm) iteration=2;// skip Therminator
 905
 906     for(int i=0; i<BINRANGE_C2global; i++){
 907
 908       C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.;
 909       if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= MomRes_C2_pp[i];
 910       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);
 911       C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.);
 912       C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2));
 913       C2ssFitting_e[i] += pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
 914       C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]);
 915       C2osFitting[i] = temp_mp->GetBinContent(i+1);
 916       if(!GeneratedSignal && !MCcase) C2osFitting[i] *= MomRes_C2_mp[i];
 917       C2osFitting_e[i] = pow(MomRes_C2_mp[i]*temp_mp->GetBinError(i+1),2);
 918       C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1));
 919       C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1));
 920       C2osFitting_e[i] += pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 921       C2osFitting_e[i] = sqrt(C2osFitting_e[i]);
 922       //
 923       C2ssSys_e[i] = pow((MomRes_C2_pp[i]-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
 924       C2ssSys_e[i] = sqrt(C2ssSys_e[i]);
 925       C2osSys_e[i] = pow((MomRes_C2_mp[i]-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 926       C2osSys_e[i] = sqrt(C2osSys_e[i]);
 927       //
 928
 929       K2SS[i] = CoulCorr2(+1, AvgQinvSS[i]);
 930       K2OS[i] = CoulCorr2(-1, AvgQinvOS[i]);
 931       //K2SS[i] = 1;
 932       //K2OS[i] = 1;
 933       //
 934       K2SS_e[i] = sqrt(pow((K2SS[i]-1)*0.02,2) + pow((K2SS[i]-Gamov(+1, AvgQinvSS[i]))*0.02,2));
 935       K2OS_e[i] = sqrt(pow((K2OS[i]-1)*0.02,2) + pow((K2OS[i]-Gamov(-1, AvgQinvOS[i]))*0.02,2));
 936       //K2SS_e[i] = 0.001;
 937       //K2OS_e[i] = 0.001;
 938
 939       //
 940       if(iteration<2){
 941         C2ssFitting[i] = C2Therm_ss->GetBinContent(i+1);// Therminator fit
 942         C2ssFitting_e[i] = C2Therm_ss->GetBinError(i+1);// Therminator fit
 943         C2osFitting[i] = C2Therm_os->GetBinContent(i+1);// Therminator fit
 944         C2osFitting_e[i] = C2Therm_os->GetBinError(i+1);// Therminator fit
 945         K2SS_e[i] = 0.;
 946         K2OS_e[i] = 0.;
 947       }
 948
 949     }
 950
 951
 952
 953     C2ssFitting[0]=-1000; C2osFitting[0]=-1000;
 954     C2ssFitting_e[0]=1000; C2osFitting_e[0]=1000;
 955     K2SS_e[0]=1000; K2OS_e[0]=1000;
 956
 957
 958
 959     par[0] = 1.0; par[1]=0.5; par[2]=0.5; par[3]=9.2; par[4] = .1; par[5] = .2; par[6] = .2; par[7] = 0.; par[8] = 0.; par[9] = 0.; par[10] = 0.01;
 960     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;
 961     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;
 962     maxVal[0] = 1.1; 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.1; maxVal[10]=1.;
 963     parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh";
 964     parName[5] = "kappa_3"; parName[6] = "kappa_4"; parName[7] = "kappa_5"; parName[8]="kappa_6"; parName[9]="Norm_2"; parName[10]="P_{coh}";
 965
 966     MyMinuit.DefineParameter(10, parName[10].c_str(), 0., stepSize[10], minVal[10], maxVal[10]); MyMinuit.FixParameter(10);// extra parameter
 967
 968     for (int i=0; i<npar; i++){
 969       MyMinuit.DefineParameter(i, parName[i].c_str(), par[i], stepSize[i], minVal[i], maxVal[i]);
 970     }
 971     //MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
 972     //MyMinuit.DefineParameter(1, parName[1].c_str(), .68, stepSize[1], minVal[1], maxVal[1]); MyMinuit.FixParameter(1);// lambda
 973     //MyMinuit.DefineParameter(2, parName[2].c_str(), 0.9999, stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
 974
 975     if(IncludeG==kFALSE || iteration<2) {
 976       MyMinuit.DefineParameter(2, parName[2].c_str(), 0., stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
 977       MyMinuit.DefineParameter(4, parName[4].c_str(), 0., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
 978     }else{
 979       Int_t err=0;
 980       Double_t tmp[1];
 981       tmp[0] = 2+1;
 982       MyMinuit.mnexcm( "RELease", tmp,  1,  err );// G
 983       tmp[0] = 4+1;
 984       MyMinuit.mnexcm( "RELease", tmp,  1,  err );// Rcoh
 985     }
 986
 987     //MyMinuit.DefineParameter(3, parName[3].c_str(), 10.88, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
 988     //MyMinuit.DefineParameter(4, parName[4].c_str(), 5., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
 989     MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
 990     MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
 991     //
 992     if(!IncludeEW){
 993       if(!IncludeEWfromTherm){
 994         // kappa3=0.24, kappa4=0.16 for testing
 995         MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
 996         MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
 997         MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
 998         MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
 999       }else{
1000         if(iteration==0){
1001           MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1002           MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1003         }else if(iteration==1){
1004           Int_t err=0;
1005           Double_t tmp[1];
1006           tmp[0] = 5+1;
1007           MyMinuit.mnexcm( "RELease", tmp,  1,  err );// k3
1008           tmp[0] = 6+1;
1009           MyMinuit.mnexcm( "RELease", tmp,  1,  err );// k4
1010         }else{// iteration=2
1011           MyMinuit.DefineParameter(5, parName[5].c_str(), ThermEW[0], stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1012           MyMinuit.DefineParameter(6, parName[6].c_str(), ThermEW[1], stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1013           MyMinuit.DefineParameter(7, parName[7].c_str(), ThermEW[2], stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
1014           MyMinuit.DefineParameter(8, parName[8].c_str(), ThermEW[3], stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
1015         }
1016       }// IncludeEWfromTherm
1017     }// IncludeEW
1018
1019     if(IncludeSS==kFALSE){
1020       MyMinuit.DefineParameter(3, parName[3].c_str(), 9.1, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
1021       MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
1022       MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
1023       MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
1024     }
1025     if(IncludeOS==kFALSE){
1026       MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2
1027     }
1028
1029
1030     int ierflg=0;
1031     double arglist[10];
1032     //arglist[0]=2;// improve Minimization Strategy (1 is default)
1033     //MyMinuit.mnexcm("SET STR",arglist,1,ierflg);
1034     //arglist[0] = 0;
1035     //MyMinuit.mnexcm("SCAN", arglist,1,ierflg);
1036     arglist[0] = 5000;
1037     MyMinuit.mnexcm("MIGRAD", arglist ,1,ierflg);
1038     // Do the minimization!
1039     cout<<"Start C2 Global fit"<<endl;
1040     MyMinuit.Migrad();// generally the best minimization algorithm
1041     cout<<"End C2 Global fit"<<endl;
1042
1043     for (int i=0; i<npar; i++){
1044       MyMinuit.GetParameter(i,OutputPar[i],OutputPar_e[i]);
1045     }
1046
1047     cout<<"Chi2/NDF = "<<Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars())<<"   Chi^2="<<Chi2_C2global<<endl;
1048     if(iteration==2) {
1049       ChiSquaredNDF->Fill(1, Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars()));
1050     }if(iteration==1) {
1051       ThermEW[0]=OutputPar[5]; ThermEW[1]=OutputPar[6]; ThermEW[2]=OutputPar[7]; ThermEW[3]=OutputPar[8];
1052       fitC2ssTherm = new TF1("fitC2ssTherm",C2ssFitFunction, 0.005,0.2, npar);
1053       for(int i=0; i<npar; i++) {
1054         fitC2ssTherm->FixParameter(i,OutputPar[i]);
1055         fitC2ssTherm->SetParError(i,OutputPar_e[i]);
1056       }
1057       fitC2osTherm = new TF1("fitC2osTherm",C2osFitFunction, 0.005,0.2, npar);
1058       for(int i=0; i<npar; i++) {
1059         fitC2osTherm->FixParameter(i,OutputPar[i]);
1060         fitC2osTherm->SetParError(i,OutputPar_e[i]);
1061       }
1062       fitC2osTherm->SetLineColor(4);
1063     }
1064     if(iteration==0){
1065       fitC2ssThermGaus = new TF1("fitC2ssThermGaus",C2ssFitFunction, 0.005,0.2, npar);
1066       for(int i=0; i<npar; i++) {
1067         fitC2ssThermGaus->FixParameter(i,OutputPar[i]);
1068         fitC2ssThermGaus->SetParError(i,OutputPar_e[i]);
1069       }
1070       fitC2osThermGaus = new TF1("fitC2osThermGaus",C2osFitFunction, 0.005,0.2, npar);
1071       for(int i=0; i<npar; i++) {
1072         fitC2osThermGaus->FixParameter(i,OutputPar[i]);
1073         fitC2osThermGaus->SetParError(i,OutputPar_e[i]);
1074       }
1075       fitC2osThermGaus->SetLineColor(4);
1076     }
1077
1078     iteration++;
1079
1080   }// iteration loop
1081
1082
1083   TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone();
1084   TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone();
1085   TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone();
1086   TH1D *C2_os_Momsys=(TH1D*)Two_ex_clone_mp->Clone();
1087   TH1D *C2_ss_Ksys=(TH1D*)Two_ex_clone_mm->Clone();
1088   TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone();
1089   TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone();
1090   TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone();
1091
1092   for(int i=0; i<BINRANGE_C2global; i++){
1093     C2_ss->SetBinContent(i+1, C2ssFitting[i]);
1094     C2_os->SetBinContent(i+1, C2osFitting[i]);
1095     C2_ss_Momsys->SetBinContent(i+1, C2ssFitting[i]);
1096     C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]);
1097     C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]);
1098     C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]);
1099     double staterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
1100     double staterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
1101     C2_ss->SetBinError(i+1, staterror_ss);
1102     C2_os->SetBinError(i+1, staterror_os);
1103     C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]);
1104     C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]);
1105     C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]);
1106     C2_os_Ksys->SetBinError(i+1, OutputPar[1]*K2OS_e[i]);
1107     //
1108     K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0);
1109     K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0);
1110   }
1111
1112   C2_ss_Momsys->SetMarkerSize(0);
1113   C2_ss_Momsys->SetFillStyle(1000);
1114   C2_ss_Momsys->SetFillColor(kRed-10);
1115   C2_os_Momsys->SetMarkerSize(0);
1116   C2_os_Momsys->SetFillStyle(1000);
1117   C2_os_Momsys->SetFillColor(kBlue-10);
1118   C2_ss_Ksys->SetMarkerSize(0);
1119   C2_ss_Ksys->SetFillStyle(3004);
1120   C2_ss_Ksys->SetFillColor(kRed);
1121   C2_os_Ksys->SetMarkerSize(0);
1122   C2_os_Ksys->SetFillStyle(3004);
1123   C2_os_Ksys->SetFillColor(kBlue);
1124
1125
1126
1127   C2_ss->GetXaxis()->SetRangeUser(0,0.098);
1128   C2_ss->GetYaxis()->SetRangeUser(0.98,1.3);
1129   C2_ss->GetXaxis()->SetTitleOffset(.8);
1130   C2_ss->GetYaxis()->SetTitleOffset(.8);
1131   C2_ss->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1132   C2_ss->GetXaxis()->SetTitleSize(0.05);
1133   C2_ss->GetYaxis()->SetTitleSize(0.05);
1134   C2_os->GetXaxis()->SetRangeUser(0,0.098);
1135   C2_os->GetYaxis()->SetRangeUser(0.98,1.3);
1136   C2_os->GetXaxis()->SetTitleOffset(.8);
1137   C2_os->GetYaxis()->SetTitleOffset(.8);
1138   C2_os->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
1139   C2_os->GetXaxis()->SetTitleSize(0.05);
1140   C2_os->GetYaxis()->SetTitleSize(0.05);
1141
1142   C2_ss->SetMarkerSize(1.5);
1143   C2_os->SetMarkerSize(1.5);
1144   C2_os->SetMarkerStyle(25);
1145   C2_os->SetMarkerColor(4);
1146
1147   TF1 *fitC2ss = new TF1("fitC2ss",C2ssFitFunction, 0.005,0.2, npar);
1148   TF1 *fitC2os = new TF1("fitC2os",C2osFitFunction, 0.005,0.2, npar);
1149   for(int i=0; i<npar; i++) {
1150     fitC2ss->FixParameter(i,OutputPar[i]);
1151     fitC2ss->SetParError(i,OutputPar_e[i]);
1152   }
1153   for(int i=0; i<npar; i++) {
1154     fitC2os->FixParameter(i,OutputPar[i]);
1155     fitC2os->SetParError(i,OutputPar_e[i]);
1156   }
1157
1158   if(!MCcase){
1159     C2_ss->DrawCopy();
1160     legend1->AddEntry(C2_ss,"same-charge","p");
1161     C2_os->DrawCopy("same");
1162     legend1->AddEntry(C2_os,"opp-charge","p");
1163     C2_ss_Momsys->DrawCopy("E2 same");
1164     C2_os_Momsys->DrawCopy("E2 same");
1165     C2_ss_Ksys->DrawCopy("E2 same");
1166     C2_os_Ksys->DrawCopy("E2 same");
1167     C2_ss->DrawCopy("same");
1168     C2_os->DrawCopy("same");
1169     fitC2ss->DrawCopy("same");
1170     fitC2os->SetLineColor(4);
1171     fitC2os->DrawCopy("same");
1172     MJ_bkg_ss->SetLineColor(1);
1173     MJ_bkg_ss->Draw("same");
1174     legend1->AddEntry(MJ_bkg_ss,"Non-femto bkg","l");
1175     legend1->Draw("same");
1176   }
1177
1178   /*
1179   ////////// C2 systematics
1180   // C2 -- base, M0, Pos B (0.035 TTC for all)
1181   //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};
1182   //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};
1183   // --, M0,
1184   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};
1185   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};
1186   // -+, M0,
1187   //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};
1188   //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};
1189   // --, M0, kt4
1190   //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};
1191   //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};
1192   // --, M0, kt5
1193   //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};
1194   //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};
1195   // --, M0, kt6
1196   //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};
1197   //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};
1198
1199   cout<<"C2 values:"<<endl;
1200   for(int ii=0; ii<20; ii++){
1201     cout<<Two_ex_clone_mm->GetBinContent(ii+1)<<", ";
1202   }
1203   cout<<endl;
1204   cout<<"C2 errors:"<<endl;
1205   for(int ii=0; ii<20; ii++){
1206     cout<<Two_ex_clone_mm->GetBinError(ii+1)<<", ";
1207   }
1208   cout<<endl;
1209   TH1D *C2ss_Sys=(TH1D*)Two_ex_clone_mm->Clone();
1210   for(int ii=1; ii<20; ii++){
1211     if(C2ss_base[ii] ==0) {
1212       C2ss_Sys->SetBinContent(ii+1, 100.);
1213       C2ss_Sys->SetBinError(ii+1, 100.);
1214       continue;
1215     }
1216     C2ss_Sys->SetBinContent(ii+1, (C2ss_base[ii]-C2ss_Sys->GetBinContent(ii+1))/C2ss_base[ii]);
1217     C2ss_Sys->SetBinError(ii+1, sqrt(fabs(pow(C2ss_Sys->GetBinError(ii+1),2) - pow(C2ss_base_e[ii],2))));
1218   }
1219   C2ss_Sys->GetXaxis()->SetRangeUser(0,0.099);
1220   C2ss_Sys->GetYaxis()->SetTitle("(C_{2}^{s}-C_{2}^{v})/C_{2}^{s}");
1221   C2ss_Sys->GetYaxis()->SetTitleOffset(1.3);
1222   C2ss_Sys->SetMinimum(-0.01);
1223   C2ss_Sys->SetMaximum(0.01);
1224   C2ss_Sys->Draw();
1225   TF1 *C2lineSys=new TF1("C2lineSys","pol0",0,0.099);
1226   //C2ss_Sys->Fit(C2lineSys,"MEQ","",0,0.099);
1227   */
1228
1229   // Momentum resolution uncorrected C2
1230   /*C2ssRaw->Draw("same");
1231   C2osRaw->Draw("same");
1232   legend1->AddEntry(C2ssRaw,"same-charge, Raw","p");
1233   legend1->AddEntry(C2osRaw,"opp-charge, Raw","p");
1234   legend1->Draw("same");
1235   */
1236
1237   // FSI corrected C2+-
1238   /*TH1D *C2_os_FSIcorrected=(TH1D*)C2_os->Clone();
1239   for(int ii=2; ii<20; ii++){
1240     double value = (C2_os_FSIcorrected->GetBinContent(ii) - (1-OutputPar[1]))/(OutputPar[1]*K2OS[ii-1]);
1241     C2_os_FSIcorrected->SetBinContent(ii,value);
1242   }
1243   C2_os_FSIcorrected->SetMinimum(0.95);
1244   C2_os_FSIcorrected->SetMarkerStyle(24);
1245   C2_os_FSIcorrected->Draw();
1246   */
1247
1248
1249   /*Two_ex_clone_mm->SetMarkerStyle(25);
1250   Two_ex_clone_mp->SetMarkerStyle(25);
1251   Two_ex_clone_mm->SetMarkerColor(2);
1252   Two_ex_clone_mp->SetMarkerColor(4);
1253   Two_ex_clone_mm->Draw("same");
1254   Two_ex_clone_mp->Draw("same");
1255   legend1->AddEntry(C2_ss,"same-charge; Momentum Resolution Corrected","p");
1256   legend1->AddEntry(C2_os,"opp-charge; Momentum Resolution Corrected","p");
1257   legend1->AddEntry(Two_ex_clone_mm,"same-charge; Raw","p");
1258   legend1->AddEntry(Two_ex_clone_mp,"opp-charge; Raw","p");
1259   legend1->Draw("same");
1260   */
1261   //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1262   //TGraph *gr2 = (TGraph*)MyMinuit.Contour(15,1,2);
1263   //gr2->GetXaxis()->SetTitle("lambda");
1264   //gr2->GetYaxis()->SetTitle("G");
1265   //gr2->SetTitle("");
1266   //gr2->SetFillColor(kYellow);
1267   //gr2->Draw("alf");
1268   //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1269   //MyMinuit.SetErrorDef(1);
1270   //TGraph *gr1 = (TGraph*)MyMinuit.Contour(15,1,2);
1271   //gr1->SetFillColor(kGreen);//38
1272   //gr1->Draw("lf");
1273
1274
1275
1276   ///////////////////////////
1277   // STAR comparison
1278   /* 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};
1279   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};
1280   TH1D *Two_star_mm=(TH1D*)Two_ex_clone_mm->Clone();
1281   Two_star_mm->Add(Two_star_mm,-1);
1282   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]);}
1283   Two_star_mm->SetMarkerColor(4);
1284   Two_star_mm->SetLineColor(4);
1285   //Two_star_mm->Draw("same");
1286   //legend1->AddEntry(Two_star_mm,"--, 200 GeV","p");
1287
1288   //Two_ex_clone_mp->Multiply(Two_ex_clone_pp);
1289   //Two_ex_clone_mp->Draw();
1290   //Two_ex_clone_mm->Multiply(Two_ex_clone_mp);
1291   //Two_ex_clone_mm->Draw();
1292   //Two_ex_clone_pp->Draw("same");
1293
1294   //legend1->Draw("same");
1295   */
1296
1297
1298
1299
1300
1301   /////////////////////////////////////////////////////////////////////////
1302   // 3-d fitting (out,side,long).  Not used for paper.
1303   /*
1304   TString *name1 = new TString("Explicit2_Charge1_0_Charge2_0_SC_0_M_");
1305   *name1 += Mbin;
1306   TString *name2 = new TString(name1->Data());
1307   TString *name3 = new TString(name1->Data());
1308   name1->Append("_ED_0_Term_1_osl_b2");// b1 (0.2<kt<0.3). b2 (0.6<kt<0.7)
1309   name2->Append("_ED_0_Term_2_osl_b2");
1310   name3->Append("_ED_0_Term_1_osl_b2_QW");
1311
1312
1313   TH3D *num_osl = (TH3D*)MyList->FindObject(name1->Data());
1314   TH3D *den_osl = (TH3D*)MyList->FindObject(name2->Data());
1315   den_osl->Scale(num_osl->Integral(28,40, 28,40, 28,40)/den_osl->Integral(28,40, 28,40, 28,40));
1316   TH3D *num_osl_QW = (TH3D*)MyList->FindObject(name3->Data());
1317
1318   for(int i=0; i<BINS_OSL; i++){
1319     for(int j=0; j<BINS_OSL; j++){
1320       for(int k=0; k<BINS_OSL; k++){
1321         if(num_osl->GetBinContent(i+1,j+1,k+1) < 1 || den_osl->GetBinContent(i+1,j+1,k+1) < 1) continue;
1322
1323         avg_q[i][j][k] = num_osl_QW->GetBinContent(i+1,j+1,k+1)/num_osl->GetBinContent(i+1,j+1,k+1);
1324         int QinvBin = int((avg_q[i][j][k])/0.005);
1325         if(QinvBin >=20) QinvBin=19;
1326         if(MomRes_term1_pp[QinvBin] ==0) continue;
1327         if(MomRes_term2_pp[QinvBin] ==0) continue;
1328         //
1329         num_osl->SetBinContent(i+1,j+1,k+1, num_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term1_pp[QinvBin]);
1330         den_osl->SetBinContent(i+1,j+1,k+1, den_osl->GetBinContent(i+1,j+1,k+1)*MomRes_term2_pp[QinvBin]);
1331         //
1332         A[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1333         B[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1334         //
1335         A_e[i][j][k] = num_osl->GetBinContent(i+1,j+1,k+1);
1336         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);
1337         A_e[i][j][k] = sqrt(A_e[i][j][k]);
1338         B_e[i][j][k] = den_osl->GetBinContent(i+1,j+1,k+1);
1339         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);
1340         B_e[i][j][k] = sqrt(B_e[i][j][k]);
1341         //
1342         K_OSL[i][j][k] = CoulCorr2(+1, avg_q[i][j][k]);
1343         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));
1344         //K_OSL_e[i][j][k] = 0;
1345       }
1346     }
1347   }
1348
1349
1350   const int npar_osl=6;//5
1351   TMinuit MyMinuit_osl(npar_osl);
1352   MyMinuit_osl.SetFCN(fcnOSL);
1353   double OutputPar_osl[npar_osl]={0};
1354   double OutputPar_osl_e[npar_osl]={0};
1355
1356   double par_osl[npar_osl];               // the start values
1357   double stepSize_osl[npar_osl];          // step sizes
1358   double minVal_osl[npar_osl];            // minimum bound on parameter
1359   double maxVal_osl[npar_osl];            // maximum bound on parameter
1360   string parName_osl[npar_osl];
1361
1362   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.;
1363   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;
1364   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.;
1365   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.;
1366   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";
1367
1368   for (int i=0; i<npar_osl; i++){
1369     MyMinuit_osl.DefineParameter(i, parName_osl[i].c_str(), par_osl[i], stepSize_osl[i], minVal_osl[i], maxVal_osl[i]);
1370   }
1371   MyMinuit_osl.FixParameter(2);// G
1372   //MyMinuit.FixParameter(1);// lambda
1373   cout<<"here!!"<<endl;
1374
1375   // Do the minimization!
1376   cout<<"Start Three-d fit"<<endl;
1377   MyMinuit_osl.Migrad();// Minuit's best minimization algorithm
1378   cout<<"End Three-d fit"<<endl;
1379   cout<<"Chi2/NDF = "<<Chi2_OSL/(NFitPoints_OSL-MyMinuit_osl.GetNumFreePars())<<endl;
1380
1381   for (int i=0; i<npar_osl; i++){
1382     MyMinuit_osl.GetParameter(i,OutputPar_osl[i],OutputPar_osl_e[i]);
1383   }
1384
1385
1386   TF3 *fit3d = new TF3("fit3d",OSLfitfunction, 0,0.2, 0,0.2, 0,0.2, npar_osl);
1387   for(int i=0; i<npar_osl; i++) {fit3d->FixParameter(i,OutputPar_osl[i]);}
1388   Int_t BinsOSL = num_osl->GetNbinsX();
1389   double LimitOSL = num_osl->GetXaxis()->GetBinUpEdge(BinsOSL);
1390   TH3D *den_timesFit = new TH3D("den_timesFit","",BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL, BinsOSL,0,LimitOSL);
1391   for(int i=0; i<BinsOSL; i++){
1392     for(int j=0; j<BinsOSL; j++){
1393       for(int k=0; k<BinsOSL; k++){
1394
1395         double qout=(i+.5)*0.005;
1396         double qside=(j+.5)*0.005;
1397         double qlong=(k+.5)*0.005;
1398         den_timesFit->SetBinContent(i+1,j+1,k+1, fit3d->Eval(qout,qside,qlong)*den_osl->GetBinContent(i+1,j+1,k+1));
1399         den_timesFit->SetBinError(i+1,j+1,k+1, 0);
1400       }
1401     }
1402   }
1403
1404   int binL=1, binH=4;
1405   TH1D *num_pro=(TH1D*)num_osl->ProjectionX("num_pro",binL,binH, binL,binH);
1406   TH1D *den_pro=(TH1D*)den_osl->ProjectionX("den_pro",binL,binH, binL,binH);
1407   TH1D *dentimesFit_pro=(TH1D*)den_timesFit->ProjectionX("dentimesFit_pro",binL,binH, binL,binH);
1408   num_pro->Sumw2();
1409   den_pro->Sumw2();
1410   num_pro->Divide(den_pro);
1411   num_pro->SetMarkerStyle(20);
1412   num_pro->SetTitle("");
1413   num_pro->GetXaxis()->SetTitle("q_{out}");
1414   num_pro->GetYaxis()->SetTitle("C_{2}");
1415   num_pro->SetMinimum(0.97);
1416   num_pro->SetMaximum(1.48);
1417   num_pro->DrawCopy();
1418   //
1419   dentimesFit_pro->Divide(den_pro);
1420   dentimesFit_pro->SetLineColor(2);
1421   dentimesFit_pro->DrawCopy("same");
1422   //
1423   */
1424   //MyMinuit.SetErrorDef(4); //note 4 and not 2!
1425   //TGraph *gr2 = (TGraph*)MyMinuit.Contour(10,1,2);
1426   //gr2->SetFillColor(kYellow);
1427   //gr2->Draw("alf");
1428   //Get contour for parameter 0 versus parameter 2 for ERRDEF=1
1429   //MyMinuit.SetErrorDef(1);
1430   //TGraph *gr1 = (TGraph*)MyMinuit.Contour(10,1,2);
1431   //gr1->SetFillColor(kGreen);//38
1432   //gr1->Draw("lf");
1433
1434
1435
1436   //////////////////////////////////////////////////////////////////////////////////////
1437
1438
1439   // To visualize the Qcut and Norm Regions
1440   //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
1441   //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);
1442   //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);
1443   //for(int bin=1; bin<=20; bin++) QcutRegion->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1444   //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1445   //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
1446   //QcutRegion->SetFillColor(4);
1447   //NormRegion1->SetFillColor(2);
1448   //NormRegion2->SetFillColor(3);
1449   //Two_ex[0][0][0][0]->Draw();
1450   //QcutRegion->Draw("same");
1451   //NormRegion1->Draw("same");
1452   //NormRegion2->Draw("same");
1453
1454
1455
1456
1457   TCanvas *can2 = new TCanvas("can2", "can2",800,0,900,900);//800,0,600,900
1458   can2->SetHighLightColor(2);
1459   gStyle->SetOptFit(0111);
1460   can2->SetFillColor(10);
1461   can2->SetBorderMode(0);
1462   can2->SetBorderSize(2);
1463   can2->SetGridx();
1464   can2->SetGridy();
1465   can2->SetFrameFillColor(0);
1466   can2->SetFrameBorderMode(0);
1467   can2->SetFrameBorderMode(0);
1468   can2->cd();
1469   gPad->SetLeftMargin(0.14);
1470   gPad->SetRightMargin(0.02);
1471
1472   TLegend *legend2 = new TLegend(.4,.67,1.,.87,NULL,"brNDC");
1473   legend2->SetBorderSize(1);
1474   legend2->SetTextSize(.06);// small .03; large .06
1475   legend2->SetFillColor(0);
1476
1477   /////////////////////////////////////////////
1478   TH3D *C3QS_3d = new TH3D("C3QS_3d","",20,0,.98, 20,0,.1, 20,0,.1);
1479   TH3D *Combinatorics_3d = new TH3D("Combinatorics_3d","",20,0,.1, 20,0,.1, 20,0,.1);
1480   //
1481   const int Q3BINS = 20;
1482   TH1D *num_withRS = new TH1D("num_withRS","",Q3BINS,0,0.2);
1483   TH1D *num_withGRS = new TH1D("num_withGRS","",Q3BINS,0,0.2);
1484   TH1D *num_withTM = new TH1D("num_withTM","",Q3BINS,0,0.2);
1485   TH1D *Cterm1 = new TH1D("Cterm1","",Q3BINS,0,0.2);
1486   TH1D *Cterm1_noMRC = new TH1D("Cterm1_noMRC","",Q3BINS,0,0.2);
1487   TH1D *Cterm2 = new TH1D("Cterm2","",Q3BINS,0,0.2);
1488   TH1D *Cterm3 = new TH1D("Cterm3","",Q3BINS,0,0.2);
1489   TH1D *Cterm4 = new TH1D("Cterm4","",Q3BINS,0,0.2);
1490   TH1D *num_QS = new TH1D("num_QS","",Q3BINS,0,0.2);
1491   TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,0.2);
1492   TH1D *Combinatorics_1d_noMRC = new TH1D("Combinatorics_1d_noMRC","",Q3BINS,0,0.2);
1493   TH1D *Coul_Riverside = new TH1D("Coul_Riverside","",Q3BINS,0,0.2);
1494   TH1D *Coul_GRiverside = new TH1D("Coul_GRiverside","",Q3BINS,0,0.2);
1495   TH1D *Coul_Omega0 = new TH1D("Coul_Omega0","",Q3BINS,0,0.2);
1496   TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,0.2);
1497   TH1D *c3_hist_STAR = new TH1D("c3_hist_STAR","",Q3BINS,0,0.2);
1498   TProfile *MomRes_2d = new TProfile("MomRes_2d","",Q3BINS,0,0.2, 0,20.,"");
1499   TProfile *MomRes_3d = new TProfile("MomRes_3d","",Q3BINS,0,0.2, 0,20.,"");
1500   TH1D *r3_num_Q3 = new TH1D("r3_num_Q3","",Q3BINS,0,0.2);
1501   TH1D *r3_den_Q3 = new TH1D("r3_den_Q3","",Q3BINS,0,0.2);
1502   r3_num_Q3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1503   r3_num_Q3->GetYaxis()->SetTitle("r_{3}");
1504   r3_num_Q3->GetYaxis()->SetTitleOffset(1.3);
1505   r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
1506   r3_num_Q3->SetMinimum(0);
1507   r3_num_Q3->SetMaximum(2.4);
1508   r3_den_Q3->SetLineColor(2);
1509   r3_num_Q3->SetMarkerStyle(20);
1510   Coul_Omega0->GetXaxis()->SetRangeUser(0,0.099);
1511   Coul_Omega0->GetXaxis()->SetLabelSize(0.04);
1512   Coul_Omega0->GetYaxis()->SetLabelSize(0.04);
1513   c3_hist_STAR->GetXaxis()->SetRangeUser(0,0.099);
1514   c3_hist_STAR->SetMinimum(0.8); c3_hist_STAR->SetMaximum(1.02);
1515   c3_hist_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1516   c3_hist_STAR->GetYaxis()->SetTitle("c_{3}*^{++-}");
1517   c3_hist_STAR->GetYaxis()->SetTitleOffset(1.6);
1518   if(SameCharge) {Cterm1_noMRC->Sumw2(); Combinatorics_1d_noMRC->Sumw2();}
1519   //
1520   double num_QS_e[Q3BINS]={0};
1521   double c3_e[Q3BINS]={0};
1522   double r3_e[Q3BINS]={0};
1523
1524   // CB=Charge Bin. 0 means pi-, 1 means pi+
1525   int CB1=0, CB2=0, CB3=0;
1526   int CP12=1, CP13=1, CP23=1;
1527   if(CHARGE==-1) {
1528     if(SameCharge) {CB1=0; CB2=0; CB3=0;}
1529     else {CB1=0; CB2=0; CB3=1; CP12=+1, CP13=-1, CP23=-1;}
1530   }else {
1531     if(SameCharge) {CB1=1; CB2=1; CB3=1;}
1532     else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;}
1533   }
1534
1535   // SC = species combination.  Always 0 meaning pi-pi-pi.
1536   int SCBin=0;
1537
1538   //
1539   ReadCoulCorrections(RValue, bValue, 10);// switch to full kt range, 10.
1540   //ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation (STAR method testing)
1541   TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",20,0,0.2);
1542   TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",20,0,0.2);
1543   //
1544   double value_num[2]={0};
1545   double value_num_e[2]={0};
1546   double value_den[2]={0};
1547   double N3_QS[2]={0};
1548   double N3_QS_e[2]={0};
1549   double OutTriplets=0, InTriplets=0;
1550   for(int i=2; i<=20; i++){// bin number
1551     //double Q12_m = (i-0.5)*(0.005);// geometric center
1552     double Q12 = AvgQinvSS[i-1]; if(!SameCharge && CHARGE==+1) {Q12 = AvgQinvOS[i-1];}// true center
1553     int Q12bin = int(Q12/0.005)+1;
1554     //
1555     for(int j=2; j<=20; j++){// bin number
1556       //double Q13_m = (j-0.5)*(0.005);// geometric center
1557       double Q13 = AvgQinvSS[j-1]; if(!SameCharge) {Q13 = AvgQinvOS[j-1];}// true center
1558       int Q13bin = int(Q13/0.005)+1;
1559       //
1560       for(int k=2; k<=20; k++){// bin number
1561         //double Q23_m = (k-0.5)*(0.005);// geometric center
1562         double Q23 = AvgQinvSS[k-1]; if(!SameCharge && CHARGE==-1) {Q23 = AvgQinvOS[k-1];}// true center
1563         int Q23bin = int(Q23/0.005)+1;
1564         //
1565         //if(fabs(i-j)>3 || fabs(i-k)>3 || fabs(j-k)>3) continue;// testing
1566
1567         double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
1568         int Q3bin = Cterm1->GetXaxis()->FindBin(Q3);
1569
1570         if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
1571         if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
1572
1573
1574         //
1575         double G_12 = Gamov(CP12, Q12);// point-source Coulomb correlation
1576         double G_13 = Gamov(CP13, Q13);
1577         double G_23 = Gamov(CP23, Q23);
1578         double K2_12 = CoulCorr2(CP12, Q12);// finite-source Coulomb+Strong correlation from Therminator.
1579         double K2_13 = CoulCorr2(CP13, Q13);
1580         double K2_23 = CoulCorr2(CP23, Q23);
1581         double K3 = 1.;// 3-body Coulomb+Strong correlation
1582         if(GRS) {// GRS approach
1583           K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
1584           if(CHARGE==+1 && !SameCharge) K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
1585         }else {
1586           K3 = CoulCorrOmega0(SameCharge, Q12, Q13, Q23);
1587           if(CHARGE==+1 && !SameCharge) K3 = CoulCorrOmega0(SameCharge, Q23, Q12, Q13);
1588         }
1589         if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;}
1590         if(K3==0) continue;
1591
1592         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
1593         double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
1594         double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
1595         double TERM4=Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event
1596         double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
1597         if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
1598         if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
1599
1600         if(Q3>0.08 && Q3<0.09){// just for testing
1601           if(Q12>0.08 || Q13>0.08 || Q23>0.08) OutTriplets++;
1602           else InTriplets++;
1603         }
1604
1605         // apply momentum resolution correction
1606         if(!MCcase && !GeneratedSignal){
1607           if(SameCharge){
1608             // 3d momentum resolution corrections
1609             TERM1 *= (MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1610             TERM2 *= (MomRes3d[0][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1611             TERM3 *= (MomRes3d[0][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1612             TERM4 *= (MomRes3d[0][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1613             TERM5 *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1614             // Triple produce of 1-d momentum resolution corrections (less accurate).
1615             /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1];
1616             TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];
1617             TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term2_pp[k-1];
1618             TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term1_pp[k-1];
1619             TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_pp[j-1]*MomRes_term2_pp[k-1];*/
1620             //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_pp[j-1]*MomRes_term1_pp[k-1]);
1621             //MomRes_3d->Fill(Q3, MomRes3d[0][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
1622           }else{
1623             if(CHARGE==-1){// pi-pi-pi+
1624               TERM1 *= (MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1625               TERM2 *= (MomRes3d[1][1]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1626               TERM3 *= (MomRes3d[1][2]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1627               TERM4 *= (MomRes3d[1][3]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1628               TERM5 *= (MomRes3d[1][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;
1629               /*TERM1 *= MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1];
1630               TERM2 *= MomRes_term1_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];
1631               TERM3 *= MomRes_term2_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_mp[k-1];
1632               TERM4 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_mp[k-1];
1633               TERM5 *= MomRes_term2_pp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_mp[k-1];*/
1634             }else {// pi-pi+pi+
1635               TERM1 *= (MomRes3d[1][0]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1636               TERM2 *= (MomRes3d[1][3]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1637               TERM3 *= (MomRes3d[1][2]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1638               TERM4 *= (MomRes3d[1][1]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1639               TERM5 *= (MomRes3d[1][4]->GetBinContent(Q23bin, Q13bin, Q12bin)-1)*MRCShift+1;
1640               /*TERM1 *= MomRes_term1_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_pp[k-1];
1641               TERM2 *= MomRes_term1_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];
1642               TERM3 *= MomRes_term2_mp[i-1]*MomRes_term1_mp[j-1]*MomRes_term2_pp[k-1];
1643               TERM4 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term1_pp[k-1];
1644               TERM5 *= MomRes_term2_mp[i-1]*MomRes_term2_mp[j-1]*MomRes_term2_pp[k-1];*/
1645             }
1646             //MomRes_2d->Fill(Q3, MomRes_term1_pp[i-1]*MomRes_term1_mp[j-1]*MomRes_term1_mp[k-1]);// always treats 12 as ss pair
1647             //MomRes_3d->Fill(Q3, MomRes3d[1][0]->GetBinContent(Q12bin, Q13bin, Q23bin));
1648           }
1649         }
1650         //
1651
1652         for(int LamType=0; LamType<2; LamType++){
1653
1654           if(LamType==1){TwoFrac=TwoFracr3; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);}// r3 assignment
1655           else {
1656             TwoFrac = OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);// Assign to C2 global fit values found previously
1657           }
1658
1659           // Purify.  Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
1660           N3_QS[LamType] = TERM1;
1661           N3_QS[LamType] -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1662           N3_QS[LamType] -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
1663           N3_QS[LamType] /= ThreeFrac;
1664           N3_QS[LamType] /= K3;
1665
1666           if(LamType==0) num_QS->Fill(Q3, N3_QS[LamType]);
1667
1668           // Isolate 3-pion cumulant
1669           value_num[LamType] = N3_QS[LamType];
1670           value_num[LamType] -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac;
1671           value_num[LamType] -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac;
1672           value_num[LamType] -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac;
1673           value_num[LamType] += 2*TERM5;
1674
1675
1676           // r3 denominator
1677           if(LamType==1 && SameCharge) {
1678             value_den[LamType] = PionNorm[CB1]->GetBinContent(i,j,k);// Raw r3 denominator
1679             if(!MCcase){
1680               // Momentum Resolution Correction Systematic variation. Only important when MRCShift != 1.0.
1681               double denMRC1 = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1682               denMRC1 *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1683               denMRC1 *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1684               double denMRC2 = (C2ssRaw->GetBinContent(i)*((MomRes_C2_pp[i-1]-1)*MRCShift+1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1685               denMRC2 *= (C2ssRaw->GetBinContent(j)*((MomRes_C2_pp[j-1]-1)*MRCShift+1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1686               denMRC2 *= (C2ssRaw->GetBinContent(k)*((MomRes_C2_pp[k-1]-1)*MRCShift+1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1687               // Non-femto background correction (Mini-Jet).
1688               double denMJ = (C2ssRaw->GetBinContent(i)*MomRes_C2_pp[i-1] - (MJ_bkg_ss->Eval(Q12)-1) - TwoFrac*K2_12 - (1-TwoFrac))/(TwoFrac*K2_12);
1689               denMJ *= (C2ssRaw->GetBinContent(j)*MomRes_C2_pp[j-1] - (MJ_bkg_ss->Eval(Q13)-1) - TwoFrac*K2_13 - (1-TwoFrac))/(TwoFrac*K2_13);
1690               denMJ *= (C2ssRaw->GetBinContent(k)*MomRes_C2_pp[k-1] - (MJ_bkg_ss->Eval(Q23)-1) - TwoFrac*K2_23 - (1-TwoFrac))/(TwoFrac*K2_23);
1691
1692               double den_ratio = sqrt(fabs(denMRC2))/sqrt(fabs(denMRC1));
1693               value_den[LamType] *= den_ratio;// apply Momentum Resolution correction systematic variation if any.
1694               double den_ratioMJ = sqrt(fabs(denMJ))/sqrt(fabs(denMRC1));
1695               if(IncludeMJcorrection) value_den[LamType] *= den_ratioMJ;// Non-femto bkg correction
1696
1697               //value_den[LamType] -= TPNRejects->GetBinContent(i,j,k);// Testing for 0-5% only to estimate the effect when C2^QS < 1.0
1698               value_den[LamType] *= (MomRes3d[0][4]->GetBinContent(Q12bin, Q13bin, Q23bin)-1)*MRCShift+1;// additional momentum resolution correction
1699             }
1700           }
1701
1702           // errors
1703           N3_QS_e[LamType] = TERM1;
1704           N3_QS_e[LamType] += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2);
1705           N3_QS_e[LamType] += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2));
1706           N3_QS_e[LamType] /= pow(ThreeFrac,2);
1707           N3_QS_e[LamType] /= pow(K3,2);
1708           if(LamType==0) num_QS_e[Q3bin-1] += N3_QS_e[LamType];
1709           // errors
1710           value_num_e[LamType] = N3_QS_e[LamType];
1711           value_num_e[LamType] += (pow(1/K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(TERM5),2));
1712           value_num_e[LamType] += (pow(1/K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(TERM5),2));
1713           value_num_e[LamType] += (pow(1/K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(TERM5),2));
1714           value_num_e[LamType] += pow(2*sqrt(TERM5),2);
1715           if(LamType==0) c3_e[Q3bin-1] += value_num_e[LamType] + TERM5;// add baseline stat error
1716           else r3_e[Q3bin-1] += value_num_e[LamType];
1717         }
1718
1719         // Fill histograms
1720         c3_hist->Fill(Q3, value_num[0] + TERM5);// for cumulant correlation function
1721         C3QS_3d->SetBinContent(i,j,k, N3_QS[0]);
1722         C3QS_3d->SetBinError(i,j,k, N3_QS_e[0]);
1723         //
1724         Coul_Riverside->Fill(Q3, G_12*G_13*G_23*TERM5);
1725         Coul_GRiverside->Fill(Q3, TERM5*CoulCorrGRS(SameCharge, Q12, Q13, Q23));
1726         Coul_Omega0->Fill(Q3, K3*TERM5);
1727         num_withGRS->Fill(Q3, N3_QS[0]);
1728         Cterm1->Fill(Q3, TERM1);
1729         Cterm2->Fill(Q3, TERM2);
1730         Cterm3->Fill(Q3, TERM3);
1731         Cterm4->Fill(Q3, TERM4);
1732         Combinatorics_1d->Fill(Q3, TERM5);
1733         Combinatorics_3d->Fill(Q12,Q13,Q23, TERM5);
1734         r3_num_Q3->Fill(Q3, value_num[1]);
1735         r3_den_Q3->Fill(Q3, value_den[1]);
1736         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));
1737         Cterm1_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k));
1738         Combinatorics_1d_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1739         double cumulant_STAR = Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k)/(K2_12*K2_13*K2_23);
1740         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k)/(K2_12);
1741         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k)/(K2_13);
1742         cumulant_STAR -= Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k)/(K2_23);
1743         c3_hist_STAR->Fill(Q3_m, cumulant_STAR + 3*Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1744
1745         ///////////////////////////////////////////////////////////
1746         // Edgeworth 3-pion expection.  Not important for r3.
1747         //double radius_exp = (11-(MomResCentBin-1))/FmToGeV;
1748         //TwoFrac=0.68; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1749         double radius_exp = (OutputPar[3])/FmToGeV;
1750         TwoFrac=OutputPar[1]; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1751         double arg12 = Q12*radius_exp;
1752         double arg13 = Q13*radius_exp;
1753         double arg23 = Q23*radius_exp;
1754         /*Float_t EW12 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1755         EW12 += 0.16/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1756         Float_t EW13 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1757         EW13 += 0.16/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1758         Float_t EW23 = 1 + 0.24/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1759         EW23 += 0.16/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);*/
1760
1761         Float_t EW12 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1762         EW12 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1763         Float_t EW13 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1764         EW13 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1765         Float_t EW23 = 1 + OutputPar[5]/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1766         EW23 += OutputPar[6]/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
1767         //
1768         double cumulant_exp=0, term1_exp=0;
1769         if(SameCharge) {
1770           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;
1771           cumulant_exp += 2*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5;
1772           term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1773           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;
1774           term1_exp += ThreeFrac*K3*cumulant_exp;
1775           //term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1776         }else {
1777           cumulant_exp = (exp(-pow(radius_exp*Q12,2))*pow(EW12,2))*TERM5 + TERM5;
1778           term1_exp = ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1779           term1_exp += TwoFrac*(1-OneFrac)*(K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)) + K2_13 + K2_23)*TERM5;
1780           term1_exp += ThreeFrac*K3*cumulant_exp;
1781           term1_exp = ((1-TwoFrac) + TwoFrac*K2_12*(1+exp(-pow(radius_exp*Q12,2))*pow(EW12,2)))*TERM5;
1782           //term1_exp = ((1-TwoFrac) + TwoFrac*K2_13)*TERM5;
1783         }
1784
1785         GenSignalExpected_num->Fill(Q3, term1_exp);
1786         GenSignalExpected_den->Fill(Q3, TERM5);
1787         ///////////////////////////////////////////////////////////
1788
1789       }
1790     }
1791   }
1792
1793   cout<<"OutTriplets: "<<OutTriplets<<"   InTriplets: "<<InTriplets<<endl;
1794   ////////////////////////////
1795
1796   // Intermediate steps
1797   num_withRS->Sumw2();
1798   num_withGRS->Sumw2();
1799   num_withTM->Sumw2();
1800   Cterm1->Sumw2();
1801   Cterm2->Sumw2();
1802   Cterm3->Sumw2();
1803   Cterm4->Sumw2();
1804   Combinatorics_1d->Sumw2();
1805   Combinatorics_3d->Sumw2();
1806   r3_num_Q3->Sumw2();
1807   r3_den_Q3->Sumw2();
1808   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]));}
1809   for(int i=0; i<Q3BINS; i++) {c3_hist_STAR->SetBinError(i+1, sqrt(c3_e[i]));}// approximate error
1810   num_withRS->Divide(Combinatorics_1d);
1811   num_withGRS->Divide(Combinatorics_1d);
1812   num_withTM->Divide(Combinatorics_1d);
1813   for(int q3bin=1; q3bin<=r3_num_Q3->GetNbinsX(); q3bin++){
1814     r3_num_Q3->SetBinError(q3bin, sqrt(r3_e[q3bin-1]));
1815   }
1816   Cterm1->Divide(Combinatorics_1d);
1817   Cterm1_noMRC->Divide(Combinatorics_1d_noMRC);
1818   Cterm2->Divide(Combinatorics_1d);
1819   Cterm3->Divide(Combinatorics_1d);
1820   Cterm4->Divide(Combinatorics_1d);
1821   c3_hist->Divide(Combinatorics_1d);
1822   c3_hist_STAR->Divide(Combinatorics_1d_noMRC);
1823   num_QS->Divide(Combinatorics_1d);
1824   r3_num_Q3->Divide(r3_den_Q3);
1825   GenSignalExpected_num->Sumw2();
1826   GenSignalExpected_num->Divide(GenSignalExpected_den);
1827
1828   //
1829   //
1830   Coul_Riverside->Divide(Combinatorics_1d);
1831   Coul_GRiverside->Divide(Combinatorics_1d);
1832   Coul_Omega0->Divide(Combinatorics_1d);
1833   for(int ii=1; ii<=Coul_Omega0->GetNbinsX(); ii++){
1834     Coul_Riverside->SetBinError(ii,0.000001);
1835     Coul_GRiverside->SetBinError(ii,0.000001);
1836     Coul_Omega0->SetBinError(ii,0.000001);
1837   }
1838   ////////////////////////////
1839
1840
1841
1842   ///////////////////////////////////////////////////////////////////////////////////////////////////
1843
1844   Coul_Riverside->SetLineColor(1);
1845   Coul_GRiverside->SetLineColor(2);
1846   Coul_Omega0->SetLineColor(4);
1847
1848   //Coul_Riverside->Draw();
1849   //legend1->AddEntry(Coul_Riverside,"Riverside","l");
1850   //Coul_GRiverside->Draw("same");
1851   //legend1->AddEntry(Coul_GRiverside,"Generalized Riverside","l");
1852   //Coul_Omega0->Draw("same");
1853   //legend1->AddEntry(Coul_Omega0,"Omega0","l");
1854
1855
1856   Cterm1->SetMarkerStyle(20);
1857   Cterm1->SetMarkerColor(4);
1858   Cterm1->SetLineColor(4);
1859   Cterm1->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1860   Cterm1->GetYaxis()->SetTitle("C_{3}");
1861   //Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}");
1862   Cterm1->SetMinimum(0.97);
1863   Cterm1->SetMaximum(1.7);// 6.1 for same-charge
1864   Cterm1->GetXaxis()->SetRangeUser(0,.095);
1865   Cterm1->GetYaxis()->SetTitleOffset(1.4);
1866   Cterm1->Draw();
1867   legend2->AddEntry(Cterm1,"C_{3} raw","p");
1868   //
1869   Cterm2->SetMarkerStyle(20);
1870   Cterm2->SetMarkerColor(7);
1871
1872   Cterm3->SetMarkerStyle(25);
1873   Cterm3->SetMarkerColor(8);
1874   Cterm4->SetMarkerStyle(26);
1875   Cterm4->SetMarkerColor(9);
1876   //Cterm2->Draw("same");
1877   //Cterm3->Draw("same");
1878   //Cterm4->Draw("same");
1879   //legend2->AddEntry(Cterm1,"++-","p");
1880
1881
1882   if(MCcase){
1883     double C3points_HIJING_mmm[10]={0, 0.961834, 1.01827, 0.999387, 1.00202, 1.00081, 1.00082, 1.00037, 0.999074, 0.999099};
1884     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};
1885     double C3points_HIJING_ppp[10]={0, 1.13015, 1.00623, 1.00536, 1.00293, 0.999964, 1.00116, 1.0007, 1.00037, 1.00105};
1886     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};
1887     TH1D *C3_HIJING_mmm=(TH1D*)Cterm1->Clone();
1888     TH1D *C3_HIJING_ppp=(TH1D*)Cterm1->Clone();
1889     for(int i=0; i<10; i++){
1890       C3_HIJING_mmm->SetBinContent(i+1, C3points_HIJING_mmm[i]);
1891       C3_HIJING_mmm->SetBinError(i+1, C3points_HIJING_mmm_e[i]);
1892       C3_HIJING_ppp->SetBinContent(i+1, C3points_HIJING_ppp[i]);
1893       C3_HIJING_ppp->SetBinError(i+1, C3points_HIJING_ppp_e[i]);
1894     }
1895     C3_HIJING_mmm->SetMarkerColor(2);
1896     C3_HIJING_mmm->SetLineColor(2);
1897     C3_HIJING_ppp->SetMarkerColor(4);
1898     C3_HIJING_ppp->SetLineColor(4);
1899     //C3_HIJING_mmm->Draw("same");
1900     //C3_HIJING_ppp->Draw("same");
1901     //legend2->AddEntry(C3_HIJING_mmm,"---","p");
1902     //legend2->AddEntry(C3_HIJING_ppp,"+++","p");
1903   }
1904
1905   num_QS->SetMarkerStyle(24);
1906   num_QS->SetMarkerColor(4);
1907   num_QS->SetLineColor(4);
1908   num_QS->GetXaxis()->SetTitle("Q_{3}");
1909   num_QS->GetYaxis()->SetTitle("C_{3}^{QS}");
1910   num_QS->GetXaxis()->SetRangeUser(0,.12);
1911   num_QS->SetMaximum(6);
1912   //num_QS->Draw("same");
1913   //legend2->AddEntry(num_QS,"C_{3}^{QS}","p");
1914
1915   c3_hist->GetXaxis()->SetTitle("Q_{3}");
1916   c3_hist->GetYaxis()->SetTitle("c_{3}^{++-}");
1917   c3_hist->GetXaxis()->SetRangeUser(0,.15);
1918   c3_hist->SetMarkerStyle(25);
1919   c3_hist->SetMarkerColor(2);
1920   c3_hist->SetLineColor(2);
1921   c3_hist->SetMaximum(3);
1922   c3_hist->SetMinimum(.9);
1923   if(!MCcase) c3_hist->Draw("same");
1924   //legend2->AddEntry(c3_hist,"#font[12]{c}_{3} (cumulant correlation)","p");
1925
1926
1927   GenSignalExpected_num->SetMarkerStyle(20);
1928   //GenSignalExpected_num->Draw("same");
1929   //legend2->AddEntry(GenSignalExpected_num,"#kappa_{3}=0.24, #kappa_{4}=0.16, #lambda=0.68, R=6 fm","p");
1930   //legend2->AddEntry(GenSignalExpected_num,"Edeworth expectation (fully chaotic)","p");
1931
1932   //MomRes_2d->SetMarkerStyle(20);
1933   //MomRes_3d->SetMarkerStyle(20);
1934   //MomRes_3d->SetMarkerColor(4);
1935   //MomRes_2d->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
1936   //MomRes_2d->GetYaxis()->SetTitle("< MRC >");
1937   //MomRes_2d->SetTitle("");
1938   //MomRes_2d->Draw();
1939   //legend2->AddEntry(MomRes_2d,"2D: Triple pair product","p");
1940   //MomRes_3d->Draw("same");
1941   //legend2->AddEntry(MomRes_3d,"3D","p");
1942
1943   //legend2->Draw("same");
1944
1945
1946   /*
1947   ////////// C3 systematics
1948   // C3 --- base, M0, (0.035 TTC for all)
1949   //double C3_base[10]={0, 1.63072, 1.6096, 1.43731, 1.28205, 1.17777, 1.11973, 1.07932, 1.05459, 1.04029};
1950   //double C3_base_e[10]={0, 0.022528, 0.00387504, 0.00115667, 0.000423799, 0.000238973, 0.000143993, 9.71502e-05, 7.02007e-05, 5.31267e-05};
1951   // C3 --- base, M0, Pos B (0.035 TTC for all)
1952   //double C3_base[10]={0, 1.62564, 1.60461, 1.438, 1.28234, 1.17827, 1.12009, 1.07953, 1.05474, 1.04037};
1953   //double C3_base_e[10]={0, 0.0239233, 0.00409002, 0.0012215, 0.000446701, 0.000251769, 0.000151651, 0.000102284, 7.38993e-05, 5.59212e-05};
1954   // C3 --+ base, M0, (0.035 TTC for all)
1955   //double C3_base[10]={0, 1.66016, 1.41961, 1.25229, 1.16459, 1.11254, 1.08012, 1.05768, 1.04265, 1.0332};
1956   //double C3_base_e[10]={0, 0.00539779, 0.00111398, 0.000387926, 0.000192906, 0.00011428, 7.24765e-05, 5.09126e-05, 3.76421e-05, 2.87533e-05};
1957
1958   // C3 --- base, M0, (New Standard: 0.04 TTC )
1959   //double C3_base[10]={0, 1.63903, 1.60254, 1.43381, 1.2794, 1.17603, 1.11806, 1.07806, 1.05345, 1.03936};
1960   //double C3_base_e[10]={0, 0.0322796, 0.00438361, 0.00122249, 0.000424557, 0.000233965, 0.000139058, 9.28136e-05, 6.66159e-05, 5.01816e-05};
1961   // C3 --- base, M1, (New Standard: 0.04 TTC )
1962   //double C3_base[10]={0, 1.6127, 1.65561, 1.49508, 1.33093, 1.21458, 1.14708, 1.099, 1.06864, 1.05064};
1963   //double C3_base_e[10]={0, 0.0425447, 0.0061176, 0.00172948, 0.000600699, 0.000329342, 0.000194832, 0.000129427, 9.25599e-05, 6.95395e-05};
1964   // C3 --- base, M2, (New Standard: 0.04 TTC )
1965   //double C3_base[10]={0, 1.6509, 1.71863, 1.54652, 1.38092, 1.25226, 1.17549, 1.12068, 1.08408, 1.06236};
1966   //double C3_base_e[10]={0, 0.0981473, 0.0143699, 0.00404612, 0.00141189, 0.000770764, 0.000453944, 0.000300452, 0.000214068, 0.000160448};
1967   // C3 --- base, M9, (New Standard: 0.04 TTC )
1968   //double C3_base[10]={0, 2.41982, 2.18303, 1.93205, 1.80399, 1.60955, 1.49305, 1.38565, 1.29873, 1.23626};
1969   //double C3_base_e[10]={0, 1.60227, 0.177274, 0.0501455, 0.018456, 0.00998147, 0.00583719, 0.00379465, 0.00264116, 0.0019391};
1970   //
1971   // C3 --+ base, M0, (New Standard: 0.04 TTC )
1972   //double C3_base[10]={0, 1.66087, 1.41943, 1.25081, 1.16313, 1.11143, 1.07917, 1.05701, 1.04215, 1.0328};
1973   //double C3_base_e[10]={0, 0.00584743, 0.00111278, 0.000374009, 0.00018315, 0.000107523, 6.78669e-05, 4.75116e-05, 3.50489e-05, 2.67328e-05};
1974   //
1975   // HIJING C3 --- base, M0
1976   //double C3_base[10]={0, 0.970005, 1.00655, 1.00352, 1.00291, 1.00071, 1.0002, 0.999524, 0.999404, 0.999397};
1977   //double C3_base_e[10]={0, 0.050845, 0.0099602, 0.00334862, 0.00138008, 0.000841743, 0.000531776, 0.000371712, 0.000274716, 0.00021};
1978   // HIJING C3 --- base, M1
1979   //double C3_base[10]={0, 1.03657, 1.00199, 0.997984, 1.00067, 1.0006, 0.999901, 0.999967, 0.999792, 0.999777};
1980   //double C3_base_e[10]={0, 0.0634232, 0.0117204, 0.0039446, 0.00163131, 0.000996638, 0.000629662, 0.000440266, 0.00032534, 0.000249};
1981   // HIJING C3 --- base, M2
1982   //double C3_base[10]={0, 1.34345, 1.04226, 1.0278, 0.99582, 1.00554, 1.00296, 1.00057, 1.00271, 1.00152};
1983   //double C3_base_e[10]={0, 0.363559, 0.0503531, 0.0170117, 0.00679185, 0.00419035, 0.00264603, 0.00184587, 0.00136663, 0.00104772};
1984   // HIJING C3 --- base, M3
1985   double C3_base[10]={0, 0.890897, 1.05222, 1.00461, 1.01364, 0.998981, 1.00225, 1.00305, 1.00235, 1.00043};
1986   double C3_base_e[10]={0, 0.297124, 0.0604397, 0.0195066, 0.00812906, 0.00490835, 0.00310751, 0.00217408, 0.00160575, 0.00123065};
1987   TH1D *C3baseHisto=(TH1D*)Cterm1->Clone();
1988   for(int ii=0; ii<10; ii++){
1989     C3baseHisto->SetBinContent(ii+1, C3_base[ii]);
1990     C3baseHisto->SetBinError(ii+1, C3_base_e[ii]);
1991   }
1992
1993   cout<<"C3 values:"<<endl;
1994   for(int ii=0; ii<10; ii++){
1995     cout<<Cterm1->GetBinContent(ii+1)<<", ";
1996   }
1997   cout<<endl;
1998   cout<<"C3 errors:"<<endl;
1999   for(int ii=0; ii<10; ii++){
2000     cout<<Cterm1->GetBinError(ii+1)<<", ";
2001   }
2002   cout<<endl;
2003   TH1D *C3_Sys=(TH1D*)Cterm1->Clone();
2004   for(int ii=1; ii<10; ii++){
2005     if(C3_base[ii] ==0) {
2006       C3_Sys->SetBinContent(ii+1, 100.);
2007       C3_Sys->SetBinError(ii+1, 100.);
2008       continue;
2009     }
2010     C3_Sys->SetBinContent(ii+1, (C3_base[ii]-C3_Sys->GetBinContent(ii+1))/C3_base[ii]);
2011     C3_Sys->SetBinError(ii+1, sqrt(fabs(pow(C3_Sys->GetBinError(ii+1),2) - pow(C3_base_e[ii],2))));
2012   }
2013   gStyle->SetOptFit(111);
2014   C3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2015   C3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2016   C3_Sys->GetYaxis()->SetTitleOffset(2);
2017   C3_Sys->SetMinimum(-0.01);
2018   C3_Sys->SetMaximum(0.01);
2019   //C3_Sys->Draw();
2020   TF1 *C3lineSys=new TF1("C3lineSys","pol3",0,0.099);
2021   //C3_Sys->Fit(C3lineSys,"MEQ","",0,0.099);
2022
2023   if(MCcase){
2024     // Plotting +++ added with --- for HIJING
2025     TLegend *legendC3Hijing = new TLegend(.5,.15,.9,.3,NULL,"brNDC");
2026     legendC3Hijing->SetBorderSize(1);
2027     legendC3Hijing->SetTextSize(.03);// small .03; large .06
2028     legendC3Hijing->SetFillColor(0);
2029
2030     Cterm1->Add(C3baseHisto);
2031     Cterm1->Scale(0.5);
2032     Cterm1->SetMinimum(0.9);
2033     Cterm1->SetMaximum(1.1);
2034     Cterm1->Draw();
2035     legendC3Hijing->AddEntry(Cterm1,"same-charge, HIJING","p");
2036     legendC3Hijing->Draw("same");
2037   }
2038   */
2039   /*
2040   ////////// c3 systematics
2041   // c3 --- base, M0, (New Standard: 0.04 TTC )
2042   double c3_base[10]={0, 1.86014, 1.47533, 1.23733, 1.09944, 1.04145, 1.01693, 1.00715, 1.00253, 1.00111};
2043   double c3_base_e[10]={0, 0.104645, 0.0120917, 0.00333303, 0.00118126, 0.0006692, 0.000405246, 0.000274163, 0.000198507, 0.000150258};
2044
2045   cout<<"c3 values:"<<endl;
2046   for(int ii=0; ii<10; ii++){
2047     cout<<c3_hist->GetBinContent(ii+1)<<", ";
2048   }
2049   cout<<endl;
2050   cout<<"c3 errors:"<<endl;
2051   for(int ii=0; ii<10; ii++){
2052     cout<<c3_hist->GetBinError(ii+1)<<", ";
2053   }
2054   cout<<endl;
2055   TH1D *c3_Sys=(TH1D*)c3_hist->Clone();
2056   for(int ii=1; ii<10; ii++){
2057     if(c3_base[ii] ==0) {
2058       c3_Sys->SetBinContent(ii+1, 100.);
2059       c3_Sys->SetBinError(ii+1, 100.);
2060       continue;
2061     }
2062     c3_Sys->SetBinContent(ii+1, (c3_base[ii]-c3_Sys->GetBinContent(ii+1))/c3_base[ii]);
2063     c3_Sys->SetBinError(ii+1, sqrt(fabs(pow(c3_Sys->GetBinError(ii+1),2) - pow(c3_base_e[ii],2))));
2064   }
2065   gStyle->SetOptFit(111);
2066   c3_Sys->GetXaxis()->SetRangeUser(0,0.099);
2067   c3_Sys->GetYaxis()->SetTitle("(C_{3}^{t1}-C_{3}^{t2})/C_{3}^{t1}");
2068   c3_Sys->GetYaxis()->SetTitleOffset(2);
2069   c3_Sys->SetMinimum(-0.01);
2070   c3_Sys->SetMaximum(0.01);
2071   c3_Sys->Draw();
2072   TF1 *c3lineSys=new TF1("c3lineSys","pol3",0,0.099);
2073   c3_Sys->Fit(c3lineSys,"MEQ","",0,0.099);
2074   */
2075
2076
2077
2078   /*TPad *pad1 = new TPad("pad1","pad1",0.0,0.0,1.,1.);
2079   gPad->SetGridx(0);
2080   gPad->SetGridy(0);
2081   gPad->SetTickx();
2082   gPad->SetTicky();
2083   pad1->SetTopMargin(0.02);//0.05
2084   //pad1->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2085   pad1->SetRightMargin(0.01);//1e-2
2086   pad1->SetBottomMargin(0.06);//0.12
2087   pad1->Divide(1,2,0,0);
2088   pad1->Draw();
2089   pad1->cd(1);
2090   gPad->SetLeftMargin(0.14);
2091   gPad->SetRightMargin(0.02);
2092   if(SameCharge){
2093     Coul_Omega0->SetMinimum(0.48);
2094     Coul_Omega0->SetMaximum(1.01);
2095   }else{
2096     Coul_Omega0->SetMinimum(0.99);
2097     Coul_Omega0->SetMaximum(1.32);
2098   }
2099   Coul_Omega0->SetMarkerStyle(20);
2100   Coul_Omega0->SetMarkerColor(4);
2101   Coul_GRiverside->SetMarkerStyle(25);
2102   Coul_GRiverside->SetMarkerColor(2);
2103   Coul_Riverside->SetMarkerStyle(23);
2104   Coul_Omega0->Draw();
2105   //Coul_Riverside->Draw("same");
2106   Coul_GRiverside->Draw("same");
2107   legend2->AddEntry(Coul_Omega0,"#Omega_{0}","p");
2108   legend2->AddEntry(Coul_GRiverside,"Generalized Riverside","p");
2109   //legend2->AddEntry(Coul_Riverside,"Riverside","p");
2110   legend2->Draw("same");
2111   TLatex *K3Label = new TLatex(-0.011,0.92,"K_{3}");// -0.011,0.92 (ss), 1.26 (os)
2112   K3Label->SetTextSize(0.08);
2113   K3Label->SetTextAngle(90);
2114   K3Label->Draw();
2115   //
2116   pad1->cd(2);
2117   gPad->SetLeftMargin(0.14);
2118   gPad->SetRightMargin(0.02);
2119   gPad->SetBottomMargin(0.13);
2120   TH1D *K3_compOmega0 = (TH1D*)Coul_Omega0->Clone();
2121   TH1D *K3_compGRS = (TH1D*)Coul_GRiverside->Clone();
2122   for(int ii=0; ii<K3_compOmega0->GetNbinsX(); ii++){
2123     K3_compOmega0->SetBinContent(ii+1, K3_compOmega0->GetBinContent(ii+1)-1.0);
2124     K3_compGRS->SetBinContent(ii+1, K3_compGRS->GetBinContent(ii+1)-1.0);
2125   }
2126   K3_compOmega0->Add(K3_compGRS,-1);
2127   K3_compOmega0->Divide(K3_compGRS);
2128   K3_compOmega0->SetMarkerStyle(20);
2129   K3_compOmega0->SetMarkerColor(1);
2130   K3_compOmega0->SetLineColor(1);
2131   K3_compOmega0->SetMinimum(-0.12);// -0.021
2132   K3_compOmega0->SetMaximum(0.12);// 0.021
2133   K3_compOmega0->SetBinContent(1,-100);
2134   K3_compOmega0->Draw();
2135   TLatex *RatioLabel = new TLatex(-.011,-0.05,"#Delta (K_{3}-1)/(K_{3}(#Omega_{0})-1)");// -0.011,0.04
2136   RatioLabel->SetTextSize(0.08);
2137   RatioLabel->SetTextAngle(90);
2138   RatioLabel->Draw();
2139   TLatex *Q3Label = new TLatex(.065,-0.147,"Q_{3} (GeV/c)");//.065,-0.147
2140   Q3Label->SetTextSize(0.08);
2141   Q3Label->Draw();
2142   */
2143
2144
2145   /////////////////////////////
2146   // 3d C3QS fitting
2147   /*
2148   for(int i=0; i<BINRANGE_C2global; i++){
2149     for(int j=0; j<BINRANGE_C2global; j++){
2150       for(int k=0; k<BINRANGE_C2global; k++){
2151         C3ssFitting[i][j][k] = C3QS_3d->GetBinContent(i+1,j+1,k+1);
2152         C3ssFitting_e[i][j][k] = C3QS_3d->GetBinError(i+1,j+1,k+1);
2153       }
2154     }
2155   }
2156
2157   TF3 *C3ssFitResult=new TF3("C3ssFitResult",C3ssFitFunction,0,0.2,0,0.2,0,0.2,npar);
2158   for(int i=0; i<npar; i++) {
2159     if(i<=5) C3ssFitResult->FixParameter(i, OutputPar[i]);
2160     else C3ssFitResult->FixParameter(i, OutputPar[i]-OutputPar_e[i]);
2161   }
2162
2163
2164   const int NbinsC3_1d=Q3BINS;
2165   TH1D *C3fit_1d=new TH1D("C3fit_1d","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2166   TH1D *C3fit_1d_den=new TH1D("C3fit_1d_den","",num_QS->GetNbinsX(),num_QS->GetXaxis()->GetBinLowEdge(1),num_QS->GetXaxis()->GetBinUpEdge(num_QS->GetNbinsX()));
2167   double C3_err[NbinsC3_1d]={0};
2168
2169
2170   for(int i=0; i<C3QS_3d->GetNbinsX(); i++){// bin number
2171     for(int j=0; j<C3QS_3d->GetNbinsY(); j++){// bin number
2172       for(int k=0; k<C3QS_3d->GetNbinsZ(); k++){// bin number
2173
2174         if(i==0 || j==0 || k==0) continue;
2175         //if(i==1 || j==1 || k==1) continue;
2176         double Q3 = sqrt(pow(AvgQinvSS[i],2) + pow(AvgQinvSS[j],2) + pow(AvgQinvSS[k],2));
2177         C3fit_1d->Fill(Q3, C3ssFitResult->Eval(AvgQinvSS[i],AvgQinvSS[j],AvgQinvSS[k])*Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2178         C3fit_1d_den->Fill(Q3, Combinatorics_3d->GetBinContent(i+1,j+1,k+1));
2179       }
2180     }
2181   }
2182
2183   C3fit_1d->Divide(C3fit_1d_den);
2184   for(int q3bin=1; q3bin<=num_QS->GetNbinsX(); q3bin++) C3fit_1d->SetBinError(q3bin, 0.001);// non-zero to display marker
2185   C3fit_1d->SetMarkerStyle(22);
2186   C3fit_1d->SetLineColor(4);
2187   //C3fit_1d->Draw("same");
2188
2189   TF1 *lineC3 = new TF1("lineC3","6",0,10);
2190   lineC3->SetLineColor(4);
2191   lineC3->SetLineStyle(2);
2192   //lineC3->Draw("same");
2193   //legend2->AddEntry(lineC3,"Chaotic Limit C_{3}^{QS}","l");
2194
2195   TF1 *linec3 = new TF1("linec3","3",0,10);
2196   linec3->SetLineColor(2);
2197   linec3->SetLineStyle(2);
2198   //linec3->Draw("same");
2199   //legend2->AddEntry(linec3,"Chaotic Limit #font[12]{c}_{3}","l");
2200   */
2201
2202   //legend2->Draw("same");
2203
2204
2205   /////////////////////////////////////////////////////////////////
2206   /////////////////////////////////////////////////////////////////
2207
2208   TCanvas *can3 = new TCanvas("can3", "can3",1600,0,700,700);
2209   can3->SetHighLightColor(2);
2210   can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
2211   //gStyle->SetOptFit(0111);
2212   can3->SetFillColor(10);
2213   can3->SetBorderMode(0);
2214   can3->SetBorderSize(2);
2215   can3->SetGridx();
2216   can3->SetGridy();
2217   can3->SetFrameFillColor(0);
2218   can3->SetFrameBorderMode(0);
2219   can3->SetFrameBorderMode(0);
2220   can3->cd();
2221   TPad *pad3 = new TPad("pad3","pad3",0.0,0.0,1.,1.);
2222   gPad->SetGridx(1);
2223   gPad->SetGridy(1);
2224   gPad->SetTickx();
2225   gPad->SetTicky();
2226   pad3->SetTopMargin(0.02);//0.05
2227   //pad3->SetLeftMargin(0.13);//.14 for wide title, .10 for narrow title, 0.08 for DeltaQ
2228   pad3->SetRightMargin(0.01);//1e-2
2229   pad3->SetBottomMargin(0.06);//0.12
2230   pad3->Draw();
2231   pad3->cd(1);
2232   gPad->SetLeftMargin(0.14);
2233   gPad->SetRightMargin(0.02);
2234
2235   TF1 *QuarticFit = new TF1("QuarticFit","[0]*(1-[1]*pow(x,4))",0,.1);
2236   QuarticFit->SetParameter(0,2); QuarticFit->SetParameter(1,0);
2237   QuarticFit->SetLineColor(4);
2238   QuarticFit->SetParName(0,"I(Q^{4})"); QuarticFit->SetParName(1,"a(Q^{4})");
2239   TF1 *QuadraticFit = new TF1("QuadraticFit","[0]*(1-[1]*pow(x,2))",0,.1);
2240   QuadraticFit->SetParameter(0,2); QuadraticFit->SetParameter(1,0);
2241   QuadraticFit->SetParName(0,"I(Q^{2})"); QuadraticFit->SetParName(1,"a(Q^{2})");
2242
2243
2244   if(SameCharge){
2245
2246     r3_num_Q3->SetMinimum(0); r3_num_Q3->SetMaximum(2.5);
2247     r3_num_Q3->GetXaxis()->SetRangeUser(0.01,0.09);
2248     r3_num_Q3->Draw();
2249     /*
2250     r3_num_Q3->Fit(QuarticFit,"IME","",0,0.1);
2251     gPad->Update();
2252     TPaveStats *Quartic_stats=(TPaveStats*)r3_num_Q3->FindObject("stats");
2253     Quartic_stats->SetX1NDC(0.2);
2254     Quartic_stats->SetX2NDC(0.5);
2255     Quartic_stats->SetY1NDC(.8);
2256     Quartic_stats->SetY2NDC(.9);
2257
2258     TH1D *r3_clone=(TH1D*)r3_num_Q3->Clone();
2259     r3_clone->Fit(QuadraticFit,"IME","",0,0.1);
2260     gPad->Update();
2261     TPaveStats *Quadratic_stats=(TPaveStats*)r3_clone->FindObject("stats");
2262     Quadratic_stats->SetX1NDC(0.55);
2263     Quadratic_stats->SetX2NDC(0.85);
2264     Quadratic_stats->SetY1NDC(.8);
2265     Quadratic_stats->SetY2NDC(.9);
2266
2267     QuarticFit->Draw("same");
2268     QuadraticFit->Draw("same");
2269     Quartic_stats->Draw("same");
2270     Quadratic_stats->Draw("same");
2271     */
2272     /*
2273     //////////
2274     // r3 Sys
2275     cout<<"r3 values:"<<endl;
2276     for(int ii=0; ii<10; ii++){
2277       cout<<r3_num_Q3->GetBinContent(ii+1)<<", ";
2278     }
2279     cout<<endl;
2280     cout<<"r3 errors:"<<endl;
2281     for(int ii=0; ii<10; ii++){
2282       cout<<r3_num_Q3->GetBinError(ii+1)<<", ";
2283     }
2284     cout<<endl;
2285
2286     // M0,pi- base values with TTC=0.035 (old standard)
2287     //double r3base[10]={0, 2.02584, 1.82659, 1.85384, 1.84294, 1.82431, 1.72205, 1.5982, 1.17379, 0.881535};
2288     //double r3base_e[10]={0, 0.15115, 0.038311, 0.0231936, 0.0209217, 0.0291552, 0.0406093, 0.0624645, 0.0933104, 0.122683};
2289     // M0,pi+ base values with TTC=0.035 (old standard)
2290     //double r3base[10]={0, 1.84738, 1.81805, 1.80057, 1.82563, 1.76585, 1.6749, 1.37454, 1.37558, 0.74562};
2291     //double r3base_e[10]={0, 0.147285, 0.0381688, 0.0231643, 0.0209571, 0.0292255, 0.0407515, 0.0627228, 0.0937758, 0.123392};
2292     // M1,pi- base values with TTC=0.035 (old standard)
2293     //double r3base[10]={0, 1.72641, 1.82547, 1.79747, 1.84141, 1.83557, 1.7994, 1.67511, 1.44238, 1.03004};
2294     //double r3base_e[10]={0, 0.192385, 0.0471262, 0.0270385, 0.0229849, 0.0302041, 0.0397806, 0.0595701, 0.0902909, 0.123976};
2295     // M2,pi- base values with TTC=0.035 (old standard)
2296     //double r3base[10]={0, 1.61726, 1.58973, 1.7834, 1.8914, 1.75445, 1.76511, 1.81367, 1.55617, 1.17797};
2297     //double r3base_e[10]={0, 0.417924, 0.0973277, 0.0527296, 0.0421998, 0.0522792, 0.065144, 0.0923891, 0.135867, 0.190967};
2298     // M9,pi- base values with TTC=0.035 (old standard)
2299     //double r3base[10]={0, 2.51235, 2.68042, 1.85333, 1.66794, 1.39832, 1.717, 1.5337, 1.80686, 1.42286};
2300     //double r3base_e[10]={0, 3.09403, 0.650924, 0.272827, 0.160562, 0.145494, 0.137767, 0.146212, 0.164592, 0.189126};
2301
2302     // M0 PosB, pi- base values with TTC=0.035 (old standard)
2303     //double r3base[10]={0, 1.81361, 1.84779, 1.81693, 1.8252, 1.78366, 1.71822, 1.41243, 1.31656, 0.719842};
2304     //double r3base_e[10]={0, 0.158656, 0.0405, 0.0244075, 0.0219925, 0.0306193, 0.0426494, 0.0655837, 0.0980259, 0.129011};
2305     // M1 PosB, pi- base values with TTC=0.035 (old standard)
2306     //double r3base[10]={0, 1.7208, 1.79016, 1.76858, 1.85008, 1.80416, 1.67274, 1.6501, 1.45098, 1.1052};
2307     //double r3base_e[10]={0, 0.202391, 0.0496228, 0.0282103, 0.0239023, 0.0313499, 0.041234, 0.061679, 0.0934039, 0.128249};
2308
2309     // M1 PID1p5 special case (M0 MomRes used),pi- base values with TTC=0.035 (old standard)
2310     //double r3base[10]={0, 1.77035, 1.82471, 1.77934, 1.80615, 1.77819, 1.71105, 1.60505, 1.44563, 1.09473};
2311     //double r3base_e[10]={0, 0.193326, 0.0471058, 0.0269928, 0.0229581, 0.0301777, 0.0397553, 0.0595472, 0.0902768, 0.123968};
2312
2313     // M0,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2314     //double r3base[10]={0, 1.93236, 1.77265, 1.71324, 1.63426, 1.56167, 1.38348, 1.25064, 0.897615, 0.691928};
2315     //double r3base_e[10]={0, 0.221991, 0.0433103, 0.023234, 0.0187678, 0.0243522, 0.0318668, 0.0459759, 0.066614, 0.0873847};
2316     // M1,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2317     //double r3base[10]={0, 1.46815, 1.76361, 1.68282, 1.6334, 1.56039, 1.44628, 1.27343, 1.01443, 0.691937};
2318     //double r3base_e[10]={0, 0.277164, 0.053134, 0.0271048, 0.0207043, 0.0253691, 0.0315642, 0.0443283, 0.0641466, 0.0858159};
2319     // M2,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2320     //double r3base[10]={0, 1.43331, 1.72759, 1.71155, 1.68957, 1.53846, 1.45209, 1.40698, 1.16202, 1.04672};
2321     //double r3base_e[10]={0, 0.594776, 0.109832, 0.0531013, 0.0383514, 0.0444564, 0.0525963, 0.0705969, 0.0998679, 0.135742};
2322     // M9,pi- base values with TTC=0.04 (New standard, r3 Interp fix, GRS)
2323     //double r3base[10]={0, 6.98181, 2.42845, 1.71962, 1.57641, 1.23102, 1.49543, 1.48063, 1.50938, 1.27714};
2324     //double r3base_e[10]={0, 7.09254, 0.731841, 0.277336, 0.152199, 0.130775, 0.118981, 0.121864, 0.132875, 0.148409};
2325     //
2326     // M0,pi- base values with 3sigma
2327     //double r3base[10]={0, 1.71114, 1.72142, 1.74749, 1.69154, 1.60411, 1.47805, 1.27561, 1.02355, 0.769853};
2328     //double r3base_e[10]={0, 0.251731, 0.0398122, 0.0196354, 0.0153904, 0.0197609, 0.0258499, 0.0377931, 0.0560878, 0.0774703};
2329     // M1,pi- base values with 3sigma
2330     //double r3base[10]={0, 2.02471, 1.78737, 1.68681, 1.69639, 1.62634, 1.57337, 1.42229, 1.28352, 0.958121};
2331     //double r3base_e[10]={0, 0.314475, 0.048507, 0.022855, 0.0169939, 0.0207322, 0.0260424, 0.0375715, 0.0569927, 0.0830874};
2332     // M2,pi- base values with 3sigma
2333     //double r3base[10]={0, 1.50807, 1.75138, 1.71876, 1.76974, 1.6322, 1.5883, 1.5079, 1.4551, 1.49369};
2334     //double r3base_e[10]={0, 0.658863, 0.0979719, 0.0436736, 0.0308217, 0.035654, 0.0425953, 0.0585648, 0.0865817, 0.128056};
2335     // M9,pi- base values with 3sigma
2336     //double r3base[10]={0, 1.17509, 1.5117, 1.98348, 1.63032, 1.57195, 1.48828, 1.48218, 1.53689, 1.39533};
2337     //double r3base_e[10]={0, 4.89181, 0.60212, 0.230712, 0.11983, 0.102686, 0.0937373, 0.0971281, 0.107949, 0.123923};
2338     //
2339     // M0,pi-,EDbin=0, Femto Plus+Minus
2340     //double r3base[10]={0, 1.77055, 1.75926, 1.7711, 1.67159, 1.57851, 1.41187, 1.11004, 0.870851, 0.442019};
2341     //double r3base_e[10]={0, 0.0688956, 0.0232182, 0.0156323, 0.0155991, 0.0225749, 0.0324584, 0.0516114, 0.0797203, 0.112197};
2342     // M1,pi-,EDbin=0, Femto Plus+Minus
2343     //double r3base[10]={0, 1.83169, 1.79042, 1.7275, 1.67116, 1.6032, 1.5423, 1.30426, 1.13238, 0.656902};
2344     //double r3base_e[10]={0, 0.0891247, 0.0286553, 0.0182713, 0.0171956, 0.0235918, 0.0324304, 0.0510863, 0.0823637, 0.125496};
2345     // M2,pi-,EDbin=0, Femto Plus+Minus
2346     //double r3base[10]={0, 1.60769, 1.81565, 1.75909, 1.74026, 1.62865, 1.58101, 1.50583, 1.3503, 1.27367};
2347     //double r3base_e[10]={0, 0.188422, 0.0581482, 0.0348698, 0.0309363, 0.0401393, 0.0519907, 0.0771148, 0.120583, 0.187679};
2348     // M3,pi-,EDbin=0, Femto Plus+Minus
2349     double r3base[10]={0, 1.52778, 1.81442, 1.78982, 1.69978, 1.70934, 1.59516, 1.56605, 1.45676, 1.17468};
2350     double r3base_e[10]={0, 0.288329, 0.0882351, 0.0510576, 0.0432015, 0.0536294, 0.0667149, 0.0950229, 0.143021, 0.217322};
2351     // M8,pi-,EDbin=0, Femto Plus+Minus
2352     //double r3base[10]={0, 2.27921, 2.40794, 1.72409, 1.66853, 1.70966, 1.52681, 1.64425, 1.18747, 1.35988};
2353     //double r3base_e[10]={0, 1.08315, 0.289559, 0.136951, 0.0919134, 0.0900869, 0.0901754, 0.101887, 0.121289, 0.148063};
2354     // M9,pi-,EDbin=0, Femto Plus+Minus
2355     //double r3base[10]={0, 2.1271, 2.12837, 1.95071, 1.4832, 1.50308, 1.41068, 1.29826, 1.10739, 0.755899};
2356     //double r3base_e[10]={0, 1.52325, 0.385705, 0.181328, 0.116165, 0.10896, 0.105705, 0.115163, 0.131994, 0.156299};
2357
2358     //
2359     // M9,pi+ (No MRC)
2360     double r3base_noMRC[10]={0, 4.16425, 0.429681, 1.56506, 1.64596, 1.73785, 1.57181, 1.51971, 1.59096, 1.54403};
2361     double r3base_noMRC_e[10]={0, 2.00855, 0.437928, 0.187872, 0.1087, 0.097367, 0.0893871, 0.0933358, 0.103603, 0.11839};
2362     // M9,pi+ (MRC normal)
2363     double r3base_MRC[10]={0, 3.84986, 0.558354, 1.60733, 1.67855, 1.75362, 1.59637, 1.54972, 1.62045, 1.57307};
2364     double r3base_MRC_e[10]={0, 1.93428, 0.402333, 0.171664, 0.0996222, 0.0903509, 0.0840615, 0.0889731, 0.100023, 0.115704};
2365     // M9,pi+ (MRC 50% larger)
2366     double r3base_MRC_over[10]={0, 3.70355, 0.614999, 1.62641, 1.69371, 1.76209, 1.60931, 1.56565, 1.63663, 1.58935};
2367     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};
2368
2369
2370
2371     TH1D *r3_noMRC=(TH1D*)r3_num_Q3->Clone();
2372     TH1D *r3_MRC=(TH1D*)r3_num_Q3->Clone();
2373     TH1D *r3_MRC_over=(TH1D*)r3_num_Q3->Clone();
2374
2375     TH1D *r3Sys=(TH1D*)r3_num_Q3->Clone();
2376     for(int ii=1; ii<10; ii++){
2377       double Stat_RelativeError = sqrt(fabs(pow(r3_num_Q3->GetBinError(ii+1),2) - pow(r3base_e[ii],2)));// + for independent stat's
2378       r3Sys->SetBinContent(ii+1, (r3base[ii]-r3_num_Q3->GetBinContent(ii+1))/r3base[ii]);
2379       r3Sys->SetBinError(ii+1, Stat_RelativeError);
2380       //
2381       r3_noMRC->SetBinContent(ii+1, r3base_noMRC[ii]);
2382       r3_MRC->SetBinContent(ii+1, r3base_MRC[ii]);
2383       r3_MRC_over->SetBinContent(ii+1, r3base_MRC_over[ii]);
2384       r3_noMRC->SetBinError(ii+1, r3base_noMRC_e[ii]);
2385       r3_MRC->SetBinError(ii+1, r3base_MRC_e[ii]);
2386       r3_MRC_over->SetBinError(ii+1, r3base_MRC_over_e[ii]);
2387     }
2388     r3_noMRC->SetMarkerStyle(20); r3_noMRC->SetMarkerColor(1); r3_noMRC->SetLineColor(1);
2389     r3_MRC->SetMarkerStyle(23); r3_MRC->SetMarkerColor(4); r3_MRC->SetLineColor(4);
2390     r3_MRC_over->SetMarkerStyle(24); r3_MRC_over->SetMarkerColor(2); r3_MRC_over->SetLineColor(2);
2391
2392     //r3_MRC->Draw();
2393     //r3_noMRC->Draw("same");
2394     //r3_MRC_over->Draw("same");
2395     TLegend *legendMRC = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2396     legendMRC->SetBorderSize(1);
2397     legendMRC->SetTextSize(.03);// small .03; large .06
2398     legendMRC->SetFillColor(0);
2399     //legendMRC->AddEntry(r3_MRC,"Momentum Resolution Corrected","p");
2400     //legendMRC->AddEntry(r3_noMRC,"No Correction","p");
2401     //legendMRC->AddEntry(r3_MRC_over,"50% Increased Correction","p");
2402     //legendMRC->Draw("same");
2403
2404     r3Sys->GetYaxis()->SetTitle("(r_{3}^{t1}-r_{3}^{t2})/r_{3}^{t1}");
2405     r3Sys->GetYaxis()->SetTitleOffset(1.8);
2406     r3Sys->SetMinimum(-0.1);
2407     r3Sys->SetMaximum(0.1);
2408     r3Sys->GetXaxis()->SetRangeUser(0,0.099);
2409     r3Sys->Draw();
2410     TF1 *lineFit=new TF1("lineFit","pol0",0,0.02);
2411     gStyle->SetOptFit(111);
2412     r3Sys->Fit(lineFit,"MEQ","",0,0.1);
2413     */
2414
2415
2416     /*
2417     ////////////////////////////////////////////////////////////////////////
2418     // The STAR method
2419     ReadCoulCorrections(0, 5, 2, 10);// Change to Gaussian R=5 fm calculation
2420
2421     TLegend *legendSTAR = new TLegend(.45,.15,.95,.3,NULL,"brNDC");
2422     legendSTAR->SetBorderSize(1);
2423     legendSTAR->SetTextSize(.03);// small .03; large .06
2424     legendSTAR->SetFillColor(0);
2425
2426     TH1D *r3_STAR = new TH1D("r3_STAR","",Q3BINS,0,0.2);
2427     TH1D *r3_den_STAR = new TH1D("r3_den_STAR","",Q3BINS,0,0.2);
2428     double r3_STAR_num_e[20]={0};
2429     double r3_STAR_den_e[20]={0};
2430
2431     for(int i=2; i<=20; i++){// bin number
2432       double Q12 = (i-0.5)*(0.005);
2433       int Q12bin = int(Q12/0.005)+1;
2434       //
2435       for(int j=2; j<=20; j++){// bin number
2436         double Q13 = (j-0.5)*(0.005);
2437         int Q13bin = int(Q13/0.005)+1;
2438         //
2439         for(int k=2; k<=20; k++){// bin number
2440           double Q23 = (k-0.5)*(0.005);
2441           int Q23bin = int(Q23/0.005)+1;
2442           //
2443           double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
2444           int Q3bin = Cterm1_noMRC->GetXaxis()->FindBin(Q3);
2445           if(Q3bin>19) continue;
2446
2447           if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;
2448           if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;
2449
2450           double K2_12 = CoulCorr2(CP12, Q12);
2451           double K2_13 = CoulCorr2(CP13, Q13);
2452           double K2_23 = CoulCorr2(CP23, Q23);
2453
2454           double r3num = Cterm1_noMRC->GetBinContent(Q3bin)/(K2_12*K2_13*K2_23) - 1;
2455           r3num -= C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1;
2456           r3num -= C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1;
2457           r3num -= C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1;
2458           double r3den = (C2ssRaw->GetBinContent(Q12bin)/(K2_12) - 1);
2459           r3den *= (C2ssRaw->GetBinContent(Q13bin)/(K2_13) - 1);
2460           r3den *= (C2ssRaw->GetBinContent(Q23bin)/(K2_23) - 1);
2461           if(r3den<0) continue;
2462           r3den = sqrt(r3den);
2463           //
2464           r3_STAR_num_e[Q3bin-1] += pow(Cterm1_noMRC->GetBinError(Q3bin)/(K2_12*K2_13*K2_23),2);
2465           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q12bin)/(K2_12),2);
2466           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q13bin)/(K2_13),2);
2467           r3_STAR_num_e[Q3bin-1] += pow(C2ssRaw->GetBinError(Q23bin)/(K2_23),2);
2468           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);
2469           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);
2470           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);
2471           //
2472           r3_STAR->Fill(Q3, r3num);
2473           r3_den_STAR->Fill(Q3, r3den);
2474         }
2475       }
2476     }
2477
2478
2479     for(int i=0; i<20; i++){
2480       r3_STAR->SetBinError(i+1, sqrt( r3_STAR_num_e[i]));
2481       r3_den_STAR->SetBinError(i+1, sqrt( r3_STAR_den_e[i]));
2482     }
2483     r3_STAR->Divide(r3_den_STAR);
2484     r3_STAR->SetMaximum(5);
2485     r3_STAR->SetMinimum(-10);
2486     r3_STAR->GetXaxis()->SetRangeUser(0,0.099);
2487     r3_STAR->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
2488     r3_STAR->GetYaxis()->SetTitle("r_{3}*");
2489     r3_STAR->SetLineColor(2); r3_STAR->SetMarkerColor(2);
2490     r3_STAR->SetMarkerStyle(kFullStar);
2491     r3_STAR->SetMarkerSize(1.2);
2492     r3_STAR->Draw();
2493     r3_STAR->Fit(QuarticFit,"IME","",0.01,0.1);
2494     gPad->Update();
2495     TPaveStats *Quartic_stats_STAR=(TPaveStats*)r3_STAR->FindObject("stats");
2496     Quartic_stats_STAR->SetX1NDC(0.2);
2497     Quartic_stats_STAR->SetX2NDC(0.5);
2498     Quartic_stats_STAR->SetY1NDC(.8);
2499     Quartic_stats_STAR->SetY2NDC(.9);
2500
2501     TH1D *r3_clone_STAR=(TH1D*)r3_STAR->Clone();
2502     r3_STAR->Fit(QuadraticFit,"IME","",0,0.1);
2503     gPad->Update();
2504     TPaveStats *Quadratic_stats_STAR=(TPaveStats*)r3_clone_STAR->FindObject("stats");
2505     Quadratic_stats_STAR->SetX1NDC(0.55);
2506     Quadratic_stats_STAR->SetX2NDC(0.85);
2507     Quadratic_stats_STAR->SetY1NDC(.8);
2508     Quadratic_stats_STAR->SetY2NDC(.9);
2509
2510     QuarticFit->Draw("same");
2511     QuadraticFit->Draw("same");
2512     Quartic_stats_STAR->Draw("same");
2513     Quadratic_stats_STAR->Draw("same");
2514     //Cterm1_noMRC->Draw();
2515
2516     /////////////////////////////////////////////////////////////////////////
2517     if(SameCharge) r3_num_Q3->Draw("same");
2518     legendSTAR->AddEntry(r3_STAR,"STAR method","p");
2519     legendSTAR->AddEntry(r3_num_Q3,"ALICE method","p");
2520     legendSTAR->Draw("same");
2521     */
2522   }// SameCharge
2523
2524   //c3_hist_STAR->Draw();
2525   //Cterm1_noMRC->Draw();
2526
2527   //r3_num_Q3->SetMarkerStyle(20);
2528   //r3_num_Q3->GetXaxis()->SetRangeUser(0,0.1);
2529   //r3_num_Q3->SetMinimum(-7);
2530   //r3_num_Q3->SetMaximum(2.2);
2531
2532
2533   /*
2534     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};
2535   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};
2536   TH1D *C3_star=(TH1D*)Cterm1->Clone();
2537   C3_star->Add(C3_star,-1);
2538   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]);}
2539   C3_star->SetMarkerColor(4);
2540   C3_star->SetLineColor(4);
2541   //C3_star->Draw("same");
2542   //legend1->AddEntry(C3_star,"---; 200 GeV","p");
2543
2544   //legend1->Draw("same");
2545
2546   //TGraph *gr_MomRes_term1=new TGraph(10,
2547   */
2548
2549
2550
2551
2552
2553   //////////////////////////////////////////////////////////////////////////////////////
2554
2555
2556
2557
2558
2559
2560   if(SaveToFile){
2561     TString *savefilename = new TString("OutFiles/OutFile");
2562     if(SameCharge) savefilename->Append("SC");
2563     else savefilename->Append("MC");
2564     //
2565     if(CHARGE==1) savefilename->Append("Pos");
2566     else savefilename->Append("Neg");
2567     //
2568     if(IncludeG) savefilename->Append("YesG");
2569     else savefilename->Append("NoG");
2570     //
2571     if(IncludeEWfromTherm) savefilename->Append("EWfromTherm");
2572     if(IncludeEW) savefilename->Append("EW");
2573     if(GRS) savefilename->Append("GRS");
2574     else savefilename->Append("Omega0");
2575
2576     if(EDbin==0) savefilename->Append("Kt3_1");
2577     else savefilename->Append("Kt3_2");
2578
2579     savefilename->Append("_Kt");
2580     *savefilename += Ktbin;
2581     savefilename->Append("_M");
2582     *savefilename += Mbin;
2583     savefilename->Append(".root");
2584     TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
2585     can1->Write("can1");
2586     C2_ss->Write("C2_ss");
2587     C2_os->Write("C2_os");
2588     C2_ss_Momsys->Write("C2_ss_Momsys");
2589     C2_os_Momsys->Write("C2_os_Momsys");
2590     C2_ss_Ksys->Write("C2_ss_Ksys");
2591     C2_os_Ksys->Write("C2_os_Ksys");
2592     fitC2ss->Write("fitC2ss");
2593     fitC2os->Write("fitC2os");
2594     if(IncludeEWfromTherm) {
2595       fitC2ssTherm->Write("fitC2ssTherm");
2596       fitC2osTherm->Write("fitC2osTherm");
2597       fitC2ssThermGaus->Write("fitC2ssThermGaus");
2598       fitC2osThermGaus->Write("fitC2osThermGaus");
2599     }
2600     C2Therm_ss->Write("C2Therm_ss");
2601     C2Therm_os->Write("C2Therm_os");
2602     K2_ss->Write("K2_ss");
2603     K2_os->Write("K2_os");
2604     Cterm1->Write("C3");
2605     Coul_Omega0->Write("Coul_Omega0");
2606     Coul_GRiverside->Write("Coul_GRiverside");
2607     GenSignalExpected_num->Write("C3_EWexpectation");
2608     num_QS->Write("num_QS");
2609     c3_hist->Write("c3");
2610     Combinatorics_1d->Write("Combinatorics_1d");
2611     ChiSquaredNDF->Write("ChiSquaredNDF");
2612     r3_num_Q3->Write("r3_Q3");
2613
2614     //
2615     savefile->Close();
2616   }
2617
2618
2619 }
2620
2621 void ReadCoulCorrections(int RVal, int bVal, int kt){
2622   ///////////////////////
2623
2624   TString *fname;
2625   if(FileSetting!=6 && SourceType==0) fname = new TString("KFile.root");
2626   else fname = new TString("KFile_Gauss.root");
2627
2628   TFile *File=new TFile(fname->Data(),"READ");
2629   if(bVal!=2 && bVal!=3 && bVal!=5 && bVal!=7 && bVal!=8 && bVal!=9) cout<<"Therminator bVal not acceptable in 2-particle Coulomb read"<<endl;
2630
2631   if(kt==10){// kt integrated
2632     TH2D *tempT_ss = (TH2D*)File->Get("K2ssT");
2633     TH2D *tempT_os = (TH2D*)File->Get("K2osT");
2634     CoulCorr2SS = (TH1D*)tempT_ss->ProjectionY("CoulCorr2SS",bBin, bBin);
2635     CoulCorr2OS = (TH1D*)tempT_os->ProjectionY("CoulCorr2OS",bBin, bBin);
2636   }else{
2637     if(kt < 1 || kt > 6) cout<<"kt bin out of range in 2-particle Coulomb read"<<endl;
2638     TH3D *tempT3_ss = (TH3D*)File->Get("K2ssT_kt");
2639     TH3D *tempT3_os = (TH3D*)File->Get("K2osT_kt");
2640     CoulCorr2SS = (TH1D*)tempT3_ss->ProjectionZ("CoulCorr2SS",bBin, bBin, kt,kt);
2641     CoulCorr2OS = (TH1D*)tempT3_os->ProjectionZ("CoulCorr2OS",bBin, bBin, kt,kt);
2642   }
2643
2644   CoulCorr2SS->SetDirectory(0);
2645   CoulCorr2OS->SetDirectory(0);
2646   File->Close();
2647 }
2648 double CoulCorr2(int chargeproduct, double Q2){// returns K2
2649   int indexL=0;
2650   int indexH=0;
2651   double slope=0;
2652
2653
2654   indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2655   indexH = indexL+1;
2656
2657   if(indexH >= Nlines) return 1.0;
2658   if(chargeproduct==1){
2659     slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
2660     slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
2661     return slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1);
2662   }else {
2663     slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1));
2664     slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1));
2665     return slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
2666   }
2667
2668   /*
2669   int Qbin=CoulCorr2SS->GetXaxis()->FindBin(Q2);
2670   if(chargeproduct==1) return CoulCorr2SS->GetBinContent(Qbin);
2671   else return CoulCorr2OS->GetBinContent(Qbin);
2672   */
2673 }
2674
2675 //----------------------------------------------------------------------
2676
2677
2678 //________________________________________________________________________
2679 void fcnC2_Global(int& npar, double* deriv, double& f, double par[], int flag){
2680
2681   double qinvSS=0;
2682
2683   double Rch=par[3]/FmToGeV;
2684   double Rcoh=par[4]/FmToGeV;
2685   double Dp=0;
2686   double CSS=0, COS=0;
2687   double SumChi2=0;
2688   double EW=0;
2689   double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2690   NFitPoints_C2global=0;
2691   if(LinkN) par[9]=par[0];// Link N factors
2692
2693   for(int i=1; i<BINRANGE_C2global; i++){
2694
2695     qinvSS = AvgQinvSS[i];
2696
2697     if(!GofP) Dp=fabs(par[2])/(1-fabs(par[2]));// p independent
2698     else Dp = fabs(par[2])/(1-fabs(par[2]));
2699
2700     //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));
2701     //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));
2702     //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);
2703     //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);
2704     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]-qinvSS/2.,2) - 2*MT/Temp);
2705     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(par[10]+qinvSS/2.,2) - 2*MT/Temp);
2706     //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2707     //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2708     double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2709     double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2710
2711     if(!GofP) {Dp1=Dp; Dp2=Dp;}
2712
2713     //
2714     double arg=qinvSS*Rch;
2715     EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2716     EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2717     EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2718     EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2719     //
2720     double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
2721     double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
2722     CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
2723     if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2724     else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
2725     //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
2726
2727     CSS *= par[1]*K2SS[i];
2728     CSS += 1-par[1];
2729     CSS *= par[0];
2730     //
2731     COS = 1;
2732     if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2733     COS *= par[1]*K2OS[i];
2734     COS += 1-par[1];
2735     COS *= par[9];// different Norm
2736     //
2737     if(C2ssFitting[i] > 0){
2738       if(IncludeSS) {
2739         double errorSS = sqrt(pow( (CSS/par[0] - (1-par[1]))/K2SS[i] * K2SS_e[i] * par[0],2) + pow(C2ssFitting_e[i],2));
2740         //double errorSS = C2ssFitting_e[i];
2741         SumChi2 += pow((CSS - C2ssFitting[i])/errorSS,2);
2742         NFitPoints_C2global++;
2743       }
2744     }
2745     if(IncludeOS) {
2746       double errorOS = sqrt(pow( (COS/par[9] - (1-par[1]))/K2OS[i] * K2OS_e[i] * par[9],2) + pow(C2osFitting_e[i],2));
2747       //double errorOS = C2osFitting_e[i];
2748       SumChi2 += pow((COS - C2osFitting[i])/errorOS,2);
2749       NFitPoints_C2global++;
2750     }
2751   }
2752
2753
2754   f = SumChi2;
2755   Chi2_C2global = f;
2756
2757 }
2758 //________________________________________________________________________
2759 double C2ssFitFunction(Double_t *x, Double_t *par)
2760 {
2761   double Rch=par[3]/FmToGeV;
2762   double Rcoh=par[4]/FmToGeV;
2763   double Dp=0;
2764   int qbin = int(fabs(x[0]/0.005));
2765   if(qbin >= BINRANGE_C2global) return 1.0;
2766   double qinvSS = AvgQinvSS[qbin];
2767   double MT = sqrt(pow(0.25 + 0.1*Ktbin_GofP,2) + pow(masspiC,2));
2768
2769   if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
2770   else Dp = fabs(par[2])/(1-fabs(par[2]));
2771
2772   //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));
2773   //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));
2774   //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);
2775   //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);
2776   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2777   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*pow(Rcoh,2)*pow(qinvSS/2.,2) - 2*MT/Temp);
2778   double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2779   double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvSS/2.,2));
2780
2781   if(!GofP) {Dp1=Dp; Dp2=Dp;}
2782   double arg=qinvSS*Rch;
2783   double EW = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2784   EW += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2785   EW += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2786   EW += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2787   double Gaus_coh = exp(-pow(Rcoh*qinvSS,2)/2.);
2788   double Gaus_ch = exp(-pow(Rch*qinvSS,2)/2.) * EW;
2789   double CSS = 1 + pow(1 + Dp*Gaus_coh/Gaus_ch,2)/((1+Dp1)*(1+Dp2)) * exp(-pow(Rch*qinvSS,2))*pow(EW,2);
2790   if(ChargeConstraint) CSS -= 4/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2791   else CSS -= pow(Gaus_coh*Dp,2)/((1+Dp1)*(1+Dp2));
2792   //else CSS -= Dp1/(1+Dp1) * Dp2/(1+Dp2);
2793   CSS *= par[1]*K2SS[qbin];
2794   CSS += 1-par[1];
2795   CSS *= par[0];
2796
2797   //cout<<qinvSS<<"  "<<Dp1/(1+Dp1) * Dp2/(1+Dp2)<<endl;
2798
2799   return CSS;
2800 }
2801 //________________________________________________________________________
2802 double C2osFitFunction(Double_t *x, Double_t *par)
2803 {
2804   if(LinkN) par[9]=par[0];// Link N factors
2805   double Rcoh=par[4]/FmToGeV;
2806   double Dp=0;
2807   int qbin = int(fabs(x[0]/0.005));
2808   if(qbin >= BINRANGE_C2global) return 1.0;
2809   double qinvOS = AvgQinvOS[qbin];
2810
2811   if(!GofP) Dp = fabs(par[2])/(1-fabs(par[2]));// p independent
2812   else Dp = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(AvgP[Ktbin_GofP-1][qbin],2));
2813   //Dp = fabs(par[2])/(1-fabs(par[2]));
2814   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));
2815   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));
2816   //double Dp1 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
2817   //double Dp2 = fabs(par[2])/(1-fabs(par[2])) * exp(-2*(pow(Rcoh,2)-pow(RT,2))*pow(qinvOS/2.,2));
2818
2819   if(!GofP) {Dp1=Dp; Dp2=Dp;}
2820   double COS = 1;
2821   if(ChargeConstraint && GofP) COS += 1/5.* Dp1/(1+Dp1) * Dp2/(1+Dp2);
2822   COS *= par[1]*K2OS[qbin];
2823   COS += 1-par[1];
2824   COS *= par[9];
2825   return COS;
2826 }
2827 //________________________________________________________________________
2828 void fcnC3_Global(int& npar, double* deriv, double& f, double par[], int flag){
2829
2830   double qinvSS_12=0;
2831   double qinvSS_13=0;
2832   double qinvSS_23=0;
2833
2834   //double D=fabs(par[2])/(1-fabs(par[2]));
2835   double Rch=par[3]/FmToGeV;
2836   //double Rcoh=par[4]/FmToGeV;
2837   double CSS=0;
2838   double SumChi2=0;
2839   double EW_12=0, EW_13=0, EW_23=0;
2840   NFitPoints_C3global=0;
2841
2842   for(int i=0; i<BINRANGE_C2global; i++){
2843     qinvSS_12 = AvgQinvSS[i];
2844     double arg_12=qinvSS_12*Rch;
2845     EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
2846     EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
2847     EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
2848     EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
2849
2850     for(int j=0; j<BINRANGE_C2global; j++){
2851       qinvSS_13 = AvgQinvSS[j];
2852       double arg_13=qinvSS_13*Rch;
2853       EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
2854       EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
2855       EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
2856       EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
2857
2858       for(int k=0; k<BINRANGE_C2global; k++){
2859         qinvSS_23 = AvgQinvSS[k];
2860         double arg_23=qinvSS_23*Rch;
2861         EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
2862         EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
2863         EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
2864         EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
2865         //
2866         if(i==0 || j==0 || k==0) continue;
2867         CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
2868         CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
2869         CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
2870         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)));
2871         //
2872         if(IncludeSS) {
2873           /*double errorSS = sqrt(pow(C3ssFitting_e[i],2));
2874           if(errorSS > 0){
2875             SumChi2 += pow((CSS - C3ssFitting[i])/errorSS,2);
2876             NFitPoints_C3global++;
2877             }*/
2878         }
2879         if(IncludeOS) {
2880           //double errorOS = sqrt(pow(C3osFitting_e[i],2));
2881           //SumChi2 += pow((COS - C3osFitting[i])/errorOS,2);
2882           //NFitPoints_C3global++;
2883         }
2884
2885
2886       }// k
2887     }// j
2888   }// i
2889
2890   f = SumChi2;
2891   Chi2_C3global = f;
2892
2893 }
2894 //______________________________________________________________________
2895 /*void C3ssFitFunction(Double_t *x, Double_t *par){
2896
2897   double qinvSS_12=0, qinvOS_12=0;
2898   double qinvSS_13=0, qinvOS_13=0;
2899   double qinvSS_23=0, qinvOS_23=0;
2900
2901   double D=fabs(par[2])/(1-fabs(par[2]));
2902   double Rch=par[3]/FmToGeV;
2903   double Rcoh=par[4]/FmToGeV;
2904   double CSS=0, COS=0;
2905   double term1=0, term2=0;
2906   double SumChi2=0;
2907   double EW_12=0, EW_13=0, EW_23=0;
2908   NFitPoints_C3global=0;
2909
2910
2911   int qbin_12 = x[0]*1000/5;
2912   qinvSS_12 = AvgQinvSS[qbin_12];
2913   qinvOS_12 = AvgQinvOS[qbin_12];
2914   double arg_12=qinvSS_12*Rch;
2915   EW_12 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_12,3) - 12.*pow(arg_12,1));
2916   EW_12 += par[6]/(24.*pow(2,2))*(16.*pow(arg_12,4) -48.*pow(arg_12,2) + 12);
2917   EW_12 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_12,5) - 160.*pow(arg_12,3) + 120*arg_12);
2918   EW_12 += par[8]/(720.*pow(2,3))*(64.*pow(arg_12,6) - 480.*pow(arg_12,4) + 720.*pow(arg_12,2) - 120);
2919   //
2920   int qbin_13 = x[1]*1000/5;
2921   qinvSS_13 = AvgQinvSS[qbin_13];
2922   qinvOS_13 = AvgQinvOS[qbin_13];
2923   double arg_13=qinvSS_13*Rch;
2924   EW_13 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_13,3) - 12.*pow(arg_13,1));
2925   EW_13 += par[6]/(24.*pow(2,2))*(16.*pow(arg_13,4) -48.*pow(arg_13,2) + 12);
2926   EW_13 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_13,5) - 160.*pow(arg_13,3) + 120*arg_13);
2927   EW_13 += par[8]/(720.*pow(2,3))*(64.*pow(arg_13,6) - 480.*pow(arg_13,4) + 720.*pow(arg_13,2) - 120);
2928   //
2929   int qbin_23 = x[2]*1000/5;
2930   qinvSS_23 = AvgQinvSS[qbin_23];
2931   qinvOS_23 = AvgQinvOS[qbin_23];
2932   double arg_23=qinvSS_23*Rch;
2933   EW_23 = 1 + par[5]/(6.*pow(2,1.5))*(8.*pow(arg_23,3) - 12.*pow(arg_23,1));
2934   EW_23 += par[6]/(24.*pow(2,2))*(16.*pow(arg_23,4) -48.*pow(arg_23,2) + 12);
2935   EW_23 += par[7]/(120.*pow(2,2.5))*(32.*pow(arg_23,5) - 160.*pow(arg_23,3) + 120*arg_23);
2936   EW_23 += par[8]/(720.*pow(2,3))*(64.*pow(arg_23,6) - 480.*pow(arg_23,4) + 720.*pow(arg_23,2) - 120);
2937   //
2938   if(qbin_12==0 || qbin_13==0 || qbin_23==0) return 0;
2939   CSS = 1 + pow(EW_12,2)*exp(-pow(Rch*qinvSS_12,2));
2940   CSS += pow(EW_13,2)*exp(-pow(Rch*qinvSS_13,2));
2941   CSS += pow(EW_23,2)*exp(-pow(Rch*qinvSS_23,2));
2942   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)));
2943   //
2944
2945   return CSS;
2946   }*/
2947 //----------------------------------------------------------------------
2948 void fcnOSL(int& npar, double* deriv, double& f, double par[], int flag){
2949
2950   double qout=0, qside=0, qlong=0;
2951   double Gaus_ch=0;
2952   double C=0;
2953   //double lnL=0;
2954   double SumChi2=0;
2955   NFitPoints_OSL=0;
2956
2957   for(int i=0; i<BINS_OSL-10; i++){// out, 0 to BINS_OSL-10
2958     for(int j=0; j<BINS_OSL-10; j++){// side, 0 to BINS_OSL-10
2959       for(int k=0; k<BINS_OSL-10; k++){// long, 0 to BINS_OSL-10
2960
2961         if(B[i][j][k] < 1. ) continue;
2962
2963         qout = (i+0.5)*0.005;
2964         qside = (j+0.5)*0.005;
2965         qlong = (k+0.5)*0.005;
2966
2967         Gaus_ch = exp(-pow(par[3]*qout/FmToGeV,2) - pow(par[4]*qside/FmToGeV,2) - pow(par[5]*qlong/FmToGeV,2));
2968
2969         //C = par[0]*(1 + par[1]*(K_OSL[i][j][k]-1) + K_OSL[i][j][k]*par[1]*Gaus_ch);
2970         C = par[0]*( (1-par[1]) + par[1]*K_OSL[i][j][k]*(1+Gaus_ch));
2971
2972         if(C < 0) continue;
2973
2974         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);
2975         error += pow((C/par[0] - (1-par[1]))/K_OSL[i][j][k] * K_OSL_e[i][j][k]*par[0],2);
2976         error = sqrt(error);
2977         SumChi2 += pow( (C - A[i][j][k]/B[i][j][k])/error,2);
2978
2979         /*
2980         if(A[i][j][k] >= 1) term1 = C*(A[i][j][k]+B[i][j][k])/(A[i][j][k]*(C+1));
2981         else term1 = 0;
2982         term2 = (A[i][j][k]+B[i][j][k])/(B[i][j][k]*(C+1));
2983
2984         if(term1 > 0.0 && term2 > 0.0){
2985           lnL += A[i][j][k]*log(term1) + B[i][j][k]*log(term2);
2986         }else if(term1==0 && term2 > 0.0){
2987           lnL += B[i][j][k]*log(term2);
2988         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
2989         */
2990         NFitPoints_OSL++;
2991       }
2992     }
2993   }
2994   f = SumChi2;
2995   //f = -2.0*lnL;
2996   Chi2_OSL = f;
2997
2998 }
2999 //________________________________________________________________________
3000 double OSLfitfunction(Double_t *x, Double_t *par)
3001 {
3002   int bin_out = int(1000*(x[0])/5.);
3003   int bin_side = int(1000*(x[1])/5.);
3004   int bin_long = int(1000*(x[2])/5.);
3005
3006   double K = CoulCorr2(+1, avg_q[bin_out][bin_side][bin_long]);
3007   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));
3008   return par[0]*(1 + par[1]*(K-1) + K*par[1]*Gaus_ch);
3009
3010 }
3011 //__________________________________________________________________________
3012
3013 void fcn_3(int& npar, double* deriv, double& f, double par[], int flag){
3014
3015   double q12=0, q13=0, q23=0;
3016   double K12=0, K13=0, K23=0, K123=0;
3017   double G1=0, G2=0, G3=0, G4=0;
3018   double EW12=0, EW13=0, EW23=0;
3019   double C=0;
3020   double term1=0, term2=0;
3021   double lnL=0;
3022
3023   // start from 2-4 MeV bins
3024   for(int i=1; i<BINRANGE_3-10; i++){// q12
3025     for(int j=1; j<BINRANGE_3-10; j++){// q13
3026       for(int k=1; k<BINRANGE_3-10; k++){// q23
3027
3028         if(B_3[i][j][k] < 1. ) continue;
3029
3030         q12 = (i+0.5)*0.002;
3031         q13 = (j+0.5)*0.002;
3032         q23 = (k+0.5)*0.002;
3033
3034         //K12 = 1/CoulCorr_CsorgoMate(2, q12);
3035         //K13 = 1/CoulCorr_CsorgoMate(2, q13);
3036         //K23 = 1/CoulCorr_CsorgoMate(2, q23);
3037         //K12 = 1/Coul_pipi_data(q12);
3038         //K13 = 1/Coul_pipi_data(q13);
3039         //K23 = 1/Coul_pipi_data(q23);
3040         K12=1.0; K13=1.0; K23=1.0;
3041         K123 = K12*K13*K23;
3042         //
3043         //K123 = 1 - pow(par[1],3)*(1-K12*K13*K23);
3044         //K12 = 1 - pow(par[1],3)*(1-K12);
3045         //K13 = 1 - pow(par[1],3)*(1-K13);
3046         //K23 = 1 - pow(par[1],3)*(1-K23);
3047         //
3048         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);
3049         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);
3050         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);
3051         //EW12=1; EW13=1; EW23=1;
3052         //
3053         G1 = K12*exp(-pow(par[2]*q12/FmToGeV,2)/2.)*EW12;
3054         G1 += K13*exp(-pow(par[2]*q13/FmToGeV,2)/2.)*EW13;
3055         G1 += K23*exp(-pow(par[2]*q23/FmToGeV,2)/2.)*EW23;
3056         //
3057         G2 = K12*exp(-pow(par[2]*q12/FmToGeV,2))*pow(EW12,2);
3058         G2 += K13*exp(-pow(par[2]*q13/FmToGeV,2))*pow(EW13,2);
3059         G2 += K23*exp(-pow(par[2]*q23/FmToGeV,2))*pow(EW23,2);
3060         //
3061         G3 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q23*q23)/2.)*EW12*EW23;
3062         G3 += exp(-pow(par[2]/FmToGeV,2)*(q13*q13 + q23*q23)/2.)*EW13*EW23;
3063         G3 += exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13)/2.)*EW12*EW13;
3064         G3 *= K123;
3065         //
3066         G4 = exp(-pow(par[2]/FmToGeV,2)*(q12*q12 + q13*q13 + q23*q23)/2.);
3067         G4 *= EW12*EW13*EW23;
3068         G4 *= K123;
3069         //
3070         C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3071         //C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3072         C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3073         C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3074         C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3075         C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3076         C *= par[0];// Norm
3077
3078
3079
3080         if(C < 0) continue;
3081
3082         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));
3083         else term1 = 0;
3084         term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1));
3085
3086         if(term1 > 0.0 && term2 > 0.0){
3087           lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2);
3088         }else if(term1==0 && term2 > 0.0){
3089           lnL += B_3[i][j][k]*log(term2);
3090         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
3091
3092       }
3093     }
3094   }
3095
3096   f = -2.0*lnL;
3097
3098 }
3099 //________________________________________________________________________
3100 double D3fitfunction_3(Double_t *x, Double_t *par)
3101 {
3102   double K12=CoulCorr2(+1, x[0]);
3103   double K13=CoulCorr2(+1, x[1]);
3104   double K23=CoulCorr2(+1, x[2]);
3105   K12=1.0; K13=1.0; K23=1.0;
3106   double K123=K12*K13*K23;
3107   //
3108   //K123 = 1 - pow(par[1],3)*(1 - K123);
3109   //K12 = 1 - pow(par[1],3)*(1-K12);
3110   //K13 = 1 - pow(par[1],3)*(1-K13);
3111   //K23 = 1 - pow(par[1],3)*(1-K23);
3112   //
3113   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);
3114   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);
3115   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);
3116   //EW12=1; EW13=1; EW23=1;
3117   //
3118   double G1=0, G2=0, G3=0, G4=0;
3119
3120   G1 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2)/2.)*EW12;
3121   G1 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2)/2.)*EW13;
3122   G1 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2)/2.)*EW23;
3123   //
3124   G2 = K12*exp(-pow(par[2]*x[0]/FmToGeV,2))*pow(EW12,2);
3125   G2 += K13*exp(-pow(par[2]*x[1]/FmToGeV,2))*pow(EW13,2);
3126   G2 += K23*exp(-pow(par[2]*x[2]/FmToGeV,2))*pow(EW23,2);
3127   //
3128   G3 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[2]*x[2])/2.)*EW12*EW23;
3129   G3 += exp(-pow(par[2]/FmToGeV,2)*(x[1]*x[1] + x[2]*x[2])/2.)*EW13*EW23;
3130   G3 += exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1])/2.)*EW12*EW13;
3131   G3 *= K123;
3132   //
3133   G4 = exp(-pow(par[2]/FmToGeV,2)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/2.);
3134   G4 *= EW12*EW13*EW23;
3135   G4 *= K123;
3136   //
3137   double C = 1 - pow(par[1],3)*(1-K123) - pow(par[1],2)*((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3138   //double C = 1 - (1-K123) - ((1-K12) + (1-K13) + (1-K23));// pure Coulomb
3139   C += pow(par[1],2)*2*par[3]*(1-par[3])*G1;// 2-mixed chaos/coherence
3140   C += pow(par[1],2)*pow(par[3],2)*G2;// 2-pure chaos
3141   C += pow(par[1],3)*2*pow(par[3],2)*(1-par[3])*G3;// 3-mixed chaos/coherence
3142   C += pow(par[1],3)*2*pow(par[3],3)*G4;// 3-pure chaos
3143   C *= par[0];// Norm
3144
3145   return C;
3146 }
3147
3148 void ReadCoulCorrections_Omega0(){
3149   // read in 3d 3-particle coulomb correlations = K3
3150
3151   TFile *coulfile;
3152   if(FileSetting!=6 && SourceType==0) coulfile = new TFile("KFile.root","READ");
3153   else coulfile = new TFile("KFile_Gauss.root","READ");
3154
3155   TString *name=new TString("K3ss_");
3156   *name += bBin-1;
3157   CoulCorrOmega0SS = (TH3D*)coulfile->Get(name->Data());
3158   CoulCorrOmega0SS->SetDirectory(0);
3159   name=new TString("K3os_");
3160   *name += bBin-1;
3161   CoulCorrOmega0OS = (TH3D*)coulfile->Get(name->Data());
3162   CoulCorrOmega0OS->SetDirectory(0);
3163   coulfile->Close();
3164
3165 }
3166 double CoulCorrGRS(bool SC, double Q_12, double Q_13, double Q_23){
3167   /*int Q12bin = CoulCorr2SS->GetXaxis()->FindBin(Q_12);
3168   int Q13bin = CoulCorr2SS->GetXaxis()->FindBin(Q_13);
3169   int Q23bin = CoulCorr2SS->GetXaxis()->FindBin(Q_23);*/
3170   int index12L = int(fabs(Q_12 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3171   int index12H = index12L+1;
3172   int index13L = int(fabs(Q_13 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3173   int index13H = index13L+1;
3174   int index23L = int(fabs(Q_23 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
3175   int index23H = index23L+1;
3176
3177   if(Tricubic){
3178     // Tricubic Interpolation
3179     double arr[4][4][4]={{{0}}};
3180     for(int x=0; x<4; x++){
3181       for(int y=0; y<4; y++){
3182         for(int z=0; z<4; z++){
3183           if(SC){
3184             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2SS->GetBinContent((index23L)+y)*CoulCorr2SS->GetBinContent((index13L)+z);
3185           }else{
3186             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2OS->GetBinContent((index23L)+y)*CoulCorr2OS->GetBinContent((index13L)+z);
3187           }
3188
3189         }
3190       }
3191     }
3192     return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3193   }else{
3194     // Trilinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3195     //
3196     double xd = (Q_12-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3197     xd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
3198     double yd = (Q_13-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3199     yd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index13H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
3200     double zd = (Q_23-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3201     zd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index23H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
3202     //
3203     if(SC){
3204       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3205       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3206       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
3207       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
3208       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3209       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3210       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
3211       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
3212       //
3213       double c0 = c00*(1-yd) + c10*yd;
3214       double c1 = c01*(1-yd) + c11*yd;
3215       return (c0*(1-zd) + c1*zd);
3216     }else{
3217       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3218       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3219       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
3220       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
3221       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3222       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3223       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
3224       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
3225       //
3226       double c0 = c00*(1-yd) + c10*yd;
3227       double c1 = c01*(1-yd) + c11*yd;
3228       return (c0*(1-zd) + c1*zd);
3229     }
3230   }
3231
3232 }
3233 double CoulCorrOmega0(bool SC, double Q_12, double Q_13, double Q_23){
3234   int Q12bin = CoulCorrOmega0SS->GetXaxis()->FindBin(Q_12);
3235   int Q13bin = CoulCorrOmega0SS->GetZaxis()->FindBin(Q_13);
3236   int Q23bin = CoulCorrOmega0SS->GetYaxis()->FindBin(Q_23);
3237   int index12L = int(fabs(Q_12 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3238   int index12H = index12L+1;
3239   int index13L = int(fabs(Q_13 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3240   int index13H = index13L+1;
3241   int index23L = int(fabs(Q_23 - CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorrOmega0SS->GetXaxis()->GetBinWidth(1)));
3242   int index23H = index23L+1;
3243
3244   if(SC){
3245     if(CoulCorrOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
3246
3247     if(Tricubic){
3248       //////////////////////////
3249       // Tricubic Interpolation
3250       double arr[4][4][4]={{{0}}};
3251       for(int x=0; x<4; x++){
3252         for(int y=0; y<4; y++){
3253           for(int z=0; z<4; z++){
3254             arr[x][y][z] = CoulCorrOmega0SS->GetBinContent((index12L)+x, (index23L)+y, (index13L)+z);
3255           }
3256         }
3257       }
3258       return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3259     }else{
3260       ///////////////////////////
3261       // Trilinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3262       //
3263       double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3264       xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3265       double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3266       yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3267       double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3268       zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3269
3270       // interpolate along x
3271       double c00=0, c10=0, c01=0, c11=0;
3272       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;
3273       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;
3274       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;
3275       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;
3276       //
3277       double c0=0, c1=0;
3278       if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
3279       if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
3280
3281       if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
3282       else {
3283         cout<<"Not all Omega0 bins well defined.  Use GRS instead"<<endl;
3284         return CoulCorrGRS(kTRUE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
3285       }
3286     }
3287
3288   }else {// mixed charge. Q12bin always designated as the same-charge pair
3289     if(CoulCorrOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) {cout<<"Problematic Omega0 bin"<<endl;}
3290
3291     if(Tricubic){
3292       //////////////////////////
3293       // Tricubic Interpolation
3294       double arr[4][4][4]={{{0}}};
3295       for(int x=0; x<4; x++){
3296         for(int y=0; y<4; y++){
3297           for(int z=0; z<4; z++){
3298             arr[x][y][z] = CoulCorrOmega0OS->GetBinContent((Q12bin)+x, (Q23bin)+y, (Q13bin)+z);
3299           }
3300         }
3301       }
3302       return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
3303     }else{
3304       ///////////////////////////
3305       // TriLinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
3306       //
3307       double xd = (Q_12-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3308       xd /= (CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorrOmega0SS->GetXaxis()->GetBinCenter(index12L+1));
3309       double yd = (Q_23-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3310       yd /= (CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23H+1)-CoulCorrOmega0SS->GetYaxis()->GetBinCenter(index23L+1));
3311       double zd = (Q_13-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3312       zd /= (CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13H+1)-CoulCorrOmega0SS->GetZaxis()->GetBinCenter(index13L+1));
3313       // interpolate along x
3314       double c00=0, c10=0, c01=0, c11=0;
3315       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;
3316       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;
3317       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;
3318       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;
3319       //
3320       double c0=0, c1=0;
3321       if(c00>0 && c10>0) c0 = c00*(1-yd) + c10*yd;
3322       if(c01>0 && c11>0) c1 = c01*(1-yd) + c11*yd;
3323
3324       if(c0>0 && c1>0) return (c0*(1-zd) + c1*zd);
3325       else {
3326         cout<<"Not all Omega0 bins well defined.  Use GRS instead"<<endl;
3327         return CoulCorrGRS(kFALSE, Q_12, Q_13, Q_23);// if cell not well defined use GRS
3328       }
3329     }
3330   }
3331
3332 }
3333
3334 double Gamov(int chargeProduct, double qinv){
3335
3336   double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
3337
3338   return arg/(exp(arg)-1);
3339 }
3340 /*double LednickyCorr(double qinv){
3341
3342   int indexL = int(fabs(1000*qinv-Lednicky_qinv[0])/(2*2.027027));// starts from 2.027 MeV
3343   int indexH = indexL+1;
3344   if( (indexL >= 73)) return 1.0;
3345   double slope = (Lednicky_CoulStrong[indexL] - Lednicky_CoulStrong[indexH])/(Lednicky_qinv[indexL]-Lednicky_qinv[indexH]);
3346   return (slope*(1000*qinv - Lednicky_qinv[indexL]) + Lednicky_CoulStrong[indexL]);
3347
3348 }
3349 void ReadLednickyFile(int Rvalue){
3350
3351   TString *filename=new TString("../Strong_FSI/converted_CoulStrong_");
3352   *filename += Rvalue;
3353   filename->Append("fm.root");
3354   TFile *Ledfile = new TFile(filename->Data(),"READ");
3355   TH1F *LedHisto = (TH1F*)Ledfile->Get("h27");
3356   for(int i=0; i<74; i++){
3357     Lednicky_qinv[i] = 2.*(LedHisto->GetXaxis()->GetBinLowEdge(i+1) + LedHisto->GetXaxis()->GetBinUpEdge(i+1))/2.;
3358     Lednicky_CoulStrong[i] = LedHisto->GetBinContent(i+1);
3359   }
3360   Ledfile->Close();
3361   }*/
3362 void ReadMomResFile(int rvalue, double lambda){
3363
3364   for(int i=0; i<40; i++){
3365     MomRes_C2_pp[i] = 1.;
3366     MomRes_C2_mp[i] = 1.;
3367     MomRes_term1_pp[i] = 1.;
3368     MomRes_term2_pp[i] = 1.;
3369     MomRes_term1_mp[i] = 1.;
3370     MomRes_term2_mp[i] = 1.;
3371     //
3372     MomResDev_C2_pp[i] = 1.;
3373     MomResDev_C2_mp[i] = 1.;
3374     MomResDev_term1_pp[i] = 1.;
3375     MomResDev_term2_pp[i] = 1.;
3376   }
3377
3378
3379   TString *momresfilename = new TString("MomResFile");
3380   momresfilename->Append("_Offline_");
3381   if(FileSetting==3) momresfilename->Append("PID1p5");
3382   else momresfilename->Append("11h");
3383   momresfilename->Append(".root");// no corresponding file for 10h
3384
3385   TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
3386
3387   TH2D *MomResWeights_pp = (TH2D*)MomResFile->Get("MomResHisto_pp");
3388   TH2D *MomResWeights_mp = (TH2D*)MomResFile->Get("MomResHisto_mp");
3389   TH2D *MomResWeights_pp_term1 = (TH2D*)MomResFile->Get("MomResHisto_pp_term1");
3390   TH2D *MomResWeights_pp_term2 = (TH2D*)MomResFile->Get("MomResHisto_pp_term2");
3391   TH2D *MomResWeights_mp_term1 = (TH2D*)MomResFile->Get("MomResHisto_mp_term1");
3392   TH2D *MomResWeights_mp_term2 = (TH2D*)MomResFile->Get("MomResHisto_mp_term2");
3393   //
3394   //
3395   TString *names3[2][5];// SC/MC, term#
3396   TString *names1D[2][5];// SC/MC, term#
3397   TString *names3_AvgK3[2];// SC/MC
3398   for(int ChProd=0; ChProd<2; ChProd++){
3399     if(ChProd==0) names3_AvgK3[ChProd] = new TString("AvgK3_SC_M");
3400     else names3_AvgK3[ChProd] = new TString("AvgK3_MC_M");
3401     *names3_AvgK3[ChProd] += MomResCentBin-1;
3402     //AvgK3[ChProd] = (TH3D*)MomResFile->Get(names3_AvgK3[ChProd]->Data());
3403     //AvgK3[ChProd]->SetDirectory(0);
3404     //
3405     for(int term=0; term<5; term++){
3406       //
3407       if(ChProd==0) {names3[ChProd][term] = new TString("MomResHisto3_SC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_SC_term");}
3408       else {names3[ChProd][term] = new TString("MomResHisto3_MC_term"); names1D[ChProd][term] = new TString("MomResHisto1D_MC_term");}
3409       *names3[ChProd][term] += term+1;
3410       *names1D[ChProd][term] += term+1;
3411       names3[ChProd][term]->Append("_M");
3412       names1D[ChProd][term]->Append("_M");
3413       *names3[ChProd][term] += MomResCentBin-1;
3414       *names1D[ChProd][term] += MomResCentBin-1;
3415       MomRes3d[ChProd][term] = (TH3D*)MomResFile->Get(names3[ChProd][term]->Data());
3416       MomRes1d[ChProd][term] = (TH1D*)MomResFile->Get(names1D[ChProd][term]->Data());
3417       MomRes3d[ChProd][term]->SetDirectory(0);
3418       MomRes1d[ChProd][term]->SetDirectory(0);
3419
3420     }
3421   }
3422
3423
3424   Int_t LambdaRbin;
3425   if(rvalue<=5) LambdaRbin = 1 + int(float(lambda-0.5+0.0001)/0.02);
3426   else LambdaRbin = 1 + 16*(rvalue-5) + int(float(lambda-0.5+0.0001)/0.02);
3427   Int_t LambdaRbinDev;
3428   if(rvalue<=5) LambdaRbinDev = 1 + 16*(1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3429   else LambdaRbinDev = 1 + 16*(rvalue-5 - 1) + int(float((lambda+0.04)-0.5+0.0001)/0.02);
3430
3431   for(int i=0; i<40; i++){
3432     MomRes_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbin, i+1);
3433     MomRes_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbin, i+1);
3434     MomRes_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbin, i+1);
3435     MomRes_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbin, i+1);
3436     MomRes_term1_mp[i] = MomResWeights_mp_term1->GetBinContent(LambdaRbin, i+1);
3437     MomRes_term2_mp[i] = MomResWeights_mp_term2->GetBinContent(LambdaRbin, i+1);
3438     //
3439     MomResDev_C2_pp[i] = MomResWeights_pp->GetBinContent(LambdaRbinDev, i+1);
3440     MomResDev_C2_mp[i] = MomResWeights_mp->GetBinContent(LambdaRbinDev, i+1);
3441     MomResDev_term1_pp[i] = MomResWeights_pp_term1->GetBinContent(LambdaRbinDev, i+1);
3442     MomResDev_term2_pp[i] = MomResWeights_pp_term2->GetBinContent(LambdaRbinDev, i+1);
3443   }
3444
3445   MomResFile->Close();
3446
3447 }
3448
3449 void fcn_r3_3d(int& npar, double* deriv, double& f, double par[], int flag){
3450
3451   double r3=0;
3452   double Chi2sum=0;
3453
3454   for(int i=fitlimitLowBin-1; i<fitlimitHighBin; i++){
3455     double q12=(i+0.5)*0.005;
3456     for(int j=fitlimitLowBin-1; j<fitlimitHighBin; j++){
3457       double q13=(j+0.5)*0.005;
3458       for(int k=fitlimitLowBin-1; k<fitlimitHighBin; k++){
3459         double q23=(k+0.5)*0.005;
3460
3461         //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);
3462         //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);
3463         r3 = par[0] - par[1]*(q12*q13 + q12*q23 + q13*q23) - par[2]*pow(q12*q13 + q12*q23 + q13*q23,2);
3464         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);
3465
3466       }
3467     }
3468   }
3469
3470   f = Chi2sum;
3471 }
3472
3473 double r3_fitfunction(Double_t *x, Double_t *par){
3474
3475   //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);
3476   //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);
3477   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);
3478
3479   return r3;
3480
3481 }
3482
3483 double cubicInterpolate (double p[4], double x) {
3484   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
3485 }
3486 double bicubicInterpolate (double p[4][4], double x, double y) {
3487         double arr[4];
3488         arr[0] = cubicInterpolate(p[0], y);
3489         arr[1] = cubicInterpolate(p[1], y);
3490         arr[2] = cubicInterpolate(p[2], y);
3491         arr[3] = cubicInterpolate(p[3], y);
3492         return cubicInterpolate(arr, x);
3493 }
3494 double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
3495         double arr[4];
3496         arr[0] = bicubicInterpolate(p[0], y, z);
3497         arr[1] = bicubicInterpolate(p[1], y, z);
3498         arr[2] = bicubicInterpolate(p[2], y, z);
3499         arr[3] = bicubicInterpolate(p[3], y, z);
3500         return cubicInterpolate(arr, x);
3501 }