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

  TString plotVal("DLcosBY");
  TString title_D0("cosTheta_B-Dl (D0)");
  TString title_Dch("cosTheta_B-Dl (Dch)");

  // 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 = 40;
  double xmin = -10;
  double xmax = 5;


  TH1D *Dlnu_p_D0 = new TH1D("Dlnu_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *Dstar_p_D0 = new TH1D("Dstar_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *D0star_p_D0 = new TH1D("D0star_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *D1_p_D0 = new TH1D("D1_p_D0",title_D0,nbin, xmin,xmax); 
  TH1D *D1p_p_D0 = new TH1D("D1p_p_D0",title_D0,nbin, xmin,xmax);  
  TH1D *D2star_p_D0 = new TH1D("D2star_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *DPi_p_D0 = new TH1D("DPi_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *DstarPi_p_D0 = new TH1D("DstarPi_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *DiffB_p_D0 = new TH1D("DiffB_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *CascL_p_D0 = new TH1D("CascL_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *LMisID_p_D0 = new TH1D("LMisID_p_D0",title_D0,nbin, xmin,xmax);
  TH1D *Comb_p_D0 = new TH1D("Comb_p_D0",title_D0,nbin, xmin,xmax);

  TH1D *D0star_p_Dch = new TH1D("D0star_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *D1_p_Dch = new TH1D("D1_p_Dch",title_Dch,nbin, xmin,xmax); 
  TH1D *D1p_p_Dch = new TH1D("D1p_p_Dch",title_Dch,nbin, xmin,xmax);  
  TH1D *D2star_p_Dch = new TH1D("D2star_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DPi_p_Dch = new TH1D("DPi_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DstarPi_p_Dch = new TH1D("DstarPi_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DiffB_p_Dch = new TH1D("DiffB_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *CascL_p_Dch = new TH1D("CascL_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *LMisID_p_Dch = new TH1D("LMisID_p_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Comb_p_Dch = new TH1D("Comb_p_Dch",title_Dch,nbin, xmin,xmax);

  TH1D *Dlnu_s_D0 = new TH1D("Dlnu_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *Dstar_s_D0 = new TH1D("Dstar_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *D0star_s_D0 = new TH1D("D0star_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *D1_s_D0 = new TH1D("D1_s_D0",title_D0,nbin, xmin,xmax); 
  TH1D *D1p_s_D0 = new TH1D("D1p_s_D0",title_D0,nbin, xmin,xmax);  
  TH1D *D2star_s_D0 = new TH1D("D2star_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *DPi_s_D0 = new TH1D("DPi_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *DstarPi_s_D0 = new TH1D("DstarPi_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *DiffB_s_D0 = new TH1D("DiffB_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *CascL_s_D0 = new TH1D("CascL_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *LMisID_s_D0 = new TH1D("LMisID_s_D0",title_D0,nbin, xmin,xmax);
  TH1D *Comb_s_D0 = new TH1D("Comb_s_D0",title_D0,nbin, xmin,xmax);

  TH1D *D0star_s_Dch = new TH1D("D0star_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *D1_s_Dch = new TH1D("D1_s_Dch",title_Dch,nbin, xmin,xmax); 
  TH1D *D1p_s_Dch = new TH1D("D1p_s_Dch",title_Dch,nbin, xmin,xmax);  
  TH1D *D2star_s_Dch = new TH1D("D2star_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DPi_s_Dch = new TH1D("DPi_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DstarPi_s_Dch = new TH1D("DstarPi_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DiffB_s_Dch = new TH1D("DiffB_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *CascL_s_Dch = new TH1D("CascL_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *LMisID_s_Dch = new TH1D("LMisID_s_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Comb_s_Dch = new TH1D("Comb_s_Dch",title_Dch,nbin, xmin,xmax);

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

  TH1D *D0star_Dch = new TH1D("D0star_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *D1_Dch = new TH1D("D1_Dch",title_Dch,nbin, xmin,xmax); 
  TH1D *D1p_Dch = new TH1D("D1p_Dch",title_Dch,nbin, xmin,xmax);  
  TH1D *D2star_Dch = new TH1D("D2star_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DPi_Dch = new TH1D("DPi_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DstarPi_Dch = new TH1D("DstarPi_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *DiffB_Dch = new TH1D("DiffB_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *CascL_Dch = new TH1D("CascL_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *LMisID_Dch = new TH1D("LMisID_Dch",title_Dch,nbin, xmin,xmax);
  TH1D *Comb_Dch = new TH1D("Comb_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

  // 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_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_Dlnu);
  treeMC->Project("Dstar_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_Dstar);
  treeMC->Project("D0star_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_D0star);
  treeMC->Project("D1_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_D1);
  treeMC->Project("D1p_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_D1p);
  treeMC->Project("D2star_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_D2star);
  treeMC->Project("DPi_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_DPi);
  treeMC->Project("DstarPi_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_DstarPi);
  treeMC->Project("DiffB_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_DiffB);
  treeMC->Project("CascL_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_CascL);
  treeMC->Project("LMisID_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_LMisID);
  treeMC->Project("Comb_p_D0",plotVal,*cut_D0&&*cut_D0Peak&&*cut_Comb);

  treeMC->Project("D0star_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_D0star);
  treeMC->Project("D1_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_D1);
  treeMC->Project("D1p_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_D1p);
  treeMC->Project("D2star_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_D2star);
  treeMC->Project("DPi_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_DPi);
  treeMC->Project("DstarPi_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_DstarPi);
  treeMC->Project("DiffB_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_DiffB);
  treeMC->Project("CascL_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_CascL);
  treeMC->Project("LMisID_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_LMisID);
  treeMC->Project("Comb_p_Dch",plotVal,*cut_Dch&&*cut_DchPeak&&*cut_Comb);

  treeMC->Project("Dlnu_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_Dlnu);
  treeMC->Project("Dstar_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_Dstar);
  treeMC->Project("D0star_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_D0star);
  treeMC->Project("D1_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_D1);
  treeMC->Project("D1p_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_D1p);
  treeMC->Project("D2star_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_D2star);
  treeMC->Project("DPi_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_DPi);
  treeMC->Project("DstarPi_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_DstarPi);
  treeMC->Project("DiffB_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_DiffB);
  treeMC->Project("CascL_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_CascL);
  treeMC->Project("LMisID_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_LMisID);
  treeMC->Project("Comb_s_D0",plotVal,*cut_D0&&*cut_D0Side&&*cut_Comb);

  treeMC->Project("D0star_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_D0star);
  treeMC->Project("D1_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_D1);
  treeMC->Project("D1p_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_D1p);
  treeMC->Project("D2star_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_D2star);
  treeMC->Project("DPi_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_DPi);
  treeMC->Project("DstarPi_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_DstarPi);
  treeMC->Project("DiffB_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_DiffB);
  treeMC->Project("CascL_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_CascL);
  treeMC->Project("LMisID_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_LMisID);
  treeMC->Project("Comb_s_Dch",plotVal,*cut_Dch&&*cut_DchSide&&*cut_Comb);

  // sideband subtraction
  Dlnu_D0->Add(Dlnu_p_D0,Dlnu_s_D0,1,-1);
  Dstar_D0->Add(Dstar_p_D0,Dstar_s_D0,1,-1);
  D0star_D0->Add(D0star_p_D0,D0star_s_D0,1,-1);
  D1_D0->Add(D1_p_D0,D1_s_D0,1,-1);
  D1p_D0->Add(D1p_p_D0,D1p_s_D0,1,-1);
  D2star_D0->Add(D2star_p_D0,D2star_s_D0,1,-1);
  DPi_D0->Add(DPi_p_D0,DPi_s_D0,1,-1);
  DstarPi_D0->Add(DstarPi_p_D0,DstarPi_s_D0,1,-1);
  DiffB_D0->Add(DiffB_p_D0,DiffB_s_D0,1,-1);
  CascL_D0->Add(CascL_p_D0,CascL_s_D0,1,-1);
  LMisID_D0->Add(LMisID_p_D0,LMisID_s_D0,1,-1);
  Comb_D0->Add(Comb_p_D0,Comb_s_D0,1,-1);

  D0star_Dch->Add(D0star_p_Dch,D0star_s_Dch,1,-1);
  D1_Dch->Add(D1_p_Dch,D1_s_Dch,1,-1);
  D1p_Dch->Add(D1p_p_Dch,D1p_s_Dch,1,-1);
  D2star_Dch->Add(D2star_p_Dch,D2star_s_Dch,1,-1);
  DPi_Dch->Add(DPi_p_Dch,DPi_s_Dch,1,-1);
  DstarPi_Dch->Add(DstarPi_p_Dch,DstarPi_s_Dch,1,-1);
  DiffB_Dch->Add(DiffB_p_Dch,DiffB_s_Dch,1,-1);
  CascL_Dch->Add(CascL_p_Dch,CascL_s_Dch,1,-1);
  LMisID_Dch->Add(LMisID_p_Dch,LMisID_s_Dch,1,-1);
  Comb_Dch->Add(Comb_p_Dch,Comb_s_Dch,1,-1);


  // 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();

}