void histMk_1D_DMass_BpBm() {

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

  gROOT->Reset();

  // Root macro to produce D mass plots
  //#######################################################
  // Change these values to meet your condition
  //#######################################################
 
  // Output histogram name
  TString outFileNameBase("histo_1D_DMass_BpBm");
  TString outHistoFileName(outFileNameBase+".root");
  TString outTxtFileName(outFileNameBase+".txt");

  // Info of input data/MC root files
  TString baseDir("/home1x/OtherMounts/hep03/khamano/R16b/Run3");
  TString nameBase("SP-1235-BToDlnu-Run3-R16b");
  int nStart = 1;
  int nStop  = 1193; //1193

  // Name of the root tree in the input root files
  TString ntpName("ntp2");

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

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

  TChain * chainMC = new TChain(ntpName);
  for (Int_t i=nStart; i<=nStop; i++) {
  chainMC->Add(baseDir+"/"+nameBase+"-"+i+".root");
  }
  TTree * treeMC = chainMC;


  // Define histograms
  int nbin = 70;
  double xmin = 1.8;
  double xmax = 1.94;


  TH1D *Dlnu_D0 = new TH1D("Dlnu_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *Dstar_D0 = new TH1D("Dstar_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *D0star_D0 = new TH1D("D0star_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *D1_D0 = new TH1D("D1_D0","D mass (D0)",nbin, xmin,xmax); 
  TH1D *D1p_D0 = new TH1D("D1p_D0","D mass (D0)",nbin, xmin,xmax);  
  TH1D *D2star_D0 = new TH1D("D2star_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *DPi_D0 = new TH1D("DPi_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *DstarPi_D0 = new TH1D("DstarPi_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *DiffB_D0 = new TH1D("DiffB_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *CascL_D0 = new TH1D("CascL_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *LMisID_D0 = new TH1D("LMisID_D0","D mass (D0)",nbin, xmin,xmax);
  TH1D *Comb_D0 = new TH1D("Comb_D0","D mass (D0)",nbin, xmin,xmax);

  TH1D *D0star_Dch = new TH1D("D0star_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *D1_Dch = new TH1D("D1_Dch","D mass (Dch)",nbin, xmin,xmax); 
  TH1D *D1p_Dch = new TH1D("D1p_Dch","D mass (Dch)",nbin, xmin,xmax);  
  TH1D *D2star_Dch = new TH1D("D2star_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *DPi_Dch = new TH1D("DPi_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *DstarPi_Dch = new TH1D("DstarPi_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *DiffB_Dch = new TH1D("DiffB_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *CascL_Dch = new TH1D("CascL_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *LMisID_Dch = new TH1D("LMisID_Dch","D mass (Dch)",nbin, xmin,xmax);
  TH1D *Comb_Dch = new TH1D("Comb_Dch","D mass (Dch)",nbin, xmin,xmax);


  // Define cuts ####################
  // Cuts on D decay modes:
  TCut * cut_D0 = new TCut("DDecayMode==1"); //D0 -> Kpi
  TCut * cut_Dch = new TCut("DDecayMode==5"); //Dch -> Kpipi

  // Cobinatrial background
  TCut * cut_Comb = new TCut("noExtraKorPi==0 || mcKaonMissID==1 || mcPiMissID==1 || mcKPifromSameD==0");

  // LMisID background
  TCut * cut_noComb = new TCut("noExtraKorPi==1 && mcKaonMissID==0 && mcPiMissID==0 && mcKPifromSameD==1");
  TCut * cut_LMisID = new TCut(*cut_noComb && "mcLepMissID==1");

  // CascL background 
  TCut * cut_noLMisID = new TCut(*cut_noComb && "mcLepMissID==0");
  TCut * cut_CascL = new TCut(*cut_noLMisID && "mcDirectLepFromB==0");

  // DiffB background
  TCut * cut_noCascL = new TCut(*cut_noLMisID && "mcDirectLepFromB==1");
  TCut * cut_DiffB = new TCut(*cut_noCascL && "trueLepFromSameB==0");

  // Signal mode
  TCut *cut_Signal = new TCut(*cut_noCascL && "trueLepFromSameB==1");
  TCut *cut_Dlnu = new TCut(*cut_Signal && "trueBDecMod==1140 || trueBDecMod==2130");
  TCut *cut_Dstar = new TCut(*cut_Signal && "trueBDecMod==1240 || trueBDecMod==2230");
  TCut *cut_D2star = new TCut(*cut_Signal && "trueBDecMod==1340 || trueBDecMod==2330");
  TCut *cut_D0star = new TCut(*cut_Signal && "trueBDecMod==1440 || trueBDecMod==2430");
  TCut *cut_D1 = new TCut(*cut_Signal && "trueBDecMod==1540 || trueBDecMod==2530");
  TCut *cut_D1p = new TCut(*cut_Signal && "trueBDecMod==1640 || trueBDecMod==2630");
  TCut *cut_DPi = new TCut(*cut_Signal && "trueBDecMod==1121 || trueBDecMod==1143 || trueBDecMod==2112 || trueBDecMod==2133");
  TCut *cut_DstarPi = new TCut(*cut_Signal && "trueBDecMod==1221 || trueBDecMod==1243 || trueBDecMod==2212 || trueBDecMod==2233");


  // Project trees
  treeMC->Project("Dlnu_D0","DMass",*cut_D0&&*cut_Dlnu);
  treeMC->Project("Dstar_D0","DMass",*cut_D0&&*cut_Dstar);
  treeMC->Project("D0star_D0","DMass",*cut_D0&&*cut_D0star);
  treeMC->Project("D1_D0","DMass",*cut_D0&&*cut_D1);
  treeMC->Project("D1p_D0","DMass",*cut_D0&&*cut_D1p);
  treeMC->Project("D2star_D0","DMass",*cut_D0&&*cut_D2star);
  treeMC->Project("DPi_D0","DMass",*cut_D0&&*cut_DPi);
  treeMC->Project("DstarPi_D0","DMass",*cut_D0&&*cut_DstarPi);
  treeMC->Project("DiffB_D0","DMass",*cut_D0&&*cut_DiffB);
  treeMC->Project("CascL_D0","DMass",*cut_D0&&*cut_CascL);
  treeMC->Project("LMisID_D0","DMass",*cut_D0&&*cut_LMisID);
  treeMC->Project("Comb_D0","DMass",*cut_D0&&*cut_Comb);

  treeMC->Project("D0star_Dch","DMass",*cut_Dch&&*cut_D0star);
  treeMC->Project("D1_Dch","DMass",*cut_Dch&&*cut_D1);
  treeMC->Project("D1p_Dch","DMass",*cut_Dch&&*cut_D1p);
  treeMC->Project("D2star_Dch","DMass",*cut_Dch&&*cut_D2star);
  treeMC->Project("DPi_Dch","DMass",*cut_Dch&&*cut_DPi);
  treeMC->Project("DstarPi_Dch","DMass",*cut_Dch&&*cut_DstarPi);
  treeMC->Project("DiffB_Dch","DMass",*cut_Dch&&*cut_DiffB);
  treeMC->Project("CascL_Dch","DMass",*cut_Dch&&*cut_CascL);
  treeMC->Project("LMisID_Dch","DMass",*cut_Dch&&*cut_LMisID);
  treeMC->Project("Comb_Dch","DMass",*cut_Dch&&*cut_Comb);


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

  Dlnu_D0->Write();
  Dstar_D0->Write();
  D0star_D0->Write();
  D1_D0->Write();
  D1p_D0->Write();
  D2star_D0->Write();
  DPi_D0->Write();
  DstarPi_D0->Write();
  DiffB_D0->Write();
  CascL_D0->Write();
  LMisID_D0->Write();
  Comb_D0->Write();

  D0star_Dch->Write();
  D1_Dch->Write();
  D1p_Dch->Write();
  D2star_Dch->Write();
  DPi_Dch->Write();
  DstarPi_Dch->Write();
  DiffB_Dch->Write();
  CascL_Dch->Write();
  LMisID_Dch->Write();
  Comb_Dch->Write();

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


  //#######################################################################
  // Open output txt file and save parameters. #####
  ofstream outFile(outTxtFileName); //open output file
  if (!outFile) {
    cout << "Cannot open file " << outFile << endl;
  }
  outFile.setf(ios::fixed);
  outFile.setf(ios::showpoint);
  outFile << "outFileNameBase = " << outFileNameBase << endl;
  outFile.close();

}