#include <iostream>
#include <fstream>

//gROOT->Reset();


//########## chenge these values ##########

// Info to create histo/array for Dlnu
// Run1, Run2, Run3, Run4 or all Run1-4 if run==0
int run = 3;

// Info of input histograms
TString runNum("Run3"); //Run1,Run2, Run3, Run4 or RunAll
TString corrEff("noCorr"); // "noCorr" or "effCorr"
TString corrDstarFF("oldFF"); //"newFF","newFFRhoMax","newFFR1Max"..,or "oldFF"
TString DlnuFF("rho0000"); // Dlnu FF slope of input histogram, "rho0000"
TString Slope("New");  // Dlnu FF slope parametrization, "Old" or "New"
TString corrBF("oldBF"); // D->Xenu BF correction, "oldBF" or "newBF"
//TString histDir("histogramsFit"); //directory of histograms, "histogramsFit"

//#########################################

// Setup file including parameters
#include "fit_3DMC_setup.h"

extern void fcnBToD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);

// Declear common variables
static ncount = 0;

double num_3D[nFiles][nbinsx][nbinsy][nbinsz];
double err_3D[nFiles][nbinsx][nbinsy][nbinsz];
enum { nBins = 2000 };  //number of bins of slope : 2000
double lowEnd[nBins];
double hiEnd[nBins];
double RN[nBins];
double RN_err[nBins];
int showChisq = 0;  // printout details of chisquare velue.
int showChisq2 = 0;  // printout details of chisquare velue.
int printChisq = 0;  // printout chi-sq values
int printRN = 0;  // printout nomalization factor values
int showBins = 1;  // print number of bins used in the fit

static Double_t vstrt[19] = {0.0210, 0.0560,
			     0.0020, 0.0056, 0.0037, 0.0037,
			     0.0060, 0.0020,
			     1.35,
			     1, 1, 1, 1,
			     1, 1, 1, 1,
			     1, 1};
/*static Double_t vstrt[19] = {0.02627, 0.05813,
			     0.00813, 0.000001, 0.00065, 0.00100,
			     0.000001, 0.00817,
			     1.5, //1.471,
			     0.9625, 0.9310, 0.2284, 0.9955,
			     0.7073, 0.7538, 0.9108, 1.1780,
			     1.0407, 0.9718};*/
static Double_t stp[19] = {0.001, 0.001, 
			   0.001, 0.001, 0.001, 0.001, 
			   0.001, 0.001,
			   0.2,
			   0.01, 0.01, 0.01, 0.01,
			   0.01, 0.01, 0.01, 0.01,
			   0.01, 0.01};

int chisqMode = 0;  // 0 = all, 1 = D0 only, 2 = D+ only
int doMigrad = 0;  // do MIGRAD
int doMinos = 0;  // do MINOS
int printInfo = 1; // printout chisq and RN when "CALL FCN"
int printDetailInfo = 0; // printout details of chisq when "CALL FCN"
int printDetailInfo2 = 0; // printout calculation of chisq when "CALL FCN"
int printInfoSBS = 0; // printout chisq and RN when "MIGRAD"
int fixSlope = 2; // fix FF slope. If = 2: use single histogram for Dlnu
int fixDlnu = 0;
int fixDstar = 0;
int fixD0star = 0;
double ratio_D0star = 0.35714; //=20/56  Br_D0star = ratio_D0star*Br_D1
int fixD1 = 0;
double Br_D1 = 0.0056; //if fixD1 = 2
int fixD1p = 0;
double ratio_D1p = 0.66071; //=37/56  Br_D1p = ratio_D1p*Br_D1
int fixD2star = 0;
double ratio_D2star = 1.85; //=37/20  Br_D2star = ratio_D2star*Br_D0star
int fixDPi = 0;
double Br_DPi = 0.0060;
int fixDstarPi = 0;
double ratio_DstarPi = 1.24074; //=67/54  Br_DstarPi = ratio_DstarPi*Br_DPi
int fixBkgd = 0; // fix background coeeficients
int fixFp0 = 0;  // fix f_+0
int fixCD0p = 0;  // fix c_D_0+
int fixT0p = 1;  // fix t_0+ = 1 otherwise t_0+ = 1/t_p0
int nCandMin = 25;  // minimum number of events in one bin
double inclBF = 0.104;  //inclusive b->clnu BF


//###########################################
//########## Main ###########################
//###########################################

void fit_3DMC() {

  cout << "start" << endl;
  // Open input files and get 3D histos
  TFile *file[nFiles];
  TH3D * inHist[nFiles];
  for (int l=0; l<nFiles; l++) {
    file[l] = new TFile(inFile[l]);
    //file[l]->ls();
    inHist[l] = (TH3D *)file[l]->Get("sub;1");
  }

  // Readout 3D histos into arrays
  cout << "reading histos into arrays" << endl;
  for (int l=0; l<nFiles; l++) {
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  num_3D[l][i][j][k] = inHist[l]->GetBinContent(i+1,j+1,k+1);
	  err_3D[l][i][j][k] = inHist[l]->GetBinError(i+1,j+1,k+1);
	}
      }
    }
  }
  /*for (int l=0; l<nFiles; l++) {
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  cout << "num_3D[0]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << num_3D[0][i][j][k] << endl;
	  cout << "err_3D[0]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << err_3D[0][i][j][k] << endl;
	  cout << "num_3D[28]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << num_3D[0][i][j][k] << endl;
	  cout << "err_3D[28]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << err_3D[0][i][j][k] << endl;
	}
      }
    }
    }*/

  // Readin integration table (RN table) ########
  TString inputTable("integFF_Dlnu_oldSlope.txt");
  if (Slope == "New") inputTable = "integFF_Dlnu_newSlope.txt"; 

  cout << "reading table : " << inputTable << endl;
  ifstream inputTableFile(inputTable);  //open input file
  if (!inputTableFile) {
    cout << "Cannot open file : " << inputTable << endl;
  }
  for (Int_t j=0;j<nBins;j++) {
    inputTableFile >> lowEnd[j];
    inputTableFile >> hiEnd[j];
    inputTableFile >> RN[j];
    inputTableFile >> RN_err[j];
  }
  inputTableFile.close();
  /*
  for (Int_t j=0;j<10;j++) {
    cout << lowEnd[j] << "  " << hiEnd[j] << "  "
	 << RN[j] << "  " << RN_err[j] << endl;
  }
  for (Int_t j=1000;j<1010;j++) {
    cout << lowEnd[j] << "  " << hiEnd[j] << "  "
	 << RN[j] << "  " << RN_err[j] << endl;
  }
  */

  //######## Fitting ######################
  //#######################################
  
  static Int_t nprm[19] = {0, 1, 
			   2, 3, 4, 5, 
			   6, 7, 
			   8, 
			   9, 10, 11, 12, 
			   13, 14, 15, 16, 
			   17, 18};

  static Double_t p0=0;
  static Double_t p1=1;
  static Double_t p2=2;
  static Double_t p3=3;
  static Double_t p4=4;
  static Double_t p5=5;
  static Double_t p6=6;
  static Double_t p7=7;
  static Double_t p8=8;
  static Double_t p9=9;
  static Double_t p10=10;
  static Double_t p11=11;
  static Double_t p12=12;
  static Double_t p13=13;
  static Double_t p14=14;
  static Double_t p15=15;
  static Double_t p16=16;
  static Double_t p17=17;
  static Double_t p18=18;
  static Double_t p19=19;
  static Double_t p20=20;
  Int_t ierflg = 0;

  TStopwatch timer;
   
  //initialize TMinuit with a maximum of 19 params
  TMinuit *gMinuit = new TMinuit(19);
  cout<<"Starting timer"<<endl;
  timer.Start();
  gMinuit->mninit(5,6,7);
  gMinuit->SetFCN(fcnBToD);

  // Set parameters
  gMinuit->mnparm(nprm[0], "B_B1",  vstrt[0], stp[0], 0,1,ierflg);
  gMinuit->mnparm(nprm[1], "B_B2",  vstrt[1], stp[1], 0,1,ierflg);
  gMinuit->mnparm(nprm[2], "B_B31", vstrt[2], stp[2], 0,1,ierflg);
  gMinuit->mnparm(nprm[3], "B_B32", vstrt[3], stp[3], 0,1,ierflg);
  gMinuit->mnparm(nprm[4], "B_B33", vstrt[4], stp[4], 0,1,ierflg);
  gMinuit->mnparm(nprm[5], "B_B34", vstrt[5], stp[5], 0,1,ierflg);
  gMinuit->mnparm(nprm[6], "B_B41", vstrt[6], stp[6], 0,1,ierflg);
  gMinuit->mnparm(nprm[7], "B_B42", vstrt[7], stp[7], 0,1,ierflg);
  gMinuit->mnparm(nprm[8], "FF slope", vstrt[8], stp[8], -1,4,ierflg);
  gMinuit->mnparm(nprm[9], "Uncorr(D0)", vstrt[9], stp[9], 0,2,ierflg);
  gMinuit->mnparm(nprm[10], "CaslcL(D0)", vstrt[10], stp[10], 0,2,ierflg);
  gMinuit->mnparm(nprm[11], "FakeL(D0)", vstrt[11], stp[11], 0,2,ierflg);
  gMinuit->mnparm(nprm[12], "ContD(D0)", vstrt[12], stp[12], 0,2,ierflg);
  gMinuit->mnparm(nprm[13], "Uncorr(Dch)", vstrt[13], stp[13], 0,2,ierflg);
  gMinuit->mnparm(nprm[14], "CaslcL(Dch)", vstrt[14], stp[14], 0,2,ierflg);
  gMinuit->mnparm(nprm[15], "FakeL(Dch)", vstrt[15], stp[15], 0,2,ierflg);
  gMinuit->mnparm(nprm[16], "ContD(Dch)", vstrt[16], stp[16], 0,2,ierflg);
  gMinuit->mnparm(nprm[17], "f_+0", vstrt[17], stp[13], 0,2,ierflg);
  gMinuit->mnparm(nprm[18], "c_D+0", vstrt[18], stp[15], 0,2,ierflg);
  if (ierflg) {
    Printf(" UNABLE TO DEFINE PARAMETER NO.");
    return ierflg;
  }
  
  // Perform fitting ############
  if (printInfo) {printChisq = 1; printRN = 1;}
  if (printDetailInfo) showChisq = 1;
  if (printDetailInfo2) showChisq2 = 1;
  gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg);
  printChisq = 0;
  printRN = 0;
  showChisq = 0;
  showChisq2 = 0;

  gMinuit->mnexcm("SET PRINT", &p0 ,1,ierflg);

  if (fixSlope) gMinuit->mnexcm("FIX", &p9 ,1,ierflg);
  if (fixDlnu) gMinuit->mnexcm("FIX", &p1 ,1,ierflg);
  if (fixDstar) gMinuit->mnexcm("FIX", &p2 ,1,ierflg);
  if (fixD0star) gMinuit->mnexcm("FIX", &p3 ,1,ierflg);
  if (fixD1) gMinuit->mnexcm("FIX", &p4 ,1,ierflg);
  if (fixD1p) gMinuit->mnexcm("FIX", &p5 ,1,ierflg);
  if (fixD2star) gMinuit->mnexcm("FIX", &p6 ,1,ierflg);
  if (fixDPi) gMinuit->mnexcm("FIX", &p7 ,1,ierflg);
  if (fixDstarPi) gMinuit->mnexcm("FIX", &p8 ,1,ierflg);
  if (fixBkgd) {
    gMinuit->mnexcm("FIX", &p10 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p11 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p12 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p13 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p14 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p15 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p16 ,1,ierflg);
    gMinuit->mnexcm("FIX", &p17 ,1,ierflg);
  }
  if (fixFp0) gMinuit->mnexcm("FIX", &p18,1,ierflg);
  if (fixCD0p) gMinuit->mnexcm("FIX", &p19,1,ierflg);

  if (doMigrad) {
    if (printInfoSBS) {printChisq = 1; printRN = 1;}
    gMinuit->mnexcm("MIGRAD", &p0 ,0,ierflg);
    gMinuit->mnexcm("SHOW COVARIANCE", &p0 ,0,ierflg);
    printChisq = 0;
    printRN = 0;
  }
  if (doMinos) gMinuit->mnexcm("MINOS",  &p0 ,0,ierflg);

  if (printInfo) {printChisq = 1; printRN = 1;}
  //if (printDetailInfo) showChisq = 1;
  //if (printDetailInfo2) showChisq2 = 1;
  gMinuit->mnexcm("CALL FCN", &p2 , 1,ierflg);
  printChisq = 0;
  printRN = 0;
  showChisq = 0;
  showChisq2 = 0;

  cout<<"Time at the end of job = "<<timer.CpuTime()<<" seconds"<<endl;

  // Now Calcualte branching fractions from fitted parameters
  // Get parameters
  Double_t BF1(0), BF1_err(0);
  Double_t BF2(0), BF2_err(0);
  Double_t BF31(0), BF31_err(0);
  Double_t BF32(0), BF32_err(0);
  Double_t BF33(0), BF33_err(0);
  Double_t BF34(0), BF34_err(0);
  Double_t BF41(0), BF41_err(0);
  Double_t BF42(0), BF42_err(0);
  Double_t slope(0), slopeErr(0);
  Double_t BkgD01(0), BkgD01_err(0);
  Double_t BkgD02(0), BkgD02_err(0);
  Double_t BkgD03(0), BkgD03_err(0);
  Double_t BkgD04(0), BkgD04_err(0);
  Double_t BkgDch1(0), BkgDch1_err(0);
  Double_t BkgDch2(0), BkgDch2_err(0);
  Double_t BkgDch3(0), BkgDch3_err(0);
  Double_t BkgDch4(0), BkgDch4_err(0);
  Double_t F_p0(0), F_p0_err(0);
  Double_t RDp0(0), RDp0_err(0);
  gMinuit->GetParameter(0, BF1, BF1_err);
  gMinuit->GetParameter(1, BF2, BF2_err);
  gMinuit->GetParameter(2, BF31, BF31_err);
  gMinuit->GetParameter(3, BF32, BF32_err);
  gMinuit->GetParameter(4, BF33, BF33_err);
  gMinuit->GetParameter(5, BF34, BF34_err);
  gMinuit->GetParameter(6, BF41, BF41_err);
  gMinuit->GetParameter(7, BF42, BF42_err);
  gMinuit->GetParameter(8, slope, slopeErr);
  gMinuit->GetParameter(9, BkgD01, BkgD01_err);
  gMinuit->GetParameter(10, BkgD02, BkgD02_err);
  gMinuit->GetParameter(11, BkgD03, BkgD03_err);
  gMinuit->GetParameter(12, BkgD04, BkgD04_err);
  gMinuit->GetParameter(13, BkgDch1, BkgDch1_err);
  gMinuit->GetParameter(14, BkgDch2, BkgDch2_err);
  gMinuit->GetParameter(15, BkgDch3, BkgDch3_err);
  gMinuit->GetParameter(16, BkgDch4, BkgDch4_err);
  gMinuit->GetParameter(17, F_p0, F_p0_err);
  gMinuit->GetParameter(18, RDp0, RDp0_err);
  cout << "BF1 = " << BF1 << endl;
  cout << "BF2 = " << BF2 << endl;
  cout << "BF31 = " << BF31 << endl;
  cout << "BF32 = " << BF32 << endl;
  cout << "BF33 = " << BF33 << endl;
  cout << "BF34 = " << BF34 << endl;
  cout << "BF41 = " << BF41 << endl;
  cout << "BF42 = " << BF42 << endl;
  cout << "BF1_err = " << BF1_err << endl;
  cout << "BF2_err = " << BF2_err << endl;
  cout << "BF31_err = " << BF31_err << endl;
  cout << "BF32_err = " << BF32_err << endl;
  cout << "BF33_err = " << BF33_err << endl;
  cout << "BF34_err = " << BF34_err << endl;
  cout << "BF41_err = " << BF41_err << endl;
  cout << "BF42_err = " << BF42_err << endl;
  cout << "ratio_D0star = " << ratio_D0star << endl;
  cout << "ratio_D1p = " << ratio_D1p << endl;
  cout << "ratio_D2star = " << ratio_D2star << endl;
  cout << "ratio_DstarPi = " << ratio_DstarPi << endl;
  if (fixD0star)  BF31 = ratio_D0star * BF32;
  if (fixD1p)  BF33 = ratio_D1p * BF32;
  if (fixD2star)  BF34 = ratio_D2star * BF31;
  if (fixDstarPi)  BF42 = ratio_DstarPi * BF41;
  cout << "BF1 = " << BF1 << endl;
  cout << "BF2 = " << BF2 << endl;
  cout << "BF31 = " << BF31 << endl;
  cout << "BF32 = " << BF32 << endl;
  cout << "BF33 = " << BF33 << endl;
  cout << "BF34 = " << BF34 << endl;
  cout << "BF41 = " << BF41 << endl;
  cout << "BF42 = " << BF42 << endl;

  // Calcurate BFs
  double TotFrac 
    = BF1 + BF2 + BF31 + BF32 + BF33 + BF34 + BF41*1.5 + BF42*1.5;
  double denom = TotFrac;
  cout << "TotFrac = " << TotFrac << endl;
  double BF_B1(0), err_B1(0);
  double BF_B2(0), err_B2(0);
  double BF_B31(0), err_B31(0);
  double BF_B32(0), err_B32(0);
  double BF_B33(0), err_B33(0);
  double BF_B34(0), err_B34(0);
  double BF_B41(0), err_B41(0);
  double BF_B42(0), err_B42(0);
  double BF_B41b(0);
  double BF_B42b(0);
  cout << "*************************" << endl;
  cout << "** Branching Fractions **" << endl;
  cout << "*************************" << endl;
  BF_B1 = BF1 * inclBF / denom;
  err_B1 = BF1_err * inclBF / denom;
  double percent = BF1_err * 100 / BF1;
  cout << "B(B->Dlnu) = " << BF_B1 << " +- " << err_B1 
       << " (" << percent << " %)" << endl;
  BF_B2 = BF2 * inclBF / denom;
  err_B2 = BF2_err * inclBF / denom;
  double percent = BF2_err * 100 / BF2;
  cout << "B(B->D*lnu) = " << BF_B2 << " +- " << err_B2 
       << " (" << percent << " %)" << endl;
  BF_B31 = BF31 * inclBF / denom;
  err_B31 = BF31_err * inclBF / denom;
  double percent = BF31_err * 100 / BF31;
  cout << "B(B->D0*lnu) = " << BF_B31 << " +- " << err_B31 
       << " (" << percent << " %)" << endl;
  BF_B32 = BF32 * inclBF / denom;
  err_B32 = BF32_err * inclBF / denom;
  double percent = BF32_err * 100 / BF32;
  cout << "B(B->D1lnu) = " << BF_B32 << " +- " << err_B32 
       << " (" << percent << " %)" << endl;
  BF_B33 = BF33 * inclBF / denom;
  err_B33 = BF33_err * inclBF / denom;
  double percent = BF33_err * 100 / BF33;
  cout << "B(B->D1'lnu) = " << BF_B33 << " +- " << err_B33 
       << " (" << percent << " %)" << endl;
  BF_B34 = BF34 * inclBF / denom;
  err_B34 = BF34_err * inclBF / denom;
  double percent = BF34_err * 100 / BF34;
  cout << "B(B->D2*lnu) = " << BF_B34 << " +- " << err_B34 
       << " (" << percent << " %)" << endl;
  BF_B41 = BF41 * inclBF / denom;
  err_B41 = BF41_err * inclBF / denom;
  double percent = BF41_err * 100 / BF41;
  cout << "B(B->DPi+lnu) = " << BF_B41 << " +- " << err_B41 
       << " (" << percent << " %)" << endl;
  BF_B41b = BF_B41/2;
  cout << "B(B->DPi0lnu) = " << BF_B41b << endl;
  BF_B42 = BF42 * inclBF / denom;
  err_B42 = BF42_err * inclBF / denom;
  double percent = BF42_err * 100 / BF42;
  cout << "B(B->D*Pi+lnu) = " << BF_B42 << " +- " << err_B42 
       << " (" << percent << " %)" << endl;
  BF_B42b = BF_B42/2;
  cout << "B(B->D*Pi0lnu) = " << BF_B42b << endl;
  double BRSum 
    = BF_B1 + BF_B2 + BF_B31 + BF_B32 + BF_B33 + BF_B34
    + BF_B41 + BF_B41b + BF_B42 + BF_B42b;
  cout << "Sum of BR's = " << BRSum << endl;

  double percent(0);
  if (slope>0) percent = slopeErr * 100 / slope;
  cout << "slope = " << slope << " +- " << slopeErr 
       << " (" << percent << " %)" << endl;

  double percent(0);
  if (BkgD01>0) percent = BkgD01_err * 100 / BkgD01;
  cout << "BkgD0(Uncorr) = " << BkgD01 << " +- " << BkgD01_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgD02>0) percent = BkgD02_err * 100 / BkgD02;
  cout << "BkgD0(Casc L) = " << BkgD02 << " +- " << BkgD02_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgD03>0) percent = BkgD03_err * 100 / BkgD03;
  cout << "BkgD0(Fake L) = " << BkgD03 << " +- " << BkgD03_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgD04>0) percent = BkgD04_err * 100 / BkgD04;
  cout << "BkgD0(Cont D) = " << BkgD04 << " +- " << BkgD04_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgDch1>0) percent = BkgDch1_err * 100 / BkgDch1;
  cout << "BkgDch(Uncorr) = " << BkgDch1 << " +- " << BkgDch1_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgDch2>0) percent = BkgDch2_err * 100 / BkgDch2;
  cout << "BkgDch(Casc L) = " << BkgDch2 << " +- " << BkgDch2_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgDch3>0) percent = BkgDch3_err * 100 / BkgDch3;
  cout << "BkgDch(Fake L) = " << BkgDch3 << " +- " << BkgDch3_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (BkgDch4>0) percent = BkgDch4_err * 100 / BkgDch4;
  cout << "BkgDch(Cont D) = " << BkgDch4 << " +- " << BkgDch4_err 
       << " (" << percent << " %)" << endl;

  double percent(0);
  if (F_p0>0) percent = F_p0_err * 100 / F_p0;
  cout << "f_+0 = f+-/f00 = " << F_p0 << " +- " << F_p0_err 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (RDp0>0) percent = RDp0_err * 100 / RDp0;
  cout << "c_D+0 = D+-/D00 = " << RDp0 << " +- " << RDp0_err 
       << " (" << percent << " %)" << endl;
  
}



//############################################################
//########## Function to calculate chi-sq ####################
//############################################################

void fcnBToD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag)
{
  //cout << "gin = " << gin << endl;
  //cout << "x[0] = " << x[0] << endl;
  //cout << "x[1] = " << x[1] << endl;
  //cout << "x[9] = " << x[9] << endl;

  // Branching fractions
  Double_t BrB1, BrB2;
  Double_t BrB31, BrB32, BrB33, BrB34;
  Double_t BrB41, BrB42;
  BrB1 = x[0];
  BrB2 = x[1];
  if (fixD0star == 2) BrB31 = x[3]*ratio_D0star;
  else BrB31 = x[2];
  if (fixD1 == 2)   BrB32 = Br_D1;
  else BrB32 = x[3];
  if (fixD1p == 2)   BrB33 = x[3]*ratio_D1p;
  else BrB33 = x[4];
  if (fixD2star == 2)   BrB34 = x[2]*ratio_D2star;
  else BrB34 = x[5];
  if (fixDPi == 2)   BrB41 = Br_DPi;
  else BrB41 = x[6];
  if (fixDstarPi == 2)   BrB42 = x[6]*ratio_DstarPi;
  else BrB42 = x[7];

  // FF slope
  Double_t rho_D_sq;
  rho_D_sq = x[8];

  // Peaking background
  Double_t Bkg1_D0, Bkg2_D0, Bkg3_D0, Bkg4_D0;
  Bkg1_D0 = x[9];
  Bkg2_D0 = x[10];
  Bkg3_D0 = x[11];
  Bkg4_D0 = x[12];
  Double_t Bkg1_Dch, Bkg2_Dch, Bkg3_Dch, Bkg4_Dch;
  Bkg1_Dch = x[13];
  Bkg2_Dch = x[14];
  Bkg3_Dch = x[15];
  Bkg4_Dch = x[16];
  
  // ratio f_+- / f_00 and t_+- / t_00
  Double_t fp0, t0p;
  fp0 = x[17];
  if (fixT0p) t0p = 1;
  else t0p = 1/t_p0;

  // ratio D+/D0
  Double_t cDp0;
  cDp0 = x[18];

  // Coefficients ################
  double c[nFiles];
  // D0KPi
  c[0] = 1; //OnPeak
  c[1] = 1; //OffPeak
  c[2] = 1; //uds_Cont
  c[3] = 1; //tautau_Cont
  c[4] = 1; //ccbar_Cont
  c[5] = Bkg4_D0; //ccbar_RealD
  c[6] = 1 * fp0; //BpBm_Comb
  c[7] = 1;  //B0B0_Comb
  c[8] = Bkg1_D0 * fp0; //BpBm_DiffB
  c[9] = Bkg1_D0; //B0B0_DiffB
  c[10] = Bkg2_D0 * fp0; //BpBm_CascL
  c[11] = Bkg2_D0; //B0B0_CascL
  c[12] = Bkg3_D0 * fp0;  //BpBm_LMisID
  c[13] = Bkg3_D0; //B0B0_LMisID

  double B_new = fp0*(B_Ddstar11 / 2 + B_Ddstar12 * (1/2 + B_Dstarp1)) * BrB31;
  double B_old = f_p0_SP6*(B_Ddstar11_SP6 / 2 + B_Ddstar12_SP6 * (1/2 + B_Dstarp1_SP6)) * B_B31_SP6;
  c[14] = B_new / B_old; //BpBm_D0star
  
  double B_new = t0p*(B_Ddstar11 + B_Ddstar12 * (1 + B_Dstarp1 / 2)) * BrB31;
  double B_old = (B_Ddstar11_SP6 + B_Ddstar12_SP6 * (1 + B_Dstarp1_SP6 / 2)) * B_B31_SP6;
  c[15] = B_new / B_old; //B0B0_D0star

  double B_new = fp0*(B_Ddstar21 / 2 + B_Ddstar22 * (1/2 + B_Dstarp1)) * BrB32;
  double B_old = f_p0_SP6*(B_Ddstar21_SP6 / 2 + B_Ddstar22_SP6 * (1/2 + B_Dstarp1_SP6)) * B_B32_SP6;
  c[16] = B_new / B_old; //BpBm_D1

  double B_new = t0p*(B_Ddstar21 + B_Ddstar22 * (1 + B_Dstarp1 / 2)) * BrB32;
  double B_old = (B_Ddstar21_SP6 + B_Ddstar22_SP6 * (1 + B_Dstarp1_SP6 / 2)) * B_B32_SP6;
  c[17] = B_new / B_old;  //B0B0_D1

  double B_new = fp0*(B_Ddstar31 / 2 + B_Ddstar32 * (1/2 + B_Dstarp1)) * BrB33;
  double B_old = f_p0_SP6*(B_Ddstar31_SP6 / 2 + B_Ddstar32_SP6 * (1/2 + B_Dstarp1_SP6)) * B_B33_SP6;
  c[18] = B_new / B_old; //BpBm_D1prime

  double B_new = t0p*(B_Ddstar31 + B_Ddstar32 * (1 + B_Dstarp1 / 2)) * BrB33;
  double B_old = (B_Ddstar31_SP6 + B_Ddstar32_SP6 * (1 + B_Dstarp1_SP6 / 2)) * B_B33_SP6;
  c[19] = B_new / B_old; //B0B0_D1prime

  double B_new = fp0*(B_Ddstar41 / 2 + B_Ddstar42 * (1/2 + B_Dstarp1)) * BrB34;
  double B_old = f_p0_SP6*(B_Ddstar41_SP6 / 2 + B_Ddstar42_SP6 * (1/2 + B_Dstarp1_SP6)) * B_B34_SP6;
  c[20] = B_new / B_old; //BpBm_D2star

  double B_new = t0p*(B_Ddstar41 + B_Ddstar42 * (1 + B_Dstarp1 / 2)) * BrB34;
  double B_old = (B_Ddstar41_SP6 + B_Ddstar42_SP6 * (1 + B_Dstarp1_SP6 / 2)) * B_B34_SP6;
  c[21] = B_new / B_old; //B0B0_D2star

  double B_new = (fp0 / 2) * BrB41;
  double B_old = (f_p0_SP6 / 2) * B_B41_SP6;
  c[22] = B_new / B_old; //BpBm_DPi

  double B_new = t0p * BrB41;
  double B_old = (1 / t_p0_SP6) * B_B41_SP6;
  c[23] = B_new / B_old; //B0B0_DPi

  double B_new = fp0 * (1/2 + B_Dstarp1) * BrB42;
  double B_old = f_p0_SP6 * (1/2 + B_Dstarp1_SP6) * B_B42_SP6;
  c[24] = B_new / B_old; //BpBm_DstarPi

  double B_new = (1 + B_Dstarp1 / 2) * t0p * BrB42;
  double B_old = ((1 + B_Dstarp1_SP6 / 2) / t_p0_SP6) * B_B42_SP6;
  c[25] = B_new / B_old; //B0B0_DstarPi

  double B_new = fp0 * BrB2;
  double B_old = f_p0_SP6 * B_B2_SP6;
  c[26] = B_new / B_old; //BpBm_Dstar

  double B_new = B_Dstarp1 * t0p * BrB2;
  double B_old = (B_Dstarp1_SP6 / t_p0_SP6) * B_B2_SP6;
  c[27] = B_new / B_old; //B0B0_Dstar

  double B_new = fp0 * BrB1;
  double B_old = f_p0_SP6 * B_B1_SP6;
  c[28] = B_new / B_old; //BpBm_Dlnu

  c[29] = 1;  //BpBm_Dlnu_0th
  c[30] = 1;  //BpBm_Dlnu_1st
  c[31] = 1;  //BpBm_Dlnu_2nd

  // DchKPiPi
  c[32] = 1;  //OnPeak
  c[33] = 1; //OffPeak
  c[34] = cDp0; //uds_Cont
  c[35] = cDp0; //tautau_Cont
  c[36] = cDp0; //ccbar_Cont
  c[37] = cDp0 * Bkg4_Dch; //ccbar_RealD
  c[38] = cDp0 * fp0; //BpBm_Comb
  c[39] = cDp0; //B0B0_Comb
  c[40] = cDp0 * Bkg1_Dch * fp0; //BpBm_DiffB
  c[41] = cDp0 * Bkg1_Dch; //B0B0_DiffB
  c[42] = cDp0 * Bkg2_Dch * fp0; //BpBm_CascL
  c[43] = cDp0 * Bkg2_Dch; //B0B0_CascL
  c[44] = cDp0 * Bkg3_Dch * fp0; //BpBm_LMisID
  c[45] = cDp0 * Bkg3_Dch; //B0B0_LMisID

  double B_new = fp0 * (B_Ddstar11 + B_Ddstar12 *( B_Dstarp1 + B_Dstar02)) * BrB31;
  double B_old = f_p0_SP6 * (B_Ddstar11_SP6 + B_Ddstar12_SP6 *( B_Dstarp1_SP6 + B_Dstar02_SP6)) * B_B31_SP6;
  c[46] = cDp0 * B_new / B_old; //BpBm_D0star

  double B_new = ((B_Ddstar11 + B_Ddstar12 * B_Dstarp1) / 2) * t0p * BrB31;
  double B_old = ((B_Ddstar11_SP6 + B_Ddstar12_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6)) * B_B31_SP6;
  c[47] = cDp0 * B_new / B_old; //B0B0_D0star

  double B_new = fp0 * (B_Ddstar21 + B_Ddstar22 *( B_Dstarp1 + B_Dstar02)) * BrB32;
  double B_old = f_p0_SP6 * (B_Ddstar21_SP6 + B_Ddstar22_SP6 *( B_Dstarp1_SP6 + B_Dstar02_SP6)) * B_B32_SP6;
  c[48] = cDp0 * B_new / B_old; //BpBm_D1

  double B_new = ((B_Ddstar21 + B_Ddstar22 * B_Dstarp1) / 2) * t0p * BrB32;
  double B_old = ((B_Ddstar21_SP6 + B_Ddstar22_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6)) * B_B32_SP6;
  c[49] = cDp0 * B_new / B_old; //B0B0_D1

  double B_new = fp0 * (B_Ddstar31 + B_Ddstar32 *( B_Dstarp1 + B_Dstar02)) * BrB33;
  double B_old = f_p0_SP6 * (B_Ddstar31_SP6 + B_Ddstar32_SP6 *( B_Dstarp1_SP6 + B_Dstar02_SP6)) * B_B33_SP6;
  c[50] = cDp0 * B_new / B_old; //BpBm_D1prime

  double B_new = ((B_Ddstar31 + B_Ddstar32 * B_Dstarp1) / 2) * t0p * BrB33;
  double B_old = ((B_Ddstar31_SP6 + B_Ddstar32_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6)) * B_B33_SP6;
  c[51] = cDp0 * B_new / B_old; //B0B0_D1prime

  double B_new = fp0 * (B_Ddstar41 + B_Ddstar42 *( B_Dstarp1 + B_Dstar02)) * BrB34;
  double B_old = f_p0_SP6 * (B_Ddstar41_SP6 + B_Ddstar42_SP6 *( B_Dstarp1_SP6 + B_Dstar02_SP6)) * B_B34_SP6;
  c[52] = cDp0 * B_new / B_old; //BpBm_D2star

  double B_new = ((B_Ddstar41 + B_Ddstar42 * B_Dstarp1) / 2) * t0p * BrB34;
  double B_old = ((B_Ddstar41_SP6 + B_Ddstar42_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6)) * B_B34_SP6;
  c[53] = cDp0 * B_new / B_old; //B0B0_D2star

  double B_new = fp0 * BrB41;
  double B_old = f_p0_SP6 * B_B41_SP6;
  c[54] = cDp0 * B_new / B_old; //BpBm_DPi

  double B_new = t0p * BrB41;
  double B_old = (1 / t_p0_SP6) * B_B41_SP6;
  c[55] = cDp0 * B_new / B_old; //B0B0_DPi

  double B_new = fp0 * (B_Dstarp2 + B_Dstarp3) * BrB42;
  double B_old = f_p0_SP6 * (B_Dstarp2_SP6 + B_Dstarp3_SP6) * B_B42_SP6;
  c[56] = cDp0 * B_new / B_old; //BpBm_DstarPi

  double B_new = (B_Dstarp2 + B_Dstarp3) * t0p * BrB42;
  double B_old = ((B_Dstarp2_SP6 + B_Dstarp3_SP6) / t_p0_SP6) * B_B42_SP6;
  c[57] = cDp0 * B_new / B_old; //B0B0_DstarPi

  double B_new = (B_Dstarp2 + B_Dstarp3) * t0p * BrB2;
  double B_old = ((B_Dstarp2_SP6 + B_Dstarp3_SP6) / t_p0_SP6) * B_B2_SP6;
  c[58] = cDp0 * B_new / B_old; //B0B0_Dstar

  double B_new = t0p * BrB1;
  double B_old = (1 / t_p0_SP6) * B_B1_SP6;
  c[59] = cDp0 * B_new / B_old; //B0B0_Dlnu

  c[60] = 1;  //B0B0_Dlnu_0th
  c[61] = 1;  //B0B0_Dlnu_1st
  c[62] = 1;  //B0B0_Dlnu_2nd


  //###########################################################
  // Make Dlnu array #################
  double num_sub_D0KPi[nbinsx][nbinsy][nbinsz];
  double var_sub_D0KPi[nbinsx][nbinsy][nbinsz];
  double err_sub_D0KPi[nbinsx][nbinsy][nbinsz];
  double num_sub_DchKPiPi[nbinsx][nbinsy][nbinsz];
  double var_sub_DchKPiPi[nbinsx][nbinsy][nbinsz];
  double err_sub_DchKPiPi[nbinsx][nbinsy][nbinsz];
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	num_sub_D0KPi[i][j][k] = 0;
	var_sub_D0KPi[i][j][k] = 0;
	err_sub_D0KPi[i][j][k] = 0;
	num_sub_DchKPiPi[i][j][k] = 0;
	var_sub_DchKPiPi[i][j][k] = 0;
	err_sub_DchKPiPi[i][j][k] = 0;
      }
    }
  }


  int ll0_D0 = 29;  //Dlnu_0th_D0
  int ll1_D0 = 30;  //Dlnu_1st_D0
  int ll2_D0 = 31;  //Dlnu_2nd_D0
  int ll0_Dch = 60;  //Dlnu_0th_Dch
  int ll1_Dch = 61;  //Dlnu_1st_Dch
  int ll2_Dch = 62;  //Dlnu_2nd_Dch
  if (fixSlope != 2) {
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  num_sub_D0KPi[i][j][k] 
	    = num_3D[ll0_D0][i][j][k] 
	    - num_3D[ll1_D0][i][j][k]*rho_D_sq
	    + num_3D[ll2_D0][i][j][k]*rho_D_sq*rho_D_sq;
	  var_sub_D0KPi[i][j][k] 
	    = err_3D[ll0_D0][i][j][k]*err_3D[ll0_D0][i][j][k]
	    - err_3D[ll1_D0][i][j][k]*err_3D[ll1_D0][i][j][k]*rho_D_sq
	    + err_3D[ll2_D0][i][j][k]*err_3D[ll2_D0][i][j][k]*rho_D_sq*rho_D_sq;
	  err_sub_D0KPi[i][j][k] = sqrt(var_sub_D0KPi[i][j][k]);
	  
	  num_sub_DchKPiPi[i][j][k]
	    = num_3D[ll0_Dch][i][j][k] 
	    - num_3D[ll1_Dch][i][j][k]*rho_D_sq
	    + num_3D[ll2_Dch][i][j][k]*rho_D_sq*rho_D_sq;
	  var_sub_DchKPiPi[i][j][k] 
	    = err_3D[ll0_Dch][i][j][k]*err_3D[ll0_Dch][i][j][k]
	    - err_3D[ll1_Dch][i][j][k]*err_3D[ll1_Dch][i][j][k]*rho_D_sq
	    + err_3D[ll2_Dch][i][j][k]*err_3D[ll2_Dch][i][j][k]*rho_D_sq*rho_D_sq;
	  err_sub_DchKPiPi[i][j][k] = sqrt(var_sub_DchKPiPi[i][j][k]);
	  
	  /*	  if (showChisq2) {
		  cout << "num_3D["<<ll0_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		  << num_3D[ll0_D0][i][j][k] << endl;
		  cout << "err_3D["<<ll0_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		  << err_3D[ll0_D0][i][j][k] << endl;
	    cout << "num_3D["<<ll1_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[ll1_D0][i][j][k] << endl;
	    cout << "err_3D["<<ll1_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[ll1_D0][i][j][k] << endl;
	    cout << "num_3D["<<ll2_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[ll2_D0][i][j][k] << endl;
	    cout << "err_3D["<<ll2_D0<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[ll2_D0][i][j][k] << endl;
	    cout << "num_sub_D0KPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_sub_D0KPi[i][j][k] << endl;
	    cout << "var_sub_D0KPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << var_sub_D0KPi[i][j][k] << endl;
	    cout << "err_sub_D0KPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_sub_D0KPi[i][j][k] << endl;
	    
	    cout << "num_3D["<<ll0_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[ll0_Dch][i][j][k] << endl;
	    cout << "err_3D["<<ll0_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[ll0_Dch][i][j][k] << endl;
	    cout << "num_3D["<<ll1_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[ll1_Dch][i][j][k] << endl;
	    cout << "err_3D["<<ll1_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[ll1_Dch][i][j][k] << endl;
	    cout << "num_3D["<<ll2_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[ll2_Dch][i][j][k] << endl;
	    cout << "err_3D["<<ll2_Dch<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[ll2_Dch][i][j][k] << endl;
	    cout << "num_sub_DchKPiPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_sub_DchKPiPi[i][j][k] << endl;
	    cout << "var_sub_DchKPiPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << var_sub_DchKPiPi[i][j][k] << endl;
	    cout << "err_sub_DchKPiPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_sub_DchKPiPi[i][j][k] << endl;
		 } */
	}
      }
    }
  }


  // Get normalization factor ################
  
  double RNM = 1;
  double RNM_err = 0;
  if (fixSlope != 2) {
    int jj(0);
    for (Int_t j=0;j<nBins;j++) {    
      if (rho_D_sq > lowEnd[j] && rho_D_sq <= hiEnd[j+1]) {
	RNM = RN[j];
	RNM_err = RN_err[j];
	jj = j;
	break;
      }
    }
    //cout << "rho_D_sq = " << rho_D_sq << endl;
    //cout << "jj = " << jj << endl;
    //cout << "lowEnd[jj] = " << lowEnd[jj] << endl;
    //cout << "hiEnd[jj+1] = " << hiEnd[jj+1] << endl;
    //cout << "RN[jj] = " << RN[jj] << endl;
    //cout << "RN_err[jj] = " << RN_err[jj] << endl;
    if (printRN) {
      cout << "rho_D_sq = " << rho_D_sq << endl;
      cout << "RNM = " << RNM << " +- " << RNM_err << endl;
    }
  }


  // Calculate chi-square #######################################
  int nBinsChi(0);
  int nBinsChi_D0(0);
  int nBinsChi_Dch(0);
  Double_t chisq(0);
  Double_t chisq_D0(0);
  Double_t chisq_Dch(0);
  Double_t nMC_D0[nbinsx][nbinsy][nbinsz];
  Double_t nMC_Dch[nbinsx][nbinsy][nbinsz];
  Double_t nCand_D0[nbinsx][nbinsy][nbinsz];
  Double_t nCand_Dch[nbinsx][nbinsy][nbinsz];
  Double_t varMC_D0[nbinsx][nbinsy][nbinsz];
  Double_t varMC_Dch[nbinsx][nbinsy][nbinsz];
  Double_t var_D0[nbinsx][nbinsy][nbinsz];
  Double_t var_Dch[nbinsx][nbinsy][nbinsz];
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	nMC_D0[i][j][k] = 0;
	nMC_Dch[i][j][k] = 0;
	nCand_D0[i][j][k] = 0;
	nCand_Dch[i][j][k] = 0;
	varMC_D0[i][j][k] = 0;
	varMC_Dch[i][j][k] = 0;
	var_D0[i][j][k] = 0;
	var_Dch[i][j][k] = 0;
      }
    }
  }

  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	if (showChisq) {
	  cout << "[i][j][k] = ["<<i<<"]["<<j<<"]["<<k<<"]" << endl;
	}

	// D0KPi
	int llOn = 0;  //for OnPeak
	int llDl = 28;  //for Dlnu
	for (int l=2; l<=27; l++) {
	  nMC_D0[i][j][k] += c[l]*num_3D[l][i][j][k];
	  double sigm1 = c[l]*err_3D[l][i][j][k];
	  varMC_D0[i][j][k] += sigm1 * sigm1;
	  /*if (showChisq2) {
	    cout << "c["<<l<<"] = " << c[l] << endl;
	    cout << "num_3D["<<l<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[l][i][j][k] << endl;
	    cout << "err_3D["<<l<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[l][i][j][k] << endl;
	    cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << nMC_D0[i][j][k] << endl;
	    cout << "sigm1 = " << sigm1 << endl;
	    cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << varMC_D0[i][j][k] << endl;
		 } */
	}

	double sigm2 = 0;
	if (fixSlope == 2) {
	  nMC_D0[i][j][k] += c[llDl]*num_3D[llDl][i][j][k];
	  sigm2 = c[llDl]*err_3D[llDl][i][j][k];
	  /*if (showChisq2) {
	    cout << "c["<<llDl<<"] = " << c[llDl] << endl;
	    cout << "num_3D["<<llDl<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[llDl][i][j][k] << endl;
	    cout << "err_3D["<<llDl<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[llDl][i][j][k] << endl;
		 }*/
	} else {
	  nMC_D0[i][j][k] += c[llDl]*num_sub_D0KPi[i][j][k]*RNM;
	  sigm2 = c[llDl]*err_sub_D0KPi[i][j][k]*RNM;
	  /*if (showChisq2) {
	    cout << "c["<<llDl<<"] = " << c[llDl] << endl;
	    cout << "num_sub_D0KPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_sub_D0KPi[i][j][k] << endl;
	    cout << "err_sub_D0KPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_sub_D0KPi[i][j][k] << endl;
		 }*/
	}
	varMC_D0[i][j][k] += sigm2 * sigm2;
	if (showChisq2) {
	  cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << nMC_D0[i][j][k] << endl;
	  cout <<  "sigm2 = " << sigm2 << endl;
	  cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << varMC_D0[i][j][k] << endl;
	}

	// Calculate chi-sq for D0KPi
	if (showChisq2) {
	  cout << "num_3D[" << llOn << "]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << num_3D[llOn][i][j][k] << endl;
	  cout << "err_3D["<<llOn<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << err_3D[llOn][i][j][k] << endl;
	}
	if (nMC_D0[i][j][k] >= nCandMin && varMC_D0[i][j][k]>0) {
	  nBinsChi_D0++;
	  nCand_D0[i][j][k] = num_3D[llOn][i][j][k] - nMC_D0[i][j][k];
	  double var_OnPeak = err_3D[llOn][i][j][k] * err_3D[llOn][i][j][k];
	  var_D0[i][j][k] = var_OnPeak + varMC_D0[i][j][k];
	  double numChisq = nCand_D0[i][j][k] * nCand_D0[i][j][k];
	  double thisChisq = 0;
	  if (var_D0[i][j][k]>0) thisChisq = numChisq / var_D0[i][j][k];
	  chisq_D0 += thisChisq;
	  if (showChisq2) {
	    cout <<  "var_OnPeak = " << var_OnPeak << endl;
	    cout <<  "numChisq = " << numChisq << endl;
	  }
	  if (showChisq) {
	    cout << "nCand_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << nCand_D0[i][j][k] << endl;
	    cout << "var_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << var_D0[i][j][k] << endl;
	    cout << "thisChisq = " << thisChisq << endl;
	    cout << "chisq_D0 = " << chisq_D0 << endl;
	  }
	}

	// DchKPiPi
	int llOn = 32;  //for OnPeak
	int llDl = 59;  //for Dlnu
	for (int l=34; l<=58; l++) {
	  nMC_Dch[i][j][k] += c[l]*num_3D[l][i][j][k];
	  double sigm1 = c[l]*err_3D[l][i][j][k];
	  varMC_Dch[i][j][k] += sigm1 * sigm1;
	  /*if (showChisq2) {
	    cout << "c["<<l<<"] = " << c[l] << endl;
	    cout << "num_3D["<<l<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[l][i][j][k] << endl;
	    cout << "err_3D["<<l<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[l][i][j][k] << endl;
	    cout << "nMC_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << nMC_Dch[i][j][k] << endl;
	    cout << "sigm1 = " << sigm1 << endl;
	    cout << "varMC_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << varMC_Dch[i][j][k] << endl;
		 }*/
	}

	double sigm2 = 0;
	if (fixSlope == 2) {
	  nMC_Dch[i][j][k] += c[llDl]*num_3D[llDl][i][j][k];
	  sigm2 = c[llDl]*err_3D[llDl][i][j][k];
	  /*if (showChisq2) {
	    cout << "c["<<llDl<<"] = " << c[llDl] << endl;
	    cout << "num_3D["<<llDl<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_3D[llDl][i][j][k] << endl;
	    cout << "err_3D["<<llDl<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_3D[llDl][i][j][k] << endl;
		 }*/
	} else {
	  nMC_Dch[i][j][k] += c[llDl]*num_sub_DchKPiPi[i][j][k]*RNM;
	  sigm2 = c[llDl]*err_sub_DchKPiPi[i][j][k]*RNM;
	  /*if (showChisq2) {
	    cout << "c["<<llDl<<"] = " << c[llDl] << endl;
	    cout << "num_sub_DchKPiPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << num_sub_DchKPiPi[i][j][k] << endl;
	    cout << "err_sub_DchKPiPi["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << err_sub_DchKPiPi[i][j][k] << endl;
		 }*/
	}
	varMC_Dch[i][j][k] += sigm2 * sigm2;
	if (showChisq2) {
	  cout << "nMC_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << nMC_Dch[i][j][k] << endl;
	  cout <<  "sigm2 = " << sigm2 << endl;
	  cout << "varMC_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << varMC_Dch[i][j][k] << endl;
	}
	
	// Calculate chi-sq for DchKPiPi
	if (showChisq2) {
	  cout << "num_3D[" << llOn << "]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << num_3D[llOn][i][j][k] << endl;
	  cout << "err_3D["<<llOn<<"]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	       << err_3D[llOn][i][j][k] << endl;
	}
	if (nMC_Dch[i][j][k] >= nCandMin && varMC_Dch[i][j][k] > 0) {
	  nBinsChi_Dch++;
	  nCand_Dch[i][j][k] = num_3D[llOn][i][j][k] - nMC_Dch[i][j][k];
	  double var_OnPeak = err_3D[llOn][i][j][k] * err_3D[llOn][i][j][k];
	  var_Dch[i][j][k] = var_OnPeak + varMC_Dch[i][j][k];
	  double numChisq = nCand_Dch[i][j][k] * nCand_Dch[i][j][k];
	  double thisChisq = 0;
	  if (var_Dch[i][j][k]>0) thisChisq = numChisq / var_Dch[i][j][k];
	  chisq_Dch += thisChisq;
	  if (showChisq2) {
	    cout <<  "var_OnPeak = " << var_OnPeak << endl;
	    cout <<  "numChisq = " << numChisq << endl;
	  }
	  if (showChisq) {
	    cout << "nCand_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << nCand_Dch[i][j][k] << endl;
	    cout << "var_Dch["<<i<<"]["<<j<<"]["<<k<<"] = " 
		 << var_Dch[i][j][k] << endl;
	    cout << "thisChisq = " << thisChisq << endl;
	    cout << "chisq_Dch = " << chisq_Dch << endl;
	  }
	}
      }
    }
  } 

  if (chisqMode == 1) {
    nBinsChi = nBinsChi_D0;
    chisq = chisq_D0;
  } else if (chisqMode == 2) {
    nBinsChi = nBinsChi_Dch;
    chisq = chisq_Dch;
  } else {
    nBinsChi = nBinsChi_D0 + nBinsChi_Dch;
    chisq = chisq_D0 + chisq_Dch;
  }

  chisq = chisq + (fp0 - 1)*(fp0 - 1)/0.001 + (cDp0 - 1)*(cDp0 - 1)/0.001;

  if (showBins) {
    cout << "Number of bins into chisq = " << nBinsChi << endl;
    cout << "                 from D0  = " << nBinsChi_D0 << endl;
    cout << "                 from Dch = " << nBinsChi_Dch << endl;
    showBins = 0;
  }

  if (printChisq) {
    cout << "chisq_D0 = " << chisq_D0 << endl;
    cout << "chisq_Dch = " << chisq_Dch << endl;
    cout << "chisq = " << chisq << endl;
    
  }


  f = chisq;
  ncount++;

}