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

void QAvertex(const Char_t *fdata, const Char_t *fmc)
{
  style();

  TFile *fdtin = TFile::Open(fdata);
  TList *ldtin = (TList *)fdtin->Get("clist");
  TH2 *hdtin = (TH2 *)ldtin->FindObject("zv");
  TH1 *hdt = (TH1 *)ldtin->FindObject("zvNoSel");
  SetHistoStyle(hdt, 20, kRed+1);
  hdt->Scale(1. / hdt->Integral());

  TH1 *hdt0010 = hdtin->ProjectionX("hdt0010", 1, 4);
  SetHistoStyle(hdt0010, 20, kRed+1);
  hdt0010->Scale(1. / hdt0010->Integral());

  TH1 *hdt7080 = hdtin->ProjectionX("hdt7080", 11, 11);
  SetHistoStyle(hdt7080, 25, kAzure-3);
  hdt7080->Scale(1. / hdt7080->Integral());
  
  TFile *fmcin = TFile::Open(fmc);
  TList *lmcin = (TList *)fmcin->Get("clist");
  TH1 *hmc = (TH1 *)lmcin->FindObject("zvNoSel");
  SetHistoStyle(hmc, 25, kAzure-3);
  hmc->Scale(1. / hmc->Integral());

  TCanvas *c = new TCanvas("cVertex", "cVertex", 800, 800);
  TH1 * hfr = c->DrawFrame(-20., 0., 20., 0.1);
  hfr->SetTitle(";#it{z}_{vtx};");
  hdt0010->Draw("same");
  hdt7080->Draw("same");
  TLegend *legend = new TLegend(0.20, 0.18+0.60, 0.50, 0.30+0.60);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(hdt0010, "0-10%", "p");
  legend->AddEntry(hdt7080, "70-80%", "p");
  legend->Draw("same");
  c->SaveAs(canvasPrefix+"vertex.pdf");
  
  TCanvas *c1 = new TCanvas("cVertexDataMC", "cVertexDataMC", 800, 800);
  hfr = c1->DrawFrame(-20., 0., 20., 0.1);
  hfr->SetTitle(";#it{z}_{vtx};");
  hdt->Draw("same");
  hmc->Draw("same");
  legend = new TLegend(0.20, 0.18+0.60, 0.50, 0.30+0.60);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(hdt, "data", "p");
  legend->AddEntry(hmc, "Monte Carlo", "p");
  legend->Draw("same");
  c1->SaveAs(canvasPrefix+"vertexDataMC.pdf");
  
  //return 0;  
}

// Extract and merge histograms in range [from, to) (to is not included)
//
TH1 *merge(const string &path, TFile **input, const int &from, const int &to,
        const bool &do_normalize = false)
{
    TH1 *result = 0;
    for(int i = from; to > i; ++i)
    {
        TH1 *hist = get(path, input[i], i);
        if (!hist)
        {
            cerr << "failed to extract: " << path << endl;

            continue;
        }

        if (!result)
            result = dynamic_cast<TH1 *>(hist->Clone());
        else
            result->Add(hist);
    }

    if (do_normalize
            && result
            && result->GetEntries())
    {
        result->Scale(1. / result->Integral());
    }

    return result;
}

///-----------------------------------------------------------------------------
Double_t effSigma(RooAbsPdf *pdf, RooRealVar *obs, Int_t nbins)
{
  TH1 *hist = pdf->createHistogram(obs->GetName(), nbins);
  hist->Scale(nbins);

  return effSigma( hist);
}

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

void QAcentrality(const Char_t *fdata)
{
  style();

  TFile *fin = TFile::Open(fdata);
  TList *lin = (TList *)fin->Get("clist");
  lin->ls();
  TH1 *hin = (TH1 *)lin->FindObject("EvCentrDist");
  Float_t sum = 1.2 * hin->Integral(hin->FindBin(0.1), hin->FindBin(79.9));
  hin->Scale(1. / sum);
  SetHistoStyle(hin, 20, kRed+1);
  TCanvas *c = new TCanvas("cQAcentrality", "cQAcentrality", 800, 800);
  TH1 * hfr = c->DrawFrame(0., 0.005, 100., 0.015);
  hfr->SetTitle(";centrality percentile;events");
  hin->Draw("same");
  c->SaveAs(canvasPrefix+"centrality.pdf");

  TH2 *hinv0 = (TH2 *)lin->FindObject("V0");
  TCanvas *cv0 = new TCanvas("cQAcentralityV0", "cQAcentralityV0", 800, 800);
  cv0->SetLogx();
  cv0->SetLogz();
  //  TH1 * hfrv0 = cv0->DrawFrame(100., -0.5, 50000., 10.5);
  // DrawBinLabelsY(hfrv0, kTRUE);
  // hfrv0->SetTitle(";V0 signal;");
  //hinv0->Draw("same,col");
  hinv0->Draw("col");
  cv0->SaveAs(canvasPrefix+"centralityV0.pdf");
}

void compareDYTemplates(TString baseURL="~/scratch0/top-newjec/syst_plotter.root")
{
  TString ch[]={"ee","mumu"};
  TString categs[]={"eq1jets",""};
  for(size_t ich=0; ich<2; ich++)
    {
      for(size_t icat=0; icat<2; icat++)
	{
	  TObjArray lowMet  = getDistributionFromPlotter(ch[ich]+"_"+categs[icat]+"lowmetdilarccosine",baseURL);
	  TH1 *lowMetH = (TH1 *) ((TList *)lowMet.At(3))->At(3);
	  formatPlot(lowMetH,1,1,1,24,0,false,false,1,1,1);
	  lowMetH->SetTitle("E_{T}^{miss}<30 GeV/c^{2}");
	  lowMetH->Scale(1./lowMetH->Integral());

	  TObjArray highMet = getDistributionFromPlotter(ch[ich]+"_"+categs[icat]+"dilarccosine",baseURL);
	  TH1 *highMetH = (TH1 *)  ((TList *)highMet.At(3))->At(3);
	  formatPlot(highMetH,1,1,1,20,0,false,false,1,1,1);
	  highMetH->SetTitle("E_{T}^{miss}>30 GeV/c^{2}");
	  highMetH->Scale(1./highMetH->Integral());

	  TString channelTitle(ich==0 ? "ee" : "#mu#mu"); 
	  if(categs[icat]=="eq1jets") channelTitle += "+ 1 jet";
	  else                        channelTitle += "+ #geq 2 jets";

	  //draw
	  TString plot(ch[ich]+categs[icat]+"_anglecomparison");
	  TCanvas *cnv = getNewCanvas(plot+"c",plot+"c",false);
	  cnv->Clear();
	  cnv->SetWindowSize(600,600);
	  cnv->cd(); 
	  TList *mc = new TList; mc->Add(lowMetH);
	  TList *data = new TList; data->Add(highMetH);
	  TList *spimpose = new TList;
	  TLegend *leg=showPlotsAndMCtoDataComparison(cnv,*mc,*spimpose,*data,false);
	  TPad *p=(TPad *)cnv->cd(1);
	  formatForCmsPublic(p,leg,"CMS simulation, " + channelTitle, 4);
	  cnv->SaveAs(plot+".C");
	  cnv->SaveAs(plot+".pdf");
	  cnv->SaveAs(plot+".png");
	}
    }
}

void rescaleBoundaryHists(std::string infile, int numSamples=-1){

  TFile* f = new TFile(infile.c_str(), "UPDATE");
  TDirectory* dir = 0;

  TIter dir_it(f->GetListOfKeys());
  TKey* dir_k;
  while ((dir_k = (TKey *)dir_it())) {
    if (TString(dir_k->GetClassName()) != "TDirectoryFile") continue;
    std::string dir_name = std::string(dir_k->GetTitle());
    if(dir_name == "") continue;
    dir = (TDirectory*)dir_k->ReadObj();
    if(dir == 0) continue;
    TIter hist_it(dir->GetListOfKeys(), kIterBackward);
    TKey* hist_k;
    while ((hist_k = (TKey *)hist_it())) {
      std::string hist_name = (hist_k->GetTitle());
      if (hist_name.find("_HI") != std::string::npos || hist_name.find("_LOW") != std::string::npos || hist_name.find("h_n_mt2bins") != std::string::npos) {
        TH1* h = (TH1*)hist_k->ReadObj();
        if(numSamples==-1)
            h->Scale(1.0/h->GetEntries());
        else
            h->Scale(1.0/numSamples);
        dir->cd();
        h->Write("",TObject::kOverwrite);
      }
    }
  }

  delete dir;

  gDirectory->GetList()->Delete();
  
  f->Write("",TObject::kOverwrite);
  f->Close();
  delete f;

}

TH1* makeCDF(TH1* h) {
	TString sName(TString(h->GetName())+TString("_CDF"));
	TString sTitle(TString(h->GetTitle())+TString(" CDF"));
	TH1* hOut = (TH1*) h->Clone(sName);
	hOut->SetTitle(sTitle);
	hOut->Reset();
	double cdf = 0;
	for (int ibin=0; ibin < h->GetNbinsX()+2; ++ibin) {
		cdf += h->GetBinContent(ibin);
		hOut->SetBinContent(ibin,cdf);
	}
	hOut->Scale(1.0/(h->Integral(0,h->GetNbinsX()+1)));
	return hOut;
}

void plotter::draw_delta_rel(TH1* hist_, TH1* result_, TString file_name){
  TH1* hist = (TH1*) hist_->Clone("hist");
  TH1* result = (TH1*) result_->Clone("result");
  hist->Divide(result);
  hist->Scale(100);
  TCanvas *c= new TCanvas("Particle Level","",600,600);
  gPad->SetLeftMargin(0.15);
  hist->SetTitle(file_name);
  hist->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  hist->GetYaxis()->SetTitle("relative uncertainty [%]");
  hist->GetYaxis()->SetTitleOffset(1.5);
  hist->GetYaxis()->SetNdivisions(505);
  hist->SetFillColor(810);
  hist->SetLineColor(810);
  hist->Draw("HIST");
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

//________________________________________________________
Int_t GFOverlay::AddHistsAt(const TObjArray &hists, const TObjArray &legends, Int_t layer,Int_t pos)
{
  // hists and legends must have same length, but might have gaps...
  // return number of hists found and added

  Int_t nHists = 0;
  for (Int_t iHist = 0; iHist < hists.GetEntriesFast(); ++iHist) {
    TH1 *hist = static_cast<TH1*>(hists[iHist]);
    if (!hist) continue;

    if (fNormalise && hist->GetEntries()) {
      hist->Scale(1./hist->GetEntries());
    }

    fHistMan->AddHistSame(hist, layer, pos, (legends[iHist] ? legends[iHist]->GetName() : 0));
    ++nHists;
  }

  return nHists;
}

std::shared_ptr<TH1> getHistogram(const std::string& name, const std::vector<Input>& inputs, int type) {

  TH1* h = nullptr;

  for (const auto& input: inputs) {
    if (input.type != type)
      continue;

    TH1* f = static_cast<TH1*>(input.file->Get(name.c_str()));
    f->Scale(input.cross_section / (input.generated_events * input.top_pt_weight));

    if (! h) {
      h = static_cast<TH1*>(f->Clone());
      h->SetDirectory(NULL);
    } else
      h->Add(f);
  }

  return std::shared_ptr<TH1>(h);
}

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

// Called just after the main event loop
// Can be used to write things out, dump a summary etc
// Return non-zero to indicate a problem
int IslandAmplitude::AfterLastEntry(TGlobalData* gData,const TSetupData *setup){

  // Print extra info if we're debugging this module:
  if(Debug()){
     cout<<"-----IslandAmplitude::AfterLastEntry(): I'm debugging!"<<endl;
  }

  double run_norm = fAmpNorm->Integral(0,-1);
  for(mapSH_t::iterator it = fAmpHist.begin(); it != fAmpHist.end(); ++it)
    {
      TH1F* h = it->second;
      TObject* obj = h->Clone((std::string(h->GetName()) + "_RunNorm").c_str());
      TH1* hn = static_cast<TH1*>(obj);
      hn->SetTitle((std::string(h->GetTitle()) + " (run normalized)").c_str());
      hn->Scale(1.0/run_norm);
    }

  for(mapSH_t::iterator it = fAmpHistNorm.begin(); it != fAmpHistNorm.end(); ++it)
    it->second->Scale(1.0/fNProcessed);

  return 0;
}

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