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

TH1* getHistogram(TFile* inputFile, const std::string& channel, double massPoint, 
		  const std::string& directory, const std::string& histogramName, double metResolution)
{
  std::string process = "";
  if   ( massPoint <  95. ) process = "ZToTauTau";
  else                      process = Form("HToTauTau_M-%1.0f", massPoint);
  
  std::string metResolution_label = "";
  if ( metResolution > 0. ) metResolution_label = Form("pfMEtRes%1.0f", metResolution);
  else metResolution_label = "pfMEtResMC";

  std::vector<TH1*> histograms;
  std::string directory_process = 
    //Form("DQMData/%s/%s/%s/%s/plotEntryType1", process.data(), channel.data(), metResolution_label.data(), directory.data());
    Form("DQMData/%s/%s/%s/plotEntryType1", process.data(), channel.data(), directory.data());
  histograms.push_back(getHistogram(inputFile, directory_process, histogramName));
  TH1* histogramSum = NULL;
  for ( std::vector<TH1*>::const_iterator histogram = histograms.begin();
	histogram != histograms.end(); ++histogram ) {
    if ( !histogramSum ) {
      std::string histogramSumName = std::string((*histogram)->GetName()).append("_summed");
      histogramSum = (TH1*)(*histogram)->Clone(histogramSumName.data());
    } else {
      histogramSum->Add(*histogram);
    }
  }

  assert(histogramSum);

  if ( !histogramSum->GetSumw2N() ) histogramSum->Sumw2();
  if ( histogramSum->Integral() > 0. ) histogramSum->Scale(1./histogramSum->Integral());
  
  return histogramSum;
}

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* loadHistogram(TFile* inputFile, const std::string& directory, const std::string& histogramName)
{
  std::string histogramName_full = Form("%s/%s", directory.data(), histogramName.data());
  TH1* histogram = dynamic_cast<TH1*>(inputFile->Get(histogramName_full.data()));
  if ( !histogram ) {
    std::cerr << "Failed to load histogram = " << histogramName_full << " from file = " << inputFile->GetName() << " !!" << std::endl;
    assert(0);
  }
  if ( !histogram->GetSumw2N() ) histogram->Sumw2();
  histogram->Rebin(4);
  return histogram;
}

TH1* loadHistogram(TFile* inputFile, const std::string& histogramName)
{
  std::cout << "loading histogram = " << histogramName << " from file = " << inputFile->GetName() << std::endl;
  TH1* histogram = dynamic_cast<TH1*>(inputFile->Get(histogramName.data()));
  std::cout << "histogram = " << histogram << std::endl;
  if ( !histogram ) {
    std::cerr << "Failed to load histogram = " << histogramName << " from file = " << inputFile->GetName() << " --> skipping !!" << std::endl;
    return 0;
  }
  if ( !histogram->GetSumw2N() ) histogram->Sumw2();
  if ( dynamic_cast<TH2*>(histogram) ) histogram = linearizeHistogram(histogram);
  //else if ( histogram->GetNbinsX() >= 100 ) histogram->Rebin(5);
  return histogram;
}

TH1* getHistogram(TFile* inputFile, const TString& directory, const TString& histogramName)
{  
  TString histogramName_full = TString("DQMData").Append("/").Append(directory);
  if ( !histogramName_full.EndsWith("/") ) histogramName_full.Append("/");
  histogramName_full.Append(histogramName);

  TH1* histogram = (TH1*)inputFile->Get(histogramName_full.Data());

  if ( histogram && !histogram->GetSumw2N() ) histogram->Sumw2();
  else if ( !histogram) 
    std::cerr << "Failed to load histogram = " << histogramName_full << " from file = " << inputFile->GetName() << " !!" << std::endl;

  return histogram;
}

TH1* getHistogram(TFile* inputFile, const TString& dqmDirectory, const TString& meName)
{  
  TString histogramName = dqmDirectory;
  if ( histogramName.Length() > 0 && !histogramName.EndsWith("/") ) histogramName.Append("/");
  histogramName.Append(meName);

  TH1* histogram = (TH1*)inputFile->Get(histogramName.Data());
  std::cout << "histogramName = " << histogramName.Data() << ": histogram = " << histogram;
  if ( histogram ) std::cout << ", integral = " << histogram->Integral();
  std::cout << std::endl; 

  if ( !histogram->GetSumw2N() ) histogram->Sumw2();

  if ( histogram->GetDimension() == 1 ) histogram->Rebin(5);

  return histogram;
}

TH1* loadHistogram(TFile* inputFile, const std::string& directory, const std::string& histogramName)
{
  std::string histogramName_full = Form("%s/%s", directory.data(), histogramName.data());
  TH1* histogram = dynamic_cast<TH1*>(inputFile->Get(histogramName_full.data()));
  if ( !histogram ) {
    std::cerr << "Failed to load histogram = " << histogramName_full << " from file = " << inputFile->GetName() << " !!" << std::endl;
    assert(0);
  }
  if ( !histogram->GetSumw2N() ) histogram->Sumw2();

  //int numBins = histogram->GetNbinsX();
  //for ( int iBin = 0; iBin <= (numBins + 1); ++iBin ) {
  //  histogram->SetBinError(iBin, 0.);
  //}
  //histogram->SetBinContent(0, 0.);
  //histogram->SetBinContent(numBins + 1, 0.);
  return histogram;
}

//.........这里部分代码省略.........
  histogramRares_density->SetFillColor(kBlue - 8); 
  legend->AddEntry(histogramRares_density, "Rares", "f");

  THStack* histogramStack_density = new THStack("stack", "");
  histogramStack_density->Add(histogramRares_density);
  histogramStack_density->Add(histogramEWK_density);
  histogramStack_density->Add(histogramTTV_density);
  histogramStack_density->Add(histogramTT_density);
  histogramStack_density->Draw("histsame");
  
  histogramBgrUncertainty_density->SetFillColor(kBlack);
  histogramBgrUncertainty_density->SetFillStyle(3344);    
  histogramBgrUncertainty_density->Draw("e2same");
  legend->AddEntry(histogramBgrUncertainty_density, "Uncertainty", "f");

  histogramTTH_density->SetLineWidth(2);
  histogramTTH_density->SetLineStyle(1);
  histogramTTH_density->SetLineColor(kBlue);
  histogramTTH_density->Draw("histsame");
  
  histogramData_density->Draw("epsame");
  histogramData_density->Draw("axissame");
  
  legend->Draw();
  
  addLabel_CMS_luminosity(0.2050, 0.9225, 0.6850);
  
  canvas->cd();
  bottomPad->Draw();
  bottomPad->cd();
  
  TH1* histogramRatio = (TH1*)histogramData->Clone("histogramRatio");
  histogramRatio->Reset();
  if ( !histogramRatio->GetSumw2N() ) histogramRatio->Sumw2();
  checkCompatibleBinning(histogramRatio, histogramBgrSum);
  histogramRatio->Divide(histogramData, histogramBgrSum);
  int numBins_bottom = histogramRatio->GetNbinsX();
  for ( int iBin = 1; iBin <= numBins_bottom; ++iBin ) {
    double binContent = histogramRatio->GetBinContent(iBin);
    if ( histogramData && histogramData->GetBinContent(iBin) >= 0. ) histogramRatio->SetBinContent(iBin, binContent - 1.0);
    else histogramRatio->SetBinContent(iBin, -10.);
  }
  histogramRatio->SetTitle("");
  histogramRatio->SetStats(false);
  histogramRatio->SetMinimum(-0.50);
  histogramRatio->SetMaximum(+0.50);
  histogramRatio->SetMarkerStyle(histogramData_density->GetMarkerStyle());
  histogramRatio->SetMarkerSize(histogramData_density->GetMarkerSize());
  histogramRatio->SetMarkerColor(histogramData_density->GetMarkerColor());
  histogramRatio->SetLineColor(histogramData_density->GetLineColor());
  histogramRatio->Draw("ep");
  
  TAxis* xAxis_bottom = histogramRatio->GetXaxis();
  xAxis_bottom->SetTitle(xAxis_top->GetTitle());
  xAxis_bottom->SetLabelColor(1);
  xAxis_bottom->SetTitleColor(1);
  xAxis_bottom->SetTitleOffset(1.20);
  xAxis_bottom->SetTitleSize(0.13);
  xAxis_bottom->SetLabelOffset(0.02);
  xAxis_bottom->SetLabelSize(0.10);
  xAxis_bottom->SetTickLength(0.055);
  
  TAxis* yAxis_bottom = histogramRatio->GetYaxis();
  yAxis_bottom->SetTitle("#frac{Data - Simulation}{Simulation}");
  yAxis_bottom->SetTitleOffset(0.80);
  yAxis_bottom->SetNdivisions(505);

//.........这里部分代码省略.........
    legend1->AddEntry(histogram_ttW_density, "t#bar{t}W", "f");
    legend1->Draw();
    TLegend* legend2 = new TLegend(0.6600, legend_y0, 0.9350, 0.9250, NULL, "brNDC");
    legend2->SetFillStyle(0);
    legend2->SetBorderSize(0);
    legend2->SetFillColor(10);
        legend2->SetTextSize(0.050); 
    legend2->AddEntry(histogram_EWK_density, "Electroweak", "f");
    legend2->AddEntry(histogram_Rares_density, "Rares", "f");
    legend2->AddEntry(histogram_fakes_density, "Fakes", "f");    
    if ( histogramErr_mc ) legend2->AddEntry(histogramErr_mc_density, "Uncertainty", "f");
    legend2->Draw();
  }

  //addLabel_CMS_luminosity(0.2100, 0.9700, 0.6350);
  addLabel_CMS_preliminary(0.2100, 0.9700, 0.6350);

  TPaveText* label_category = 0;
  if ( showLegend ) label_category = new TPaveText(0.6600, legend_y0 - 0.0550, 0.9350, legend_y0, "NDC");
  else label_category = new TPaveText(0.2350, 0.8500, 0.5150, 0.9100, "NDC");
  label_category->SetTextAlign(13);
  label_category->AddText(label.data());
  label_category->SetTextSize(0.055);
  label_category->SetTextColor(1);
  label_category->SetFillStyle(0);
  label_category->SetBorderSize(0);
  label_category->Draw();

  canvas->cd();
  bottomPad->Draw();
  bottomPad->cd();
 
  TH1* histogramRatio = (TH1*)histogram_data_density->Clone("histogramRatio");
  if ( !histogramRatio->GetSumw2N() ) histogramRatio->Sumw2();
  histogramRatio->SetTitle("");
  histogramRatio->SetStats(false);
  histogramRatio->SetMinimum(-0.99);
  histogramRatio->SetMaximum(+0.99);
  histogramRatio->SetMarkerColor(histogram_data_density->GetMarkerColor());
  histogramRatio->SetMarkerStyle(histogram_data_density->GetMarkerStyle());
  histogramRatio->SetMarkerSize(histogram_data_density->GetMarkerSize());
  histogramRatio->SetLineColor(histogram_data_density->GetLineColor());

  TH1* histogramRatioUncertainty = (TH1*)histogram_data_density->Clone("histogramRatioUncertainty");
  if ( !histogramRatioUncertainty->GetSumw2N() ) histogramRatioUncertainty->Sumw2();
  histogramRatioUncertainty->SetMarkerColor(10);
  histogramRatioUncertainty->SetMarkerSize(0);
  setStyle_uncertainty(histogramRatioUncertainty);

  int numBins_bottom = histogramRatio->GetNbinsX();
  for ( int iBin = 1; iBin <= numBins_bottom; ++iBin ) {
    double binContent_data = histogram_data_density->GetBinContent(iBin);
    double binError_data = histogram_data_density->GetBinError(iBin);
    double binContent_mc = 0;
    double binError_mc = 0;
    if ( histogramSum_mc && histogramErr_mc ) {
      binContent_mc = histogramSum_mc_density->GetBinContent(iBin);
      binError_mc = histogramErr_mc_density->GetBinError(iBin);
    } else {
      TList* histograms = histogramStack_mc->GetHists();
      TIter nextHistogram(histograms);
      double binError2_mc = 0.;
      while ( TH1* histogram_density = dynamic_cast<TH1*>(nextHistogram()) ) {
        binContent_mc += histogram_density->GetBinContent(iBin);
        binError2_mc += square(histogram_density->GetBinError(iBin));
      }