#include <fstream> 
#include <string>
#include <vector>
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include <cmath>
#include <complex>

struct dPointsEffs_t {
  double val;
  double err;
  double val_2;
  double err_2;
};

dPointsEffs_t getEff1530(int bNumber, double rangeLow, double rangeHigh);
dPointsEffs_t getEff1820(int bNumber, int kmLamNumber, double rangeLow, double rangeHigh);
dPointsEffs_t getEff2030(int bNumber, int kmLamNumber, double rangeLow, double rangeHigh);

// Global mass constants for the final states (in GeV)
const double m_xi  = 1.31486; // Xi-
const double m_pi  = 0.13957; // pi+
const double m_lam = 1.11568; // Lambda
const double m_k   = 0.49368; // K-

// Relativistic two-body phase space factor
double RelativisticRho(double m, double m1, double m2) {
    if (m <= (m1 + m2)) return 0.0; // Return 0 below threshold for real width
    double term1 = 1.0 - (std::pow(m1 + m2, 2) / std::pow(m, 2));
    double term2 = 1.0 - (std::pow(m1 - m2, 2) / std::pow(m, 2));
    return std::sqrt(term1 * term2);
}

// Complex-valued phase space to handle analytic continuation below threshold
std::complex<double> ComplexRho(double m, double m1, double m2) {
    double term1 = 1.0 - (std::pow(m1 + m2, 2) / std::pow(m, 2));
    double term2 = 1.0 - (std::pow(m1 - m2, 2) / std::pow(m, 2));
    // scimath-like behavior using std::complex sqrt
    return std::sqrt(std::complex<double>(term1 * term2, 0.0));
}
// Common denominator calculation shared by both channels
std::complex<double> GetFlatteDenominator(double m, double m0, double g1, double g2) {
    std::complex<double> rho1 = ComplexRho(m, m_xi, m_pi);
    std::complex<double> rho2 = ComplexRho(m, m_lam, m_k);
    
    std::complex<double> i(0.0, 1.0);
    return std::pow(m0, 2) - std::pow(m, 2) - i * (g1 * rho1 + g2 * rho2);
}

// Channel 1 Intensity: Xi + pi
double Intensity_XiPi(double m, double A, double m0, double g1, double g2) {
    double rho1 = RelativisticRho(m, m_xi, m_pi);
    std::complex<double> denom = GetFlatteDenominator(m, m0, g1, g2);
    
    double numerator = std::pow(A, 2) * g1 * rho1;
    return numerator / std::norm(denom); // std::norm returns |z|^2
}
// Channel 2 Intensity: Lambda + K
double Intensity_LamK(double m, double A, double m0, double g1, double g2) {
    double rho2 = RelativisticRho(m, m_xi, m_pi); // Note: Keep channel 1 open if m > threshold
    double rho2_actual = RelativisticRho(m, m_lam, m_k); 
    std::complex<double> denom = GetFlatteDenominator(m, m0, g1, g2);
    
    // The numerator evaluates to zero below 1.610 GeV automatically via rho2_actual
    double numerator = std::pow(A, 2) * g2 * rho2_actual;
    return numerator / std::norm(denom);
}



class  MyFitFun1 {
public:
  double Evaluate(double *xPtr, double *parPtr) {
    double fitVal = 0.0;
    double xVal = *xPtr;

    double smearVal = 0.00926;
    //double smearVal = 0.0071;

    smearVal = parPtr[0];
    bool kmLamType = false;
    if (smearVal < 0) kmLamType = true;
    smearVal = fabs(smearVal);

    smearVal = (11.7*xVal - 7.8)/1000;
    if (kmLamType == true) smearVal = (16.8*xVal - 23)/1000;
    
    double mag1   = fabs(parPtr[1]);
    double center1 = parPtr[2];
    double gamma1  = parPtr[3];
    double bw1 = TMath::Voigt(xVal-center1,smearVal,gamma1)*mag1;

    double mag2   = fabs(parPtr[4]);
    double center2 = parPtr[5];
    double gamma2  = parPtr[6];
    double bw2 = TMath::Voigt(xVal-center2,smearVal,gamma2)*mag2;

    double mag3   = fabs(parPtr[7]);
    double center3 = parPtr[8];
    double gamma3  = parPtr[9];
    double bw3 = TMath::Voigt(xVal-center3,smearVal,gamma3)*mag3;

    double mag4   = fabs(parPtr[10]);
    double center4 = parPtr[11];
    double gamma4  = parPtr[12];
    double bw4 = TMath::Voigt(xVal-center4,smearVal,gamma4)*mag4;

    double mag5   = fabs(parPtr[13]);
    double center5 = parPtr[14];
    double gamma5  = parPtr[15];
    double bw5 = TMath::Voigt(xVal-center5,smearVal,gamma5)*mag5;

    double mag6   = fabs(parPtr[16]);
    double center6 = parPtr[17];
    double gamma6  = parPtr[18];
    double bw6 = TMath::Voigt(xVal-center6,smearVal,gamma6)*mag6;

    double mag7   = fabs(parPtr[19]);
    double center7 = parPtr[20];
    double gamma7  = parPtr[21];
    double bw7 = TMath::Voigt(xVal-center7,smearVal,gamma7)*mag7;

    double mag8   = fabs(parPtr[22]);
    double center8 = parPtr[23];
    double gamma8  = parPtr[24];
    double bw8 = TMath::Voigt(xVal-center8,smearVal,gamma8)*mag8;

    int bkgPar0 = 25;
    double bkg = parPtr[bkgPar0] + parPtr[bkgPar0+1]*xVal + parPtr[bkgPar0+2]*pow(xVal,2);
    double sCenter1 = parPtr[bkgPar0+3];
    double sWidth1 = parPtr[bkgPar0+4];
    double sFun1 = 1.0/(exp((-xVal + sCenter1)*sWidth1) + 1);

    double bkgDalitz = (parPtr[bkgPar0] + parPtr[bkgPar0+1]*xVal +
			parPtr[bkgPar0+2]*pow(xVal,2) + parPtr[bkgPar0+3]*pow(xVal,3))*parPtr[bkgPar0+4];

    ////FLATTE////
    //Some constants
    double mK = 0.493677;
    double mLam = 1.115683;
    double mCasM = 1.32170;
    double mPi0 = 0.1349768;
    
    int flattePar0 = bkgPar0+5;

    double flatteMag = parPtr[flattePar0];
    double gCasPi    = parPtr[flattePar0+1];
    //double gKLam     = gCasPi*parPtr[flattePar0+2];
    double gKLam     = parPtr[flattePar0+2];
    double centerFlatte = parPtr[flattePar0+3];

    double GammaKLam  = gKLam*sqrt(fabs(0.25*xVal*xVal - 0.25*(mK+mLam)*(mK+mLam)));
    double ECasCM = (xVal*xVal - mPi0*mPi0 + mCasM*mCasM)/(2*xVal);
    double pCasCMSq = ECasCM*ECasCM - mCasM*mCasM;
    double pCasCM = 0;
    if (pCasCMSq > 0) pCasCM =  sqrt(pCasCMSq);
    double GammaCasPi = gCasPi*pCasCM;

    //KLam threshold
    bool xValAboveKLam = true;
    if (xVal < (mK + mLam)) xValAboveKLam = false;
    double denomRealSQRT = centerFlatte*centerFlatte - xVal*xVal;
    double denomImSQRT  = centerFlatte*GammaCasPi + centerFlatte*GammaKLam;
    if (xValAboveKLam == false){
      double denomRealSQRT = centerFlatte*centerFlatte - xVal*xVal + centerFlatte*GammaKLam;
      double denomImSQRT  = centerFlatte*GammaCasPi;
    }
    double denomFlatte = denomRealSQRT*denomRealSQRT + denomImSQRT*denomImSQRT;
    double flatteDist = flatteMag*GammaCasPi/denomFlatte;

    if (kmLamType == true)flatteDist = flatteMag*GammaKLam/denomFlatte;
    if (kmLamType == true && xValAboveKLam == false) flatteDist = 0; 


    //////////
    flatteDist                        = Intensity_XiPi(xVal,flatteMag,centerFlatte,gCasPi,gKLam);
    if (kmLamType == true) flatteDist = Intensity_LamK(xVal,flatteMag,centerFlatte,gCasPi,gKLam);
    
    //////////
    double magF   = flatteMag;
    double centerF = centerFlatte;
    double gammaF  = gCasPi;

    //smearVal = 0;
    //centerF = 1.6;
    //magF = 5;
    //gammaF = 0.02;
    //cout<<"centerF = "<<centerF<<endl;
    //cout<<"flattePar0+3 = "<<flattePar0+3<<endl;
    //flatteDist = TMath::Voigt(xVal-centerF,smearVal,gammaF)*magF;
    //flatteDist = 20;
    //////////
    
    //sFun1 = 1;
    double bwTot = bw1 + bw2 + bw3 + bw4 + bw5 + bw6 + bw7 + bw8;
    //fitVal = bwTot + sFun1*bkg;
    fitVal = bwTot + bkgDalitz + flatteDist;

    int res9Par0 = flattePar0+4;
    double mag9   = fabs(parPtr[res9Par0]);
    double center9 = parPtr[res9Par0+1];
    double gamma9  = parPtr[res9Par0+2];
    double bw9 = TMath::Voigt(xVal-center9,smearVal,gamma9)*mag9;
    fitVal = fitVal + bw9;

    //if (xVal>2.011 && xVal<2.103) cout<<"res9Par0+1 = "<<res9Par0+1<<endl;
    
    return fitVal;
  }
};

//Numbering isobar,fit-histo



//SimFit
int nParj01 = 37;
int iparj01[37] = {0, //smear
  1, 
  2, //Shared (c1)
  3, //Shared (Gamma1)
  4, 
  5, //Shared (c2)
  6, //Shared (Gamma2)
  7,
  8, //Shared (c3)
  9, //Shared (Gamma3)
  10,
  11, //Shared (c4)
  12, //Shared (Gamma4)
  13,
  14, //Shared (c5)
  15, //Shared (Gamma5)
  16,
  17, //Shared (c6)
  18, //Shared (Gamma6)
  19,
  20, //Shared (c7)
  21, //Shared (Gamma7)
  22,
  23, //Shared (c8)
  24, //Shared (Gamma8)
  25,
  26,
  27,
  28,
  29,
  30,
  31, //Flatte Shared: gCasPi
  32, //Flatte Shared: gKLam
  33,  //Flatte Shared: center
  34,
  35, //Shared (c9)
  36  //Shared (Gamma9)
};
int parj02Base = 37;
int nParj02 = 37;
int iparj02[37] = { parj02Base,
  parj02Base+1, 
  2, //Shared (c1)
  3, //Shared (Gamma1)
  parj02Base+2, 
  5, //Shared (c2)
  6, //Shared (Gamma2)
  parj02Base+3, 
  8, //Shared (c3)
  9, //Shared (Gamma3)
  parj02Base+4, 
  11, //Shared (c4)
  12, //Shared (Gamma4)
  parj02Base+5, 
  14, //Shared (c5)
  15, //Shared (Gamma5)
  parj02Base+6, 
  17, //Shared (c6)
  18, //Shared (Gamma6)
  parj02Base+7, 
  20, //Shared (c7)
  21, //Shared (Gamma7)
  parj02Base+8, 
  23, //Shared (c8)
  24, //Shared (Gamma8)
  parj02Base+9, 
  parj02Base+10, 
  parj02Base+11, 
  parj02Base+12, 
  parj02Base+13,
  parj02Base+14,
  31, //Flatte Shared: gCasPi
  32, //Flatte Shared: gKLam
  33,  //Flatte Shared: center
  parj02Base+15,
  35, //Shared (c9)
  36  //Shared (Gamma9)
};
// Create the GlobalCHi2 structure                                                           \
                                                                                     
struct GlobalChi2 {
  GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,
               ROOT::Math::IMultiGenFunction & f2) :
    fChi2_1(&f1), fChi2_2(&f2) {}

 
  double operator() (const double *par) const {
    double p1[nParj01];
    for (int i = 0; i < nParj01; ++i) p1[i] = par[iparj01[i] ];

    double p2[nParj02];
    for (int i = 0; i < nParj02; ++i) p2[i] = par[iparj02[i] ];

    //double p3[nParj03];
    //for (int i = 0; i < nParj03; ++i) p3[i] = par[iparj03[i] ];

    //double p4[nParj04];
    //for (int i = 0; i < nParj04; ++i) p4[i] = par[iparj04[i] ];

    //return (*fChi2_1)(p1) + (*fChi2_2)(p2) + (*fChi2_3)(p3) + (*fChi2_4)(p4);
    //return (*fChi2_1)(p1) + (*fChi2_2)(p2) + (*fChi2_3)(p3);
    return (*fChi2_1)(p1) + (*fChi2_2)(p2);
  }

  const  ROOT::Math::IMultiGenFunction * fChi2_1;
  const  ROOT::Math::IMultiGenFunction * fChi2_2;
  //const  ROOT::Math::IMultiGenFunction * fChi2_3;
  //const  ROOT::Math::IMultiGenFunction * fChi2_4;
};



void simFitBranchesPlus(int bNumber = -1,int rangeType=-1){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  int kmLamNumber = 13;
  //int bNumber = -1;

  
  //fitType 0-> old 1620, 1-> new 1620
  
  //if (rfHybridType == 1) simPars[10] = 0.0092;
  //if (rfHybridType == 1) simPars[15] = 0.0124;
  //if (rfHybridType == 1) simPars[20] = 0.0111;

  //Setup the output file
  char  fId[120];
  char  hId[120];


  sprintf(fId,".//cm1StarSp17Sp18Fa182023BestCLStage3_final.root");
  TFile *inFileData = new TFile(fId);

  sprintf(fId,".//cm1StarKmLamBestCLSp17Sp18Fa182020Stage3_all.root");
  TFile *inFileData2 = new TFile(fId);

  sprintf(fId,".//cm1StarMCBestCLStage3pExpo6BKGdalitzRP1.root");
  TFile *inFileData3 = new TFile(fId);

  sprintf(fId,".//kmLamMCBestCLStage3BKGdalitz.root");
  TFile *inFileData4 = new TFile(fId);


  TH1D *hj01;
  TH1D *hj02;
  //TH1D *hj03;
  //TH1D *hj04;

  TH1D *hDalitz;
  TH1D *hDalitz2;

  //inFileData->GetObject("hMassXiStarCut2MCl6V4B0123",hj01);
  inFileData->GetObject( "hMassXiStarCut2MCl6V4B0123SBS",hj01);
  inFileData3->GetObject("hMassXiStarCut2MCl6V4B0123",hDalitz);

  cout<<"T:1"<<endl;
  
  sprintf(hId,"hMassXiStarCl4V4SBsubCut%d",kmLamNumber);inFileData2->GetObject(hId,hj02);
  sprintf(hId,"hMassXiStarCl4V4Cut%d",kmLamNumber); inFileData4->GetObject(hId,hDalitz2);

  //cout<<"kmLam hId = "<<hId<<endl;
  
  //inFileData2->GetObject("hMassXiStarCl4V4SBsubCut9",hj02);
  //inFileData2->GetObject("hMassXiStarCl4V4Cut2",hj02);

  //inFileData2->GetObject("hMassXiStarCl4V4SBsubCut9",hj02);

  //inFileData->GetObject("hMassXiStarCut1MCl6V4B2",hj03);      
  //inFileData->GetObject("hMassXiStarCut1MCl6V4B3",hj04);      


 if (bNumber >-1 && bNumber <4) {
   sprintf(hId,"hMassXiStarCut2MCl6V4B%d",bNumber);inFileData->GetObject(hId,hj01);
   //sprintf(hId,"hMassXiStarCut2MCl6V4B%dSBS",bNumber);inFileData->GetObject(hId,hj01);
    sprintf(hId,"hMassXiStarCut2MCl6V4B%d",bNumber);inFileData3->GetObject(hId,hDalitz);
    //sprintf(hId,"hMassXiStarCl4V4Cut%dB%d",kmLamNumber,bNumber);inFileData2->GetObject(hId,hj02);
    //sprintf(hId,"hMassXiStarCl4V4Cut%dB%d",kmLamNumber,bNumber); inFileData4->GetObject(hId,hDalitz2);
    sprintf(hId,"hMassXiStarCl4V4SBsubCut%dB%d",kmLamNumber,bNumber);inFileData2->GetObject(hId,hj02);
    sprintf(hId,"hMassXiStarCl4V4SBsubCut%dB%d",kmLamNumber,bNumber); inFileData4->GetObject(hId,hDalitz2);
  } else {
   sprintf(hId,"hMassXiStarCut2MCl6V4B0123");inFileData->GetObject(hId,hj01);
   //sprintf(hId,"hMassXiStarCut2MCl6V4B0123SBS");inFileData->GetObject(hId,hj01);
    sprintf(hId,"hMassXiStarCut2MCl6V4B0123");inFileData3->GetObject(hId,hDalitz);
    //sprintf(hId,"hMassXiStarCl4V4Cut%d",kmLamNumber);inFileData2->GetObject(hId,hj02);
    //sprintf(hId,"hMassXiStarCl4V4Cut%d",kmLamNumber); inFileData4->GetObject(hId,hDalitz2);
    sprintf(hId,"hMassXiStarCl4V4SBsubCut%d",kmLamNumber);inFileData2->GetObject(hId,hj02);
    sprintf(hId,"hMassXiStarCl4V4SBsubCut%d",kmLamNumber); inFileData4->GetObject(hId,hDalitz2);
  }


  
 TH1D *tagged_lumi;
 TH1D *tagged_lumi2;

  inFileData->GetObject("tagged_lumi",tagged_lumi);
  inFileData2->GetObject("tagged_lumi",tagged_lumi2);

  double e1 = 6.2;
  double e2 = 7.6;
  double e3 = 8.3;
  double e4 = 9.0;
  double e5 = 10.0;
  double e6 = 11.4;

  int lBinE0 = tagged_lumi->FindBin(e1+0.0001);
  int uBinE0 = tagged_lumi->FindBin(e2-0.0001);
  int lBinE1 = tagged_lumi->FindBin(e2+0.0001);
  int uBinE1 = tagged_lumi->FindBin(e3-0.0001);
  int lBinE2 = tagged_lumi->FindBin(e3+0.0001);
  int uBinE2 = tagged_lumi->FindBin(e4-0.0001);
  int lBinE3 = tagged_lumi->FindBin(e4+0.0001);
  int uBinE3 = tagged_lumi->FindBin(e5-0.0001);

  int lBinE0123 = tagged_lumi->FindBin(e1+0.0001);
  int uBinE0123 = tagged_lumi->FindBin(e5-0.0001);

  double lumiVal = tagged_lumi->Integral(lBinE0123,uBinE0123);
  if (bNumber == 0) lumiVal = tagged_lumi->Integral(lBinE0,uBinE0);
  if (bNumber == 1) lumiVal = tagged_lumi->Integral(lBinE1,uBinE1);
  if (bNumber == 2) lumiVal = tagged_lumi->Integral(lBinE2,uBinE2);
  if (bNumber == 3) lumiVal = tagged_lumi->Integral(lBinE3,uBinE3);

  //1530 efficiencies
  double rangeLow1530 = 1.51, rangeHigh1530 = 1.56;
  //rangeLow1530 = 1.97; rangeHigh1530 = 2.08;
  dPointsEffs_t effData1530;
  effData1530 = getEff1530(bNumber,rangeLow1530,rangeHigh1530);
  double eff1530      = effData1530.val;
  double eff1530Err   = effData1530.err;

  //1620 efficiencies
  //TEMP
  double rangeLow1620 = 1.52, rangeHigh1620 = 1.7;
  dPointsEffs_t effData1620;
  //effData1530 = getEff1530(bNumber,rangeLow1530,rangeHigh1530);
  double eff1620      = effData1530.val;
  double eff1620Err   = effData1530.err;

  //1690 efficiencies
  //TEMP
  double rangeLow1690 = 1.65, rangeHigh1690 = 1.75;
  dPointsEffs_t effData1690;
  //effData1530 = getEff1530(bNumber,rangeLow1530,rangeHigh1530);
  double eff1690      = effData1530.val;
  double eff1690Err   = effData1530.err;

  //1820 efficiencies
  double rangeLow1820 = 1.76, rangeHigh1820 = 1.89;
  dPointsEffs_t effData1820;
  effData1820 = getEff1820(bNumber,kmLamNumber,rangeLow1820,rangeHigh1820);
  double eff1820      = effData1820.val;
  double eff1820Err   = effData1820.err;
  double eff1820_2    = effData1820.val_2;
  double eff1820_2Err = effData1820.err_2;


  //2030 efficiencies
  double rangeLow2030 = 1.99, rangeHigh2030 = 2.07;
  dPointsEffs_t effData2030;
  effData2030 = getEff2030(bNumber,kmLamNumber,rangeLow2030,rangeHigh2030);
  double eff2030      = effData2030.val;
  double eff2030Err   = effData2030.err;
  double eff2030_2    = effData2030.val_2;
  double eff2030_2Err = effData2030.err_2;


  
  //2120 efficiencies
  //double rangeLow2120 = 2.08, rangeHigh2120 = 2.18;
  double rangeLow2120 = 1.6, rangeHigh2120 = 2.6;
  //TEMP
  double eff2120      = effData2030.val;
  double eff2120Err   = effData2030.err;
  double eff2120_2    = effData2030.val_2;
  double eff2120_2Err = effData2030.err_2;
  
  cout<<"####   1530   ####"<<endl;

  cout<<"eff1530   = "<<eff1530<<" +/- "<<eff1530Err<<endl;
  
  cout<<"####   2030   ####"<<endl;

  cout<<"eff2030   = "<<eff2030<<" +/- "<<eff2030Err<<endl;
  cout<<"eff2030_2 = "<<eff2030_2<<" +/- "<<eff2030_2Err<<endl;
  
  cout<<"##################"<<endl;
  //return;
  
  double rangeLow  = 1.45;
  double rangeLow2 = 1.6;
  double rangeHigh = 2.1;

  if (rangeType == 1) {
    //rangeLow  = 1.5;
    //rangeLow2  = 1.5;
    //rangeHigh = 1.73;
  }
  if (rangeType == 2) {
    rangeHigh = 2.6;
  }
  if (rangeType == 3) {
    //rangeLow = 2.25;
    rangeHigh = 2.2;
  }
  
  //rangeLow =  1.57;
  //rangeHigh = 1.67;
  
  bool reBin = false;
  if (reBin) {
    hj01->Rebin(2);
    hj02->Rebin(2);
    //hj03->Rebin(2);
    //hj04->Rebin(2);
  }
  //hDalitz->Rebin(2);

  int lineWidth = 1;
  hDalitz->Scale(0.25);
  hDalitz->SetLineColor(2);
  hDalitz->SetLineWidth(lineWidth);
  hDalitz->GetXaxis()->SetRangeUser(rangeLow,rangeHigh);
  hDalitz->Fit("pol3");

  hDalitz->GetFunction("pol3")->SetLineColor(kRed);

  double p0_Dalitz = hDalitz->GetFunction("pol3")->GetParameter(0);
  double p1_Dalitz = hDalitz->GetFunction("pol3")->GetParameter(1);
  double p2_Dalitz = hDalitz->GetFunction("pol3")->GetParameter(2);
  double p3_Dalitz = hDalitz->GetFunction("pol3")->GetParameter(3);

  //hDalitz2->Rebin(2);
  hDalitz2->Scale(0.25);
  hDalitz2->SetLineColor(2);
  hDalitz2->SetLineWidth(lineWidth);
  hDalitz2->GetXaxis()->SetRangeUser(rangeLow,rangeHigh);
  hDalitz2->Fit("pol3");

  hDalitz2->GetFunction("pol3")->SetLineColor(kRed);

  double p0_Dalitz2 = hDalitz2->GetFunction("pol3")->GetParameter(0);
  double p1_Dalitz2 = hDalitz2->GetFunction("pol3")->GetParameter(1);
  double p2_Dalitz2 = hDalitz2->GetFunction("pol3")->GetParameter(2);
  double p3_Dalitz2 = hDalitz2->GetFunction("pol3")->GetParameter(3);
  
  //TH1D *hj04; inFileData->GetObject("hMassXiStarCut2MCl6V4",hj04);

  //hj01->Draw();






  MyFitFun1 * fptr1 = new MyFitFun1();
  int nPar = 37;
  TF1 *fj01 = new TF1("fj01",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02 = new TF1("fj02",fptr1,&MyFitFun1::Evaluate,rangeLow2,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj03 = new TF1("fj03",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj04 = new TF1("fj04",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  fj01->SetLineWidth(lineWidth);
  fj02->SetLineWidth(lineWidth);
  
  //TF1 *fj01a = new TF1("fj01a",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj02a = new TF1("fj02a",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj03a = new TF1("fj03a",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj04a = new TF1("fj04a",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  //TF1 *fj01b = new TF1("fj01b",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj02b = new TF1("fj02b",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj03b = new TF1("fj03b",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj04b = new TF1("fj04b",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  //TF1 *fj01c = new TF1("fj01c",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj02c = new TF1("fj02c",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj03c = new TF1("fj03c",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  //TF1 *fj04c = new TF1("fj04c",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Flatte = new TF1("fj01Flatte",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Flatte = new TF1("fj02Flatte",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Bkg1530 = new TF1("fj01Bkg1530",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Bkg1820 = new TF1("fj01Bkg1820",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Bkg1820 = new TF1("fj02Bkg1820",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Bkg1620 = new TF1("fj01Bkg1620",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Bkg1620 = new TF1("fj02Bkg1620",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Bkg1690 = new TF1("fj01Bkg1690",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Bkg1690 = new TF1("fj02Bkg1690",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Bkg2030 = new TF1("fj01Bkg2030",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Bkg2030 = new TF1("fj02Bkg2030",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Bkg2120 = new TF1("fj01Bkg2120",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Bkg2120 = new TF1("fj02Bkg2120",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  
  TF1 *fj01Peak1 = new TF1("fj01Peak1",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");


  TF1 *fj01Peak2 = new TF1("fj01Peak2",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak2 = new TF1("fj02Peak2",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Peak3 = new TF1("fj01Peak3",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak3 = new TF1("fj02Peak3",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");

  TF1 *fj01Peak4 = new TF1("fj01Peak4",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak4 = new TF1("fj02Peak4",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Peak5 = new TF1("fj01Peak5",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak5 = new TF1("fj02Peak5",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Peak6 = new TF1("fj01Peak6",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak6 = new TF1("fj02Peak6",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Peak7 = new TF1("fj01Peak7",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak7 = new TF1("fj02Peak7",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Peak8 = new TF1("fj01Peak8",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak8 = new TF1("fj02Peak8",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01Peak9 = new TF1("fj01Peak9",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02Peak9 = new TF1("fj02Peak9",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  
  TF1 *fj01back = new TF1("fj01back",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj02back = new TF1("fj02back",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj03back = new TF1("fj03back",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");
  TF1 *fj04back = new TF1("fj04back",fptr1,&MyFitFun1::Evaluate,rangeLow,rangeHigh,nPar,"MyFitFun1","Evaluate");


  ROOT::Math::WrappedMultiTF1 wfj01(*fj01,1);
  ROOT::Math::WrappedMultiTF1 wfj02(*fj02,1);
  //ROOT::Math::WrappedMultiTF1 wfj03(*fj03,1);
  //ROOT::Math::WrappedMultiTF1 wfj04(*fj04,1);

  ROOT::Fit::DataOptions opt;

  ROOT::Fit::DataRange rangej01;
  rangej01.SetRange(rangeLow,rangeHigh);
  ROOT::Fit::BinData dataj01(opt,rangej01);
  ROOT::Fit::FillData(dataj01, hj01);

  ROOT::Fit::DataRange rangej02;
  rangej02.SetRange(rangeLow,rangeHigh);
  ROOT::Fit::BinData dataj02(opt,rangej02);
  ROOT::Fit::FillData(dataj02, hj02);

  //ROOT::Fit::DataRange rangej03;
  //rangej03.SetRange(rangeLow,rangeHigh);
  //ROOT::Fit::BinData dataj03(opt,rangej03);
  //ROOT::Fit::FillData(dataj03, hj03);

  //ROOT::Fit::DataRange rangej04;
  //rangej04.SetRange(rangeLow,rangeHigh);
  //ROOT::Fit::BinData dataj04(opt,rangej04);
  //ROOT::Fit::FillData(dataj04, hj04);

  ROOT::Fit::Chi2Function chi2_j01(dataj01, wfj01);
  ROOT::Fit::Chi2Function chi2_j02(dataj02, wfj02);
  //ROOT::Fit::Chi2Function chi2_j03(dataj03, wfj03);
  //ROOT::Fit::Chi2Function chi2_j04(dataj04, wfj04);

  //GlobalChi2 globalChi2(chi2_j01, chi2_j02, chi2_j03, chi2_j04);
  //GlobalChi2 globalChi2(chi2_j01, chi2_j02, chi2_j03);
  GlobalChi2 globalChi2(chi2_j01, chi2_j02);

 
  ROOT::Fit::Fitter fitter;


const int Npar = 53;//parj02Base+15+1+1

double sf1 = 0.01;
//Xi(1530)
double a1f1 = 1;
double c1f1 = 1.5352;
double w1f1 = 0.0099;
//Xi(1620)
double a2f1 = 0;
double c2f1 = 1.600;
double w2f1 = 0.032;
//Xi(1690)
double a3f1 = 1;
double c3f1 = 1.685;
double w3f1 = 0.010;
//Xi(1750)
double a4f1 = 1;
double c4f1 = 1.755;
double w4f1 = 0.01;
//Xi(1820)
double a5f1 = 1;
double c5f1 = 1.823;
double w5f1 = 0.024;
//Xi(1870) 
double a6f1 = 0;
double c6f1 = 1.870;
double w6f1 = 0.01;
//Xi(1950)
double a7f1 = 1;
double c7f1 = 1.920;
double w7f1 = 0.015;
//Xi(2030)
double a8f1 = 1;
double c8f1 = 2.025;
double w8f1 = 0.02;
//Xi(2120)
double a9f1 = 0;
double c9f1 = 2.120;
double w9f1 = 0.02;
//BKG1
//double b0f1 = 10;
//double b1f1 = 0;
//double b2f1 = 0;
//double b3f1 = 1.51;
//double b4f1 = 90.0;
 //BKG1-Dalitz
double b0f1 = p0_Dalitz;
double b1f1 = p1_Dalitz;
double b2f1 = p2_Dalitz;
double b3f1 = p3_Dalitz;
double b4f1 = 1.57;
//Flatte
 double flatte0f1 = 1;
 double flatte1f1 = 0.15;//gCasPi
 double flatte2f1 = 0.15;//gKmLam
 double flatte3f1 = 1.600;
 
//f2
 double sf2 = -0.0075;
 double a1f2 = 0;//1530 FIX
 double a2f2 = 0;//1620 FIX
 double a3f2 = 2;//1690
 double a4f2 = 0;//1750
 double a5f2 = 0;//1820
 double a6f2 = 0;//1870
 double a7f2 = 0;//1950
 double a8f2 = 0;//2030
 double a9f2 = 0;//2120
 double flatte0f2 = 1;
//BKG2
//double b0f2 = 0;
//double b1f2 = 0;//FIX
//double b2f2 = 0;//FIX
//double b3f2 = 1.23;
//double b4f2 = 200;
//BKG2
//double b0f2 = 0;
//double b1f2 = 0;
//double b2f2 = 0;
//double b3f2 = 0;
//double b4f2 = 1;
double b0f2 = p0_Dalitz2;
double b1f2 = p1_Dalitz2;
double b2f2 = p2_Dalitz2;
double b3f2 = p3_Dalitz2;
double b4f2 = 2.87;

 a1f1 = 0;a2f1 = 0;a3f1 = 0;a4f1 = 0;a5f1 = 0;a6f1 = 0;a7f1 = 0;a8f1 = 0;
 int flatteFit = 1;
 if (flatteFit == 1){
   a2f1 = 0;
   a2f2 = 0;
 }else{
   flatte0f1 = 0;
   flatte0f2 = 0;
 }

 //b4f1 = 1.38387;
 //b4f2 = 4.06721; 
 bool lockBKG = false;
 if(lockBKG==true){
   //b4f1 = 1.30474;
   //b4f2 = 3.58957;

   b4f1 = 1.38387;
   b4f2 = 4.06721; 
 }
 
 double simPars[Npar] = {sf1, //0
   a1f1,c1f1,w1f1, //1-3
   a2f1,c2f1,w2f1, //4-6
   a3f1,c3f1,w3f1, //7-9
   a4f1,c4f1,w4f1, //10-12
   a5f1,c5f1,w5f1, //13-15
   a6f1,c6f1,w6f1, //16-18
   a7f1,c7f1,w7f1, //19-21
   a8f1,c8f1,w8f1, //22-24
   b0f1,b1f1,b2f1,b3f1,b4f1,//25-29
   flatte0f1,flatte1f1,flatte2f1,flatte3f1, //30-33
   a9f1,c9f1,w9f1, //34-36
   sf2, //37
   a1f2,a2f2,a3f2,a4f2,a5f2,a6f2,a7f2,a8f2, //38-45
   b0f2,b1f2,b2f2,b3f2,b4f2, //46-50
   flatte0f2, //51
   a9f2 //52
 };

  
  // create before the parameter settings in order to fix or set range on them                                                          
  fitter.Config().SetParamsSettings(Npar,simPars);

  fitter.Config().ParSettings(0).SetLimits(0.005,0.01);//cout<<"T:0"<<endl;
  //1530
  fitter.Config().ParSettings(1).SetLimits(0,1000);//cout<<"T:1"<<endl;
  //fitter.Config().ParSettings(1).Fix();
  fitter.Config().ParSettings(2).SetLimits(c1f1-0.001,c1f1+0.001);//cout<<"T:2"<<endl;
  fitter.Config().ParSettings(3).SetLimits(w1f1-0.0019,w1f1+0.0017);//cout<<"T:3"<<endl;
  //1620
  if (flatteFit == 1){
    fitter.Config().ParSettings(4).Fix();
  }else{
    fitter.Config().ParSettings(4).SetLimits(0,1000);//cout<<"T:4"<<endl;
  }
  fitter.Config().ParSettings(5).SetLimits(c2f1-0.001,c2f1+0.001);//cout<<"T:5"<<endl;
  fitter.Config().ParSettings(6).SetLimits(w2f1-0.009,w2f1+0.008);//cout<<"T:6"<<endl;
  //1690
  fitter.Config().ParSettings(7).SetLimits(0,1000);//cout<<"T:7"<<endl;
  //fitter.Config().ParSettings(7).Fix();
  fitter.Config().ParSettings(8).SetLimits(c3f1-0.05,c3f1+0.020);//cout<<"T:8"<<endl;
  fitter.Config().ParSettings(9).SetLimits(w3f1-0.015,w3f1+0.015);//cout<<"T:9"<<endl;
  //1750
  //fitter.Config().ParSettings(10).Fix();
  fitter.Config().ParSettings(10).SetLimits(0,1000);//cout<<"T:10"<<endl;
  fitter.Config().ParSettings(11).SetLimits(c4f1-0.005,c4f1+0.005);//cout<<"T:11"<<endl;
  fitter.Config().ParSettings(12).SetLimits(w4f1-0.010,w4f1+0.005);//cout<<"T:12"<<endl;
  //1820
  //fitter.Config().ParSettings(13).Fix();
  fitter.Config().ParSettings(13).SetLimits(0,1000);//cout<<"T:13"<<endl;
  fitter.Config().ParSettings(14).SetLimits(c5f1-0.005,c5f1+0.005);//cout<<"T:14"<<endl;
  //fitter.Config().ParSettings(15).SetLimits(w5f1-0.010,w5f1+0.015);//cout<<"T:15"<<endl;
  fitter.Config().ParSettings(15).SetLimits(w5f1-0.005,w5f1+0.005);//cout<<"T:15"<<endl;
  //1870 
  //fitter.Config().ParSettings(16).Fix();
  fitter.Config().ParSettings(16).SetLimits(0,1000);//cout<<"T:16"<<endl;
  fitter.Config().ParSettings(17).SetLimits(c6f1-0.01,c6f1+0.01);//cout<<"T:17"<<endl;
  fitter.Config().ParSettings(18).SetLimits(w6f1-0.010,w6f1+0.005);//cout<<"T:18"<<endl;
  //1950
  //fitter.Config().ParSettings(19).Fix();
  fitter.Config().ParSettings(19).SetLimits(0,1000);//cout<<"T:19"<<endl;
  fitter.Config().ParSettings(20).SetLimits(c7f1-0.010,c7f1+0.010);//cout<<"T:20"<<endl;
  fitter.Config().ParSettings(21).SetLimits(w7f1-0.005,w7f1+0.005);//cout<<"T:21"<<endl;
  //2030
  //fitter.Config().ParSettings(22).Fix();
  fitter.Config().ParSettings(22).SetLimits(0,1000);//cout<<"T:22"<<endl;
  fitter.Config().ParSettings(23).SetLimits(c8f1-0.005,c8f1+0.005);//cout<<"T:23"<<endl;
  fitter.Config().ParSettings(24).SetLimits(w8f1-0.005,w8f1+0.007);//cout<<"T:24"<<endl;
  //BKG
  //fitter.Config().ParSettings(25).SetLimits(0,1000);
  //fitter.Config().ParSettings(28).SetLimits(1.48,1.52);
  //fitter.Config().ParSettings(29).SetLimits(80.0,150.0);
  //fitter.Config().ParSettings(28).SetLimits(1.28,1.30);
  //fitter.Config().ParSettings(29).SetLimits(80.0,150.0);
  if (lockBKG==true) fitter.Config().ParSettings(29).Fix();   
  //BKG-Dalitz
  fitter.Config().ParSettings(25).Fix();//cout<<"T:25"<<endl;
  fitter.Config().ParSettings(26).Fix();//cout<<"T:26"<<endl;
  fitter.Config().ParSettings(27).Fix();//cout<<"T:27"<<endl;
  fitter.Config().ParSettings(28).Fix();//cout<<"T:28"<<endl;
  //fitter.Config().ParSettings(29).Fix();
  //fitter.Config().ParSettings(29).SetLimits(1.56,1.58);
  //Flatte
  if (flatteFit == 1){
    fitter.Config().ParSettings(30).SetLimits(0.0,100); //cout<<"T:29"<<endl;  //MAG 
  }else{
    fitter.Config().ParSettings(30).Fix(); //MAG
  }

  //Flatte
  //fitter.Config().ParSettings(31).Fix();//gCasPi
  fitter.Config().ParSettings(31).SetLimits(0.1,0.4);//cout<<"T:30"<<endl;  //gCasPi
  fitter.Config().ParSettings(32).SetLimits(0.01,0.2);//cout<<"T:31"<<endl;  //gKmLam
  //fitter.Config().ParSettings(32).Fix();//gKmLam
  fitter.Config().ParSettings(33).SetLimits(1.6,1.640);//cout<<"T:32"<<endl;//center

  //2120
  //fitter.Config().ParSettings(34).Fix();
  bool fit2120 = false;
  if (rangeType == 2) fit2120 = true;
  if (fit2120 == true) {
    fitter.Config().ParSettings(34).SetLimits(0,1000);//cout<<"T:33"<<endl;
  }else{
    fitter.Config().ParSettings(34).Fix();
  }
  //fitter.Config().ParSettings(35).SetLimits(c9f1-0.015,c9f1+0.015);//cout<<"T:34"<<endl;
  //fitter.Config().ParSettings(36).SetLimits(0.013,w9f1+0.012);//cout<<"T:35"<<endl;

  fitter.Config().ParSettings(35).SetLimits(2.1,2.15);//cout<<"T:34"<<endl;
  //fitter.Config().ParSettings(36).SetLimits(0.005,0.35);//cout<<"T:35"<<endl;

  
  //F2
  fitter.Config().ParSettings(37).SetLimits(-0.010,-0.005);//SMEAR FOR F2
  fitter.Config().ParSettings(38).Fix();//Xi(1530) FOR F2
  //fitter.Config().ParSettings(39).Fix();
  //fitter.Config().ParSettings(40).Fix();//Xi(1690) FOR F2 
  //fitter.Config().ParSettings(41).Fix();//Xi(1750) FOR F2
  //fitter.Config().ParSettings(42).Fix();//Xi(1820) FOR F2 
  //fitter.Config().ParSettings(43).Fix();//Xi(1870) FOR F2
  //fitter.Config().ParSettings(44).Fix();//Xi(1950) FOR F2
  //fitter.Config().ParSettings(45).Fix();//Xi(2030) FOR F2 

  //Xi(1620) FOR F2
  if (flatteFit == 1){
    fitter.Config().ParSettings(39).Fix();
  }else{
    fitter.Config().ParSettings(39).SetLimits(0,1000);
  }
  
  //fitter.Config().ParSettings(39).SetLimits(0,1000);//Xi(1620) FOR F2
  fitter.Config().ParSettings(40).SetLimits(0,1000);//Xi(1690) FOR F2
  fitter.Config().ParSettings(41).SetLimits(0,1000);//Xi(1750) FOR F2
  fitter.Config().ParSettings(42).SetLimits(0,1000);//Xi(1820) FOR F2
  fitter.Config().ParSettings(43).SetLimits(0,1000);//Xi(1950) FOR F2 
  fitter.Config().ParSettings(44).SetLimits(0,1000);//Xi(1970) FOR F2
  fitter.Config().ParSettings(45).SetLimits(0,1000);//Xi(2030) FOR F2
  if (fit2120 == true) {
    fitter.Config().ParSettings(52).SetLimits(0,1000);//Xi(2120) FOR F2
  }else{
    fitter.Config().ParSettings(52).Fix();//Xi(2120) FOR F2
  }


  //fitter.Config().ParSettings(46).SetLimits(0,1000);//b0f2
  fitter.Config().ParSettings(46).Fix();//b0f2
  fitter.Config().ParSettings(47).Fix();//b1f2
  fitter.Config().ParSettings(48).Fix();//b2f2
  fitter.Config().ParSettings(49).Fix();//b3f2
  if (lockBKG==true) fitter.Config().ParSettings(50).Fix();   
  //fitter.Config().ParSettings(50).SetLimits(2.8,2.94);

  if (flatteFit == 1){//MAG
    fitter.Config().ParSettings(51).SetLimits(0.0,200);
  }else{
    fitter.Config().ParSettings(51).Fix(); 
  }

  //fitter.Config().ParSettings(51).Fix();
  
  //fitter.Config().ParSettings(49).SetLimits(1.24,1.28);
  //fitter.Config().ParSettings(50).SetLimits(100.0,250.0);
  
  //fitter.Config().ParSettings(0).SetLimits(0,100000);
  //fitter.Config().ParSettings(0).SetStepSize(1);
  //fitter.Config().ParSettings(1).SetLimits(1.535-0.01,1.535+0.01);
  //fitter.Config().ParSettings(1).SetStepSize(0.0001);// 0.1 MeV
  //fitter.Config().ParSettings(2).SetLimits(0.0099-0.0099,0.0099+0.0017+0.01);
  //fitter.Config().ParSettings(2).SetStepSize(0.0001);// 0.1 MeV

  //fitter.Config().ParSettings(3).SetLimits(0,100000);
  //fitter.Config().ParSettings(3).SetStepSize(1);
  //fitter.Config().ParSettings(4).SetLimits(1.56,1.64);
  //fitter.Config().ParSettings(4).SetStepSize(0.0001);// 0.1 MeV
  //fitter.Config().ParSettings(5).SetLimits(0.00001,0.05);
  //fitter.Config().ParSettings(5).SetStepSize(0.0001);// 0.1 MeV

  //fitter.Config().ParSettings(6).SetLimits(0.0,100000);
  //fitter.Config().ParSettings(6).SetStepSize(0.0001);
  //fitter.Config().ParSettings(7).SetLimits(1.65,1.7+0.001);

  //fitter.Config().ParSettings(7).SetLimits(1.65,1.7);
  //fitter.Config().ParSettings(7).SetStepSize(0.0001);// 0.1 MeV
  //fitter.Config().ParSettings(8).SetLimits(0.00001,0.0106+0.005);
  //fitter.Config().ParSettings(8).SetStepSize(0.0001);// 0.1 MeV

  //fitter.Config().ParSettings(9).Fix(); //p0
  //fitter.Config().ParSettings(10).Fix();//p1

  //fitter.Config().ParSettings(11).SetLimits(0,100000);
  //fitter.Config().ParSettings(11).SetStepSize(1);
  //fitter.Config().ParSettings(12).SetLimits(0,100000);
  //fitter.Config().ParSettings(12).SetStepSize(1);
  //fitter.Config().ParSettings(13).SetLimits(0,100000);
  //fitter.Config().ParSettings(13).SetStepSize(1);
  //fitter.Config().ParSettings(14).Fix();//p0
  //fitter.Config().ParSettings(15).Fix();//p1

  //fitter.Config().ParSettings(16).SetLimits(0,100000);
  //fitter.Config().ParSettings(16).SetStepSize(1);
  //fitter.Config().ParSettings(17).SetLimits(0,100000);
  //fitter.Config().ParSettings(17).SetStepSize(1);
  //fitter.Config().ParSettings(18).SetLimits(0,100000);
  //fitter.Config().ParSettings(18).SetStepSize(1);
  //fitter.Config().ParSettings(19).Fix();//p0
  //fitter.Config().ParSettings(20).Fix();//p1

  //fitter.Config().ParSettings(21).SetLimits(0,100000);
  //fitter.Config().ParSettings(21).SetStepSize(1);
  //fitter.Config().ParSettings(22).SetLimits(0,100000);
  //fitter.Config().ParSettings(22).SetStepSize(1);
  //fitter.Config().ParSettings(23).SetLimits(0,100000);
  //fitter.Config().ParSettings(23).SetStepSize(1);
  //fitter.Config().ParSettings(24).Fix();//p0
  //fitter.Config().ParSettings(25).Fix();//p1

  fitter.Config().MinimizerOptions().SetPrintLevel(0);
  fitter.Config().SetMinimizer("Minuit","Migrad");
  //fitter.FitFCN(Npar,globalChi2,0,dataj01.Size()+dataj02.Size()+dataj03.Size()+dataj04.Size(),true);
  //fitter.FitFCN(Npar,globalChi2,0,dataj01.Size()+dataj02.Size()+dataj03.Size(),true);
  fitter.FitFCN(Npar,globalChi2,0,dataj01.Size()+dataj02.Size(),true);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  fj01->SetFitResult( result, iparj01);
  fj02->SetFitResult( result, iparj02);
  //fj03->SetFitResult( result, iparj03);
  //fj04->SetFitResult( result, iparj04);

  //fj01a->SetFitResult( result, iparj01);
  //fj02a->SetFitResult( result, iparj02);
  //fj03a->SetFitResult( result, iparj03);
  //fj04a->SetFitResult( result, iparj04);

  //fj01b->SetFitResult( result, iparj01);
  //fj02b->SetFitResult( result, iparj02);
  //fj03b->SetFitResult( result, iparj03);
  //fj04b->SetFitResult( result, iparj04);

  //fj01c->SetFitResult( result, iparj01);
  //fj02c->SetFitResult( result, iparj02);
  //fj03c->SetFitResult( result, iparj03);
  //fj04c->SetFitResult( result, iparj04);

  fj01back->SetFitResult( result, iparj01);
  fj02back->SetFitResult( result, iparj02);

  fj01Flatte->SetFitResult( result, iparj01);
  fj02Flatte->SetFitResult( result, iparj02);

  fj01Bkg1530->SetFitResult( result, iparj01);
  fj01Bkg1530->SetLineWidth(lineWidth);
  
  fj01Bkg1620->SetFitResult( result, iparj01);
  fj02Bkg1620->SetFitResult( result, iparj02);

  fj01Bkg1690->SetFitResult( result, iparj01);
  fj02Bkg1690->SetFitResult( result, iparj02);

  fj01Bkg1820->SetFitResult( result, iparj01);
  fj02Bkg1820->SetFitResult( result, iparj02);

  fj01Bkg2030->SetFitResult( result, iparj01);
  fj02Bkg2030->SetFitResult( result, iparj02);

  fj01Bkg2120->SetFitResult( result, iparj01);
  fj02Bkg2120->SetFitResult( result, iparj02);

  
  fj01Peak1->SetFitResult( result, iparj01);


  fj01Peak2->SetFitResult( result, iparj01);
  fj02Peak2->SetFitResult( result, iparj02);
  
  fj01Peak3->SetFitResult( result, iparj01);
  fj02Peak3->SetFitResult( result, iparj02);
  
  fj01Peak4->SetFitResult( result, iparj01);
  fj02Peak4->SetFitResult( result, iparj02);
  
  fj01Peak5->SetFitResult( result, iparj01);
  fj02Peak5->SetFitResult( result, iparj02);
  
  fj01Peak6->SetFitResult( result, iparj01);
  fj02Peak6->SetFitResult( result, iparj02);
  
  fj01Peak7->SetFitResult( result, iparj01);
  fj02Peak7->SetFitResult( result, iparj02);
  
  fj01Peak8->SetFitResult( result, iparj01);
  fj02Peak8->SetFitResult( result, iparj02);
  
  fj01Peak9->SetFitResult( result, iparj01);
  fj02Peak9->SetFitResult( result, iparj02);
  
  //fj03back->SetFitResult( result, iparj03);
  //fj04back->SetFitResult( result, iparj04);

  //fj01->SetParameter(9,0);
  //fj01->SetParameter(10,0);

  //RED
  //fj01a->SetParameter(3,0);fj01a->SetParameter(6,0);fj01a->SetParameter(9,0);
  //fj02a->SetParameter(3,0);fj02a->SetParameter(6,0);fj02a->SetParameter(9,0);
  //fj03a->SetParameter(3,0);fj03a->SetParameter(6,0);fj03a->SetParameter(9,0);
  //fj04a->SetParameter(3,0);fj04a->SetParameter(6,0);fj04a->SetParameter(9,0);

  //GREEN
  //fj01b->SetParameter(0,0);fj01b->SetParameter(6,0);fj01b->SetParameter(9,0);
  //fj02b->SetParameter(0,0);fj02b->SetParameter(6,0);fj02b->SetParameter(9,0);
  //fj03b->SetParameter(0,0);fj03b->SetParameter(6,0);fj03b->SetParameter(9,0);
  //fj04b->SetParameter(0,0);fj04b->SetParameter(6,0);fj04b->SetParameter(9,0);

  //BLUE
  //fj01c->SetParameter(0,0);fj01c->SetParameter(3,0);fj01c->SetParameter(9,0);
  //fj02c->SetParameter(0,0);fj02c->SetParameter(3,0);fj02c->SetParameter(9,0);
  //fj03c->SetParameter(0,0);fj03c->SetParameter(3,0);fj03c->SetParameter(9,0);
  ///fj04c->SetParameter(0,0);fj04c->SetParameter(3,0);fj04c->SetParameter(9,0);

  //Background
  fj01back->SetParameter(1,0);
  fj01back->SetParameter(4,0);
  fj01back->SetParameter(7,0);
  fj01back->SetParameter(10,0);
  fj01back->SetParameter(13,0);
  fj01back->SetParameter(16,0);
  fj01back->SetParameter(19,0);
  fj01back->SetParameter(22,0);
  fj01back->SetParameter(30,0);
  fj01back->SetParameter(34,0);

  fj02back->SetParameter(1,0);
  fj02back->SetParameter(4,0);
  fj02back->SetParameter(7,0);
  fj02back->SetParameter(10,0);
  fj02back->SetParameter(13,0);
  fj02back->SetParameter(16,0);
  fj02back->SetParameter(19,0);
  fj02back->SetParameter(22,0);
  fj02back->SetParameter(30,0);
  fj02back->SetParameter(34,0);

  //Flatte
  fj01Flatte->SetParameter(1,0);
  fj01Flatte->SetParameter(4,0);
  fj01Flatte->SetParameter(7,0);
  fj01Flatte->SetParameter(10,0);
  fj01Flatte->SetParameter(13,0);
  fj01Flatte->SetParameter(16,0);
  fj01Flatte->SetParameter(19,0);
  fj01Flatte->SetParameter(22,0);
  fj01Flatte->SetParameter(29,0);
  fj01Flatte->SetParameter(34,0);

  fj02Flatte->SetParameter(1,0);
  fj02Flatte->SetParameter(4,0);
  fj02Flatte->SetParameter(7,0);
  fj02Flatte->SetParameter(10,0);
  fj02Flatte->SetParameter(13,0);
  fj02Flatte->SetParameter(16,0);
  fj02Flatte->SetParameter(19,0);
  fj02Flatte->SetParameter(22,0);
  fj02Flatte->SetParameter(29,0);
  fj02Flatte->SetParameter(34,0);

  //Bkg1530
  fj01Bkg1530->SetParameter(1,0);
  
  //Bkg1620
  fj01Bkg1620->SetParameter(4,0);
  fj02Bkg1620->SetParameter(4,0);

  //Bkg1690
  fj01Bkg1690->SetParameter(4,0);
  fj02Bkg1690->SetParameter(4,0);

  //Bkg1820
  fj01Bkg1820->SetParameter(13,0);
  fj02Bkg1820->SetParameter(13,0);

  //Bkg2030
  fj01Bkg2030->SetParameter(22,0);
  fj02Bkg2030->SetParameter(22,0);

  //Bkg2120
  fj01Bkg2120->SetParameter(34,0);
  fj02Bkg2120->SetParameter(34,0);

  //Peak1
  //fj01Peak1->SetParameter(1,0);
  fj01Peak1->SetParameter(4,0);
  fj01Peak1->SetParameter(7,0);
  fj01Peak1->SetParameter(10,0);
  fj01Peak1->SetParameter(13,0);
  fj01Peak1->SetParameter(16,0);
  fj01Peak1->SetParameter(19,0);
  fj01Peak1->SetParameter(22,0);
  fj01Peak1->SetParameter(29,0);
  fj01Peak1->SetParameter(30,0);
  fj01Peak1->SetParameter(34,0);

  //Peak3
  fj01Peak3->SetParameter(1,0);
  fj01Peak3->SetParameter(4,0);
  //fj01Peak3->SetParameter(7,0);
  fj01Peak3->SetParameter(10,0);
  fj01Peak3->SetParameter(13,0);
  fj01Peak3->SetParameter(16,0);
  fj01Peak3->SetParameter(19,0);
  fj01Peak3->SetParameter(22,0);
  fj01Peak3->SetParameter(29,0);
  fj01Peak3->SetParameter(30,0);
  fj01Peak3->SetParameter(34,0);

  fj02Peak3->SetParameter(1,0);
  fj02Peak3->SetParameter(4,0);
  //fj02Peak3->SetParameter(7,0);
  fj02Peak3->SetParameter(10,0);
  fj02Peak3->SetParameter(13,0);
  fj02Peak3->SetParameter(16,0);
  fj02Peak3->SetParameter(19,0);
  fj02Peak3->SetParameter(22,0);
  fj02Peak3->SetParameter(29,0);
  fj02Peak3->SetParameter(30,0);
  fj02Peak3->SetParameter(34,0);

  //Peak5
  fj01Peak5->SetParameter(1,0);
  fj01Peak5->SetParameter(4,0);
  fj01Peak5->SetParameter(7,0);
  fj01Peak5->SetParameter(10,0);
  //fj01Peak5->SetParameter(13,0);
  fj01Peak5->SetParameter(16,0);
  fj01Peak5->SetParameter(19,0);
  fj01Peak5->SetParameter(22,0);
  fj01Peak5->SetParameter(29,0);
  fj01Peak5->SetParameter(30,0);
  fj01Peak5->SetParameter(34,0);

  fj02Peak5->SetParameter(1,0);
  fj02Peak5->SetParameter(4,0);
  fj02Peak5->SetParameter(7,0);
  fj02Peak5->SetParameter(10,0);
  //fj02Peak5->SetParameter(13,0);
  fj02Peak5->SetParameter(16,0);
  fj02Peak5->SetParameter(19,0);
  fj02Peak5->SetParameter(22,0);
  fj02Peak5->SetParameter(29,0);
  fj02Peak5->SetParameter(30,0);
  fj02Peak5->SetParameter(34,0);
  
  //Peak8
  fj01Peak8->SetParameter(1,0);
  fj01Peak8->SetParameter(4,0);
  fj01Peak8->SetParameter(7,0);
  fj01Peak8->SetParameter(10,0);
  fj01Peak8->SetParameter(13,0);
  fj01Peak8->SetParameter(16,0);
  fj01Peak8->SetParameter(19,0);
  //fj01Peak8->SetParameter(22,0);
  fj01Peak8->SetParameter(29,0);
  fj01Peak8->SetParameter(30,0);
  fj01Peak8->SetParameter(34,0);

  fj02Peak8->SetParameter(1,0);
  fj02Peak8->SetParameter(4,0);
  fj02Peak8->SetParameter(7,0);
  fj02Peak8->SetParameter(10,0);
  fj02Peak8->SetParameter(13,0);
  fj02Peak8->SetParameter(16,0);
  fj02Peak8->SetParameter(19,0);
  //fj02Peak8->SetParameter(22,0);
  fj02Peak8->SetParameter(29,0);
  fj02Peak8->SetParameter(30,0);
  fj02Peak8->SetParameter(34,0);

  //Peak9
  fj01Peak9->SetParameter(1,0);
  fj01Peak9->SetParameter(4,0);
  fj01Peak9->SetParameter(7,0);
  fj01Peak9->SetParameter(10,0);
  fj01Peak9->SetParameter(13,0);
  fj01Peak9->SetParameter(16,0);
  fj01Peak9->SetParameter(19,0);
  fj01Peak9->SetParameter(22,0);
  fj01Peak9->SetParameter(29,0);
  fj01Peak9->SetParameter(30,0);
  //fj01Peak9->SetParameter(34,0);

  fj02Peak9->SetParameter(1,0);
  fj02Peak9->SetParameter(4,0);
  fj02Peak9->SetParameter(7,0);
  fj02Peak9->SetParameter(10,0);
  fj02Peak9->SetParameter(13,0);
  fj02Peak9->SetParameter(16,0);
  fj02Peak9->SetParameter(19,0);
  fj02Peak9->SetParameter(22,0);
  fj02Peak9->SetParameter(29,0);
  fj02Peak9->SetParameter(30,0);
  //fj02Peak9->SetParameter(34,0);
  

  //fj02back->SetParameter(0,0);
  //fj03back->SetParameter(0,0);
  //fj04back->SetParameter(0,0);

  //fj01a->SetLineColor(2);
  //fj02a->SetLineColor(2);
  //fj03a->SetLineColor(2);
  //fj04a->SetLineColor(2);

  //fj01b->SetLineColor(3);
  //fj02b->SetLineColor(3);
  //fj03b->SetLineColor(3);
  //fj04b->SetLineColor(3);

  //fj01c->SetLineColor(4);
  //fj02c->SetLineColor(4);
  //fj03c->SetLineColor(4);
  //fj04c->SetLineColor(4);

  //hj01->SetMinimum(0);
  hj01->SetLineColor(1);
  hj01->SetMarkerStyle(21);
  hj01->SetMarkerColor(1);
  hj01->GetXaxis()->SetTitle("Mass(#Xi^{-}#pi^{0})/GeV");
  hj01->GetYaxis()->SetTitle("Yield");
  hj01->GetXaxis()->SetTitleSize(0.06);
  hj01->GetYaxis()->SetTitleSize(0.06);
  hj01->GetXaxis()->SetTitleOffset(0.8);
  hj01->GetYaxis()->SetTitleOffset(0.4);

  //hj02->SetMinimum(0);
  hj02->SetLineColor(1);
  hj02->SetMarkerStyle(21);
  hj02->SetMarkerColor(1);
  hj02->GetXaxis()->SetTitle("Mass(#Lambda#pi^{-})/GeV");
  hj02->GetYaxis()->SetTitle("Yield");
  hj02->GetXaxis()->SetTitleSize(0.06);
  hj02->GetYaxis()->SetTitleSize(0.06);
  hj02->GetXaxis()->SetTitleOffset(0.8);
  hj02->GetYaxis()->SetTitleOffset(0.4);

  //hj03->SetMinimum(0);
  //hj03->SetLineColor(1);
  //hj03->SetMarkerStyle(21);
  //hj03->SetMarkerColor(1);
  //hj03->GetXaxis()->SetTitle("Mass(#Xi^{-} #pi^{0})/GeV");
  //hj03->GetYaxis()->SetTitle("Yield");
  //hj03->GetXaxis()->SetTitleSize(0.06);
  //hj03->GetYaxis()->SetTitleSize(0.06);
  //hj03->GetXaxis()->SetTitleOffset(0.8);
  //hj03->GetYaxis()->SetTitleOffset(0.4);

  ////hj04->SetMinimum(0);
  //hj04->SetLineColor(1);
  //hj04->SetMarkerStyle(21);
  //hj04->SetMarkerColor(1);
  //hj04->GetXaxis()->SetTitle("Mass(#Xi^{-} #pi^{0})/GeV");
  //hj04->GetYaxis()->SetTitle("Yield");
  //hj04->GetXaxis()->SetTitleSize(0.05);
  //hj04->GetYaxis()->SetTitleSize(0.05);
  //hj04->GetYaxis()->SetTitleOffset(0.4);


  int iShow = 1;

  TCanvas * c1 = new TCanvas("Simfit","AllSpFull",
                             10,10,800,600);
  c1->Divide(1,2);
  c1->cd(1);
  //gStyle->SetOptFit(1111);
  if (rangeType == 1) {
    rangeHigh = 1.73;
  }
  hj01->GetXaxis()->SetRangeUser(rangeLow,rangeHigh);
  if (rangeType !=1) hj01->SetMinimum(-10);
  hj01->SetLineWidth(lineWidth);
  double maxHj01 = hj01->GetMaximum();
  hj01->SetMaximum(1.25*maxHj01);


  hj01->Draw("e1");


  fj01->SetLineWidth(lineWidth);
  fj01->Draw("same");

  double scaleFacDalitz = fj01back->GetParameter(29);
  hDalitz->GetFunction("pol3")->SetBit(TF1::kNotDraw);
  hDalitz->Scale(scaleFacDalitz);
  //hDalitz->Draw("same");

  if (iShow == 0){
    //fj01a->Draw("same");
    //fj01b->Draw("same");
    //fj01c->Draw("same");
  }


  if (iShow == 0) {
    //legend->AddEntry(fj01a,"#Xi(1530)","l");
    //legend->AddEntry(fj01b,"#Xi(1620)","l");
    //legend->AddEntry(fj01c,"#Xi(1690)","l");
  }
  if (iShow == 1) {
    //legend->AddEntry(fj01back,"Background","l");
    //legend->AddEntry(fj01,"Signal","l");

  }
  //legend->Draw();

  c1->cd(2);
  if (rangeType !=1) hj02->SetMinimum(-10);
  //gStyle->SetOptFit(1111);
  hj02->GetXaxis()->SetRangeUser(rangeLow,rangeHigh);
  hj02->Draw("e1");

  fj02->GetXaxis()->SetRangeUser(1.6,3.0);
  fj02->Draw("same");

  double scaleFacDalitz2 = fj02back->GetParameter(29);

  hDalitz2->GetFunction("pol3")->SetBit(TF1::kNotDraw);
  hDalitz2->Scale(scaleFacDalitz2);
  //hDalitz2->Draw("same");
  cout<<"scaleFacDalitz2 = "<<scaleFacDalitz2<<endl;

  if (iShow == 0){
    //fj02a->Draw("same");
    //fj02b->Draw("same");
    //fj02c->Draw("same");
  }


  double C_Xi_1530_val = fj01->GetParameter(1);
  double C_Xi_1530_err = fj01->GetParError(1);
  double W_Xi_1530_val = fj01->GetParameter(2);
  double W_Xi_1530_err = fj01->GetParError(2);
  cout<<""<<endl;
  cout<<"-------------------"<<endl;
  cout<<"PDG: Xi(1530) Center = 1535.0     +/- 0.6"<<endl;
  cout<<"Fit: Xi(1530) Center = "<<1000*C_Xi_1530_val<<" +/- "<<1000*C_Xi_1530_err<<endl;
  cout<<""<<endl;
  cout<<"PDG: Xi(1530) Width = 9.9 +1.7/-1.9"<<endl;
  cout<<"Fit: Xi(1530) Width = "<<1000*W_Xi_1530_val<<" +/- "<<1000*W_Xi_1530_err<<endl;

  double C_Xi_1620_val = fj01->GetParameter(4);
  double C_Xi_1620_err = fj01->GetParError(4);
  double W_Xi_1620_val = fj01->GetParameter(5);
  double W_Xi_1620_err = fj01->GetParError(5);
  cout<<""<<endl;
  cout<<"-------------------"<<endl;
  cout<<"PDG: Xi(1620) Center = ~ 1620   +/- ?"<<endl;
  cout<<"Fit: Xi(1620) Center = "<<1000*C_Xi_1620_val<<" +/- "<<1000*C_Xi_1620_err<<endl;
  cout<<""<<endl;
  cout<<"PDG: Xi(1620) Width = 32 +8/-9 "<<endl;
  cout<<"Fit: Xi(1620) Width = "<<1000*W_Xi_1620_val<<" +/- "<<1000*W_Xi_1620_err<<endl;

  double C_Xi_1690_val = fj01->GetParameter(7);
  double C_Xi_1690_err = fj01->GetParError(7);
  double W_Xi_1690_val = fj01->GetParameter(8);
  double W_Xi_1690_err = fj01->GetParError(8);
  cout<<""<<endl;
  cout<<"-------------------"<<endl;
  cout<<"PDG: Xi(1690) Center = 1690     +/- 10"<<endl;
  cout<<"Fit: Xi(1690) Center = "<<1000*C_Xi_1690_val<<" +/- "<<1000*C_Xi_1690_err<<endl;
  cout<<""<<endl;
  cout<<"PDG: Xi(1690) Width = 20 +/- 15"<<endl;
  cout<<"Fit: Xi(1690) Width = "<<1000*W_Xi_1690_val<<" +/- "<<1000*W_Xi_1690_err<<endl;

  double smearVal = 0.00926;
  
  //Xi 1530
  double sigma_Xi_1530 = sqrt(pow(W_Xi_1530_val/2.35,2)+pow(smearVal,2));
  double low_Xi_1530 = C_Xi_1530_val - 3*sigma_Xi_1530;
  double high_Xi_1530 = C_Xi_1530_val + 3*sigma_Xi_1530;
  int low_Xi_1530_Bin = hj01->GetXaxis()->FindBin(low_Xi_1530) + 1;
  int high_Xi_1530_Bin = hj01->GetXaxis()->FindBin(high_Xi_1530) - 1;

  //TH1D * hYBin0 = (TH1D*) hj01->Clone("hYBin0");
  //TH1D * hYBin1 = (TH1D*) hj02->Clone("hYBin1");
  //TH1D * hYBin2 = (TH1D*) hj03->Clone("hYBin2");
  //TH1D * hYBin3 = (TH1D*) hj04->Clone("hYBin3");

  TH1D * hBackBin0 = (TH1D*) hj01->Clone("hBackBin0");
  //TH1D * hBackBin1 = (TH1D*) hj02->Clone("hBackBin1");
  //TH1D * hBackBin2 = (TH1D*) hj03->Clone("hBackBin2");
  //TH1D * hBackBin3 = (TH1D*) hj04->Clone("hBackBin3");
  hBackBin0->Eval(fj01back);
  //hBackBin1->Eval(fj02back);
  //hBackBin2->Eval(fj03back);
  //hBackBin3->Eval(fj04back);
  //hYBin0->Add(hYBin0,hBackBin0,1.0,-1.0);
  //hYBin1->Add(hYBin1,hBackBin1,1.0,-1.0);
  //hYBin2->Add(hYBin2,hBackBin2,1.0,-1.0);
  //hYBin3->Add(hYBin3,hBackBin3,1.0,-1.0);

  //B0
  //double errVal_Xi_1530_B0;
  //double yieldVal_Xi_1530_B0 = hYBin0->IntegralAndError(low_Xi_1530_Bin,high_Xi_1530_Bin,errVal_Xi_1530_B0);

  //B1
  //double errVal_Xi_1530_B1;
  //double yieldVal_Xi_1530_B1 = hYBin1->IntegralAndError(low_Xi_1530_Bin,high_Xi_1530_Bin,errVal_Xi_1530_B1);

  //B2
  //double errVal_Xi_1530_B2;
  //double yieldVal_Xi_1530_B2 = hYBin2->IntegralAndError(low_Xi_1530_Bin,high_Xi_1530_Bin,errVal_Xi_1530_B2);

  //B3
  //double errVal_Xi_1530_B3;
  //double yieldVal_Xi_1530_B3 = hYBin3->IntegralAndError(low_Xi_1530_Bin,high_Xi_1530_Bin,errVal_Xi_1530_B3);

  //hYBin0->GetXaxis()->SetRange(low_Xi_1530_Bin,high_Xi_1530_Bin);
  //hYBin1->GetXaxis()->SetRange(low_Xi_1530_Bin,high_Xi_1530_Bin);
  //hYBin2->GetXaxis()->SetRange(low_Xi_1530_Bin,high_Xi_1530_Bin);
  //hYBin3->GetXaxis()->SetRange(low_Xi_1530_Bin,high_Xi_1530_Bin);
  fj01back->SetLineColor(1);
  fj02back->SetLineColor(1);
  fj01back->SetLineWidth(lineWidth);
  fj02back->SetLineWidth(lineWidth);
  fj01back->SetLineStyle(2);
  fj02back->SetLineStyle(2);
  fj01Flatte->SetLineColor(30);
  fj01Flatte->SetLineWidth(2);

  fj02Flatte->GetXaxis()->SetRangeUser(1.6,3.0);
  fj02Flatte->SetLineColor(30);
  fj02Flatte->SetLineWidth(2);

  fj01Peak1->SetLineColor(2);
  //fj01Peak1->SetLineWidth(lineWidth);
  fj01Peak1->SetLineWidth(2);
 
  fj01Peak3->SetLineColor(9);
  fj02Peak3->SetLineColor(9);
  fj01Peak3->SetLineWidth(2);
  fj02Peak3->SetLineWidth(2);

  fj01Peak5->SetLineColor(3);
  fj02Peak5->SetLineColor(3);
  fj01Peak5->SetLineWidth(2);
  fj02Peak5->SetLineWidth(2);
  
  fj01Peak8->SetLineColor(4);
  fj02Peak8->SetLineColor(4);
  fj01Peak8->SetLineWidth(2);
  fj02Peak8->SetLineWidth(2);

  fj01Peak9->SetLineColor(7);
  fj02Peak9->SetLineColor(7);
  fj01Peak9->SetLineWidth(2);
  fj02Peak9->SetLineWidth(2);
  //fj03back->SetLineColor(2);
  //fj04back->SetLineColor(2);
  //if (iShow == 1) {
  c1->cd(1);
  fj01back->Draw("same");


  int peak1Yes = 1;
  int flatteYes = 1;
  int peak2Yes = 0;
  int peak3Yes = 0;
  int peak5Yes = 0;
  int peak8Yes = 0;
  int peak9Yes = 0; if (fit2120 == true) peak9Yes = 1; 
  if (rangeType == 2){
    peak5Yes = 0;
    peak8Yes = 0;
  }

  
  fj01Peak1->GetXaxis()->SetRangeUser(1.4,1.7);
  fj01Flatte->GetXaxis()->SetRangeUser(1.4,1.9);
  fj01Peak3->GetXaxis()->SetRangeUser(1.6,1.75);
  fj01Peak5->GetXaxis()->SetRangeUser(1.7,2.0);
  fj01Peak8->GetXaxis()->SetRangeUser(1.7,2.2);
  fj01Peak9->GetXaxis()->SetRangeUser(1.6,2.6);

  fj02Flatte->GetXaxis()->SetRangeUser(1.4,1.9);
  fj02Peak3->GetXaxis()->SetRangeUser(1.6,1.75);
  fj02Peak5->GetXaxis()->SetRangeUser(1.7,2.0);
  fj02Peak8->GetXaxis()->SetRangeUser(1.7,2.2);
  fj02Peak9->GetXaxis()->SetRangeUser(1.6,2.6);


  if (flatteYes == 1) fj01Flatte->Draw("same");
  if (peak1Yes == 1) fj01Peak1->Draw("same");
  if (peak3Yes == 1) fj01Peak3->Draw("same");
  if (peak5Yes == 1) fj01Peak5->Draw("same");
  if (peak8Yes == 1) fj01Peak8->Draw("same");
  if (peak9Yes == 1) fj01Peak9->Draw("same");


  
  //1530 yield extraction
  TH1D * hj01Back1530 = (TH1D*) hj01->Clone("hj01Back1530");
  TH1D * hj01Signal1530 = (TH1D*) hj01->Clone("hj01Signal1530");
  hj01Back1530->Eval(fj01Bkg1530);
  //hj01Back1530->Sumw2();
  hj01Signal1530->Add(hj01Signal1530,hj01Back1530,1,-1);
  //fj01Bkg1530->Draw("same");
  hj01Signal1530->SetLineColor(2);
  hj01Signal1530->SetMarkerColor(2);
  hj01Signal1530->SetFillColor(2);
  hj01Signal1530->SetFillStyle(3001);
  hj01Signal1530->GetXaxis()->SetRangeUser(rangeLow1530,rangeHigh1530);
  int lBin1530 = hj01Signal1530->FindBin(rangeLow1530);
  int uBin1530 = hj01Signal1530->FindBin(rangeHigh1530);
  if (bNumber == 1) lBin1530 = lBin1530-1;
  if (bNumber == 1) uBin1530 = uBin1530+1;
  if (bNumber > 1) lBin1530 = lBin1530+1;
  rangeLow1530 =  hj01Signal1530->GetBinCenter(lBin1530);
  rangeHigh1530 =  hj01Signal1530->GetBinCenter(uBin1530);
  hj01Signal1530->GetXaxis()->SetRangeUser(rangeLow1530,rangeHigh1530);


  double err1530;
  double yield1530 = hj01Signal1530->IntegralAndError(lBin1530,uBin1530,err1530,"");
  //cout<<"yield1530, err/yield = "<<yield1530<<" +/- "<<err1530<<" , "<<err1530/yield1530<<endl;
  hj01Signal1530->SetLineWidth(lineWidth);
  if (peak1Yes == 1) hj01Signal1530->Draw("same");
  if (peak1Yes == 1) hj01Signal1530->Draw("samehist");
  
  
  double normedYield1530 = yield1530/lumiVal;
  //cout<<"normedYield1530 = "<<normedYield1530<<endl;
  
  double cs1530 = normedYield1530*3.0/(1000*eff1530);
  double cs1530Err = sqrt(pow(cs1530*(err1530/yield1530),2) + pow(cs1530*eff1530Err/eff1530,2));
  cout<<"### 1530 ###\n"<<endl;
  cout<<"\t\tcs1530 (includes factor of 3) = "<<cs1530<<" , "<<cs1530Err<<endl;
  cout<<"############\n"<<endl;

  //1620 yield extraction
  TH1D * hj01Back1620 = (TH1D*) hj01->Clone("hj01Back1620");
  TH1D * hj01Signal1620 = (TH1D*) hj01->Clone("hj01Signal1620");
  hj01Back1620->Eval(fj01Bkg1620);
  //hj01Back1620->Sumw2();
  hj01Signal1620->Add(hj01Signal1620,hj01Back1620,1,-1);
  fj01Bkg1620->SetLineWidth(lineWidth);
  fj01Bkg1620->Draw("same");

  hj01Signal1620->SetLineColor(30);
  hj01Signal1620->SetMarkerColor(30);
  hj01Signal1620->SetFillColor(30);
  hj01Signal1620->SetFillStyle(3001);
  hj01Signal1620->GetXaxis()->SetRangeUser(rangeLow1620,rangeHigh1620);
  int lBin1620 = hj01Signal1620->FindBin(rangeLow1620);
  int uBin1620 = hj01Signal1620->FindBin(rangeHigh1620);
  double err1620;
  double yield1620 = hj01Signal1620->IntegralAndError(lBin1620,uBin1620,err1620,"");
  //cout<<"yield1620, err/yield = "<<yield1620<<" +/- "<<err1620<<" , "<<err1620/yield1620<<endl;
  if (peak2Yes == 1) hj01Signal1620->Draw("same");
  if (peak2Yes == 1) hj01Signal1620->Draw("samehist");

  //1820 yield extraction
  TH1D * hj01Back1820 = (TH1D*) hj01->Clone("hj01Back1820");
  TH1D * hj01Signal1820 = (TH1D*) hj01->Clone("hj01Signal1820");
  hj01Back1820->Eval(fj01Bkg1820);
  //hj01Back1820->Sumw2();
  hj01Signal1820->Add(hj01Signal1820,hj01Back1820,1,-1);
  fj01Bkg1820->SetLineWidth(lineWidth);
  if (peak5Yes == 1) fj01Bkg1820->Draw("same");
  hj01Signal1820->SetLineColor(3);
  hj01Signal1820->SetMarkerColor(3);
  hj01Signal1820->SetFillColor(3);
  hj01Signal1820->SetFillStyle(3001);

  //double rangeLow1820 = 1.78;
  //double rangeHigh1820 = 1.89;
  hj01Signal1820->GetXaxis()->SetRangeUser(rangeLow1820,rangeHigh1820);

  
  int lBin1820 = hj01Signal1820->FindBin(rangeLow1820);
  int uBin1820 = hj01Signal1820->FindBin(rangeHigh1820);
  double err1820;
  double yield1820 = hj01Signal1820->IntegralAndError(lBin1820,uBin1820,err1820,"");
  cout<<"yield1820, err/yield = "<<yield1820<<" +/- "<<err1820<<" , "<<err1820/yield1820<<endl;
  hj01Signal1820->SetLineWidth(lineWidth);
  if (peak5Yes == 1) hj01Signal1820->Draw("same");
  if (peak5Yes == 1) hj01Signal1820->Draw("samehist");

  double normedYield1820 = yield1820/lumiVal;
  //cout<<"normedYield1820 = "<<normedYield1820<<endl;

  
  //2030 yield extraction
  TH1D * hj01Back2030 = (TH1D*) hj01->Clone("hj01Back2030");
  TH1D * hj01Signal2030 = (TH1D*) hj01->Clone("hj01Signal2030");
  hj01Back2030->Eval(fj01Bkg2030);
  hj01Signal2030->Add(hj01Signal2030,hj01Back2030,1,-1);
  fj01Bkg2030->SetLineWidth(lineWidth);
  if (peak8Yes == 1)  fj01Bkg2030->Draw("same");
  hj01Signal2030->SetLineColor(4);
  hj01Signal2030->SetMarkerColor(4);
  hj01Signal2030->SetFillColor(4);
  hj01Signal2030->SetFillStyle(3001);
  hj01Signal2030->GetXaxis()->SetRangeUser(rangeLow2030,rangeHigh2030);

  int lBin2030 = hj01Signal2030->FindBin(rangeLow2030);
  int uBin2030 = hj01Signal2030->FindBin(rangeHigh2030);
  double err2030;
  double yield2030 = hj01Signal2030->IntegralAndError(lBin2030,uBin2030,err2030,"");
  //cout<<"\tyield2030, err/yield = "<<yield2030<<" +/- "<<err2030<<" , "<<err2030/yield2030<<endl;
  if (peak8Yes == 1) hj01Signal2030->Draw("same");
  if (peak8Yes == 1) hj01Signal2030->Draw("samehist");

  double normedYield2030 = yield2030/lumiVal;
  cout<<"normedYield2030 = "<<normedYield2030<<endl;

  double cs2030 = normedYield2030*3.0/(1000*eff2030);
  double cs2030Err = sqrt(pow(cs2030*(err2030/yield2030),2) + pow(cs2030*eff2030Err/eff2030,2));


  //2120 yield extraction
  TH1D * hj01Back2120 = (TH1D*) hj01->Clone("hj01Back2120");
  TH1D * hj01Signal2120 = (TH1D*) hj01->Clone("hj01Signal2120");
  hj01Back2120->Eval(fj01Bkg2120);
  hj01Signal2120->Add(hj01Signal2120,hj01Back2120,1,-1);
  fj01Bkg2120->SetLineWidth(lineWidth);
  if (peak9Yes == 1)  fj01Bkg2120->Draw("same");
  hj01Signal2120->SetLineColor(7);
  hj01Signal2120->SetMarkerColor(7);
  hj01Signal2120->SetFillColor(7);
  hj01Signal2120->SetFillStyle(3001);
  hj01Signal2120->GetXaxis()->SetRangeUser(rangeLow2120,rangeHigh2120);

  int lBin2120 = hj01Signal2120->FindBin(rangeLow2120);
  int uBin2120 = hj01Signal2120->FindBin(rangeHigh2120);
  double err2120;
  double yield2120 = hj01Signal2120->IntegralAndError(lBin2120,uBin2120,err2120,"");
  //cout<<"\tyield2120, err/yield = "<<yield2120<<" +/- "<<err2120<<" , "<<err2120/yield2120<<endl;
  //if (peak9Yes == 1) hj01Signal2120->Draw("same");
  //if (peak9Yes == 1) hj01Signal2120->Draw("samehist");

  double normedYield2120 = yield2120/lumiVal;
  cout<<"normedYield2120 = "<<normedYield2120<<endl;

  Double_t xl1=0.13, yl1 = 0.8, xl2=xl1+0.15, yl2=yl1+0.05;
  TLegend *legend2 = new TLegend(xl1,yl1,xl2,yl2);
  legend2->SetFillColor(0);
  legend2->SetLineColor(0);
  legend2->SetTextSize(0.06);
  legend2->AddEntry(fj01,"Full fit","l");
  legend2->Draw();

  xl1=xl1+0.25;
  xl2=xl2+0.25;
  TLegend *legend = new TLegend(xl1,yl1,xl2,yl2);
  legend->SetFillColor(0);
  legend->SetLineColor(0);
  legend->SetTextSize(0.06);
  legend->AddEntry(fj01back,"Phase space background","l");
  legend->Draw();

  //peak1Yes = 1;
  //peak5Yes = 1;
  //peak8Yes = 1;

  int peakTotYes = peak1Yes+flatteYes+peak3Yes+peak5Yes+peak8Yes+peak9Yes;
  
  xl1 = 0.7;
  xl2=xl1+0.15;
  yl1 = 0.47;
  yl2=yl1+0.4;

  double deltaY = yl2-yl1;
  deltaY *= peakTotYes/6.0;
  yl1 = yl2 - deltaY;
  
  TLegend *legendB = new TLegend(xl1,yl1,xl2,yl2);
  legendB->SetFillColor(0);
  legendB->SetLineColor(0);
  legendB->SetTextSize(0.06);
  if (peak1Yes==1) legendB->AddEntry(hj01Signal1530,"#Xi(1530)","f");
  if (flatteYes==1) legendB->AddEntry(fj01Flatte,"#Xi(1620) Flatte","l");
  if (peak3Yes==1)  legendB->AddEntry(fj01Peak3,"#Xi(1690)","l");
  if (peak5Yes==1)  legendB->AddEntry(hj01Signal1820,"#Xi(1820)","f");
  if (peak8Yes==1)  legendB->AddEntry(hj01Signal2030,"#Xi(2030)","f");
  if (peak9Yes==1)  legendB->AddEntry(fj01Peak9,"#Xi(2120)","l");
  legendB->Draw();



  c1->cd(2);




  fj02back->GetXaxis()->SetRangeUser(1.6,3.0);
  //fj02Peak5->GetXaxis()->SetRangeUser(1.6,3.0);
  //fj02Peak8->GetXaxis()->SetRangeUser(1.6,3.0);
  //fj02Peak9->GetXaxis()->SetRangeUser(1.6,3.0);
  //fj02Flatte->GetXaxis()->SetRangeUser(1.6,3.0);

  fj02back->Draw("same");
  if (flatteYes == 1) fj02Flatte->Draw("same");
  if (peak3Yes == 1) fj02Peak3->Draw("same");
  if (peak5Yes == 1) fj02Peak5->Draw("same");
  if (peak8Yes == 1) fj02Peak8->Draw("same");
  if (peak9Yes == 1) fj02Peak9->Draw("same");

  //1820 yield extraction
  TH1D * hj02Back1820 = (TH1D*) hj02->Clone("hj02Back1820");
  TH1D * hj02Signal1820 = (TH1D*) hj02->Clone("hj02Signal1820");
  hj02Back1820->Eval(fj02Bkg1820);
  //hj02Back1820->Sumw2();
  hj02Signal1820->Add(hj02Signal1820,hj02Back1820,1,-1);
  fj02Bkg1820->GetXaxis()->SetRangeUser(1.6,3.0);
  fj02Bkg1820->SetLineWidth(lineWidth);
  if (peak5Yes == 1) fj02Bkg1820->Draw("same");
  hj02Signal1820->SetLineColor(3);
  hj02Signal1820->SetMarkerColor(3);
  hj02Signal1820->SetFillColor(3);
  hj02Signal1820->SetFillStyle(3001);

  //RESET
  //rangeLow1820 = 1.76;
  //rangeHigh1820 = 1.91;
  lBin1820 = hj01Signal1820->FindBin(rangeLow1820);
  uBin1820 = hj01Signal1820->FindBin(rangeHigh1820);
  hj02Signal1820->GetXaxis()->SetRangeUser(rangeLow1820,rangeHigh1820);

  double err1820_2;
  double yield1820_2 = hj02Signal1820->IntegralAndError(lBin1820,uBin1820,err1820_2,"");
  //cout<<"\tyield1820_2, err/yield = "<<yield1820_2<<" +/- "<<err1820_2<<" , "<<err1820_2/yield1820_2<<endl;

  if (peak5Yes == 1) hj02Signal1820->Draw("same");
  if (peak5Yes == 1)  hj02Signal1820->Draw("samehist");

  double normedYield1820_2 = yield1820_2/lumiVal;
  //cout<<"normedYield1820_2 = "<<normedYield1820_2<<endl;

  // TEMP
  //double eff1820 = eff2030;
  // double eff1820Err = eff2030Err;
  //double eff1820_2 = eff2030_2;
  //double eff1820_2Err = eff2030_2Err;

  int binDalitz1820 = hDalitz->FindBin(1.82);
  int binDalitz2030 = hDalitz->FindBin(2.03);
  double d1820   = hDalitz->GetBinContent(binDalitz1820);
  double d1820_2 = hDalitz2->GetBinContent(binDalitz1820);

  double d2030 = hDalitz->GetBinContent(binDalitz2030);
  double d2030_2 = hDalitz2->GetBinContent(binDalitz2030);

  double effRat1820 = d1820/d1820_2;
  double effRat2030 = d2030/d2030_2;

  cout<<"effRat1820, effRat2030 = "<<effRat1820<<" , "<<effRat2030<<endl;
  //Taking the 2030 as point of unity. Assuming Xi pi at 1820 statys static
  //Now Drop the 1820 efficiency for Lambda K
  eff1820_2 = eff2030_2/(effRat1820/effRat2030);
  eff1820_2Err = eff2030_2Err/(effRat1820/effRat2030);
  
  double cs1820 = 3*normedYield1820/(1000*eff1820);
  double cs1820Err = sqrt(pow(cs1820*(err1820/yield1820),2) + pow(cs1820*eff1820Err/eff1820,2));

  cout<<"cs1820, err1820, yield1820, eff1820Err, eff1820 = "<<cs1820<<" , "<<err1820<<" , "<<eff1820Err<<" , "<< eff1820<<endl;
  
  double cs1820_2 = normedYield1820_2/(1000*eff1820_2);
  double cs1820_2Err = sqrt(pow(cs1820_2*err1820_2/yield1820_2,2) + pow(cs1820_2*eff1820_2Err/eff1820_2,2));

  double ratio1820 = cs1820/cs1820_2;
  double ratio1820Err = cs1820/cs1820_2*sqrt(pow(cs1820Err/cs1820,2)+pow(cs1820_2Err/cs1820_2,2));

  
  cout<<"### 1820 ###\n"<<endl;
  cout<<"\t\tEST: cs1820 (includes factor of 3) = "<<cs1820<<" +/- "<<cs1820Err<<endl;
  cout<<"\t\tEST: cs1820_2 = "<<cs1820_2<<" +/- "<<cs1820_2Err<<endl;
  cout<<"EST: ratio1820 Xi pi / k Lammbda = "<<ratio1820<<" +/- "<<ratio1820Err<<endl;
  cout<<"############\n"<<endl;  

  
  //2030 yield extraction
  TH1D * hj02Back2030 = (TH1D*) hj02->Clone("hj02Back2030");
  TH1D * hj02Signal2030 = (TH1D*) hj02->Clone("hj02Signal2030");
  hj02Back2030->Eval(fj02Bkg2030);
  //hj02Back2030->Sumw2();
  hj02Signal2030->Add(hj02Signal2030,hj02Back2030,1,-1);
  fj02Bkg2030->GetXaxis()->SetRangeUser(1.6,3.0);
  fj02Bkg2030->SetLineWidth(lineWidth);
  if (peak8Yes == 1) fj02Bkg2030->Draw("same");
  hj02Signal2030->SetLineColor(4);
  hj02Signal2030->SetMarkerColor(4);
  hj02Signal2030->SetFillColor(4);
  hj02Signal2030->SetFillStyle(3001);
  hj02Signal2030->GetXaxis()->SetRangeUser(rangeLow2030,rangeHigh2030);
  //int lBin2030 = hj02Signal2030->FindBin(1.97);                                                                                       
  //int uBin2030 = hj02Signal2030->FindBin(2.08);                                                                                       
  double err2030_2;
  double yield2030_2 = hj02Signal2030->IntegralAndError(lBin2030,uBin2030,err2030_2,"");
  //cout<<"\tyield2030_2, err/yield = "<<yield2030_2<<" +/- "<<err2030_2<<" , "<<err2030_2/yield2030_2<<endl;
  if (peak8Yes == 1) hj02Signal2030->Draw("same");
  if (peak8Yes == 1) hj02Signal2030->Draw("samehist");

 double normedYield2030_2 = yield2030_2/lumiVal;
 //cout<<"normedYield2030_2 = "<<normedYield2030_2<<endl;

 double cs2030_2 = normedYield2030_2/(1000*eff2030_2);
 double cs2030_2Err = sqrt(pow(cs2030_2*(err2030_2/yield2030_2),2) + pow(cs2030_2*eff2030_2Err/eff2030_2,2));

  double ratio2030 = cs2030/cs2030_2;
  double ratio2030Err = cs2030/cs2030_2*sqrt(pow(cs2030Err/cs2030,2)+pow(cs2030_2Err/cs2030_2,2));
  
  cout<<"### 2030 ###\n"<<endl;
  cout<<"\t\tcs2030 (includes factor of 3) = "<<cs2030<<" , "<<cs2030Err<<endl;
  cout<<"\t\tcs2030_2 = "<<cs2030_2<<" , "<<cs2030_2Err<<endl;
  cout<<"ratio2030 Xi pi / k Lammbda = "<<ratio2030<<" +/- "<<ratio2030Err<<endl;
  cout<<"############\n"<<endl;


  //2120 yield extraction
  TH1D * hj02Back2120 = (TH1D*) hj02->Clone("hj02Back2120");
  TH1D * hj02Signal2120 = (TH1D*) hj02->Clone("hj02Signal2120");
  hj02Back2120->Eval(fj02Bkg2120);
  //hj02Back2120->Sumw2();
  hj02Signal2120->Add(hj02Signal2120,hj02Back2120,1,-1);
  fj02Bkg2120->GetXaxis()->SetRangeUser(1.6,3.0);
  fj02Bkg2120->SetLineWidth(lineWidth);
  if (peak9Yes == 1) fj02Bkg2120->Draw("same");
  hj02Signal2120->SetLineColor(7);
  hj02Signal2120->SetMarkerColor(7);
  hj02Signal2120->SetFillColor(7);
  hj02Signal2120->SetFillStyle(3001);
  hj02Signal2120->GetXaxis()->SetRangeUser(rangeLow2120,rangeHigh2120);
  //int lBin2120 = hj02Signal2120->FindBin(1.97);                                                                                       
  //int uBin2120 = hj02Signal2120->FindBin(2.08);                                                                                       
  double err2120_2;
  double yield2120_2 = hj02Signal2120->IntegralAndError(lBin2120,uBin2120,err2120_2,"");
  //cout<<"\tyield2120_2, err/yield = "<<yield2120_2<<" +/- "<<err2120_2<<" , "<<err2120_2/yield2120_2<<endl;
  //if (peak9Yes == 1) hj02Signal2120->Draw("same");
  //if (peak9Yes == 1) hj02Signal2120->Draw("samehist");



  
  //hYBin1->SetFillStyle(3001);
    //hYBin1->SetFillColor(3);
    //hYBin1->Draw("samehist");
    //c1->cd(3);
    //fj03back->Draw("same");
    //hYBin2->SetFillStyle(3001);
    //hYBin2->SetFillColor(3);
    //hYBin2->Draw("samehist");
    //c1->cd(4);
    //fj04back->Draw("same");
    //hYBin3->SetFillStyle(3001);
    //hYBin3->SetFillColor(3);
    //hYBin3->Draw("samehist");
  //}

  //cout<<""<<endl;
  //cout<<"yield_Xi_1530 B0 = "<<yieldVal_Xi_1530_B0<<" +/- "<<errVal_Xi_1530_B0<<endl;
  //cout<<"yield_Xi_1530 B1 = "<<yieldVal_Xi_1530_B1<<" +/- "<<errVal_Xi_1530_B1<<endl;
  //cout<<"yield_Xi_1530 B2 = "<<yieldVal_Xi_1530_B2<<" +/- "<<errVal_Xi_1530_B2<<endl;
  //cout<<"yield_Xi_1530 B3 = "<<yieldVal_Xi_1530_B3<<" +/- "<<errVal_Xi_1530_B3<<endl;

  //TH1D *hYield = new TH1D("hYield","",4,-0.5,3.5); 

  //hYield->SetBinContent(1,yieldVal_Xi_1530_B0);
  //hYield->SetBinError(1,errVal_Xi_1530_B0);

  //hYield->SetBinContent(2,yieldVal_Xi_1530_B1);
  //hYield->SetBinError(2,errVal_Xi_1530_B1);

  //hYield->SetBinContent(3,yieldVal_Xi_1530_B2);
  //hYield->SetBinError(3,errVal_Xi_1530_B2);

  //hYield->SetBinContent(4,yieldVal_Xi_1530_B3);
  //hYield->SetBinError(4,errVal_Xi_1530_B3);

  //outFile->cd();
  //hj01->Write();
  //hj02->Write();
  //hj03->Write();
  //hj04->Write();
  //hYield->Write();
  //outFile->Write();

}

dPointsEffs_t getEff1530(int bNumber, double rangeLow, double rangeHigh){  
  char  fId[120];
  char  hId[120];

  //Efficiency                                                                                                                                                                 
  sprintf(fId,".//cm1StarMCBestCLStage3pExpo6R1530E.root");
  TFile *inFileMC = new TFile(fId);
  TH1D *hXStarMC;
  if (bNumber >-1 && bNumber <4) {
    sprintf(hId,"hMassXiStarCut2MCl6V4B%d",bNumber);inFileMC->GetObject(hId,hXStarMC);
  } else {
    sprintf(hId,"hMassXiStarCut2MCl6V4B0123");inFileMC->GetObject(hId,hXStarMC);
  }
  int lBin1530 = hXStarMC->FindBin(rangeLow);
  int uBin1530 = hXStarMC->FindBin(rangeHigh);
  ////////

  //cout<<"file seen & histogram = "<<fId<<" , "<<hId<<endl;
  
  //double errMCSeen1530;
  //double mcSeen1530 = hXStarMC->IntegralAndError(lBin1530,uBin1530,errMCSeen1530,"");
  double mcSeen1530 = hXStarMC->Integral();

  //THROWN
  TFile *inFileMCThrown = new TFile(".//mcOutThrownCasStar_R1530EAll.root");
  TTree *mass_Cm1Star_Tree = (TTree*)inFileMCThrown->Get("mass_Cm1Star_Tree");
  char  tString[120];
  double e1 = 6.2;
  double e2 = 7.6;
  double e3 = 8.3;
  double e4 = 9.0;
  double e5 = 10.0;
  double e6 = 11.4;
  sprintf(tString,"eGamma>%f && eGamma<%f",e1,e5);
  if (bNumber == 0) sprintf(tString,"eGamma>%f && eGamma<%f",e1,e2);
  if (bNumber == 1) sprintf(tString,"eGamma>%f && eGamma<%f",e2,e3);
  if (bNumber == 2) sprintf(tString,"eGamma>%f && eGamma<%f",e3,e4);
  if (bNumber == 3) sprintf(tString,"eGamma>%f && eGamma<%f",e4,e5);
  double mcThrown1530 = mass_Cm1Star_Tree->Draw("eGamma",tString,"goff");

  //cout<<"Thrown File: & command "<<".//mcOutThrownCasStar_R1530EAll.root"<<" , "<<"mass_Cm1Star_Tree->Draw(\"eGamma\",\""<< tString <<"\",\"goff\")"<<endl;

  //cout<<"mcThrown1530 = "<<mcThrown1530<<endl;
  
  double eff1530 = mcSeen1530/mcThrown1530;

  double prob1530 = eff1530;
  double eff1530Err = sqrt(eff1530*(1-eff1530)/mcThrown1530);
  /////
  
  dPointsEffs_t effData;
  
  effData.val = eff1530;
  effData.err = eff1530Err;
  effData.val_2 = -1;
  effData.err_2 = -1; 
  
  return effData;
}

dPointsEffs_t getEff1820(int bNumber, int kmLamNumber, double rangeLow, double rangeHigh){  
  char  fId[120];
  char  hId[120];

  //Efficiency                                                                                                                                                                 
  sprintf(fId,".//cm1StarMCBestCLStage3pExpo6R1820A.root");
  TFile *inFileMC = new TFile(fId);
  TH1D *hXStarMC;
  if (bNumber >-1 && bNumber <4) {
    sprintf(hId,"hMassXiStarCut2MCl6V4B%dSBS",bNumber);inFileMC->GetObject(hId,hXStarMC);
  } else {
    sprintf(hId,"hMassXiStarCut2MCl6V4B0123SBS");inFileMC->GetObject(hId,hXStarMC);
  }
  int lBin1820 = hXStarMC->FindBin(rangeLow);
  int uBin1820 = hXStarMC->FindBin(rangeHigh);
  ////////

  //double errMCSeen1820;
  //double mcSeen1820 = hXStarMC->IntegralAndError(lBin1820,uBin1820,errMCSeen1820,"");
  double mcSeen1820 = hXStarMC->Integral();


  sprintf(fId,".//kmLamMCBestCLStage3R1820A.root");
  TFile *inFileMC_2 = new TFile(fId);
  TH1D *hXStarMC_2;     inFileMC_2->GetObject("",hXStarMC_2);
  if (bNumber >-1 && bNumber <4) {
  sprintf(hId,"hMassXiStarCl4V4SBsubCut%dB%d",kmLamNumber,bNumber); inFileMC_2->GetObject(hId,hXStarMC_2);
  } else {
  sprintf(hId,"hMassXiStarCl4V4SBsubCut%d",kmLamNumber); inFileMC_2->GetObject(hId,hXStarMC_2);
  }
  //double errMCSeen1820_2;
  //double mcSeen1820_2 = hXStarMC_2->IntegralAndError(lBin1820,uBin1820,errMCSeen1820_2,"");
  double mcSeen1820_2 = hXStarMC_2->Integral();

  //THROWN
  char  tString[120];
  double e1 = 6.2;
  double e2 = 7.6;
  double e3 = 8.3;
  double e4 = 9.0;
  double e5 = 10.0;
  double e6 = 11.4;
  sprintf(tString,"eGamma>%f && eGamma<%f",e1,e5);
  if (bNumber == 0) sprintf(tString,"eGamma>%f && eGamma<%f",e1,e2);
  if (bNumber == 1) sprintf(tString,"eGamma>%f && eGamma<%f",e2,e3);
  if (bNumber == 2) sprintf(tString,"eGamma>%f && eGamma<%f",e3,e4);
  if (bNumber == 3) sprintf(tString,"eGamma>%f && eGamma<%f",e4,e5);
  
  TFile *inFileMCThrown = new TFile(".//mcOutThrownCasStar_R1820AAll.root");
  TTree *mass_Cm1Star_Tree = (TTree*)inFileMCThrown->Get("mass_Cm1Star_Tree");
  double mcThrown1820 = mass_Cm1Star_Tree->Draw("eGamma",tString,"goff");
  double eff1820 = mcSeen1820/mcThrown1820;

  double prob1820 = eff1820;
  double eff1820Err = sqrt(eff1820*(1-eff1820)/mcThrown1820);
  
  TFile *inFileMCThrown_2 = new TFile(".//mcOutThrownKmLam_R1820AAll.root");
  TTree *mass_kmLam_Tree = (TTree*)inFileMCThrown_2->Get("mass_kmLam_Tree");
  double mcThrown1820_2 = mass_kmLam_Tree->Draw("eGamma",tString,"goff");
  double eff1820_2 = mcSeen1820_2/mcThrown1820_2;

  double prob1820_2 = eff1820_2;
  double eff1820_2Err = sqrt(eff1820_2*(1-eff1820_2)/mcThrown1820_2);
  //////////////

  //cout<<"TEST: eff1820 "<<eff1820<<" +/- "<<eff1820Err<<endl;
  
  dPointsEffs_t effData;
  
  effData.val = eff1820;
  effData.err = eff1820Err;
  effData.val_2 = eff1820_2;
  effData.err_2 = eff1820_2Err;
  
  return effData;
}

/// 2030 ///

dPointsEffs_t getEff2030(int bNumber, int kmLamNumber, double rangeLow, double rangeHigh){  
  char  fId[120];
  char  hId[120];

  //Efficiency                                                                                                                                                                 
  sprintf(fId,".//cm1StarMCBestCLStage3pExpo6R2030B.root");
  TFile *inFileMC = new TFile(fId);
  TH1D *hXStarMC;
  if (bNumber >-1 && bNumber <4) {
    sprintf(hId,"hMassXiStarCut2MCl6V4B%dSBS",bNumber);inFileMC->GetObject(hId,hXStarMC);
  } else {
    sprintf(hId,"hMassXiStarCut2MCl6V4B0123SBS");inFileMC->GetObject(hId,hXStarMC);
  }
  int lBin2030 = hXStarMC->FindBin(rangeLow);
  int uBin2030 = hXStarMC->FindBin(rangeHigh);
  ////////

  //double errMCSeen2030;
  //double mcSeen2030 = hXStarMC->IntegralAndError(lBin2030,uBin2030,errMCSeen2030,"");
  double mcSeen2030 = hXStarMC->Integral();


  sprintf(fId,".//kmLamMCBestCLStage3R2030ZB.root");
  TFile *inFileMC_2 = new TFile(fId);
  TH1D *hXStarMC_2;     inFileMC_2->GetObject("",hXStarMC_2);
  if (bNumber >-1 && bNumber <4) {
  sprintf(hId,"hMassXiStarCl4V4SBsubCut%dB%d",kmLamNumber,bNumber); inFileMC_2->GetObject(hId,hXStarMC_2);
  } else {
  sprintf(hId,"hMassXiStarCl4V4SBsubCut%d",kmLamNumber); inFileMC_2->GetObject(hId,hXStarMC_2);
  }
  //double errMCSeen2030_2;
  //double mcSeen2030_2 = hXStarMC_2->IntegralAndError(lBin2030,uBin2030,errMCSeen2030_2,"");
  double mcSeen2030_2 = hXStarMC_2->Integral();

  //THROWN
  char  tString[120];
  double e1 = 6.2;
  double e2 = 7.6;
  double e3 = 8.3;
  double e4 = 9.0;
  double e5 = 10.0;
  double e6 = 11.4;
  sprintf(tString,"eGamma>%f && eGamma<%f",e1,e5);
  if (bNumber == 0) sprintf(tString,"eGamma>%f && eGamma<%f",e1,e2);
  if (bNumber == 1) sprintf(tString,"eGamma>%f && eGamma<%f",e2,e3);
  if (bNumber == 2) sprintf(tString,"eGamma>%f && eGamma<%f",e3,e4);
  if (bNumber == 3) sprintf(tString,"eGamma>%f && eGamma<%f",e4,e5);
  
  TFile *inFileMCThrown = new TFile(".//mcOutThrownCasStar_R2030BAll.root");
  TTree *mass_Cm1Star_Tree = (TTree*)inFileMCThrown->Get("mass_Cm1Star_Tree");
  double mcThrown2030 = mass_Cm1Star_Tree->Draw("eGamma",tString,"goff");
  double eff2030 = mcSeen2030/mcThrown2030;

  double prob2030 = eff2030;
  double eff2030Err = sqrt(eff2030*(1-eff2030)/mcThrown2030);
  
  TFile *inFileMCThrown_2 = new TFile(".//mcOutThrownKmLam_R2030ZBAll.root");
  TTree *mass_kmLam_Tree = (TTree*)inFileMCThrown_2->Get("mass_kmLam_Tree");
  double mcThrown2030_2 = mass_kmLam_Tree->Draw("eGamma",tString,"goff");
  double eff2030_2 = mcSeen2030_2/mcThrown2030_2;

  double prob2030_2 = eff2030_2;
  double eff2030_2Err = sqrt(eff2030_2*(1-eff2030_2)/mcThrown2030_2);
  //////////////

  //cout<<"TEST: eff2030 "<<eff2030<<" +/- "<<eff2030Err<<endl;
  
  dPointsEffs_t effData;
  
  effData.val = eff2030;
  effData.err = eff2030Err;
  effData.val_2 = eff2030_2;
  effData.err_2 = eff2030_2Err;
  
  return effData;
}
