#include <iostream>
#include <fstream>

gROOT->Reset();

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

void histMk_3D_B0B0_Dlnu_DchKPiPi() {

  // Root macro to produce 3D histogram of B0B0 MC.
  // B0B0 MC is luminosity nomalized to OnPeak data.
  //###########################################################
  // Change these values to meet your condition
  //###########################################################
 
  TString RunNum("RunAll");
  if (run == 1) RunNum = "Run1";
  if (run == 2) RunNum = "Run2";
  if (run == 3) RunNum = "Run3";
  if (run == 4) RunNum = "Run4";

  TString MCCorr("noCorr");
  if (correctTrk) MCCorr ="effCorr";

  // Name base of output files
  TString outFileNameBase("histo_3D_"+RunNum+"_"+MCCorr);

  if (newSlope) {
    TString Slope("New");
  } else {
    TString Slope("Old");
  }

  // Info of output files
  TString out3DHistoName(outFileNameBase+"_B0B0_Dlnu_"+DlnuFF+Slope+"_DchKPiPi.root");
  TString outInfoFileName(outFileNameBase+"_B0B0_Dlnu_"+DlnuFF+Slope+"_DchKPiPi.txt");
  TString histTitle_peak("B0B0 to Dlnu (DchKPipi) D mass peak");
  TString histTitle_side("B0B0 to Dlnu (DchKPipi) D mass sideband");
  TString histTitle_sub("B0B0 to Dlnu (DchKPipi) sideband subtracted");

  // B decay mode selection
  int selBDecMod1 = 3120;
  int selBDecMod2 = 4110;
  int selBDecMod3 = -1;
  int selBDecMod4 = -1;

  // D decay mode selection
  int DDecMod_3D = DDecMod_DchKPiPi;

  // D mass range
  float DMass1 = DMass_DchKPiPi1; // lower edge of sideband
  float DMass2 = DMass_DchKPiPi2; // lower edge of D mass peak
  float DMass3 = DMass_DchKPiPi3; // upper edge of D mass peak
  float DMass4 = DMass_DchKPiPi4; // upper edge of sideband

  // Input files
  TString inFileName_Run1(baseDir_Run1+"/"+inFile_Run1_B0B0);
  int nStart_Run1 = nStart_Run1_B0B0;
  int nStop_Run1 = nStop_Run1_B0B0;
  TString inFileName_Run2(baseDir_Run2+"/"+inFile_Run2_B0B0);
  int nStart_Run2 = nStart_Run2_B0B0;
  int nStop_Run2 = nStop_Run2_B0B0;
  TString inFileName_Run3(baseDir_Run3+"/"+inFile_Run3_B0B0);
  int nStart_Run3 = nStart_Run3_B0B0;
  int nStop_Run3 = nStop_Run3_B0B0;
  TString inFileName_Run4(baseDir_Run4+"/"+inFile_Run4_B0B0);
  int nStart_Run4 = nStart_Run4_B0B0;
  int nStop_Run4 = nStop_Run4_B0B0;

  // Luminosity normalization factors
  double c_Run1 = c_Run1_B0B0;
  double c_Run2 = c_Run2_B0B0;
  double c_Run3 = c_Run3_B0B0;
  double c_Run4 = c_Run4_B0B0;
  double c_All = c_All_B0B0;


  //###########################################################
  // Get normalization factor
  
  double RN = 1;
  double mB = 5.2790;
  double mD = 1.8646;
  TRandom *rdm1 = new TRandom();
  TRandom2 *rdm2 = new TRandom2();
  TRandom *rdm3 = new TRandom();
  TRandom3 *rdm4 = new TRandom3();
  int N = 100000;
  int N_old = 0;
  int N_new = 0;
  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++;

    // Get new FF
    double y2 = (rdm3->Rndm(i))*2; 
    double trueW = (rdm4->Rndm(i))*0.592 + 1;

    if (newSlope) {
      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_new > 0 && N_old > 0) {
    RN_err = sqrt((double)(N_old + N_new)/(N_old*N_new));
  }
  cout << "N_old = " << N_old << endl;
  cout << "N_new = " << N_new << endl;
  cout << "RN = " << RN << endl;
  cout << "RN_err = " << RN_err << endl;


  //#######################################################################
  // Create D mass peak and sideband arrays, readin data and fill the arrays.

  cout << "starat B0B0ToDlnu\n";
  cout << "start chain\n";
  if (run == 1) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run1; i<=nStop_Run1; i++) {
      chain_3D->Add(inFileName_Run1+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run1;
  } else if (run == 2) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run2; i<=nStop_Run2; i++) {
      chain_3D->Add(inFileName_Run2+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run2;
  } else if (run == 3) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run3; i<=nStop_Run3; i++) {
      chain_3D->Add(inFileName_Run3+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run3;
  } else if (run == 4) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run4; i<=nStop_Run4; i++) {
      chain_3D->Add(inFileName_Run4+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run4;
  } else {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run1; i<=nStop_Run1; i++) {
      chain_3D->Add(inFileName_Run1+"-"+i+".root");
    }
    for (Int_t i=nStart_Run2; i<=nStop_Run2; i++) {
      chain_3D->Add(inFileName_Run2+"-"+i+".root");
    }
    for (Int_t i=nStart_Run3; i<=nStop_Run3; i++) {
      chain_3D->Add(inFileName_Run3+"-"+i+".root");
    }
    for (Int_t i=nStart_Run4; i<=nStop_Run4; i++) {
      chain_3D->Add(inFileName_Run4+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_All;
  }





  // Define and initialise arrays for 3D histo
  double num_peak[nbinsx][nbinsy][nbinsz];
  double num_side[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_peak[i][j][k] = 0;
	num_side[i][j][k] = 0;
      }
    }
  }

  cout << "read tree\n";
  // variables to read out
  Int_t nDl_3D(0);
  Int_t DDecayMode_3D[200];
  Float_t DMass_3D[200];
  Int_t noExtraKorPi_3D[200], mcKaonMissID_3D[200];
  Int_t mcPiMissID_3D[200], mcKPifromSameD_3D[200];
  Int_t mcLepMissID_3D[200], mcDirectLepFromB_3D[200];
  Int_t trueLepFromSameB_3D[200];
  Int_t trueBDecMod_3D[200];
  double trkEffCorr_3D[200], PIDCorr_3D[200];
  double mcW_3D[200];
  Float_t varx_3D[200], vary_3D[200], varz_3D[200];
  
  // Set reference between the above variables and root-tuple entries
  tree_3D->SetBranchAddress("nDl", &nDl_3D);
  tree_3D->SetBranchAddress("DDecayMode", &DDecayMode_3D);
  tree_3D->SetBranchAddress("DMass", &DMass_3D);
  tree_3D->SetBranchAddress("noExtraKorPi", &noExtraKorPi_3D);
  tree_3D->SetBranchAddress("mcKaonMissID", &mcKaonMissID_3D);
  tree_3D->SetBranchAddress("mcPiMissID", &mcPiMissID_3D);
  tree_3D->SetBranchAddress("mcKPifromSameD", &mcKPifromSameD_3D);
  tree_3D->SetBranchAddress("mcLepMissID", &mcLepMissID_3D);
  tree_3D->SetBranchAddress("mcDirectLepFromB", &mcDirectLepFromB_3D);
  tree_3D->SetBranchAddress("trueBDecMod", &trueBDecMod_3D);
  tree_3D->SetBranchAddress("trueLepFromSameB", &trueLepFromSameB_3D);
  tree_3D->SetBranchAddress("trkEffCorr", &trkEffCorr_3D);
  tree_3D->SetBranchAddress("PIDCorr", &PIDCorr_3D);
  tree_3D->SetBranchAddress("mcW", &mcW_3D);
  tree_3D->SetBranchAddress(varx, &varx_3D);
  tree_3D->SetBranchAddress(vary, &vary_3D);
  tree_3D->SetBranchAddress(varz, &varz_3D);

  // Read out root-tuple entries and fill arrays
  Int_t nevent_3D =(Int_t)tree_3D->GetEntries(); // number of events
  cout << "number of 3D events = " << nevent_3D << endl;
  int numCand(0), numSelCand(0);
  int numCand_DDecMod(0), numCand_Index(0);
  for (Int_t i=0;i<nevent_3D;i++) { // loop over events
    if (i==1) cout << "3D event 1" << endl;
    if (i==100) cout << "3D event 100" << endl;
    if (i==1000) cout << "3D event 1000" << endl;
    if (i==10000) cout << "3D event 10000" << endl;
    if (i==100000) cout << "3D event 100000" << endl;
    if (i==1000000) cout << "3D event 1000000" << endl;
    if (i==2000000) cout << "3D event 2000000" << endl;
    if (i==3000000) cout << "3D event 3000000" << endl;
    if (i==4000000) cout << "3D event 4000000" << endl;
    
    tree_3D->GetEntry(i); //read complete event in memory

    for (Int_t j=0;j<nDl_3D;j++) { // loop over Dl candidates
      int DDecMod = DDecayMode_3D[j];
      float DhadMass = DMass_3D[j];
      int noExtKorPi = noExtraKorPi_3D[j];
      int kaonMissID = mcKaonMissID_3D[j];
      int piMissID = mcPiMissID_3D[j];
      int KPiFromSameD = mcKPifromSameD_3D[j];
      int lepMissID = mcLepMissID_3D[j];
      int directLepFromB = mcDirectLepFromB_3D[j];
      int lepFromSameB = trueLepFromSameB_3D[j];
      int mcBDecMod = trueBDecMod_3D[j];
      double trkCorr = trkEffCorr_3D[j] * trkCorrModif;
      double PidCorr = PIDCorr_3D[j] * PidCorrModif;
      double trueW = mcW_3D[j];
      float xvar = varx_3D[j];
      float yvar = vary_3D[j];
      float zvar = varz_3D[j];

      numCand++;
      if (DDecMod != DDecMod_3D) continue;
      numCand_DDecMod++;
      if (noExtKorPi == 0) continue;
      if (kaonMissID == 1) continue;
      if (piMissID == 1) continue;
      if (KPiFromSameD == 0) continue;
      if (lepMissID == 1) continue;
      if (directLepFromB == 0) continue;
      if (lepFromSameB == 0) continue;
      int goodBDecMod(0);
      if (mcBDecMod == selBDecMod1 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod2 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod3 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod4 )  goodBDecMod=1;
      if (goodBDecMod != 1) continue;

      // Get indeces of this candidate
      //cout << "getting indecies" << endl;
      int xIndex(-1);
      for (Int_t k=0;k<nbinsx;k++) {
	if (xbins[k]<  xvar &&  xvar <= xbins[k+1]) {
	  xIndex = k;
	}
      }
      int yIndex(-1);
      for (Int_t k=0;k<nbinsy;k++) {
	if (ybins[k]< yvar && yvar <= ybins[k+1]) {
	  yIndex = k;
	}
      }
      int zIndex(-1);
      for (Int_t k=0;k<nbinsz;k++) {
	if (zbins[k]< zvar && zvar <= zbins[k+1]) {
	  zIndex = k;
	}
      }

      numCand_Index++;

      if (xIndex>=0 && yIndex>=0 && zIndex>=0) {
	numSelCand++;

	// Get FF correction factor
	double W_D = 1;
	
	if (newDlnuFF) {
	  // Get q-square
	  double mB = 5.2790;
	  double mD = 1.8646;
	  double q2 = mB*mB + mD*mD - 2*mB*mD*trueW;
	  //cout << "q2 = " << q2 << endl;
	  
	  // 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.87;
	  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;
	  //cout << "fpf = " << fpf << endl;
	  
	  // Get R
	  double mcR = 2*sqrt(mB*mD) / (mB+mD);

	  // Get oldFF
	  double oldFF = mcR * fpf;
	  //cout << "oldFF = " << oldFF << endl;
	  
	  // Get newFF
	  if (newSlope) {
	    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);
	  }
	  
	  // Get weight
	  double W_D_root = newFF / oldFF; 
	  W_D = W_D_root * W_D_root * RN;
	}
	  
	if (DMass2<DhadMass && DhadMass<DMass3) {
	  if (correctTrk) {
	    num_peak[xIndex][yIndex][zIndex] 
	      += (trkCorr * PidCorr)*W_D;
	  } else {
	    num_peak[xIndex][yIndex][zIndex] += W_D;
	  }
	}
	if ((DMass1<DhadMass && DhadMass<DMass2) 
	    ||(DMass3<DhadMass && DhadMass<DMass4)){
	  if (correctTrk) {
	    num_side[xIndex][yIndex][zIndex] 
	      += (trkCorr * PidCorr)*W_D;
	  } else {
	    num_side[xIndex][yIndex][zIndex] += W_D;
	  }
	}
      }
      
    } //end of Dl loop
  } //end of event loop
  cout << "number of processed candidates = " << numCand << endl;
  cout << "number of candidates after D decay mode selection = " 
       << numCand_DDecMod << endl;
  cout << "number of candidates before index selection = " 
       << numCand_Index << endl;
  cout << "number of selected candidates = " << numSelCand << endl;
  

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

  // Define and fill histos
  TH3D *peak 
    = new TH3D("peak", histTitle_peak,
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
  TH3D *side 
    = new TH3D("side", histTitle_side, 
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
  TH3D *sub 
    = new TH3D("sub", histTitle_sub, 
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
 
  // Now fill the histograms
  cout << "filling histograms\n";
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	peak->SetBinContent(i+1,j+1,k+1,num_peak[i][j][k]);
        side->SetBinContent(i+1,j+1,k+1,num_side[i][j][k]);

	double peak3D = conv * num_peak[i][j][k];
	double side3D = conv * num_side[i][j][k];
	double sub3D = peak3D - side3D;
	double var3D = num_peak[i][j][k] + num_side[i][j][k];
	double err3D = conv * sqrt(var3D);
        sub->SetBinContent(i+1,j+1,k+1,sub3D);
        sub->SetBinError(i+1,j+1,k+1,err3D);
      }
    }
  }

  //#############################################################
  // Open output file and write histograms in it. ######
  cout <<"saving histograms to " << out3DHistoName << endl;
  TFile *f = new TFile(out3DHistoName, "Recreate");

  peak->Write();
  side->Write();
  sub->Write();

  f->Write();
  f->Close();

  //#######################################################################
  // Open output file and save parameters. #####
  cout << "dumping parameters to " << outInfoFileName << endl;
  ofstream outFile(outInfoFileName); //open output file
  if (!outFile) {
    cout << "Cannot open file " << outFile << endl;
  }
  outFile.setf(ios::fixed);
  outFile.setf(ios::showpoint);
  outFile << "Input files" << endl;
  outFile << "  run = " << run << endl;
  if (run == 1) {
    outFile << "  File name base = " << inFileName_Run1 << endl;
    outFile << "  start number = " << nStart_Run1 << endl;
    outFile << "  stop number = " << nStop_Run1 << endl;
  } else if (run == 2) {
    outFile << "  File name base = " << inFileName_Run2 << endl;
    outFile << "  start number = " << nStart_Run2 << endl;
    outFile << "  stop number = " << nStop_Run2 << endl;
  } else if (run == 3) {
    outFile << "  File name base = " << inFileName_Run3 << endl;
    outFile << "  start number = " << nStart_Run3 << endl;
    outFile << "  stop number = " << nStop_Run3 << endl;
  } else if (run == 4) {
    outFile << "  File name base = " << inFileName_Run4 << endl;
    outFile << "  start number = " << nStart_Run4 << endl;
    outFile << "  stop number = " << nStop_Run4 << endl;
  } else {
    outFile << "  File name base = " << inFileName_Run1 << endl;
    outFile << "  start number = " << nStart_Run1 << endl;
    outFile << "  stop number = " << nStop_Run1 << endl;
    outFile << "  File name base = " << inFileName_Run2 << endl;
    outFile << "  start number = " << nStart_Run2 << endl;
    outFile << "  stop number = " << nStop_Run2 << endl;
    outFile << "  File name base = " << inFileName_Run3 << endl;
    outFile << "  start number = " << nStart_Run3 << endl;
    outFile << "  stop number = " << nStop_Run3 << endl;
    outFile << "  File name base = " << inFileName_Run4 << endl;
    outFile << "  start number = " << nStart_Run4 << endl;
    outFile << "  stop number = " << nStop_Run4 << endl;
  }
  outFile << "################" << endl;
  outFile << "Output files" << endl;
  outFile << "  rho_D_sq = " << rho_D_sq << endl;
  outFile << "  histograms : " << out3DHistoName << endl;
  outFile << "  information : " << outInfoFileName << endl;
  outFile << "################" << endl;
  outFile << "Normalization info" << endl;
  outFile << "  N_old = " << N_old << endl;
  outFile << "  N_new = " << N_new << endl;
  outFile << "  RN = " << RN << endl;
  outFile << "  RN_err = " << RN_err << endl;
  outFile << "################" << endl;
  outFile << "Event number info" << endl;
  outFile << "  number of events processed = " << nevent_3D << endl;
  outFile << "  number of candidates processed = " << numCand << endl;
  outFile << "  number of candidates after D decay mode selection = " 
	  << numCand_DDecMod << endl;
  outFile << "  number of candidates before index selection = " 
	  << numCand_Index << endl;
  outFile << "  number of candidates selected = " << numSelCand << endl;
  outFile << "################" << endl;
  outFile << "Luminosity normalization factor" << endl;
  outFile << "  conv = " << conv << endl;
  outFile << "B decay mode selection" << endl;
  outFile << "  selBDecMod1 = " << selBDecMod1 << endl;
  outFile << "  selBDecMod2 = " << selBDecMod2 << endl;
  outFile << "  selBDecMod3 = " << selBDecMod3 << endl;
  outFile << "  selBDecMod4 = " << selBDecMod4 << endl;
  outFile << "D decay mode selection" << endl;
  outFile << "  DDecMod_3D = " << DDecMod_3D << endl;
  outFile << "D mass cut" << endl;
  outFile << "  DMass1 = " << DMass1 << endl;
  outFile << "  DMass2 = " << DMass2 << endl;
  outFile << "  DMass3 = " << DMass3 << endl;
  outFile << "  DMass4 = " << DMass4 << endl;
  outFile << "3 variables" << endl;
  outFile << "  Variable1 = " << varx << endl;
  outFile << "  Variable2 = " << vary << endl;
  outFile << "  Variable3 = " << varz << endl;
  outFile << "################" << endl;
  outFile << "Do or do not tracking and PID correction" << endl;
  outFile << "  correctTrk = " << correctTrk << endl;
  outFile << "track and PID correction modificatin factors" << endl;
  outFile << "modified correction = Corr * Modif" << endl;
  outFile << "  trkCorrModif = " << trkCorrModif << endl;
  outFile << "  PidCorrModif = " << PidCorrModif << endl;
  outFile << "FF weight W = [1 - rho_D_sq*(w-1)]^2" << endl;
  outFile << "  slope rho_D_sq = " << rho_D_sq << endl;
  outFile.close();

}