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