C++ TH1::GetBinError方法代码示例

void dominik()
{
  TH1* matHistogramRoman = static_cast<TH1*>(extractObjectFromFile("lyRoman.root", "lightYieldProjectionY")->At(0));
  TList* objects = extractObjectFromFile("c.root", "chargeBins");
  TH1* matHistogramDominik = new TH1D("matHistogramDominik", ";channel;light yield / pixels", 512, -0.5, 512-0.5);
  int sipmIt = 0;
  for (int i = 0; i < objects->GetSize(); ++i) {
    TH1* h = static_cast<TH1*>(objects->At(i));
    if (h->GetLineColor() == 8) {
      for (int bin = 1; bin <= 128; ++bin) {
        matHistogramDominik->SetBinContent(512 - (sipmIt * 128 + bin - 1), h->GetBinContent(bin));
        if (h->GetBinError(bin) > 0)
          matHistogramDominik->SetBinError(512 - (sipmIt * 128 + bin - 1), h->GetBinError(bin));
      }
      ++sipmIt;
    }
  }
  TCanvas* c = new TCanvas;
  c->Divide(1, 2);
  c->cd(1);
  matHistogramDominik->Draw();
  matHistogramRoman->Draw("SAME");
  c->cd(2);
  TH1* h = static_cast<TH1*>(matHistogramDominik->Clone());
  h->Add(matHistogramRoman, -1);
  h->Draw();
}

// ===  FUNCTION  ============================================================
//         Name:  TPlot::GetEffHist
//  Description:  
// ===========================================================================
TH1* TPlot::GetEffHist(std::string Eff, std::string det, std::string algo)
{

  TH1* hNum = GetHist1D(listEff[Eff].first, det, algo);

  TH1* hDem = GetHist1D(listEff[Eff].second, det, algo);

  TH1* temp = (TH1*)hNum->Clone(Eff.c_str());

  assert(hNum->GetNbinsX() == hDem->GetNbinsX());
  for (int i = 0; i < hNum->GetNbinsX(); ++i)
  {
    double val = hNum->GetBinContent(i) / hDem->GetBinContent(i);
    double valerr = val * sqrt( pow(hNum->GetBinError(i)/hNum->GetBinContent(i), 2) +
        pow(hDem->GetBinError(i)/hDem->GetBinContent(i), 2) );
    if (isnan(val)) { val = 0; valerr = 0; }
    temp->SetBinContent(i, val);
    temp->SetBinError(i, valerr);
    std::cout << " bin " <<i <<" val " << val << std::endl;
  }

  //temp->Divide(hDem);
  temp->GetYaxis()->SetTitle("Efficiency");

  return temp;
}       // -----  end of function TPlot::GetEffHist  -----

void fit() {
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      fitmmp(h);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

void ScaleTo10pb( TH1& hist) {

  int bins = hist.GetNbinsX();
  double error = 0;

  for(int i=1; i<= bins; i++) {
    error = 8.0 * hist.GetBinError(i);
    // if( error > 1.0 ) error = 1.0;
    hist.SetBinError(i, error);
  }

}

RooHistN::RooHistN(const TH1 &data, Double_t nominalBinWidth, Double_t nSigma, RooAbsData::ErrorType etype, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the contents of the specified TH1 object
  // which may have fixed or variable bin widths. Error bars are
  // calculated using Poisson statistics. Prints a warning and rounds
  // any bins with non-integer contents. Use the optional parameter to
  // specify the confidence level in units of sigma to use for
  // calculating error bars. The nominal bin width specifies the
  // default used by addBin(), and is used to set the relative
  // normalization of bins with different widths. If not set, the
  // nominal bin width is calculated as range/nbins.

  initialize();
  // copy the input histogram's name and title
  SetName(data.GetName());
  SetTitle(data.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  // TH1::GetYaxis() is not const (why!?)
  setYAxisLabel(const_cast<TH1&>(data).GetYaxis()->GetTitle());
  
  // initialize our contents from the input histogram's contents
  Int_t nbin= data.GetNbinsX();
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data.GetBinCenter(bin);
    Stat_t y= data.GetBinContent(bin);
    Stat_t dy = data.GetBinError(bin) ;
    if (etype==RooAbsData::Poisson) {
      addBin(x,roundBin(y),data.GetBinWidth(bin),xErrorFrac);
    } else {
      addBinWithError(x,y,dy,dy,data.GetBinWidth(bin),xErrorFrac);
    }
  }
  // add over/underflow bins to our event count
  _entries+= data.GetBinContent(0) + data.GetBinContent(nbin+1);
}

void Template::makeResidualsControlPlot(const string& tag, unsigned int rebin)
/*****************************************************************/
{
    if(numberOfDimensions()>0 && m_template->GetNbinsX()%rebin!=0) return;
    if(numberOfDimensions()>1 && m_template->GetNbinsY()%rebin!=0) return;
    if(numberOfDimensions()>2 && m_template->GetNbinsZ()%rebin!=0) return;

    stringstream cpName, cpRawName, resMapName, resDistName, relErrDistName;
    cpName << m_name << "_cp";
    cpRawName << m_name << "_rawcp";
    resMapName << m_name << "_resmap_" << tag << "_rebin" << rebin;
    resDistName << m_name << "_resdist_" << tag << "_rebin" << rebin;
    relErrDistName << m_name << "_relerrdist_" << tag << "_rebin" << rebin;
    TH1* cpTmp = NULL;
    TH1* cpRawTmp = NULL;
    if(rebin==1)
    {
        cpTmp = dynamic_cast<TH1*>(m_template->Clone(cpName.str().c_str()));
        cpRawTmp = dynamic_cast<TH1*>(m_rawTemplate->Clone(cpRawName.str().c_str()));
    }
    else
    {
        if(numberOfDimensions()==2)
        {
            cpTmp = dynamic_cast<TH2F*>(m_template)->Rebin2D(rebin, rebin, cpName.str().c_str());
            cpRawTmp = dynamic_cast<TH2F*>(m_rawTemplate)->Rebin2D(rebin, rebin, cpName.str().c_str());
        }
        else if(numberOfDimensions()==3)
        {
            cpTmp = dynamic_cast<TH3F*>(m_template)->Rebin3D(rebin, rebin, rebin, cpName.str().c_str());
            cpRawTmp = dynamic_cast<TH3F*>(m_rawTemplate)->Rebin3D(rebin, rebin, rebin, cpName.str().c_str());
        }
    }
    TH1* resMap = dynamic_cast<TH1*>(cpTmp->Clone(resMapName.str().c_str()));
    TH1D* resDist = new TH1D(resDistName.str().c_str(), resDistName.str().c_str(), 30, -3, 3);
    resDist->StatOverflows();
    TH1D* relErrDist = new TH1D(relErrDistName.str().c_str(), relErrDistName.str().c_str(), 200, -1, 1);
    relErrDist->StatOverflows();
    unsigned int nbins1 = cpTmp->GetNbinsX();
    unsigned int nbins2 = cpTmp->GetNbinsY();
    unsigned int nbins3 = cpTmp->GetNbinsZ();
    for(unsigned int b1=1;b1<=nbins1;b1++)
    {
        for(unsigned int b2=1;b2<=nbins2;b2++)
        {
            for(unsigned int b3=1;b3<=nbins3;b3++)
            {
                double tmpValue = cpTmp->GetBinContent(b1,b2,b3);
                double tmpRawValue = cpRawTmp->GetBinContent(b1,b2,b3);
                double tmpRawError = cpRawTmp->GetBinError(b1,b2,b3);
                if(tmpValue>0. && tmpRawValue>0. && tmpRawError>0.)
                {
                    double res = (tmpRawValue-tmpValue)/tmpRawError;
                    double relErr = (tmpRawValue-tmpValue)/tmpRawValue;
                    if(numberOfDimensions()==2)
                    {
                        resMap->SetBinContent(b1,b2,b3, res);
                        resMap->SetBinError(b1,b2,b3, 0.);
                    }
                    resDist->Fill(res);
                    relErrDist->Fill(relErr);
                }
            }
        }
    }
    stringstream plotMapName, plotDistName, plotErrDistName;
    plotMapName << "control_" << getName() << "_resMap" << "_" << tag << "_rebin" << rebin;
    plotDistName << "control_" << getName() << "_resDist" << "_" << tag << "_rebin" << rebin;
    plotErrDistName << "control_" << getName() << "_relErrDist" << "_" << tag << "_rebin" << rebin;
    if(numberOfDimensions()==2)
    {
        TCanvas* c = new TCanvas(plotMapName.str().c_str(),plotMapName.str().c_str(), 700,700);
        resMap->SetContour(99);
        resMap->SetAxisRange(-3., 3., "z");
        resMap->Draw("color z");
        resMap->SetXTitle(getVariable(0).c_str());
        resMap->SetYTitle(getVariable(1).c_str());
        addControlPlot(c);
    }
    TCanvas* c2 = new TCanvas(plotDistName.str().c_str(),plotDistName.str().c_str(), 700,700);
    resDist->SetLineColor(kBlack);
    resDist->SetLineWidth(2);
    resDist->SetMarkerColor(kBlack);
    resDist->SetMarkerStyle(20);
    resDist->SetXTitle("(raw-template)/error_{raw}");
    resDist->Draw();
    addControlPlot(c2);

    TCanvas* c3 = new TCanvas(plotErrDistName.str().c_str(),plotErrDistName.str().c_str(), 700,700);
    relErrDist->SetLineColor(kBlack);
    relErrDist->SetLineWidth(2);
    relErrDist->SetMarkerColor(kBlack);
    relErrDist->SetMarkerStyle(20);
    relErrDist->SetXTitle("(raw-template)/raw");
    relErrDist->Draw();
    addControlPlot(c3);

    if(cpTmp) delete cpTmp;
    if(cpRawTmp) delete cpRawTmp;
}

void Fit(TString SIGNAL,TString HISTO,double scaleSGN)
{
  gROOT->ForceStyle();
  TFile *fDat = TFile::Open("Histo_flatTree_data_tmva"+SIGNAL+".root");
  TFile *fBkg = TFile::Open("Histo_flatTree_qcd_weights_tmva"+SIGNAL+".root");
  TFile *fSgn = TFile::Open("Histo_flatTree_"+SIGNAL+"_weights_tmva"+SIGNAL+".root");
  
  TH1 *hDat = (TH1*)fDat->Get(HISTO);
  TH1 *hBkgRaw = (TH1*)fBkg->Get(HISTO);
  TH1 *hSgn = (TH1*)fSgn->Get(HISTO);
  TH1 *hDat_JESlo = (TH1*)fDat->Get(HISTO+"_JESlo");
  TH1 *hBkgRaw_JESlo = (TH1*)fBkg->Get(HISTO+"_JESlo");
  TH1 *hSgn_JESlo = (TH1*)fSgn->Get(HISTO+"_JESlo");
  TH1 *hDat_JESup = (TH1*)fDat->Get(HISTO+"_JESup");
  TH1 *hBkgRaw_JESup = (TH1*)fBkg->Get(HISTO+"_JESup");
  TH1 *hSgn_JESup = (TH1*)fSgn->Get(HISTO+"_JESup");
  
  TH1F *hBkg = (TH1F*)hBkgRaw->Clone("Bkg");
  TH1F *hBkg_JESlo = (TH1F*)hBkgRaw_JESlo->Clone("Bkg_JESlo");
  TH1F *hBkg_JESup = (TH1F*)hBkgRaw_JESup->Clone("Bkg_JESup");
  
  hBkg->Smooth(2);
  hBkg_JESlo->Smooth(2);
  hBkg_JESup->Smooth(2);
  hSgn->Smooth(2);
  hSgn_JESlo->Smooth(2);
  hSgn_JESup->Smooth(2);
  
  double lumi = 4967;
  hBkg->Scale(lumi);
  hBkg_JESlo->Scale(lumi);
  hBkg_JESup->Scale(lumi);
  double k_factor = hDat->Integral()/hBkg->Integral();
  double k_factor_JESlo = hDat->Integral()/hBkg_JESlo->Integral();
  double k_factor_JESup = hDat->Integral()/hBkg_JESup->Integral();
  hBkg->Scale(k_factor);
  cout<<"Signal entries = "<<hSgn->GetEntries()<<endl;
  hSgn->Scale(lumi/scaleSGN);
  hBkg_JESlo->Scale(k_factor_JESlo);
  hSgn_JESlo->Scale(lumi/scaleSGN);
  hBkg_JESup->Scale(k_factor_JESup);
  hSgn_JESup->Scale(lumi/scaleSGN);
  hSgn_JESlo->Scale(hSgn->Integral()/hSgn_JESlo->Integral());
  hSgn_JESup->Scale(hSgn->Integral()/hSgn_JESup->Integral());
  
  TH1 *hBkg_STATlo = (TH1*)hBkg->Clone(HISTO+"_STATlo");
  TH1 *hSgn_STATlo = (TH1*)hSgn->Clone(HISTO+"_STATlo");
  TH1 *hBkg_STATup = (TH1*)hBkg->Clone(HISTO+"_STATup");
  TH1 *hSgn_STATup = (TH1*)hSgn->Clone(HISTO+"_STATup");
  
  float y1,e1;
  for(int i=0;i<hBkg->GetNbinsX();i++) {
    y1 = hBkg->GetBinContent(i+1);
    e1 = hBkg->GetBinError(i+1);
    hBkg_STATlo->SetBinContent(i+1,y1-e1);
    hBkg_STATup->SetBinContent(i+1,y1+e1);
    y1 = hSgn->GetBinContent(i+1);
    e1 = hSgn->GetBinError(i+1);
    hSgn_STATlo->SetBinContent(i+1,y1-e1);
    hSgn_STATup->SetBinContent(i+1,y1+e1);
  }
  hBkg_STATlo->Scale(hBkg->Integral()/hBkg_STATlo->Integral());
  hBkg_STATup->Scale(hBkg->Integral()/hBkg_STATup->Integral());
  hSgn_STATlo->Scale(hSgn->Integral()/hSgn_STATlo->Integral());
  hSgn_STATup->Scale(hSgn->Integral()/hSgn_STATup->Integral());
  
  double xMIN = hBkg->GetBinLowEdge(1);
  double xMAX = hBkg->GetBinLowEdge(hBkg->GetNbinsX()+1);
  double xMIN2 = hDat->GetBinLowEdge(hDat->FindFirstBinAbove(0.5));
  double xMAX2 = hDat->GetBinLowEdge(hDat->FindLastBinAbove(0.5)+1);
  RooRealVar x("x","x",xMIN2,xMAX2);
  
  RooDataHist data("data","dataset with x",x,hDat);
  RooDataHist bkg("qcd","bkg with x",x,hBkg);
  RooDataHist sgn("signal","sgn with x",x,hSgn);
  
  RooHistPdf bkgPDF("bkgPDF","bkgPDF",x,bkg,0);
  RooHistPdf sgnPDF("sgnPDF","sgnPDF",x,sgn,0);
  
  RooRealVar f("f","f",0,0.,1.);
  
  RooAddPdf model("model","model",RooArgList(sgnPDF,bkgPDF),RooArgList(f));
  
  RooFitResult* r = model.fitTo(data,Save());
  r->Print("v");
  double N = hDat->Integral();
  double B = hBkg->Integral();
  double S = hSgn->Integral();
  double m = f.getVal();
  double e = f.getError();
  cout<<"k-factor = "<<k_factor<<endl;
  cout<<N<<" "<<B<<" "<<S<<endl;
  cout<<"Total cross section =       "<<N/lumi<<" pb"<<endl;
  cout<<"Model cross section =       "<<S/lumi<<" pb"<<endl;
  cout<<"Fitted signal strength =    "<<m*N/S<<endl;
  cout<<"Fitted signal error =       "<<e*N/S<<endl;
  double p = 0.95;
  double xup = (N/S)*(m+sqrt(2.)*e*TMath::ErfInverse((1-p)*TMath::Erf(m/e)+p));
  cout<<"Bayesian Upper limit =      "<<xup<<endl;
  RooPlot* frame1 = x.frame();  
//.........这里部分代码省略.........

/** 
 * Make a GraphSysErr object
 * 
 * @param d  Directory
 * @param c1 Least centrality 
 * @param c2 Largest centrality 
 * 
 * @return Newly created GraphSysErr
 */
TObject* MakeGSE(TDirectory* d, Double_t c1, Double_t c2)
{
  if (!gROOT->GetClass("GraphSysErr")) return 0;
  TString  bin; bin.Form("%03dd%02d_%03dd%02d",
			 Int_t(c1), Int_t(c1*100)%100, 
			 Int_t(c2), Int_t(c2*100)%100);
  TString sub(bin); sub.Prepend("cent"); sub.Append("/dndeta");
  TString nme(bin); nme.Prepend("CENT_");
  TH1*    g = GetH1(d, sub);
  if (!g) return 0;
  
  Color_t  col     = g->GetMarkerColor();
  // Double_t bg      = (1-c1/100)*(2-0.1)+0.1;
  // Double_t c       = TMath::Power(c1/100,2)*(6.2-0.4)+0.4;
  Double_t bg      = CSysEval(c2, 0.02, 0.001);
  Double_t c       = CSysEval(c2, 0.005, 0.075);
  GraphSysErr* gse = new GraphSysErr(g->GetNbinsX());
  gse->SetName(nme);
  gse->SetTitle(Form("%5.1f - %5.1f%%", c1, c2));
  gse->SetKey("author", "PREGHENELLA : 20150");
  gse->SetKey("title", "dNch/deta in PbPb at 5023 GeV");
  gse->SetKey("obskey", "DN/DETARAP");
  gse->SetKey("reackey", "PB PB --> CHARGED X");
  gse->SetKey("laboratory", "CERN");
  gse->SetKey("accelerator", "LHC");
  gse->SetKey("detector", "TRACKLETS");
  gse->SetKey("reference", "ALICE-AN-2830");
  gse->AddQualifier("CENTRALITY IN PCT", Form("%.1f TO %.1f",c1,c2));
  gse->AddQualifier("SQRT(S)/NUCLEON IN GEV", "5023");
  gse->SetXTitle("ETARAP");
  gse->SetYTitle("DN/DETARAP");
  gse->SetMarkerStyle(g->GetMarkerStyle());
  gse->SetMarkerSize(g->GetMarkerSize());
  gse->SetDataOption(GraphSysErr::kNoTick);
  gse->SetMarkerColor(col);
  gse->SetLineColor(col);
  gse->SetFillColor(col);
  gse->SetSumFillColor(col);
  gse->SetSumLineColor(col);
  gse->SetSumOption(GraphSysErr::kBox);
  gse->SetCommonSumFillColor(col);
  gse->SetCommonSumLineColor(col);
  gse->SetCommonSumOption(GraphSysErr::kBox);
  MakeCommon(gse, "Particle composition",	0.01,   col);
  MakeCommon(gse, "Weak decay",			0.01,   col);
  MakeCommon(gse, "pT extrapolation",		0.02,   col);
  MakeCommon(gse, "EG dependence",		0.02,   col);
  MakeCommon(gse, "Background subrtaction",	bg,     col);
  MakeCommon(gse, "Centrality",			c,      col);
  Int_t acc = MakeP2P(gse, "Acceptance", col);

  Int_t j =  0;
  for (Int_t i = 1; i <= g->GetNbinsX(); i++) {
    Double_t eta  = g->GetXaxis()->GetBinCenter(i);
    Double_t eEta = g->GetXaxis()->GetBinWidth(i)/2;
    Double_t xo   = TMath::Abs(eta)+eEta;
    if (xo > 2) continue;
    Double_t ea  = 0.02*TMath::Power(xo/2,2);
    gse->SetPoint(j, eta, g->GetBinContent(i));
    gse->SetPointError(j, eEta, eEta);
    gse->SetStatError(j, g->GetBinError(i),g->GetBinError(i));
    gse->SetSysError(acc, j, eEta, eEta, ea/100, ea/100);
    j++;
  }
  return gse;
    
  
}

TList* readHistos(vector<string> inputRootFiles, bool fillNameStrings){
  TList* thelistHistos = new TList();
  for (int j=0; j< inputRootFiles.size(); j++) 
    {
      int numMax = 0;
      //cout <<"inputProcess: "<<inputProcess.at(j)<<endl;
      // Access input ROOT file (can access over secure shell)
      cout <<"inputRootFiles.at(j) "<<inputRootFiles.at(j)<<endl; 
      TFile *rootTFile = new TFile((inputRootFiles.at(j)).c_str()); // open root file
      TDirectory *current_sourcedir = gDirectory;
      TIter nextkey( current_sourcedir->GetListOfKeys() );
      TKey *key;
      
      // loop over keys(ROOT objects) within the root file
      while ((key = (TKey*)nextkey())) 
        {
	  TObject *obj = key->ReadObj();
	  
	  TH1 *histoObj = (TH1*)obj;
	  TH1* h1 = 0; //points to the new copied hist we will make
          int totalbins = 0;
	  // rename histogram
          if ( obj->IsA()->InheritsFrom( "TH1" ) ) 
            {
              string histname = obj->GetName();
	      size_t found = histname.find("lumi");
              string theProcess = "lumi";

              string new_theProcess = "_" + theProcess;
              histname.erase(found-1, new_theProcess.size());
              string histName = histname;
              cout << histName << endl;
              for (int o = 0; o < inputHistoName.size(); o++)
                {   
                  if ( histName == inputHistoName.at(o) )
		    {
                      int numBins = histoObj->GetXaxis()->GetNbins();
		      double lowerEdge = histoObj->GetXaxis()->GetXmin();
		      double upperEdge = histoObj->GetXaxis()->GetXmax();
		      
		      h1 = new TH1F (histName.c_str(), histName.c_str(),
			             numBins, lowerEdge, upperEdge);
                      
                      TH1F* tmp2 = static_cast<TH1F*>(h1);
                      totalbins = tmp2->GetSize();

                      for (int k=0; k <= totalbins; k++)
                        {
                          h1->SetBinContent(k,(histoObj->GetBinContent(k)));
                          h1->SetBinError(k,(histoObj->GetBinError(k)));
                        }

                      h1 = Rebinh(h1);
                      //cout<<"h1 "<<h1->GetName()<<"   Integral  "<<h1->Integral()<<endl; 
		      thelistHistos->Add(h1);
		      if(fillNameStrings) nHistList++;
		      if(fillNameStrings) theHistNameStrings.push_back(histName);
		      numMax++;
		    } // close if loop on selected histograms
		} // close for loop on selected histograms
	    } // close if loop InheritsFrom("TH1")
	} // close while loop
      // If first input ROOT file (Data), store total number of histograms
      if (j == 0)
          numHistos = numMax;
      // Else exit the code if a ROOT file has a different number of histograms
      else
        {
          if (numMax != numHistos)
            {
              cerr << "ERROR: Input Root Files DO NOT have the same number"
                   << " of histograms." << endl;
              exit (1);
            }
        }
    } // close for loop over root files

    return thelistHistos;
}

void getHistoInfo ()
{
 
  listHistos = new TList();    // Histograms normalized to one
  nHistList = 0;
  
  // loop over root files
  for (int j=0; j< inputRootFile.size(); j++) 
    {
      int numMax = 0;

      // Access input ROOT file (can access over secure shell)
      
      TFile *rootTFile = new TFile((inputRootFile.at(j)).c_str()); // open root file
      std::cout <<"inputRootFile.at(j)).c_str()  "<<inputRootFile.at(j).c_str()<<std::endl;
      TDirectory *current_sourcedir = gDirectory;
      TIter nextkey( current_sourcedir->GetListOfKeys() );
      TKey *key;
      
      // loop over keys(ROOT objects) within the root file
      while ((key = (TKey*)nextkey())) 
        {
	  TObject *obj = key->ReadObj();
	  
	  TH1 *histoObj = (TH1*)obj;
	  TH1* h1 = 0; //points to the new copied hist we will make
	  int totalbins = 0;
	  // rename histogram
          if ( obj->IsA()->InheritsFrom( "TH1" ) ) 
            {

              string histname = obj->GetName();
              
              size_t found = histname.find("lumi");
              string theProcess = "lumi";

              string new_theProcess = "_" + theProcess;
              histname.erase(found-1, new_theProcess.size());
              string histName = histname;
 
	      int numBins = histoObj->GetXaxis()->GetNbins();
	      double lowerEdge = histoObj->GetXaxis()->GetXmin();
	      double upperEdge = histoObj->GetXaxis()->GetXmax();

	      if ( obj->IsA()->InheritsFrom( "TH2" ) ) 
		{
		  int numBinsY = histoObj->GetYaxis()->GetNbins();
		  double lowerEdgeY = histoObj->GetYaxis()->GetXmin();
		  double upperEdgeY = histoObj->GetYaxis()->GetXmax();
		  h1 = new TH2F (histName.c_str(), histName.c_str(),
				 numBins, lowerEdge, upperEdge,
				 numBinsY, lowerEdgeY, upperEdgeY);
	      
		  //here (and below) GetSize() will allow us to
		  //iterate over 1 or 2d hists.  We pick up overflow
		  //bins too, but that's fine
		  TH2F* tmp1 = static_cast<TH2F*>(h1);
		  totalbins = tmp1->GetSize();
		} else {
	      h1 = new TH1F (histName.c_str(), histName.c_str(),
				   numBins, lowerEdge, upperEdge);
	      TH1F* tmp2 = static_cast<TH1F*>(h1);
	      totalbins = tmp2->GetSize();
	      }

	      // Apply scale factor?
	      if (inputApplyScaleFactor.at(j) == "True")
		{
		  for (int k=0; k <= totalbins; k++)
		    {
		      double scf = inputScaleFactor.at(j);
		      h1->SetBinContent(k,(histoObj->GetBinContent(k)*scf));
		      h1->SetBinError(k,(histoObj->GetBinError(k))*scf);
		    }
		} else
		{
		  for (int k=0; k <= totalbins; k++)
		    {
		      h1->SetBinContent(k,(histoObj->GetBinContent(k)));
		      h1->SetBinError(k,(histoObj->GetBinError(k)));
		    }
		}
              h1 = Rebinh(h1);
	      listHistos->Add(h1);
	      nHistList++;
	      theHistNameStrings.push_back(histName);
	      numMax++;
	    } // close if loop InheritsFrom("TH1")
	} // close while loop
      // If first input ROOT file (Data), store total number of histograms
      if (j == 0)
        numHistos = numMax;
      // Else exit the code if a ROOT file has a different number of histograms
      else
        {
          if (numMax != numHistos)
            {
              cerr << "ERROR: Input Root Files DO NOT have the same number"
                   << " of histograms." << endl;
              exit (1);
//.........这里部分代码省略.........

void
Example(Bool_t fit=true)
{
  // Load the class - if not already done 
  if (!gROOT->GetClass("GraphSysErr"))
    gROOT->LoadMacro("GraphSysErr.C+g");

  // Adjust size along X of common errors 
  gStyle->SetErrorX(.2);
  // Adjust size of hat, cap, ... - depends on canvas size!
  gStyle->SetEndErrorSize(10);

  // Make our object 
  GraphSysErr* gse = new GraphSysErr("foo", "Gaussian");
  // Draw data with-out ticks 
  gse->SetDataOption(GraphSysErr::kNoTick);
  gse->SetXTitle("X");
  gse->SetYTitle("Y");
  // Set some key/value pairs 
  gse->SetKey("laboratory", "The Center");
  gse->SetKey("accelerator", "Mega Collider");
  gse->SetKey("detector", "Huge Experiment");
  gse->SetKey("author", "Christensen");
  gse->SetKey("reference","Jour.All.Things A1,999");
  gse->SetKey("doi","9999-9999-9999-9999");
  gse->SetKey("abstract", "The data");
  gse->SetKey("location", "In the paper");
  gse->SetKey("reackey", "graviton -> tachyons");
  gse->SetKey("obskey", "GUT");
  // Adding qualifiers 
  gse->AddQualifier("question", "Life, universe, and everything");
  
  // Two sources of common errors one relative, one absolue
  UInt_t cm1 = gse->DefineCommon("Common 0.05", false, .05);
  UInt_t cm2 = gse->DefineCommon("Common 10%", true, .1);
  
  // Two sources of point-to-point errors, one relative, one absolute
  UInt_t pp1 = gse->DeclarePoint2Point("Point-to-Point 0.1-0.2", true);
  UInt_t pp2 = gse->DeclarePoint2Point("Point-to-Point 5-10%", false);
  
  // Set options on summed errors (in case of option COMBINED)
  gse->SetSumLineColor(kRed+2);
  gse->SetSumLineWidth(2);
  gse->SetSumTitle("All errors");
  gse->SetSumOption(GraphSysErr::kHat);
  
  // Set attributes of common errors 
  gse->SetSysFillColor(cm1, kRed+2);
  gse->SetSysFillStyle(cm1, 3001);
  gse->SetSysLineColor(cm1, kRed+2);
  gse->SetSysFillColor(cm2, kCyan+2);
  gse->SetSysFillStyle(cm2, 3001);
  gse->SetSysOption(cm1, GraphSysErr::kBox);
  gse->SetSysOption(cm2, GraphSysErr::kRect);
  
  // Set attributes of other errors 
  gse->SetSysLineColor(pp1, kBlue+2);
  gse->SetSysLineWidth(pp1, 2);
  gse->SetSysLineColor(pp2, kGreen+2);
  gse->SetSysLineWidth(pp2, 3);
  gse->SetSysOption(pp1, GraphSysErr::kBar);
  gse->SetSysOption(pp2, GraphSysErr::kHat);
  
  // Fill a histogram with a Guassian random deviate 
  TH1* h = new TH1F("h", "h", 30, -3, 3);
  h->Sumw2();
  h->SetDirectory(0);
  h->FillRandom("gaus",1000);
  h->Scale(1./1000, "width");
  
  // Fill in the data points 
  for (Int_t i = 0; i < h->GetNbinsX(); i++) { 
    Int_t    bin = i+1;
    Double_t x   = h->GetXaxis()->GetBinCenter(bin);
    Double_t y   = h->GetBinContent(bin);
    Double_t sta = h->GetBinError(bin);
    Double_t w   = h->GetXaxis()->GetBinWidth(bin);
    
    // Set data 
    gse->SetPoint(i, x, y);
    gse->SetPointError(i, w/2, w/2);
    gse->SetStatError(i, sta);
    
    // Set point-to-point errors 
    gse->SetSysError(pp1, i, 0., gRandom->Uniform(0.1, 0.2));
    gse->SetSysError(pp2, i, 0., 0., 
		     gRandom->Uniform(0.05, 0.1),
		     gRandom->Uniform(0.05, 0.1));
  } 
  // Remove temporary histogram
  delete h;

  // Build our canvas 
  TCanvas* c = new TCanvas("c","c", 1400, 1000);
  c->SetFillColor(0);
  c->SetFillStyle(0);
  c->SetTopMargin(0.01);
  c->SetRightMargin(0.01);
  
  // Draw or fit (and draw) a Guassian to the data
//.........这里部分代码省略.........

void DrawTwoInPad(TVirtualPad* p,
                  Int_t        sub,
                  TH1*         h1,
                  TH1*         h2,
		  Bool_t       ratio,
		  Bool_t       logy=false,
		  Bool_t       legend=false)
{
  TVirtualPad* pp = p->cd(sub);
  pp->SetRightMargin(0.02);
  pp->SetLeftMargin(0.10);
  TVirtualPad* ppp = pp;
  if (ratio) {
    pp->Divide(1,2,0,0);
    ppp = pp->cd(1);				
    ppp->SetRightMargin(0.02);
  }
  if (logy) ppp->SetLogy();
  TH1* hs[] = { h1, h2, 0 };
  if (h1->GetMaximum() < h2->GetMaximum()) {
    hs[0] = h2;
    hs[1] = h1;
  }
  TH1** ph = hs;
  Double_t size = (ratio ? 0.1 : 0.05);
  Double_t off  = (ratio ? 0.6 : 0.5);
  h1->SetFillStyle(3004);
  h2->SetFillStyle(3005);
  while (*ph) {
    TString opt("hist");
    if (ph != hs) opt.Append(" same");

    TH1* copy = (*ph)->DrawCopy(opt);
    copy->GetXaxis()->SetLabelSize(2*size);
    copy->GetYaxis()->SetLabelSize(size);
    copy->GetYaxis()->SetTitleSize(size);
    copy->GetYaxis()->SetTitleOffset(off);
    copy->SetYTitle(copy->GetTitle());
    copy->SetTitle("");
    copy->SetDirectory(0);
    ph++;
  }
  TString s1 = h1->GetYaxis()->GetTitle();
  TString s2 = h2->GetYaxis()->GetTitle();
  
  if (legend) { 
    TLegend* l = new TLegend(0.6, 0.1, 0.9, 0.9);
    l->SetBorderSize(0);
    TLegendEntry* e = l->AddEntry("dummy", s1, "lf");
    l->SetFillColor(kWhite);
    e->SetFillColor(kBlack);
    e->SetFillStyle(h1->GetFillStyle());
    e = l->AddEntry("dummy", s2, "lf");
    e->SetFillColor(kBlack);
    e->SetFillStyle(h2->GetFillStyle());
    l->Draw();
  }
  if (!ratio) return;
  ppp = pp->cd(2);
  ppp->SetRightMargin(0.02);
  TH1* r = static_cast<TH1*>(h1->Clone(Form("ratio%s", h1->GetName())));
  r->SetDirectory(0);
  r->SetTitle("");
  r->GetXaxis()->SetLabelSize(size);
  r->GetYaxis()->SetLabelSize(size);
  r->GetYaxis()->SetTitleSize(0.9*size);
  r->GetYaxis()->SetTitleOffset(0.9*off);
  r->SetMarkerStyle(20);
  r->SetMarkerColor(h1->GetFillColor()+1);
  r->SetFillStyle(3007);
  r->SetYTitle(Form("#frac{%s}{%s}", s1.Data(), s2.Data()));

  // r->Add(h2, -1);
  // r->Divide(h1);
  if (!r->IsA()->InheritsFrom(TProfile::Class())) {
    r->GetSumw2()->Set(0); // r->Sumw2(false);
    h2->GetSumw2()->Set(0); // h2->Sumw2(false);
  }
  r->Divide(h2);
  Printf("%s", r->GetName());
  for (UShort_t bin = 1; bin <= r->GetNbinsX(); bin++) {
    Printf("  bin # %2d: Diff=%g+/-%g", bin, r->GetBinContent(bin),
	   r->GetBinError(bin));
    r->SetBinError(bin, 0);
  }
  r->GetSumw2()->Set(0); //r->Sumw2(false);
  r->SetMarkerSize(4);
  r->SetMaximum(r->GetMaximum()*1.2);
  r->SetMinimum(r->GetMinimum()*0.8);
  r->Draw("hist text30");  
  p->Modified();
  p->Update();
  p->cd();
}

void makeStack(TString myVar, TString myCut, TString myName, TString myAxisNameX, TString myAxisNameY, 
               vector<const Sample*>& listOfSignals, vector<const Sample*>& listOfSamples, vector<const Sample*> listOfDatasets, 
               TString inFileName,
               bool isBlind, bool isLog, bool drawSignal, bool drawLegend,
               int nBins, float xLow, float xHigh,
               float* xlowVec)
{
  // prepare the input file
  TFile* infile = new TFile(inFileName, "READ"); 
  infile -> cd();
  
  // prepare the stack
  THStack *hs = new THStack("hs","");
  // prepare the histos pointers
  TH1F*   hist[20];
  // prepare the tree pointers
  TTree*  tree[20];
  // prepare the legend
  TLegend* leg = new TLegend(.7485,.7225,.9597,.9604);
  leg->SetFillColor(0);
  // prepare the colors
  Int_t col[20] = {46,2,12,5,3,4,9,7,47,49,49,50,51,52,53,54,55,56,57,58};
  // prepare the cut
  if (isBlind) myCut += "*(phoMetDeltaPhi < 2.9)";        
  // prepare the Y axis lable
  if (xlowVec != 0) myAxisNameY = "Events/" + myAxisNameY;
  else {
    float binWidth = (xHigh-xLow)/nBins;
    TString tempString;
    tempString.Form("%.2f ",binWidth); 
    myAxisNameY = "Events/" + tempString + myAxisNameY;
  }
  // prepare the legend strings
  vector<TString> theLegends;
  
  // loop through the datasets and produce the plots
  TH1F* hdata;
  TH1F* hsignal;
  //prepare data and signal histos
  if (xlowVec != 0) hdata   = new TH1F("hdata","",nBins,xlowVec);
  else hdata = new TH1F("hdata","",nBins,xLow,xHigh);
  if (xlowVec != 0) hsignal = new TH1F("hsignal","",nBins,xlowVec);
  else hsignal = new TH1F("hsignal","",nBins,xLow,xHigh);

  TTree*  treedata[20];
  for (UInt_t iDatas=0; iDatas < listOfDatasets.size(); iDatas++) {
    //get the tree
    treedata[iDatas] = (TTree*) infile -> Get(listOfDatasets.at(iDatas)->Name()->Data());

    //fill the histogram
    if ( iDatas == 0 ) treedata[iDatas] -> Draw(myVar + " >> hdata","evt_weight*kf_weight*pu_weight" + myCut);
    else treedata[iDatas] -> Draw(myVar + " >>+ hdata","evt_weight*kf_weight*pu_weight" + myCut);
    
    if ( isBlind && iDatas == 0 ) leg -> AddEntry(hdata, "DATA (19.8 fb^{-1})", "pl");    
    
  }//end loop on datasets
  if (xlowVec != 0) {
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinError  (iBin,hdata->GetBinError(iBin)/hdata->GetBinWidth(iBin));
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinContent(iBin,hdata->GetBinContent(iBin)/hdata->GetBinWidth(iBin));
  }

  TTree*  treesignal[20];
  for (UInt_t iSignal=0; iSignal < listOfSignals.size(); iSignal++) {
    //get the tree
    treesignal[iSignal] = (TTree*) infile -> Get(listOfSignals.at(iSignal)->Name()->Data());

    //fill the histogram
    TString thisScale = Form("%f *", *(listOfSignals.at(iSignal)->Scale()));
    if ( iSignal == 0 ) treesignal[iSignal] -> Draw(myVar + " >> hsignal",thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    else treesignal[iSignal] -> Draw(myVar + " >>+ hsignal",thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    
    if ( drawSignal && iSignal == 0 ) leg -> AddEntry(hsignal, "Signal", "l");    
    
  }//end loop on signals
  if (xlowVec != 0) {
    for (int iBin = 1; iBin <= nBins; iBin++) hsignal->SetBinError  (iBin,hsignal->GetBinError(iBin)/hsignal->GetBinWidth(iBin));
    for (int iBin = 1; iBin <= nBins; iBin++) hsignal->SetBinContent(iBin,hsignal->GetBinContent(iBin)/hsignal->GetBinWidth(iBin));
  }
  hsignal -> SetLineColor(49);
  hsignal -> SetLineWidth(4.0);
       
  int theHistCounter = 0;
  // loop through the samples and produce the plots
  for (UInt_t iSample=0; iSample < listOfSamples.size(); iSample++) {

    //determine if the histo is first of the series
    bool isFirstOfSerie = (*listOfSamples.at(iSample)->Legend()).CompareTo(" ");
    bool isLastOfSerie = false;
    if (iSample == listOfSamples.size() - 1) isLastOfSerie = true;
    if (iSample < listOfSamples.size() - 1 && (*listOfSamples.at(iSample+1)->Legend()).CompareTo(" ") != 0) isLastOfSerie = true;
    
    //get the tree
    tree[iSample] = (TTree*) infile -> Get(listOfSamples.at(iSample)->Name()->Data());
    //if sample first of the list create a new histogram
    if (isFirstOfSerie) {
       TString thisHistName = "h_" + *(listOfSamples.at(iSample)->Name());
       //variable bin histo
       if (xlowVec != 0) hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xlowVec);
       //fixed bin histo
       else hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xLow,xHigh);
//.........这里部分代码省略.........

void Process(TString fname, TString oldFolder, int toMass, int fromMass)
{




  std::string channels[] = {"data_obs", "WH", "TT",  "WjLF", "WjHF", "ZjLF", "ZjHF" , "VV" , "s_Top"};
  std::string systs[] = {"eff_b", "fake_b", "res_j", "scale_j" , "stat" };

  kount = 0;  

  gROOT->SetStyle("Plain");
  setTDRStyle();

  TFile * file = new TFile(fname.Data(), "READ");
  std::cout << "reading " << fname.Data() << std::endl;
  TString outname(massS[toMass]);
  outname.Append("_June8.root");
  fname.ReplaceAll(".root",outname.Data());
 
///BEGIN CHECK RBIN
  int addRightBin = 0, addRightBinm1 = 0 , rightBinNo = 0;
  float a,overflow;
  RooWorkspace *mytempWS =  (RooWorkspace*) file->Get(oldFolder.Data());
  RooRealVar BDT("CMS_vhbb_BDT_Wln", "BDT", -1, 1);

  float Signal = 0;
  float Background = 0;
  float Backgroundm1 = 0;
  float Data = 0;

  RooDataHist* tempRooDataHistNomS = (RooDataHist*)  mytempWS->data(channels[1].c_str());
  TH1 *tempHistNomS = tempRooDataHistNomS->createHistogram(channels[1].c_str(),BDT,RooFit::Binning(bins));
  tempHistNomS->Rebin(rebin);


  TH1 *tempHistNomD = tempRooDataHistNomS->createHistogram(channels[0].c_str(),BDT,RooFit::Binning(bins));
  tempHistNomD->Rebin(rebin);



  for(int i = 1; i <= tempHistNomS->GetNbinsX(); i++)
  { 
    if(tempHistNomS->GetBinContent(i) > 0)
      { rightBinNo = i; Signal = tempHistNomS->GetBinContent(i); Data = tempHistNomS->GetBinError(i); }
  }

  for(int i = 2; i < 9; i++)
  {
   RooDataHist* tempRooDataHistNom = (RooDataHist*)  mytempWS->data(channels[i].c_str());
   TH1 *tempHistNom = tempRooDataHistNom->createHistogram(channels[i].c_str(),BDT,RooFit::Binning(bins));
   tempHistNom->Rebin(rebin);
   Background += tempHistNom->GetBinContent(rightBinNo);
   

   //   if(tempHistNom->GetBinContent(rightBinNo-1) == 0)
   // {
   // Backgroundm1 = 0;  
   //std::cout << "ARGHGHGHGHGH bkg is 0 still at left" << std::endl;
   //   std::cin >> ;
   // }

   Backgroundm1+=  tempHistNom->GetBinContent(rightBinNo-1); 

   overflow = tempHistNom->GetBinContent(rightBinNo+1) ;
   
   if(tempHistNom->GetBinContent(rightBinNo+1) > 0)
   {
   std::cout << "ARGHGHGHGHGH overflow at right" << std::endl;
   //   std::cin >> ;
   }
   a+= overflow;
  }

  if( (Background ==0) ) addRightBin = 1;
else  addRightBin = 0;


std::cout << "################# folder" << oldFolder << std::endl;
 std::cout<< "################# CHECK RBIN:: right bin n  " << rightBinNo << " signal: " << Signal << " bkg: " << Background << " bkgm1: " <<  Backgroundm1  << " Data:  " << Data <<  " at right there is an overflow of: "<< a << std::endl;
std::cout << "########################### CHECK RBIN:: REBINNING: " << addRightBin << std::endl;  
 if ( (Backgroundm1 == 0)  ) 
   { addRightBinm1 =1 ;
     std::cout << "ARGHGHGHGHGH " << Backgroundm1 << " at left" << std::endl;
    
   }

std::cout << "########################### need to rebin further? " << addRightBinm1  << std::endl;  
///END CHECK RBIN



  TFile * outfile = new TFile(fname.Data(), "RECREATE");
  using namespace RooFit;
  RooWorkspace *myWS = new RooWorkspace(oldFolder.Data(),oldFolder.Data());
  myWS->factory("CMS_vhbb_BDT_Wln[-1.,1.]"); ///NEW VARIABLE NAME HERE 

  TString oldFolder2(oldFolder.Data());
  oldFolder2.Append("2");
  RooWorkspace *myWS2 = new RooWorkspace(oldFolder2.Data(),oldFolder2.Data());
//.........这里部分代码省略.........