void calc_SBRatio_D0Kpi2() {

#include <iostream.h>
#include <fstream.h>

  gROOT->Reset();

  // Root macro to calculate uncertainty of signal
  //#######################################################
  // Change these values to meet your condition
  //#######################################################
 
  // Info of input data/MC root files
  TString baseDir("/home1x/OtherMounts/hep03/khamano/R16b/Run3");
  //TString baseDir("/home1x/OtherMounts/hep03/khamano/R16b/kCut1");
  TString nameBase_BpBm("SP-1235-BToDlnu-Run3-R16b");
  int nStart_BpBm = 1;
  int nStop_BpBm  = 1193; //218 //87
  TString nameBase_B0B0bar("SP-1237-BToDlnu-Run3-R16b");
  int nStart_B0B0bar = 10000;
  int nStop_B0B0bar  = 855; //220 //88
  TString nameBase_ccbar("SP-1005-BToDlnu-Run3-R16b");
  int nStart_ccbar = 1000;
  int nStop_ccbar  = 208; //208 //83
  TString nameBase_uds("SP-998-BToDlnu-Run3-R16b");
  int nStart_uds = 1000;
  int nStop_uds  = 143; //143 //58
  TString nameBase_tautau("SP-3429-BToDlnu-Run3-R16b");
  int nStart_tautau = 100;
  int nStop_tautau  = 11; //11 //5
  // Name of the root tree in the input root files
  TString ntpName("ntp2");


  // Info of output files
  TString outFileName("SBRatio2_D0Kpi_Run3_test.txt");
  TString outHistoFileName("histo2_nDl_D0Kpi_Run3_test.root");

  // D mass peak region
  float DMassMin = 1.840; //1.855; //1.845; //for D0
  float DMassMax = 1.888; //1.875; //1.885;

  // D mass off peak (side band)region
  float DSideMin1 = 1.816;  //for Dch
  float DSideMax1 = 1.840; 
  float DSideMin2 = 1.888;  //for Dch
  float DSideMax2 = 1.912; 

  // Select D decay mode
  int setDDecMod = 1; // D0Kpi
  //int setDDecMod = 2; // D0K3pi
  //int setDDecMod = 3; // D0Kpipi0
  //int setDDecMod = 4; // D0Kspipi

  // Vertex cut
  float DprobMin = 0.01; //0.001
  float BprobMin = 0.01; //0.001

  // Thrust cut
  float cosThrMax = 0.80;

  // Select a bin
  float DMomMin = 0.2;
  float DMomMax = 2.8;
  float LMomMin = 0.8;
  float LMomMax = 2.8;
  float cosBYMin = -10;
  float cosBYMax = 5;

  // uncertainty on background
  float rho = 0;

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

  // Define output historgams
  int nbin = 5;
  float xmin = 0;
  float xmax = 5;

  TH1D *TH1nDl_total = new TH1D("TH1nDl_total","nDl_total",nbin, xmin,xmax);
  TH1D *TH1nDl_sig = new TH1D("TH1nDl_sig","nDl_sig",nbin, xmin,xmax);
  TH1D *TH1nDl_bkg = new TH1D("TH1nDl_bkg","nDl_bkg",nbin, xmin,xmax);

  // Chain input root files and define trees.
  cout << "start chain\n";

  TChain * chain_MC = new TChain(ntpName);
  for (Int_t i=nStart_BpBm; i<=nStop_BpBm; i++) {
  chain_MC->Add(baseDir+"/"+nameBase_BpBm+"-"+i+".root");
  }
  for (Int_t i=nStart_B0B0bar; i<=nStop_B0B0bar; i++) {
  chain_MC->Add(baseDir+"/"+nameBase_B0B0bar+"-"+i+".root");
  }
  for (Int_t i=nStart_ccbar; i<=nStop_ccbar; i++) {
  chain_MC->Add(baseDir+"/"+nameBase_ccbar+"-"+i+".root");
  }
  for (Int_t i=nStart_uds; i<=nStop_uds; i++) {
  chain_MC->Add(baseDir+"/"+nameBase_uds+"-"+i+".root");
  }
  for (Int_t i=nStart_tautau; i<=nStop_tautau; i++) {
  chain_MC->Add(baseDir+"/"+nameBase_tautau+"-"+i+".root");
  }
  TTree * tree_MC = chain_MC;

  int numDl(0);
  int DDecMod[200], DType[200];
  double fitDprob[200], fitBprob[200];
  float cosThr[200];
  int isABEvent[200], isDirectFromB[200], isFromSameB[200];
  int BDecMod[200], mcDDecMod[200];
  int noExtKPi[200], KMissID[200], PiMissID[200], KPiSameD[200];
  float Dmass[200];
  float pD[200], pLep[200], cosBY[200];
  float pAddLep[200];
  tree_MC->SetBranchAddress("nDl", &numDl);
  tree_MC->SetBranchAddress("DDecayMode", &DDecMod);
  tree_MC->SetBranchAddress("DhadType", &DType);
  tree_MC->SetBranchAddress("FitDprob", &fitDprob);
  tree_MC->SetBranchAddress("FitBprob", &fitBprob);
  tree_MC->SetBranchAddress("cosDlnonDl", &cosThr);
  tree_MC->SetBranchAddress("trueLepFromB", &isABEvent);
  tree_MC->SetBranchAddress("mcDirectLepFromB", &isDirectFromB);
  tree_MC->SetBranchAddress("trueLepFromSameB", &isFromSameB);
  tree_MC->SetBranchAddress("trueBDecMod", &BDecMod);
  tree_MC->SetBranchAddress("trueDDecMod", &mcDDecMod);
  tree_MC->SetBranchAddress("noExtraKorPi", &noExtKPi);
  tree_MC->SetBranchAddress("mcKaonMissID", &KMissID);
  tree_MC->SetBranchAddress("mcPiMissID", &PiMissID);
  tree_MC->SetBranchAddress("mcKPifromSameD", &KPiSameD);
  tree_MC->SetBranchAddress("DMass", &Dmass);
  tree_MC->SetBranchAddress("DMom", &pD);
  tree_MC->SetBranchAddress("DlLeptMom", &pLep);
  tree_MC->SetBranchAddress("DLcosBY", &cosBY);
  tree_MC->SetBranchAddress("AddLepMom", &pAddLep);

  int nSig_Peak(0), nBkg_Peak(0);
  int nBchD_Peak(0), nBchDstar_Peak(0), nBchD2star_Peak(0), nBchD0star_Peak(0);
  int nBchD1_Peak(0), nBchD1H_Peak(0), nBchDpi_Peak(0), nBchDstarpi_Peak(0);
  int nB0D_Peak(0), nB0Dstar_Peak(0), nB0D2star_Peak(0), nB0D0star_Peak(0);
  int nB0D1_Peak(0), nB0D1H_Peak(0), nB0Dpi_Peak(0), nB0Dstarpi_Peak(0);
  int nNonB_Peak(0), nNonB_realD_Peak(0), nLepMissID_Peak(0);
  int nNotDirectB_Peak(0), nNotFrmSameB_Peak(0);
  int nKMissID_Peak(0), nPiMissID_Peak(0), nKPiNotSameD_Peak(0);
  int nExtraKPi_Peak(0), nOther_Peak(0);

  int nSig_Side(0), nBkg_Side(0);
  int nBchD_Side(0), nBchDstar_Side(0), nBchD2star_Side(0), nBchD0star_Side(0);
  int nBchD1_Side(0), nBchD1H_Side(0), nBchDpi_Side(0), nBchDstarpi_Side(0);
  int nB0D_Side(0), nB0Dstar_Side(0), nB0D2star_Side(0), nB0D0star_Side(0);
  int nB0D1_Side(0), nB0D1H_Side(0), nB0Dpi_Side(0), nB0Dstarpi_Side(0);
  int nNonB_Side(0), nNonB_realD_Side(0), nLepMissID_Side(0);
  int nNotDirectB_Side(0), nNotFrmSameB_Side(0);
  int nKMissID_Side(0), nPiMissID_Side(0), nKPiNotSameD_Side(0);
  int nExtraKPi_Side(0), nOther_Side(0);

  Int_t nEvent =(Int_t)tree_MC->GetEntries(); // get number of events
  cout << "nEvent = " << nEvent << endl;
  int nSelEvent(0);
  for (Int_t i=0;i<nEvent;i++) { // loop over events
    if (i==1) cout << "event 1" << endl;
    if (i==100) cout << "event 100" << endl;
    if (i==1000) cout << "event 1000" << endl;
    if (i==10000) cout << "event 10000" << endl;
    if (i==100000) cout << "event 100000" << endl;
    if (i==1000000) cout << "event 1000000" << endl;
    if (i==2000000) cout << "event 2000000" << endl;
    if (i==3000000) cout << "event 3000000" << endl;
    if (i==4000000) cout << "event 4000000" << endl;

    tree_MC->GetEntry(i); //read complete event in memory
    int nDl_total(0), nDl_sig(0), nDl_bkg(0);
    int nDl_side(0);

    for (Int_t j=0;j<numDl;j++) { // loop over Dl candidates
      int DDecMod1 = DDecMod[j];
      int DType1 = DType[j];
      int BDecMod1 = BDecMod[j];
      double fitDprob1 = fitDprob[j]; 
      double fitBprob1 = fitBprob[j]; 
      float cosThr1 = cosThr[j];
      //cout << "cosThr1 = " << cosThr1 << endl;
      int isABEvent1 = isABEvent[j];
      int isDirectFromB1 = isDirectFromB[j];
      int isFromSameB1 = isFromSameB[j];
      int mcDDecMod1 = mcDDecMod[j];
      int noExtKPi1 = noExtKPi[j];
      int KMissID1 = KMissID[j];
      int PiMissID1 = PiMissID[j];
      int KPiSameD1 = KPiSameD[j];
      float Dmass1 = Dmass[j];
      float pD1 = pD[j];
      float pLep1 = pLep[j];
      float cosBY1 = cosBY[j];
      float pAddLep1 = pAddLep[j];

      // Select a bin
      if ((pD1<DMomMin || DMomMax<pD1) || (pLep1<LMomMin || LMomMax<pLep1) 
	  || (cosBY1<cosBYMin || cosBYMax<cosBY1))  continue;

      // D decay mode selection
      if (DDecMod1 != setDDecMod) continue;
      if (DType1 == 1) continue;

      // Vertex cut
      if (fitDprob1 < DprobMin) continue;
      if (fitBprob1 < BprobMin) continue;
      
      // Thrust cut
      if (cosThr1<-cosThrMax || cosThr1>cosThrMax) continue;
      
      // DMass cut
      if (Dmass1<DSideMin1 || DSideMax2<Dmass1) continue;

      // AddLep cut
      //if (pAddLep1<0.1) continue;

      // DMass peak region:
      if (DMassMin<Dmass1 && Dmass1<DMassMax) {

	// Sum up each signal mode and background
	if (isABEvent1 == 0) {
	  if (KMissID1 == 0 && PiMissID1 == 0 && KPiSameD1 == 1
	      && noExtKPi1 == 1) {nNonB_realD_Peak++; nDl_bkg++;}
	  else {nNonB_Peak++; nDl_bkg++;}
	}
	if (isABEvent1 == 1) {
	  if (KMissID1 == 1) {nKMissID_Peak++;  nDl_bkg++;}
	  if (KMissID1 == 0) {
	    if (PiMissID1 == 1) {nPiMissID_Peak++; nDl_bkg++;}
	    if (PiMissID1 == 0) {
	      if (KPiSameD1 == 0) {nKPiNotSameD_Peak++; nDl_bkg++;}
	      if (KPiSameD1 == 1) {
		if (noExtKPi1 == 0) {nExtraKPi_Peak++; nDl_bkg++;}
		if (noExtKPi1 == 1) {
		  if (BDecMod1 == 0) {nLepMissID_Peak++; nDl_bkg++;}
		  if (BDecMod1 != 0) {
		    if (isDirectFromB1 == 0) {nNotDirectB_Peak++; nDl_bkg++;}
		    if (isDirectFromB1 == 1) {
		      if (isFromSameB1 == 0) {nNotFrmSameB_Peak++; nDl_bkg++;}
		      if (isFromSameB1 == 1) {
			if (BDecMod1==1140 || BDecMod1==2130) {
			  nBchD_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1240 || BDecMod1==2230) {
			  nBchDstar_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1340 || BDecMod1==2330) {
			  nBchD2star_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1440 || BDecMod1==2430) {
			  nBchD0star_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1540 || BDecMod1==2530) {
			  nBchD1_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1640 || BDecMod1==2630) {
			  nBchD1H_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==1121 || BDecMod1==1143 
			    || BDecMod1==2112 || BDecMod1==2133) {
			  nBchDpi_Peak++; 
			  nDl_sig++;
			} 
			if (BDecMod1==1221 || BDecMod1==1243 
			    || BDecMod1==2212 || BDecMod1==2233) {
			  nBchDstarpi_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3120 || BDecMod1==4110) {
			  nB0D_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3220 || BDecMod1==4210) {
			  nB0Dstar_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3320 || BDecMod1==4310) {
			  nB0D2star_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3420 || BDecMod1==4410) {
			  nB0D0star_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3520 || BDecMod1==4510) {
			  nB0D1_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3620 || BDecMod1==4610) {
			  nB0D1H_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3123 || BDecMod1==3142
			    || BDecMod1==4113 || BDecMod1==4131) {
			  nB0Dpi_Peak++; 
			  nDl_sig++;
			}
			if (BDecMod1==3223 || BDecMod1==3242 
			    || BDecMod1==4213 || BDecMod1==4231) {
			  nB0Dstarpi_Peak++; 
			  nDl_sig++;
			}
			if (nDl_sig==0)  {
			  nOther_Peak++; 
			  nDl_bkg++;
			}
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}
      }

      // DMass side band:
      if (Dmass1<DMassMin || DMassMax<Dmass1) {

	// Sum up each signal mode and background
	if (isABEvent1 == 0) {
	  if (KMissID1 == 0 && PiMissID1 == 0 && KPiSameD1 == 1
	      && noExtKPi1 == 1) {nNonB_realD_Side++;}
	  else {nNonB_Side++;}
	}
	if (isABEvent1 == 1) {
	  if (KMissID1 == 1) nKMissID_Side++;
	  if (KMissID1 == 0) {
	    if (PiMissID1 == 1) nPiMissID_Side++;
	    if (PiMissID1 == 0) {
	      if (KPiSameD1 == 0) nKPiNotSameD_Side++;
	      if (KPiSameD1 == 1) {
		if (noExtKPi1 == 0) nExtraKPi_Side++;
		if (noExtKPi1 == 1) {
		  if (BDecMod1 == 0) nLepMissID_Side++;
		  if (BDecMod1 != 0) {
		    if (isDirectFromB1 == 0) nNotDirectB_Side++;
		    if (isDirectFromB1 == 1) {
		      if (isFromSameB1 == 0) nNotFrmSameB_Side++;
		      if (isFromSameB1 == 1) {
			if (BDecMod1==1140 || BDecMod1==2130) {
			  nBchD_Side++;
			  nDl_side++;
			}
			if (BDecMod1==1240 || BDecMod1==2230) {
			  nBchDstar_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==1340 || BDecMod1==2330) {
			  nBchD2star_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==1440 || BDecMod1==2430) {
			  nBchD0star_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==1540 || BDecMod1==2530) {
			  nBchD1_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==1640 || BDecMod1==2630) {
			  nBchD1H_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==1121 || BDecMod1==1143 
			    || BDecMod1==2112 || BDecMod1==2133) {
			  nBchDpi_Side++;
			  nDl_side++;
			}
			if (BDecMod1==1221 || BDecMod1==1243 
			    || BDecMod1==2212 || BDecMod1==2233) {
			  nBchDstarpi_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3120 || BDecMod1==4110) {
			  nB0D_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3220 || BDecMod1==4210) {
			  nB0Dstar_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3320 || BDecMod1==4310) {
			  nB0D2star_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3420 || BDecMod1==4410) {
			  nB0D0star_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3520 || BDecMod1==4510) {
			  nB0D1_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3620 || BDecMod1==4610) {
			  nB0D1H_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3123 || BDecMod1==3142
			    || BDecMod1==4113 || BDecMod1==4131) {
			  nB0Dpi_Side++; 
			  nDl_side++;
			}
			if (BDecMod1==3223 || BDecMod1==3242 
			    || BDecMod1==4213 || BDecMod1==4231) {
			  nB0Dstarpi_Side++; 
			  nDl_side++;
			}
			if (nDl_side==0)   nOther_Side++; 
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}
      }

    } //end of Dl loop

    nDl_total = nDl_sig + nDl_bkg;
    if (nDl_total>0) nSelEvent++;
    if (nDl_total>0) TH1nDl_total->Fill(nDl_total);
    if (nDl_sig>0)  TH1nDl_sig->Fill(nDl_sig);
    if (nDl_bkg>0)  TH1nDl_bkg->Fill(nDl_bkg);

  } //end of event loop

  nSig_Peak = nBchD_Peak + nBchDstar_Peak + nBchD2star_Peak 
    + nBchD0star_Peak + nBchD1_Peak + nBchD1H_Peak 
    + nBchDpi_Peak + nBchDstarpi_Peak 
    + nB0D_Peak + nB0Dstar_Peak + nB0D2star_Peak 
    + nB0D0star_Peak + nB0D1_Peak
    + nB0D1H_Peak + nB0Dpi_Peak + nB0Dstarpi_Peak;
  nBkg_Peak = nNonB_Peak + nNonB_realD_Peak + nLepMissID_Peak 
    + nNotDirectB_Peak + nNotFrmSameB_Peak 
    + nKMissID_Peak + nPiMissID_Peak
    + nKPiNotSameD_Peak + nExtraKPi_Peak + nOther_Peak;

  nSig_Side = nBchD_Side + nBchDstar_Side + nBchD2star_Side 
    + nBchD0star_Side + nBchD1_Side + nBchD1H_Side 
    + nBchDpi_Side + nBchDstarpi_Side 
    + nB0D_Side + nB0Dstar_Side + nB0D2star_Side 
    + nB0D0star_Side + nB0D1_Side
    + nB0D1H_Side + nB0Dpi_Side + nB0Dstarpi_Side;
  nBkg_Side = nNonB_Side + nNonB_realD_Side + nLepMissID_Side 
    + nNotDirectB_Side + nNotFrmSameB_Side 
    + nKMissID_Side + nPiMissID_Side
    + nKPiNotSameD_Side + nExtraKPi_Side + nOther_Side;

  float deltaS2 = nSig_Peak + nBkg_Peak + rho*rho*nBkg_Peak*nBkg_Peak;
  float Sig = nSig_Peak;
  float SBRatio(0);
  if (Sig>0) SBRatio = sqrt(deltaS2)/Sig;

  float deltaS2_sub = nSig_Peak + nBkg_Peak + nSig_Side + nBkg_Side
    + rho*rho*nBkg_Peak*nBkg_Peak;
  float Sig_sub = nSig_Peak - nSig_Side;
  float SBRatio_sub(0);
  if (Sig_sub>0) SBRatio_sub = sqrt(deltaS2_sub)/Sig_sub;

  cout << "nSig_Peak = " << nSig_Peak << endl;
  cout << "nBkg_Peak = " << nBkg_Peak << endl;
  cout << "SBRatio = " << SBRatio << endl;

  // Open output file and save results. #####
  ofstream outFile(outFileName); //open output file
  if (!outFile) {
    cout << "Cannot open file " << outFile << endl;
  }
  outFile.setf(ios::fixed);
  outFile.setf(ios::showpoint);
  outFile << "D mass peak region" << endl;
  outFile << "number of events = " << nEvent << endl;
  outFile << "number of selected events = " << nSelEvent << endl;
  outFile << "number of Bp->Dlnu = " << nBchD_Peak << endl;
  outFile << "number of Bp->D*lnu = " << nBchDstar_Peak << endl;
  outFile << "number of Bp->D0**lnu = " << nBchD0star_Peak << endl;
  outFile << "number of Bp->D1**lnu = " << nBchD1_Peak << endl;
  outFile << "number of Bp->D1H**lnu = " << nBchD1H_Peak << endl;
  outFile << "number of Bp->D2**lnu = " << nBchD2star_Peak << endl;
  outFile << "number of Bp->DPilnu = " << nBchDpi_Peak << endl;
  outFile << "number of Bp->D*Pilnu = " << nBchDstarpi_Peak << endl;
  outFile << "number of B0->Dlnu = " << nB0D_Peak << endl;
  outFile << "number of B0->D*lnu = " << nB0Dstar_Peak << endl;
  outFile << "number of B0->D0**lnu = " << nB0D0star_Peak << endl;
  outFile << "number of B0->D1**lnu = " << nB0D1_Peak << endl;
  outFile << "number of B0->D1H**lnu = " << nB0D1H_Peak << endl;
  outFile << "number of B0->D2**lnu = " << nB0D2star_Peak << endl;
  outFile << "number of B0->DPilnu = " << nB0Dpi_Peak << endl;
  outFile << "number of B0->D*Pilnu = " << nB0Dstarpi_Peak << endl;
  outFile << "number of signal = " << nSig_Peak << endl;
  outFile << "number of non-B without real D = " << nNonB_Peak << endl;
  outFile << "number of non-B with real D = " << nNonB_realD_Peak << endl;
  outFile << "number of kaon miss ID = " << nKMissID_Peak << endl;
  outFile << "number of pi miss ID = " << nPiMissID_Peak << endl;
  outFile << "number of K and Pi not from same D = " 
	  << nKPiNotSameD_Peak << endl;
  outFile << "number of extra K or Pi from same D = " 
	  << nExtraKPi_Peak << endl;
  outFile << "number of lepton miss ID = " << nLepMissID_Peak << endl;
  outFile << "number of lepton not direct from B = " 
	  << nNotDirectB_Peak << endl;
  outFile << "number of lepton and D not from same B = " 
	  << nNotFrmSameB_Peak << endl;
  outFile << "number of other background = " << nOther_Peak << endl;
  outFile << "number of background = " << nBkg_Peak << endl;
  outFile << "Delta_S / S = " << SBRatio << endl;
  outFile << "Delta_S / S (sideband subtracted) = " << SBRatio_sub << endl;
  outFile << "-----------------------------------" << endl;
  outFile << "D mass side band" << endl;
  outFile << "number of Bp->Dlnu = " << nBchD_Side << endl;
  outFile << "number of Bp->D*lnu = " << nBchDstar_Side << endl;
  outFile << "number of Bp->D0**lnu = " << nBchD0star_Side << endl;
  outFile << "number of Bp->D1**lnu = " << nBchD1_Side << endl;
  outFile << "number of Bp->D1H**lnu = " << nBchD1H_Side << endl;
  outFile << "number of Bp->D2**lnu = " << nBchD2star_Side << endl;
  outFile << "number of Bp->DPilnu = " << nBchDpi_Side << endl;
  outFile << "number of Bp->D*Pilnu = " << nBchDstarpi_Side << endl;
  outFile << "number of B0->Dlnu = " << nB0D_Side << endl;
  outFile << "number of B0->D*lnu = " << nB0Dstar_Side << endl;
  outFile << "number of B0->D0**lnu = " << nB0D0star_Side << endl;
  outFile << "number of B0->D1**lnu = " << nB0D1_Side << endl;
  outFile << "number of B0->D1H**lnu = " << nB0D1H_Side << endl;
  outFile << "number of B0->D2**lnu = " << nB0D2star_Side << endl;
  outFile << "number of B0->DPilnu = " << nB0Dpi_Side << endl;
  outFile << "number of B0->D*Pilnu = " << nB0Dstarpi_Side << endl;
  outFile << "number of signal = " << nSig_Side << endl;
  outFile << "number of non-B without real D = " << nNonB_Side << endl;
  outFile << "number of non-B with real D = " << nNonB_realD_Side << endl;
  outFile << "number of kaon miss ID = " << nKMissID_Side << endl;
  outFile << "number of pi miss ID = " << nPiMissID_Side << endl;
  outFile << "number of K and Pi not from same D = " 
	  << nKPiNotSameD_Side << endl;
  outFile << "number of extra K or Pi from same D = " 
	  << nExtraKPi_Side << endl;
  outFile << "number of lepton miss ID = " << nLepMissID_Side << endl;
  outFile << "number of lepton not direct from B = " 
	  << nNotDirectB_Side << endl;
  outFile << "number of lepton and D not from same B = " 
	  << nNotFrmSameB_Side << endl;
  outFile << "number of other background = " << nOther_Side << endl;
  outFile << "number of background = " << nBkg_Side << endl;
  outFile.close();

  // Open output file and write histograms in it. ######
  TFile *f = new TFile(outHistoFileName, "Recreate");

  TH1nDl_total->Write();
  TH1nDl_sig->Write();
  TH1nDl_bkg->Write();

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


}