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

void LATcorr::UpdateErrorbars()
{

   if(afterR)	 afterR->Sumw2();
   if(beforeR)	 beforeR->Sumw2(); 

}

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* 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* 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* fillHistogram(TTree* testTree, const std::string& varName, const std::string& selection, const std::string& frWeightName,
		   const std::string& histogramName, unsigned numBinsX, double xMin, double xMax)
{
  std::cout << "<fillHistogram>:" << std::endl;
  std::cout << " testTree = " << testTree << std::endl;
  std::cout << " varName = " << varName << std::endl;
  std::cout << " selection = " << selection << std::endl;
  std::cout << " frWeightName = " << frWeightName << std::endl;
  std::cout << " histogramName = " << histogramName << std::endl;

  TH1* histogram = new TH1F(histogramName.data(), histogramName.data(), numBinsX, xMin, xMax);
  histogram->Sumw2();

  TFile* dummyOutputFile = new TFile("dummyOutputFile.root", "RECREATE");

  TTree* selTree = ( selection != "" ) ? testTree->CopyTree(selection.data()) : testTree;
  std::cout << " selTree = " << selTree << std::endl;

  Float_t eventWeight = 1.;
  selTree->SetBranchAddress("weight", &eventWeight);

  Float_t var = 0.;
  selTree->SetBranchAddress(varName.data(), &var);
  
  Float_t frWeight = 1.;
  if ( frWeightName != "" ) {
    std::cout << "--> setting branch-address of fake-rate weight..." << std::endl;
    selTree->SetBranchAddress(frWeightName.data(), &frWeight);
  }

  int numEntries = selTree->GetEntries();
  std::cout << "--> numEntries = " << numEntries << std::endl;
  //if ( maxEntriesToProcess != -1 ) numEntries = TMath::Min(numEntries, maxEntriesToProcess);
  for ( int iEntry = 0 ; iEntry < numEntries; ++iEntry ) {
    selTree->GetEvent(iEntry);  

    //std::cout << "iEntry = " << iEntry << ": var = " << var << ", frWeight = " << frWeight << std::endl;

    if ( frWeightName != "" ) {
      if ( TMath::Abs(frWeight) < 1. ) // some entries have weight O(-100)
                                       // --> indication of technical problem with k-NearestNeighbour tree ?
	histogram->Fill(var, frWeight*eventWeight);
    } else {
      histogram->Fill(var, eventWeight);
    }
  }

  delete dummyOutputFile;

  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* GetHist(const string histname)
{
	//hists are already scaled to 10fb-1
	TH1* h = dynamic_cast<TH1*> (files[0]->Get(histname.c_str()));
	if (h == NULL)
	{
		cout << "hist " << histname << " not found in " <<  "!" << endl;
		assert (false);
	}
	TH1* hist = dynamic_cast<TH1*> (h->Clone());
	hist->Sumw2();
	hist->SetLineWidth(2);

	return hist;
}

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

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

vector<TH1*> GetHist(const string histname)
{
	//hists are already scaled to 10fb-1
	vector<TH1*> hists;
	for (int i=0; i< nBins; ++i)
	{
		TH1* h = dynamic_cast<TH1*> (files[i]->Get(histname.c_str()));
		if (h == NULL)
		{
			cout << "hist " << histname << " not found in file # " << i <<  " !" << endl;
			assert (false);
		}
		TH1* hist = dynamic_cast<TH1*> (h->Clone());
		hist->Sumw2();
		hist->SetLineWidth(2);
		hists.push_back(hist);
	}

	return hists;
}

  // Do the loop here, so that we can use options like "errors"
  void Draw( const TString & xTitle = "", const TString & yTitle = "", const bool errors = false ) {

    // Create a new THStack so that it handle tha maximum
    // THStack stack(name_, title_);
    THStack * stack = new THStack(name_, title_);

    int colorIndex = 0;
    if( !(histoList_.empty()) ) {
      std::vector<TH1*>::iterator histoIter = histoList_.begin();
      for( ; histoIter != histoList_.end(); ++histoIter, ++colorIndex ) {
        TH1 * histo = *histoIter;
        if(errors) histo->Sumw2();
        // histo->SetNormFactor(1);
        if( colorIndex < 4 ) histo->SetLineColor(colors_[colorIndex]);
        else histo->SetLineColor(colorIndex);
        // Draw and get the maximum value
        TString normalizedHistoName(histo->GetName());
        TH1 * normalizedHisto = (TH1*)histo->Clone(normalizedHistoName+"clone");
        normalizedHisto->Scale(1/normalizedHisto->Integral());
        stack->Add(normalizedHisto);
      }
      // Take the maximum of all the drawed histograms
      // First we need to draw the histogram, or getAxis() will return 0... (see root code...)
      canvas_->Draw();
      canvas_->cd();
      stack->Draw("nostack");
      stack->GetYaxis()->SetTitleOffset(1.2);
      stack->GetYaxis()->SetTitle(yTitle);
      stack->GetXaxis()->SetTitle(xTitle);
      stack->GetXaxis()->SetTitleColor(kBlack);
      stack->Draw("nostack");
      legend_->Draw("same");

      canvas_->Update();
      canvas_->Draw();
      canvas_->ForceUpdate();
      //canvas_->Print("test.pdf");
      canvas_->Write();

    }
  }

TH1* getMonitorElement(TFile* inputFile, const TString& dqmDirectory, const char* dqmSubDirectory, const TString& meName)
{
  TString meName_full = TString("DQMData").Append("/");
  if ( dqmDirectory != "") meName_full.Append(dqmDirectory).Append("/");
  meName_full.Append(dqmSubDirectory).Append("/").Append(meName);
  std::cout << "meName_full = " << meName_full << std::endl;
  
  TH1* me = (TH1*)inputFile->Get(meName_full);
  std::cout << "me = " << me <<  std::endl;
  
  //if ( !me->GetSumw2() ) me->Sumw2();
  me->Sumw2();

  me->Rebin(2);

  me->Scale(1./me->Integral());

  me->SetMaximum(1.);
  me->SetStats(false);

  return me;
}

  TH1 *Draw(const plotVar_t& pv, const TCut& cut, const TCut& cutSQ ) {
    cout << "\tDrawing " << pv.plotvar << " for sample = " << info_.samplename << " ... ";
    TString hname = TString("th1")+ pv.outfile + Form("%d",info_.index);
    TH1 *histo = new TH1D(hname, hname, pv.NBINS, pv.MINRange, pv.MAXRange);
    assert(histo);
    histo->Sumw2();
    assert(tree_);
    cout << tree_->Draw(pv.plotvar+TString(">>")+hname, cut, "goff") << " events" << endl;

    if (strlen((const char *)cutSQ)) {
      hname = TString("th1") + pv.outfile + Form("%d",info_.index) + TString("SQ");
      TH1 *histoSQ = new TH1D(hname, hname, pv.NBINS, pv.MINRange, pv.MAXRange);
      tree_->Draw(pv.plotvar+TString(">>")+hname, cutSQ, "goff");
      for(int hi=1;hi<=pv.NBINS;hi++) {
	histo->SetBinError(hi,sqrt(histoSQ->GetBinContent(hi)));
      }
      delete histoSQ;
    }
    if (info_.nMCevents)
      histo->Scale(info_.xsecpblumi*info_.otherscale/(double)info_.nMCevents);

    return histo;
  }