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