void calc_SBRatio_D0Kpi_lepMissID() {

#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/test");
  TString baseDir("/home1x/OtherMounts/hep03/khamano/R16b/thrustCut");
  TString nameBase_BpBm("SP-1235-BToDlnu-Run3-R16b");
  int nStart_BpBm = 1;
  int nStop_BpBm  = 218; //218 //87
  TString nameBase_B0B0bar("SP-1237-BToDlnu-Run3-R16b");
  int nStart_B0B0bar = 1;
  int nStop_B0B0bar  = 220; //220 //88
  TString nameBase_ccbar("SP-1005-BToDlnu-Run3-R16b");
  int nStart_ccbar = 1;
  int nStop_ccbar  = 208; //208 //83
  TString nameBase_uds("SP-998-BToDlnu-Run3-R16b");
  int nStart_uds = 1;
  int nStop_uds  = 143; //143 //58
  TString nameBase_tautau("SP-3429-BToDlnu-Run3-R16b");
  int nStart_tautau = 1;
  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("SBRatio_D0Kpi_lepMissID.txt");

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

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

  // 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];
  int lepType[200];
  float pD[200], pLep[200], cosBY[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("DlLepType", &lepType);

  int nSig_Elec(0), nBkg_Elec(0);
  int nLepMissID_NonB_Elec(0), nLepMissID_DmissReco_Elec(0); 
  int nLepMissID_BmissReco_Elec(0);
  int nNonB_Elec(0);
  int nNotDirectB_Elec(0), nNotFrmSameB_Elec(0);
  int nKMissID_Elec(0), nPiMissID_Elec(0), nKPiNotSameD_Elec(0);
  int nExtraKPi_Elec(0);

  int nSig_Mu(0), nBkg_Mu(0);
  int nLepMissID_NonB_Mu(0), nLepMissID_DmissReco_Mu(0); 
  int nLepMissID_BmissReco_Mu(0);
  int nNonB_Mu(0);
  int nNotDirectB_Mu(0), nNotFrmSameB_Mu(0);
  int nKMissID_Mu(0), nPiMissID_Mu(0), nKPiNotSameD_Mu(0);
  int nExtraKPi_Mu(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];
      int lepType1 = lepType[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;

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

	// Lepton flavour cut
	if (lepType1==1) {
	  // Sum up each signal mode and background
	  if (isABEvent1 == 0) {
	    nNonB_Elec++; 
	    if (BDecMod1 != 0) {
	      nLepMissID_NonB_Elec++;
	    }
	  }
	  if (isABEvent1 == 1) {
	    if (KMissID1 == 1) {
	      nKMissID_Elec++; 
	      if (BDecMod1 != 0) {
		nLepMissID_DmissReco_Elec++;
	      }
	    }
	    if (KMissID1 == 0) {
	      if (PiMissID1 == 1) {
		nPiMissID_Elec++; 
		if (BDecMod1 != 0) {
		  nLepMissID_DmissReco_Elec++;
		}
	      }
	      if (PiMissID1 == 0) {
		if (KPiSameD1 == 0) {
		  nKPiNotSameD_Elec++; 
		  if (BDecMod1 != 0) {
		    nLepMissID_DmissReco_Elec++;
		  }
		}
		if (KPiSameD1 == 1) {
		  if (noExtKPi1 == 0) {
		    nExtraKPi_Elec++;
		    if (BDecMod1 != 0) {
		      nLepMissID_DmissReco_Elec++;
		    }
		  }
		  if (noExtKPi1 == 1) {
		    if (BDecMod1 == 0) {
		      nLepMissID_BmissReco_Elec++; 
		    }
		    if (BDecMod1 != 0) {
		      if (isDirectFromB1 == 0) nNotDirectB_Elec++;
		      if (isDirectFromB1 == 1) {
			if (isFromSameB1 == 0) nNotFrmSameB_Elec++;
			if (isFromSameB1 == 1) {
			  nSig_Elec++;
			}
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}

 	// Lepton flavour cut
	if (lepType1==2) {
	  // Sum up each signal mode and background
	  if (isABEvent1 == 0) {
	    nNonB_Mu++; 
	    if (BDecMod1 != 0) {
	      nLepMissID_NonB_Mu++;
	    }
	  }
	  if (isABEvent1 == 1) {
	    if (KMissID1 == 1) {
	      nKMissID_Mu++; 
	      if (BDecMod1 != 0) {
		nLepMissID_DmissReco_Mu++;
	      }
	    }
	    if (KMissID1 == 0) {
	      if (PiMissID1 == 1) {
		nPiMissID_Mu++; 
		if (BDecMod1 != 0) {
		  nLepMissID_DmissReco_Mu++;
		}
	      }
	      if (PiMissID1 == 0) {
		if (KPiSameD1 == 0) {
		  nKPiNotSameD_Mu++; 
		  if (BDecMod1 != 0) {
		    nLepMissID_DmissReco_Mu++;
		  }
		}
		if (KPiSameD1 == 1) {
		  if (noExtKPi1 == 0) {
		    nExtraKPi_Mu++;
		    if (BDecMod1 != 0) {
		      nLepMissID_DmissReco_Mu++;
		    }
		  }
		  if (noExtKPi1 == 1) {
		    if (BDecMod1 == 0) {
		      nLepMissID_BmissReco_Mu++; 
		    }
		    if (BDecMod1 != 0) {
		      if (isDirectFromB1 == 0) nNotDirectB_Mu++;
		      if (isDirectFromB1 == 1) {
			if (isFromSameB1 == 0) nNotFrmSameB_Mu++;
			if (isFromSameB1 == 1) {
			  nSig_Mu++;
			}
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}
      }
    } //end of Dl loop
  } //end of event loop

  nBkg_Elec = nNonB_Elec + nLepMissID_BmissReco_Elec + nNotDirectB_Elec 
    + nNotFrmSameB_Elec + nKMissID_Elec + nPiMissID_Elec
    + nKPiNotSameD_Elec + nExtraKPi_Elec;

  nBkg_Mu = nNonB_Mu + nLepMissID_BmissReco_Mu + nNotDirectB_Mu 
    + nNotFrmSameB_Mu + nKMissID_Mu + nPiMissID_Mu
    + nKPiNotSameD_Mu + nExtraKPi_Mu;

  float deltaS2_Elec = nSig_Elec + nBkg_Elec + rho*rho*nBkg_Elec*nBkg_Elec;
  float Sig_Elec = nSig_Elec;
  float SBRatio_Elec(0);
  if (Sig_Elec>0) SBRatio_Elec = sqrt(deltaS2_Elec)/Sig_Elec;

  cout << "nSig_Elec = " << nSig_Elec << endl;
  cout << "nBkg_Elec = " << nBkg_Elec << endl;
  cout << "SBRatio_Elec = " << SBRatio_Elec << endl;

  float deltaS2_Mu = nSig_Mu + nBkg_Mu + rho*rho*nBkg_Mu*nBkg_Mu;
  float Sig_Mu = nSig_Mu;
  float SBRatio_Mu(0);
  if (Sig_Mu>0) SBRatio_Mu = sqrt(deltaS2_Mu)/Sig_Mu;

  cout << "nSig_Mu = " << nSig_Mu << endl;
  cout << "nBkg_Mu = " << nBkg_Mu << endl;
  cout << "SBRatio_Mu = " << SBRatio_Mu << 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 << "Electron" << endl;
  outFile << "number of events = " << nEvent << endl;
  outFile << "number of signal = " << nSig_Elec << endl;
  outFile << "number of non-B = " << nNonB_Elec << endl;
  outFile << "number of kaon miss ID = " << nKMissID_Elec << endl;
  outFile << "number of pi miss ID = " << nPiMissID_Elec << endl;
  outFile << "number of K and Pi not from same D = " 
	  << nKPiNotSameD_Elec << endl;
  outFile << "number of extra K or Pi from same D = " 
	  << nExtraKPi_Elec << endl;
  outFile << "number of lepton miss ID in non-B = " 
	  << nLepMissID_NonB_Elec << endl;
  outFile << "number of lepton miss ID in DmissReco = " 
	  << nLepMissID_DmissReco_Elec << endl;
  outFile << "number of lepton miss ID in BmissReco = " 
	  << nLepMissID_BmissReco_Elec << endl;
  outFile << "number of lepton not direct from B = " 
	  << nNotDirectB_Elec << endl;
  outFile << "number of lepton and D not from same B = " 
	  << nNotFrmSameB_Elec << endl;
  outFile << "number of background = " << nBkg_Elec << endl;
  outFile << "Delta_S / S = " << SBRatio_Elec << endl;
  outFile << "-----------------------------------" << endl;
  outFile << "Muon" << endl;
  outFile << "number of signal = " << nSig_Mu << endl;
  outFile << "number of non-B = " << nNonB_Mu << endl;
  outFile << "number of kaon miss ID = " << nKMissID_Mu << endl;
  outFile << "number of pi miss ID = " << nPiMissID_Mu << endl;
  outFile << "number of K and Pi not from same D = " 
	  << nKPiNotSameD_Mu << endl;
  outFile << "number of extra K or Pi from same D = " 
	  << nExtraKPi_Mu << endl;
  outFile << "number of lepton miss ID in non-B = " 
	  << nLepMissID_NonB_Mu << endl;
  outFile << "number of lepton miss ID in DmissReco = " 
	  << nLepMissID_DmissReco_Mu << endl;
  outFile << "number of lepton miss ID in BmissReco = " 
	  << nLepMissID_BmissReco_Mu << endl;
  outFile << "number of lepton not direct from B = " 
	  << nNotDirectB_Mu << endl;
  outFile << "number of lepton and D not from same B = " 
	  << nNotFrmSameB_Mu << endl;
  outFile << "number of background = " << nBkg_Mu << endl;
  outFile << "Delta_S / S = " << SBRatio_Mu << endl;
  outFile.close();


}