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

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

// PValue for finding a local p-value as observed in 'bin' or worse 
void ToyExperiments::printGlobalPValueOfLocalFluctuation(unsigned int bin, unsigned int nExperiments) const {
  std::cout << "Determining global p-value for observed fluctuation in bin " << bin << " from " << nExperiments << " toy experiments ...  " << std::flush;

  // Find the predicted yields that correspond
  // to the local p-value 'localPValue'
  std::vector<unsigned int> limitYields = yields(localPValue(bin,observedYields_.at(bin)));

  TH1* hIsAbovePValue = new TH1D("hIsAbovePValue","",2,0,2);

  const double minCorr = findMinValidRandomNumberForCorrelatedUncertainties();

  for(unsigned int p = 0; p < nExperiments; ++p) {
    bool isAbovePValue = false;
    double rCorr = rand_->Gaus(0.,1.);
    while( rCorr <= minCorr ) {
      rCorr = rand_->Gaus(0.,1.);
    }
    for(unsigned int b = 0; b < Parameters::nBins(); ++b) {
      double prediction = -1.;
      bool negativePrediction = true;
      while( negativePrediction ) {
	double rUncorr = rand_->Gaus(0.,1.);
	double uncorrVar = rUncorr * uncorrelatedUncerts_.at(b);
	double corrVar   = rCorr   * correlatedUncerts_.at(b);
	prediction = meanPredictions_.at(b) + uncorrVar + corrVar;
	if( prediction >= 0. ) {
	  negativePrediction = false;
	}
      }
      double predictedYield = rand_->Poisson(prediction);
      if( predictedYield >= limitYields.at(b) ) {
	isAbovePValue = true;
	break;
      }      
    }
    if( isAbovePValue ) {
      hIsAbovePValue->Fill(1);
    } else {
      hIsAbovePValue->Fill(0);
    }
  }
  std::cout << "ok" << std::endl;

  double lpv      = localPValue(bin,observedYields_.at(bin));
  double gpUncorr = 1. - pow(1.-lpv,Parameters::nBins());
  double gpCorr   = hIsAbovePValue->Integral(2,2)/hIsAbovePValue->Integral(1,2);

  std::cout << "\n\n----- Global p-value for observed fluctuation in bin " << bin << " -----" << std::endl;
  std::cout << "  local p-value                           : " << lpv << " (" << TMath::NormQuantile(1.-lpv) << "sig)" << std::endl;
  std::cout << "  global p-value (without correlations)   : " << gpUncorr  << " (" << TMath::NormQuantile(1.-gpUncorr) << "sig)" << std::endl;
  std::cout << "  global p-value (including correlations) : " << gpCorr  << " (" << TMath::NormQuantile(1.-gpCorr) << "sig)" << std::endl;
}

void cutFlow(bool eff = false)
{
    using namespace std;
    
    double integral, preInt = 1.0;
    
    vector<pair<string,TFile*> > fnames;
    vector<pair<string,string> > hnames;

    hnames.push_back(pair<string,string>("none","hNu/cut0_none/mWR"));
    hnames.push_back(pair<string,string>("LLJJ Pt","hNu/cut1_LLJJpt/mWR"));
    hnames.push_back(pair<string,string>("trig","hNu/cut2_TrigMatches/mWR"));
    hnames.push_back(pair<string,string>("vertex","hNu/cut3_Vertex/mWR"));
    hnames.push_back(pair<string,string>("mu1 pt","hNu/cut4_Mu1HighPt/mWR"));
    hnames.push_back(pair<string,string>("Mll","hNu/cut5_diLmass/mWR"));
    hnames.push_back(pair<string,string>("MWR","hNu/cut6_mWRmass/mWR"));

    fnames.push_back(pair<string,TFile*>("ttbar powheg",   new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_TTTo2L2Nu2B_7TeV-powheg-pythia6.root")));
    fnames.push_back(pair<string,TFile*>("ttbar madgraph", new TFile("/local/cms/user/dahmes/wr2011/bgMC/Summer11/aug30/ttbar-PFJets.root")));
    fnames.push_back(pair<string,TFile*>("Z+Jets sherpa",  new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_DYToLL_M-50_7TeV-sherpa.root")));
    fnames.push_back(pair<string,TFile*>("W+Jets sherpa",  new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_WToLNu_7TeV-sherpa.root")));
    fnames.push_back(pair<string,TFile*>("ZZ",             new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_ZZ_TuneZ2_7TeV_pythia6_tauola.root")));
    fnames.push_back(pair<string,TFile*>("WZ",             new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_WZ_TuneZ2_7TeV_pythia6_tauola.root")));
    fnames.push_back(pair<string,TFile*>("WW",             new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_WW_TuneZ2_7TeV_pythia6_tauola.root")));
    fnames.push_back(pair<string,TFile*>("tW",             new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_T_TuneZ2_tW-channel-DR_7TeV-powheg-tauola.root")));
    fnames.push_back(pair<string,TFile*>("tbarW",          new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_Tbar_TuneZ2_tW-channel-DR_7TeV-powheg-tauola.root")));
    
	printf("%8s", "cut");
	for(vector<pair<string,TFile*> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
	{
		printf(eff?" & %15s        ":" & %15s", i->first.c_str());
	}
	printf(" \\\\ \\hline\n");
	for(vector<pair<string,string> >::const_iterator j = hnames.begin(); j != hnames.end(); j++)
	{
		printf("%8s", j->first.c_str());
		for(vector<pair<string,TFile*> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
		{
			TH1* h = (TH1*)i->second->Get(j->second.c_str());
			integral = h->Integral(0, h->GetNbinsX()+1);
			
			printf((j == hnames.begin() && eff)?" & %15.0f        ":" & %15.0f", integral);
			if(j != hnames.begin() && eff) 
			{
				TH1* hpre = (TH1*)i->second->Get((j-1)->second.c_str());
				preInt = hpre->Integral(0, h->GetNbinsX()+1);
				printf("  (%4.2f)", integral/preInt);
			}
		} 
		printf(" \\\\ \\hline\n");
	}
}

void fadc_fit_heights(std::string fadc_id,std::string hist_type)
{
  /*****************************************************************/
  // Prepare the canvas
  gStyle->SetOptFit(1111);
  TCanvas *AlCapCanvas = (TCanvas *) gROOT->GetListOfCanvases()->At(0);
  AlCapCanvas->Clear();
  AlCapCanvas->Divide(4,2);

  //  gROOT->ProcessLine(".L common/get_histogram.C"); // get_histogram() is called here
  /*****************************************************************/
  const int n_channels = 8;
  std::string bank_names[n_channels] = {"Na", "Nb", "Nc", "Nd", "Ne", "Nf", "Ng", "Nh"};
  std::string name;

  double mean = 0;
  double rms = 0;

  for (int iChn = 0; iChn < n_channels; iChn++) {
    name=bank_names[iChn]+fadc_id;
    TH1* hist = get_histogram(name, hist_type);
    mean = hist->GetMean();
    rms = hist->GetRMS();

      AlCapCanvas->cd(iChn+1);

    if (hist->Integral()!=0) {
      TF1* gaus = new TF1("gaus","gaus",mean-rms,mean+rms);
      hist->Fit("gaus","Q");
      hist->GetXaxis()->SetRange(mean-2*rms,mean+2*rms);
      hist->Draw();
    }
  }
}

void makeplotTwo(TH1& hist1, TH1& hist2, const char* plotname, int logy) {

  hist1.SetLineColor(4);
  hist1.SetMarkerColor(4);
  gStyle->SetOptStat(0);
  TCanvas* can = new TCanvas( plotname, plotname, 500, 500);
  hist1.Draw( );
  hist2.Draw( "same" );
  TLegend *leg = new TLegend(0.55,0.8,0.89,0.92);
  leg->AddEntry( &hist1,"Reconstructed","LP");
  leg->AddEntry( &hist2,"Generated","L");
  leg->SetFillColor(0);
  leg->SetLineColor(0);
  leg->Draw();
  can->SetLogy( logy );
  if(logy==2) hist1.GetYaxis()->SetMoreLogLabels();

  std::string plot(cmEnergy);
  plot.append(plotname);
  can->SaveAs( (plot+".eps").c_str() );
  can->SaveAs( (plot+".gif").c_str() );
  can->SaveAs( (plot+".root").c_str() );
  // delete can;
  cout << hist1.Integral() << endl;
}

// show the histogram in first slot, try a Gaussian fit with given parameters
void PRadHistCanvas::UpdateHist(int index, TObject *tob, int range_min, int range_max)
{
    --index;
    if(index < 0 || index >= canvases.size())
        return;

    canvases[index]->cd();
    canvases[index]->SetGrid();
    gPad->SetLogy();

    TH1 *hist = (TH1*)tob;

    hist->GetXaxis()->SetRangeUser(hist->FindFirstBinAbove(0,1) - 10,
                                   hist->FindLastBinAbove(0,1) + 10);

    hist->GetXaxis()->SetLabelSize(HIST_LABEL_SIZE);
    hist->GetYaxis()->SetLabelSize(HIST_LABEL_SIZE);

    // try to fit gaussian in certain range
    if(range_max > range_min
       && hist->Integral(range_min, range_max + 1) > 0)
    {
        TF1 *fit = new TF1("", "gaus", range_min, range_max);
        fit->SetLineColor(kRed);
        fit->SetLineWidth(2);
        hist->Fit(fit,"qlR");
    }

    hist->SetFillColor(fillColors[index]);
    hist->Draw();
    canvases[index]->Refresh();
}

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

void makeplotThree(TH1& hist1, TH1& hist2, TH1& hist3, const char* plotname, int log) {

  hist1.SetLineColor(4);
  hist1.SetMarkerColor(4);
  hist3.SetLineColor(2);
  hist3.SetMarkerColor(2);

  gStyle->SetOptStat(0);
  TCanvas* can = new TCanvas( plotname, plotname, 500, 500);
  hist1.Draw( );
  hist2.Draw( "same" );
  hist3.Draw( "same" );
  TLegend *leg = new TLegend(0.48,0.7,0.89,0.92);
  leg->AddEntry( &hist1,"CaloJet","LP");
  leg->AddEntry( &hist3,"PF Jet","LP");
  leg->AddEntry( &hist2,"GenJet","L");
  leg->SetFillColor(0);
  leg->SetLineColor(0);
  leg->Draw();
  can->SetLogy( log );
  if(log==2) {
    can->SetLogy( 1 );
    hist1.GetYaxis()->SetMoreLogLabels();
  }

  std::string plot("ratio-");
  plot.append(plotname);
  can->SaveAs( (plot+".eps").c_str() );
  can->SaveAs( (plot+".gif").c_str() );
  can->SaveAs( (plot+".root").c_str() );
  // delete can;
  cout << hist1.Integral() << endl;
}

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 GetSigEfficiency(TString path, TString tag){
  
  
  setTDRStyle();
  //gStyle->SetPalette(1);
  
  TFile * file = new TFile(path);
  
  TH1* hnsig =   (TH1F*)file->Get(("eventcutflow"));
  float nsig = float(hnsig->GetBinContent(2)); 
  
  TString hist = tag + "MassRegion/h_Nelectrons_" + tag + "MassRegion";
  TH1* h =  (TH1*)file->Get(hist.Data());
  cout << h << endl;
  cout << h->Integral() << endl;
  cout << nsig << endl;
  cout << "Total efficiency  " << tag << " = " << 100* (h->Integral() / nsig) << endl;
}

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

TH1* getMonitorElement(TFile* inputFile, const TString& dqmDirectoryName, const char* processName, const TString& meName)
{
  TString meName_full = TString("DQMData/tauFakeRate/harvested").Append("/").Append(processName).Append("/").Append(dqmDirectoryName).Append("/").Append(meName);
  //TString meName_full = TString("DQMData").Append("/").Append(dqmDirectoryName).Append("/").Append(meName);
  std::cout << "meName_full = " << meName_full << std::endl;
  
  TH1* me = (TH1*)inputFile->Get(meName_full);
  std::cout << "me = " << me <<  std::endl;
  std::cout << "integral = " << me->Integral() << std::endl;
  
  return me;
}

double IntegralStack(THStack* stack) {

  double val=0;
  TList* list = stack->GetHists();
  TIter it(list, true);
  TObject* obj=0;
  while( (obj = it.Next()) ) {
    TH1* h = dynamic_cast<TH1*>(obj);
    val += h->Integral();
  }
  return val;
}

TEST_F(EvalHistMethods, CreateHistogram1D) {
  evaluator->SetNormalizationBuffer(norm);
  evaluator->SetParameterBuffer(params);

  TH1* hist = evaluator->CreateHistogram();

  EXPECT_EQ(2, hist->GetNbinsX());
  ASSERT_FLOAT_EQ(1.0, hist->Integral("width"));
  ASSERT_FLOAT_EQ(1.6, hist->GetBinContent(hist->FindBin(0.25)));
  ASSERT_FLOAT_EQ(0.4, hist->GetBinContent(hist->FindBin(0.75)));

  delete hist;
}