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

  /** 
   * 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);
  }

TH1* getWeightedAverage(TH1* h1, TH1* h2, char* histName) {


  TH1* h = h1->Clone(histName);

  int NBinsMax = h1->GetNbinsX() +2;

  for(int i=0; i < NBinsMax; ++i) {

    double x1 = h1->GetBinContent(i);
    double e1 = h1->GetBinError(i);
    double x2 = h2->GetBinContent(i);
    double e2 = h2->GetBinError(i);

    double x = 0.0;
    double e = 0.0;

    if(x1 < 0.1) {x = x2; e = e2; }
    else if(x2 < 0.1) {x = x1; e = e1; }
    else {
      x = (x1/(e1*e1) + x2/(e2*e2)) / (1.0/(e1*e1) + 1.0/(e2*e2));
      e = 1.0 / sqrt(1.0/(e1*e1) + 1.0/(e2*e2));
    }

    h->SetBinContent(i, x);
    h->SetBinError(i, e);
  }

  return h;
}

TH1* GenerateConsistentHist(TH1* hbg, TH1* hprod)
{
    TH1* hpprod = 0;
    if (hbg!=NULL && hprod!=NULL)
    {
        hpprod = dynamic_cast<TH1*>(hbg->Clone(hprod->GetName()));
        if (hpprod!=NULL)
	{
            Reset(hpprod);
            Double_t binctr=0, x=0, ex=0;
//            cout << hpprod->GetName() << endl;
            for(Int_t i=1; i<=hpprod->GetNbinsX(); i++)
	    {
                binctr = hbg->GetBinCenter(i);
                x  = hprod->GetBinContent(hprod->GetXaxis()->FindBin(binctr));
                ex = hprod->GetBinError(hprod->GetXaxis()->FindBin(binctr));
                if (x!=0 && ex/x > 1.0001)
                    cout << "bin:" << i << "\tcontent:" << x << "\terr:" << ex
                            << " \terr/content:" << ex/x << endl;
                hpprod->SetBinContent(i,x);
                hpprod->SetBinError(i,ex);
	    }

	}
    }
    return hpprod;
}

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

TH1* subtractHistograms(const std::string& newHistogramName, const TH1* histogram_central, const TH1* histogram_shift, int mode)
{
  const TH1* histogramMinuend = 0;
  const TH1* histogramSubtrahend = 0;
  if ( compIntegral(histogram_central) >= compIntegral(histogram_shift) ) {
    histogramMinuend    = histogram_central;
    histogramSubtrahend = histogram_shift;
  } else {
    histogramMinuend    = histogram_shift;
    histogramSubtrahend = histogram_central;
  }
  if ( !(histogramMinuend && histogramSubtrahend) ) return 0;
  TH1* newHistogram = (TH1*)histogramMinuend->Clone(newHistogramName.data());
  newHistogram->Reset();
  if ( !newHistogram->GetSumw2N() ) newHistogram->Sumw2();
  int numBins = newHistogram->GetNbinsX();
  for ( int iBin = 0; iBin <= (numBins + 1); ++iBin ) {
    double newBinContent = histogramMinuend->GetBinContent(iBin) - histogramSubtrahend->GetBinContent(iBin);
    double newBinError2 = square(histogramMinuend->GetBinError(iBin)) + square(histogramSubtrahend->GetBinError(iBin));
    if ( mode == kRelative ) {
      if ( histogram_central->GetBinContent(iBin) > 0. ) {
	newBinContent /= histogram_central->GetBinContent(iBin);
	newBinError2 /= square(histogram_central->GetBinContent(iBin));
      } else {
	newBinContent = 0.;
	newBinError2 = 0.;
      }
    }
    newHistogram->SetBinContent(iBin, newBinContent);
    assert(newBinError2 >= 0.);
    newHistogram->SetBinError(iBin, TMath::Sqrt(newBinError2));
  }
  return newHistogram;
}

TH1 *
YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth)
{
  /* find highest edge in histo */
  Int_t binhi;
  Double_t hi;
  for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
    if (h->GetBinContent(ibin) != 0.) {
      binhi = ibin + 1;
      hi = h->GetBinLowEdge(ibin + 1);
      break;
    }
  }
  if(max<hi) {
	Printf("Warning! You should probably set a higher max value (Max = %f, hi = %f)", max, hi);
    return 0x0;
  }
  Int_t nbins = (max - hi) / binwidth;
 if(nbins<1)
	return 0x0;	 
  TH1 *hhi = new TH1F("hhi", "", nbins, hi, max);
  
  /* integrate function in histogram bins */
  Double_t cont, err, width;
  for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
    width = hhi->GetBinWidth(ibin + 1);
    cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
    err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
    hhi->SetBinContent(ibin + 1, cont / width);
    hhi->SetBinError(ibin + 1, err / width);
  }

  return hhi;
}

TH1 *
YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth)
{
  /* find lowest edge in histo */
  Int_t binlo;
  Double_t lo;
  for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
    if (h->GetBinContent(ibin) != 0.) {
      binlo = ibin;
      lo = h->GetBinLowEdge(ibin);
      break;
    }
  }
  
  Int_t nbins = (lo - min) / binwidth;
  if(nbins<1)
	return 0x0;		
  TH1 *hlo = new TH1F("hlo", "", nbins, min, lo);
  
  /* integrate function in histogram bins */
  Double_t cont, err, width;
  for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
    width = hlo->GetBinWidth(ibin + 1);
    cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
    err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
    hlo->SetBinContent(ibin + 1, cont / width);
    hlo->SetBinError(ibin + 1, err / width);
  }

  return hlo;
}

  //____________________________________________________________________
  void Run(const char* newName,        const char* oldName,
	   const char* newTitle="New", const char* oldTitle="Old")
  {
    TFile* newFile = TFile::Open(newName,"READ");
    TFile* oldFile = TFile::Open(oldName,"READ");
    if (!newFile || !oldFile) return;

    TH1* newCent = GetH1(newFile, "realCent");
    TH1* oldCent = GetH1(oldFile, "realCent");
    if (!newCent || !oldCent) return;

    TString  t; t.Form("#it{R}=#frac{%s}{%s}", newTitle, oldTitle);
    TCanvas* c     = new TCanvas("c", t, 1200, 800);
    c->SetTopMargin(0.01);
    c->SetRightMargin(0.20);
    fLegend = new TLegend(1-c->GetRightMargin(),
			  c->GetBottomMargin(),
			  1, 1-c->GetTopMargin(),
			  t);
    fLegend->SetFillStyle(0);
    fLegend->SetBorderSize(0);
    THStack* stack = new THStack("ratios","");
			       
    fMin = +1e6;
    fMax = -1e6;
    TH1* one = 0;
    for (Int_t i = newCent->GetNbinsX(); i--;) {
      Double_t c1 = newCent->GetXaxis()->GetBinLowEdge(i+1);
      Double_t c2 = newCent->GetXaxis()->GetBinUpEdge(i+1);
      Info("", "c1=%f c2=%f", c1, c2);
      TH1*     r  = One(newFile, oldFile, c1, c2);    
      if (!r) continue;
      if (!one) {
	one = static_cast<TH1*>(r->Clone("one"));
	one->SetDirectory(0);
	one->Reset();
	for (Int_t j = 1; j <= one->GetNbinsX(); j++) {
	  one->SetBinContent(j,1);
	  one->SetBinError  (j,0);
	}
      }
      // r->Add(one, i-1);
      // r->Scale(TMath::Power(10,i));
      stack->Add(r);
    }
    stack->Draw("nostack");
    stack->SetMinimum(0.95*fMin);
    stack->SetMaximum(1.05*fMax);
    stack->GetHistogram()->SetXTitle("#eta");
    stack->GetHistogram()->SetYTitle("#it{R}");
    fLegend->Draw();
    c->Modified();
    c->Update();
    c->cd();
    c->SaveAs(Form("%sover%s.png", newTitle, oldTitle));
  }  

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);
  }

}

makeMuonTriggerEfficiencyLUT()
{
  TH1* histoLUT = new TH2F("muonTriggerEfficiencyCorrection", "muonTriggerEfficiencyCorrection", 
			   1, -2.5, +2.5, 1, 0., 10000.);

  histoLUT->SetBinContent(1, 1, 0.98);
  histoLUT->SetBinError(1, 1, 0.02);

  TFile* outputFile = new TFile("../data/muonTriggerEfficiencyCorrection.root", "RECREATE");
  histoLUT->Write();
  delete outputFile;
}

TH1* compHistogramErr(const TH1* histogram)
{
  if ( !histogram ) return 0;
  std::string histogramErrName = Form("%sErr", histogram->GetName());
  TH1* histogramErr = (TH1*)histogram->Clone();
  int numBins = histogram->GetNbinsX();
  for ( int iBin = 1; iBin <= numBins; ++iBin ) {
    double binError = histogram->GetBinError(iBin);
    histogramErr->SetBinContent(iBin, binError);
    histogramErr->SetBinError(iBin, 0.);
  }
  return histogramErr;
}

TH1* divideHistogramByBinWidth(TH1* histogram)
{
  std::string histogramDensityName = Form("%s_density", histogram->GetName());
  TH1* histogramDensity = (TH1*)histogram->Clone(histogramDensityName.data());
  TAxis* xAxis = histogram->GetXaxis();
  int numBins = xAxis->GetNbins();
  for ( int iBin = 1; iBin <= numBins; ++iBin ) {
    double binContent = histogram->GetBinContent(iBin);
    double binError = histogram->GetBinError(iBin);
    double binWidth = xAxis->GetBinWidth(iBin);
    histogramDensity->SetBinContent(iBin, binContent/binWidth);
    histogramDensity->SetBinError(iBin, binError/binWidth);
  }
  return histogramDensity;
}

TH1 *computeEffVsCut(TH1 *discrShape, double total)
{
	TH1 *h = discrShape->Clone(Form("%s_effVsCut", discrShape->GetName()));
	h->Sumw2();
	h->SetMaximum(1.5);
	h->SetMinimum(1e-3);

	unsigned int n = h->GetNbinsX();
	for(unsigned int bin = 1; bin <= n; bin++) {
		double efficiency = h->Integral(bin, n + 1) / total;
		double error = sqrt(efficiency * (1 - efficiency) / total);
		h->SetBinContent(bin, efficiency);
		h->SetBinError(bin, error);
	}

	return h;
}

// A function that get histogram and sets contents to 0 
// if entries are too small
TH1 * getHisto(TFile * file, const char * name, unsigned int rebin) {
  TObject * h = file->Get(name);
  if(h == 0)
    throw edm::Exception(edm::errors::Configuration) 
      << "Can't find object " << name << "\n";
  TH1 * histo = dynamic_cast<TH1*>(h);
  if(histo == 0)
    throw edm::Exception(edm::errors::Configuration) 
      << "Object " << name << " is of type " << h->ClassName() << ", not TH1\n";
  histo->Rebin(rebin);  
  for(int i = 1; i <= histo->GetNbinsX(); ++i) {
    if(histo->GetBinContent(i) < 0.1) {
      histo->SetBinContent(i, 0.0);
      histo->SetBinError(i, 0.0);
    }
  }
  return histo;
}