void histMk_1D_LepMom_OnPeak() {

#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_LepMom_OnPeak");
  TString outHistoFileName(outFileNameBase+".root");
  TString outTxtFileName(outFileNameBase+".txt");

  TString plotVal("DlLeptMom");
  TString title_D0("Lepton momentum (D0)");
  TString title_Dch("Lepton momentum (Dch)");

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

  // 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 = 40;
  double xmin = 0.8;
  double xmax = 2.8;

  TH1D *OnPeak_p_D0 = new TH1D("OnPeak_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *OnPeak_p_Dch = new TH1D("OnPeak_p_Dch",title_Dch,nbin, xmin,xmax);

  TH1D *OnPeak_s_D0 = new TH1D("OnPeak_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *OnPeak_s_Dch = new TH1D("OnPeak_s_Dch",title_Dch,nbin, xmin,xmax);

  TH1D *OnPeak_D0 = new TH1D("OnPeak_D0",title_D0,nbin, xmin,xmax);
  TH1D *OnPeak_Dch = new TH1D("OnPeak_Dch",title_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

  // Cuts on D mass:
  TCut * cut_D0Peak = new TCut("1.840 < DMass && DMass < 1.888");  //Peak
  TCut * cut_D0Side = new TCut("(1.816<DMass && DMass<1.840) || (1.888<DMass && DMass<1.912)");  //Sideband
  TCut * cut_DchPeak = new TCut("1.845 < DMass && DMass < 1.893");  //Peak
  TCut * cut_DchSide = new TCut("(1.821<DMass && DMass<1.845) || (1.893<DMass && DMass<1.917)");  //Sideband


  // Project trees
  treeMC->Project("OnPeak_p_D0",plotVal,*cut_D0&&*cut_D0Peak);
  treeMC->Project("OnPeak_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak);
  treeMC->Project("OnPeak_s_D0",plotVal,*cut_D0&&*cut_D0Side);
  treeMC->Project("OnPeak_s_Dch",plotVal,*cut_Dch&&*cut_DchSide);

  OnPeak_D0->Add(OnPeak_p_D0,OnPeak_s_D0,1,-1);
  OnPeak_Dch->Add(OnPeak_p_Dch,OnPeak_s_Dch,1,-1);



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

  OnPeak_D0->Write();
  OnPeak_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();

}