#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"
TString Slope("Old");  // Slope parametrization, "Old" or "New"
TString corrBF("oldBF"); // "oldBF" or "newBF"

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

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

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

// Chane these to meet your requirements
int requireD0Only = 0;  // D0 only
int requireDchOnly = 0;  // D+ only
int printChisq = 0;  // printout chi-sq values
int doMinos = 0;  // do MINOS
int fixSlope = 1; // fix FF slope and use single histogram for Dlnu
int fixDlnu = 0;
double Br_Dlnu = 0.0210;
int fixDstar = 0;
double Br_Dstar = 0.0560;
int fixD0star = 1;
double ratio_D0star = 0.35714; //=20/56  // Br_D0star = ratio_D0star*Br_D1
int fixD1 = 0;
double Br_D1 = 0.0056;
int fixD1p = 1;
double ratio_D1p = 0.66071; //=37/56  // Br_D1p = ratio_D1p*Br_D1
int fixD2star = 1;
double ratio_D2star = 0.66071; //=37/56  // Br_D2star = ratio_D2star*Br_D1
int fixDPi = 0;
double Br_DPi = 0.0060;
int fixDstarPi = 1;
double ratio_DstarPi = 0.33333; //=20/60 //=67/56  // Br_DstarPi = ratio_DstarPi*Br_DPi
int fixBkgd = 0; // fix background coeeficients to 1
int fixFp0 = 1;  // fix f_+0 = 1
int fixCD0p = 1;  // fix c_D_0+ = 1
int nCandMin = 5;  // minimum number of events in one bin
double inclBF = 0.100;  //inclusive b->clnu BF

// Declear common variables
static ncount = 0;
double num_3D[nFiles][nbinsx][nbinsy][nbinsz];
double err_3D[nFiles][nbinsx][nbinsy][nbinsz];
double totN_org[nFiles];
int showBins = 1;  // print number of bins used in the fit
int showTotNum = 1;  //print total number of candidates


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

void fit_3D_pdf() {

  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);
	}
      }
    }
  }

  // Get the total number of candidates ##############
  for (int l=0; l<nFiles; l++) {
    totN_org[l] = 0;
  }
  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++) {
	  totN_org[l] += c[l]*num_3D[l][i][j][k];
	}
      }
    }
  }


  //######## Fitting ######################
  //#######################################
  
  static Int_t nprm[16] = {0, 1, 
			   2, 3, 4, 5, 
			   6, 7, 
			   8, 
			   9, 10, 11, 12, 
			   13, 14, 15};
  static Double_t vstrt[16] = {0.1, 0.44,
			       0.007, 0.036, 0.023, 0.017,
			       0.033, 0.016,
			       0,
			       0.1, 0.1, 0.02, 0.05,
			       1, 1, 1};
  /* static Double_t vstrt[16] = {0.4859, 0.05569,
			       0.0020, 0.006518, 0.0037, 0.0037,
			       0.000934, 0.0020,
			       2.342,
			       0.87594, 1.09287, 0, 1.04067,
			       1.02439, 1, 0.94820};  */
  static Double_t stp[16] = {0.01, 0.01, 
			     0.01, 0.01, 0.01, 0.01, 
			     0.01, 0.01,
			     0.01,
			     0.01, 0.01, 0.01, 0.01,
			     0.01, 0.01, 0.01};
  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;
  Int_t ierflg = 0;

  TStopwatch timer;
   
  //initialize TMinuit with a maximum of 16 params
  TMinuit *gMinuit = new TMinuit(16);
  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,5,ierflg);
  gMinuit->mnparm(nprm[9], "Uncorr", vstrt[9], stp[9], 0,1,ierflg);
  gMinuit->mnparm(nprm[10], "CaslcL", vstrt[10], stp[10], 0,1,ierflg);
  gMinuit->mnparm(nprm[11], "FakeL", vstrt[11], stp[11], 0,1,ierflg);
  gMinuit->mnparm(nprm[12], "ContD", vstrt[12], stp[12], 0,1,ierflg);
  gMinuit->mnparm(nprm[13], "f_+0", vstrt[13], stp[13], 0,2,ierflg);
  //gMinuit->mnparm(nprm[14], "t_0+", vstrt[14], stp[14], 0,2,ierflg);
  gMinuit->mnparm(nprm[15], "c_D+0", vstrt[15], stp[15], 0,2,ierflg);
  if (ierflg) {
    Printf(" UNABLE TO DEFINE PARAMETER NO.");
    return ierflg;
  }
  
  // Perform fitting ############
  gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg);

  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);
  }
  if (fixFp0) gMinuit->mnexcm("FIX", &p14,1,ierflg);
  if (fixCD0p) gMinuit->mnexcm("FIX", &p16,1,ierflg);

  gMinuit->mnexcm("MIGRAD", &p0 ,0,ierflg);
  if (doMinos) gMinuit->mnexcm("MINOS",  &p0 ,0,ierflg);
  gMinuit->mnexcm("SHOW COVARIANCE", &p0 ,0,ierflg);

  gMinuit->mnexcm("CALL FCN", &p2 , 1,ierflg);

  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 Bkg1(0), Bkg1_err(0);
  Double_t Bkg2(0), Bkg2_err(0);
  Double_t Bkg3(0), Bkg3_err(0);
  Double_t Bkg4(0), Bkg4_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, Bkg1, Bkg1_err);
  gMinuit->GetParameter(10, Bkg2, Bkg2_err);
  gMinuit->GetParameter(11, Bkg3, Bkg3_err);
  gMinuit->GetParameter(12, Bkg4, Bkg4_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 << "Bkg1 = " << Bkg1 << endl;
  cout << "Bkg2 = " << Bkg2 << endl;
  cout << "Bkg3 = " << Bkg3 << endl;
  cout << "Bkg4 = " << Bkg4 << 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 << "Bkg1_err = " << Bkg1_err << endl;
  cout << "Bkg2_err = " << Bkg2_err << endl;
  cout << "Bkg3_err = " << Bkg3_err << endl;
  cout << "Bkg4_err = " << Bkg4_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 * BF32;
  if (fixDstarPi)  BF42 = ratio_DstarPi * BF32;
  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 << "Bkg1 = " << Bkg1 << endl;
  cout << "Bkg2 = " << Bkg2 << endl;
  cout << "Bkg3 = " << Bkg3 << endl;
  cout << "Bkg4 = " << Bkg4 << endl;

  // Calcurate BFs
  double TotFrac 
    = BF1 + BF2 + BF31 + BF32 + BF33 + BF34 + BF41 + BF42 
    + Bkg1 + Bkg2 + Bkg3 + Bkg4;
  double TotBkgd = Bkg1 + Bkg2 + Bkg3 + Bkg4;
  double denom = TotFrac * (1 - TotBkgd);
  cout << "TotFrac = " << TotFrac << endl;
  cout << "TotBkgd = " << TotBkgd << endl;
  cout << "inclBF = " << inclBF << 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);
  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 << "BF_B1 = " << BF_B1 << " +- " << err_B1 
       << " (" << percent << " %)" << endl;
  BF_B2 = BF2 * inclBF / denom;
  err_B2 = BF2_err * inclBF / denom;
  double percent = BF2_err * 100 / BF2;
  cout << "BF_B2 = " << BF_B2 << " +- " << err_B2 
       << " (" << percent << " %)" << endl;
  BF_B31 = BF31 * inclBF / denom;
  err_B31 = BF31_err * inclBF / denom;
  double percent = BF31_err * 100 / BF31;
  cout << "BF_B31 = " << BF_B31 << " +- " << err_B31 
       << " (" << percent << " %)" << endl;
  BF_B32 = BF32 * inclBF / denom;
  err_B32 = BF32_err * inclBF / denom;
  double percent = BF32_err * 100 / BF32;
  cout << "BF_B32 = " << BF_B32 << " +- " << err_B32 
       << " (" << percent << " %)" << endl;
  BF_B33 = BF33 * inclBF / denom;
  err_B33 = BF33_err * inclBF / denom;
  double percent = BF33_err * 100 / BF33;
  cout << "BF_B33 = " << BF_B33 << " +- " << err_B33 
       << " (" << percent << " %)" << endl;
  BF_B34 = BF34 * inclBF / denom;
  err_B34 = BF34_err * inclBF / denom;
  double percent = BF34_err * 100 / BF34;
  cout << "BF_B34 = " << BF_B34 << " +- " << err_B34 
       << " (" << percent << " %)" << endl;
  BF_B41 = BF41 * inclBF / denom;
  err_B41 = BF41_err * inclBF / denom;
  double percent = BF41_err * 100 / BF41;
  cout << "BF_B41 = " << BF_B41 << " +- " << err_B41 
       << " (" << percent << " %)" << endl;
  BF_B42 = BF42 * inclBF / denom;
  err_B42 = BF42_err * inclBF / denom;
  double percent = BF42_err * 100 / BF42;
  cout << "BF_B42 = " << BF_B42 << " +- " << err_B42 
       << " (" << percent << " %)" << endl;
  double percent(0);
  if (slope>0) percent = slopeErr * 100 / slope;
  cout << "slope = " << slope << " +- " << slopeErr 
       << " (" << percent << " %)" << endl;
  
}



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

void fcnBToD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag)
{
  enum { nPdf = 13 };

  // Branching fractions
  double BF[nPdf];
  BF[0] = 1;

  if (fixDlnu)  BF[5] = Br_Dlnu;
  else BF[5] = x[0];
  if (fixDstar)   BF[6] = Br_Dstar;
  else BF[6] = x[1];
  if (fixD0star) BF[7] = x[3]*ratio_D0star;
  else BF[7] = x[2];
  if (fixD1)   BF[8] = Br_D1;
  else BF[8] = x[3];
  if (fixD1p)   BF[9] = x[3]*ratio_D1p;
  else BF[9] = x[4];
  if (fixD2star)  BF[10] = x[3]*ratio_D2star;
  else BF[10] = x[5];
  if (fixDPi)  BF[11] = Br_DPi;
  else BF[11] = x[6];
  if (fixDstarPi)  BF[12] = x[6]*ratio_DstarPi;
  else BF[12] = x[7];

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

  // Peaking background
  BF[1] = x[9];
  BF[2] = x[10];
  BF[3] = x[11];
  BF[4] = x[12];
  
  // ratio f_+- / f_00 and t_+- / t_00
  Double_t fp0, t0p;
  fp0 = x[13];
  //t0p = x[14];
  t0p = 1/t_p0;

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

  // Coefficients
  double c[nFiles];
  c[0] = 1;  //OnPeak_D0KPi
  c[1] = 1;  //OffPeak_D0KPi

  c[2] = 1; //BpBm_Comb_D0KPi
  c[3] = 1; //B0B0_Comb_D0KPi
  c[4] = 1; //ccbar_Cont_D0KPi
  c[5] = 1; //uds_Cont_D0KPi
  c[6] = 1; //tautau_Cont_D0KPi

  c[7] = fp0/f_p0_SP6; //BpBm_DiffB_D0KPi
  c[8] = 1;            //B0B0_DiffB_D0KPi
  c[9] = fp0/f_p0_SP6; //BpBm_CascL_D0KPi
  c[10] = 1;            //B0B0_CascL_D0KPi
  c[11] = fp0/f_p0_SP6; //BpBm_LMisID_D0KPi
  c[12] = 1;            //B0B0_LMisID_D0KPi
  c[13] = 1;            //ccbar_RealD_D0KPi

  double B_new = fp0;
  double B_old = f_p0_SP6;
  c[14] = B_new / B_old; //BpBm_Dlnu_D0KPi

  double B_new = B_Dstarp1 * t0p;
  double B_old = (B_Dstarp1_SP6 / t_p0_SP6);
  c[15] = B_new / B_old; //BpBm_Dstar_D0KPi
  
  double B_new = fp0;
  double B_old = f_p0_SP6;
  c[16] = B_new / B_old; //B0B0_Dstar_D0KPi

  double B_new = fp0*(B_Ddstar11 / 2 + B_Ddstar12 * (1/2 + B_Dstarp1));
  double B_old 
    = f_p0_SP6*(B_Ddstar11_SP6 / 2 + B_Ddstar12_SP6 * (1/2 + B_Dstarp1_SP6));
  c[17] = B_new / B_old; //BpBm_D0star_D0KPi

  double B_new = t0p*(B_Ddstar11 + B_Ddstar12 * (1 + B_Dstarp1 / 2));
  double B_old = (B_Ddstar11_SP6 + B_Ddstar12_SP6 * (1 + B_Dstarp1_SP6 / 2));
  c[18] = B_new / B_old; //B0B0_D0star_D0KPi

  double B_new = fp0*(B_Ddstar21 / 2 + B_Ddstar22 * (1/2 + B_Dstarp1));
  double B_old 
    = f_p0_SP6*(B_Ddstar21_SP6 / 2 + B_Ddstar22_SP6 * (1/2 + B_Dstarp1_SP6));
  c[19] = B_new / B_old; //BpBm_D1_D0KPi

  double B_new = t0p*(B_Ddstar21 + B_Ddstar22 * (1 + B_Dstarp1 / 2));
  double B_old = (B_Ddstar21_SP6 + B_Ddstar22_SP6 * (1 + B_Dstarp1_SP6 / 2));
  c[20] = B_new / B_old; //B0B0_D1_D0KPi

  double B_new = fp0*(B_Ddstar31 / 2 + B_Ddstar32 * (1/2 + B_Dstarp1));
  double B_old 
    = f_p0_SP6*(B_Ddstar31_SP6 / 2 + B_Ddstar32_SP6 * (1/2 + B_Dstarp1_SP6));
  c[21] = B_new / B_old; //BpBm_D1prime_D0KPi

  double B_new = t0p*(B_Ddstar31 + B_Ddstar32 * (1 + B_Dstarp1 / 2));
  double B_old = (B_Ddstar31_SP6 + B_Ddstar32_SP6 * (1 + B_Dstarp1_SP6 / 2));
  c[22] = B_new / B_old; //B0B0_D1prime_D0KPi

  double B_new = fp0*(B_Ddstar41 / 2 + B_Ddstar42 * (1/2 + B_Dstarp1));
  double B_old 
    = f_p0_SP6*(B_Ddstar41_SP6 / 2 + B_Ddstar42_SP6 * (1/2 + B_Dstarp1_SP6));
  c[23] = B_new / B_old; //BpBm_D2star_D0KPi

  double B_new = t0p*(B_Ddstar41 + B_Ddstar42 * (1 + B_Dstarp1 / 2));
  double B_old = (B_Ddstar41_SP6 + B_Ddstar42_SP6 * (1 + B_Dstarp1_SP6 / 2));
  c[24] = B_new / B_old; //B0B0_D2star_D0KPi

  double B_new = (fp0 / 2);
  double B_old = (f_p0_SP6 / 2);
  c[25] = B_new / B_old; //BpBm_DPi_D0KPi

  double B_new = t0p;
  double B_old = (1 / t_p0_SP6);
  c[26] = B_new / B_old; //B0B0_DPi_D0KPi

  double B_new = fp0 * (1/2 + B_Dstarp1);
  double B_old = f_p0_SP6 * (1/2 + B_Dstarp1_SP6);
  c[27] = B_new / B_old; //BpBm_DstarPi_D0KPi

  double B_new = (1 + B_Dstarp1 / 2) * t0p;
  double B_old = ((1 + B_Dstarp1_SP6 / 2) / t_p0_SP6);
  c[28] = B_new / B_old; //B0B0_DstarPi_D0KPi


  c[29] = 1; //OnPeak_DchKPiPi
  c[30] = 1; //OffPeak_DchKPiPi

  c[31] = cDp0 * 1; //BpBm_Comb_DchKPiPi
  c[32] = cDp0 * 1; //B0B0_Comb_DchKPiPi
  c[33] = cDp0 * 1; //ccbar_Cont_DchKPiPi
  c[34] = cDp0 * 1; //uds_Cont_DchKPiP
  c[35] = cDp0 * 1; //tautau_Cont_DchKPiPi

  c[36] = cDp0 * fp0/f_p0_SP6; //BpBm_DiffB_DchKPiPi
  c[37] = cDp0;                //B0B0_DiffB_DchKPiPi    
  c[38] = cDp0 * fp0/f_p0_SP6; //BpBm_CascL_DchKPiPi
  c[39] = cDp0;                //B0B0_CascL_DchKPiPi
  c[40] = cDp0 * fp0/f_p0_SP6; //BpBm_LMisID_DchKPiP
  c[41] = cDp0;                //B0B0_LMisID_DchKPiPi
  c[42] = cDp0;                //ccbar_RealD_DchKPiPi

  double B_new = t0p;
  double B_old = (1 / t_p0_SP6);
  c[43] = cDp0 * B_new / B_old; //B0B0_Dlnu_DchKPiPi

  double B_new = (B_Dstarp2 + B_Dstarp3) * t0p;
  double B_old = ((B_Dstarp2_SP6 + B_Dstarp3_SP6) / t_p0_SP6);
  c[44] = cDp0 * B_new / B_old; //B0B0_Dstar_DchKPiPi

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

  double B_new = ((B_Ddstar11 + B_Ddstar12 * B_Dstarp1) / 2) * t0p;
  double B_old 
    = ((B_Ddstar11_SP6 + B_Ddstar12_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6));
  c[46] = cDp0 * B_new / B_old;  //B0B0_D0star_DchKPiPi

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

  double B_new = ((B_Ddstar21 + B_Ddstar22 * B_Dstarp1) / 2) * t0p;
  double B_old 
    = ((B_Ddstar21_SP6 + B_Ddstar22_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6));
  c[48] = cDp0 * B_new / B_old; //B0B0_D1_DchKPiPi

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

  double B_new = ((B_Ddstar31 + B_Ddstar32 * B_Dstarp1) / 2) * t0p;
  double B_old 
    = ((B_Ddstar31_SP6 + B_Ddstar32_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6));
  c[50] = cDp0 * B_new / B_old; //B0B0_D1prime_DchKPiPi

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

  double B_new = ((B_Ddstar41 + B_Ddstar42 * B_Dstarp1) / 2) * t0p;
  double B_old 
    = ((B_Ddstar41_SP6 + B_Ddstar42_SP6 * B_Dstarp1_SP6) / (2 * t_p0_SP6));
  c[52] = cDp0 * B_new / B_old; //B0B0_D2star_DchKPiPi

  double B_new = fp0;
  double B_old = f_p0_SP6;
  c[53] = cDp0 * B_new / B_old; //BpBm_DPi_DchKPiPi

  double B_new = t0p;
  double B_old = (1 / t_p0_SP6);
  c[54] = cDp0 * B_new / B_old; //B0B0_DPi_DchKPiPi

  double B_new = fp0 * (B_Dstarp2 + B_Dstarp3);
  double B_old = f_p0_SP6 * (B_Dstarp2_SP6 + B_Dstarp3_SP6);
  c[55] = cDp0 * B_new / B_old; //BpBm_DstarPi_DchKPiPi

  double B_new = (B_Dstarp2 + B_Dstarp3) * t0p;
  double B_old = ((B_Dstarp2_SP6 + B_Dstarp3_SP6) / t_p0_SP6);
  c[56] = cDp0 * B_new / B_old; //B0B0_DstarPi_DchKPiPi


  c[57]=c[58]=c[59]=1;
  c[60]=c[61]=c[62]=1;

  // Get the total number of candidates ##############
  double totN[nFiles];
  for (int l=0; l<nFiles; l++) {
    totN[l] = 0;
  }
  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++) {
	  totN[l] += c[l]*num_3D[l][i][j][k];
	}
      }
    }
  }

  // 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];

  double totN_Dlnu_D0(0);
  double totN_Dlnu_Dch(0);
  if (!fixSlope) {
    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[57][i][j][k] - num_3D[58][i][j][k]*rho_D_sq
	    + num_3D[59][i][j][k]*rho_D_sq*rho_D_sq;
	  var_sub_D0KPi[i][j][k] 
	    = err_3D[57][i][j][k]*err_3D[57][i][j][k]
	    - err_3D[58][i][j][k]*err_3D[58][i][j][k]*rho_D_sq
	    + err_3D[59][i][j][k]*err_3D[59][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[60][i][j][k] - num_3D[61][i][j][k]*rho_D_sq
	    + num_3D[62][i][j][k]*rho_D_sq*rho_D_sq;
	  var_sub_DchKPiPi[i][j][k] 
	    = err_3D[60][i][j][k]*err_3D[60][i][j][k]
	    - err_3D[61][i][j][k]*err_3D[61][i][j][k]*rho_D_sq
	    + err_3D[62][i][j][k]*err_3D[62][i][j][k]*rho_D_sq*rho_D_sq;
	  err_sub_DchKPiPi[i][j][k] = sqrt(var_sub_DchKPiPi[i][j][k]);

	  //cout << "num_3D[57]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << num_3D[57][i][j][k] << endl;
	  //cout << "err_3D[57]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << err_3D[57][i][j][k] << endl;
	  //cout << "num_3D[58]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << num_3D[58][i][j][k] << endl;
	  //cout << "err_3D[58]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << err_3D[58][i][j][k] << endl;
	  //cout << "num_3D[59]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << num_3D[59][i][j][k] << endl;
	  //cout << "err_3D[59]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << err_3D[59][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;
	}
      }
    }

    // Get total numbers
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  totN_Dlnu_D0 += c[14]*num_sub_D0KPi[i][j][k];
	  totN_Dlnu_Dch += c[43]*num_sub_DchKPiPi[i][j][k];
	}
      }
    }

  }

  if (showTotNum) {
    for (int l=0; l<nFiles; l++) {
      cout << "totN[" << l << "] = " << totN[l] << endl;
    }
    cout << "totN_Dlnu_D0 = " << totN_Dlnu_D0 << endl;
    cout << "totN_Dlnu_Dch = " << totN_Dlnu_Dch << endl;
    double totN_MC_D0(0), totN_MC_Dch(0);
    for (int l=7; l<=28; l++) {
      totN_MC_D0 += totN[l];
    }
    for (int l=36; l<=56; l++) {
      totN_MC_Dch += totN[l];
    }
    cout << "totN_MC_D0 = " << totN_MC_D0 << endl;
    cout << "totN_MC_Dch = " << totN_MC_Dch << endl;

    cout << "OnPeak_D0 = " << totN[0] << endl;
    cout << "DiffB_D0 = " << totN[7] + totN[8] 
	 << " (" << (totN[7] + totN[8])*100/totN_MC_D0 << " %)" << endl;
    cout << "CascL_D0 = " << totN[9] + totN[10] 
	 << " (" << (totN[9] + totN[10])*100/totN_MC_D0 << " %)" << endl;
    cout << "LMisID_D0 = " << totN[11] + totN[12] 
	 << " (" << (totN[11] + totN[12])*100/totN_MC_D0 << " %)" << endl;
    cout << "realD_D0 = " << totN[13]  
	 << " (" << totN[13]*100/totN_MC_D0 << " %)" << endl;
    if (fixSlope) {
      cout << "Dlnu_D0 = " << totN[14]  
	   << " (" << totN[14]*100/totN_MC_D0 << " %)" << endl;
    } else {
      cout << "Dlnu_D0 = " << totN_Dlnu_D0  
	   << " (" << totN_Dlnu_D0*100/totN_MC_D0 << " %)" << endl;
    }
    cout << "Dstar_D0 = " << totN[15] + totN[16] 
	 << " (" << (totN[15] + totN[16])*100/totN_MC_D0 << " %)" << endl;
    cout << "D0star_D0 = " << totN[17] + totN[18] 
	 << " (" << (totN[17] + totN[18])*100/totN_MC_D0 << " %)" << endl;
    cout << "D1_D0 = " << totN[19] + totN[20] 
	 << " (" << (totN[19] + totN[20])*100/totN_MC_D0 << " %)" << endl;
    cout << "D1prime_D0 = " << totN[21] + totN[22] 
	 << " (" << (totN[21] + totN[22])*100/totN_MC_D0 << " %)" << endl;
    cout << "D2star_D0 = " << totN[23] + totN[24] 
	 << " (" << (totN[23] + totN[24])*100/totN_MC_D0 << " %)" << endl;
    cout << "DPi_D0 = " << totN[25] + totN[26] 
	 << " (" << (totN[25] + totN[26])*100/totN_MC_D0 << " %)" << endl;
    cout << "DstarPi_D0 = " << totN[27] + totN[28] 
	 << " (" << (totN[27] + totN[28])*100/totN_MC_D0 << " %)" << endl;
    cout << "OnPeak_Dch = " << totN[29] << endl;
    cout << "DiffB_Dch = " << totN[36] + totN[37] 
	 << " (" << (totN[36] + totN[37])*100/totN_MC_Dch << " %)" << endl;
    cout << "CascL_Dch = " << totN[38] + totN[39] 
	 << " (" << (totN[38] + totN[39])*100/totN_MC_Dch << " %)" << endl;
    cout << "LMisID_Dch = " << totN[40] + totN[41] 
	 << " (" << (totN[40] + totN[41])*100/totN_MC_Dch << " %)" << endl;
    cout << "realD_Dch = " << totN[42]  
	 << " (" << totN[42]*100/totN_MC_Dch << " %)" << endl;
    cout << "Dlnu_org_Dch = " << totN[43]  
	 << " (" << totN[43]*100/totN_MC_Dch << " %)" << endl;
    cout << "Dlnu_new_Dch = " << totN_Dlnu_Dch  
	 << " (" << totN_Dlnu_Dch*100/totN_MC_Dch << " %)" << endl;
    cout << "Dstar_Dch = " << totN[44]
	 << " (" << totN[44]*100/totN_MC_Dch << " %)" << endl;
    cout << "D0star_Dch = " << totN[45] + totN[46] 
	 << " (" << (totN[45] + totN[46])*100/totN_MC_Dch << " %)" << endl;
    cout << "D1_Dch = " << totN[47] + totN[48] 
	 << " (" << (totN[47] + totN[48])*100/totN_MC_Dch << " %)" << endl;
    cout << "D1prime_Dch = " << totN[49] + totN[50] 
	 << " (" << (totN[49] + totN[50])*100/totN_MC_Dch << " %)" << endl;
    cout << "D2star_Dch = " << totN[51] + totN[52] 
	 << " (" << (totN[51] + totN[52])*100/totN_MC_Dch << " %)" << endl;
    cout << "DPi_Dch = " << totN[53] + totN[54] 
	 << " (" << (totN[53] + totN[54])*100/totN_MC_Dch << " %)" << endl;
    cout << "DstarPi_Dch = " << totN[55] + totN[56] 
	 << " (" << (totN[55] + totN[56])*100/totN_MC_Dch << " %)" << endl;
    showTotNum = 0;
  }


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

  // Make PDF
  Double_t PDF_D0[nPdf][nbinsx][nbinsy][nbinsz];
  Double_t var_D0[nPdf][nbinsx][nbinsy][nbinsz];
  Double_t PDF_Dch[nPdf][nbinsx][nbinsy][nbinsz];
  Double_t var_Dch[nPdf][nbinsx][nbinsy][nbinsz];

  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {

	// Data[0]
	double num = c[0]*num_3D[0][i][j][k];
	double den = totN[0];
	double num2 = c[0]*err_3D[0][i][j][k]*c[0]*err_3D[0][i][j][k];
	for (int l=2; l<=6; l++) {
	  num -= c[l]*num_3D[l][i][j][k];
	  den -= totN[l];
	  num2 += c[l]*err_3D[l][i][j][k]*c[l]*err_3D[l][i][j][k];
	}
	PDF_D0[0][i][j][k] = num/den;
	double den2 = den*den;
	var_D0[0][i][j][k] = num2/den2;

	// DiffB[1], CascL[2], FakeL[3]
	for (int p=1; p<=3; p++) {
	  int p1 = 5 + 2*p;
	  int p2 = 6 + 2*p;
	  double num = c[p1]*num_3D[p1][i][j][k] + c[p2]*num_3D[p2][i][j][k];
	  double den = totN[p1] + totN[p2];
	  PDF_D0[p][i][j][k] = num/den;
	  double num2 
	    = c[p1]*err_3D[p1][i][j][k]*c[p1]*err_3D[p1][i][j][k]
	    + c[p2]*err_3D[p2][i][j][k]*c[p2]*err_3D[p2][i][j][k];
	  double den2 = den*den;
	  var_D0[p][i][j][k] = num2/den2;
	}

	// RealD[4]
	double num = c[13]*num_3D[13][i][j][k];
	double den = totN[13];
	PDF_D0[4][i][j][k] = num/den;
	double num2 
	  = c[13]*err_3D[13][i][j][k]*c[13]*err_3D[13][i][j][k];
	double den2 = den*den;
	var_D0[4][i][j][k] = num2/den2;

	// Dlnu[5]
	if (fixSlope) {
	  double num = c[14]*num_3D[14][i][j][k];
	  double den = totN[14];
	  PDF_D0[5][i][j][k] = num/den;
	  double num2 
	    = c[14]*err_3D[14][i][j][k]*c[14]*err_3D[14][i][j][k];
	  double den2 = den*den;
	  var_D0[5][i][j][k] = num2/den2;
	} else {
	  double num = c[14]*num_sub_D0KPi[i][j][k];
	  double den = totN_Dlnu_D0;
	  PDF_D0[5][i][j][k] = num/den;
	  double num2
	    = c[14]*err_sub_D0KPi[i][j][k]*c[14]*err_sub_D0KPi[i][j][k];
	  double den2 = den*den;
	  var_D0[5][i][j][k] = num2/den2;
	}

	// Dstar[6], D0star[7], D1[8], D1prime[9], D2star[10], 
	// DPi[11],DstarPi[12] 
	for (int p=6; p<=12; p++) {
	  int p1 = 3 + 2*p;
	  int p2 = 4 + 2*p;
	  double num = c[p1]*num_3D[p1][i][j][k] + c[p2]*num_3D[p2][i][j][k];
	  double den = totN[p1] + totN[p2];
	  PDF_D0[p][i][j][k] = num/den;
	  double num2 
	    = c[p1]*err_3D[p1][i][j][k]*c[p1]*err_3D[p1][i][j][k]
	    + c[p2]*err_3D[p2][i][j][k]*c[p2]*err_3D[p2][i][j][k];
	  double den2 = den*den;
	  var_D0[p][i][j][k] = num2/den2;
	}
      }
    }
  } 

  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {

	// Data[0]
	double num = c[29]*num_3D[29][i][j][k]; 
	double den = totN[29];
	double num2
	  = c[29]*err_3D[29][i][j][k]*c[29]*err_3D[29][i][j][k];
	for (int l=31; l<=35; l++) {
	  num -= c[l]*num_3D[l][i][j][k];
	  den -= totN[l];
	  num2 += c[l]*err_3D[l][i][j][k]*c[l]*err_3D[l][i][j][k];
	}
	PDF_Dch[0][i][j][k] = num/den;
	double den2 = den*den;
	var_Dch[0][i][j][k] = num2/den2;

	// DiffB[1], CascL[2], FakeL[3]
	for (int p=1; p<=3; p++) {
	  int p1 = 34 + 2*p;
	  int p2 = 35 + 2*p;
	  double num = c[p1]*num_3D[p1][i][j][k] + c[p2]*num_3D[p2][i][j][k];
	  double den = totN[p1] + totN[p2];
	  PDF_Dch[p][i][j][k] = num/den;
	  double num2 
	    = c[p1]*err_3D[p1][i][j][k]*c[p1]*err_3D[p1][i][j][k]
	    + c[p2]*err_3D[p2][i][j][k]*c[p2]*err_3D[p2][i][j][k];
	  double den2 = den*den;
	  var_Dch[p][i][j][k] = num2/den2;
	}

	// RealD[4]
	double num = c[42]*num_3D[42][i][j][k];
	double den = totN[42];
	PDF_Dch[4][i][j][k] = num/den;
	double num2 
	  = c[42]*err_3D[42][i][j][k]*c[42]*err_3D[42][i][j][k];
	double den2 = den*den;
	var_Dch[4][i][j][k] = num2/den2;

	// Dlnu[5]
	if (fixSlope) {
	  double num = c[43]*num_3D[43][i][j][k];
	  double den = totN[43];
	  PDF_Dch[5][i][j][k] = num/den;
	  double num2 
	    = c[43]*err_3D[43][i][j][k]*c[43]*err_3D[43][i][j][k];
	  double den2 = den*den;
	  var_Dch[5][i][j][k] = num2/den2;
	} else {
	  double num = c[43]*num_sub_DchKPiPi[i][j][k];
	  double den = totN_Dlnu_Dch;
	  PDF_Dch[5][i][j][k] = num/den;
	  double num2
	    = c[43]*err_sub_DchKPiPi[i][j][k]*c[43]*err_sub_DchKPiPi[i][j][k];
	  double den2 = den*den;
	  var_Dch[5][i][j][k] = num2/den2;
	}

	// Dstar[6]
	double num = c[44]*num_3D[44][i][j][k];
	double den = totN[44];
	PDF_Dch[6][i][j][k] = num/den;
	double num2 
	  = c[44]*err_3D[44][i][j][k]*c[44]*err_3D[44][i][j][k];
	double den2 = den*den;
	var_Dch[6][i][j][k] = num2/den2;

	// D0star[7], D1[8], D1prime[9], D2star[10], 
	// DPi[11],DstarPi[12] 
	for (int p=7; p<=12; p++) {
	  int p1 = 31 + 2*p;
	  int p2 = 32 + 2*p;
	  double num = c[p1]*num_3D[p1][i][j][k] + c[p2]*num_3D[p2][i][j][k];
	  double den = totN[p1] + totN[p2];
	  PDF_Dch[p][i][j][k] = num/den;
	  double num2 
	    = c[p1]*err_3D[p1][i][j][k]*c[p1]*err_3D[p1][i][j][k]
	    + c[p2]*err_3D[p2][i][j][k]*c[p2]*err_3D[p2][i][j][k];
	  double den2 = den*den;
	  var_Dch[p][i][j][k] = num2/den2;
	}
      }
    }
  } 

  // 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);
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	if (num_3D[0][i][j][k] >= nCandMin) {
	  double pMC_tot(0);
	  double varMC_tot(0);
	  for (int l=1; l<nPdf; l++) {
	    pMC_tot += BF[l]*PDF_D0[l][i][j][k];
	    varMC_tot += BF[l]*BF[l]*var_D0[l][i][j][k];
	  }
	  double num = PDF_D0[0][i][j][k] - pMC_tot;
	  double num2 = num * num;
	  double den = var_D0[0][i][j][k] + varMC_tot;
	  if (den>0) {
	    nBinsChi_D0++;
	    chisq_D0 += num2/den;
	  }
	  //cout << "chisq_D0 = " << chisq_D0 << endl;
	}
      }
    } 
  }
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	if (num_3D[29][i][j][k] >= nCandMin) {
	  double pMC_tot(0);
	  double varMC_tot(0);
	  for (int l=1; l<nPdf; l++) {
	    pMC_tot += BF[l]*PDF_Dch[l][i][j][k];
	    varMC_tot += BF[l]*BF[l]*var_Dch[l][i][j][k];
	  }
	  double num = PDF_Dch[0][i][j][k] - pMC_tot;
	  double num2 = num * num;
	  double den = var_Dch[0][i][j][k] + varMC_tot;
	  if (den>0) {
	    nBinsChi_Dch++;
	    chisq_Dch += num2/den;
	  }
	  //cout << "chisq_Dch = " << chisq_Dch << endl;
	}
      }
    } 
  }

  if (requireD0Only) {
    nBinsChi = nBinsChi_D0;
    chisq = chisq_D0;
  } else if (requireDchOnly) {
    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 D+ = " << 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++;

}