void histMk_1D_DMom_cc() {

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

  TString plotVal("DMom");
  TString title_D0("D momentum (D0)");
  TString title_Dch("D momentum (Dch)");

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

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


  TH1D *RealD_p_D0 = new TH1D("RealD_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *Cont_p_D0 = new TH1D("Cont_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *RealD_p_Dch = new TH1D("RealD_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Cont_p_Dch = new TH1D("Cont_p_Dch",title_Dch,nbin, xmin,xmax);

  TH1D *RealD_s_D0 = new TH1D("RealD_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *Cont_s_D0 = new TH1D("Cont_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *RealD_s_Dch = new TH1D("RealD_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Cont_s_Dch = new TH1D("Cont_s_Dch",title_Dch,nbin, xmin,xmax);

  TH1D *RealD_D0 = new TH1D("RealD_D0",title_D0,nbin, xmin,xmax);
  TH1D *Cont_D0 = new TH1D("Cont_D0",title_D0,nbin, xmin,xmax);
  TH1D *RealD_Dch = new TH1D("RealD_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Cont_Dch = new TH1D("Cont_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

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

  // RealD background
  TCut * cut_RealD = new TCut("noExtraKorPi==1 && mcKaonMissID==0 && mcPiMissID==0 && mcKPifromSameD==1");


  // Project trees
  treeMC->Project("RealD_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_RealD);
  treeMC->Project("Cont_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_Cont);
  treeMC->Project("RealD_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_RealD);
  treeMC->Project("Cont_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_Cont);

  treeMC->Project("RealD_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_RealD);
  treeMC->Project("Cont_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_Cont);
  treeMC->Project("RealD_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_RealD);
  treeMC->Project("Cont_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_Cont);

  // sideband subtraction
  RealD_D0->Add(RealD_p_D0,RealD_s_D0,1,-1);
  Cont_D0->Add(Cont_p_D0,Cont_s_D0,1,-1);
  RealD_Dch->Add(RealD_p_Dch,RealD_s_Dch,1,-1);
  Cont_Dch->Add(Cont_p_Dch,Cont_s_Dch,1,-1);


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

  RealD_D0->Write();
  Cont_D0->Write();

  RealD_Dch->Write();
  Cont_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();

}