#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("effCorr"); // "noCorr" or "effCorr"
TString corrDstarFF("newFF"); //"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("newBF"); // D->Xenu BF correction, "oldBF" or "newBF"
TString histDir("histogramsFit"); //directory of histograms, "histogramsFit"

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

// Setup file including parameters
#include "fit_3D_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];
int N_old = 0;
int showChisq = 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

int N = 100000;  // for MC integration.
int requireD0Only = 0;  // D0 only
int requireDchOnly = 0;  // D+ only
int doMigrad = 1;  // do MIGRAD
int doMinos = 1;  // do MINOS
int printInfo = 1; // printout chisq and RN when "CALL FCN"
int printDetailInfo = 0; // printout details of chisq when "MIGRAD"
int printInfoSBS = 0; // printout chisq and RN when "MIGRAD"
int fixSlope = 0; // 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 = 0;
double ratio_D0star = 0.35714; //=20/56  Br_D0star = ratio_D0star*Br_D1
int fixD1 = 0;
double Br_D1 = 0.0056;
int fixD1p = 0;
double ratio_D1p = 0.66071; //=37/56  Br_D1p = ratio_D1p*Br_D1
int fixD2star = 0;
double ratio_D2star = 0.66071; //=37/56  Br_D2star = ratio_D2star*Br_D1
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 to 1
int fixFp0 = 0;  // fix f_+0 = 1
int fixCD0p = 0;  // fix c_D_0+ = 1
int fixT0p = 0;  // fix t_0+ = 1 otherwise t_0+ = 1/t_p0
int nCandMin = 5;  // minimum number of events in one bin
double inclBF = 0.103;  //inclusive b->clnu BF


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

void fit_3D_integ() {

  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 integratio of old FF

  if (!fixSlope) {
    double mB = 5.2790;
    double mD = 1.8646;
    TRandom *rdm1 = new TRandom();
    TRandom2 *rdm2 = new TRandom2();
    for (Int_t i=0;i<N;i++) {
      double y1 = (rdm1->Rndm(i))*2; 
      double q2 = (rdm2->Rndm(i))*11.66;
      
      // Get Fplus
      double mass=mD;
      double mtb, mbb;
      double msd, mx,mb,nf,nfp; 
      double msq,bx2,mbx,mtx;
      double zji,cji,w,gammaji,chiji,betaji_fppfm;
      double rfppfm,rfpmfm,f3fppfm,f3fpmfm,fppfm,fpmfm,al,ai,rcji,f3; 
      double mqm,msb,bb2,mup,mum,bbx2,tm,wt,r2,betaji_fpmfm;
      
      double t=q2;
      
      msb=5.2;
      msd=0.33;
      bb2=0.431*0.431;
      mbb=5.31;
      nf = 4.0;
      
      msq=1.82;
      bx2=0.45*0.45;
      mbx=0.75*2.01+0.25*1.97;
      nfp = 3.0;
      
      mtb = msb + msd;
      mtx = msq + msd;
      mb=5.279;
      mx=mass;
      
      mup=1.0/(1.0/msq+1.0/msb);
      mum=1.0/(1.0/msq-1.0/msb);
      bbx2=0.5*(bb2+bx2);
      tm=(mb-mx)*(mb-mx);
      if ( t>tm ) t=0.99*tm;
      wt=1.0+(tm-t)/(2.0*mbb*mbx);
      
      mqm = 0.1;
      
      double PI = 3.1415;
      double lqcd2 = 0.04;
      double nflav = 4;
      double As0=0.6;
      if ( mqm > 0.6 ) {
	if ( mqm < 1.85 ) {
	  nflav = 3.0;}
	
	As0 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( mqm*mqm/lqcd2);
      }
      
      double As2=0.6;
      if ( msq > 0.6 ) {
	if ( msq < 1.85 ) {
	  nflav = 3.0;}
	
	As2 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msq*msq/lqcd2);
      }
      r2=3.0/(4.0*msb*msq)+3*msd*msd/(2*mbb*mbx*bbx2) + 
	(16.0/(mbb*mbx*(33.0-2.0*nfp)))*
	log(As0/As2);
      
      f3 = sqrt(mtx/mtb)*pow(sqrt(bx2*bb2)/bbx2,1.5) /
	((1.0+r2*(tm-t)/12.0)*(1.0+r2*(tm-t)/12.0));
      
      w = 1.0 + (( tm - t ) / ( 2.0* mb * mx ));
      rcji = ( 1/sqrt(w*w -1 ))*log( w + sqrt( w*w -1 ));
      al = (8.0 / ( 33.0 - 2.0*nfp ))*(w*rcji -1.0 );
      ai = -1.0* ( 6.0/( 33.0 - 2.0*nf));  
      
      double As1=0.6;
      if ( msb > 0.6 ) {
	if ( msb < 1.85 ) {
	  nflav = 3.0;}
	
	As1 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msb*msb/lqcd2);
      }
    
      
      cji = pow(( As1 / As2 ),ai);
      
      zji = msq / msb;
      
      gammaji = 2+((2.0*zji)/(1-zji))*log(zji);
      gammaji = -1.0*gammaji;
      
      chiji = -1.0 - ( gammaji / ( 1- zji ));
      betaji_fppfm = gammaji - (2.0/3.0)*chiji;
      betaji_fpmfm = gammaji + (2.0/3.0)*chiji;
      
      double As3=0.6;
      double m3=sqrt(msb*msq);
      
      if ( m3 > 0.6 ) {
	if ( msq < 1.85 ) {
	  nflav = 3.0;}
	
	As3 = 12.0*PI / ( 33.0 - 2.0*nflav) /
	  log( m3*m3/lqcd2);
      }
      rfppfm = cji *(1.0 + betaji_fppfm*As3/PI);
      
      rfpmfm = cji *(1.0 + betaji_fpmfm*As3/PI);
      
      f3fppfm = f3*pow(( mbb / mtb ),-0.5)*pow((mbx/mtx),0.5);
      f3fpmfm = f3*pow(( mbb / mtb ),0.5)*pow((mbx/mtx),-0.5);
      fppfm = f3fppfm* rfppfm * ( 2.0 - 
				  ( ( mtx/msq)*(1- ( (msd*msq*bb2)
						     /(2.0*mup*mtx*bbx2)))));
      fpmfm = f3fpmfm* rfpmfm * ( mtb/msq) * ( 1 - 
					       ( ( msd*msq*bb2)/
						 ( 2.0*mup*mtx*bbx2)));
      
      double fpf = (fppfm + fpmfm)/2.0;
      
      // Get R
      double mcR = 2*sqrt(mB*mD) / (mB+mD);
      
      // Get old FF
      double oldFF = mcR * fpf;
      double oldFF2 = oldFF*oldFF;
      double rto1 = (mB*mB + mD*mD - t)/(2*mB*mD);
      double rto2 = rto1*rto1;
      double rto3 = rto2 - 1;
      double rto4 = sqrt(rto3*rto3*rto3);
      double oldInt =  rto4*oldFF2;
      if (y1 <= oldInt) N_old++;
    }
  }
  cout << "N_old = " << N_old << 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 vstrt[19] = {0.0210, 0.0560,
			       0.0020, 0.0056, 0.0037, 0.0037,
			       0.0060, 0.0020,
			       0.8,
			       1, 1, 1, 1,
			       1, 1, 1, 1,
			       1, 1};
  static Double_t stp[19] = {0.0001, 0.0001, 
			     0.0001, 0.0001, 0.0001, 0.0001, 
			     0.0001, 0.0001,
			     0.0001,
			     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;
  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,5,ierflg);
  gMinuit->mnparm(nprm[9], "Uncorr", 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;
  gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg);
  printChisq = 0;
  printRN = 0;
  showChisq = 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;
  gMinuit->mnexcm("CALL FCN", &p2 , 1,ierflg);
  printChisq = 0;
  printRN = 0;
  showChisq = 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 * 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;

  // 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;
  if (fixDlnu)  BrB1 = Br_Dlnu;
  else BrB1 = x[0];
  if (fixDstar)   BrB2 = Br_Dstar;
  else BrB2 = x[1];
  if (fixD0star) BrB31 = x[3]*ratio_D0star;
  else BrB31 = x[2];
  if (fixD1)   BrB32 = Br_D1;
  else BrB32 = x[3];
  if (fixD1p)   BrB33 = x[3]*ratio_D1p;
  else BrB33 = x[4];
  if (fixD2star)   BrB34 = x[3]*ratio_D2star;
  else BrB34 = x[5];
  if (fixDPi)   BrB41 = Br_DPi;
  else BrB41 = x[6];
  if (fixDstarPi)   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];
  c[0] = 1;
  c[28] = 1;
  c[2] = 1;
  c[29] = 1;

  c[14] = Bkg1_D0 * fp0;
  c[15] = Bkg1_D0;
  c[16] = Bkg2_D0 * fp0;
  c[17] = Bkg2_D0;
  c[18] = Bkg3_D0 * fp0;
  c[19] = Bkg3_D0;
  c[20] = 1;
  c[21] = 1;
  c[22] = Bkg4_D0;
  c[23] = 1;
  c[24] = 1;
  c[25] = 1;

  c[42] = cDp0 * Bkg1_Dch * fp0;
  c[43] = cDp0 * Bkg1_Dch;
  c[44] = cDp0 * Bkg2_Dch * fp0;
  c[45] = cDp0 * Bkg2_Dch;
  c[46] = cDp0 * Bkg3_Dch * fp0;
  c[47] = cDp0 * Bkg3_Dch;
  c[48] = cDp0;
  c[49] = cDp0;
  c[50] = cDp0 * Bkg4_Dch;
  c[51] = cDp0;
  c[52] = cDp0;
  c[53] = cDp0;

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

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

  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;
  
  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[2] = B_new / B_old;

  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[4] = B_new / B_old;

  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[6] = B_new / B_old;
;
  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[8] = B_new / B_old;

  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[3] = B_new / B_old;

  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[5] = B_new / B_old;

  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[7] = B_new / B_old;

  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[9] = B_new / B_old;


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

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

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

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

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

  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[54] = cDp0 * B_new / B_old;

  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[30] = cDp0 * B_new / B_old;

  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[32] = cDp0 * B_new / B_old;

  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[34] = cDp0 * B_new / B_old;

  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[36] = cDp0 * B_new / B_old;

  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[31] = cDp0 * B_new / B_old;

  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[33] = cDp0 * B_new / B_old;

  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[35] = cDp0 * B_new / B_old;

  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[37] = cDp0 * B_new / B_old;

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

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

  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[40] = cDp0 * B_new / B_old;

  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[41] = cDp0 * B_new / B_old;


  //###########################################################
  // Get normalization factor
  
  double RN = 1;
  if (!fixSlope) {
    TRandom *rdm3 = new TRandom();
    TRandom3 *rdm4 = new TRandom3();
    int N_new = 0;
    for (Int_t i=0;i<N;i++) {
     
      // Get new FF
      double y2 = (rdm3->Rndm(i))*2; 
      double trueW = (rdm4->Rndm(i))*0.592 + 1;
      
      if (Slope == "New") {
	double trueZ 
	  = (sqrt(trueW+1) - sqrt(2))/(sqrt(trueW+1) + sqrt(2));
	double N01_new = 8*rho_D_sq*trueZ;
	double N02_new = (51*rho_D_sq - 10)*trueZ*trueZ;
	double N03_new = (252*rho_D_sq - 84)*trueZ*trueZ*trueZ;
	double newFF = 1 - N01_new + N02_new - N03_new;
      } else {
	double newFF = 1 - rho_D_sq*(trueW - 1);
      }
      double newFF2 = newFF*newFF;
      double rto5 = trueW*trueW - 1;
      double rto6 = sqrt(rto5*rto5*rto5);
      double newInt = rto6*newFF2;
      if (y2 <= newInt) N_new++;
    }
    
    if (N_new > 0) {
      RN = (double) N_old / N_new;
    }
    double RN_err = 0;
    if (N_old > 0 && N_new > 0) {
      double err1 = (double) N_old + N_new;
      double err2 = (double) N_old * N_new;
      double err3 = err1/err2;
      RN_err = RN * sqrt(err3);
    }
    //cout << "N_old = " << N_old << endl;
    //cout << "N_new = " << N_new << endl;
    //cout << "RN = " << RN << endl;
    //cout << "RN_err = " << RN_err << endl;
    if (printRN) cout << "RN = " << RN << " +- " << RN_err << endl;
  }


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


  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]*RN - num_3D[58][i][j][k]*rho_D_sq*RN
	    + num_3D[59][i][j][k]*rho_D_sq*rho_D_sq*RN;
	  num_sub_DchKPiPi[i][j][k]
	    = num_3D[60][i][j][k]*RN - num_3D[61][i][j][k]*rho_D_sq*RN
	    + num_3D[62][i][j][k]*rho_D_sq*rho_D_sq*RN;
	  var_sub_D0KPi[i][j][k] 
	    = err_3D[57][i][j][k]*err_3D[57][i][j][k]*RN
	    - err_3D[58][i][j][k]*err_3D[58][i][j][k]*rho_D_sq*RN
	    + err_3D[59][i][j][k]*err_3D[59][i][j][k]*rho_D_sq*rho_D_sq*RN;
	  var_sub_DchKPiPi[i][j][k] 
	    = err_3D[60][i][j][k]*err_3D[60][i][j][k]*RN
	    - err_3D[61][i][j][k]*err_3D[61][i][j][k]*rho_D_sq*RN
	    + err_3D[62][i][j][k]*err_3D[62][i][j][k]*rho_D_sq*rho_D_sq*RN;
	  err_sub_D0KPi[i][j][k] = sqrt(var_sub_D0KPi[i][j][k]);
	  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;
	}
      }
    }
    double totN_D0KPi(0), totN_DchKPiPi(0);
    double totN_org_D0KPi(0), totN_org_DchKPiPi(0);
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  totN_D0KPi += num_sub_D0KPi[i][j][k];
	  totN_DchKPiPi += num_sub_DchKPiPi[i][j][k];
	  totN_org_D0KPi += num_3D[55][i][j][k];
	  totN_org_DchKPiPi += num_3D[56][i][j][k];
	}
      }
    }
    if (totN_D0KPi > 0 && totN_DchKPiPi > 0) {
      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_sub_D0KPi[i][j][k] * totN_org_D0KPi/totN_D0KPi;
	    var_sub_D0KPi[i][j][k] 
	      = num_sub_D0KPi[i][j][k] 
	      * totN_org_D0KPi*totN_org_D0KPi/(totN_D0KPi*totN_D0KPi);
	    num_sub_DchKPiPi[i][j][k] 
	      = num_sub_DchKPiPi[i][j][k] * totN_org_DchKPiPi/totN_DchKPiPi;
	    var_sub_DchKPiPi[i][j][k] 
	      = num_sub_DchKPiPi[i][j][k] 
	      * totN_org_DchKPiPi*totN_org_DchKPiPi/(totN_DchKPiPi*totN_DchKPiPi);
	  }
	}
      }
    }
    

  }



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

  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;
	}
	//cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << nMC_D0[i][j][k] << endl;
	//cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << varMC_D0[i][j][k] << endl;

	for (int l=2; l<=13; l++) {
	  //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;
	  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;
	  //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;
	}
	//cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << nMC_D0[i][j][k] << endl;
	//cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << varMC_D0[i][j][k] << endl;

	if (fixSlope) {
	  nMC_D0[i][j][k] += c[55]*num_3D[55][i][j][k];
	  double sigm1 = c[55]*err_3D[55][i][j][k];
	} else {
	  nMC_D0[i][j][k] += c[55]*num_sub_D0KPi[i][j][k];
	  double sigm1 = c[55]*err_sub_D0KPi[i][j][k];
	}
	varMC_D0[i][j][k] += sigm1 * sigm1;
	//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;
	//cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << nMC_D0[i][j][k] << endl;
	//cout <<  "c[55] = " << c[55] << endl;
	//cout <<  "sigm1 = " << sigm1 << endl;
	//cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << varMC_D0[i][j][k] << endl;

	for (int l=14; 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;
	  //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;
	}
	//cout << "nMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << nMC_D0[i][j][k] << endl;
	//cout << "varMC_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << varMC_D0[i][j][k] << endl;
       	//cout << "num_3D[0]["<<i<<"]["<<j<<"]["<<k<<"] = " 
	//     << num_3D[0][i][j][k] << endl;

	if (num_3D[0][i][j][k] >= nCandMin && varMC_D0[i][j][k]>0) {
	  nBinsChi_D0++;
	  nCand_D0[i][j][k] = num_3D[0][i][j][k] - nMC_D0[i][j][k];
	  //cout << "nCand_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << nCand_D0[i][j][k] << endl;
	  var_D0[i][j][k] 
	    = err_3D[0][i][j][k] * err_3D[0][i][j][k] + varMC_D0[i][j][k]; 
	  //cout << "var_D0["<<i<<"]["<<j<<"]["<<k<<"] = " 
	  //     << var_D0[i][j][k] << endl;
	  chisq_D0 += (nCand_D0[i][j][k] * nCand_D0[i][j][k] / var_D0[i][j][k]);
	  if (showChisq) cout << "chisq_D0 = " << chisq_D0 << endl;
	}
	//cout << "chisq_D0 = " << chisq_D0 << endl;

	for (int l=30; l<=41; 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 (fixSlope) {
	  nMC_Dch[i][j][k] += c[56]*num_3D[56][i][j][k];
	  double sigm1 = c[56]*err_3D[56][i][j][k];
	} else {
	  nMC_Dch[i][j][k] += c[56]*num_sub_DchKPiPi[i][j][k];
	  double sigm1 = c[56]*err_sub_DchKPiPi[i][j][k];
	}
	varMC_Dch[i][j][k] += sigm1 * sigm1;
	
	for (int l=42; l<=54; 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 (num_3D[28][i][j][k] >= nCandMin && varMC_Dch[i][j][k] > 0) {
	  nBinsChi_Dch++;
	  nCand_Dch[i][j][k] = num_3D[28][i][j][k] - nMC_Dch[i][j][k];
	  var_Dch[i][j][k] 
	    = err_3D[28][i][j][k] * err_3D[28][i][j][k] + varMC_Dch[i][j][k];
	  chisq_Dch += nCand_Dch[i][j][k]*nCand_Dch[i][j][k]/var_Dch[i][j][k];
	  if (showChisq) 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;
  //chisq = chisq + (fp0 - 1)*(fp0 - 1)/0.001 + (cDp0 - 1)*(cDp0 - 1)/0.001
  //  + (BrB1 - 0.0215)*(BrB1 - 0.0215)/(0.0022*0.0022);

  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++;

}