PWGCF/FEMTOSCOPY/macros/Plot_PDCumulants_TPR.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <assert.h>
   7
   8 #include "TVector2.h"
   9 #include "TFile.h"
  10 #include "TString.h"
  11 #include "TF1.h"
  12 #include "TF3.h"
  13 #include "TH1.h"
  14 #include "TH2.h"
  15 #include "TH3.h"
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TMath.h"
  19 #include "TText.h"
  20 #include "TRandom3.h"
  21 #include "TArray.h"
  22 #include "TLegend.h"
  23 #include "TStyle.h"
  24 #include "TMinuit.h"
  25 #include "TCanvas.h"
  26 #include "TLatex.h"
  27 #include "TPaveStats.h"
  28 #include "TASImage.h"
  29 #include "TSpline.h"
  30
  31
  32 #define BohrR 1963.6885 // Bohr Radius for pions
  33 #define FmToGeV 0.19733 // conversion to fm
  34 #define PI 3.1415926
  35 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  36 double kappa3 = 0.1; // 0.16 (default), 0.05 or 0.3 as variations
  37 double kappa4 = 0.5; // 0.4 (default),
  38
  39 using namespace std;
  40
  41 int CollisionType_def=1;// 0=PbPb, 1=pPb, 2=pp
  42 bool MCcase_def=0;// MC data?
  43 int CHARGE_def=-1;// -1 or +1: + or - pions for same-charge case, --+ or ++- for mixed-charge case
  44 bool SameCharge_def=kTRUE;// 3-pion same-charge?
  45 bool AddCC=kTRUE;
  46 bool Gaussian=1;// Gaussian or Exponential?
  47 bool UseC2Bkg=1;// use Pythia/DPMJET bkg?
  48 //
  49 bool MuonCorrection=1;// apply Muon correction?
  50 bool IncludeExpansion=kTRUE;// Include EdgeWorth coefficients?
  51 bool FixExpansionAvg=1;
  52 int Mbin_def=18;// 0-19 (0=1050-2000 pions, 19=0-5 pions)
  53 int EDbin_def=0;// 0-2: Kt3 bin
  54 int Ktbin_def=1;// 1-6, 10=full range
  55 double MRCShift=1.0;// 1.0=full Momentum Resolution Correction. 1.1 for 10% systematic increase
  56 bool FitRangeShift=0;// 30% reduction in Fit Range
  57 //
  58 //
  59 bool SaveToFile_def=kFALSE;// Save outputs to file?
  60 int SourceType=0;// 0=Therminator, 1=Therminator with 50fm cut (keep at 0 for default), 2=Gaussian
  61 bool ConstantFSI=kFALSE;// Constant FSI's for each kt bin?
  62 bool GofP=kFALSE;// Include momentum dependence of coherent fraction?
  63 bool ChargeConstraint=kFALSE;// Include Charge Constraint for coherent states?
  64 bool LinkN=kFALSE;// Make N(++)=N(+-)?
  65 bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
  66 bool Tricubic=kFALSE;// Tricubic or Trilinear interpolation?  Tricubic was found to be unstable.
  67 bool IncludeSS=kTRUE;// Include same-charge two-pion correlations in the fit?
  68 bool IncludeOS=kFALSE;// Include mixed-charge two-pion correlations in the fit?
  69 //
  70 //
  71 double ExpPower=2;// default Gaussian(2), Exponential(1)
  72 //
  73 //
  74 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)
  75 const int Sbins_3=1;// 3-particle Species bins. 1=Pion-Pion-Pion only. max=10
  76 const int fitlimitLowBin = 3;
  77 const int fitlimitHighBin = 14;// max = 20 (100MeV)
  78 bool LHC11h=kTRUE;// always kTRUE
  79 bool ExplicitLoop=kFALSE;// always kFALSE
  80 bool PairCut=kTRUE;// always kTRUE
  81 int FSIindex;
  82 int Ktlowbin;
  83 int Kthighbin;
  84 int MomResCentBin;// momentum resolution centrality bin
  85 float TwoFrac;// Lambda parameter
  86 float OneFrac;// Lambda^{1/2}
  87 float ThreeFrac;// Lambda^{3/2}
  88 float Q2Limit;
  89 float Q3Limit;
  90
  91 double lambda_PM = 0.7;
  92
  93
  94 TH1D *CoulCorr2SS;
  95 TH1D *CoulCorr2OS;
  96 TH1D *CoulCorr2SS_2;// for interpolation in transition from Therminator to Gauss K factors
  97 TH1D *CoulCorr2OS_2;// for interpolation in transition from Therminator to Gauss K factors
  98 //
  99
 100 void ReadCoulCorrections(int);
 101 void ReadMomResFile(int);
 102 void ReadMuonCorrections(int,int);
 103 double CoulCorr2(int, double);
 104 double CoulCorrGRS(bool, double, double, double);
 105 double Dfitfunction_c3(Double_t *x, Double_t *par);
 106 double Gamov(int, double);
 107 double C2ssFitFunction(Double_t *x, Double_t *par);
 108 double C2osFitFunction(Double_t *x, Double_t *par);
 109 double cubicInterpolate(double[4], double);
 110 double bicubicInterpolate(double[4][4], double, double);
 111 double tricubicInterpolate(double[4][4][4], double, double, double);
 112 void DrawALICELogo(Bool_t, Float_t, Float_t, Float_t, Float_t);
 113
 114 //
 115 void fcnC2_Global(int&, double*, double&, double[], int);
 116
 117 const int BINRANGE_C2global=400;// 40
 118 double C2ssFitting[BINRANGE_C2global];
 119 double C2osFitting[BINRANGE_C2global];
 120 double C2ssFitting_e[BINRANGE_C2global];
 121 double C2osFitting_e[BINRANGE_C2global];
 122 double K2SS[BINRANGE_C2global];
 123 double K2OS[BINRANGE_C2global];
 124 double K2SS_e[BINRANGE_C2global];
 125 double K2OS_e[BINRANGE_C2global];
 126 double Chi2_C2global;
 127 double NFitPoints_C2global;
 128 double Chi2_C3global;
 129 double NFitPoints_C3global;
 130
 131 const int BINS_OSL=40;// out,side,long bins
 132 double A[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 133 double B[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 134 double A_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 135 double B_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 136 double avg_q[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 137 double K_OSL[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 138 double K_OSL_e[BINS_OSL][BINS_OSL][BINS_OSL];// out,side,long
 139 double Chi2_OSL;
 140 double NFitPoints_OSL;
 141
 142 const int BINRANGE_3=40;// q12,q13,q23
 143 int BINLIMIT_3=20;
 144 double A_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 145 double A_3_e[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 146 double B_3[BINRANGE_3][BINRANGE_3][BINRANGE_3];// q12,q13,q23
 147 double BinCenters[400];
 148 double BinWidthQ2;
 149 double Chi2_c3;
 150 double NFitPoints_c3;
 151 void fcn_c3(int&, double*, double&, double[], int);
 152
 153 double C3_base[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 154 double N_term1[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 155 double N_term2[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 156 double N_term3[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 157 double N_term4[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 158 double N_term5[BINRANGE_3][BINRANGE_3][BINRANGE_3];
 159
 160
 161 TH3D *MomRes3d[2][5];// SC/MC, term#
 162 TH1D *MomRes1d[2][5];// SC/MC, term#
 163 TH1D *MomResC2[2];// SC/MC
 164 int MRC2index;// index for C2 MRC
 165
 166 double AvgP[6][20];// kt bin, qinv bin
 167
 168 TH1D *C2muonCorrection_SC;
 169 TH1D *C2muonCorrection_MC;
 170 TH1D *C3muonCorrection;
 171
 172 TF1 *StrongSC;// same-charge pion strong FSI
 173 TF1 *MixedChargeSysFit;// mixed-charge 3-pion cumulant residue obtained from Plot_plotsTPR.C
 174 TF1 *BkgMCFit;// Pythia or DPMJET fit
 175 //
 176
 177
 178 void Plot_PDCumulants_TPR(bool SaveToFile=SaveToFile_def, int CollisionType=CollisionType_def, bool MCcase=MCcase_def, bool SameCharge=SameCharge_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
 179
 180   EDbin_def=EDbin;
 181   Ktbin_def=Ktbin;
 182   CollisionType_def=CollisionType;
 183   SaveToFile_def=SaveToFile;
 184   MCcase_def=MCcase;
 185   CHARGE_def=CHARGE;
 186   SameCharge_def=SameCharge;// 3-pion same-charge
 187   Mbin_def=Mbin;
 188   //
 189
 190   // PbPb then pPb then pp
 191   //double Nch_means[3][20]={{1828,1407,1071,896,757,619,496,402,324,252,171,117,83,63,50,36,   36,36,36,36},  {71,71,71,71,71,71,71,71,71,71,71,71,   71,57,45,32,23,16,9.7,5},   {47,47,47,47,47,47,47,47,47,47,47,47,47,   47,37,27,20,15,8.6,4.6}};
 192
 193   kappa3 = 0.1;
 194   kappa4 = 0.5;
 195   //kappa3 = 0.05 + 0.25*(Mbin/19.);// linear dependence of kappa3
 196
 197   //double MainFitParams[8]={0};
 198   //
 199   //double RefMainFitParams[8]={3.58014,  2.84635,  0.652247,  0.831958,  3.13505,  2.33717,  0.761147,  0.909211};
 200
 201   if(Gaussian) ExpPower=2.;
 202   else ExpPower=1.;
 203
 204   // old way v5 and before
 205   /*if(Mbin==0) FSIindex = 0;// 0-5% Therm
 206   else if(Mbin==1) FSIindex = 1;// 0-10% Therm
 207   else if(Mbin<=3) FSIindex = 2;// 10-20% Therm
 208   else if(Mbin<=5) FSIindex = 3;// 20-30% Therm
 209   else if(Mbin<=7) FSIindex = 4;// 30-40% Therm
 210   else if(Mbin<=9) FSIindex = 5;// 40-50% Therm
 211   else if(Mbin<=11) FSIindex = 6;// Gauss R=3.5 fm
 212   else if(Mbin<=14) FSIindex = 7;// Gauss R=3.0 fm
 213   else if(Mbin<=16) FSIindex = 8;// Gauss R=2.5 fm
 214   else FSIindex = 9;// Gauss R=2.0 fm*/
 215
 216   if(CollisionType==0){
 217     if(Mbin==0) FSIindex = 0;// 0-5% Therm
 218     else if(Mbin==1) FSIindex = 1;// 0-10% Therm
 219     else if(Mbin<=3) FSIindex = 2;// 10-20% Therm
 220     else if(Mbin<=5) FSIindex = 3;// 20-30% Therm
 221     else if(Mbin<=7) FSIindex = 4;// 30-40% Therm
 222     else if(Mbin<=9) FSIindex = 5;// 40-50% Therm
 223     else FSIindex = 6;// EW R=4.0 fm
 224     // next bin is 3.0 fm which is used for interpolation
 225   }else{
 226     if(Mbin<=14) FSIindex = 8;// EW R=2.5 fm
 227     else if(Mbin<=16) FSIindex = 9;// EW R=2.0 fm
 228     else FSIindex = 10;// EW R=1.7 fm
 229     // last bin is 1.25 fm which is used for interpolation
 230   }
 231
 232   // chooses lambda=0.7 for R=10,8,6,4,2
 233   if(Mbin<=2) MRC2index=138;
 234   else if(Mbin<=5) MRC2index=106;
 235   else if(Mbin<=9) MRC2index=74;
 236   else if(Mbin<=13) MRC2index=42;// was 9 by mistake, now 13 (corrected on 22/10/2013)
 237   else MRC2index=10;
 238
 239
 240
 241   TwoFrac = lambda_PM;
 242   OneFrac = sqrt(TwoFrac);
 243   ThreeFrac = pow(TwoFrac,3/2.);
 244
 245   if(CollisionType==0 && Mbin<10) BINLIMIT_3=20;
 246   else BINLIMIT_3=40;
 247
 248   Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
 249   Kthighbin=(Ktbin)*2+4;
 250   if(Ktbin==2) {cout<<"Kthighbin increased"<<endl; Kthighbin = 20;}// to make <KT3> comparible to <kT>
 251
 252
 253   BkgMCFit=new TF1("BkgMCFit","[0]*([1] + [2]*exp(-pow([3]*x,2)))",0,2.0);
 254   BkgMCFit->SetParameter(0,1);
 255   BkgMCFit->SetParameter(1,1.001);
 256   BkgMCFit->SetParameter(2,0.04102);
 257   BkgMCFit->SetParameter(3,1.874);
 258   BkgMCFit->SetLineColor(1);
 259   if(!MCcase){
 260     // kt bin(0.2<kt<0.3, 0.3<kt<1.0), Mult bin(low to high)
 261     float paramsPP_0[2][10]={{1.00014,1.00019,9.99733e-01,9.99379e-01,9.99488e-01, 9.99488e-01,9.99488e-01,9.99488e-01,9.99488e-01,9.99488e-01},{9.97267e-01,9.98918e-01,9.99117e-01,9.98878e-01,9.99457e-01, 9.99457e-01,9.99457e-01,9.99457e-01,9.99457e-01,9.99457e-01}};
 262     float paramsPP_1[2][10]={{1.00198,1.00156,1.00109,1.00068,1.00079, 1.00079,1.00079,1.00079,1.00079,1.00079},
 263                              {9.98091e-01,9.98893e-01,9.99782e-01,9.99595e-01,1.00037, 1.00037,1.00037,1.00037,1.00037,1.00037}};
 264     float paramsPP_2[2][10]={{4.51532e-02,4.30415e-02,4.10306e-02,3.42739e-02,2.95607e-02, 2.95607e-02,2.95607e-02,2.95607e-02,2.95607e-02,2.95607e-02},
 265                              {9.89882e-02,1.00657e-01,8.93560e-02,7.55225e-02,6.52689e-02, 6.52689e-02,6.52689e-02,6.52689e-02,6.52689e-02,6.52689e-02}};
 266     float paramsPP_3[2][10]={{1.57172,1.91807,1.87353,1.88549,2.22286, 2.22286,2.22286,2.22286,2.22286,2.22286},
 267                              {1.49834,1.65562,1.66501,1.75097,1.78597, 1.78597,1.78597,1.78597,1.78597,1.78597}};
 268     //////////////
 269     float paramsPPb_0[2][10]={{9.99479e-01,9.98772e-01,9.99295e-01,9.99532e-01,9.99430e-01,9.99642e-01,9.99424e-01, 9.99424e-01,9.99424e-01,9.99424e-01},
 270                               {9.91941e-01,9.96569e-01,9.97183e-01,9.98106e-01,9.98591e-01,9.98717e-01,9.98517e-01, 9.98517e-01,9.98517e-01,9.98517e-01}};
 271     float paramsPPb_1[2][10]={{1.00093,1.00023,1.00068,1.00082,1.00077,1.00087,1.00069, 1.00069,1.00069,1.00069},
 272                               {9.93727e-01,9.98118e-01,9.98630e-01,9.99584e-01,1.00005,1.00017,9.99972e-01, 9.99972e-01,9.99972e-01,9.99972e-01}};
 273     float paramsPPb_2[2][10]={{3.05502e-02,2.21034e-02,1.39907e-02,9.82480e-03,6.48672e-03,5.37259e-03,4.54498e-03, 4.54498e-03,4.54498e-03,4.54498e-03},
 274                               {8.86336e-02,5.53244e-02,3.59877e-02,2.57145e-02,1.82072e-02,1.44253e-02,1.28923e-02, 1.28923e-02,1.28923e-02, 1.28923e-02}};
 275     float paramsPPb_3[2][10]={{1.69217,1.51289,1.63134,1.71567,1.71918,1.78037,2.10736, 2.10736,2.10736,2.10736},
 276                               {1.21996,1.34570,1.35660,1.37801,1.39951,1.39093,1.33747, 1.33747,1.33747,1.33747}};
 277
 278     int ktindexBkg = Ktbin - 1;
 279     if(ktindexBkg > 1) ktindexBkg = 1;
 280     if(CollisionType==2){
 281       BkgMCFit->FixParameter(0, paramsPP_0[ktindexBkg][19-Mbin]);
 282       BkgMCFit->FixParameter(1, paramsPP_1[ktindexBkg][19-Mbin]);
 283       BkgMCFit->FixParameter(2, paramsPP_2[ktindexBkg][19-Mbin]);
 284       BkgMCFit->FixParameter(3, paramsPP_3[ktindexBkg][19-Mbin]);
 285     }else if(CollisionType==1){
 286       BkgMCFit->FixParameter(0, paramsPPb_0[ktindexBkg][19-Mbin]);
 287       BkgMCFit->FixParameter(1, paramsPPb_1[ktindexBkg][19-Mbin]);
 288       BkgMCFit->FixParameter(2, paramsPPb_2[ktindexBkg][19-Mbin]);
 289       BkgMCFit->FixParameter(3, paramsPPb_3[ktindexBkg][19-Mbin]);
 290     }else{// no Bkg for PbPb
 291       BkgMCFit->FixParameter(0, 1);
 292       BkgMCFit->FixParameter(1, 1);
 293       BkgMCFit->FixParameter(2, 0);
 294       BkgMCFit->FixParameter(3, 1);
 295     }
 296
 297   }
 298   // Core-Halo modeling of single-particle and triple-particle core fraction
 299   //float AvgN[10]={1.99266, 3.97789, 6.4624, 8.94042, 11.4194, 13.8987, 16.385, 18.8756, 21.3691, 26.0742};// 13c (avg total pion mult)/2.  2.0sigma
 300   //
 301
 302   // bin centers from QS+FSI
 303   double BinCenterPbPbCentral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509};
 304   double BinCenterPbPbPeripheral[40]={0.00206385, 0.00818515, 0.0129022, 0.0177584, 0.0226881, 0.027647, 0.032622, 0.0376015, 0.042588, 0.0475767, 0.0525692, 0.0575625, 0.0625569, 0.0675511, 0.0725471, 0.0775436, 0.0825399, 0.0875364, 0.0925339, 0.0975321, 0.102529, 0.107527, 0.112525, 0.117523, 0.122522, 0.12752, 0.132519, 0.137518, 0.142516, 0.147515, 0.152514, 0.157513, 0.162513, 0.167512, 0.172511, 0.177511, 0.18251, 0.187509, 0.192509, 0.197509};
 305   double BinCenterpPbAndpp[40]={0.00359275, 0.016357, 0.0257109, 0.035445, 0.045297, 0.0552251, 0.0651888, 0.0751397, 0.0851088, 0.0951108, 0.105084, 0.115079, 0.12507, 0.135059, 0.145053, 0.155049, 0.16505, 0.175038, 0.185039, 0.195034, 0.205023, 0.215027, 0.225024, 0.235023, 0.245011, 0.255017, 0.265017, 0.275021, 0.285021, 0.295017, 0.305018, 0.315018, 0.325013, 0.335011, 0.345016, 0.355019, 0.365012, 0.375016, 0.385017, 0.395016};
 306   if(CollisionType==0){
 307     Q3Limit = 0.1 + 0.1*Mbin/16.;
 308     for(int i=0; i<40; i++){
 309       if(Mbin<10) BinCenters[i] = BinCenterPbPbCentral[i];
 310       else BinCenters[i] = BinCenterPbPbPeripheral[i];
 311     }
 312     BinWidthQ2 = 0.005;
 313   }else{
 314     Q3Limit = 0.3 + 0.2*fabs(Mbin-10)/9.;
 315     for(int i=0; i<40; i++){
 316       BinCenters[i] = BinCenterpPbAndpp[i];
 317     }
 318     BinWidthQ2 = 0.01;
 319   }
 320   if(FitRangeShift) Q3Limit -= 0.3*Q3Limit;
 321   Q2Limit = Q3Limit/sqrt(2.);
 322   //Q2Limit = 1.2;
 323   //cout<<Q2Limit<<"  "<<Q3Limit<<endl;
 324
 325   // extend BinCenters for high q
 326   for(int index=40; index<400; index++){
 327     if(CollisionType==0) BinCenters[index] = (index+0.5)*(0.005);
 328     else BinCenters[index] = (index+0.5)*(0.010);
 329   }
 330   // Set 0's to 3-particle fit arrays
 331   for(int i=1; i<=BINLIMIT_3; i++){// bin number
 332     for(int j=1; j<=BINLIMIT_3; j++){// bin number
 333       for(int k=1; k<=BINLIMIT_3; k++){// bin number
 334         A_3[i-1][j-1][k-1]=0;
 335         A_3_e[i-1][j-1][k-1]=0;
 336         B_3[i-1][j-1][k-1]=0;
 337       }
 338     }
 339   }
 340
 341   // same-charge pion strong FSI (obtained from ratio of CoulStrong to Coul at R=7 Gaussian from Lednicky's code)
 342   StrongSC=new TF1("StrongSC","[0]+[1]*exp(-pow([2]*x,2))",0,1000);
 343   StrongSC->FixParameter(0,9.99192e-01);
 344   StrongSC->FixParameter(1,-8.01113e-03);
 345   StrongSC->FixParameter(2,5.35912e-02);
 346   // mixed-charge 3-pion cumulant residue obtained from Plot_plotsTPR.C
 347   MixedChargeSysFit=new TF1("MixedChargeSysFit","[0]+[1]*exp(-pow([2]*x/0.19733,2))",0,.5);
 348   MixedChargeSysFit->FixParameter(0,1);
 349   MixedChargeSysFit->FixParameter(1,.06);
 350   MixedChargeSysFit->FixParameter(2,5.5 - 6*(Mbin/19.));
 351
 352
 353   gStyle->SetOptStat(0);
 354   gStyle->SetOptDate(0);
 355   gStyle->SetOptFit(0111);
 356
 357   TFile *_file0;
 358   if(CollisionType==0){// PbPb
 359     if(MCcase) {
 360       if(Mbin<=1){
 361         _file0 = new TFile("Results/PDC_12a17a_R10_2.root","READ");
 362       }else{
 363         _file0 = new TFile("Results/PDC_12a17e_R10.root","READ");
 364       }
 365     }else{
 366       //_file0 = new TFile("Results/PDC_10h_11h_0to50_50to100.root","READ");
 367       //_file0 = new TFile("Results/PDC_10h_11h_0to50_50to100_3Ktbins.root","READ");
 368       //_file0 = new TFile("Results/PDC_10h_2Kt3bins.root","READ");
 369       //_file0 = new TFile("Results/PDC_10h_NewNorm.root","READ");
 370       //_file0 = new TFile("Results/PDC_10h_dEta0p03dPhi0p03.root","READ");
 371       //_file0 = new TFile("Results/PDC_10h_noPID.root","READ");
 372       //_file0 = new TFile("Results/PDC_10h_1percentCentral.root","READ");
 373       //_file0 = new TFile("Results/PDC_10h_11h_2Kt3bins.root","READ");// v5 and before
 374       //_file0 = new TFile("Results/PDC_10h_dEta0p025dPhi0p055.root","READ");
 375       //_file0 = new TFile("Results/PDC_11h_3Kt3bins_FB5and7overlap.root","READ");
 376       //_file0 = new TFile("Results/PDC_10h_11h_V0binning.root","READ");
 377       _file0 = new TFile("Results/PDC_10h_11h_NclsFix.root","READ");// standard
 378     }
 379   }else if(CollisionType==1){// pPb
 380     if(!MCcase){
 381       //_file0 = new TFile("Results/PDC_13c_kMB.root","READ");
 382       //_file0 = new TFile("Results/PDC_13bc.root","READ");
 383       //_file0 = new TFile("Results/PDC_13bc_3Ktbins.root","READ");
 384       //_file0 = new TFile("Results/PDC_13bc_kAnyINT_kINT7.root","READ");
 385       //_file0 = new TFile("Results/PDC_13bc_2Kt3bins.root","READ");
 386       //_file0 = new TFile("Results/PDC_13bc_kINT7.root","READ");// paper v1
 387       //_file0 = new TFile("Results/PDC_13c_posEta.root","READ");// eta sub-range
 388       //_file0 = new TFile("Results/PDC_13c_negEta.root","READ");// eta sub-range
 389       //_file0 = new TFile("Results/PDC_13c_0p4to0p8eta.root","READ");// eta sub-range
 390       //_file0 = new TFile("Results/PDC_13c_neg0p4toneg0p8eta.root","READ");// eta sub-range
 391       //_file0 = new TFile("Results/PDC_13bc_dEta0p025dPhi0p055_kINT7_kHighMult.root","READ");
 392       //_file0 = new TFile("Results/PDC_13bc_FB5and7overlap_V0binning.root","READ");
 393       //_file0 = new TFile("Results/PDC_13bc_kAnyINT_V0binning.root","READ");
 394       //_file0 = new TFile("Results/PDC_13bcd_kHighMult.root","READ");
 395       //_file0 = new TFile("Results/PDC_13bcd_kINT7.root","READ");
 396       //_file0 = new TFile("Results/PDC_13f_kHighMult.root","READ");
 397       //_file0 = new TFile("Results/PDC_13e_kHighMult.root","READ");
 398       //_file0 = new TFile("Results/PDC_13bcd_kINT7_plus_kHighMult.root","READ");// v5 and before
 399       //_file0 = new TFile("Results/PDC_13bcde_kINT7_NclsFix.root","READ");
 400       _file0 = new TFile("Results/PDC_13bcde_kINT7_kHighMult_NclsFix.root","READ");// standard
 401     }else{
 402       _file0 = new TFile("Results/PDC_13b2_efix_p1234_R2_2Ktbins.root","READ");// standard (gen level)
 403       //_file0 = new TFile("Results/PDC_13b2_efix_p1234_R2_RecLevel.root","READ");// Reconstructed tracks
 404       //_file0 = new TFile("Results/PDC_13b2_p1_noTTC.root","READ");
 405     }
 406   }else{// pp
 407     if(!MCcase){
 408       //_file0 = new TFile("Results/PDC_10cde.root","READ");
 409       //_file0 = new TFile("Results/PDC_10cde_3Ktbins.root","READ");
 410       //_file0 = new TFile("Results/PDC_10e_kHighMult.root","READ");
 411       //_file0 = new TFile("Results/PDC_10e_kAnyINT.root","READ");
 412       //_file0 = new TFile("Results/PDC_10cde_kAnyINT.root","READ");
 413       //_file0 = new TFile("Results/PDC_10cde_kAnyINT_Plus_kHighMult.root","READ");
 414       //_file0 = new TFile("Results/K0sRange_10cde.root","READ");// for K0s peak removal
 415       //_file0 = new TFile("Results/PDC_10cde_2Kt3bins.root","READ");
 416       //_file0 = new TFile("Results/PDC_10e_kHighMult.root","READ");
 417       //_file0 = new TFile("Results/PDC_10cde_kMB.root","READ");
 418       //_file0 = new TFile("Results/PDC_10cde_kMB_plus_kHighMult_2Kt3bins.root","READ");// v2 paper
 419       //_file0 = new TFile("Results/PDC_10cde_kMB_plus_kHighMult_AOD147.root","READ");// v5 and before
 420       //_file0 = new TFile("Results/PDC_10cde_FB5and7overlap.root","READ");
 421       //_file0 = new TFile("Results/PDC_10d_noTTC.root","READ");
 422       //_file0 = new TFile("Results/PDC_10cde_kMB.root","READ");
 423       _file0 = new TFile("Results/PDC_10cde_kMB_kHighMult_NclsFix.root","READ");// standard
 424     }else{
 425       _file0 = new TFile("Results/PDC_10f6a_R2_2Ktbins.root","READ");// standard
 426       //_file0 = new TFile("Results/PDC_10f6a_R10_genLevel_2Ktbins.root","READ");
 427     }
 428   }
 429
 430   ReadCoulCorrections(FSIindex);
 431   ReadMomResFile(Mbin);
 432   ReadMuonCorrections(CollisionType,Mbin);
 433   //
 434   /////////////////////////////////////////////////////////
 435
 436
 437   double NormQcutLow;
 438   double NormQcutHigh;
 439   if(CollisionType==0) {// PbPb
 440     if(Mbin<10){
 441       NormQcutLow = 0.15;// was 0.15
 442       NormQcutHigh = 0.175;// was 0.175
 443     }else{
 444       NormQcutLow = 0.3;// was 0.3
 445       NormQcutHigh = 0.35;// was 0.35
 446     }
 447   }else{// pPb and pp
 448     NormQcutLow = 1.0;// was 1.0
 449     NormQcutHigh = 1.2;// was 1.2
 450   }
 451
 452
 453   TList *MyList;
 454   TDirectoryFile *tdir;
 455   if(!MCcase) tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputThreePionRadiiAnalysis.root");
 456
 457
 458   if(CollisionType==0){
 459     if(!MCcase){
 460       if(Mbin<6) MyList=(TList*)tdir->Get("ThreePionRadiiOutput_1");
 461       else MyList=(TList*)tdir->Get("ThreePionRadiiOutput_2");
 462     }else{
 463       MyList=(TList*)_file0->Get("MyList");
 464     }
 465   }else {
 466     if(!MCcase) {
 467       if(CollisionType==2 && Mbin<15) MyList=(TList*)tdir->Get("ThreePionRadiiOutput_2");
 468       else MyList=(TList*)tdir->Get("ThreePionRadiiOutput_1");
 469     }else MyList=(TList*)_file0->Get("MyList");
 470   }
 471   _file0->Close();
 472
 473   TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
 474   //
 475   cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
 476
 477   // below lines determine fractional cross-sections without correction for efficiency
 478   //TH1F *fMultDist3 = (TH1F*)MyList->FindObject("fMultDist3");// total event count before track cuts
 479   //for(int mm=19; mm>=10; mm--){
 480   //cout<<"fracton of cross-section = "<<Events->GetBinContent(mm+1)/fMultDist3->GetEntries()<<endl;
 481   //}
 482
 483   TF1 *Unity = new TF1("Unity","1",0,10);
 484   Unity->SetLineStyle(2); Unity->SetLineColor(1);
 485
 486   // Explicit Loop Histos
 487   TH2D *Two_ex_2d[2][2][Sbins_2][2];
 488   TH1D *Two_ex[2][2][Sbins_2][2];
 489
 490   // Normalizations
 491   double NormH_pc[5]={0};
 492
 493   double norm_ex_2[6][2]={{0}};
 494
 495   // 3d method histograms
 496   TH3D *Three_3d[2][2][2][Sbins_3][5];
 497   TH1D *Three_1d[2][2][2][Sbins_3][5];
 498   ///////////////////////////////////
 499   // Create Histograms
 500
 501   // 2-particle
 502   for(int c1=0; c1<2; c1++){
 503     for(int c2=0; c2<2; c2++){
 504       for(int sc=0; sc<Sbins_2; sc++){
 505         for(int term=0; term<2; term++){
 506
 507           TString *name2_ex = new TString("Explicit2_Charge1_");
 508           *name2_ex += c1;
 509           name2_ex->Append("_Charge2_");
 510           *name2_ex += c2;
 511           name2_ex->Append("_SC_");
 512           *name2_ex += sc;
 513           name2_ex->Append("_M_");
 514           *name2_ex += Mbin;
 515           name2_ex->Append("_ED_");
 516           *name2_ex += 0;// EDbin
 517           name2_ex->Append("_Term_");
 518           *name2_ex += term+1;
 519
 520           if(sc==0 || sc==3 || sc==5){
 521             if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
 522           }
 523
 524           Two_ex_2d[c1][c2][sc][term] = (TH2D*)MyList->FindObject(name2_ex->Data());
 525           Two_ex_2d[c1][c2][sc][term]->Sumw2();
 526           TString *proname = new TString(name2_ex->Data());
 527           proname->Append("_pro");
 528
 529           if(Ktbin==10) {Ktlowbin=1; Kthighbin=Two_ex_2d[c1][c2][sc][term]->GetNbinsX();}//full kt interval
 530           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)
 531
 532           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));
 533           Two_ex[c1][c2][sc][term]->Scale(norm_ex_2[sc][0]/norm_ex_2[sc][term]);
 534
 535           Two_ex[c1][c2][sc][term]->SetMarkerStyle(20);
 536           Two_ex[c1][c2][sc][term]->GetXaxis()->SetTitle("#font[12]{q}_{inv}");
 537           Two_ex[c1][c2][sc][term]->GetYaxis()->SetTitle("#font[12]{C}_{2}");
 538           Two_ex[c1][c2][sc][term]->SetTitle("");
 539
 540         }// term
 541       }// sc
 542
 543       // 3-particle
 544       for(int c3=0; c3<2; c3++){
 545         for(int sc=0; sc<Sbins_3; sc++){
 546           for(int term=0; term<5; term++){
 547
 548             TString *name3_ex = new TString("Explicit3_Charge1_");
 549             *name3_ex += c1;
 550             name3_ex->Append("_Charge2_");
 551             *name3_ex += c2;
 552             name3_ex->Append("_Charge3_");
 553             *name3_ex += c3;
 554             name3_ex->Append("_SC_");
 555             *name3_ex += sc;
 556             name3_ex->Append("_M_");
 557             *name3_ex += Mbin;
 558             name3_ex->Append("_ED_");
 559             *name3_ex += EDbin;
 560             name3_ex->Append("_Term_");
 561             *name3_ex += term+1;
 562
 563
 564             TString *name3_pc = new TString("PairCut3_Charge1_");
 565             *name3_pc += c1;
 566             name3_pc->Append("_Charge2_");
 567             *name3_pc += c2;
 568             name3_pc->Append("_Charge3_");
 569             *name3_pc += c3;
 570             name3_pc->Append("_SC_");
 571             *name3_pc += sc;
 572             name3_pc->Append("_M_");
 573             *name3_pc += Mbin;
 574             name3_pc->Append("_ED_");
 575             *name3_pc += EDbin;
 576             name3_pc->Append("_Term_");
 577             *name3_pc += term+1;
 578
 579             ///////////////////////////////////////
 580             // skip degenerate histograms
 581             if(sc==0 || sc==6 || sc==9){// Identical species
 582               if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
 583               if( (c1+c2+c3)==2) {if(c1!=0) continue;}
 584             }else if(sc!=5){
 585               if( (c1+c2)==1) {if(c1!=0) continue;}
 586             }else {}// do nothing for pi-k-p case
 587
 588             /////////////////////////////////////////
 589
 590
 591             if(PairCut){
 592
 593               TString *nameNorm=new TString("PairCut3_Charge1_");
 594               *nameNorm += c1;
 595               nameNorm->Append("_Charge2_");
 596               *nameNorm += c2;
 597               nameNorm->Append("_Charge3_");
 598               *nameNorm += c3;
 599               nameNorm->Append("_SC_");
 600               *nameNorm += sc;
 601               nameNorm->Append("_M_");
 602               *nameNorm += Mbin;
 603               nameNorm->Append("_ED_");
 604               *nameNorm += 0;// EDbin
 605               nameNorm->Append("_Term_");
 606               *nameNorm += term+1;
 607               nameNorm->Append("_Norm");
 608               NormH_pc[term] = ((TH1D*)MyList->FindObject(nameNorm->Data()))->Integral();
 609
 610               if(NormH_pc[term] > 0){
 611
 612                 if(sc<=2) {
 613
 614                   TString *name_3d = new TString(name3_pc->Data());
 615                   name_3d->Append("_3D");
 616                   Three_3d[c1][c2][c3][sc][term] = (TH3D*)MyList->FindObject(name_3d->Data());
 617                   Three_3d[c1][c2][c3][sc][term]->Sumw2();
 618                   Three_3d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]);
 619                   TString *name_Q3 = new TString(name3_pc->Data());
 620                   name_Q3->Append("_Q3");
 621                   Three_1d[c1][c2][c3][sc][term] = (TH1D*)MyList->FindObject(name_Q3->Data());
 622                   Three_1d[c1][c2][c3][sc][term]->Sumw2();
 623                   Three_1d[c1][c2][c3][sc][term]->Scale(NormH_pc[0]/NormH_pc[term]);
 624                   //cout<<"Scale factor = "<<NormH_pc[0]/NormH_pc[term]<<endl;
 625
 626
 627                 }
 628
 629               }else{
 630                 cout<<"normalization = 0.  Skipping this SC.  SC="<<sc<<"  c1="<<c1<<"  c2="<<c2<<"  c3="<<c3<<endl;
 631               }
 632
 633
 634
 635             }// pair cut
 636           }// term
 637
 638         }// SC
 639
 640
 641       }// c3
 642
 643     }// c2
 644   }// c1
 645
 646
 647   cout<<"Done getting Histograms"<<endl;
 648
 649
 650
 651   TCanvas *can1 = new TCanvas("can1", "can1",11,53,800,800);
 652   can1->SetHighLightColor(2);
 653   //can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 654   gStyle->SetOptFit(0111);
 655   can1->SetFillColor(10);
 656   can1->SetBorderMode(0);
 657   can1->SetBorderSize(2);
 658   can1->SetFrameFillColor(0);
 659   can1->SetFrameBorderMode(0);
 660   can1->SetFrameBorderMode(0);
 661   gPad->SetRightMargin(0.025); gPad->SetLeftMargin(0.1); gPad->SetTopMargin(0.02);
 662   gPad->SetGridx(0);
 663   gPad->SetGridy(0);
 664   TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
 665   legend1->SetBorderSize(1);
 666   legend1->SetTextSize(.04);// small .03; large .036
 667   legend1->SetFillColor(0);
 668
 669   TLatex *Specif_System = new TLatex(0.23, 0.9,"ALICE p-Pb #sqrt{#font[12]{s}_{NN}}=5.02 TeV, #LT#font[12]{N}_{ch}#GT = 9.8 #pm 0.5");
 670   TLatex *Specif_KT3 = new TLatex(0.23, 0.8,"0.16<#font[12]{K}_{T,3}<0.3 GeV/#font[12]{c}");
 671   TLatex *Specif_kT = new TLatex(0.23, 0.8,"0.2<#font[12]{k}_{T}<0.3 GeV/#font[12]{c}");
 672   TLatex *Specif_FSI = new TLatex(0.23, 0.7,"FSI uncorrected");
 673   TLatex *Specif_Disclaimer = new TLatex(0.23, 0.6,"Statistical errors only");
 674   Specif_FSI->SetTextFont(42); Specif_System->SetTextFont(42); Specif_KT3->SetTextFont(42); Specif_Disclaimer->SetTextFont(42);
 675   Specif_kT->SetTextFont(42);
 676   Specif_FSI->SetTextSize(0.04); Specif_System->SetTextSize(0.04); Specif_KT3->SetTextSize(0.04); Specif_Disclaimer->SetTextSize(0.04);
 677   Specif_kT->SetTextSize(0.04);
 678   Specif_FSI->SetNDC(); Specif_System->SetNDC(); Specif_KT3->SetNDC(); Specif_Disclaimer->SetNDC();
 679   Specif_kT->SetNDC();
 680
 681   /////////////////////////
 682   // 2-particle legend
 683   // 0 = pi-pi
 684   // 1 = pi-k
 685   // 2 = pi-p
 686   // 3 = k-k
 687   // 4 = k-p
 688   // 5 = p-p
 689   /////////////////////////
 690   TH1D *Two_pipi_mp = (TH1D*)Two_ex[0][1][0][0]->Clone();
 691   Two_pipi_mp->Divide(Two_ex[0][1][0][1]);
 692
 693   // Normalization range counting.  Just to evaluate required normalization width in qinv.
 694   //cout<<Two_ex[0][0][0][0]->Integral(Two_ex[0][0][0][0]->GetXaxis()->FindBin(0), Two_ex[0][0][0][0]->GetXaxis()->FindBin(0.1))<<endl;
 695   //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;
 696   //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;
 697
 698
 699
 700   const int SCOI_2=0;
 701
 702   TH1D *Two_ex_clone_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 703   Two_ex_clone_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 704   TH1D *Two_ex_clone_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 705   Two_ex_clone_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 706   TH1D *Two_ex_clone_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 707   Two_ex_clone_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 708
 709
 710   Two_ex_clone_mm->GetYaxis()->SetTitle("#font[12]{C}_{2}");
 711   Two_ex_clone_mm->SetTitle("");
 712   Two_ex_clone_mp->GetYaxis()->SetTitle("#font[12]{C}_{2}");
 713   Two_ex_clone_mm->SetMarkerColor(2);
 714   Two_ex_clone_mm->SetLineColor(2);
 715   Two_ex_clone_mp->SetMarkerColor(1);
 716   Two_ex_clone_pp->SetMarkerColor(4);
 717   Two_ex_clone_pp->SetLineColor(4);
 718   Two_ex_clone_mm->GetXaxis()->SetRangeUser(0,0.6);
 719   Two_ex_clone_mm->SetMinimum(0.95);
 720   Two_ex_clone_mm->SetMaximum(2.4);
 721
 722
 723   if(MCcase){
 724     Two_ex_clone_mp->SetMarkerColor(4);
 725     Two_ex_clone_mp->SetLineColor(4);
 726     Two_ex_clone_mm->Add(Two_ex_clone_pp);
 727     Two_ex_clone_mm->Scale(0.5);
 728     Two_ex_clone_mm->GetYaxis()->SetTitle("#font[12]{C}_{2}^{MC}");
 729     Two_ex_clone_mm->GetYaxis()->SetTitleOffset(1.3);
 730     Two_ex_clone_mm->SetMinimum(0.95);
 731     Two_ex_clone_mm->SetMaximum(1.3);
 732     Two_ex_clone_mm->DrawCopy();
 733     gStyle->SetOptFit(1111);
 734     Two_ex_clone_mm->Fit(BkgMCFit,"","",0,1.8);
 735     //BkgMCFit->Draw("same");
 736     Two_ex_clone_mp->DrawCopy("same");
 737     legend1->AddEntry(Two_ex_clone_mm,"same-charge","p");
 738     //legend1->AddEntry(Two_ex_clone_pp,"++","p");
 739     legend1->AddEntry(Two_ex_clone_mp,"mixed-charge","p");
 740     legend1->Draw("same");
 741   }
 742
 743   //for(int i=0; i<Two_ex_clone_mm->GetNbinsX(); i++){
 744   //cout<<Two_ex_clone_mm->GetBinError(i+1)<<", ";
 745   //}
 746
 747   /////////////////////////////////////////////////////
 748   // Global fitting C2os and C2ss
 749   double C2ssSys_e[BINRANGE_C2global]={0};
 750   double C2osSys_e[BINRANGE_C2global]={0};
 751   //
 752   const int npar=10;// 10
 753   TMinuit MyMinuit(npar);
 754   MyMinuit.SetFCN(fcnC2_Global);
 755   double arglist_C2 = 0; int dummy=0;
 756   MyMinuit.mnexcm("SET NOWarnings",&arglist_C2,1, dummy);
 757   arglist_C2 = -1;
 758   MyMinuit.mnexcm("SET PRint",&arglist_C2,1, dummy);
 759   double OutputPar[npar]={0};
 760   double OutputPar_e[npar]={0};
 761
 762   double par[npar];               // the start values
 763   double stepSize[npar];          // step sizes
 764   double minVal[npar];            // minimum bound on parameter
 765   double maxVal[npar];            // maximum bound on parameter
 766   string parName[npar];
 767
 768   TH1D *temp_mm=(TH1D*)Two_ex[0][0][SCOI_2][0]->Clone();
 769   temp_mm->Divide(Two_ex[0][0][SCOI_2][1]);
 770   TH1D *temp_mp=(TH1D*)Two_ex[0][1][SCOI_2][0]->Clone();
 771   temp_mp->Divide(Two_ex[0][1][SCOI_2][1]);
 772   TH1D *temp_pp=(TH1D*)Two_ex[1][1][SCOI_2][0]->Clone();
 773   temp_pp->Divide(Two_ex[1][1][SCOI_2][1]);
 774   TH1D *C2ssRaw=(TH1D*)temp_mm->Clone();// MRC uncorrected
 775   TH1D *C2osRaw=(TH1D*)temp_mp->Clone();// MRC uncorrected
 776   C2ssRaw->SetMarkerStyle(24);
 777   C2osRaw->SetMarkerStyle(21);//21
 778   C2ssRaw->SetMarkerColor(2);
 779   C2osRaw->SetMarkerColor(4);
 780
 781   for(int i=0; i<BINRANGE_C2global; i++){
 782     if(i >= temp_mm->GetNbinsX()) continue;// bin limit
 783
 784     C2ssFitting[i] = (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.;
 785     double MRC_SC = 1.0;
 786     double MRC_MC = 1.0;
 787     if(MomResC2[0]->GetBinContent(i+1) > 0) MRC_SC = (MomResC2[0]->GetBinContent(i+1)-1)*MRCShift + 1;
 788     if(MomResC2[1]->GetBinContent(i+1) > 0) MRC_MC = (MomResC2[1]->GetBinContent(i+1)-1)*MRCShift + 1;
 789     if(!GeneratedSignal && !MCcase) C2ssFitting[i] *= MRC_SC;
 790     C2ssFitting_e[i] = pow(MRC_SC*sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2);
 791     C2ssRaw->SetBinContent(i+1, (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.);
 792     C2ssRaw->SetBinError(i+1, pow(sqrt(pow(temp_mm->GetBinError(i+1),2) + pow(temp_pp->GetBinError(i+1),2))/sqrt(2.),2));
 793     C2ssFitting_e[i] = sqrt(C2ssFitting_e[i]);
 794     C2osFitting[i] = temp_mp->GetBinContent(i+1);
 795     if(!GeneratedSignal && !MCcase) C2osFitting[i] *= MRC_MC;
 796     C2osFitting_e[i] = pow(MRC_MC*temp_mp->GetBinError(i+1),2);
 797     C2osRaw->SetBinContent(i+1, temp_mp->GetBinContent(i+1));
 798     C2osRaw->SetBinError(i+1, temp_mp->GetBinError(i+1));
 799     C2osFitting_e[i] += pow((MRC_MC-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 800     C2osFitting_e[i] = sqrt(C2osFitting_e[i]);
 801     //
 802     C2ssSys_e[i] = pow((MRC_SC-1)*0.1 * (temp_mm->GetBinContent(i+1) + temp_pp->GetBinContent(i+1))/2.,2);
 803     C2ssSys_e[i] = sqrt(C2ssSys_e[i]);
 804     C2osSys_e[i] = pow((MRC_MC-1)*0.1 * temp_mp->GetBinContent(i+1),2);
 805     C2osSys_e[i] = sqrt(C2osSys_e[i]);
 806     //
 807     if(UseC2Bkg) C2ssFitting[i] /= BkgMCFit->Eval(BinCenters[i]);
 808     //
 809     // method with undilution for plotting purposes
 810     //C2ssFitting[i] = (C2ssFitting[i] - (1-lambda_PM)) / (CoulCorr2(+1, BinCenters[i]) * lambda_PM);
 811     //
 812     K2SS[i] = CoulCorr2(+1, BinCenters[i]);
 813     K2OS[i] = CoulCorr2(-1, BinCenters[i]);
 814     //K2SS[i] = 1;
 815     //K2OS[i] = 1;
 816     //
 817     K2SS_e[i] = 0.0;
 818     K2OS_e[i] = 0.0;
 819   }
 820
 821
 822
 823
 824
 825   C2ssFitting[0]=-1000; C2osFitting[0]=-1000;
 826   C2ssFitting_e[0]=1000; C2osFitting_e[0]=1000;
 827   K2SS_e[0]=1000; K2OS_e[0]=1000;
 828   int ierflg=0;
 829   double arglist[10];
 830   //arglist[0]=2;// improve Minimization Strategy (1 is default)
 831   //MyMinuit.mnexcm("SET STR",arglist,1,ierflg);
 832   //arglist[0] = 0;
 833   //MyMinuit.mnexcm("SCAN", arglist,1,ierflg);
 834   arglist[0] = 5000;// 5000
 835   MyMinuit.mnexcm("MIGRAD", arglist ,1,ierflg);
 836
 837   TF1 *fitC2ss_Base = new TF1("fitC2ss_Base",C2ssFitFunction, 0.005,2, npar);//0.2
 838   TF1 *fitC2ss_Expan = new TF1("fitC2ss_Expan",C2ssFitFunction, 0.005,2, npar);//0.2
 839   fitC2ss_Base->SetLineStyle(2);
 840   TH1D *fitC2ss_h = new TH1D("fitC2ss_h","",Two_ex_clone_mm->GetNbinsX(),Two_ex_clone_mm->GetXaxis()->GetBinLowEdge(1), Two_ex_clone_mm->GetXaxis()->GetBinUpEdge(Two_ex_clone_mm->GetNbinsX()));
 841
 842   for(int ft=0; ft<2; ft++){// Gaussian or EW
 843     if(ft==1 && !IncludeExpansion) continue;
 844     par[0] = 1.0; par[1]=0.7; par[2]=0.5; par[3]=7.2; par[4] = .1; par[5] = .0; par[6] = .0; par[7] = 0.; par[8] = 0.; par[9] = 0.;// was par[1]=0.5
 845     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;
 846     minVal[0] = 0.955; minVal[1] = 0.5; minVal[2] = 0.; minVal[3] = 0.1; minVal[4] = 0.001; minVal[5] = -10.; minVal[6] = -10.; minVal[7] = -10.; minVal[8]=-10; minVal[9] = 0.995;// minVal[1] was 0.2
 847     maxVal[0] = 1.1; maxVal[1] = 1.3; maxVal[2] = 0.99; maxVal[3] = 15.; maxVal[4] = 2.; maxVal[5] = 10.; maxVal[6] = 10.; maxVal[7] = 10.; maxVal[8]=10.; maxVal[9]=1.1;// maxVal[1] was 1.0
 848     if(Gaussian==kFALSE) maxVal[1]=2.0;
 849
 850     parName[0] = "Norm"; parName[1] = "#Lambda"; parName[2] = "G"; parName[3] = "Rch"; parName[4] = "Rcoh";
 851     parName[5] = "coeff_3"; parName[6] = "coeff_4"; parName[7] = "coeff_5"; parName[8]="coeff_6"; parName[9]="Norm_2";
 852
 853     for (int i=0; i<npar; i++){
 854       MyMinuit.DefineParameter(i, parName[i].c_str(), par[i], stepSize[i], minVal[i], maxVal[i]);
 855     }
 856
 857     //MyMinuit.DefineParameter(1, parName[1].c_str(), 0.7, stepSize[1], minVal[1], maxVal[1]); MyMinuit.FixParameter(1);// lambda
 858     //MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
 859     MyMinuit.DefineParameter(2, parName[2].c_str(), 0., stepSize[2], minVal[2], maxVal[2]); MyMinuit.FixParameter(2);// G
 860     MyMinuit.DefineParameter(4, parName[4].c_str(), 0., stepSize[4], minVal[4], maxVal[4]); MyMinuit.FixParameter(4);// Rcoh
 861     MyMinuit.DefineParameter(7, parName[7].c_str(), 0, stepSize[7], minVal[7], maxVal[7]); MyMinuit.FixParameter(7);// k5
 862     MyMinuit.DefineParameter(8, parName[8].c_str(), 0, stepSize[8], minVal[8], maxVal[8]); MyMinuit.FixParameter(8);// k6
 863
 864     //
 865     if(ft==0){
 866       MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
 867       MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
 868     }else{// IncludeExpansion
 869       if(FixExpansionAvg){
 870         MyMinuit.DefineParameter(5, parName[5].c_str(), kappa3, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
 871         MyMinuit.DefineParameter(6, parName[6].c_str(), kappa4, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
 872
 873       }else{
 874         Int_t err=0;
 875         Double_t tmp[1];
 876         tmp[0] = 5+1;
 877         MyMinuit.mnexcm( "RELease", tmp,  1,  err );// kappa3
 878         tmp[0] = 6+1;
 879         MyMinuit.mnexcm( "RELease", tmp,  1,  err );// kappa4
 880         //
 881       }
 882     }
 883
 884     if(IncludeSS==kFALSE){
 885       MyMinuit.DefineParameter(3, parName[3].c_str(), 9.1, stepSize[3], minVal[3], maxVal[3]); MyMinuit.FixParameter(3);// Rch
 886       MyMinuit.DefineParameter(0, parName[0].c_str(), .998, stepSize[0], minVal[0], maxVal[0]); MyMinuit.FixParameter(0);// N
 887       MyMinuit.DefineParameter(5, parName[5].c_str(), 0, stepSize[5], minVal[5], maxVal[5]); MyMinuit.FixParameter(5);// k3
 888       MyMinuit.DefineParameter(6, parName[6].c_str(), 0, stepSize[6], minVal[6], maxVal[6]); MyMinuit.FixParameter(6);// k4
 889     }
 890     if(IncludeOS==kFALSE){
 891       MyMinuit.DefineParameter(9, parName[9].c_str(), 1.002, stepSize[9], minVal[9], maxVal[9]); MyMinuit.FixParameter(9);// N_2
 892     }
 893
 894
 895     // Do the minimization!
 896     if(!MCcase){
 897       cout<<"Start C2 Global fit"<<endl;
 898       MyMinuit.Migrad();// generally the best minimization algorithm
 899       cout<<"End C2 Global fit"<<endl;
 900     }
 901     for (int i=0; i<npar; i++){
 902       MyMinuit.GetParameter(i,OutputPar[i],OutputPar_e[i]);
 903     }
 904     MyMinuit.mnexcm("SHOw PARameters", &arglist_C2, 1, ierflg);
 905     cout<<"C2 fit: Chi2/NDF = "<<Chi2_C2global/(NFitPoints_C2global - MyMinuit.GetNumFreePars())<<"   Chi^2="<<Chi2_C2global<<endl;
 906
 907     if(ft==0){// Gaussian or Exponential
 908       for(int i=0; i<npar; i++) {
 909         fitC2ss_Base->FixParameter(i,OutputPar[i]);
 910         fitC2ss_Base->SetParError(i,OutputPar_e[i]);
 911       }
 912
 913     }else{// Edgeworth or Laguerre
 914       for(int i=0; i<npar; i++) {
 915         fitC2ss_Expan->FixParameter(i,OutputPar[i]);
 916         fitC2ss_Expan->SetParError(i,OutputPar_e[i]);
 917       }
 918       for(int bin=1; bin<=Two_ex_clone_mm->GetNbinsX(); bin++){
 919         double qinv = Two_ex_clone_mm->GetXaxis()->GetBinCenter(bin);
 920         if(!MCcase) fitC2ss_h->SetBinContent(bin, fitC2ss_Expan->Eval(qinv));// uncomment to show expansion
 921       }
 922
 923       double lambdaStar=0, lambdaStar_e=0;
 924       if(Gaussian){
 925         lambdaStar = OutputPar[1]*pow(1 + OutputPar[6]/8.,2);
 926         lambdaStar_e = pow(OutputPar_e[1]*(1 + OutputPar[6]/8.),2);
 927         lambdaStar_e = sqrt(lambdaStar_e);
 928       }else{
 929         lambdaStar = OutputPar[1]*pow(1 - OutputPar[5] + OutputPar[6],2);
 930         lambdaStar_e = pow(OutputPar_e[1]*(1 - OutputPar[5] + OutputPar[6]),2);
 931         lambdaStar_e = sqrt(lambdaStar_e);
 932       }
 933       cout<<"lambda* = "<<lambdaStar<<" +- "<<lambdaStar_e<<endl;
 934     }
 935     //if(ft==0) {MainFitParams[0]=OutputPar[3]; MainFitParams[2]=OutputPar[1];}
 936     //else {MainFitParams[4]=OutputPar[3]; MainFitParams[6]=OutputPar[1];}
 937   }// ft
 938
 939
 940   TH1D *C2_ss=(TH1D*)Two_ex_clone_mm->Clone();
 941   TH1D *C2_os=(TH1D*)Two_ex_clone_mp->Clone();
 942   TH1D *C2_ss_Momsys=(TH1D*)Two_ex_clone_mm->Clone();
 943   TH1D *C2_os_Momsys=(TH1D*)Two_ex_clone_mp->Clone();
 944   TH1D *C2_ss_Ksys=(TH1D*)Two_ex_clone_mm->Clone();
 945   TH1D *C2_os_Ksys=(TH1D*)Two_ex_clone_mp->Clone();
 946   TH1D *K2_ss = (TH1D*)Two_ex_clone_mm->Clone();
 947   TH1D *K2_os = (TH1D*)Two_ex_clone_mp->Clone();
 948
 949   for(int i=0; i<BINRANGE_C2global; i++){
 950     C2_ss->SetBinContent(i+1, C2ssFitting[i]);
 951     C2_os->SetBinContent(i+1, C2osFitting[i]);
 952     C2_ss_Momsys->SetBinContent(i+1, C2ssFitting[i]);
 953     C2_os_Momsys->SetBinContent(i+1, C2osFitting[i]);
 954     C2_ss_Ksys->SetBinContent(i+1, C2ssFitting[i]);
 955     C2_os_Ksys->SetBinContent(i+1, C2osFitting[i]);
 956     double Toterror_ss = sqrt(fabs(pow(C2ssFitting_e[i],2)-pow(C2ssSys_e[i],2)));
 957     double Toterror_os = sqrt(fabs(pow(C2osFitting_e[i],2)-pow(C2osSys_e[i],2)));
 958     C2_ss->SetBinError(i+1, Toterror_ss);
 959     C2_os->SetBinError(i+1, Toterror_os);
 960     C2_ss_Momsys->SetBinError(i+1, C2ssSys_e[i]);
 961     C2_os_Momsys->SetBinError(i+1, C2osSys_e[i]);
 962     C2_ss_Ksys->SetBinError(i+1, OutputPar[1]*K2SS_e[i]);
 963     C2_os_Ksys->SetBinError(i+1, OutputPar[1]*K2OS_e[i]);
 964     //
 965     K2_ss->SetBinContent(i+1, K2SS[i]); K2_ss->SetBinError(i+1, 0);
 966     K2_os->SetBinContent(i+1, K2OS[i]); K2_os->SetBinError(i+1, 0);
 967   }
 968
 969   C2_ss_Momsys->SetMarkerSize(0);
 970   C2_ss_Momsys->SetFillStyle(1000);
 971   C2_ss_Momsys->SetFillColor(kRed-10);
 972   C2_os_Momsys->SetMarkerSize(0);
 973   C2_os_Momsys->SetFillStyle(1000);
 974   C2_os_Momsys->SetFillColor(kBlue-10);
 975   C2_ss_Ksys->SetMarkerSize(0);
 976   C2_ss_Ksys->SetFillStyle(3004);
 977   C2_ss_Ksys->SetFillColor(kRed);
 978   C2_os_Ksys->SetMarkerSize(0);
 979   C2_os_Ksys->SetFillStyle(3004);
 980   C2_os_Ksys->SetFillColor(kBlue);
 981
 982
 983   C2_ss->GetXaxis()->SetRangeUser(0,Q2Limit);// 0,0.098
 984   C2_ss->GetYaxis()->SetRangeUser(0.95,2.4);//0.98,1.3
 985   C2_ss->GetXaxis()->SetTitleOffset(.8);
 986   C2_ss->GetYaxis()->SetTitleOffset(.85);
 987   C2_ss->GetXaxis()->SetTitle("#font[12]{q} (GeV/c)");
 988   C2_ss->GetYaxis()->SetTitle("#font[12]{C}_{2}^{#pm#pm}");
 989   C2_ss->GetXaxis()->SetTitleSize(0.05);
 990   C2_ss->GetYaxis()->SetTitleSize(0.05);
 991   C2_os->GetXaxis()->SetRangeUser(0,0.6);
 992   C2_os->GetYaxis()->SetRangeUser(0.98,1.3);
 993   C2_os->GetXaxis()->SetTitleOffset(.8);
 994   C2_os->GetYaxis()->SetTitleOffset(.8);
 995   C2_os->GetXaxis()->SetTitle("#font[12]{q} (GeV/c)");
 996   C2_os->GetXaxis()->SetTitleSize(0.05);
 997   C2_os->GetYaxis()->SetTitleSize(0.05);
 998
 999   //C2_ss->SetMarkerSize(1.5);
1000   //C2_os->SetMarkerSize(1.5);
1001   C2_os->SetMarkerStyle(25);
1002   C2_os->SetMarkerColor(4);
1003
1004
1005   fitC2ss_h->SetLineWidth(2);
1006   fitC2ss_h->SetLineColor(2);
1007
1008
1009   if(!MCcase){
1010
1011     C2_ss->DrawCopy();
1012     legend1->AddEntry(C2_ss,"same-charge","p");
1013     //C2_os->DrawCopy("same");
1014     legend1->AddEntry(C2_os,"mixed-charge","p");
1015
1016     Specif_System->Draw("same");
1017     Specif_kT->Draw("same");
1018     Specif_FSI->Draw("same");
1019     Specif_Disclaimer->Draw("same");
1020
1021     TLatex *Stats1 = new TLatex(.5,.55, "#lambda = 0.41 #pm 0.004");
1022     Stats1->SetNDC();
1023     Stats1->SetTextSize(0.06);
1024     //Stats1->Draw("same");
1025     TLatex *Stats2 = new TLatex(.5,.45, "R_{inv} = 1.59 #pm 0.01 fm");
1026     Stats2->SetNDC();
1027     Stats2->SetTextSize(0.06);
1028     //Stats2->Draw("same");
1029     /*TF1 *BkgFitC2 = new TF1("BkgFitC2","1+[0]*exp(-pow([1]*x/0.19733,2))",0,1);
1030     BkgFitC2->SetParameter(0,0.08);
1031     BkgFitC2->SetParameter(1,0.5);
1032     BkgFitC2->SetLineColor(1);
1033     C2_ss->Fit(BkgFitC2,"IME","",0.2,0.8);
1034     BkgFitC2->Draw("same");*/
1035
1036     Unity->Draw("same");
1037     //fitC2ss_h->Draw("C same");
1038     //fitC2ss_Base->Draw("C same");
1039     //fitC2ss_Expan->Draw("C same");
1040     //legend1->AddEntry(fitC2ss_Base,"Gaussian","l");
1041     //legend1->AddEntry(fitC2ss_Expan,"Edgeworth","l");
1042     //legend1->Draw("same");
1043   }
1044
1045
1046
1047   cout<<"============================================"<<endl;
1048   cout<<"Start 3-pion section"<<endl;
1049
1050   TCanvas *can2 = new TCanvas("can2", "can2",800,0,800,800);//800,0,600,900
1051   can2->SetHighLightColor(2);
1052   gStyle->SetOptFit(0111);
1053   can2->SetFillColor(10);
1054   can2->SetBorderMode(0);
1055   can2->SetBorderSize(2);
1056   can2->SetGridx();
1057   can2->SetGridy();
1058   can2->SetFrameFillColor(0);
1059   can2->SetFrameBorderMode(0);
1060   can2->SetFrameBorderMode(0);
1061   can2->cd();
1062   gPad->SetRightMargin(0.02);
1063
1064   TLegend *legend2 = new TLegend(.58,.55, .97,.65,NULL,"brNDC");
1065   legend2->SetBorderSize(1);
1066   legend2->SetTextSize(.03);// small .03; large .06
1067   legend2->SetFillColor(0);
1068
1069   /////////////////////////////////////////////
1070   TH3D *C3QS_3d = new TH3D("C3QS_3d","",BINRANGE_3,0,.4, BINRANGE_3,0,.4, BINRANGE_3,0,.4);
1071   TH3D *Combinatorics_3d = new TH3D("Combinatorics_3d","",BINRANGE_3,0,.4, BINRANGE_3,0,.4, BINRANGE_3,0,.4);
1072
1073   //
1074   const float Q3HistoLimit=0.5;
1075   const int Q3BINS = int((Q3HistoLimit+0.0001)/0.01);
1076   TH1D *num_withRS = new TH1D("num_withRS","",Q3BINS,0,Q3HistoLimit);
1077   TH1D *num_withGRS = new TH1D("num_withGRS","",Q3BINS,0,Q3HistoLimit);
1078   TH1D *num_withTM = new TH1D("num_withTM","",Q3BINS,0,Q3HistoLimit);
1079   TH1D *Cterm1 = new TH1D("Cterm1","",Q3BINS,0,Q3HistoLimit);
1080   TH1D *Cterm1_noMRC = new TH1D("Cterm1_noMRC","",Q3BINS,0,Q3HistoLimit);
1081   TH1D *Cterm2 = new TH1D("Cterm2","",Q3BINS,0,Q3HistoLimit);
1082   TH1D *Cterm3 = new TH1D("Cterm3","",Q3BINS,0,Q3HistoLimit);
1083   TH1D *Cterm4 = new TH1D("Cterm4","",Q3BINS,0,Q3HistoLimit);
1084   TH1D *num_QS = new TH1D("num_QS","",Q3BINS,0,Q3HistoLimit);
1085   TH1D *Combinatorics_1d = new TH1D("Combinatorics_1d","",Q3BINS,0,Q3HistoLimit);
1086   TH1D *Combinatorics_1d_noMRC = new TH1D("Combinatorics_1d_noMRC","",Q3BINS,0,Q3HistoLimit);
1087   TH1D *Coul_Riverside = new TH1D("Coul_Riverside","",Q3BINS,0,Q3HistoLimit);
1088   TH1D *Coul_GRiverside = new TH1D("Coul_GRiverside","",Q3BINS,0,Q3HistoLimit);
1089   TH1D *c3_hist = new TH1D("c3_hist","",Q3BINS,0,Q3HistoLimit);
1090
1091   if(SameCharge) {Cterm1_noMRC->Sumw2(); Combinatorics_1d_noMRC->Sumw2();}
1092   //
1093   double num_QS_e[Q3BINS];
1094   double c3_e[Q3BINS];
1095
1096   for(int ii=0; ii<Q3BINS; ii++){
1097     num_QS_e[ii]=0.;
1098     c3_e[ii]=0.;
1099   }
1100
1101   // CB=Charge Bin. 0 means pi-, 1 means pi+
1102   int CB1=0, CB2=0, CB3=0;
1103   int CP12=1, CP13=1, CP23=1;
1104   int CB1_2=0, CB2_2=0, CB3_2=0;
1105   if(CHARGE==-1) {
1106     if(SameCharge) {CB1=0; CB2=0; CB3=0;}
1107     else {CB1=0; CB2=0; CB3=1; CP12=+1, CP13=-1, CP23=-1;}
1108   }else {
1109     if(SameCharge) {CB1=1; CB2=1; CB3=1;}
1110     else {CB1=0; CB2=1; CB3=1; CP12=-1, CP13=-1, CP23=+1;}
1111   }
1112   if(AddCC){
1113     if(CHARGE==-1) {
1114       if(SameCharge) {CB1_2=1; CB2_2=1; CB3_2=1;}
1115       else {CB1_2=0; CB2_2=1; CB3_2=1;};
1116     }else {
1117       if(SameCharge) {CB1_2=0; CB2_2=0; CB3_2=0;}
1118       else {CB1_2=0; CB2_2=0; CB3_2=1;}
1119     }
1120   }
1121
1122
1123   // SC = species combination.  Always 0 meaning pi-pi-pi.
1124   int SCBin=0;
1125   //
1126
1127   TH1D *GenSignalExpected_num=new TH1D("GenSignalExpected_num","",Q3BINS,0,Q3HistoLimit);
1128   TH1D *GenSignalExpected_den=new TH1D("GenSignalExpected_den","",Q3BINS,0,Q3HistoLimit);
1129   //
1130   double value_num;
1131   double value_num_e;
1132   double N3_QS;
1133   double N3_QS_e;
1134
1135   for(int i=2; i<=BINLIMIT_3; i++){// bin number
1136     double Q12 = BinCenters[i-1];// true center
1137     //int Q12bin = int(Q12/0.01)+1;
1138     //
1139     for(int j=2; j<=BINLIMIT_3; j++){// bin number
1140       double Q13 = BinCenters[j-1];// true center
1141       //int Q13bin = int(Q13/0.01)+1;
1142       //
1143       for(int k=2; k<=BINLIMIT_3; k++){// bin number
1144         double Q23 = BinCenters[k-1];// true center
1145         //int Q23bin = int(Q23/0.01)+1;
1146         //
1147
1148         double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
1149         int Q3bin = Cterm1->GetXaxis()->FindBin(Q3);
1150
1151         if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
1152         if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
1153
1154
1155         //
1156         double G_12 = Gamov(CP12, Q12);// point-source Coulomb correlation
1157         double G_13 = Gamov(CP13, Q13);
1158         double G_23 = Gamov(CP23, Q23);
1159         double K2_12 = CoulCorr2(CP12, Q12);// finite-source Coulomb+Strong correlation from Therminator.
1160         double K2_13 = CoulCorr2(CP13, Q13);
1161         double K2_23 = CoulCorr2(CP23, Q23);
1162         double K3 = 1.;// 3-body Coulomb+Strong correlation
1163         if(SameCharge || CHARGE==-1) K3 = CoulCorrGRS(SameCharge, Q12, Q13, Q23);
1164         else K3 = CoulCorrGRS(SameCharge, Q23, Q12, Q13);
1165
1166         if(MCcase && !GeneratedSignal) { K2_12=1.; K2_13=1.; K2_23=1.; K3=1.;}
1167         if(K3==0) continue;
1168
1169         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
1170         double TERM2=Three_3d[CB1][CB2][CB3][SCBin][1]->GetBinContent(i,j,k);// N2*N1. 1 and 2 from same-event
1171         double TERM3=Three_3d[CB1][CB2][CB3][SCBin][2]->GetBinContent(i,j,k);// N2*N1. 1 and 3 from same-event
1172         double TERM4=Three_3d[CB1][CB2][CB3][SCBin][3]->GetBinContent(i,j,k);// N2*N1. 2 and 3 from same-event
1173         double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
1174         if(AddCC){
1175           if(CHARGE==-1){// base is (SC,MC,MC)
1176             TERM1 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][0]->GetBinContent(j,k,i);
1177             TERM2 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][1]->GetBinContent(j,k,i);
1178             TERM3 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][2]->GetBinContent(j,k,i);
1179             TERM4 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][3]->GetBinContent(j,k,i);
1180             TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(j,k,i);
1181           }else {// base is (MC,MC,SC)
1182             TERM1 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][0]->GetBinContent(k,i,j);
1183             TERM2 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][1]->GetBinContent(k,i,j);
1184             TERM3 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][2]->GetBinContent(k,i,j);
1185             TERM4 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][3]->GetBinContent(k,i,j);
1186             TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(k,i,j);
1187           }
1188         }
1189
1190         if(TERM1==0 && TERM2==0 && TERM3==0 && TERM4==0 && TERM5==0) continue;
1191         if(TERM1==0 || TERM2==0 || TERM3==0 || TERM4==0 || TERM5==0) continue;
1192
1193         double muonCorr_C3=1.0, muonCorr_C2_12=1.0, muonCorr_C2_13=1.0, muonCorr_C2_23=1.0;
1194         if(MuonCorrection && !MCcase){
1195           if(SameCharge){
1196             muonCorr_C3 = C3muonCorrection->GetBinContent(C3muonCorrection->GetXaxis()->FindBin(Q3));
1197             muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q12));
1198             muonCorr_C2_13 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q13));
1199             muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q23));
1200           }else{
1201             muonCorr_C3 = C3muonCorrection->GetBinContent(C3muonCorrection->GetXaxis()->FindBin(Q3));
1202             if(CHARGE==-1){// pi-pi-pi+
1203               muonCorr_C2_12 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q12));
1204               muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q13));
1205               muonCorr_C2_23 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q23));
1206             }else{// pi-pi+pi+
1207               muonCorr_C2_12 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q12));
1208               muonCorr_C2_13 = C2muonCorrection_MC->GetBinContent(C2muonCorrection_MC->GetXaxis()->FindBin(Q13));
1209               muonCorr_C2_23 = C2muonCorrection_SC->GetBinContent(C2muonCorrection_SC->GetXaxis()->FindBin(Q23));
1210             }
1211           }
1212         }
1213
1214         // apply momentum resolution correction
1215         if(!MCcase && !GeneratedSignal){
1216           int MRC_Q3bin = int(Q3/0.01) + 1;
1217           if(SameCharge){
1218             TERM1 *= (MomRes1d[0][0]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1219             TERM2 *= (MomRes1d[0][1]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1220             TERM3 *= (MomRes1d[0][2]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1221             TERM4 *= (MomRes1d[0][3]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1222             TERM5 *= (MomRes1d[0][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1223           }else{
1224             // removed due to over-correction of 1-D MRC approximation
1225             TERM1 *= (MomRes1d[1][0]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1226             TERM2 *= (MomRes1d[1][1]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1227             TERM3 *= (MomRes1d[1][2]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1228             TERM4 *= (MomRes1d[1][3]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1229             TERM5 *= (MomRes1d[1][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1230           }
1231         }
1232         //
1233
1234         TwoFrac=lambda_PM;// 0.7
1235         OneFrac=pow(TwoFrac,.5); ThreeFrac=pow(TwoFrac,1.5);
1236
1237
1238         // Purify.  Isolate pure 3-pion QS correlations using Lambda and K3 (removes lower order correlations)
1239         N3_QS = TERM1;
1240         N3_QS -= ( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) )*TERM5;
1241         N3_QS -= (1-OneFrac)*(TERM2 + TERM3 + TERM4 - 3*(1-TwoFrac)*TERM5);
1242         N3_QS /= ThreeFrac;
1243         N3_QS /= K3;
1244         N3_QS *=  muonCorr_C3;
1245
1246         num_QS->Fill(Q3, N3_QS);
1247
1248         // Isolate 3-pion cumulant
1249         value_num = N3_QS;
1250         value_num -= (TERM2 - (1-TwoFrac)*TERM5)/K2_12/TwoFrac * muonCorr_C2_12;
1251         value_num -= (TERM3 - (1-TwoFrac)*TERM5)/K2_13/TwoFrac * muonCorr_C2_13;
1252         value_num -= (TERM4 - (1-TwoFrac)*TERM5)/K2_23/TwoFrac * muonCorr_C2_23;
1253         value_num += 2*TERM5;
1254
1255
1256
1257
1258         // errors
1259         N3_QS_e = TERM1;
1260         N3_QS_e += pow(( pow(1-OneFrac,3) +3*OneFrac*pow(1-OneFrac,2) )*sqrt(TERM5),2);
1261         N3_QS_e += (pow((1-OneFrac),2)*(TERM2 + TERM3 + TERM4) + pow((1-OneFrac)*3*(1-TwoFrac)*sqrt(TERM5),2));
1262         N3_QS_e /= pow(ThreeFrac,2);
1263         N3_QS_e /= pow(K3,2);
1264         num_QS_e[Q3bin-1] += N3_QS_e;
1265         // errors
1266         value_num_e = N3_QS_e;
1267         value_num_e += (pow(1/K2_12/TwoFrac*sqrt(TERM2),2) + pow((1-TwoFrac)/K2_12/TwoFrac*sqrt(TERM5),2));
1268         value_num_e += (pow(1/K2_13/TwoFrac*sqrt(TERM3),2) + pow((1-TwoFrac)/K2_13/TwoFrac*sqrt(TERM5),2));
1269         value_num_e += (pow(1/K2_23/TwoFrac*sqrt(TERM4),2) + pow((1-TwoFrac)/K2_23/TwoFrac*sqrt(TERM5),2));
1270         value_num_e += pow(2*sqrt(TERM5),2);
1271
1272         c3_e[Q3bin-1] += value_num_e + TERM5;// add baseline stat error
1273
1274
1275         // Fill histograms
1276         c3_hist->Fill(Q3, value_num + TERM5);// for cumulant correlation function
1277         C3QS_3d->SetBinContent(i,j,k, N3_QS);
1278         C3QS_3d->SetBinError(i,j,k, N3_QS_e);
1279         //
1280         Coul_Riverside->Fill(Q3, G_12*G_13*G_23*TERM5);
1281         Coul_GRiverside->Fill(Q3, TERM5*CoulCorrGRS(SameCharge, Q12, Q13, Q23));
1282         num_withGRS->Fill(Q3, N3_QS);
1283         Cterm1->Fill(Q3, TERM1);
1284         Cterm2->Fill(Q3, TERM2);
1285         Cterm3->Fill(Q3, TERM3);
1286         Cterm4->Fill(Q3, TERM4);
1287         Combinatorics_1d->Fill(Q3, TERM5);
1288         Combinatorics_3d->Fill(Q12,Q13,Q23, TERM5);
1289
1290         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));
1291         Cterm1_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][0]->GetBinContent(i,j,k));
1292         Combinatorics_1d_noMRC->Fill(Q3_m, Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k));
1293
1294         //
1295         A_3[i-1][j-1][k-1] = value_num + TERM5;
1296         B_3[i-1][j-1][k-1] = TERM5;
1297         A_3_e[i-1][j-1][k-1] = sqrt(value_num_e + TERM5);
1298
1299         ///////////////////////////////////////////////////////////
1300
1301       }
1302     }
1303   }
1304
1305
1306
1307   ////////////////////////////
1308
1309   // Intermediate steps
1310   num_withRS->Sumw2();
1311   num_withGRS->Sumw2();
1312   num_withTM->Sumw2();
1313   Cterm1->Sumw2();
1314   Cterm2->Sumw2();
1315   Cterm3->Sumw2();
1316   Cterm4->Sumw2();
1317   Combinatorics_1d->Sumw2();
1318   Combinatorics_3d->Sumw2();
1319   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]));}
1320
1321
1322   num_withRS->Divide(Combinatorics_1d);
1323   num_withGRS->Divide(Combinatorics_1d);
1324   num_withTM->Divide(Combinatorics_1d);
1325   Cterm1->Divide(Combinatorics_1d);
1326   Cterm1_noMRC->Divide(Combinatorics_1d_noMRC);
1327   Cterm2->Divide(Combinatorics_1d);
1328   Cterm3->Divide(Combinatorics_1d);
1329   Cterm4->Divide(Combinatorics_1d);
1330   c3_hist->Divide(Combinatorics_1d);
1331   num_QS->Divide(Combinatorics_1d);
1332
1333
1334
1335   ///////////////////////////////////////////////////////////////////////////////////////////////////
1336
1337
1338   Cterm1->SetMarkerStyle(20);
1339   Cterm1->SetMarkerColor(4);
1340   Cterm1->SetLineColor(4);
1341   Cterm1->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})");
1342   Cterm1->GetYaxis()->SetTitle("#font[12]{C}_{3}");
1343   //Cterm1->SetTitle("#pi^{-}#pi^{-}#pi^{-}");
1344   Cterm1->SetMinimum(0.97);
1345   Cterm1->SetMaximum(5.7);// 6.1 for same-charge
1346   Cterm1->GetXaxis()->SetRangeUser(0,Q3Limit);
1347   Cterm1->GetYaxis()->SetTitleOffset(1.4);
1348   Cterm1->Draw();
1349   legend2->AddEntry(Cterm1,"#font[12]{C}_{3} raw","p");
1350   //
1351   Cterm2->SetMarkerStyle(20);
1352   Cterm2->SetMarkerColor(7);
1353
1354   Cterm3->SetMarkerStyle(25);
1355   Cterm3->SetMarkerColor(8);
1356   Cterm4->SetMarkerStyle(26);
1357   Cterm4->SetMarkerColor(9);
1358   //Cterm2->Draw("same");
1359   //Cterm3->Draw("same");
1360   //Cterm4->Draw("same");
1361   //legend2->AddEntry(Cterm1,"++-","p");
1362
1363
1364   if(MCcase){
1365     double C3points_HIJING_mmm[10]={0, 0.961834, 1.01827, 0.999387, 1.00202, 1.00081, 1.00082, 1.00037, 0.999074, 0.999099};
1366     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};
1367     double C3points_HIJING_ppp[10]={0, 1.13015, 1.00623, 1.00536, 1.00293, 0.999964, 1.00116, 1.0007, 1.00037, 1.00105};
1368     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};
1369     TH1D *C3_HIJING_mmm=(TH1D*)Cterm1->Clone();
1370     TH1D *C3_HIJING_ppp=(TH1D*)Cterm1->Clone();
1371     for(int i=0; i<10; i++){
1372       C3_HIJING_mmm->SetBinContent(i+1, C3points_HIJING_mmm[i]);
1373       C3_HIJING_mmm->SetBinError(i+1, C3points_HIJING_mmm_e[i]);
1374       C3_HIJING_ppp->SetBinContent(i+1, C3points_HIJING_ppp[i]);
1375       C3_HIJING_ppp->SetBinError(i+1, C3points_HIJING_ppp_e[i]);
1376     }
1377     C3_HIJING_mmm->SetMarkerColor(2);
1378     C3_HIJING_mmm->SetLineColor(2);
1379     C3_HIJING_ppp->SetMarkerColor(4);
1380     C3_HIJING_ppp->SetLineColor(4);
1381     //C3_HIJING_mmm->Draw("same");
1382     //C3_HIJING_ppp->Draw("same");
1383     //legend2->AddEntry(C3_HIJING_mmm,"---","p");
1384     //legend2->AddEntry(C3_HIJING_ppp,"+++","p");
1385   }
1386
1387   num_QS->SetMarkerStyle(24);
1388   num_QS->SetMarkerColor(4);
1389   num_QS->SetLineColor(4);
1390   num_QS->GetXaxis()->SetTitle("#font[12]{Q}_{3}");
1391   num_QS->GetYaxis()->SetTitle("#font[12]{C}_{3}^{QS}");
1392   num_QS->GetXaxis()->SetRangeUser(0,.12);
1393   num_QS->SetMaximum(6);
1394   //num_QS->Draw("same");
1395   //legend2->AddEntry(num_QS,"C_{3}^{QS}","p");
1396
1397   c3_hist->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})");
1398   c3_hist->GetYaxis()->SetTitle("#font[12]{C}_{3} or #font[12]{#bf{c}}_{3}");
1399   c3_hist->GetYaxis()->SetTitleOffset(1.4);
1400   c3_hist->GetXaxis()->SetRangeUser(0,Q3Limit);
1401   c3_hist->GetYaxis()->SetRangeUser(0.9,4);
1402   c3_hist->SetMarkerStyle(25);
1403   c3_hist->SetMarkerColor(2);
1404   c3_hist->SetLineColor(2);
1405   c3_hist->SetMaximum(3);
1406   c3_hist->SetMinimum(.7);
1407   if(!MCcase) c3_hist->Draw("same");
1408   legend2->AddEntry(c3_hist,"#font[12]{#bf{c}}_{3} (cumulant correlation)","p");
1409   c3_hist->Draw();
1410
1411   //for(int i=1; i<=15; i++) cout<<c3_hist->GetBinContent(i)<<", ";
1412   //cout<<endl;
1413   //for(int i=1; i<=15; i++) cout<<c3_hist->GetBinError(i)<<", ";
1414   //cout<<endl;
1415
1416   const int npar_c3=5;// 10
1417   TMinuit MyMinuit_c3(npar_c3);
1418   MyMinuit_c3.SetFCN(fcn_c3);
1419   int ierflg_c3=0;
1420   double arglist_c3 = 0;
1421   MyMinuit_c3.mnexcm("SET NOWarnings",&arglist_c3,1, ierflg_c3);
1422   arglist_c3 = -1;
1423   MyMinuit_c3.mnexcm("SET PRint",&arglist_c3,1, ierflg_c3);
1424   //arglist_c3=2;// improve Minimization Strategy (1 is default)
1425   //MyMinuit_c3.mnexcm("SET STR",&arglist_c3,1,ierflg_c3);
1426   //arglist_c3 = 0;
1427   //MyMinuit_c3.mnexcm("SCAN", &arglist_c3,1,ierflg_c3);
1428   arglist_c3 = 5000;
1429   MyMinuit_c3.mnexcm("MIGRAD", &arglist_c3 ,1,ierflg_c3);
1430
1431   TF1 *c3Fit1D_Base=new TF1("c3Fit1D_Base","[0]*(1+[1]*exp(-pow([2]*x/0.19733,[3]) * [4]/2.))",0,1);
1432   c3Fit1D_Base->FixParameter(3,ExpPower);
1433   TF1 *c3Fit1D_Expan;
1434   if(Gaussian) {
1435     c3Fit1D_Base->FixParameter(4,1.0);
1436     c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)),3))",0,1);
1437   }else {
1438     c3Fit1D_Base->FixParameter(4,sqrt(3.));
1439     c3Fit1D_Expan=new TF1("c3Fit1D_Expan","[0]*(1+[1]*exp(-pow([2]*x/0.19733,1) * sqrt(3.)/2.) * pow(1 + [3]*([2]*x/sqrt(3.)/0.19733 - 1) + [4]/2*(pow([2]*x/sqrt(3.)/0.19733,2) - 4*[2]*x/sqrt(3.)/0.19733 + 2),3))",0,1);
1440   }
1441   double OutputPar_c3[npar_c3]={0};
1442   double OutputPar_c3_e[npar_c3]={0};
1443
1444   double par_c3[npar_c3];               // the start values
1445   double stepSize_c3[npar_c3];          // step sizes
1446   double minVal_c3[npar_c3];            // minimum bound on parameter
1447   double maxVal_c3[npar_c3];            // maximum bound on parameter
1448   string parName_c3[npar_c3];
1449   if(SameCharge && !MCcase){
1450     for(int ft=0; ft<2; ft++){// Gaussian, Edgeworth
1451       if(ft==1 && !IncludeExpansion) continue;
1452       par_c3[0] = 1.0; par_c3[1] = .5; par_c3[2] = 5; par_c3[3] = kappa3; par_c3[4] = kappa4;// was par_c3[3] = 0.; par_c3[4] = 0.;
1453       stepSize_c3[0] = 0.01; stepSize_c3[1] = 0.1; stepSize_c3[2] = 0.2; stepSize_c3[3] = 0.01; stepSize_c3[4] = 0.01;
1454       minVal_c3[0] = 0.95; minVal_c3[1] = 0.2; minVal_c3[2] = 0.5; minVal_c3[3] = -1; minVal_c3[4] = -1;// was minVal_c3[3] = -1; minVal_c3[4] = -1;
1455       maxVal_c3[0] = 1.1; maxVal_c3[1] = 2.5; maxVal_c3[2] = 15.; maxVal_c3[3] = +1; maxVal_c3[4] = +1;// was minVal_c3[3] = +1; minVal_c3[4] = +1;
1456       parName_c3[0] = "N"; parName_c3[1] = "#lambda_{3}"; parName_c3[2] = "R_{inv}"; parName_c3[3] = "coeff_{3}"; parName_c3[4] = "coeff_{4}";
1457       c3Fit1D_Base->SetParName(0,"N"); c3Fit1D_Base->SetParName(1,"#lambda_{3}"); c3Fit1D_Base->SetParName(2,"R_{inv}");
1458       for (int i=0; i<npar_c3; i++){
1459         MyMinuit_c3.DefineParameter(i, parName_c3[i].c_str(), par_c3[i], stepSize_c3[i], minVal_c3[i], maxVal_c3[i]);
1460       }
1461       //
1462       if(!FixExpansionAvg) {
1463         maxVal_c3[1] = 2.0;
1464         double RadiusFix;
1465         if(CollisionType==0) RadiusFix = 10 - 7*Mbin/15.;
1466         else RadiusFix = 2.5 - 1.25*(Mbin-12)/7.;
1467         MyMinuit_c3.DefineParameter(2, parName_c3[2].c_str(), RadiusFix, stepSize_c3[2], minVal_c3[2], maxVal_c3[2]); MyMinuit_c3.FixParameter(2);
1468         MyMinuit_c3.DefineParameter(1, parName_c3[1].c_str(), 2.0, stepSize_c3[1], minVal_c3[1], maxVal_c3[1]); MyMinuit_c3.FixParameter(1);
1469       }
1470       // kappa_3 and kappa_4
1471       if(ft==0){// Gaussian
1472         MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), 0., stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3);
1473         MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), 0., stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
1474       }else{// Edgeworth
1475         if(FixExpansionAvg){
1476           MyMinuit_c3.DefineParameter(3, parName_c3[3].c_str(), kappa3, stepSize_c3[3], minVal_c3[3], maxVal_c3[3]); MyMinuit_c3.FixParameter(3);
1477           MyMinuit_c3.DefineParameter(4, parName_c3[4].c_str(), kappa4, stepSize_c3[4], minVal_c3[4], maxVal_c3[4]); MyMinuit_c3.FixParameter(4);
1478         }else{
1479           Int_t err=0;
1480           Double_t tmp[1];
1481           tmp[0] = 3+1;
1482           MyMinuit_c3.mnexcm( "RELease", tmp,  1,  err );// kappa3
1483           tmp[0] = 4+1;
1484           MyMinuit_c3.mnexcm( "RELease", tmp,  1,  err );// kappa4
1485         }
1486       }
1487
1488
1489       /////////////////////////////////////////////////////////////
1490       // Do the minimization!
1491       cout<<"Start Three-d fit"<<endl;
1492       MyMinuit_c3.Migrad();// Minuit's best minimization algorithm
1493       cout<<"End Three-d fit"<<endl;
1494       /////////////////////////////////////////////////////////////
1495       MyMinuit_c3.mnexcm("SHOw PARameters", &arglist_c3, 1, ierflg);
1496       cout<<"c3 Fit: Chi2/NDF = "<<Chi2_c3/(NFitPoints_c3-MyMinuit_c3.GetNumFreePars())<<endl;
1497
1498       for(int i=0; i<npar_c3; i++){
1499         MyMinuit_c3.GetParameter(i,OutputPar_c3[i],OutputPar_c3_e[i]);
1500       }
1501       if(ft==0){
1502         c3Fit1D_Base->SetParameter(0,OutputPar_c3[0]); c3Fit1D_Base->SetParError(0,OutputPar_c3_e[0]);
1503         c3Fit1D_Base->SetParameter(1,OutputPar_c3[1]); c3Fit1D_Base->SetParError(1,OutputPar_c3_e[1]);
1504         c3Fit1D_Base->SetParameter(2,OutputPar_c3[2]); c3Fit1D_Base->SetParError(2,OutputPar_c3_e[2]);
1505         c3Fit1D_Base->SetLineStyle(2);
1506         c3Fit1D_Base->Draw("same");
1507       }else{
1508         c3Fit1D_Expan->SetParameter(0,OutputPar_c3[0]); c3Fit1D_Expan->SetParError(0,OutputPar_c3_e[0]);
1509         c3Fit1D_Expan->SetParameter(1,OutputPar_c3[1]); c3Fit1D_Expan->SetParError(1,OutputPar_c3_e[1]);
1510         c3Fit1D_Expan->SetParameter(2,OutputPar_c3[2]); c3Fit1D_Expan->SetParError(2,OutputPar_c3_e[2]);
1511         c3Fit1D_Expan->SetParameter(3,OutputPar_c3[3]); c3Fit1D_Expan->SetParError(3,OutputPar_c3_e[3]);
1512         c3Fit1D_Expan->SetParameter(4,OutputPar_c3[4]); c3Fit1D_Expan->SetParError(4,OutputPar_c3_e[4]);
1513         c3Fit1D_Expan->SetLineStyle(1);
1514         c3Fit1D_Expan->Draw("same");
1515         double lambdaStar=0, lambdaStar_e=0;
1516         if(Gaussian){
1517           lambdaStar = OutputPar_c3[1]*pow(1 + OutputPar_c3[4]/8.,3);
1518           lambdaStar_e = pow(OutputPar_c3_e[1]*pow(1 + OutputPar_c3[4]/8.,3),2);
1519           lambdaStar_e = sqrt(lambdaStar_e);
1520         }else{
1521           lambdaStar = OutputPar_c3[1]*pow(1 - OutputPar_c3[3] + OutputPar_c3[4],3);
1522           lambdaStar_e = pow(OutputPar_c3_e[1]*pow(1 - OutputPar_c3[3] + OutputPar_c3[4],3),2);
1523           lambdaStar_e = sqrt(lambdaStar_e);
1524         }
1525         cout<<"lambda*_3 = "<<lambdaStar<<" +- "<<lambdaStar_e<<endl;
1526       }
1527       //if(ft==0) {MainFitParams[1]=OutputPar_c3[2]; MainFitParams[3]=OutputPar_c3[1];}
1528       //else {MainFitParams[5]=OutputPar_c3[2]; MainFitParams[7]=OutputPar_c3[1];}
1529     }// fit type
1530   }// SC and !MC
1531
1532   // project 3D EW fit onto 1D histogram
1533   if(SameCharge && !MCcase){
1534     for(int i=2; i<=BINLIMIT_3; i++){// bin number
1535       double Q12 = BinCenters[i-1];// true center
1536       for(int j=2; j<=BINLIMIT_3; j++){// bin number
1537         double Q13 = BinCenters[j-1];// true center
1538         for(int k=2; k<=BINLIMIT_3; k++){// bin number
1539           double Q23 = BinCenters[k-1];// true center
1540           //
1541           double Q3 = sqrt(pow(Q12,2) + pow(Q13,2) + pow(Q23,2));
1542
1543           if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) continue;// not all configurations are possible
1544           if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) continue;// not all configurations are possible
1545
1546           double TERM5=Three_3d[CB1][CB2][CB3][SCBin][4]->GetBinContent(i,j,k);// N1*N1*N1. All from different events (pure combinatorics)
1547           if(TERM5==0) continue;
1548           //
1549           if(AddCC){
1550             if(CHARGE==-1){// base is (SC,MC,MC)
1551               TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(j,k,i);
1552             }else {// base is (MC,MC,SC)
1553               TERM5 += Three_3d[CB1_2][CB2_2][CB3_2][SCBin][4]->GetBinContent(k,i,j);
1554             }
1555           }
1556           //
1557           int MRC_Q3bin = int(Q3/0.01) + 1;
1558           TERM5 *= (MomRes1d[0][4]->GetBinContent(MRC_Q3bin)-1)*MRCShift+1;
1559           //
1560           double radius_exp = c3Fit1D_Expan->GetParameter(2)/FmToGeV;
1561           double arg12 = Q12*radius_exp;
1562           double arg13 = Q13*radius_exp;
1563           double arg23 = Q23*radius_exp;
1564
1565           Float_t EW12 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
1566           EW12 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
1567           Float_t EW13 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
1568           EW13 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
1569           Float_t EW23 = 1 + c3Fit1D_Expan->GetParameter(3)/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
1570           EW23 += c3Fit1D_Expan->GetParameter(4)/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
1571           //
1572           double cumulant_fitvalue=0;
1573           cumulant_fitvalue = c3Fit1D_Expan->GetParameter(0)*(c3Fit1D_Expan->GetParameter(1)*EW12*EW13*EW23*exp(-pow(radius_exp,2)/2. * (pow(Q12,2)+pow(Q13,2)+pow(Q23,2)))*TERM5 + TERM5);
1574
1575           GenSignalExpected_num->Fill(Q3, cumulant_fitvalue);
1576           GenSignalExpected_den->Fill(Q3, TERM5);
1577           ///////////////////////////////////////////////////////////
1578
1579
1580
1581         }
1582       }
1583     }
1584   }
1585   GenSignalExpected_num->Sumw2();
1586   GenSignalExpected_num->Divide(GenSignalExpected_den);
1587
1588   TSpline3 *c3Fit1D_ExpanSpline = new TSpline3(GenSignalExpected_num);
1589   if(CollisionType==0) c3Fit1D_ExpanSpline->SetLineColor(1);
1590   if(CollisionType==1) c3Fit1D_ExpanSpline->SetLineColor(2);
1591   if(CollisionType==2) c3Fit1D_ExpanSpline->SetLineColor(4);
1592   c3Fit1D_ExpanSpline->SetLineWidth(2);
1593   double xpoints[1000], ypoints[1000];
1594   bool splineOnly=kFALSE;
1595   for(int iii=0; iii<1000; iii++){
1596     xpoints[iii] = 0 + (iii+0.5)*0.001;
1597     if(!splineOnly){
1598       ypoints[iii] = c3Fit1D_Expan->Eval(xpoints[iii]);
1599       if(c3Fit1D_Expan->Eval(xpoints[iii])<2 && fabs(c3Fit1D_Expan->Eval(xpoints[iii])-c3Fit1D_ExpanSpline->Eval(xpoints[iii])) < 0.01) splineOnly=kTRUE;
1600     }
1601     else {ypoints[iii] = c3Fit1D_ExpanSpline->Eval(xpoints[iii]); splineOnly=kTRUE;}
1602   }
1603   TGraph *gr_c3Spline = new TGraph(1000,xpoints,ypoints);
1604   gr_c3Spline->SetLineWidth(2);
1605   if(CollisionType==0) gr_c3Spline->SetLineColor(1);
1606   if(CollisionType==1) gr_c3Spline->SetLineColor(2);
1607   if(CollisionType==2) gr_c3Spline->SetLineColor(4);
1608   gr_c3Spline->Draw("c same");
1609   //c3Fit1D_ExpanSpline->Draw("same");
1610
1611
1612   // uncomment to display fit box
1613   //c3_hist->GetListOfFunctions()->Add(c3Fit1D_Base);
1614
1615   // Fit Range Systematics
1616   //for(int i=0; i<8; i++){cout<<MainFitParams[i]<<",  ";}
1617   //cout<<endl;
1618   //for(int i=0; i<4; i++){cout<<int(100*(MainFitParams[i]-RefMainFitParams[i])/MainFitParams[i])<<"$\\%$ & ";}// Gaussian
1619   //cout<<endl;
1620   //for(int i=4; i<8; i++){cout<<int(100*(MainFitParams[i]-RefMainFitParams[i])/MainFitParams[i])<<"$\\%$ & ";}// Edgeworth
1621   //cout<<endl;
1622
1623
1624   /*TString *lamName=new TString("#lambda_{3}=");
1625   TString *RName=new TString("#R_{inv,3}=");
1626   *lamName += int(100*c3Fit1D_Base->GetParameter(1))/100.; lamName->Append(" #pm "); *lamName += int(100*c3Fit1D_Base->GetParError(1))/100.;
1627   *RName += round(100*c3Fit1D_Base->GetParameter(2))/100; RName->Append(" #pm "); *RName += int(100*c3Fit1D_Base->GetParError(2))/100.; RName->Append(" fm");
1628   TLatex *fitBox1=new TLatex(0.5,0.9,lamName->Data());
1629   fitBox1->SetNDC(kTRUE);
1630   fitBox1->Draw();
1631   TLatex *fitBox2=new TLatex(0.5,0.8,RName->Data());
1632   fitBox2->SetNDC(kTRUE);
1633   fitBox2->Draw();*/
1634   //TLatex *fitBox3=new TLatex(0.5,0.5,"<N_{rec,pions}>=103");
1635   //fitBox3->SetNDC(kTRUE);
1636   //fitBox3->Draw();
1637
1638     //TPaveStats *c3Stats=(TPaveStats*)c3Fit1D_Base->FindObject("stats");
1639     //c3Stats->SetX1NDC(0.15);
1640     //c3Stats->SetX2NDC(0.52);
1641     //c3Stats->SetY1NDC(.2);
1642     //c3Stats->SetY2NDC(.3);
1643     //c3Stats->Draw("same");
1644
1645
1646     // The below 1D fit method has less well defined bin centers
1647     //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,1);
1648     //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + (-0.12/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + (0.43/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1);
1649     //TF1 *c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * (1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/0.19733,3) - 12.*pow([2]*x/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/0.19733,4) -48.*pow([2]*x/0.19733,2) + 12))))",0,1);
1650     //c3Fit1D->SetLineColor(4);
1651     //c3Fit1D->SetParameter(0,1); c3Fit1D->SetParName(0,"N");
1652     //c3Fit1D->SetParameter(1,0.5); c3Fit1D->SetParName(1,"#lambda_{3}");
1653     //c3Fit1D->SetParameter(2,5); c3Fit1D->SetParName(2,"R_{inv}");
1654     //c3Fit1D->SetParLimits(2,0.5,15);
1655     //c3Fit1D->SetLineColor(1);
1656     //
1657     //c3Fit1D->FixParameter(0,OutputPar_c3[0]);
1658     //c3Fit1D->FixParameter(1,OutputPar_c3[1]);
1659     //c3Fit1D->FixParameter(2,OutputPar_c3[2]);
1660     //c3Fit1D->SetParName(3,"coeff_{3}"); c3Fit1D->SetParName(4,"coeff_{4}");
1661     //c3Fit1D->SetParameter(3,.24); c3Fit1D->SetParameter(3,.16);
1662     //c3_hist->Fit(c3Fit1D,"IMEN","",0.0,Q3Limit);
1663     //c3Fit1D->Draw("same");
1664     //cout<<"1d fit: Chi2/NDF = "<<c3Fit1D->GetChisquare()/c3Fit1D->GetNDF()<<endl;
1665     //dentimesFit_c3->DrawCopy("l same");
1666     //
1667     //c3Fit1D->FixParameter(0, fitcopy_c3->GetParameter(0));
1668     //c3Fit1D->FixParameter(1, fitcopy_c3->GetParameter(1));
1669     //c3Fit1D->FixParameter(2, fitcopy_c3->GetParameter(2));
1670
1671     //TF1 *c3Fit1D_2=new TF1("c3Fit1D_2","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.))",0,1);
1672     //c3Fit1D=new TF1("c3Fit1D","[0]*(1+[1]*exp(-pow([2]*x/0.19733,2)/2.) * pow(1 + ([3]/(6.*pow(2,1.5))*(8.*pow([2]*x/sqrt(3.)/0.19733,3) - 12.*pow([2]*x/sqrt(3.)/0.19733,1))) + ([4]/(24.*pow(2,2))*(16.*pow([2]*x/sqrt(3.)/0.19733,4) -48.*pow([2]*x/sqrt(3.)/0.19733,2) + 12)),3))",0,1);
1673     //c3Fit1D->FixParameter(0,OutputPar_c3[0]);
1674     //c3Fit1D->FixParameter(1,OutputPar_c3[1]);
1675     //c3Fit1D->FixParameter(2,OutputPar_c3[2]);
1676     //c3Fit1D->FixParameter(3,OutputPar_c3[3]);
1677     //c3Fit1D->FixParameter(4,OutputPar_c3[4]);
1678     //c3Fit1D->SetLineColor(4);
1679     //c3Fit1D->Draw("same");
1680
1681
1682   //MomRes1d[1][0]->Draw();
1683
1684   //for(int i=1; i<50; i++) cout<<c3_hist->GetBinContent(i)<<", ";
1685   //cout<<endl;
1686   //for(int i=1; i<50; i++) cout<<c3_hist->GetBinError(i)<<", ";
1687   // pp M2, pi-
1688
1689
1690   Cterm1->Draw("same");
1691   legend2->Draw("same");
1692
1693   /*if(SameCharge==kTRUE && MCcase==kFALSE){
1694     TString *savename=new TString("C3c3_SC_");
1695     if(CollisionType==0) savename->Append("PbPb_K");
1696     if(CollisionType==1) savename->Append("pPb_K");
1697     if(CollisionType==2) savename->Append("pp_K");
1698     *savename += EDbin;
1699     savename->Append("_M");
1700     *savename += Mbin;
1701     savename->Append(".eps");
1702     can2->SaveAs(savename->Data());
1703     }*/
1704
1705   //for(int ii=0; ii<10; ii++) cout<<c3_hist->GetBinContent(ii+1)<<", ";
1706   //Coul_GRiverside->Draw();
1707   //Coul_Riverside->Draw("same");
1708   /*TLegend *legend4 = new TLegend(.3,.8, .5,.95,NULL,"brNDC");
1709   legend4->SetBorderSize(0);
1710   legend4->SetTextFont(42);//42
1711   legend4->SetTextSize(.04);// small .03; large .06
1712   legend4->SetFillColor(0);
1713
1714   gPad->SetRightMargin(0.025); gPad->SetLeftMargin(0.1); gPad->SetTopMargin(0.02);
1715   TwoFrac=lambda_PM; OneFrac=pow(TwoFrac,.5), ThreeFrac=pow(TwoFrac,1.5);
1716   gPad->SetGridx(0);
1717   gPad->SetGridy(0);
1718
1719   TH1D *c3_Extended = (TH1D*)Three_1d[CB1][CB2][CB3][SCBin][0]->Clone();
1720   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], -( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) ));
1721   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][1], -(1-OneFrac));
1722   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][2], -(1-OneFrac));
1723   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][3], -(1-OneFrac));
1724   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], (1-OneFrac)*3*(1-TwoFrac));
1725   c3_Extended->Scale(1/ThreeFrac);
1726   //
1727   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][1], -1/TwoFrac);
1728   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][2], -1/TwoFrac);
1729   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][3], -1/TwoFrac);
1730   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], 3*(1-TwoFrac)/TwoFrac);
1731   c3_Extended->Add(Three_1d[CB1][CB2][CB3][SCBin][4], 3);
1732   //
1733   c3_Extended->Divide(Three_1d[CB1][CB2][CB3][SCBin][4]);
1734   c3_Extended->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/#font[12]{c})");
1735   c3_Extended->GetYaxis()->SetTitle("#font[12]{#bf{c}}_{3}^{#pm#pm#pm}");
1736   c3_Extended->SetMinimum(0.9); c3_Extended->SetMaximum(2.0);
1737   c3_Extended->SetMarkerStyle(24);
1738   c3_Extended->SetMarkerColor(2);
1739   c3_Extended->SetLineColor(2);
1740
1741   //
1742   TH1D *C3_Extended = (TH1D*)Three_1d[CB1][CB2][CB3][SCBin][0]->Clone();
1743   C3_Extended->Divide(Three_1d[CB1][CB2][CB3][SCBin][4]);
1744   C3_Extended->GetXaxis()->SetTitle("#font[12]{Q}_{3} (GeV/c)");
1745   C3_Extended->GetYaxis()->SetTitle("#font[12]{C}_{3}");
1746   c3_Extended->GetXaxis()->SetTitleOffset(0.83);
1747   c3_Extended->GetYaxis()->SetTitleOffset(0.8);
1748   c3_Extended->GetYaxis()->SetTitleSize(0.05);
1749   c3_Extended->GetXaxis()->SetTitleSize(0.05);
1750   C3_Extended->SetMarkerStyle(20);
1751   C3_Extended->SetMarkerColor(4);
1752   c3_Extended->GetXaxis()->SetRangeUser(0,2);
1753   c3_Extended->SetMarkerStyle(20);
1754   c3_Extended->Draw();*/
1755   //legend4->AddEntry(c3_Extended,"c_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p");
1756   //legend4->AddEntry(C3_Extended,"C_{3}, #pi^{#pm}#pi^{#pm}#pi^{#pm}, Coulomb Uncorrected","p");
1757   //legend4->Draw("same");
1758   //Unity->Draw("same");
1759
1760   //Specif_System->Draw("same");
1761   //Specif_KT3->Draw("same");
1762   //Specif_FSI->Draw("same");
1763   //Specif_Disclaimer->Draw("same");
1764
1765
1766   /*
1767   gPad->SetGridx(0); gPad->SetGridy(0);
1768   int binLow=20;
1769   int binHigh=50;
1770   //TwoFrac=0.9; OneFrac=pow(TwoFrac,0.5); ThreeFrac=pow(TwoFrac,1.5);
1771   TH1D *K0sProjection = (TH1D*)(Three_3d[0][0][1][SCBin][0]->ProjectionY("K0sProjection",binLow,binHigh,binLow,binHigh));
1772   TH1D *K0sProjection_term1 = (TH1D*)(Three_3d[0][0][1][SCBin][0]->ProjectionY("K0sProjection_term1",binLow,binHigh,binLow,binHigh));
1773   TH1D *K0sProjection_term2 = (TH1D*)(Three_3d[0][0][1][SCBin][1]->ProjectionY("K0sProjection_term2",binLow,binHigh,binLow,binHigh));
1774   TH1D *K0sProjection_term3 = (TH1D*)(Three_3d[0][0][1][SCBin][2]->ProjectionY("K0sProjection_term3",binLow,binHigh,binLow,binHigh));
1775   TH1D *K0sProjection_term4 = (TH1D*)(Three_3d[0][0][1][SCBin][3]->ProjectionY("K0sProjection_term4",binLow,binHigh,binLow,binHigh));
1776   TH1D *K0sProjection_term5 = (TH1D*)(Three_3d[0][0][1][SCBin][4]->ProjectionY("K0sProjection_term5",binLow,binHigh,binLow,binHigh));
1777   K0sProjection_term1->Add(K0sProjection_term5, -( pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2) ));
1778   K0sProjection_term1->Add(K0sProjection_term2, -(1-OneFrac));
1779   K0sProjection_term1->Add(K0sProjection_term3, -(1-OneFrac));
1780   K0sProjection_term1->Add(K0sProjection_term4, -(1-OneFrac));
1781   K0sProjection_term1->Add(K0sProjection_term5, (1-OneFrac)*3*(1-TwoFrac));
1782   K0sProjection_term1->Scale(1/ThreeFrac);
1783   //
1784   K0sProjection_term1->Add(K0sProjection_term2, -1/TwoFrac);
1785   K0sProjection_term1->Add(K0sProjection_term3, -1/TwoFrac);
1786   K0sProjection_term1->Add(K0sProjection_term4, -1/TwoFrac);
1787   K0sProjection_term1->Add(K0sProjection_term5, 3*(1-TwoFrac)/TwoFrac);
1788   K0sProjection_term1->Add(K0sProjection_term5, 3);
1789   //
1790   K0sProjection->Divide(K0sProjection_term5);
1791   K0sProjection_term1->Divide(K0sProjection_term5);
1792   K0sProjection_term1->GetXaxis()->SetRangeUser(0,.5);
1793   K0sProjection_term1->SetMarkerStyle(24);
1794   K0sProjection_term1->SetMarkerColor(4);
1795   K0sProjection->SetMarkerStyle(20);
1796   K0sProjection->SetMarkerColor(4);
1797   K0sProjection->SetMinimum(0.92);
1798   K0sProjection->GetXaxis()->SetTitle("#font[12]{q_{13}^{#pm#mp}} (GeV/#font[12]{c})");
1799   K0sProjection->GetYaxis()->SetTitle("#font[12]{C_{3}} or #font[12]{#bf{c}_{3}}");
1800   K0sProjection->GetXaxis()->SetTitleOffset(1.2); K0sProjection->GetYaxis()->SetTitleOffset(1.3);
1801
1802   K0sProjection->Draw();
1803   // Sys errors
1804   TH1D *Sys_K0sProjection=(TH1D*)K0sProjection->Clone();
1805   TH1D *Sys_K0sProjection_term1=(TH1D*)K0sProjection_term1->Clone();
1806   Sys_K0sProjection->SetMarkerSize(0); Sys_K0sProjection_term1->SetMarkerSize(0);
1807   Sys_K0sProjection->SetFillColor(kBlue-10); Sys_K0sProjection_term1->SetFillColor(kBlue-10);
1808   Sys_K0sProjection->SetFillStyle(1000); Sys_K0sProjection_term1->SetFillStyle(1000);
1809   for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) {
1810     Sys_K0sProjection->SetBinError(i+1, 0.01 * Sys_K0sProjection->GetBinContent(i+1));
1811     Sys_K0sProjection_term1->SetBinError(i+1, 0.01 * Sys_K0sProjection_term1->GetBinContent(i+1));
1812   }
1813   for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection->GetBinContent(i+1)<<", ";
1814   cout<<endl;
1815   for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection->GetBinError(i+1)<<", ";
1816   cout<<endl;
1817   for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection_term1->GetBinContent(i+1)<<", ";
1818   cout<<endl;
1819   for(int i=0; i<Sys_K0sProjection->GetNbinsX(); i++) cout<<K0sProjection_term1->GetBinError(i+1)<<", ";
1820   cout<<endl;
1821
1822   Sys_K0sProjection->Draw("E2 same");
1823   Sys_K0sProjection_term1->Draw("E2 same");
1824   K0sProjection->Draw("same");
1825   K0sProjection_term1->Draw("same");
1826
1827   legend4->AddEntry(K0sProjection,"#font[12]{C_{3}^{#pm#pm#mp}} projection","p");
1828   legend4->AddEntry(K0sProjection_term1,"#font[12]{#bf{c}_{3}^{#pm#pm#mp}} projection","p");
1829   legend4->Draw("same");
1830
1831
1832   */
1833
1834
1835   //////////////////////////////////////////////////////////////////////////////////////
1836
1837
1838
1839
1840
1841
1842   if(SaveToFile){
1843     TString *savefilename = new TString("OutFiles/OutFile");
1844     if(CollisionType==0) savefilename->Append("PbPb");
1845     else if(CollisionType==1) savefilename->Append("pPb");
1846     else savefilename->Append("pp");
1847     //
1848     if(MCcase) savefilename->Append("MonteCarlo");
1849     //
1850     if(SameCharge) savefilename->Append("SC");
1851     else savefilename->Append("MC");
1852     //
1853     if(CHARGE==1) savefilename->Append("Pos");
1854     else savefilename->Append("Neg");
1855     //
1856
1857     savefilename->Append("Kt3_");
1858     *savefilename += EDbin+1;
1859
1860     savefilename->Append("_M");
1861     *savefilename += Mbin;
1862     savefilename->Append(".root");
1863     TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
1864     can1->Write("can1");
1865     C2_ss->Write("C2_ss");
1866     C2_os->Write("C2_os");
1867     fitC2ss_h->Write("fitC2ss_h");
1868     fitC2ss_Base->Write("fitC2ss_Base");
1869     fitC2ss_Expan->Write("fitC2ss_Expan");
1870     Cterm1->Write("C3");
1871     c3_hist->Write("c3");
1872     Combinatorics_1d->Write("Combinatorics_1d");
1873     c3Fit1D_Base->Write("c3Fit1D_Base");
1874     c3Fit1D_Expan->Write("c3Fit1D_Expan");
1875     c3Fit1D_ExpanSpline->Write("c3Fit1D_ExpanSpline");
1876     gr_c3Spline->Write("gr_c3Spline");
1877     MyMinuit_c3.Write("MyMinuit_c3");
1878     //
1879     savefile->Close();
1880   }
1881
1882
1883 }
1884 void ReadCoulCorrections(int index){
1885   ///////////////////////
1886   TString *fname2 = new TString("KFile.root");
1887
1888   TFile *File=new TFile(fname2->Data(),"READ");
1889   if(index>=12) cout<<"FSI index not acceptable in 2-particle Coulomb read"<<endl;
1890   TString *nameSS = new TString("K2ss_");
1891   TString *nameOS = new TString("K2os_");
1892   *nameSS += index;
1893   *nameOS += index;
1894   TString *nameSS_2 = new TString("K2ss_");
1895   TString *nameOS_2 = new TString("K2os_");
1896   *nameSS_2 += index+1; *nameOS_2 += index+1;
1897   /*if(index<9) {*nameSS_2 += index+1; *nameOS_2 += index+1;}
1898     else {*nameSS_2 += index; *nameOS_2 += index;}*/
1899   TH1D *temp_ss = (TH1D*)File->Get(nameSS->Data());
1900   TH1D *temp_os = (TH1D*)File->Get(nameOS->Data());
1901   TH1D *temp_ss_2 = (TH1D*)File->Get(nameSS_2->Data());
1902   TH1D *temp_os_2 = (TH1D*)File->Get(nameOS_2->Data());
1903
1904   CoulCorr2SS = (TH1D*)temp_ss->Clone();
1905   CoulCorr2OS = (TH1D*)temp_os->Clone();
1906   CoulCorr2SS_2 = (TH1D*)temp_ss_2->Clone();
1907   CoulCorr2OS_2 = (TH1D*)temp_os_2->Clone();
1908
1909   //
1910   CoulCorr2SS->SetDirectory(0);
1911   CoulCorr2OS->SetDirectory(0);
1912   CoulCorr2SS_2->SetDirectory(0);
1913   CoulCorr2OS_2->SetDirectory(0);
1914
1915   File->Close();
1916 }
1917 double CoulCorr2(int chargeproduct, double Q2){// returns K2
1918   int indexL=0;
1919   int indexH=0;
1920   double slope=0;
1921   double value1=1.0, value2=1.0;
1922
1923   indexL = int(fabs(Q2 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
1924   indexH = indexL+1;
1925
1926
1927   if(indexH >= CoulCorr2SS->GetNbinsX()) {
1928     return (Gamov(chargeproduct, Q2) * CoulCorr2SS->GetBinContent(CoulCorr2SS->GetNbinsX()-1) / Gamov(chargeproduct, CoulCorr2SS->GetXaxis()->GetBinCenter(CoulCorr2SS->GetNbinsX()-1)));
1929   }
1930   // Low index value
1931   if(chargeproduct==1){
1932     slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
1933     slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
1934     value1 = slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1);
1935   }else {
1936     slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1));
1937     slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1));
1938     value1 = slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
1939   }
1940   // High index value
1941   if(chargeproduct==1){
1942     slope = (CoulCorr2SS_2->GetBinContent(indexL+1) - CoulCorr2SS_2->GetBinContent(indexH+1));
1943     slope /= (CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexH+1));
1944     value2 = slope*(Q2 - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS_2->GetBinContent(indexL+1);
1945   }else {
1946     slope = (CoulCorr2OS_2->GetBinContent(indexL+1) - CoulCorr2OS_2->GetBinContent(indexH+1));
1947     slope /= (CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexH+1));
1948     value2 = slope*(Q2 - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS_2->GetBinContent(indexL+1);
1949   }
1950
1951   if(value1>0 && value2>0){// interpolation only done between 9 and 12 Mbin.
1952     return (value1 + (value2-value1)/(2.));// always use the half-way point.
1953   }else if(value1>0){
1954     return value1;
1955   }else if(value2>0){
1956     return value2;
1957   }else return 1.0;
1958
1959
1960   //
1961   // old way with Gaussian sources for low mult systems
1962   /*if(indexH >= CoulCorr2SS->GetNbinsX()) {
1963     return (Gamov(chargeproduct, Q2) * CoulCorr2SS->GetBinContent(CoulCorr2SS->GetNbinsX()-1) / Gamov(chargeproduct, CoulCorr2SS->GetXaxis()->GetBinCenter(CoulCorr2SS->GetNbinsX()-1)));
1964   }
1965   if(chargeproduct==1){
1966     slope = (CoulCorr2SS->GetBinContent(indexL+1) - CoulCorr2SS->GetBinContent(indexH+1));
1967     slope /= (CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS->GetXaxis()->GetBinCenter(indexH+1));
1968     value1 = slope*(Q2 - CoulCorr2SS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS->GetBinContent(indexL+1);
1969   }else {
1970     slope = (CoulCorr2OS->GetBinContent(indexL+1) - CoulCorr2OS->GetBinContent(indexH+1));
1971     slope /= (CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS->GetXaxis()->GetBinCenter(indexH+1));
1972     value1 = slope*(Q2 - CoulCorr2OS->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS->GetBinContent(indexL+1);
1973   }
1974   if(Mbin_def<=9 || Mbin_def>12) return value1;
1975   else{// interpolate in K factor transition region
1976     if(chargeproduct==1){
1977       slope = (CoulCorr2SS_2->GetBinContent(indexL+1) - CoulCorr2SS_2->GetBinContent(indexH+1));
1978       slope /= (CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexH+1));
1979       value2 = slope*(Q2 - CoulCorr2SS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2SS_2->GetBinContent(indexL+1);
1980     }else {
1981       slope = (CoulCorr2OS_2->GetBinContent(indexL+1) - CoulCorr2OS_2->GetBinContent(indexH+1));
1982       slope /= (CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1) - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexH+1));
1983       value2 = slope*(Q2 - CoulCorr2OS_2->GetXaxis()->GetBinCenter(indexL+1)) + CoulCorr2OS_2->GetBinContent(indexL+1);
1984     }
1985
1986     if(value1>0 && value2>0){// interpolation only done between 9 and 12 Mbin.
1987       return (value1 + (Mbin_def-9)*(value2-value1)/(3.));
1988     }else if(value1>0){
1989       return value1;
1990     }else if(value2>0){
1991       return value2;
1992     }else return 1.0;
1993
1994   }
1995   */
1996   //
1997
1998 }
1999
2000 //----------------------------------------------------------------------
2001
2002
2003 //________________________________________________________________________
2004 void fcnC2_Global(int& npar, double* deriv, double& f, double par[], int flag){
2005
2006   double qinvSS=0;
2007
2008   double Rch=par[3]/FmToGeV;
2009   double CSS=0, COS=0;
2010   double SumChi2=0;
2011   double Expan=0;
2012   NFitPoints_C2global=0;
2013   if(LinkN) par[9]=par[0];// Link N factors
2014
2015   for(int i=1; i<BINRANGE_C2global; i++){
2016
2017     qinvSS = BinCenters[i];
2018     if(qinvSS > Q2Limit) continue;
2019
2020     //
2021     double arg=qinvSS*Rch;
2022     if(Gaussian){// Edgeworth expansion
2023       Expan = 1 + par[5]/(6.*pow(2,1.5)) * (8.*pow(arg,3) - 12.*pow(arg,1));
2024       Expan += par[6]/(24.*pow(2,2)) * (16.*pow(arg,4) -48.*pow(arg,2) + 12);
2025       Expan += par[7]/(120.*pow(2,2.5)) * (32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2026       Expan += par[8]/(720.*pow(2,3)) * (64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2027     }else{
2028       Expan = 1 + par[5] * (arg - 1);
2029       Expan += par[6]/2. * (pow(arg,2) - 4*arg + 2);
2030       Expan += par[7] * 0.;
2031       Expan += par[8] * 0.;
2032     }
2033     //Expan = 1.0;
2034     //
2035
2036     // old way without undilution
2037     /*CSS = 1 + exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2038     CSS *= par[1]*K2SS[i];
2039     CSS += 1-par[1];
2040     CSS *= par[0];*/
2041     // undilution method
2042     //CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2043     //CSS *= par[0];
2044     // undilution method with correct normalization location (was previously applied to C2^QS)
2045     CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2046     CSS *= lambda_PM * K2SS[i];
2047     CSS += 1-lambda_PM;
2048     CSS *= par[0];
2049
2050     //
2051     COS = 1;
2052     COS *= par[1]*K2OS[i];
2053     COS += 1-par[1];
2054     COS *= par[9];// different Norm
2055     //
2056     if(C2ssFitting[i] > 0){
2057       if(IncludeSS) {
2058         //double errorSS = sqrt(pow( (CSS/par[0] - (1-par[1]))/K2SS[i] * K2SS_e[i] * par[0],2) + pow(C2ssFitting_e[i],2));
2059         double errorSS = sqrt(pow(C2ssFitting_e[i],2));// new way with undilution already done in C2ssFitting[i]
2060         SumChi2 += pow((CSS - C2ssFitting[i])/errorSS,2);
2061         NFitPoints_C2global++;
2062       }
2063     }
2064     if(IncludeOS) {
2065       double errorOS = sqrt(pow( (COS/par[9] - (1-par[1]))/K2OS[i] * K2OS_e[i] * par[9],2) + pow(C2osFitting_e[i],2));
2066       SumChi2 += pow((COS - C2osFitting[i])/errorOS,2);
2067       NFitPoints_C2global++;
2068     }
2069   }
2070
2071
2072   f = SumChi2;
2073   Chi2_C2global = f;
2074
2075 }
2076 //________________________________________________________________________
2077 double C2ssFitFunction(Double_t *x, Double_t *par)
2078 {
2079   double Rch=par[3]/FmToGeV;
2080   int qbin = int(fabs(x[0]/BinWidthQ2));
2081
2082   double qinvSS = BinCenters[qbin];
2083
2084   double arg=qinvSS*Rch;
2085   double Expan = 1;
2086   if(Gaussian){// Edgeworth expansion
2087     Expan += par[5]/(6.*pow(2,1.5))*(8.*pow(arg,3) - 12.*pow(arg,1));
2088     Expan += par[6]/(24.*pow(2,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
2089     Expan += par[7]/(120.*pow(2,2.5))*(32.*pow(arg,5) - 160.*pow(arg,3) + 120*arg);
2090     Expan += par[8]/(720.*pow(2,3))*(64.*pow(arg,6) - 480.*pow(arg,4) + 720.*pow(arg,2) - 120);
2091   }else{// Laguerre
2092     Expan = 1 + par[5] * (arg - 1);
2093     Expan += par[6]/2. * (pow(arg,2) - 4*arg + 2);
2094     Expan += par[7] * 0.;
2095     Expan += par[8] * 0.;
2096   }
2097
2098   // old way without undilution
2099   /*double CSS = 1 + exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2100   double K2=1.0;
2101   if(qbin < BINRANGE_C2global) K2=K2SS[qbin];
2102   CSS *= par[1]*K2;
2103   CSS += 1-par[1];
2104   CSS *= par[0];*/
2105   // undilution method
2106   //double CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2107   //CSS *= par[0];
2108   // undilution method with correct normalization location (was previously applied to C2^QS)
2109   double CSS = 1 + par[1]*exp(-pow(Rch*qinvSS,ExpPower))*pow(Expan,2);
2110   CSS *= lambda_PM * K2SS[qbin];
2111   CSS += 1-lambda_PM;
2112   CSS *= par[0];
2113
2114   return CSS;
2115 }
2116 //________________________________________________________________________
2117 double C2osFitFunction(Double_t *x, Double_t *par)
2118 {
2119   if(LinkN) par[9]=par[0];// Link N factors
2120   int qbin = int(fabs(x[0]/BinWidthQ2));
2121
2122   //double qinvOS = BinCenters[qbin];
2123
2124   double COS = 1;
2125   double K2=1.0;
2126   if(qbin < BINRANGE_C2global) K2=K2OS[qbin];
2127   COS *= par[1]*K2;
2128   COS += 1-par[1];
2129   COS *= par[9];
2130   return COS;
2131 }
2132 //__________________________________________________________________________
2133
2134 void fcn_c3(int& npar, double* deriv, double& f, double par[], int flag){
2135
2136   double q12=0, q13=0, q23=0;
2137   double Expan12=0, Expan13=0, Expan23=0;
2138   double C=0;
2139   double Rch=par[2]/FmToGeV;
2140   double SumChi2=0;
2141   //double lnL=0, term1=0, term2=0;
2142   NFitPoints_c3=0;
2143
2144   for(int i=0; i<=BINLIMIT_3; i++){// q12
2145     for(int j=0; j<=BINLIMIT_3; j++){// q13
2146       for(int k=0; k<=BINLIMIT_3; k++){// q23
2147
2148         if(B_3[i][j][k] == 0) continue;
2149         if(A_3[i][j][k] == 0) continue;
2150         if(A_3_e[i][j][k] == 0) continue;
2151
2152         q12 = BinCenters[i];
2153         q13 = BinCenters[j];
2154         q23 = BinCenters[k];
2155         double q3 = sqrt(pow(q12,2)+pow(q13,2)+pow(q23,2));
2156         if(q3 > Q3Limit) continue;
2157         //
2158         double arg12 = q12*Rch;
2159         double arg13 = q13*Rch;
2160         double arg23 = q23*Rch;
2161         if(Gaussian){// Edgeworth expansion
2162           Expan12 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg12,3) - 12*pow(arg12,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg12,4) -48*pow(arg12,2) + 12);
2163           Expan13 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg13,3) - 12*pow(arg13,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg13,4) -48*pow(arg13,2) + 12);
2164           Expan23 = 1 + par[3]/(6.*pow(2,1.5))*(8*pow(arg23,3) - 12*pow(arg23,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg23,4) -48*pow(arg23,2) + 12);
2165         }else{// Laguerre expansion
2166           Expan12 = 1 + par[3]*(arg12 - 1) + par[4]/2.*(pow(arg12,2) - 4*arg12 + 2);
2167           Expan13 = 1 + par[3]*(arg13 - 1) + par[4]/2.*(pow(arg13,2) - 4*arg13 + 2);
2168           Expan23 = 1 + par[3]*(arg23 - 1) + par[4]/2.*(pow(arg23,2) - 4*arg23 + 2);
2169         }
2170         //
2171         C = 1 + par[1]*exp(-(pow(arg12,ExpPower)+pow(arg13,ExpPower)+pow(arg23,ExpPower))/2.)*Expan12*Expan13*Expan23;
2172         C *= par[0];// Norm
2173
2174         double error = pow(A_3_e[i][j][k]/B_3[i][j][k],2);
2175         error += pow(sqrt(B_3[i][j][k])*A_3[i][j][k]/pow(B_3[i][j][k],2),2);
2176         error = sqrt(error);
2177         SumChi2 += pow( (C - A_3[i][j][k]/B_3[i][j][k])/error,2);
2178         //
2179         /*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));
2180         else term1 = 0;
2181         term2 = (A_3[i][j][k]+B_3[i][j][k])/(B_3[i][j][k]*(C+1));
2182
2183         if(term1 > 0.0 && term2 > 0.0){
2184           lnL += A_3[i][j][k]*log(term1) + B_3[i][j][k]*log(term2);
2185         }else if(term1==0 && term2 > 0.0){
2186           lnL += B_3[i][j][k]*log(term2);
2187         }else {cout<<"WARNING -- term1 and term2 are negative"<<endl;}
2188         */
2189         float DegenerateCount=1;
2190         if(i==j && i==k) DegenerateCount=1;
2191         else if(i==j || i==k || j==k) DegenerateCount=3;
2192         else DegenerateCount=6;
2193         NFitPoints_c3 += 1/DegenerateCount;
2194
2195       }
2196     }
2197   }
2198   //f = -2.0*lnL;// log-liklihood minimization
2199   f = SumChi2;// Chi2 minimization
2200   Chi2_c3 = f;
2201
2202 }
2203 //________________________________________________________________________
2204 double Dfitfunction_c3(Double_t *x, Double_t *par)
2205 {
2206   double Rch = par[2]/FmToGeV;
2207   double arg12 = x[0]*Rch;
2208   double arg13 = x[1]*Rch;
2209   double arg23 = x[2]*Rch;
2210   double Expan12 = 1, Expan13 = 1, Expan23 = 1;
2211   if(Gaussian){// Edworth Expansion
2212     Expan12 += par[3]/(6.*pow(2,1.5))*(8*pow(arg12,3) - 12*pow(arg12,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg12,4) -48*pow(arg12,2) + 12);
2213     Expan13 += par[3]/(6.*pow(2,1.5))*(8*pow(arg13,3) - 12*pow(arg13,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg13,4) -48*pow(arg13,2) + 12);
2214     Expan23 += par[3]/(6.*pow(2,1.5))*(8*pow(arg23,3) - 12*pow(arg23,1)) + par[4]/(24.*pow(2,2))*(16*pow(arg23,4) -48*pow(arg23,2) + 12);
2215   }else{// Laguerre Expansion
2216     Expan12 = 1 + par[3]*(arg12 - 1) + par[4]/2.*(pow(arg12,2) - 4*arg12 + 2);
2217     Expan13 = 1 + par[3]*(arg13 - 1) + par[4]/2.*(pow(arg13,2) - 4*arg13 + 2);
2218     Expan23 = 1 + par[3]*(arg23 - 1) + par[4]/2.*(pow(arg23,2) - 4*arg23 + 2);
2219   }
2220
2221   //
2222   double C = 1 + par[1]*exp(-(pow(arg12,ExpPower)+pow(arg13,ExpPower)+pow(arg23,ExpPower))/2.)*Expan12*Expan13*Expan23;
2223   C *= par[0];// Norm
2224
2225   return C;
2226 }
2227 //________________________________________________________________________
2228 double CoulCorrGRS(bool SC, double Q_12, double Q_13, double Q_23){
2229
2230   int index12L = int(fabs(Q_12 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2231   int index12H = index12L+1;
2232   int index13L = int(fabs(Q_13 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2233   int index13H = index13L+1;
2234   int index23L = int(fabs(Q_23 - CoulCorr2SS->GetXaxis()->GetBinWidth(1)/2.)/(CoulCorr2SS->GetXaxis()->GetBinWidth(1)));
2235   int index23H = index23L+1;
2236
2237   if(Tricubic){
2238     // Tricubic Interpolation
2239     double arr[4][4][4]={{{0}}};
2240     for(int x=0; x<4; x++){
2241       for(int y=0; y<4; y++){
2242         for(int z=0; z<4; z++){
2243           if(SC){
2244             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2SS->GetBinContent((index23L)+y)*CoulCorr2SS->GetBinContent((index13L)+z);
2245           }else{
2246             arr[x][y][z] = CoulCorr2SS->GetBinContent((index12L)+x)*CoulCorr2OS->GetBinContent((index23L)+y)*CoulCorr2OS->GetBinContent((index13L)+z);
2247           }
2248
2249         }
2250       }
2251     }
2252     return tricubicInterpolate(arr, Q_12, Q_23, Q_13);
2253   }else{
2254     // Trilinear Interpolation.  See for instance: https://en.wikipedia.org/wiki/Trilinear_interpolation
2255     //
2256     double value1=1.0, value2=1.0;
2257     //
2258     double xd = (Q_12-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
2259     xd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index12H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index12L+1));
2260     double yd = (Q_13-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
2261     yd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index13H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index13L+1));
2262     double zd = (Q_23-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
2263     zd /= (CoulCorr2SS->GetXaxis()->GetBinCenter(index23H+1)-CoulCorr2SS->GetXaxis()->GetBinCenter(index23L+1));
2264     //
2265     if(SC){
2266       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
2267       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
2268       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*(1-xd);
2269       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23L+1)*xd;
2270       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
2271       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13L+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
2272       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*(1-xd);
2273       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2SS->GetBinContent(index13H+1)*CoulCorr2SS->GetBinContent(index23H+1)*xd;
2274       //
2275       double c0 = c00*(1-yd) + c10*yd;
2276       double c1 = c01*(1-yd) + c11*yd;
2277       value1 = (c0*(1-zd) + c1*zd);
2278     }else{
2279       double c00 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
2280       c00 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
2281       double c10 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*(1-xd);
2282       c10 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23L+1)*xd;
2283       double c01 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
2284       c01 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13L+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
2285       double c11 = CoulCorr2SS->GetBinContent(index12L+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*(1-xd);
2286       c11 += CoulCorr2SS->GetBinContent(index12H+1)*CoulCorr2OS->GetBinContent(index13H+1)*CoulCorr2OS->GetBinContent(index23H+1)*xd;
2287       //
2288       double c0 = c00*(1-yd) + c10*yd;
2289       double c1 = c01*(1-yd) + c11*yd;
2290       value1 = (c0*(1-zd) + c1*zd);
2291     }
2292
2293     //if(Mbin_def<=9 || Mbin_def>12) return value1;// old way for v5
2294     //
2295     // interpolation for K factor tansition
2296     if(SC){
2297       double c00 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*(1-xd);
2298       c00 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*xd;
2299       double c10 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*(1-xd);
2300       c10 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23L+1)*xd;
2301       double c01 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*(1-xd);
2302       c01 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13L+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*xd;
2303       double c11 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*(1-xd);
2304       c11 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2SS_2->GetBinContent(index13H+1)*CoulCorr2SS_2->GetBinContent(index23H+1)*xd;
2305       //
2306       double c0 = c00*(1-yd) + c10*yd;
2307       double c1 = c01*(1-yd) + c11*yd;
2308       value2 = (c0*(1-zd) + c1*zd);
2309     }else{
2310       double c00 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*(1-xd);
2311       c00 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*xd;
2312       double c10 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*(1-xd);
2313       c10 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23L+1)*xd;
2314       double c01 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*(1-xd);
2315       c01 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13L+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*xd;
2316       double c11 = CoulCorr2SS_2->GetBinContent(index12L+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*(1-xd);
2317       c11 += CoulCorr2SS_2->GetBinContent(index12H+1)*CoulCorr2OS_2->GetBinContent(index13H+1)*CoulCorr2OS_2->GetBinContent(index23H+1)*xd;
2318       //
2319       double c0 = c00*(1-yd) + c10*yd;
2320       double c1 = c01*(1-yd) + c11*yd;
2321       value2 = (c0*(1-zd) + c1*zd);
2322     }
2323
2324     if(value1>0 && value2>0){
2325       return (value1 + (value2-value1)/(2.));// always take the half-way point
2326       //return (value1 + (Mbin_def-9)*(value2-value1)/(3.));// old way for v5
2327     }else if(value1>0){
2328       return value1;
2329     }else if(value2>0){
2330       return value2;
2331     }else return 1.0;
2332
2333   }
2334
2335 }
2336
2337 double Gamov(int chargeProduct, double qinv){
2338
2339   double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
2340
2341   return arg/(exp(arg)-1);
2342 }
2343 void ReadMomResFile(int Mbin){
2344
2345   int MRCindex_L=0, MRCindex_H=0;
2346   float MRCindex_weight=0;
2347   if(Mbin<=2) {MRCindex_L=0; MRCindex_H=1; MRCindex_weight = Mbin/2.;}
2348   else if(Mbin<=6) {MRCindex_L=1; MRCindex_H=2; MRCindex_weight = fabs(Mbin-3)/3.;}
2349   else if(Mbin<=11) {MRCindex_L=2; MRCindex_H=3; MRCindex_weight = fabs(Mbin-7)/4.;}
2350   else {MRCindex_L=3; MRCindex_H=4; MRCindex_weight = fabs(Mbin-12)/7.;}
2351
2352
2353   TH1D *temp_L[2][5];
2354   TH1D *temp_H[2][5];
2355   TH1D *tempC2_L[2];
2356   TH1D *tempC2_H[2];
2357   TString *momresfilenameL = new TString("MomResFile_M");
2358   *momresfilenameL += MRCindex_L;
2359   momresfilenameL->Append(".root");
2360   TFile *MomResFileL = new TFile(momresfilenameL->Data(),"READ");
2361   TString *names1D_L[2][5];// SC/MC, term#
2362   for(int ChProd=0; ChProd<2; ChProd++){
2363     for(int term=0; term<5; term++){
2364       //
2365       if(ChProd==0) {names1D_L[ChProd][term] = new TString("MRC3_SC_term");}
2366       else {names1D_L[ChProd][term] = new TString("MRC3_MC_term");}
2367       *names1D_L[ChProd][term] += term+1;
2368       temp_L[ChProd][term] = (TH1D*)MomResFileL->Get(names1D_L[ChProd][term]->Data());
2369       temp_L[ChProd][term]->SetDirectory(0);
2370
2371     }
2372     TString *C2MRCname=new TString("MomResHisto_");
2373     if(ChProd==0) C2MRCname->Append("pp");
2374     else C2MRCname->Append("mp");
2375     tempC2_L[ChProd]=(TH1D*)(((TH2D*)(MomResFileL->Get(C2MRCname->Data())))->ProjectionY("C2MRCproL",MRC2index,MRC2index));
2376     tempC2_L[ChProd]->SetDirectory(0);
2377   }
2378
2379   //
2380   TString *momresfilenameH = new TString("MomResFile_M");
2381   *momresfilenameH += MRCindex_H;
2382   momresfilenameH->Append(".root");
2383   TFile *MomResFileH = new TFile(momresfilenameH->Data(),"READ");
2384   TString *names1D_H[2][5];// SC/MC, term#
2385   for(int ChProd=0; ChProd<2; ChProd++){
2386     for(int term=0; term<5; term++){
2387       //
2388       if(ChProd==0) {names1D_H[ChProd][term] = new TString("MRC3_SC_term");}
2389       else {names1D_H[ChProd][term] = new TString("MRC3_MC_term");}
2390       *names1D_H[ChProd][term] += term+1;
2391       temp_H[ChProd][term] = (TH1D*)MomResFileH->Get(names1D_H[ChProd][term]->Data());
2392       temp_H[ChProd][term]->SetDirectory(0);
2393       MomRes1d[ChProd][term] = (TH1D*)MomResFileH->Get(names1D_H[ChProd][term]->Data());
2394       MomRes1d[ChProd][term]->SetDirectory(0);
2395     }
2396     TString *C2MRCname=new TString("MomResHisto_");
2397     if(ChProd==0) C2MRCname->Append("pp");
2398     else C2MRCname->Append("mp");
2399     tempC2_H[ChProd]=(TH1D*)(((TH2D*)(MomResFileH->Get(C2MRCname->Data())))->ProjectionY("C2MRCproH",MRC2index,MRC2index));
2400     tempC2_H[ChProd]->SetDirectory(0);
2401     TString *name=new TString("MomResC2_");
2402     *name += ChProd;
2403     MomResC2[ChProd] = (TH1D*)(((TH2D*)(MomResFileH->Get(C2MRCname->Data())))->ProjectionY(name->Data(),MRC2index,MRC2index));
2404     MomResC2[ChProd]->SetDirectory(0);
2405   }
2406
2407   for(int ChProd=0; ChProd<2; ChProd++){
2408     // C3 MRC
2409     for(int term=0; term<5; term++){
2410       for(int bin=1; bin<=temp_H[ChProd][term]->GetNbinsX(); bin++){
2411         double value=1;
2412         if(temp_L[ChProd][term]->GetBinContent(bin)>0 && temp_H[ChProd][term]->GetBinContent(bin)>0){// both have entries
2413           value = temp_L[ChProd][term]->GetBinContent(bin) + MRCindex_weight * temp_H[ChProd][term]->GetBinContent(bin);
2414         }else if(temp_L[ChProd][term]->GetBinContent(bin)>0){
2415           value = temp_L[ChProd][term]->GetBinContent(bin);
2416         }else if(temp_H[ChProd][term]->GetBinContent(bin)>0){
2417           value = temp_H[ChProd][term]->GetBinContent(bin);
2418         }else value=1.0;
2419
2420         MomRes1d[ChProd][term]->SetBinContent(bin,value);
2421       }
2422     }
2423     // C2 MRC
2424     for(int bin=1; bin<=tempC2_H[ChProd]->GetNbinsX(); bin++){
2425       double value=1;
2426       if(tempC2_L[ChProd]->GetBinContent(bin)>0 && tempC2_H[ChProd]->GetBinContent(bin)>0){// both have entries
2427         value = tempC2_L[ChProd]->GetBinContent(bin)*(1-MRCindex_weight) + MRCindex_weight * tempC2_H[ChProd]->GetBinContent(bin);
2428       }else if(tempC2_L[ChProd]->GetBinContent(bin)>0){
2429         value = tempC2_L[ChProd]->GetBinContent(bin);
2430       }else if(tempC2_H[ChProd]->GetBinContent(bin)>0){
2431         value = tempC2_H[ChProd]->GetBinContent(bin);
2432       }else value=1.0;
2433       MomResC2[ChProd]->SetBinContent(bin,value);
2434     }
2435   }
2436
2437   for(int ChProd=0; ChProd<2; ChProd++){
2438     for(int term=0; term<5; term++){
2439       for(int i=0; i<MomRes1d[0][0]->GetNbinsX(); i++){
2440         if(SourceType==0 && Mbin<10){
2441           if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.1) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0);
2442         }else if(SourceType==0 && Mbin>=10){
2443           if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.2) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0);
2444         }else{
2445           if(MomRes1d[ChProd][term]->GetXaxis()->GetBinCenter(i+1) > 0.5) MomRes1d[ChProd][term]->SetBinContent(i+1, 1.0);
2446         }
2447
2448       }
2449     }
2450
2451   }
2452
2453   MomResFileL->Close();
2454   MomResFileH->Close();
2455
2456 }
2457 double cubicInterpolate (double p[4], double x) {
2458   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
2459 }
2460 double bicubicInterpolate (double p[4][4], double x, double y) {
2461         double arr[4];
2462         arr[0] = cubicInterpolate(p[0], y);
2463         arr[1] = cubicInterpolate(p[1], y);
2464         arr[2] = cubicInterpolate(p[2], y);
2465         arr[3] = cubicInterpolate(p[3], y);
2466         return cubicInterpolate(arr, x);
2467 }
2468 double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
2469         double arr[4];
2470         arr[0] = bicubicInterpolate(p[0], y, z);
2471         arr[1] = bicubicInterpolate(p[1], y, z);
2472         arr[2] = bicubicInterpolate(p[2], y, z);
2473         arr[3] = bicubicInterpolate(p[3], y, z);
2474         return cubicInterpolate(arr, x);
2475 }
2476 //____________________________________________________________________________________________________
2477 void ReadMuonCorrections(int ct, int mbin){
2478   TString *name = new TString("MuonCorrection_");
2479   if(ct==0){// PbPb
2480     if(mbin<=1) *name += 0;
2481     else if(mbin<=3) *name += 1;
2482     else if(mbin<=5) *name += 2;
2483     else if(mbin<=7) *name += 3;
2484     else if(mbin<10) *name += 4;
2485     else *name += 5;
2486   }else *name += 6;
2487   name->Append(".root");
2488   TFile *muonfile=new TFile(name->Data(),"READ");
2489   C2muonCorrection_SC = (TH1D*)muonfile->Get("C2muonCorrection_SC");
2490   C2muonCorrection_MC = (TH1D*)muonfile->Get("C2muonCorrection_MC");
2491   if(SameCharge_def) C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_SC");
2492   else C3muonCorrection = (TH1D*)muonfile->Get("C3muonCorrection_MC");
2493   //
2494   C2muonCorrection_SC->SetDirectory(0);
2495   C2muonCorrection_MC->SetDirectory(0);
2496   C3muonCorrection->SetDirectory(0);
2497   //
2498   muonfile->Close();
2499 }
2500 //____________________________________________________________________________________________________
2501 void DrawALICELogo(Bool_t prel, Float_t x1, Float_t y1, Float_t x2, Float_t y2)
2502 {
2503   // revision on July 24th or 25th 2012
2504   // correct for aspect ratio of figure plus aspect ratio of pad (coordinates are NDC!)
2505   x2 = x1 + (y2 - y1) * (466. / 523) * gPad->GetWh() * gPad->GetHNDC() / (gPad->GetWNDC() * gPad->GetWw());
2506   // Printf("%f %f %f %f", x1, x2, y1, y2);
2507
2508   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo", x1, y1, x2, y2);
2509   myPadLogo->SetLeftMargin(0);
2510   myPadLogo->SetTopMargin(0);
2511   myPadLogo->SetRightMargin(0);
2512   myPadLogo->SetBottomMargin(0);
2513   myPadLogo->Draw();
2514   myPadLogo->cd();
2515   TASImage *myAliceLogo = new TASImage((prel) ? "~/Pictures/2011-Nov-24-ALICE_PRELIMARY_logo_BLACK_small_usage_design.eps" : "~/Pictures/2011-Nov-24-ALICE_PERFORMANCE_logo_BLACK_small_usage_design.eps");
2516   myAliceLogo->Draw();
2517 }