PWGCF/FEMTOSCOPY/macros/Plot_FourPion.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <assert.h>
   7
   8 #include "TVector2.h"
   9 #include "TFile.h"
  10 #include "TString.h"
  11 #include "TF1.h"
  12 #include "TF3.h"
  13 #include "TH1.h"
  14 #include "TH2.h"
  15 #include "TH3.h"
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TMath.h"
  19 #include "TText.h"
  20 #include "TRandom3.h"
  21 #include "TArray.h"
  22 #include "TLegend.h"
  23 #include "TStyle.h"
  24 #include "TMinuit.h"
  25 #include "TCanvas.h"
  26 #include "TLatex.h"
  27 #include "TPaveStats.h"
  28 #include "TASImage.h"
  29
  30 #define BohrR 1963.6885 // Bohr Radius for pions
  31 #define FmToGeV 0.19733 // conversion to fm
  32 #define PI 3.1415926
  33 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  34
  35 using namespace std;
  36 // 0(standard), 1(fcSq=0.65), 2(fcSq=0.75), 3(B minus), 4(B plus),
  37 // 5(Linear instead of Cubic,also fc^2=0.75), 6(10h), 7(MRC 10% increase), 8(MuonCorr 92% pion-pair purity)
  38 int FileSetting=0;
  39 //
  40 bool MCcase_def=0;// MC data?
  41 int CollisionType=0;// PbPb, pPb, pp
  42 //
  43 int Mbin_def=0;// 0-9: centrality bin in widths of 5%
  44 bool SameCharge_def=1;// same-charge?
  45 int CHARGE_def=1;// -1 or +1: + or - pions for same-charge case, --+ or -++,  ---+ or -+++
  46 int MixedCharge4pionType_def = 2;// 1(---+) or 2(--++)
  47 //
  48 int EDbin_def=0;// 0 or 1: Kt3 bin
  49 int TPNbin=0;// TPN bin for r3 and r4
  50 int Gbin = int( (0) /2. ) + 5;// +5 (Rcoh=0), +55 (Rcoh=Rch)
  51 int c3FitType = 1;// EW(1), LG(2)
  52 int Ktbin_def=1;// 1(0.2-0.3),..., 6(0.7-0.8), 10 = Full Range
  53 //
  54 bool MRC=1;// Momentum Resolution Corrections?
  55 bool MuonCorrection=1;// correct for Muons?
  56 bool FSICorrection=1;// correct For Final-State-Interactions
  57 bool InterpCorrection=1;// correct for finite bin interpolation
  58 //
  59 int f_choice=0;// 0(Core/Halo), 1(40fm), 2(70fm), 3(100fm)
  60 //
  61 //
  62 bool SaveToFile_def=kFALSE;// Save outputs to file?
  63 bool GeneratedSignal=kFALSE;// Was the QS+FSI signal generated in MC?
  64 //
  65 //
  66 //
  67 //
  68
  69 int fFSIindex=0;
  70 int Ktlowbin;
  71 int Kthighbin;
  72 float TwoFrac;// Lambda parameter
  73 float OneFrac;// Lambda^{1/2}
  74 float ThreeFrac;// Lambda^{3/2}
  75 float FourFrac;// lambda^{4/2}
  76 double Qcut_pp = 0.6;// 0.6
  77 double Qcut_PbPb = 0.1;
  78 double NormQcutLow_pp = 1.0;
  79 double NormQcutHigh_pp = 1.5;
  80 double NormQcutLow_PbPb = .15;
  81 double NormQcutHigh_PbPb = .2;// was .175
  82
  83 int ThreeParticleRebin;
  84 int FourParticleRebin;
  85
  86 int RbinMRC;
  87
  88 TH1D *fFSIss[12];
  89 TH1D *fFSIos[12];
  90
  91 void SetFSICorrelations();
  92 void SetFSIindex(Float_t);
  93 Float_t FSICorrelation(Int_t, Int_t, Float_t);
  94 void SetMuonCorrections();
  95 void SetMomResCorrections();
  96 double Gamov(int, double);
  97 double cubicInterpolate(double[4], double);
  98 double bicubicInterpolate(double[4][4], double, double);
  99 double tricubicInterpolate(double[4][4][4], double, double, double);
 100 //
 101 float fact(float);
 102
 103
 104 TH1D *MRC_SC_2[2];
 105 TH1D *MRC_MC_2[2];
 106 TH1D *MRC_SC_3[3];
 107 TH1D *MRC_MC_3[4];
 108 TH1D *MRC_SC_4[5];
 109 TH1D *MRC_MC1_4[6];
 110 TH1D *MRC_MC2_4[6];
 111
 112
 113 double AvgQinvSS[30];
 114 double AvgQinvOS[30];
 115 double BinCentersQ4[20];
 116
 117 //
 118 TH1D *C2muonCorrectionSC;
 119 TH1D *C2muonCorrectionMC;
 120 TH1D *WeightmuonCorrection;
 121 TH1D *C3muonCorrectionSC[2];
 122 TH1D *C3muonCorrectionMC[3];
 123 TH1D *C4muonCorrectionSC[4];
 124 TH1D *C4muonCorrectionMC1[5];
 125 TH1D *C4muonCorrectionMC2[5];
 126
 127
 128 void Plot_FourPion(bool SaveToFile=SaveToFile_def, bool MCcase=MCcase_def, bool SameCharge=SameCharge_def, int MixedCharge4pionType=MixedCharge4pionType_def, int EDbin=EDbin_def, int CHARGE=CHARGE_def, int Mbin=Mbin_def, int Ktbin=Ktbin_def){
 129
 130   EDbin_def=EDbin;
 131   SaveToFile_def=SaveToFile;
 132   MCcase_def=MCcase;
 133   CHARGE_def=CHARGE;
 134   SameCharge_def=SameCharge;// 3-pion same-charge
 135   MixedCharge4pionType_def=MixedCharge4pionType;
 136   Mbin_def=Mbin;
 137   Ktbin_def=Ktbin;
 138   //
 139   Ktlowbin=(Ktbin)*2+3;// kt bins are 0.5 GeV/c wide (0-0.5, 0.5-1.0 ...)
 140   Kthighbin=(Ktbin)*2+4;
 141
 142   //
 143   if(FileSetting==1) TwoFrac=0.65;
 144   else if(FileSetting==2 || FileSetting==5) TwoFrac=0.75;
 145   else TwoFrac=0.7;
 146
 147   OneFrac = sqrt(TwoFrac);
 148   ThreeFrac = pow(TwoFrac, 1.5);
 149   FourFrac = pow(TwoFrac, 2.);
 150
 151   //// Core/Halo, 40fm, 70fm, 100fm
 152   float ThermShift_f33[4]={0., 0.06933, 0.01637, 0.006326};
 153   float ThermShift_f32[4]={0., -0.0185, -0.005301, -0.002286};
 154   float ThermShift_f31[4]={0., -0.01382, -0.0004682, 0.0005337};
 155   float f33prime = ThreeFrac;
 156   float f32prime = TwoFrac*(1-OneFrac);
 157   float f31prime = pow(1-OneFrac,3) + 3*OneFrac*pow(1-OneFrac,2);
 158   f33prime += ThermShift_f33[f_choice];
 159   f32prime += ThermShift_f32[f_choice];
 160   f31prime += ThermShift_f31[f_choice];
 161   float f33 = f33prime;
 162   float f32 = f32prime/TwoFrac;
 163   float f31 = f31prime - 3*((1-TwoFrac)*(1-OneFrac) + ThermShift_f32[f_choice]*(1-TwoFrac)/TwoFrac);
 164   //cout<<f33 + 3*f32 + f31<<endl;
 165
 166   //// Core/Halo, 40fm, 70fm, 100fm
 167   float ThermShift_f44[4]={0., 0.08741, 0.005284, -0.01064};
 168   float ThermShift_f43[4]={0., -0.01126, -0.001486, 0.001686};
 169   float ThermShift_f42[4]={0., -0.006466, -7.683e-05, 0.0004572};
 170   float ThermShift_f41[4]={0., -0.003556, 0.00112, 0.00115};
 171   float f44prime = FourFrac;
 172   float f43prime = ThreeFrac*(1-OneFrac);
 173   float f42prime = TwoFrac*pow(1-OneFrac,2);
 174   float f41prime = pow(1-OneFrac,4) + 4*OneFrac*pow(1-OneFrac,3);
 175   f44prime += ThermShift_f44[f_choice];
 176   f43prime += ThermShift_f43[f_choice];
 177   f42prime += ThermShift_f42[f_choice];
 178   f41prime += ThermShift_f41[f_choice];
 179   float f44 = f44prime;
 180   float f43 = f43prime/f33prime;
 181   float f42 = f42prime/TwoFrac;
 182   f42 -= 2*f43prime*f32prime/f33prime/TwoFrac;
 183   float f41 = f41prime;
 184   f41 -= 4*f43prime*f31prime/f33prime;
 185   f41 -= 6*f42prime*(1-TwoFrac)/TwoFrac;
 186   f41 += 12*f43prime/f33prime*f32prime/TwoFrac*(1-TwoFrac);
 187   //cout<<f44 + 4*f43 + 6*f42 + f41<<endl;
 188   //cout<<f44<<"  "<<f43<<"  "<<f42<<"  "<<f41<<endl;
 189   //
 190   // Core/Halo, 40fm, 70fm, 100fm
 191   //float TherminatorMod_f33[4]={1., 1.123, 1.022, 1.008};// 1.008
 192   //float TherminatorMod_f32[4]={1., 0.844, 0.933, 0.967};// .967
 193   //float TherminatorMod_f31[4]={1., 0.828, 0.982, 1.028};// 1.028
 194   /*float TherminatorMod_f44[4]={1., 1.188, 1.008, 0.985};// .985
 195   float TherminatorMod_f43[4]={1., 0.885, 0.979, 1.026};// 1.026
 196   float TherminatorMod_f42[4]={1., 0.687, 0.99, 1.08};// 1.08
 197   float TherminatorMod_f41[4]={1., 0.806, 1.33, 1.548};//1.548
 198   float f44 = FourFrac;
 199   float f43 = (1-OneFrac);
 200   float f42 = -pow(1-OneFrac,2);
 201   float f41 = -3*pow(1-OneFrac,4) - 8*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2);
 202   f44 *= TherminatorMod_f44[f_choice];
 203   f43 *= TherminatorMod_f43[f_choice];
 204   f42 *= TherminatorMod_f42[f_choice];
 205   f41 *= TherminatorMod_f41[f_choice];*/
 206   //f44prime *= TherminatorMod_f44[f_choice];
 207   //f43prime *= TherminatorMod_f43[f_choice];
 208   //f42prime *= TherminatorMod_f42[f_choice];
 209   //f41prime *= TherminatorMod_f41[f_choice];
 210   //float f44 = f44prime;
 211   //float f43 = f43prime/ThreeFrac;
 212   //float f42 = f42prime/TwoFrac - 2*pow(1-OneFrac,2);
 213   //float f41 = f41prime - 4*pow(1-OneFrac,4) - 12*OneFrac*pow(1-OneFrac,3) + 6*(1-TwoFrac)*pow(1-OneFrac,2);
 214
 215   cout<<"Mbin = "<<Mbin<<"   KT3 = "<<EDbin<<"   SameCharge = "<<SameCharge<<endl;
 216   if(!SameCharge) cout<<"4-pion MixedCharge type = "<<MixedCharge4pionType<<endl;
 217
 218   if(CollisionType==0){
 219     if(Mbin==0) {RbinMRC=10;}
 220     else if(Mbin==1) {RbinMRC=9;}
 221     else if(Mbin<=3) {RbinMRC=8;}
 222     else if(Mbin<=5) {RbinMRC=7;}
 223     else {RbinMRC=6;}
 224   }else{
 225     RbinMRC=2;
 226   }
 227
 228   if(MCcase) MuonCorrection=kFALSE;
 229
 230   if(CollisionType==0){
 231     ThreeParticleRebin=2;
 232     FourParticleRebin=3;
 233   }else{
 234     ThreeParticleRebin=2;
 235     FourParticleRebin=6;
 236   }
 237   float Cutoff_FullWeight_Q3[10]={0.065, 0.065, 0.075, 0.075, 0.075, 0.075, 0.085, 0.085, 0.085, 0.085};
 238   float Cutoff_FullWeight_Q4[10]={0.115, 0.115, 0.115, 0.130, 0.130, 0.130, 0.145, 0.145, 0.145, 0.145};
 239
 240
 241   //
 242   // avg Qinv within the 5 MeV bins (0-5, 5-10,...) for true bin mean values.  From unsmeared HIJING 0-5% with input signal
 243   double AvgQinvSS_temp[30]={0.00421164, 0.00796583, 0.0128473, 0.0178262, 0.0228416, 0.0276507, 0.0325368, 0.0376114, 0.0424707, 0.0476274, 0.0526015, 0.0575645, 0.0625478, 0.0675416, 0.0725359, 0.0775219, 0.0825521, 0.0875211, 0.0925303, 0.0975319, 0.102544, 0.10753, 0.112499, 0.117545, 0.122522, 0.127522, 0.132499, 0.137514, 0.142495, 0.147521};
 244   double AvgQinvOS_temp[30]={0.00352789, 0.00797596, 0.0128895, 0.0177198, 0.0226397, 0.0276331, 0.0326309, 0.0376471, 0.0426217, 0.047567, 0.0525659, 0.0575472, 0.0625886, 0.0675261, 0.0725543, 0.077564, 0.0825617, 0.0875465, 0.092539, 0.0975348, 0.102529, 0.107527, 0.112506, 0.117531, 0.122536, 0.1275, 0.132508, 0.137508, 0.14251, 0.147511};
 245   for(int i=0; i<30; i++){
 246     AvgQinvSS[i] = AvgQinvSS_temp[i];
 247     AvgQinvOS[i] = AvgQinvOS_temp[i];
 248   }
 249
 250
 251
 252   TH2D *TwoParticle_2d[2][2][2];// ch1,ch2,term
 253   TH1D *TwoParticle[2][2][2];// ch1,ch2,term
 254   TH2D *UnitMult_2d[2][2][2];// ch1,ch2,term
 255   TH1D *UnitMult[2][2][2];// ch1,ch2,term
 256   double norm_2[2]={0};
 257   double norm_2_UM[2]={0};
 258   //
 259   TH1D *ThreeParticle[2][2][2][5];// ch1,ch2,ch3,term
 260   TProfile *K3avg[2][2][2][4];
 261   TH2D *TPFullWeight_ThreeParticle_2D[2];// charge
 262   TH2D *TPNegFullWeight_ThreeParticle_2D[2];// charge
 263   TH1D *TPN_ThreeParticle[2];// charge
 264   TH1D *TPFullWeight_ThreeParticle[2];// charge
 265   TH1D *TPNegFullWeight_ThreeParticle[2];// charge
 266   double norm_3[5]={0};
 267   //
 268   TH1D *FourParticle[2][2][2][2][13];// ch1,ch2,ch3,ch4,term
 269   TProfile *K4avg[2][2][2][2][12];
 270   TH2D *TPFullWeight_FourParticle_2D[2];// charge
 271   TH2D *TPNegFullWeight_FourParticle_2D[2];// charge
 272   TH1D *TPFullWeight_FourParticle[2];// charge
 273   TH1D *TPNegFullWeight_FourParticle[2];// charge
 274   TH1D *TPN_FourParticle[2];// charge
 275   TH3D *FullBuildFromFits_3D;// charge
 276   TH1D *FullBuildFromFits;// charge
 277   TH3D *PartialBuildFromFits_3D;// charge
 278   TH1D *PartialBuildFromFits;// charge
 279   double norm_4[13]={0};
 280
 281
 282   gStyle->SetOptStat(0);
 283   gStyle->SetOptDate(0);
 284   gStyle->SetOptFit(0111);
 285
 286   TFile *_file0;
 287   if(CollisionType==0){// PbPb
 288     if(MCcase) {
 289       if(Mbin<=1){
 290         //_file0 = new TFile("Results/PDC_12a17a_GeneratedSignal.root","READ");
 291         //_file0 = new TFile("Results/PDC_12a17a_pTSpectrumWeight.root","READ");
 292         _file0 = new TFile("Results/PDC_12a17a_Qweights.root","READ");
 293       }
 294     }else{
 295       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0_FullRunWrongWeightsNoPadRowTTCandMCTTC_RawWeightFile.root","READ");
 296       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p16to0p25_0p25to1p0.root","READ");
 297       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT_0p2to1p0.root","READ");
 298       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_pT0p16to0p25_KT0p35.root","READ");
 299       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_1percentCentEM.root","READ");
 300       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_0p02eta0p045phi_0p03eta0p067phi.root","READ");
 301       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_c3FitBuild.root","READ");
 302       if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_LowNorm_HighNorm.root","READ");
 303       //if(FileSetting==0) _file0 = new TFile("Results/PDC_11h_extendedGweights.root","READ");// Preliminary results
 304       if(FileSetting==5) _file0 = new TFile("Results/PDC_11h_Cubic_Linear.root","READ");
 305       if(FileSetting==1 || FileSetting==2) _file0 = new TFile("Results/PDC_11h_Lam0p65_Lam0p75.root","READ");
 306       if(FileSetting==3) _file0 = new TFile("Results/PDC_11h_Cubic_Linear_Bminus.root","READ");
 307       if(FileSetting==4) _file0 = new TFile("Results/PDC_11h_Cubic_Linear_Bplus.root","READ");
 308       if(FileSetting==6) _file0 = new TFile("Results/PDC_10h_Cubic_Linear.root","READ");
 309       if(FileSetting==7 || FileSetting==8) _file0 = new TFile("Results/PDC_11h_MRC10percIncrease_Muon92percent.root","READ");
 310     }
 311   }else if(CollisionType==1){// pPb
 312     _file0 = new TFile("Results/PDC_13bc_kINT7_LowNorm_HighNorm.root","READ");
 313   }else{// pp
 314     _file0 = new TFile("Results/PDC_10bcde_kMB_LowNorm_HighNorm.root","READ");
 315   }
 316
 317   SetFSIindex(10.);
 318   SetFSICorrelations();
 319   SetMomResCorrections();
 320   SetMuonCorrections();
 321   //
 322   /////////////////////////////////////////////////////////
 323
 324
 325   double NormQcutLow;
 326   double NormQcutHigh;
 327   if(CollisionType==0) {
 328     NormQcutLow = NormQcutLow_PbPb;
 329     NormQcutHigh = NormQcutHigh_PbPb;
 330   }else {
 331     NormQcutLow = NormQcutLow_pp;
 332     NormQcutHigh = NormQcutHigh_pp;
 333   }
 334
 335
 336   TList *MyList;
 337   if(!MCcase){
 338     TDirectoryFile *tdir = (TDirectoryFile*)_file0->Get("PWGCF.outputFourPionAnalysis.root");
 339     if(FileSetting==2 || FileSetting==5 || FileSetting==8 || CollisionType!=0) MyList=(TList*)tdir->Get("FourPionOutput_2");
 340     else MyList=(TList*)tdir->Get("FourPionOutput_1");
 341     //MyList=(TList*)_file0->Get("MyList");
 342   }else{
 343     MyList=(TList*)_file0->Get("MyList");
 344   }
 345   _file0->Close();
 346
 347
 348   TH1D *Events = (TH1D*)MyList->FindObject("fEvents2");
 349   //
 350
 351   cout<<"#Events = "<<Events->Integral(Mbin+1,Mbin+1)<<endl;
 352
 353
 354
 355   ///////////////////////////////////
 356   // Get Histograms
 357
 358   // 2-particle
 359   for(int c1=0; c1<2; c1++){
 360     for(int c2=0; c2<2; c2++){
 361       for(int term=0; term<2; term++){
 362
 363         TString *name2 = new TString("TwoParticle_Charge1_");
 364         *name2 += c1;
 365         name2->Append("_Charge2_");
 366         *name2 += c2;
 367         name2->Append("_M_");
 368         *name2 += Mbin;
 369         name2->Append("_ED_");
 370         *name2 += 0;// EDbin
 371         name2->Append("_Term_");
 372         *name2 += term+1;
 373
 374         if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
 375
 376         TwoParticle_2d[c1][c2][term] = (TH2D*)MyList->FindObject(name2->Data());
 377         TwoParticle_2d[c1][c2][term]->Sumw2();
 378         TString *proname = new TString(name2->Data());
 379         proname->Append("_pro");
 380
 381         if(Ktbin==10) {Ktlowbin=1; Kthighbin=TwoParticle_2d[c1][c2][term]->GetNbinsX();}
 382         TwoParticle[c1][c2][term] = (TH1D*)TwoParticle_2d[c1][c2][term]->ProjectionY(proname->Data(),Ktlowbin,Kthighbin);// bins:5-6 (kt:0.2-0.3)
 383
 384         norm_2[term] = TwoParticle[c1][c2][term]->Integral(TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutLow),TwoParticle[c1][c2][term]->GetXaxis()->FindBin(NormQcutHigh));
 385         //cout<<"2-pion norms  "<<norm_2[term]<<endl;
 386         TwoParticle[c1][c2][term]->Scale(norm_2[0]/norm_2[term]);
 387
 388         TwoParticle[c1][c2][term]->SetMarkerStyle(20);
 389         TwoParticle[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 390         TwoParticle[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
 391         TwoParticle[c1][c2][term]->SetTitle("");
 392         //
 393
 394         TString *nameUM=new TString(name2->Data());
 395         nameUM->Append("_UnitMult");
 396         UnitMult_2d[c1][c2][term] = (TH2D*)MyList->FindObject(nameUM->Data());
 397         UnitMult_2d[c1][c2][term]->Sumw2();
 398         TString *pronameUM = new TString(nameUM->Data());
 399         pronameUM->Append("_pro");
 400         UnitMult[c1][c2][term] = (TH1D*)UnitMult_2d[c1][c2][term]->ProjectionY(pronameUM->Data(),13,13);// 11 means 1000 pions
 401         norm_2_UM[term] = UnitMult[c1][c2][term]->Integral(UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.0),UnitMult[c1][c2][term]->GetXaxis()->FindBin(1.2));
 402         UnitMult[c1][c2][term]->Scale(norm_2_UM[0]/norm_2_UM[term]);
 403         //
 404         UnitMult[c1][c2][term]->SetMarkerStyle(20);
 405         UnitMult[c1][c2][term]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 406         UnitMult[c1][c2][term]->GetYaxis()->SetTitle("C_{2}");
 407         UnitMult[c1][c2][term]->SetTitle("");
 408       }// term
 409
 410
 411       // 3-particle
 412       for(int c3=0; c3<2; c3++){
 413         for(int term=0; term<5; term++){
 414
 415           TString *name3 = new TString("ThreeParticle_Charge1_");
 416           *name3 += c1;
 417           name3->Append("_Charge2_");
 418           *name3 += c2;
 419           name3->Append("_Charge3_");
 420           *name3 += c3;
 421           name3->Append("_M_");
 422           *name3 += Mbin;
 423           name3->Append("_ED_");
 424           *name3 += EDbin;
 425           name3->Append("_Term_");
 426           *name3 += term+1;
 427
 428           TString *nameNorm3=new TString(name3->Data());
 429           nameNorm3->Append("_Norm");
 430           //
 431           TString *nameK3=new TString(name3->Data());
 432           nameK3->Append("_Kfactor");
 433           //
 434           TString *nameTPN3=new TString(name3->Data());
 435           nameTPN3->Append("_TwoPartNorm");
 436           //
 437           TString *nameNegTPN3=new TString(name3->Data());
 438           nameNegTPN3->Append("_TwoPartNegNorm");
 439           //
 440           name3->Append("_1D");
 441
 442           ///////////////////////////////////////
 443           // skip degenerate histograms
 444           if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
 445           if( (c1+c2+c3)==2) {if(c1!=0) continue;}
 446           ////////////////////////////////////////
 447
 448           norm_3[term] = ((TH1D*)MyList->FindObject(nameNorm3->Data()))->Integral();
 449           ThreeParticle[c1][c2][c3][term] = (TH1D*)MyList->FindObject(name3->Data());
 450           ThreeParticle[c1][c2][c3][term]->Sumw2();
 451           //if(c1==c2 && c1==c3) cout<<"3-pion norms  "<<norm_3[term]<<endl;
 452           ThreeParticle[c1][c2][c3][term]->Scale(norm_3[0]/norm_3[term]);
 453           ThreeParticle[c1][c2][c3][term]->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
 454           ThreeParticle[c1][c2][c3][term]->GetYaxis()->SetTitle("C_{3}");
 455           ThreeParticle[c1][c2][c3][term]->SetMarkerStyle(20);
 456           ThreeParticle[c1][c2][c3][term]->SetTitle("");
 457           //
 458           ThreeParticle[c1][c2][c3][term]->Rebin(ThreeParticleRebin);
 459
 460           //
 461           if(term<4){
 462             K3avg[c1][c2][c3][term] = (TProfile*)MyList->FindObject(nameK3->Data());
 463             K3avg[c1][c2][c3][term]->Rebin(ThreeParticleRebin);
 464             if(MCcase || !FSICorrection) {
 465               K3avg[c1][c2][c3][term]->Reset();
 466               for(int ii=1; ii<=K3avg[c1][c2][c3][term]->GetNbinsX(); ii++) {
 467                 K3avg[c1][c2][c3][term]->Fill(K3avg[c1][c2][c3][term]->GetXaxis()->GetBinCenter(ii), 1);
 468               }
 469             }
 470           }
 471           //
 472           if(term==4 && c1==c2 && c1==c3){
 473             TPFullWeight_ThreeParticle_2D[c1] = (TH2D*)MyList->FindObject(nameTPN3->Data());
 474             TPFullWeight_ThreeParticle_2D[c1]->Scale(norm_3[0]/norm_3[term]);
 475             TPFullWeight_ThreeParticle_2D[c1]->RebinY(ThreeParticleRebin);
 476             //
 477             TPNegFullWeight_ThreeParticle_2D[c1] = (TH2D*)MyList->FindObject(nameNegTPN3->Data());
 478             TPNegFullWeight_ThreeParticle_2D[c1]->Scale(norm_3[0]/norm_3[term]);
 479             TPNegFullWeight_ThreeParticle_2D[c1]->RebinY(ThreeParticleRebin);
 480             //
 481             for(int binx=1; binx<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsX(); binx++) {
 482               for(int biny=1; biny<=TPFullWeight_ThreeParticle_2D[c1]->GetNbinsY(); biny++) {
 483                 TPFullWeight_ThreeParticle_2D[c1]->SetBinError(binx, biny, 0);
 484                 if(binx!=4){
 485                   TPFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny));
 486                   TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(4, biny));
 487                   if(InterpCorrection){
 488                     double q3 = TPFullWeight_ThreeParticle_2D[c1]->GetYaxis()->GetBinCenter(biny);
 489                     double InterCorr = 1.024 - 0.2765*q3;
 490                     //if(q3<0.1) cout<<binx<<"  "<<biny<<"  "<<q3<<" "<<InterCorr<<endl;
 491                     TPFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny));
 492                     TPNegFullWeight_ThreeParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_ThreeParticle_2D[c1]->GetBinContent(binx, biny));
 493                   }
 494                 }
 495               }
 496             }
 497
 498             TString *proName=new TString(nameTPN3->Data()); TString *proNameNeg=new TString(nameNegTPN3->Data());
 499             proName->Append("_pro");  proNameNeg->Append("_pro");
 500             TPN_ThreeParticle[c1] = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin);
 501             //
 502             proName->Append("_FullWeight"); proNameNeg->Append("_FullWeight");
 503             TPFullWeight_ThreeParticle[c1] = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), Gbin, Gbin);
 504             TPNegFullWeight_ThreeParticle[c1] = (TH1D*)TPNegFullWeight_ThreeParticle_2D[c1]->ProjectionY(proNameNeg->Data(), Gbin, Gbin);
 505             proName->Append("_FullWeightDen"); proNameNeg->Append("_FullWeightDen");
 506             TH1D *tempDen = (TH1D*)TPFullWeight_ThreeParticle_2D[c1]->ProjectionY(proName->Data(), 4, 4);
 507             TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[c1]->ProjectionY(proNameNeg->Data(), 4, 4);
 508             // Add Pos with Neg weights
 509             tempDen->Add(tempDenNeg);
 510             TPFullWeight_ThreeParticle[c1]->Add(TPNegFullWeight_ThreeParticle[c1]);
 511             //
 512             //TPFullWeight_ThreeParticle[c1]->Add(tempDen);// now added above in Interp section
 513             TPFullWeight_ThreeParticle[c1]->Divide(tempDen);
 514             TPFullWeight_ThreeParticle[c1]->SetLineColor(1);
 515             //
 516           }
 517
 518         }// term 3-particle
 519
 520
 521         // 4-particle
 522         for(int c4=0; c4<2; c4++){
 523           for(int term=0; term<13; term++){
 524
 525             TString *name4 = new TString("FourParticle_Charge1_");
 526             *name4 += c1;
 527             name4->Append("_Charge2_");
 528             *name4 += c2;
 529             name4->Append("_Charge3_");
 530             *name4 += c3;
 531             name4->Append("_Charge4_");
 532             *name4 += c4;
 533             name4->Append("_M_");
 534             *name4 += Mbin;
 535             name4->Append("_ED_");
 536             *name4 += EDbin;
 537             name4->Append("_Term_");
 538             *name4 += term+1;
 539
 540             TString *nameNorm4=new TString(name4->Data());
 541             nameNorm4->Append("_Norm");
 542             //
 543             TString *nameK4=new TString(name4->Data());
 544             nameK4->Append("_Kfactor");// or Weighted
 545             //
 546             TString *nameTPN4=new TString(name4->Data());
 547             nameTPN4->Append("_TwoPartNorm");
 548             //
 549             TString *nameNegTPN4=new TString(name4->Data());
 550             nameNegTPN4->Append("_TwoPartNegNorm");
 551             //
 552             TString *nameFitBuild=new TString(name4->Data());
 553             nameFitBuild->Append("_FullBuildFromFits");
 554             //
 555             TString *namePartialFitBuild=new TString(name4->Data());
 556             namePartialFitBuild->Append("_PartialBuildFromFits");
 557             //
 558             name4->Append("_1D");
 559             ///////////////////////////////////////
 560             // skip degenerate histograms
 561             if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
 562             if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
 563             if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
 564             /////////////////////////////////////////
 565             norm_4[term] = ((TH1D*)MyList->FindObject(nameNorm4->Data()))->Integral();
 566
 567             //if( (c1+c2+c3+c4)==4) cout<<"4-pion norms  "<<norm_4[term]<<endl;
 568             if(norm_4[term] > 0){
 569               //if(c1==c2 && c1==c3 && c1==c4) cout<<term<<"  "<<norm_4[0]/norm_4[term]<<endl;
 570               FourParticle[c1][c2][c3][c4][term] = (TH1D*)MyList->FindObject(name4->Data());
 571               FourParticle[c1][c2][c3][c4][term]->Sumw2();
 572               FourParticle[c1][c2][c3][c4][term]->Scale(norm_4[0]/norm_4[term]);
 573               FourParticle[c1][c2][c3][c4][term]->GetXaxis()->SetTitle("Q_{4} (GeV/c)");
 574               FourParticle[c1][c2][c3][c4][term]->GetYaxis()->SetTitle("C_{4}");
 575               FourParticle[c1][c2][c3][c4][term]->SetMarkerStyle(20);
 576               FourParticle[c1][c2][c3][c4][term]->SetTitle("");
 577               //
 578               FourParticle[c1][c2][c3][c4][term]->Rebin(FourParticleRebin);
 579
 580             }else{
 581               cout<<"4-particle normalization = 0."<<endl;
 582             }
 583             if(term<12){
 584               K4avg[c1][c2][c3][c4][term] = (TProfile*)MyList->FindObject(nameK4->Data());
 585               K4avg[c1][c2][c3][c4][term]->Rebin(FourParticleRebin);
 586               if(MCcase || !FSICorrection) {
 587                 K4avg[c1][c2][c3][c4][term]->Reset();
 588                 for(int ii=1; ii<=K4avg[c1][c2][c3][c4][term]->GetNbinsX(); ii++) {
 589                   K4avg[c1][c2][c3][c4][term]->Fill(K4avg[c1][c2][c3][c4][term]->GetXaxis()->GetBinCenter(ii), 1);
 590                 }
 591               }
 592             }
 593             if(term==12 && c1==c2 && c1==c3 && c1==c4){
 594
 595               TPFullWeight_FourParticle_2D[c1] = (TH2D*)MyList->FindObject(nameTPN4->Data());
 596               TPFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]);
 597               TPFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin);
 598               //
 599               TPNegFullWeight_FourParticle_2D[c1] = (TH2D*)MyList->FindObject(nameNegTPN4->Data());
 600               TPNegFullWeight_FourParticle_2D[c1]->Scale(norm_4[0]/norm_4[term]);
 601               TPNegFullWeight_FourParticle_2D[c1]->RebinY(FourParticleRebin);
 602               //
 603               if(c1==0 && !MCcase){
 604                 FullBuildFromFits_3D = (TH3D*)MyList->FindObject(nameFitBuild->Data());
 605                 FullBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]);
 606                 FullBuildFromFits_3D->RebinZ(FourParticleRebin);
 607                 //
 608                 PartialBuildFromFits_3D = (TH3D*)MyList->FindObject(namePartialFitBuild->Data());
 609                 PartialBuildFromFits_3D->Scale(norm_4[0]/norm_4[term]);
 610                 PartialBuildFromFits_3D->RebinZ(FourParticleRebin);
 611               }
 612               //
 613               for(int binx=1; binx<=TPFullWeight_FourParticle_2D[c1]->GetNbinsX(); binx++) {
 614                 for(int biny=1; biny<=TPFullWeight_FourParticle_2D[c1]->GetNbinsY(); biny++) {
 615                   TPFullWeight_FourParticle_2D[c1]->SetBinError(binx, biny, 0);
 616                   if(binx!=4){
 617                     TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny));
 618                     TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny) + TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(4, biny));
 619
 620                     if(InterpCorrection){// Interpolator correction
 621                       double q4 = TPFullWeight_FourParticle_2D[c1]->GetYaxis()->GetBinCenter(biny);
 622                       double InterCorr = 1.032 - 0.2452*q4;
 623                       TPFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny));
 624                       TPNegFullWeight_FourParticle_2D[c1]->SetBinContent(binx, biny, InterCorr * TPNegFullWeight_FourParticle_2D[c1]->GetBinContent(binx, biny));
 625                     }
 626                   }
 627                 }
 628               }
 629               TString *proName=new TString(nameTPN4->Data()); TString *proNegName=new TString(nameNegTPN4->Data());
 630               TString *proNameFitBuild=new TString(nameFitBuild->Data()); TString *proNamePartialFitBuild=new TString(namePartialFitBuild->Data());
 631               proName->Append("_pro"); proNegName->Append("_pro"); proNameFitBuild->Append("_pro"); proNamePartialFitBuild->Append("_pro");
 632               TPN_FourParticle[c1] = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), TPNbin, TPNbin);
 633               //
 634               proName->Append("_FullWeight"); proNegName->Append("_FullWeight");
 635               TPFullWeight_FourParticle[c1] = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), Gbin, Gbin);
 636               TPNegFullWeight_FourParticle[c1] = (TH1D*)TPNegFullWeight_FourParticle_2D[c1]->ProjectionY(proNegName->Data(), Gbin, Gbin);
 637               proName->Append("_FullWeightDen"); proNegName->Append("_FullWeightDen");
 638               TH1D *tempDen = (TH1D*)TPFullWeight_FourParticle_2D[c1]->ProjectionY(proName->Data(), 4, 4);
 639               TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[c1]->ProjectionY(proNegName->Data(), 4, 4);
 640               //
 641               // Add Pos and Neg Weights
 642               tempDen->Add(tempDenNeg);
 643               TPFullWeight_FourParticle[c1]->Add(TPNegFullWeight_FourParticle[c1]);
 644               //
 645               //TPFullWeight_FourParticle[c1]->Add(tempDen);// now added above in Interp section
 646               TPFullWeight_FourParticle[c1]->Divide(tempDen);
 647               TPFullWeight_FourParticle[c1]->SetLineColor(1);
 648               /*TString *ErrName=new TString(nameTPN4->Data());
 649               ErrName->Append("Err");
 650               TH2D *temperr2D = (TH2D*)MyList->FindObject(ErrName->Data());
 651               TH1D *temperr = (TH1D*)temperr2D->ProjectionY("tesst",4,4);
 652               temperr->Rebin(FourParticleRebin);
 653               cout.precision(8);
 654               cout<<temperr->GetBinContent(3)<<endl;
 655               cout<<(temperr->GetBinContent(5) / tempDen->GetBinContent(5))<<"  "<<TPFullWeight_FourParticle[c1]->GetBinContent(5)<<endl;
 656               */
 657               if(c1==0 && !MCcase){
 658                 FullBuildFromFits = (TH1D*)FullBuildFromFits_3D->ProjectionZ(proNameFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin);
 659                 TH1D *tempDen2 = (TH1D*)FullBuildFromFits_3D->ProjectionZ("tempDen2", c3FitType, c3FitType, 4, 4);
 660                 tempDen2->Scale(1/100.);// It was filled 100 times with the same value
 661                 FullBuildFromFits->Add(tempDen2);
 662                 FullBuildFromFits->Divide(tempDen2);
 663                 FullBuildFromFits->SetLineColor(4);
 664                 //
 665                 PartialBuildFromFits = (TH1D*)PartialBuildFromFits_3D->ProjectionZ(proNamePartialFitBuild->Data(), c3FitType, c3FitType, Gbin, Gbin);
 666                 PartialBuildFromFits->Add(tempDen2);
 667                 PartialBuildFromFits->Divide(tempDen2);
 668                 PartialBuildFromFits->SetLineColor(2);
 669               }
 670             }
 671           }// term 4-particle
 672
 673         }// c4
 674       }// c3
 675     }// c2
 676   }// c1
 677
 678   TH1D *fTPNRejects3pion = (TH1D*)MyList->FindObject("fTPNRejects3pion2");
 679   TH1D *fTPNRejects4pion1 = (TH1D*)MyList->FindObject("fTPNRejects4pion1");
 680
 681   //TH2D *c4QSFitNum2D = (TH2D*)MyList->FindObject("fc4QSFitNum");
 682   //TH2D *c4QSFitDen2D = (TH2D*)MyList->FindObject("fc4QSFitDen");
 683   //c4QSFitNum2D->RebinY(3);
 684   //c4QSFitDen2D->RebinY(3);
 685
 686   cout<<"Done getting Histograms"<<endl;
 687
 688   TF1 *Unity=new TF1("Unity","1",0,100);
 689   Unity->SetLineStyle(2);
 690
 691
 692   int ch1_2=0,ch2_2=0;
 693   int ch1_3=0,ch2_3=0,ch3_3=0;
 694   int ch1_4=0,ch2_4=0,ch3_4=0,ch4_4=0;
 695
 696   if(SameCharge && CHARGE==+1) {ch1_2=1; ch2_2=1;   ch1_3=1; ch2_3=1; ch3_3=1;   ch1_4=1; ch2_4=1; ch3_4=1; ch4_4=1;}
 697   if(!SameCharge){
 698     ch1_2=0; ch2_2=1;
 699     //
 700     if(CHARGE==-1) { ch1_3=0; ch2_3=0; ch3_3=1; }
 701     else { ch1_3=0; ch2_3=1; ch3_3=1; }
 702     //
 703     if(CHARGE==-1 && MixedCharge4pionType==1) {ch1_4=0; ch2_4=0; ch3_4=0; ch4_4=1;}
 704     if(CHARGE==+1 && MixedCharge4pionType==1) {ch1_4=0; ch2_4=1; ch3_4=1; ch4_4=1;}
 705     if(MixedCharge4pionType==2) {ch1_4=0; ch2_4=0; ch3_4=1; ch4_4=1;}
 706   }
 707
 708   ///////////////////////////////////////////////////////////
 709   // 2-pion section
 710   cout<<"2-pion section"<<endl;
 711
 712   TCanvas *can1 = new TCanvas("can1", "can1",11,53,700,500);
 713   can1->SetHighLightColor(2);
 714   can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 715   gStyle->SetOptFit(0111);
 716   can1->SetFillColor(10);
 717   can1->SetBorderMode(0);
 718   can1->SetBorderSize(2);
 719   can1->SetGridx();
 720   can1->SetGridy();
 721   can1->SetFrameFillColor(0);
 722   can1->SetFrameBorderMode(0);
 723   can1->SetFrameBorderMode(0);
 724   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
 725
 726   TLegend *legend1 = new TLegend(.6,.67,.87,.87,NULL,"brNDC");
 727   legend1->SetBorderSize(1);
 728   legend1->SetTextSize(.04);// small .03; large .036
 729   legend1->SetFillColor(0);
 730
 731
 732   TH1D *TERM1_2=(TH1D*)TwoParticle[ch1_2][ch2_2][0]->Clone();
 733   TH1D *TERM2_2=(TH1D*)TwoParticle[ch1_2][ch2_2][1]->Clone();
 734   //
 735   if(SameCharge){
 736     TERM1_2->Multiply(MRC_SC_2[0]);
 737     TERM2_2->Multiply(MRC_SC_2[1]);
 738   }else {
 739     TERM1_2->Multiply(MRC_MC_2[0]);
 740     TERM2_2->Multiply(MRC_MC_2[1]);
 741   }
 742
 743   TH1D *C2=(TH1D*)TERM1_2->Clone();
 744   C2->Divide(TERM2_2);
 745   C2->GetXaxis()->SetRangeUser(0,0.1);
 746   C2->SetMinimum(0.98); C2->SetMaximum(2.5);
 747   /*C2->GetXaxis()->SetTitleSize(0.06);
 748   C2->GetXaxis()->SetLabelSize(0.06);
 749   C2->GetYaxis()->SetTitleSize(0.06);
 750   C2->GetYaxis()->SetLabelSize(0.06);
 751   C2->GetXaxis()->SetTitleOffset(1.05);
 752   C2->GetYaxis()->SetTitleOffset(1.1);
 753   C2->GetXaxis()->SetNdivisions(606);
 754   C2->GetYaxis()->SetNdivisions(505);*/
 755   C2->Draw();
 756
 757
 758   TH1D *C2QS=(TH1D*)TERM1_2->Clone();
 759   C2QS->Add(TERM2_2, -(1-TwoFrac));
 760   C2QS->Scale(1/TwoFrac);
 761   for(int bin=1; bin<=C2QS->GetNbinsX(); bin++) {
 762     Float_t K2 = 1.0;
 763     if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, C2QS->GetXaxis()->GetBinCenter(bin));
 764     C2QS->SetBinContent(bin, C2QS->GetBinContent(bin)/K2);
 765     C2QS->SetBinError(bin, C2QS->GetBinError(bin)/K2);
 766   }
 767   C2QS->Divide(TERM2_2);
 768   if(SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionSC);
 769   if(!SameCharge && MuonCorrection) C2QS->Multiply(C2muonCorrectionMC);
 770   C2QS->SetMarkerColor(2); C2QS->SetLineColor(2);
 771   if(!MCcase) C2QS->Draw("same");
 772
 773
 774   // To visualize the Qcut and Norm Regions
 775   //TH1D *QcutRegion = new TH1D("QcutRegion","",400,0,2);
 776   //TH1D *NormRegion1 = new TH1D("NormRegion1","",400,0,2);
 777   //TH1D *NormRegion2 = new TH1D("NormRegion2","",400,0,2);
 778   //TERM1_2->Scale(1/TERM1_2->Integral());
 779   //for(int bin=1; bin<=16; bin++) QcutRegion->SetBinContent(bin,TERM1_2->GetBinContent(bin));
 780   //for(int bin=213; bin<=220; bin++) NormRegion1->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
 781   //for(int bin=31; bin<=35; bin++) NormRegion2->SetBinContent(bin,Two_ex[0][0][0][0]->GetBinContent(bin));
 782   //QcutRegion->SetFillColor(4);
 783   //NormRegion1->SetFillColor(2);
 784   //NormRegion2->SetFillColor(3);
 785   //TERM1_2->GetYaxis()->SetTitle("Pair probability");
 786   //TERM1_2->GetYaxis()->SetTitleOffset(1.3);
 787   //TERM1_2->GetXaxis()->SetRangeUser(0,1.6);
 788   //TERM1_2->Draw();
 789   //QcutRegion->Draw("same");
 790   //cout<<TERM1_2->Integral(1,16) / TERM1_2->Integral()<<endl;
 791
 792   // Print 2-pion values
 793   /*for(int bin=1; bin<=20; bin++){
 794     cout<<C2QS->GetBinContent(bin)<<", ";
 795   }
 796   cout<<endl;
 797   for(int bin=1; bin<=20; bin++){
 798     cout<<C2QS->GetBinError(bin)<<", ";
 799   }
 800   cout<<endl;*/
 801   double C2refPoints_0p65[20]={-0.218603, 1.00087, 1.00215, 1.00209, 1.00202, 1.00215, 1.00162, 1.00129, 1.00118, 1.00076, 1.00049, 1.00021, 0.999996, 0.999733, 0.999622, 0.999554, 0.999764, 0.999679, 0.999677, 0.999641};
 802   double C2refPoints_e_0p65[20]={0.309152, 0.000473757, 0.000305599, 0.00022603, 0.000179596, 0.000149295, 0.000128009, 0.000112371, 0.00010047, 9.11363e-05, 8.37061e-05, 7.76871e-05, 7.27445e-05, 6.86217e-05, 6.5155e-05, 6.22148e-05, 5.9718e-05, 5.75613e-05, 5.57066e-05, 5.41028e-05};
 803   double C2refPoints_0p75[20]={-0.135326, 0.968257, 0.984681, 0.991465, 0.99503, 0.997249, 0.998093, 0.998645, 0.999135, 0.999172, 0.999237, 0.999213, 0.999197, 0.999095, 0.999105, 0.999123, 0.999385, 0.999366, 0.999408, 0.999411};
 804   double C2refPoints_e_0p75[20]={0.19138, 0.000422359, 0.000272193, 0.000201183, 0.000159788, 0.000132797, 0.000113845, 9.99271e-05, 8.93377e-05, 8.10334e-05, 7.44239e-05, 6.90702e-05, 6.46743e-05, 6.10077e-05, 5.79248e-05, 5.53103e-05, 5.30903e-05, 5.11725e-05, 4.95234e-05, 4.80974e-05};
 805   TH1D *C2ref_0p65=(TH1D*)C2QS->Clone();
 806   TH1D *C2ref_0p75=(TH1D*)C2QS->Clone();
 807   for(int bin=1; bin<=20; bin++){
 808     C2ref_0p65->SetBinContent(bin, C2refPoints_0p65[bin-1]);
 809     C2ref_0p65->SetBinError(bin, C2refPoints_e_0p65[bin-1]);
 810     C2ref_0p75->SetBinContent(bin, C2refPoints_0p75[bin-1]);
 811     C2ref_0p75->SetBinError(bin, C2refPoints_e_0p75[bin-1]);
 812   }
 813   C2ref_0p65->SetMarkerColor(4); C2ref_0p65->SetLineColor(4);
 814   C2ref_0p75->SetMarkerColor(6); C2ref_0p75->SetLineColor(6);
 815   //C2ref_0p65->Draw("same");
 816   //C2ref_0p75->Draw("same");
 817
 818
 819   TH1D *UnitMultC2=(TH1D*)UnitMult[0][0][0]->Clone();
 820   TH1D *UnitMultC2den=(TH1D*)UnitMult[0][0][1]->Clone();
 821   UnitMultC2->Divide(UnitMultC2den);
 822   UnitMultC2->GetXaxis()->SetRangeUser(0,0.3);
 823   UnitMultC2->SetMinimum(0.97);
 824   UnitMultC2->SetMaximum(1.25);
 825   UnitMultC2->Draw();
 826
 827   TH1D *UnitMultC2QS=(TH1D*)UnitMult[0][0][0]->Clone();
 828   UnitMultC2QS->Add(UnitMult[0][0][1], -(1-TwoFrac));
 829   UnitMultC2QS->Scale(1/TwoFrac);
 830   for(int bin=1; bin<=UnitMultC2QS->GetNbinsX(); bin++) {
 831     Float_t K2 = 1.0;
 832     if(FSICorrection) K2 = FSICorrelation(ch1_2, ch2_2, UnitMultC2QS->GetXaxis()->GetBinCenter(bin));
 833     UnitMultC2QS->SetBinContent(bin, UnitMultC2QS->GetBinContent(bin)/K2);
 834     UnitMultC2QS->SetBinError(bin, UnitMultC2QS->GetBinError(bin)/K2);
 835   }
 836   UnitMultC2QS->Divide(UnitMultC2den);
 837   UnitMultC2QS->SetMarkerColor(2); UnitMultC2QS->SetLineColor(2);
 838   UnitMultC2QS->Draw("same");
 839
 840   ///////////////////////////////////////////////////////////
 841   // 3-pion
 842   cout<<"3-pion section"<<endl;
 843
 844   TCanvas *can2 = new TCanvas("can2", "can2",600,53,700,500);
 845   can2->SetHighLightColor(2);
 846   can2->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 847   gStyle->SetOptFit(0111);
 848   can2->SetFillColor(10);
 849   can2->SetBorderMode(0);
 850   can2->SetBorderSize(2);
 851   can2->SetGridx();
 852   can2->SetGridy();
 853   can2->SetFrameFillColor(0);
 854   can2->SetFrameBorderMode(0);
 855   can2->SetFrameBorderMode(0);
 856   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
 857   TLegend *legend2 = (TLegend*)legend1->Clone();
 858   legend2->SetX1(0.45); legend2->SetX2(0.95);  legend2->SetY1(0.6); legend2->SetY2(0.95);
 859
 860   TH1D *TERM1_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][0]->Clone();
 861   TH1D *TERM2_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][1]->Clone();
 862   TH1D *TERM3_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][2]->Clone();
 863   TH1D *TERM4_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][3]->Clone();
 864   TH1D *TERM5_3=(TH1D*)ThreeParticle[ch1_3][ch2_3][ch3_3][4]->Clone();
 865   //
 866
 867   if(SameCharge){
 868     TERM1_3->Multiply(MRC_SC_3[0]);
 869     TERM2_3->Multiply(MRC_SC_3[1]);
 870     TERM3_3->Multiply(MRC_SC_3[1]);
 871     TERM4_3->Multiply(MRC_SC_3[1]);
 872     TERM5_3->Multiply(MRC_SC_3[2]);
 873   }else{
 874     if(CHARGE==-1){// --+ but MRC is stored as -++
 875       TERM1_3->Multiply(MRC_MC_3[0]);
 876       TERM2_3->Multiply(MRC_MC_3[2]);
 877       TERM3_3->Multiply(MRC_MC_3[1]);
 878       TERM4_3->Multiply(MRC_MC_3[1]);
 879       TERM5_3->Multiply(MRC_MC_3[3]);
 880     }else{
 881       TERM1_3->Multiply(MRC_MC_3[0]);
 882       TERM2_3->Multiply(MRC_MC_3[1]);
 883       TERM3_3->Multiply(MRC_MC_3[1]);
 884       TERM4_3->Multiply(MRC_MC_3[2]);
 885       TERM5_3->Multiply(MRC_MC_3[3]);
 886     }
 887   }
 888
 889
 890   TH1D *N3QS = (TH1D*)TERM1_3->Clone();
 891   N3QS->Add(TERM5_3, -f31);
 892   N3QS->Add(TERM2_3, -f32);
 893   N3QS->Add(TERM3_3, -f32);
 894   N3QS->Add(TERM4_3, -f32);
 895   N3QS->Scale(1/f33);
 896   N3QS->Multiply(K3avg[ch1_3][ch2_3][ch3_3][0]);
 897
 898
 899   TH1D *C3QS = (TH1D*)N3QS->Clone();
 900   C3QS->Divide(TERM5_3);
 901   if(SameCharge) C3QS->Multiply(C3muonCorrectionSC[0]);
 902   else C3QS->Multiply(C3muonCorrectionMC[0]);
 903
 904
 905   C3QS->SetMarkerStyle(20);
 906   C3QS->SetMarkerColor(4);
 907   C3QS->GetXaxis()->SetRangeUser(0,0.1);
 908   C3QS->SetMinimum(0.9); C3QS->SetMaximum(5.0);
 909
 910   TH1D *n3QS = (TH1D*)N3QS->Clone();
 911   TH1D *c3QS = (TH1D*)N3QS->Clone();
 912   for(int bin=1; bin<=n3QS->GetNbinsX(); bin++){
 913     if(n3QS->GetBinContent(bin)==0) continue;
 914     double MuonCorr1=1.0, MuonCorr2=1.0, MuonCorr3=1.0, MuonCorr4=1.0;
 915     if(SameCharge){
 916       MuonCorr1 = C3muonCorrectionSC[0]->GetBinContent(bin);
 917       MuonCorr2 = C3muonCorrectionSC[1]->GetBinContent(bin);
 918       MuonCorr3 = MuonCorr2;
 919       MuonCorr4 = MuonCorr2;
 920     }else if(CHARGE==-1){
 921       MuonCorr1 = C3muonCorrectionMC[0]->GetBinContent(bin);
 922       MuonCorr2 = C3muonCorrectionMC[2]->GetBinContent(bin);// --
 923       MuonCorr3 = C3muonCorrectionMC[1]->GetBinContent(bin);// -+
 924       MuonCorr4 = MuonCorr3;// -+
 925     }else{
 926       MuonCorr1 = C3muonCorrectionMC[0]->GetBinContent(bin);
 927       MuonCorr2 = C3muonCorrectionMC[1]->GetBinContent(bin);// -+
 928       MuonCorr3 = MuonCorr2;
 929       MuonCorr4 = C3muonCorrectionMC[2]->GetBinContent(bin);
 930     }
 931     double value = n3QS->GetBinContent(bin) * MuonCorr1;
 932     value -= ( TERM2_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][1]->GetBinContent(bin) / TwoFrac * MuonCorr2;
 933     value -= ( TERM3_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][2]->GetBinContent(bin) / TwoFrac * MuonCorr3;
 934     value -= ( TERM4_3->GetBinContent(bin) - (1-TwoFrac)*TERM5_3->GetBinContent(bin)) * K3avg[ch1_3][ch2_3][ch3_3][3]->GetBinContent(bin) / TwoFrac * MuonCorr4;
 935     //value -= TERM5_3->GetBinContent(bin)*(MuonCorr1 - MuonCorr2 - MuonCorr3 - MuonCorr4);
 936     value += 2*TERM5_3->GetBinContent(bin);
 937     //
 938     n3QS->SetBinContent(bin, value);
 939     c3QS->SetBinContent(bin, value + TERM5_3->GetBinContent(bin));
 940   }
 941   c3QS->Divide(TERM5_3);
 942
 943   C3QS->Draw();
 944   c3QS->Draw("same");
 945   //int lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
 946   //int highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
 947   //double SF=C3QS->Integral(lowBin, highBin); SF /= TPFullWeight_ThreeParticle[ch1_3]->Integral(lowBin, highBin);
 948   //TPFullWeight_ThreeParticle[ch1_3]->Scale(SF);
 949   //
 950   if(SameCharge) TPFullWeight_ThreeParticle[ch1_3]->Draw("same");
 951
 952   //cout<<TPFullWeight_ThreeParticle[ch1_3]->GetBinContent(2)<<endl;
 953   //TH1D *C3raw = (TH1D*)TERM1_3->Clone();
 954   //C3raw->Divide(TERM5_3);
 955   //C3raw->SetMarkerColor(4);
 956   //C3raw->Draw("same");
 957
 958   legend2->AddEntry(C3QS,"C_{3}^{QS}","p");
 959   legend2->AddEntry(c3QS,"c_{3}^{QS}{2-pion removal}","p");
 960   if(SameCharge) legend2->AddEntry(TPFullWeight_ThreeParticle[ch1_3],"C_{3}^{QS} built","l");
 961   legend2->Draw("same");
 962
 963
 964   ///////////////////////////////////////
 965   // C3QS built ratio
 966   TCanvas *can2_2 = new TCanvas("can2_2", "can2_2",1200,53,700,500);
 967   can2_2->SetHighLightColor(2);
 968   can2_2->Range(-0.7838115,-1.033258,0.7838115,1.033258);
 969   gStyle->SetOptFit(0111);
 970   can2_2->SetFillColor(10);
 971   can2_2->SetBorderMode(0);
 972   can2_2->SetBorderSize(2);
 973   can2_2->SetGridx();
 974   can2_2->SetGridy();
 975   can2_2->SetFrameFillColor(0);
 976   can2_2->SetFrameBorderMode(0);
 977   can2_2->SetFrameBorderMode(0);
 978   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
 979   TLegend *legend2_2 = (TLegend*)legend1->Clone();
 980   legend2_2->SetX1(0.25); legend2_2->SetX2(0.6);  legend2_2->SetY1(0.8); legend2_2->SetY2(0.95);
 981
 982   TH1D *C3QSratio = (TH1D*)C3QS->Clone();
 983
 984   //TH1D *Chi2NDF_PointSize_3 = new TH1D("Chi2NDF_PointSize_3","",40,-0.5,39.5);
 985   //TH1D *Chi2NDF_FullSize_3 = new TH1D("Chi2NDF_FullSize_3","",40,-0.5,39.5);
 986   TH1D *Chi2NDF_PointSize_3 = new TH1D("Chi2NDF_PointSize_3","",100,-0.5,99.5);
 987   TH1D *Chi2NDF_FullSize_3 = new TH1D("Chi2NDF_FullSize_3","",100,-0.5,99.5);
 988   Chi2NDF_PointSize_3->SetLineColor(4); Chi2NDF_FullSize_3->SetLineColor(2);
 989   Chi2NDF_PointSize_3->SetMarkerColor(4); Chi2NDF_FullSize_3->SetMarkerColor(2);
 990   Chi2NDF_PointSize_3->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_3->GetYaxis()->SetTitle("#sqrt{#chi^{2}}");
 991
 992   if(SameCharge && CollisionType==0){
 993
 994     TH1D *tempDen = (TH1D*)TPFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPFullWeight3_Den", 4, 4);
 995     TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY("TPNegFullWeight3_Den", 4, 4);
 996     tempDen->Add(tempDenNeg);// Add Pos and Neg Den
 997
 998     for(int binG=5; binG<=104; binG++){// was 44
 999       TString *proName=new TString("TPFullWeight3_");
1000       *proName += binG;
1001       TH1D *tempNum = (TH1D*)TPFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY(proName->Data(), binG, binG);
1002       proName->Append("_Neg");
1003       TH1D *tempNumNeg = (TH1D*)TPNegFullWeight_ThreeParticle_2D[ch1_3]->ProjectionY(proName->Data(), binG, binG);
1004       //
1005       // Add Pos and Neg Num
1006       tempNum->Add(tempNumNeg);
1007       //
1008       //tempNum->Add(tempDen);// now added in histogram read section
1009       tempNum->Divide(tempDen);
1010       //lowBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
1011       //highBin = C3QS->GetXaxis()->FindBin(Cutoff_FullWeight_Q3[Mbin]);
1012       //SF=C3QS->Integral(lowBin, highBin);
1013       //SF /= tempNum->Integral(lowBin, highBin);
1014       //tempNum->Scale(SF);
1015
1016       double Chi2=0;
1017       double NDF=0;
1018       for(int binQ3=3; binQ3<=3; binQ3++){
1019         if(tempNum->GetBinContent(binQ3) <=0) continue;
1020         double value = C3QS->GetBinContent(binQ3) - tempNum->GetBinContent(binQ3);
1021         double err = C3QS->GetBinError(binQ3);
1022         if(err <=0) continue;
1023
1024         Chi2 += pow(value / err,2);
1025         //
1026         NDF += 1;
1027       }
1028       //if(binG<25) {Chi2NDF_PointSize_3->SetBinContent(1 + 2*(binG-5), sqrt(Chi2)/NDF); Chi2NDF_PointSize_3->SetBinError(1 + 2*(binG-5), 0.001);}
1029       //else {Chi2NDF_FullSize_3->SetBinContent(1 + 2*(binG-25), sqrt(Chi2)/NDF); Chi2NDF_FullSize_3->SetBinError(1 + 2*(binG-25), 0.001);}
1030       if(binG<55) {Chi2NDF_PointSize_3->SetBinContent(1 + 2*(binG-5), sqrt(Chi2)/NDF); Chi2NDF_PointSize_3->SetBinError(1 + 2*(binG-5), 0.001);}
1031       else {Chi2NDF_FullSize_3->SetBinContent(1 + 2*(binG-55), sqrt(Chi2)/NDF); Chi2NDF_FullSize_3->SetBinError(1 + 2*(binG-55), 0.001);}
1032       //if(NDF>0) cout<<binG<<"  "<<sqrt(Chi2)/NDF<<endl;
1033     }
1034     Chi2NDF_PointSize_3->SetMarkerStyle(20); Chi2NDF_FullSize_3->SetMarkerStyle(20);
1035
1036     C3QSratio->Divide(TPFullWeight_ThreeParticle[ch1_3]);
1037     //
1038     C3QSratio->SetMinimum(0.94); C3QSratio->SetMaximum(1.02);
1039     C3QSratio->GetYaxis()->SetTitleOffset(1.2);
1040     C3QSratio->GetYaxis()->SetTitle("C_{3}^{QS} / C_{3}(built)");
1041     C3QSratio->Draw();
1042     Unity->Draw("same");
1043
1044     Chi2NDF_PointSize_3->Draw();
1045     Chi2NDF_FullSize_3->Draw("same");
1046     legend2_2->AddEntry(Chi2NDF_PointSize_3,"R_{coh}=0 (Point Source)","p");
1047     legend2_2->AddEntry(Chi2NDF_FullSize_3,"R_{coh}=R_{ch}","p");
1048     legend2_2->Draw("same");
1049   }
1050
1051   // r3
1052   TH1D *r3;
1053   if(SameCharge && CollisionType==0){
1054     r3 = (TH1D*)n3QS->Clone();
1055     TPN_ThreeParticle[ch1_3]->Multiply(MRC_SC_3[2]);
1056     r3->Divide(TPN_ThreeParticle[ch1_3]);
1057     r3->GetXaxis()->SetRangeUser(0,0.08);
1058     r3->SetMinimum(0.5); r3->SetMaximum(2.5);
1059     r3->GetYaxis()->SetTitle("r_{3}");
1060     //
1061     r3->Draw();
1062     //TPN_ThreeParticle[ch1_3]->Draw();
1063     //fTPNRejects3pion->SetLineColor(2);
1064     //fTPNRejects3pion->Draw("same");
1065   }
1066
1067
1068   // Print 3-pion points
1069   for(int bin=1; bin<=10; bin++){
1070     //cout<<C3QS->GetBinContent(bin)<<", ";
1071     //cout<<c3QS->GetBinContent(bin)<<", ";
1072     //cout<<TPFullWeight_ThreeParticle[ch1_3]->GetBinContent(bin)<<", ";
1073   }
1074   //cout<<endl;
1075   for(int bin=1; bin<=10; bin++){
1076     //cout<<C3QS->GetBinError(bin)<<", ";
1077     //cout<<c3QS->GetBinError(bin)<<", ";
1078   }
1079   //cout<<endl;
1080
1081
1082   ////////////////////////////////////////////////////////////////
1083   // 4-pion
1084   cout<<"4-pion section"<<endl;
1085
1086   TCanvas *can3 = new TCanvas("can3", "can3",11,600,700,500);// 60 was 600
1087   can3->SetHighLightColor(2);
1088   can3->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1089   gStyle->SetOptFit(0111);
1090   can3->SetFillColor(10);
1091   can3->SetBorderMode(0);
1092   can3->SetBorderSize(2);
1093   can3->SetGridx();
1094   can3->SetGridy();
1095   can3->SetFrameFillColor(0);
1096   can3->SetFrameBorderMode(0);
1097   can3->SetFrameBorderMode(0);
1098   gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
1099   TLegend *legend3 = (TLegend*)legend1->Clone();
1100   legend3->SetX1(0.47); legend3->SetX2(0.98);  legend3->SetY1(0.7); legend3->SetY2(0.98);
1101
1102   TH1D *TERMS_4[13];
1103   for(int index=0; index<13; index++){
1104     TERMS_4[index] = (TH1D*)FourParticle[ch1_4][ch2_4][ch3_4][ch4_4][index]->Clone();
1105     //
1106     TH1D *MRC_temp;
1107     if(SameCharge){
1108       if(index==0) MRC_temp = (TH1D*)MRC_SC_4[0]->Clone();
1109       else if(index<=4) MRC_temp = (TH1D*)MRC_SC_4[1]->Clone();
1110       else if(index<=10) MRC_temp = (TH1D*)MRC_SC_4[2]->Clone();
1111       else if(index==11) MRC_temp = (TH1D*)MRC_SC_4[3]->Clone();
1112       else MRC_temp = (TH1D*)MRC_SC_4[4]->Clone();
1113     }else if(CHARGE==-1 && MixedCharge4pionType==1){
1114       if(index==0) MRC_temp = (TH1D*)MRC_MC1_4[0]->Clone();//0
1115       else if(index!=1 && index<=4) MRC_temp = (TH1D*)MRC_MC1_4[1]->Clone();//1
1116       else if(index==1) MRC_temp = (TH1D*)MRC_MC1_4[2]->Clone();//2
1117       else if(index==7 || index==9 || index==10) MRC_temp = (TH1D*)MRC_MC1_4[3]->Clone();//3
1118       else if(index!=12) MRC_temp = (TH1D*)MRC_MC1_4[4]->Clone();//4
1119       else MRC_temp = (TH1D*)MRC_MC1_4[5]->Clone();//5
1120     }else if(CHARGE==+1 && MixedCharge4pionType==1){
1121       if(index==0) MRC_temp = (TH1D*)MRC_MC1_4[0]->Clone();
1122       else if(index<=3) MRC_temp = (TH1D*)MRC_MC1_4[1]->Clone();
1123       else if(index==4) MRC_temp = (TH1D*)MRC_MC1_4[2]->Clone();
1124       else if(index==5 || index==6 || index==7) MRC_temp = (TH1D*)MRC_MC1_4[3]->Clone();
1125       else if(index!=12) MRC_temp = (TH1D*)MRC_MC1_4[4]->Clone();
1126       else MRC_temp = (TH1D*)MRC_MC1_4[5]->Clone();
1127     }else{// --++
1128       if(index==0) MRC_temp = (TH1D*)MRC_MC2_4[0]->Clone();
1129       else if(index<=4) MRC_temp = (TH1D*)MRC_MC2_4[1]->Clone();
1130       else if(index==5 || index==10) MRC_temp = (TH1D*)MRC_MC2_4[2]->Clone();//2
1131       else if(index<=9) MRC_temp = (TH1D*)MRC_MC2_4[3]->Clone();// 3
1132       else if(index==11) MRC_temp = (TH1D*)MRC_MC2_4[4]->Clone();// 4
1133       else MRC_temp = (TH1D*)MRC_MC2_4[5]->Clone();
1134     }
1135     //
1136     TERMS_4[index]->Multiply(MRC_temp);
1137   }
1138
1139   TH1D *C4raw=(TH1D*)TERMS_4[2]->Clone();
1140   C4raw->Divide(TERMS_4[12]);
1141   C4raw->GetXaxis()->SetRangeUser(0,0.2);
1142
1143
1144   TH1D *N4QS=(TH1D*)TERMS_4[0]->Clone();
1145   //
1146   N4QS->Add(TERMS_4[1], -f43);
1147   N4QS->Add(TERMS_4[2], -f43);
1148   N4QS->Add(TERMS_4[3], -f43);
1149   N4QS->Add(TERMS_4[4], -f43);
1150   //
1151   N4QS->Add(TERMS_4[5], -f42);
1152   N4QS->Add(TERMS_4[6], -f42);
1153   N4QS->Add(TERMS_4[7], -f42);
1154   N4QS->Add(TERMS_4[8], -f42);
1155   N4QS->Add(TERMS_4[9], -f42);
1156   N4QS->Add(TERMS_4[10], -f42);
1157   //
1158   N4QS->Add(TERMS_4[12], -f41);
1159   //
1160   N4QS->Scale(1/f44);
1161   //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Scale(1.005);// K scale variation
1162   N4QS->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]);
1163   //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]);// K factorization test (MC1)
1164   //N4QS->Divide(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]);// K factorization test (MC2)
1165
1166   TH1D *C4QS=(TH1D*)N4QS->Clone();
1167   C4QS->GetXaxis()->SetRangeUser(0,0.179);
1168   C4QS->SetMinimum(0.5);
1169   C4QS->SetMarkerColor(4);
1170
1171
1172   //
1173   /*TH1D *C4QS_basic=(TH1D*)TERMS_4[0]->Clone();
1174   for(int ii=1; ii<=C4QS_basic->GetNbinsX(); ii++){
1175     double value = C4QS_basic->GetBinContent(ii) - TERMS_4[12]->GetBinContent(ii);
1176     value /= f44;
1177     value += TERMS_4[12]->GetBinContent(ii);
1178     C4QS_basic->SetBinContent(ii, value);
1179     }
1180   C4QS_basic->Divide(TERMS_4[12]);
1181   //C4QS_basic->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]);
1182   C4QS_basic->SetMarkerColor(1);
1183   //C4QS_basic->Draw("same");
1184   */
1185   //
1186
1187   TH1D *C4_pieces[12];
1188   for(int i=0; i<12; i++){
1189     C4_pieces[i]=(TH1D*)TERMS_4[i]->Clone();
1190     C4_pieces[i]->Divide(TERMS_4[12]);
1191
1192     if(i==0) C4_pieces[i]->SetMarkerColor(1);
1193     else if(i<5) C4_pieces[i]->SetMarkerColor(2);
1194     else if(i<10) C4_pieces[i]->SetMarkerColor(3);
1195     else C4_pieces[i]->SetMarkerColor(4);
1196   }
1197   C4_pieces[0]->GetXaxis()->SetRangeUser(0,0.2);
1198   C4_pieces[0]->SetMinimum(0.5); C4_pieces[0]->SetMaximum(3);
1199   //C4_pieces[0]->Draw();
1200   //C4_pieces[1]->Draw("same");
1201   //C4_pieces[5]->Draw("same");
1202   //C4_pieces[11]->Draw("same");
1203   //
1204
1205   TH1D *C4_2pairs_FromSquares = (TH1D*)TERMS_4[5]->Clone();
1206   C4_2pairs_FromSquares->Divide(TERMS_4[12]);
1207   for(int bin=1; bin<=C4_2pairs_FromSquares->GetNbinsX(); bin++){
1208     double value = C4_2pairs_FromSquares->GetBinContent(bin);
1209     value -= 1;
1210     value *= value;
1211     value += 1;
1212     C4_2pairs_FromSquares->SetBinContent(bin, value);
1213   }
1214   C4_2pairs_FromSquares->SetMarkerColor(6);
1215   //C4_2pairs_FromSquares->Draw("same");
1216   //
1217   TH1D *C4QS_2pairs=(TH1D*)TERMS_4[11]->Clone();
1218   C4QS_2pairs->Add(TERMS_4[12], -pow(1-TwoFrac,2));// N1^4
1219   C4QS_2pairs->Add(TERMS_4[5], -2*TwoFrac*(1-TwoFrac));// N2 * N1^2
1220   C4QS_2pairs->Scale(1/pow(TwoFrac,2));
1221   C4QS_2pairs->Multiply(K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]);
1222   C4QS_2pairs->Divide(TERMS_4[12]);
1223   C4QS_2pairs->SetMarkerColor(6);
1224   //C4QS_2pairs->Draw("same");
1225   //
1226
1227
1228
1229   TH1D *n4QS = (TH1D*)N4QS->Clone();
1230   TH1D *c4QS = (TH1D*)N4QS->Clone();
1231   TH1D *c4QS_RemovalStage1 = (TH1D*)N4QS->Clone();
1232   TH1D *c4QS_RemovalStage2 = (TH1D*)N4QS->Clone();
1233   TH1D *c4QS_RemovalStage3 = (TH1D*)N4QS->Clone();
1234
1235
1236   for(int bin=1; bin<=n4QS->GetNbinsX(); bin++){
1237     if(n4QS->GetBinContent(bin)==0) continue;
1238     bool exitCode=kFALSE;
1239     for(int ii=0; ii<13; ii++) {if(TERMS_4[ii]->GetBinContent(bin) < 1) exitCode=kTRUE;}
1240     if(exitCode) continue;
1241     //cout<<TERMS_4[0]->GetBinContent(bin)<<endl;
1242     //
1243     double N4QSvalue = N4QS->GetBinContent(bin);
1244     double SubtractionTerm[11] = {0};
1245     double ErrorTerm[12] = {0};
1246     ErrorTerm[0] = n4QS->GetBinError(bin);
1247     //
1248     // 3-pion terms
1249     if(SameCharge || MixedCharge4pionType==1){
1250       SubtractionTerm[0] = (TERMS_4[1]->GetBinContent(bin) - f32*TERMS_4[5]->GetBinContent(bin)  - f32*TERMS_4[6]->GetBinContent(bin)  - f32*TERMS_4[8]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;// 5 6 8
1251       SubtractionTerm[0] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]->GetBinContent(bin);
1252       ErrorTerm[1] = sqrt(TERMS_4[1]->GetBinContent(bin) + f32*f32*TERMS_4[5]->GetBinContent(bin)  + f32*f32*TERMS_4[6]->GetBinContent(bin) + f32*f32*TERMS_4[8]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1253       ErrorTerm[1] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][1]->GetBinContent(bin);
1254       //
1255       SubtractionTerm[1] = (TERMS_4[2]->GetBinContent(bin) - f32*TERMS_4[5]->GetBinContent(bin)  - f32*TERMS_4[7]->GetBinContent(bin)  - f32*TERMS_4[9]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1256       SubtractionTerm[1] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][2]->GetBinContent(bin);
1257       ErrorTerm[2] = sqrt(TERMS_4[2]->GetBinContent(bin) + f32*f32*TERMS_4[5]->GetBinContent(bin)  + f32*f32*TERMS_4[7]->GetBinContent(bin) + f32*f32*TERMS_4[9]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1258       ErrorTerm[2] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][2]->GetBinContent(bin);
1259       //
1260       SubtractionTerm[2] = (TERMS_4[3]->GetBinContent(bin) - f32*TERMS_4[6]->GetBinContent(bin)  - f32*TERMS_4[7]->GetBinContent(bin)  - f32*TERMS_4[10]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1261       SubtractionTerm[2] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][3]->GetBinContent(bin);
1262       ErrorTerm[3] = sqrt(TERMS_4[3]->GetBinContent(bin) + f32*f32*TERMS_4[6]->GetBinContent(bin)  + f32*f32*TERMS_4[7]->GetBinContent(bin) + f32*f32*TERMS_4[10]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1263       ErrorTerm[3] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][3]->GetBinContent(bin);
1264       //
1265       SubtractionTerm[3] = (TERMS_4[4]->GetBinContent(bin) - f32*TERMS_4[8]->GetBinContent(bin)  - f32*TERMS_4[9]->GetBinContent(bin)  - f32*TERMS_4[10]->GetBinContent(bin) - f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1266       SubtractionTerm[3] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][4]->GetBinContent(bin);
1267       ErrorTerm[4] = sqrt(TERMS_4[4]->GetBinContent(bin) + f32*f32*TERMS_4[8]->GetBinContent(bin)  + f32*f32*TERMS_4[9]->GetBinContent(bin) + f32*f32*TERMS_4[10]->GetBinContent(bin) + f31*f31*TERMS_4[12]->GetBinContent(bin)) / f33;
1268       ErrorTerm[4] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][4]->GetBinContent(bin);
1269       //
1270     }
1271     // 2-pion terms
1272     SubtractionTerm[4] = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin);
1273     ErrorTerm[5] = sqrt( TERMS_4[5]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac;
1274     //
1275     SubtractionTerm[5] = ( TERMS_4[6]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][6]->GetBinContent(bin) / TwoFrac;
1276     ErrorTerm[6] = sqrt( TERMS_4[6]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][6]->GetBinContent(bin) / TwoFrac;
1277     //
1278     SubtractionTerm[6] = ( TERMS_4[7]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][7]->GetBinContent(bin) / TwoFrac;
1279     ErrorTerm[7] = sqrt( TERMS_4[7]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][7]->GetBinContent(bin) / TwoFrac;
1280     //
1281     SubtractionTerm[7] = ( TERMS_4[8]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][8]->GetBinContent(bin) / TwoFrac;
1282     ErrorTerm[8] = sqrt( TERMS_4[8]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][8]->GetBinContent(bin) / TwoFrac;
1283     //
1284     SubtractionTerm[8] = ( TERMS_4[9]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][9]->GetBinContent(bin) / TwoFrac;
1285     ErrorTerm[9] = sqrt( TERMS_4[9]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][9]->GetBinContent(bin) / TwoFrac;
1286     //
1287     SubtractionTerm[9] = ( TERMS_4[10]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][10]->GetBinContent(bin) / TwoFrac;
1288     ErrorTerm[10] = sqrt( TERMS_4[10]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][10]->GetBinContent(bin) / TwoFrac;
1289     //
1290
1291     //
1292     // 2 2-pion terms: cos(q12*x12)*cos(q34*x34)
1293     SubtractionTerm[10] = TERMS_4[11]->GetBinContent(bin);// N22
1294     SubtractionTerm[10] -= 2*(1-TwoFrac) * TERMS_4[5]->GetBinContent(bin);
1295     SubtractionTerm[10] += pow(1-TwoFrac,2) * TERMS_4[12]->GetBinContent(bin);
1296     SubtractionTerm[10] /= pow(TwoFrac,2);
1297     SubtractionTerm[10] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]->GetBinContent(bin);
1298     SubtractionTerm[10] *= C4muonCorrectionSC[3]->GetBinContent(bin);
1299     double intermed = ( TERMS_4[5]->GetBinContent(bin) - (1-TwoFrac)*TERMS_4[12]->GetBinContent(bin)) / TwoFrac;
1300     intermed *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin);
1301     intermed = intermed * C4muonCorrectionSC[2]->GetBinContent(bin);
1302     SubtractionTerm[10] -= 2*intermed;
1303     SubtractionTerm[10] += 2*TERMS_4[12]->GetBinContent(bin);
1304     //
1305     ErrorTerm[11] = TERMS_4[11]->GetBinContent(bin);
1306     ErrorTerm[11] += pow(1-TwoFrac,4) * TERMS_4[12]->GetBinContent(bin);
1307     ErrorTerm[11] += pow(2*(1-TwoFrac)*TwoFrac,2) * TERMS_4[5]->GetBinContent(bin);
1308     ErrorTerm[11] = sqrt(ErrorTerm[11]);
1309     ErrorTerm[11] /= pow(TwoFrac,2);
1310     ErrorTerm[11] *= K4avg[ch1_4][ch2_4][ch3_4][ch4_4][11]->GetBinContent(bin);
1311     ErrorTerm[11] = pow(ErrorTerm[11],2);
1312     ErrorTerm[11] += pow( 2*sqrt(TERMS_4[5]->GetBinContent(bin) + pow(1-TwoFrac,2)*TERMS_4[12]->GetBinContent(bin)) * K4avg[ch1_4][ch2_4][ch3_4][ch4_4][5]->GetBinContent(bin) / TwoFrac, 2);
1313     ErrorTerm[11] += 4*TERMS_4[12]->GetBinContent(bin);
1314     ErrorTerm[11] = sqrt(ErrorTerm[11]);
1315
1316     //
1317     double FinalValue[4]={0};
1318     double FinalValue_e[4]={0};
1319     if(SameCharge) {
1320       //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1321       N4QSvalue *= C4muonCorrectionSC[0]->GetBinContent(bin);
1322       ErrorTerm[0] *= C4muonCorrectionSC[0]->GetBinContent(bin);
1323     }else if(MixedCharge4pionType==1) {
1324       //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1325       N4QSvalue *= C4muonCorrectionMC1[0]->GetBinContent(bin);
1326       ErrorTerm[0] *= C4muonCorrectionMC1[0]->GetBinContent(bin);
1327     }else {
1328       //N4QSvalue = (N4QSvalue - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[0]->GetBinContent(bin) + TERMS_4[12]->GetBinContent(bin);
1329       N4QSvalue *= C4muonCorrectionMC2[0]->GetBinContent(bin);
1330       ErrorTerm[0] *= C4muonCorrectionMC2[0]->GetBinContent(bin);
1331     }
1332     for(int ii=0; ii<4; ii++) {
1333       FinalValue[ii] = N4QSvalue;
1334       FinalValue_e[ii] = pow(ErrorTerm[0],2);
1335     }
1336
1337     if(SameCharge) {
1338       //FinalValue[0] -= 4*(SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[1]->GetBinContent(bin);
1339       //FinalValue[0] += 6*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionSC[2]->GetBinContent(bin);
1340       FinalValue[0] -= 4*(SubtractionTerm[0] * C4muonCorrectionSC[1]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1341       FinalValue[0] += 6*(SubtractionTerm[4] * C4muonCorrectionSC[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1342       FinalValue[0] -= 3*(SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin));
1343
1344       FinalValue_e[0] += 4*pow(ErrorTerm[1]*C4muonCorrectionSC[1]->GetBinContent(bin),2);
1345       FinalValue_e[0] += 6*pow(ErrorTerm[5]*C4muonCorrectionSC[2]->GetBinContent(bin),2);
1346       FinalValue_e[0] += 3*pow(ErrorTerm[11]*C4muonCorrectionSC[3]->GetBinContent(bin),2);
1347       FinalValue_e[0] += (4*C4muonCorrectionSC[1]->GetBinContent(bin) + 6*C4muonCorrectionSC[2]->GetBinContent(bin) + 3*C4muonCorrectionSC[3]->GetBinContent(bin)) * TERMS_4[12]->GetBinContent(bin);
1348       //
1349       FinalValue[1] -= 6*SubtractionTerm[4] - 6*TERMS_4[12]->GetBinContent(bin);// 2-pion removal
1350       FinalValue_e[1] += 6*pow(ErrorTerm[5],2) + 6*TERMS_4[12]->GetBinContent(bin);
1351       FinalValue[2] = FinalValue[1] - 3*SubtractionTerm[10] + 3*TERMS_4[12]->GetBinContent(bin);// further 2-pair removal
1352       FinalValue_e[2] += FinalValue_e[1] + 3*pow(ErrorTerm[11],2) + 3*TERMS_4[12]->GetBinContent(bin);
1353       FinalValue[3] = SubtractionTerm[10];
1354     }else{
1355       if(MixedCharge4pionType==1 && CHARGE==-1){
1356         //FinalValue[0] -= (SubtractionTerm[0] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin);// [2] is +++ MuonCorr
1357         FinalValue[0] -= SubtractionTerm[0] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin);// [2] is +++ MuonCorr
1358         FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1359         //
1360         FinalValue_e[0] += pow(ErrorTerm[1]*C4muonCorrectionMC1[2]->GetBinContent(bin),2);
1361         FinalValue[1] -= 3*SubtractionTerm[4] - 3*TERMS_4[12]->GetBinContent(bin);
1362         FinalValue_e[1] += 3*pow(ErrorTerm[5],2) + 3*TERMS_4[12]->GetBinContent(bin);
1363       }else if(MixedCharge4pionType==1 && CHARGE==+1){
1364         //FinalValue[0] -= (SubtractionTerm[3] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC1[2]->GetBinContent(bin);
1365         FinalValue[0] -= SubtractionTerm[3] * C4muonCorrectionMC1[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin);
1366         FinalValue[0] -= 3*(SubtractionTerm[6] * C4muonCorrectionMC1[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1367         //
1368         FinalValue_e[0] += pow(ErrorTerm[4]*C4muonCorrectionMC1[2]->GetBinContent(bin),2);
1369         FinalValue[1] -= 3*SubtractionTerm[9] - 3*TERMS_4[12]->GetBinContent(bin);
1370         FinalValue_e[1] += 3*pow(ErrorTerm[10],2) + 3*TERMS_4[12]->GetBinContent(bin);
1371       }else{// --++ case
1372         //FinalValue[0] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin);// old way
1373         FinalValue[0] -= 2*(SubtractionTerm[4] * C4muonCorrectionMC2[2]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1374         FinalValue[0] -= (SubtractionTerm[10] - TERMS_4[12]->GetBinContent(bin));
1375         FinalValue[0] -= 4*(SubtractionTerm[6] * C4muonCorrectionMC2[3]->GetBinContent(bin) - TERMS_4[12]->GetBinContent(bin));
1376         FinalValue_e[0] += 2*pow(ErrorTerm[5],2) + pow(ErrorTerm[11],2) + 2*TERMS_4[12]->GetBinContent(bin);
1377         //
1378         FinalValue[1] -= 2*(SubtractionTerm[4] - TERMS_4[12]->GetBinContent(bin)) * C4muonCorrectionMC2[2]->GetBinContent(bin);
1379         FinalValue_e[1] += 2*pow(ErrorTerm[5],2) + 2*TERMS_4[12]->GetBinContent(bin);
1380       }
1381     }
1382
1383     for(int ii=0; ii<4; ii++) FinalValue_e[ii] = sqrt(FinalValue_e[ii]);
1384     //
1385
1386     n4QS->SetBinContent(bin, FinalValue[0] - TERMS_4[12]->GetBinContent(bin));
1387     n4QS->SetBinError(bin, sqrt( pow(FinalValue_e[0],2) + TERMS_4[12]->GetBinContent(bin)));
1388     c4QS->SetBinContent(bin, FinalValue[0]);
1389     c4QS->SetBinError(bin, FinalValue_e[0]);
1390     C4QS->SetBinContent(bin, N4QSvalue);
1391     C4QS->SetBinError(bin, ErrorTerm[0]);
1392     //
1393     c4QS_RemovalStage1->SetBinContent(bin, FinalValue[1]);
1394     c4QS_RemovalStage1->SetBinError(bin, FinalValue_e[1]);
1395     //
1396     c4QS_RemovalStage2->SetBinContent(bin, FinalValue[2]);
1397     c4QS_RemovalStage2->SetBinError(bin, FinalValue_e[2]);
1398     //
1399     c4QS_RemovalStage3->SetBinContent(bin, FinalValue[3]);
1400     c4QS_RemovalStage3->SetBinError(bin, FinalValue_e[3]);
1401
1402   }
1403
1404   C4QS->Divide(TERMS_4[12]);
1405   c4QS->Divide(TERMS_4[12]);
1406   c4QS->SetMarkerColor(1); c4QS->SetLineColor(1);
1407   //
1408   c4QS_RemovalStage1->Divide(TERMS_4[12]);
1409   c4QS_RemovalStage1->SetMarkerColor(2); c4QS_RemovalStage1->SetLineColor(2);
1410   c4QS_RemovalStage2->Divide(TERMS_4[12]);
1411   c4QS_RemovalStage2->SetMarkerColor(6); c4QS_RemovalStage2->SetLineColor(6);
1412   c4QS_RemovalStage3->Divide(TERMS_4[12]);
1413   c4QS_RemovalStage3->SetMarkerColor(7); c4QS_RemovalStage3->SetLineColor(7);
1414   //
1415   //
1416
1417
1418
1419   if(SameCharge) {
1420     if(CollisionType==0) C4QS->SetMaximum(7);
1421     else C4QS->SetMaximum(12);
1422   }else if(MixedCharge4pionType==1) C4QS->SetMaximum(4);
1423   else C4QS->SetMaximum(3);
1424
1425   if(CollisionType!=0) C4QS->GetXaxis()->SetRangeUser(0,0.5);
1426
1427   C4QS->Draw();
1428   c4QS_RemovalStage1->Draw("same");
1429   if(SameCharge) c4QS_RemovalStage2->Draw("same");
1430   //c4QS_RemovalStage3->Draw("same");
1431   c4QS->Draw("same");
1432
1433
1434   //cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(9)<<endl;
1435
1436   if(SameCharge && !MCcase) {
1437     //TPFullWeight_FourParticle[ch1_4]->Draw("same");
1438     FullBuildFromFits->Draw("same");
1439     PartialBuildFromFits->Draw("same");
1440   }
1441   legend3->AddEntry(C4QS,"C_{4}^{QS}","p");
1442   legend3->AddEntry(c4QS_RemovalStage1,"c_{4}^{QS}{2-pion removal}","p");
1443   if(SameCharge) legend3->AddEntry(c4QS_RemovalStage2,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1444   if(SameCharge) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair+3-pion removal}","p");
1445   if(!SameCharge && MixedCharge4pionType==1) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+3-pion removal}","p");
1446   if(!SameCharge && MixedCharge4pionType==2) legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1447   //legend3->AddEntry(c4QS,"c_{4}^{QS}{2-pion+2-pair removal}","p");
1448
1449   if(SameCharge) legend3->AddEntry(TPFullWeight_FourParticle[ch1_4],"C_{4}^{QS} built","l");
1450   legend3->Draw("same");
1451
1452   //K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->Draw();
1453
1454
1455
1456   //TH1D *c4QSFitNum = (TH1D*)c4QSFitNum2D->ProjectionY("c4QSFitNum",5,5);
1457   //TH1D *c4QSFitDen = (TH1D*)c4QSFitDen2D->ProjectionY("c4QSFitDen",5,5);
1458   //c4QSFitNum->Divide(c4QSFitDen);
1459   //for(int bin=1; bin<20; bin++){
1460     //double fc4 =
1461     //c4QSFitNum->SetBinContent(bin, (c4QSFitNum->GetBinContent(bin) - 1.0)*
1462   //}
1463   //c4QSFitNum->Draw("same");
1464
1465   //C4QS_Syst->Draw("E2 same");
1466   //C4QS->Draw("same");
1467   //C4raw->Draw();
1468   //TH1D *testTerm = (TH1D*)TERMS_4[11]->Clone();
1469   //testTerm->Divide(TERMS_4[12]);
1470   //testTerm->SetMinimum(1); testTerm->SetMaximum(1.4); testTerm->GetXaxis()->SetRangeUser(0,0.14);
1471   //testTerm->Draw();
1472
1473   ////////////////////////////////////////////////////////////////
1474   if(SameCharge && CollisionType==0){
1475     TCanvas *can4 = new TCanvas("can4", "can4",600,600,700,500);
1476     can4->SetHighLightColor(2);
1477     can4->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1478     gStyle->SetOptFit(0111);
1479     can4->SetFillColor(10);
1480     can4->SetBorderMode(0);
1481     can4->SetBorderSize(2);
1482     can4->SetGridx();
1483     can4->SetGridy();
1484     can4->SetFrameFillColor(0);
1485     can4->SetFrameBorderMode(0);
1486     can4->SetFrameBorderMode(0);
1487     gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);
1488     TLegend *legend3_2 = (TLegend*)legend1->Clone();
1489     legend3_2->SetX1(0.45); legend3_2->SetX2(0.98);  legend3_2->SetY1(0.6); legend3_2->SetY2(0.95);
1490     //
1491     //TH1D *Chi2NDF_PointSize_4 = new TH1D("Chi2NDF_PointSize_4","",40,-0.5,39.5);
1492     //TH1D *Chi2NDF_FullSize_4 = new TH1D("Chi2NDF_FullSize_4","",40,-0.5,39.5);
1493     TH1D *Chi2NDF_PointSize_4 = new TH1D("Chi2NDF_PointSize_4","",100,-0.5,99.5);
1494     TH1D *Chi2NDF_FullSize_4 = new TH1D("Chi2NDF_FullSize_4","",100,-0.5,99.5);
1495     Chi2NDF_PointSize_4->SetLineColor(4); Chi2NDF_FullSize_4->SetLineColor(2);
1496     Chi2NDF_PointSize_4->SetMarkerColor(4); Chi2NDF_FullSize_4->SetMarkerColor(2);
1497     Chi2NDF_PointSize_4->GetXaxis()->SetTitle("Coherent fraction (%)"); Chi2NDF_PointSize_4->GetYaxis()->SetTitle("#sqrt{#chi^{2}}");
1498
1499
1500     TH1D *tempDen = (TH1D*)TPFullWeight_FourParticle_2D[ch1_4]->ProjectionY("TPFullWeight4_Den", 4, 4);
1501     TH1D *tempDenNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[ch1_4]->ProjectionY("TPNegFullWeight4_Den", 4, 4);
1502     tempDen->Add(tempDenNeg);// Add Pos and Neg Weight
1503
1504     for(int binG=5; binG<=104; binG++){// 44
1505       TString *proName=new TString("TPFullWeight4_");
1506       *proName += binG;
1507       TH1D *tempNum = (TH1D*)TPFullWeight_FourParticle_2D[ch1_4]->ProjectionY(proName->Data(), binG, binG);
1508       proName->Append("_Neg");
1509       TH1D *tempNumNeg = (TH1D*)TPNegFullWeight_FourParticle_2D[ch1_4]->ProjectionY(proName->Data(), binG, binG);
1510       //
1511       // Add Pos and Neg Weights
1512       tempNum->Add(tempNumNeg);
1513       //
1514       //tempNum->Add(tempDen);// now added in InterpCorr section
1515       tempNum->Divide(tempDen);
1516
1517
1518       double Chi2=0;
1519       double NDF=0;
1520       for(int binQ4=4; binQ4<=4; binQ4++){
1521         if(tempNum->GetBinContent(binQ4) <=0) continue;
1522         double value = C4QS->GetBinContent(binQ4) - tempNum->GetBinContent(binQ4);
1523         double err = pow(C4QS->GetBinError(binQ4),2);
1524
1525         err = sqrt(err);
1526         if(err<=0) continue;
1527
1528         Chi2 += pow(value / err,2);
1529         //
1530         NDF += 1;
1531       }
1532       //if(binG<25) {Chi2NDF_PointSize_4->SetBinContent(1 + 2*(binG-5), sqrt(fabs(Chi2))/NDF); Chi2NDF_PointSize_4->SetBinError(1 + 2*(binG-5), 0.001);}
1533       //else {Chi2NDF_FullSize_4->SetBinContent(1 + 2*(binG-25), sqrt(fabs(Chi2))/NDF); Chi2NDF_FullSize_4->SetBinError(1 + 2*(binG-25), 0.001);}
1534       if(binG<55) {Chi2NDF_PointSize_4->SetBinContent(1 + 2*(binG-5), sqrt(fabs(Chi2))/NDF); Chi2NDF_PointSize_4->SetBinError(1 + 2*(binG-5), 0.001);}
1535       else {Chi2NDF_FullSize_4->SetBinContent(1 + 2*(binG-55), sqrt(fabs(Chi2))/NDF); Chi2NDF_FullSize_4->SetBinError(1 + 2*(binG-55), 0.001);}
1536     }
1537     Chi2NDF_PointSize_4->SetMarkerStyle(20); Chi2NDF_FullSize_4->SetMarkerStyle(20);
1538
1539
1540     Chi2NDF_PointSize_4->Draw();
1541     Chi2NDF_FullSize_4->Draw("same");
1542     legend3_2->AddEntry(Chi2NDF_PointSize_4,"R_{coh}=0 (Point Source)","p");
1543     legend3_2->AddEntry(Chi2NDF_FullSize_4,"R_{coh}=R_{ch}","p");
1544     legend3_2->Draw("same");
1545
1546
1547
1548   }
1549
1550   // Print 4-pion points
1551   for(int bin=1; bin<=12; bin++){
1552     //cout<<C4QS->GetBinContent(bin)<<", ";
1553     //cout<<c4QS->GetBinContent(bin)<<", ";
1554     //cout<<TPFullWeight_FourParticle[ch1_4]->GetBinContent(bin)<<", ";
1555     //cout<<C4raw->GetBinContent(bin)<<", ";
1556     //cout<<K4avg[ch1_4][ch2_4][ch3_4][ch4_4][0]->GetBinContent(bin)<<", ";
1557   }
1558   cout<<endl;
1559   for(int bin=1; bin<=12; bin++){
1560     //cout<<c4QS->GetBinContent(bin)<<", ";
1561     //cout<<C4QS->GetBinError(bin)<<", ";
1562     //cout<<c4QS->GetBinError(bin)<<", ";
1563     //cout<<C4raw->GetBinError(bin)<<", ";
1564   }
1565   //cout<<endl;
1566   ////////////////////////////////////////////////////////////////
1567   // r4
1568
1569   TF1 *ChaoticLimit_r42 = new TF1("ChaoticLimit_r42","6",0,10);
1570   ChaoticLimit_r42->SetLineStyle(2);
1571   TH1D *r42;
1572   if(SameCharge && CollisionType==0){
1573     /*TCanvas *can5 = new TCanvas("can5", "can5",1200,600,700,500);
1574     can5->SetHighLightColor(2);
1575     can5->Range(-0.7838115,-1.033258,0.7838115,1.033258);
1576     gStyle->SetOptFit(0111);
1577     can5->SetFillColor(10);
1578     can5->SetBorderMode(0);
1579     can5->SetBorderSize(2);
1580     can5->SetGridx();
1581     can5->SetGridy();
1582     can5->SetFrameFillColor(0);
1583     can5->SetFrameBorderMode(0);
1584     can5->SetFrameBorderMode(0);
1585     gPad->SetRightMargin(0.02); gPad->SetTopMargin(0.02);*/
1586
1587     r42 = (TH1D*)n4QS->Clone();
1588     TPN_FourParticle[ch1_4]->Multiply(MRC_SC_4[4]);
1589     r42->Divide(TPN_FourParticle[ch1_4]);
1590     r42->GetXaxis()->SetRangeUser(0,0.08);
1591     r42->SetMinimum(0.5); r42->SetMaximum(14);
1592     r42->GetYaxis()->SetTitle("r_{4}{2}");
1593     //
1594     //r42->Draw();
1595     //ChaoticLimit_r42->Draw("same");
1596     //TPN_FourParticle[ch1_4]->Draw();
1597     //fTPNRejects4pion1->SetLineColor(2);
1598     //fTPNRejects4pion1->Draw("same");
1599   }
1600
1601   //double EA = exp(-pow(0.005 * 10/FmToGeV,2)/2.);
1602   //TF1 *C2mag = new TF1("C2mag","pow(1-x,2)*pow([0],2) + 2*x*(1-x)*[0]",0,1);
1603   //C2mag->FixParameter(0, EA);
1604   //TF1 *C4norm = new TF1("C4norm","( 6*pow(1-x,4)*pow([0],4) + 24*x*pow(1-x,3)*pow([0],3) + 4*(2*pow(1-x,3)*pow([0],3) + 6*x*pow(1-x,2)*pow([0],2)) + 3*pow(C2mag,2) + 6*C2mag + 1) / (6*pow(C2mag,2) + 8*pow(C2mag,1.5) + 3*pow(C2mag,2) + 6*C2mag + 1)",0,1);
1605   //C4norm->FixParameter(0, EA);
1606   //C4norm->Draw();
1607
1608   /*TF1 *C2mag = new TF1("C2mag","pow(1-[0],2)*exp(-pow(x*[1],2)/6.) + 2*[0]*(1-[0])*exp(-pow(x*[1],2)/12.)",0,0.2);
1609   C2mag->FixParameter(0, 0.25);
1610   C2mag->FixParameter(1, 10./FmToGeV);
1611   TF1 *C4norm = new TF1("C4norm","( 6*pow(1-[0],4)*exp(-pow(x*[1],2)/3.) + 24*[0]*pow(1-[0],3)*exp(-pow(x*[1],2)/4.) + 4*(2*pow(1-[0],3)*exp(-pow(x*[1],2)/4.) + 6*[0]*pow(1-[0],2)*exp(-pow(x*[1],2)/6.)) + 3*pow(C2mag,2) + 6*C2mag + 1) / (6*pow(C2mag,2) + 8*pow(C2mag,1.5) + 3*pow(C2mag,2) + 6*C2mag + 1)",0,0.2);
1612   C4norm->FixParameter(0, 0.25);
1613   C4norm->FixParameter(1, 10./FmToGeV);*/
1614   //C4norm->Draw();
1615
1616
1617   if(SaveToFile){
1618     TString *savefilename = new TString("OutFiles/OutFile");
1619     if(MCcase) savefilename->Append("MonteCarlo");
1620     //
1621     if(SameCharge) savefilename->Append("SC");
1622     else savefilename->Append("MC");
1623     //
1624     if(!SameCharge) *savefilename += MixedCharge4pionType_def;
1625     //
1626     if(CHARGE==1) savefilename->Append("_Pos_");
1627     else savefilename->Append("_Neg_");
1628     //
1629
1630     savefilename->Append("KT_");
1631     *savefilename += EDbin+1;
1632
1633     savefilename->Append("_M");
1634     *savefilename += Mbin;
1635     savefilename->Append(".root");
1636     TFile *savefile = new TFile(savefilename->Data(),"RECREATE");
1637     //
1638     C2->Write("C2");
1639     C2QS->Write("C2QS");
1640     C3QS->Write("C3QS");
1641     c3QS->Write("c3QS");
1642     C4QS->Write("C4QS");
1643     c4QS->Write("c4QS");
1644     c4QS_RemovalStage1->Write("c4QS_RemovalStage1");
1645     if(SameCharge) {
1646       r3->Write("r3");
1647       r42->Write("r42");
1648       c4QS_RemovalStage2->Write("c4QS_RemovalStage2");
1649       TPFullWeight_ThreeParticle[ch1_3]->Write("C3QS_built");
1650       TPFullWeight_FourParticle[ch1_4]->Write("C4QS_built");
1651       //
1652       TPFullWeight_ThreeParticle_2D[ch1_3]->Write("C3QS_built2D");
1653       TPFullWeight_FourParticle_2D[ch1_4]->Write("C4QS_built2D");
1654       TPNegFullWeight_ThreeParticle_2D[ch1_3]->Write("C3QS_Negbuilt2D");
1655       TPNegFullWeight_FourParticle_2D[ch1_4]->Write("C4QS_Negbuilt2D");
1656     }
1657     //
1658     savefile->Close();
1659
1660   }
1661
1662
1663 }
1664
1665 //________________________________________________________________________
1666 void SetFSICorrelations(){
1667   // read in 2-particle and 3-particle FSI correlations = K2 & K3
1668   // 2-particle input histo from file is binned in qinv.  3-particle in qinv of each pair
1669   TFile *fsifile = new TFile("KFile.root","READ");
1670   if(!fsifile->IsOpen()) {
1671     cout<<"No FSI file found!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
1672   }
1673
1674   TH1D *temphistoSS[12];
1675   TH1D *temphistoOS[12];
1676   for(Int_t MB=0; MB<12; MB++) {
1677     TString *nameK2SS = new TString("K2ss_");
1678     *nameK2SS += MB;
1679     temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
1680     //
1681     TString *nameK2OS = new TString("K2os_");
1682     *nameK2OS += MB;
1683     temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
1684     //
1685     fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
1686     fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
1687     fFSIss[MB]->SetDirectory(0);
1688     fFSIos[MB]->SetDirectory(0);
1689   }
1690   //
1691
1692   fsifile->Close();
1693
1694 }
1695
1696 double Gamov(int chargeProduct, double qinv){
1697
1698   double arg = chargeProduct*2.*PI/(BohrR*qinv/2.);
1699
1700   return arg/(exp(arg)-1);
1701 }
1702
1703 void SetMomResCorrections(){
1704
1705   TString *momresfilename = new TString("MomResFile");
1706   if(FileSetting==7) momresfilename->Append("_10percentIncrease");
1707   if(CollisionType!=0) momresfilename->Append("_ppAndpPb");
1708   momresfilename->Append(".root");
1709
1710   TFile *MomResFile = new TFile(momresfilename->Data(),"READ");
1711   TString *proName[28];
1712   for(int ii=0; ii<28; ii++){
1713     proName[ii] = new TString("MRC_pro_");
1714     *proName[ii] += ii;
1715   }
1716
1717   MRC_SC_2[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_SC_term1"))->ProjectionY(proName[0]->Data(), RbinMRC, RbinMRC));
1718   MRC_SC_2[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_SC_term2"))->ProjectionY(proName[1]->Data(), RbinMRC, RbinMRC));
1719   MRC_SC_2[0]->SetDirectory(0); MRC_SC_2[1]->SetDirectory(0);
1720   //
1721   MRC_MC_2[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_MC_term1"))->ProjectionY(proName[2]->Data(), RbinMRC, RbinMRC));
1722   MRC_MC_2[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_2_MC_term2"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC));
1723   MRC_MC_2[0]->SetDirectory(0); MRC_MC_2[1]->SetDirectory(0);
1724   //
1725   MRC_SC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term1"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
1726   MRC_SC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term2"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC));
1727   MRC_SC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_SC_term3"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC));
1728   MRC_SC_3[0]->SetDirectory(0); MRC_SC_3[1]->SetDirectory(0); MRC_SC_3[2]->SetDirectory(0);
1729   //
1730   MRC_MC_3[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term1"))->ProjectionY(proName[7]->Data(), RbinMRC, RbinMRC));
1731   MRC_MC_3[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term2"))->ProjectionY(proName[8]->Data(), RbinMRC, RbinMRC));
1732   MRC_MC_3[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term3"))->ProjectionY(proName[9]->Data(), RbinMRC, RbinMRC));
1733   MRC_MC_3[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_3_MC_term4"))->ProjectionY(proName[10]->Data(), RbinMRC, RbinMRC));
1734   MRC_MC_3[0]->SetDirectory(0); MRC_MC_3[1]->SetDirectory(0); MRC_MC_3[2]->SetDirectory(0); MRC_MC_3[3]->SetDirectory(0);
1735   //
1736   MRC_SC_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term1"))->ProjectionY(proName[11]->Data(), RbinMRC, RbinMRC));
1737   MRC_SC_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term2"))->ProjectionY(proName[12]->Data(), RbinMRC, RbinMRC));
1738   MRC_SC_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term3"))->ProjectionY(proName[13]->Data(), RbinMRC, RbinMRC));
1739   MRC_SC_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term4"))->ProjectionY(proName[14]->Data(), RbinMRC, RbinMRC));
1740   MRC_SC_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_SC_term5"))->ProjectionY(proName[15]->Data(), RbinMRC, RbinMRC));
1741   MRC_SC_4[0]->SetDirectory(0); MRC_SC_4[1]->SetDirectory(0); MRC_SC_4[2]->SetDirectory(0); MRC_SC_4[3]->SetDirectory(0);
1742   MRC_SC_4[4]->SetDirectory(0);
1743   //
1744   MRC_MC1_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term1"))->ProjectionY(proName[16]->Data(), RbinMRC, RbinMRC));
1745   MRC_MC1_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term2"))->ProjectionY(proName[17]->Data(), RbinMRC, RbinMRC));
1746   MRC_MC1_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term3"))->ProjectionY(proName[18]->Data(), RbinMRC, RbinMRC));
1747   MRC_MC1_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term4"))->ProjectionY(proName[19]->Data(), RbinMRC, RbinMRC));
1748   MRC_MC1_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term5"))->ProjectionY(proName[20]->Data(), RbinMRC, RbinMRC));
1749   MRC_MC1_4[5] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC1_term6"))->ProjectionY(proName[21]->Data(), RbinMRC, RbinMRC));
1750   MRC_MC1_4[0]->SetDirectory(0); MRC_MC1_4[1]->SetDirectory(0); MRC_MC1_4[2]->SetDirectory(0); MRC_MC1_4[3]->SetDirectory(0);
1751   MRC_MC1_4[4]->SetDirectory(0); MRC_MC1_4[5]->SetDirectory(0);
1752   //
1753   MRC_MC2_4[0] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term1"))->ProjectionY(proName[22]->Data(), RbinMRC, RbinMRC));
1754   MRC_MC2_4[1] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term2"))->ProjectionY(proName[23]->Data(), RbinMRC, RbinMRC));
1755   MRC_MC2_4[2] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term3"))->ProjectionY(proName[24]->Data(), RbinMRC, RbinMRC));
1756   MRC_MC2_4[3] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term4"))->ProjectionY(proName[25]->Data(), RbinMRC, RbinMRC));
1757   MRC_MC2_4[4] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term5"))->ProjectionY(proName[26]->Data(), RbinMRC, RbinMRC));
1758   MRC_MC2_4[5] = (TH1D*)(((TH2D*)MomResFile->Get("MRC_4_MC2_term6"))->ProjectionY(proName[27]->Data(), RbinMRC, RbinMRC));
1759   MRC_MC2_4[0]->SetDirectory(0); MRC_MC2_4[1]->SetDirectory(0); MRC_MC2_4[2]->SetDirectory(0); MRC_MC2_4[3]->SetDirectory(0);
1760   MRC_MC2_4[4]->SetDirectory(0); MRC_MC2_4[5]->SetDirectory(0);
1761
1762   if(!MRC || MCcase_def){
1763     for(int bin=1; bin<=MRC_SC_2[0]->GetNbinsX(); bin++){
1764       MRC_SC_2[0]->SetBinContent(bin, 1.0); MRC_SC_2[1]->SetBinContent(bin, 1.0);
1765       MRC_MC_2[0]->SetBinContent(bin, 1.0); MRC_MC_2[1]->SetBinContent(bin, 1.0);
1766     }
1767     for(int bin=1; bin<=MRC_SC_3[0]->GetNbinsX(); bin++){
1768       MRC_SC_3[0]->SetBinContent(bin, 1.0); MRC_SC_3[1]->SetBinContent(bin, 1.0); MRC_SC_3[2]->SetBinContent(bin, 1.0);
1769       MRC_MC_3[0]->SetBinContent(bin, 1.0); MRC_MC_3[1]->SetBinContent(bin, 1.0); MRC_MC_3[2]->SetBinContent(bin, 1.0);
1770       MRC_MC_3[3]->SetBinContent(bin, 1.0);
1771     }
1772     for(int bin=1; bin<=MRC_SC_4[0]->GetNbinsX(); bin++){
1773       MRC_SC_4[0]->SetBinContent(bin, 1.0); MRC_SC_4[1]->SetBinContent(bin, 1.0); MRC_SC_4[2]->SetBinContent(bin, 1.0);
1774       MRC_SC_4[3]->SetBinContent(bin, 1.0); MRC_SC_4[4]->SetBinContent(bin, 1.0);
1775       MRC_MC1_4[0]->SetBinContent(bin, 1.0); MRC_MC1_4[1]->SetBinContent(bin, 1.0); MRC_MC1_4[2]->SetBinContent(bin, 1.0);
1776       MRC_MC1_4[3]->SetBinContent(bin, 1.0); MRC_MC1_4[4]->SetBinContent(bin, 1.0); MRC_MC1_4[5]->SetBinContent(bin, 1.0);
1777       //
1778       MRC_MC2_4[0]->SetBinContent(bin, 1.0); MRC_MC2_4[1]->SetBinContent(bin, 1.0); MRC_MC2_4[2]->SetBinContent(bin, 1.0);
1779       MRC_MC2_4[3]->SetBinContent(bin, 1.0); MRC_MC2_4[4]->SetBinContent(bin, 1.0); MRC_MC2_4[5]->SetBinContent(bin, 1.0);
1780     }
1781   }
1782   MomResFile->Close();
1783
1784 }
1785
1786
1787 double cubicInterpolate (double p[4], double x) {
1788   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
1789 }
1790 double bicubicInterpolate (double p[4][4], double x, double y) {
1791         double arr[4];
1792         arr[0] = cubicInterpolate(p[0], y);
1793         arr[1] = cubicInterpolate(p[1], y);
1794         arr[2] = cubicInterpolate(p[2], y);
1795         arr[3] = cubicInterpolate(p[3], y);
1796         return cubicInterpolate(arr, x);
1797 }
1798 double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {
1799         double arr[4];
1800         arr[0] = bicubicInterpolate(p[0], y, z);
1801         arr[1] = bicubicInterpolate(p[1], y, z);
1802         arr[2] = bicubicInterpolate(p[2], y, z);
1803         arr[3] = bicubicInterpolate(p[3], y, z);
1804         return cubicInterpolate(arr, x);
1805 }
1806 float fact(float n){
1807   return (n < 1.00001) ? 1 : fact(n - 1) * n;
1808 }
1809 //________________________________________________________________________________________
1810 void SetMuonCorrections(){
1811   TString *name = new TString("MuonCorrection");
1812   if(FileSetting==8) name->Append("_92percent");
1813   if(CollisionType!=0) name->Append("_ppAndpPb");
1814   name->Append(".root");
1815   TFile *MuonFile=new TFile(name->Data(),"READ");
1816   TString *proName[22];
1817   for(int ii=0; ii<22; ii++){
1818     proName[ii] = new TString("MuonCorr_pro_");
1819     *proName[ii] += ii;
1820   }
1821   C2muonCorrectionSC = (TH1D*)(((TH2D*)MuonFile->Get("C2muonCorrection_SC"))->ProjectionY(proName[0]->Data(), RbinMRC, RbinMRC));
1822   C2muonCorrectionMC = (TH1D*)(((TH2D*)MuonFile->Get("C2muonCorrection_MC"))->ProjectionY(proName[1]->Data(), RbinMRC, RbinMRC));
1823   WeightmuonCorrection = (TH1D*)(((TH2D*)MuonFile->Get("WeightmuonCorrection"))->ProjectionY(proName[2]->Data(), RbinMRC, RbinMRC));
1824   C2muonCorrectionSC->SetDirectory(0); C2muonCorrectionMC->SetDirectory(0); WeightmuonCorrection->SetDirectory(0);
1825   //
1826   C3muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term1"))->ProjectionY(proName[3]->Data(), RbinMRC, RbinMRC));
1827   C3muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_SC_term2"))->ProjectionY(proName[4]->Data(), RbinMRC, RbinMRC));
1828   C3muonCorrectionMC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term1"))->ProjectionY(proName[5]->Data(), RbinMRC, RbinMRC));
1829   C3muonCorrectionMC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term2"))->ProjectionY(proName[6]->Data(), RbinMRC, RbinMRC));
1830   C3muonCorrectionMC[2] = (TH1D*)(((TH2D*)MuonFile->Get("C3muonCorrection_MC_term3"))->ProjectionY(proName[7]->Data(), RbinMRC, RbinMRC));
1831   C3muonCorrectionSC[0]->SetDirectory(0); C3muonCorrectionSC[1]->SetDirectory(0);
1832   C3muonCorrectionMC[0]->SetDirectory(0); C3muonCorrectionMC[1]->SetDirectory(0); C3muonCorrectionMC[2]->SetDirectory(0);
1833   //
1834   C4muonCorrectionSC[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term1"))->ProjectionY(proName[8]->Data(), RbinMRC, RbinMRC));
1835   C4muonCorrectionSC[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term2"))->ProjectionY(proName[9]->Data(), RbinMRC, RbinMRC));
1836   C4muonCorrectionSC[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term3"))->ProjectionY(proName[10]->Data(), RbinMRC, RbinMRC));
1837   C4muonCorrectionSC[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_SC_term4"))->ProjectionY(proName[11]->Data(), RbinMRC, RbinMRC));
1838   C4muonCorrectionSC[0]->SetDirectory(0); C4muonCorrectionSC[1]->SetDirectory(0);
1839   C4muonCorrectionSC[2]->SetDirectory(0); C4muonCorrectionSC[3]->SetDirectory(0);
1840   //
1841   C4muonCorrectionMC1[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term1"))->ProjectionY(proName[12]->Data(), RbinMRC, RbinMRC));
1842   C4muonCorrectionMC1[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term2"))->ProjectionY(proName[13]->Data(), RbinMRC, RbinMRC));
1843   C4muonCorrectionMC1[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term3"))->ProjectionY(proName[14]->Data(), RbinMRC, RbinMRC));
1844   C4muonCorrectionMC1[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term4"))->ProjectionY(proName[15]->Data(), RbinMRC, RbinMRC));
1845   C4muonCorrectionMC1[4] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC1_term5"))->ProjectionY(proName[16]->Data(), RbinMRC, RbinMRC));
1846   C4muonCorrectionMC1[0]->SetDirectory(0); C4muonCorrectionMC1[1]->SetDirectory(0);
1847   C4muonCorrectionMC1[2]->SetDirectory(0); C4muonCorrectionMC1[3]->SetDirectory(0); C4muonCorrectionMC1[4]->SetDirectory(0);
1848   //
1849   C4muonCorrectionMC2[0] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term1"))->ProjectionY(proName[17]->Data(), RbinMRC, RbinMRC));
1850   C4muonCorrectionMC2[1] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term2"))->ProjectionY(proName[18]->Data(), RbinMRC, RbinMRC));
1851   C4muonCorrectionMC2[2] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term3"))->ProjectionY(proName[19]->Data(), RbinMRC, RbinMRC));
1852   C4muonCorrectionMC2[3] = (TH1D*)(((TH2D*)MuonFile->Get("C4muonCorrection_MC2_term4"))->ProjectionY(proName[20]->Data(), RbinMRC, RbinMRC));
1853   C4muonCorrectionMC2[4] = (TH1D*)C4muonCorrectionSC[3]->Clone();
1854   C4muonCorrectionMC2[0]->SetDirectory(0); C4muonCorrectionMC2[1]->SetDirectory(0);
1855   C4muonCorrectionMC2[2]->SetDirectory(0); C4muonCorrectionMC2[3]->SetDirectory(0); C4muonCorrectionMC2[4]->SetDirectory(0);
1856   //
1857   //
1858   if(!MuonCorrection || MCcase_def){
1859     for(int bin=1; bin<=C2muonCorrectionSC->GetNbinsX(); bin++){
1860       C2muonCorrectionSC->SetBinContent(bin, 1.0);
1861       C2muonCorrectionMC->SetBinContent(bin, 1.0);
1862       WeightmuonCorrection->SetBinContent(bin, 1.0);
1863     }
1864     for(int bin=1; bin<=C3muonCorrectionSC[0]->GetNbinsX(); bin++){
1865       C3muonCorrectionSC[0]->SetBinContent(bin, 1.0); C3muonCorrectionSC[1]->SetBinContent(bin, 1.0);
1866       C3muonCorrectionMC[0]->SetBinContent(bin, 1.0); C3muonCorrectionMC[1]->SetBinContent(bin, 1.0); C3muonCorrectionMC[2]->SetBinContent(bin, 1.0);
1867     }
1868     for(int bin=1; bin<=C4muonCorrectionSC[0]->GetNbinsX(); bin++){
1869       C4muonCorrectionSC[0]->SetBinContent(bin, 1.0); C4muonCorrectionSC[1]->SetBinContent(bin, 1.0); C4muonCorrectionSC[2]->SetBinContent(bin, 1.0);
1870       C4muonCorrectionSC[3]->SetBinContent(bin, 1.0);
1871       C4muonCorrectionMC1[0]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[1]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[2]->SetBinContent(bin, 1.0);
1872       C4muonCorrectionMC1[3]->SetBinContent(bin, 1.0); C4muonCorrectionMC1[4]->SetBinContent(bin, 1.0);
1873       //
1874       C4muonCorrectionMC2[0]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[1]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[2]->SetBinContent(bin, 1.0);
1875       C4muonCorrectionMC2[3]->SetBinContent(bin, 1.0); C4muonCorrectionMC2[4]->SetBinContent(bin, 1.0);
1876     }
1877   }
1878
1879   MuonFile->Close();
1880 }
1881 //________________________________________________________________________
1882 void SetFSIindex(Float_t R){
1883   if(!MCcase_def){
1884     if(CollisionType==0){
1885       if(Mbin_def==0) fFSIindex = 0;//0-5%
1886       else if(Mbin_def==1) fFSIindex = 1;//5-10%
1887       else if(Mbin_def<=3) fFSIindex = 2;//10-20%
1888       else if(Mbin_def<=5) fFSIindex = 3;//20-30%
1889       else if(Mbin_def<=7) fFSIindex = 4;//30-40%
1890       else if(Mbin_def<=9) fFSIindex = 5;//40-50%
1891       else if(Mbin_def<=12) fFSIindex = 6;//40-50%
1892       else if(Mbin_def<=15) fFSIindex = 7;//40-50%
1893       else if(Mbin_def<=18) fFSIindex = 8;//40-50%
1894       else fFSIindex = 8;//90-100%
1895     }else fFSIindex = 9;// pp and pPb
1896   }else{// FSI binning for MC
1897     if(R>=10.) fFSIindex = 0;
1898     else if(R>=9.) fFSIindex = 1;
1899     else if(R>=8.) fFSIindex = 2;
1900     else if(R>=7.) fFSIindex = 3;
1901     else if(R>=6.) fFSIindex = 4;
1902     else if(R>=5.) fFSIindex = 5;
1903     else if(R>=4.) fFSIindex = 6;
1904     else if(R>=3.) fFSIindex = 7;
1905     else if(R>=2.) fFSIindex = 8;
1906     else fFSIindex = 9;
1907   }
1908 }
1909 //________________________________________________________________________
1910 Float_t FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
1911   // returns 2-particle Coulomb correlations = K2
1912   Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
1913   Int_t qbinH = qbinL+1;
1914   if(qbinL <= 0) return 1.0;
1915   if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) {// Scaled Gamov approximation
1916     int chargeproduct = 1;
1917     if(charge1!=charge2) {
1918       chargeproduct = -1;
1919       Float_t ScaleFac = (fFSIos[fFSIindex]->GetBinContent(fFSIos[fFSIindex]->GetNbinsX()-1) - 1);
1920       ScaleFac /= (Gamov(chargeproduct, fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(fFSIos[fFSIindex]->GetNbinsX()-1)) - 1);
1921       return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1);
1922     }else{
1923       Float_t ScaleFac = (fFSIss[fFSIindex]->GetBinContent(fFSIss[fFSIindex]->GetNbinsX()-1) - 1);
1924       ScaleFac /= (Gamov(chargeproduct, fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(fFSIss[fFSIindex]->GetNbinsX()-1)) - 1);
1925       return ( (Gamov(chargeproduct, qinv)-1)*ScaleFac + 1);
1926     }
1927   }
1928
1929   Float_t slope=0;
1930   if(charge1==charge2){
1931     slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
1932     slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
1933     return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
1934   }else {
1935     slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
1936     slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
1937     return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
1938   }
1939 }