C++ TGraphAsymmErrors类代码示例

TGraphAsymmErrors *getEfficiencyFrom(TH1F *h)
{
  if(h==0) return 0;

  // compute efficiency: N(x>cut) / N as function of x
  // we assume overflows / underflows are included
  
  TGraphAsymmErrors *effgr = new TGraphAsymmErrors;
  effgr->SetName(h->GetName()+TString("eff"));
  effgr->SetTitle(TString("#varepsilon_{")+h->GetTitle()+TString("}"));
  effgr->SetMarkerStyle(h->GetMarkerStyle());
  effgr->SetMarkerSize(1.3);
  effgr->SetMarkerColor(h->GetLineColor());
  effgr->SetLineColor(h->GetLineColor());
  effgr->SetFillStyle(0);
  int nbins=h->GetXaxis()->GetNbins();
  Double_t tot=h->Integral();
  for(int ibin=1; ibin<=nbins; ibin++)
    {
      Double_t cut = h->GetBinLowEdge(ibin);
      Double_t err;
      Double_t tot_cut=h->IntegralAndError(ibin,nbins,err);
      effgr->SetPoint(ibin-1,cut,tot_cut/tot);
      effgr->SetPointError(ibin-1,0,0,err/tot,err/tot);
    }
  
  return effgr;
}

// nameSuffix is appended to histograms name, making the graph's name
TGraphAsymmErrors* makeEffGraph (TH1F& h, TString nameSuffix = "eff")
{
  int n = 1+h.GetNbinsX();

  vector<Double_t> ps(n); // a partial sum: integral of h from -infinity to high edge of bin n, counted TH1 style
  for (int i=0; i<n; ++i) {
    ps [i] = h.GetBinContent (i) + (i>0 ? ps [i-1] : 0);
  }

  Double_t *x = new Double_t [n], *y = new Double_t [n], *exl = new Double_t [n], *exh = new Double_t [n];
  Double_t *eyl = new Double_t [n], *eyh = new Double_t [n]; // for the graph

  for (int i=0; i<n; ++i) {
    exl [i] = exh [i] = 0;
    x [i] = h.GetBinLowEdge (i+1); // low edge of nBins+1 is the right edge of h.
    Double_t eff, deff;
    calcPoisEff (ps[n-1] - ps[i], ps[i], eff, deff);
    y [i] = eff;
    makeAsymmEffErrors (eff, deff, eyl [i], eyh [i]);
  }
  TGraphAsymmErrors *gae = new TGraphAsymmErrors (n, x, y, exl, exh, eyl, eyh);
  gae->SetTitle (h.GetTitle());
  gae->GetXaxis()->SetTitle (h.GetXaxis()->GetTitle());
  gae->GetYaxis()->SetTitle ("Efficiency");
  gae->SetName (TString (h.GetName()) + nameSuffix);

  delete[] x;
  delete[] y;
  delete[] exl;
  delete[] exh;
  delete[] eyl;
  delete[] eyh;
  return gae;
}

void applyEfficiencyToTauMuons() {

    TFile* _fPtSm = new TFile("smearedMuonPtFromTau.07.09.2013.root","read");
    TFile* _fEff = new TFile("mcWEffForTauMuonStudy_07_08_2013.root","read");
    TFile* _fPtSmEff = new TFile("smearedEffAppliedMuonPtFromTau.root","recreate");

    _fPtSm->cd();
    TH1F* hPtSm = (TH1F*)_fPtSm->Get("hSmearedMuonPt");
    _fEff->cd();
    TGraphAsymmErrors* grEff = (TGraphAsymmErrors*)_fEff->Get("pEffWAccCuts");

    double* yEff = grEff->GetY();
    TH1F*  hnew = (TH1F*)hPtSm->Clone("hnew");
    ///Apply efficiency bin by bin
    for(int igr=0; igr<grEff->GetN(); ++igr) {

        double ptNominal = hPtSm->GetBinContent(igr+1);
        std::cout << "pt nominal " << ptNominal << " eff " << yEff[igr] << " = " << yEff[igr]*ptNominal << std::endl;
        hnew->SetBinContent(igr+1,yEff[igr]*ptNominal);
    }

    _fPtSmEff->cd();
    hnew->Write();
    _fPtSmEff->Write();
    _fPtSmEff->Close();
}

void
printVectorsToFile(wGraph_t *wg) // , const string& filename)
{
  if (wg->gr->InheritsFrom("TGraphAsymmErrors")) {
    TGraphAsymmErrors *agr = (TGraphAsymmErrors *)wg->gr;
    cout<<"lox\tctrx\thix\tloy\tctry\thiy\tex\tey"<<endl;
    for (int i=0; i<agr->GetN(); i++) {
      double x,y, ex, ey, lox,hix,loy,hiy;
      agr->GetPoint(i,x,y);
      ex = agr->GetErrorX(i);
      ey = agr->GetErrorY(i);
      lox = x-ex; hix = x+ex;
      loy = y-ey; hiy = y+ey;
      cout <<lox<<"\t"<<x<<"\t"<<hix<<"\t";
      cout <<loy<<"\t"<<y<<"\t"<<hiy<<"\t"<<ex<<"\t"<<ey<<endl;
    }
  } else {
    cout<<"x\ty\t"<<endl;
    for (int i=0; i<wg->gr->GetN(); i++) {
      double x,y;
      wg->gr->GetPoint(i,x,y);
      cout<<x<<"\t"<<y<<endl;
    }
  }
}                                                  // printVectorsToFile

//*************************************************************************************************
//
//*************************************************************************************************
TGraphAsymmErrors* MakeCurrentWPSigEffVsBkgEffGraph(Double_t signalEff, Double_t bkgEff, string name ) {
  //Make Met Plots
  double SigEff[1];
  double BkgEff[1];
  double SigEffErrLow[1];
  double SigEffErrHigh[1];
  double BkgEffErrLow[1];
  double BkgEffErrHigh[1];
  double NSigTotal = 0;
  double NBkgTotal = 0;
  double cutValue;

  SigEff[0] = signalEff;
  SigEffErrLow[0] = 0;
  SigEffErrHigh[0] = 0;
  BkgEff[0] = bkgEff;
  BkgEffErrLow[0] = 0;
  BkgEffErrHigh[0] = 0;

  TGraphAsymmErrors *tmpSigEffVsBkgEff = new TGraphAsymmErrors (1, BkgEff, SigEff, BkgEffErrLow, BkgEffErrHigh , SigEffErrLow, SigEffErrHigh );
  tmpSigEffVsBkgEff->SetName(name.c_str());
  tmpSigEffVsBkgEff->SetTitle("");
  tmpSigEffVsBkgEff->GetXaxis()->SetTitle("Bkg Eff");
  tmpSigEffVsBkgEff->GetYaxis()->SetTitle("Signal Eff");
  tmpSigEffVsBkgEff->GetYaxis()->SetTitleOffset(1.1);
  tmpSigEffVsBkgEff->GetXaxis()->SetTitleOffset(1.05);
  tmpSigEffVsBkgEff->SetMarkerColor(kBlack);
  tmpSigEffVsBkgEff->SetLineColor(kBlack);
  tmpSigEffVsBkgEff->SetMarkerSize(1.5);

  return tmpSigEffVsBkgEff;
}

TGraphAsymmErrors* FilterBins(std::vector<int> binsToSelect, TGraphAsymmErrors* inputGR) {
 int numbin = binsToSelect.size();
 
 TString name = Form ("%s_new",inputGR->GetName());
 TString title = Form ("%s",inputGR->GetTitle());
 
 TGraphAsymmErrors* newGR = new TGraphAsymmErrors();
 newGR -> SetName (name);
 
 for (int i=0; i< binsToSelect.size(); i++) {
  
  double X = i+0.5;
  double Y = (inputGR->GetY()) [binsToSelect.at(i)-1];
  
  double errXUp      = inputGR->GetErrorXhigh(binsToSelect.at(i)-1);
  double errXDown    = inputGR->GetErrorXlow(binsToSelect.at(i)-1);
  double errYUp      = inputGR->GetErrorYhigh(binsToSelect.at(i)-1);
  double errYDown    = inputGR->GetErrorYlow(binsToSelect.at(i)-1);
  
  newGR->SetPoint(i, X, Y);
  newGR->SetPointError(i, errXDown, errXUp, errYDown, errYUp);
  
  //     std::cout << " i = " << i << " X = " << X << " Y = " << Y << std::endl;
 }
 
 return newGR;
}

TGraph* makeGraph_response_asymmetric(std::vector<histogram_vs_X_Type>& histograms_vs_X, const std::vector<double>& y_true_array)
{
  //std::cout << "<makeGraph_response_asymmetric>:" << std::endl;

  unsigned numPoints = histograms_vs_X.size();
  assert(numPoints > 0);

  TGraphAsymmErrors* graph = new TGraphAsymmErrors(numPoints);

  for ( unsigned iPoint = 0; iPoint < numPoints; ++iPoint ) {
    double x = histograms_vs_X[iPoint].x_;
    double xErrUp = histograms_vs_X[iPoint].xErrUp_;
    double xErrDown = histograms_vs_X[iPoint].xErrDown_;

    TH1* histogram = histograms_vs_X[iPoint].histogram_;

    double max, maxErrUp, maxErrDown;
    compHistogram_yMax(histogram, max, maxErrUp, maxErrDown);
    //std::cout << "histogram = " << histogram->GetName() << ": max = " << max << " + " << maxErrUp << " - " << maxErrDown << std::endl;

    double y_true = y_true_array[iPoint];
    //std::cout << " y(true) = " << y_true << std::endl;

    double y = max/y_true;
    double yErrUp = maxErrUp/y_true;
    double yErrDown = maxErrDown/y_true;

    graph->SetPoint(iPoint, x, y);
    graph->SetPointError(iPoint, xErrDown, xErrUp, yErrDown, yErrUp);
  }

  return graph;
}

TGraphAsymmErrors* ROC( TH1* hSignal, TH1* hQCD){
  TGraphAsymmErrors* out = new TGraphAsymmErrors();
  int nbins = hSignal->GetNbinsX();
  double totalS = hSignal->GetEntries();
  double totalB = hQCD->GetEntries();
  cout << "=============================================================" << endl;
  cout << hSignal->GetTitle() << endl;
  cout << "totalS = " << totalS << "      totalB= " << totalB << endl;
  cout << "=============================================================" << endl;
  int bin;
  for(int i=0 ; i < 51; i++){
    bin=i+1;
//    bin=i;
    double Bi = hQCD->Integral(0,bin);
    double Si = hSignal->Integral(0,bin);
    double eff_qcd = Bi / totalB;
    double eff_signal = Si/ totalS; 
    double err_qcd = sqrt(Bi) / totalB;
    double soverb = eff_signal/sqrt(eff_signal+eff_qcd);
    cout.setf(ios::fixed);
    cout.precision(5);
    cout << "isolation value = "  << 0.01*(bin) << "      signal eff = " << eff_signal << "      background eff = " << eff_qcd << "      s/sqrt(s+b)= " << soverb << endl ;
    cout.precision();
    //out->SetPoint(i, eff_signal, 1-eff_qcd);
    out->SetPoint(i, eff_qcd, eff_signal);
    out->SetPointEXhigh(i,err_qcd);
    out->SetPointEXlow(i, err_qcd);
  }   
  return out;
}

void MakeNsignalEff_pt15(){
  
  setTDRStyle();
  gStyle->SetPalette(1);

  TH1* medium = makehist("PreSelection_medium_pt15");
  TH1* tight = makehist("PreSelection_tight_pt15");
  TH1* tight_dxy10= makehist("PreSelection_iso_10_10_pt15");
  TH1* tight_anal = makehist("PreSelection_pt15");


  TLegend* legendH = new TLegend(0.6, 0.7, 0.9, 0.9);
  legendH->SetFillColor(kWhite);
  legendH->SetTextSize(0.03);


  medium->GetXaxis()->SetTitle("m_{N} GeV");
  medium->GetYaxis()->SetTitle("ID efficiency");
  
  medium->SetMarkerColor(kRed);
  tight->SetMarkerColor(kRed);
  tight_dxy10->SetMarkerColor(kRed);
  tight_anal->SetMarkerColor(kRed);
  
  medium->SetMarkerStyle(20.);
  tight->SetMarkerStyle(21.);
  tight_dxy10->SetMarkerStyle(22.);
  tight_anal->SetMarkerStyle(23.);

  legendH->AddEntry(medium, "medium ID", "p");
  legendH->AddEntry(tight, "tight ID", "p");
  legendH->AddEntry(tight_dxy10, "tight+ dxy ", "p");
  legendH->AddEntry(tight_anal, "tight+ dxy+ iso ", "p");



  medium->Draw("p");
  tight->Draw("psame");
  tight_dxy10->Draw("psame");
  tight_anal->Draw("psame");


  legendH->Draw();


  TGraphAsymmErrors * g = new TGraphAsymmErrors(heff);
  g->SetLineWidth(2.0);
  g->SetMarkerSize(2.);
  //  g->Draw( "9pXsame" );
  
  
  CMS_lumi( c1, 2, 11 );
  c1->Update();
  c1->RedrawAxis();
  
  
  c1->SaveAs(("/home/jalmond/WebPlots/PreApproval/SignalPlots/SignalEff_presel_med_tight_pt15.pdf" ));
}

void MakeNPunziLowMass_lowmasscuts(){
  
  setTDRStyle();
  gStyle->SetPalette(1);

  TH1F* tight_anal1 = makehist("PreSelection_lowmass");
  TH1F* tight_anal2 = makehist("PreSelection_lowmass2");
  TH1F* tight_anal3 = makehist("PreSelection_lowmass3");


  TLegend* legendH = new TLegend(0.6, 0.7, 0.9, 0.9);
  legendH->SetFillColor(kWhite);
  legendH->SetTextSize(0.03);

  
  tight_anal3->GetXaxis()->SetTitle("m_{N} GeV");
  tight_anal3->GetYaxis()->SetTitle("ID efficiency");
 
  tight_anal1->SetMarkerColor(kRed);
  tight_anal1->SetMarkerStyle(20.);
  tight_anal2->SetMarkerColor(kRed);
  tight_anal2->SetMarkerStyle(21.);
  tight_anal3->SetMarkerColor(kRed);
  tight_anal3->SetMarkerStyle(22.);






  legendH->AddEntry(tight_anal1, "LowMass", "p");
  legendH->AddEntry(tight_anal2, "LowMass + m(ee) < 60", "p");
  legendH->AddEntry(tight_anal3, "LowMass + m(eejj) < 155", "p");

  tight_anal3->GetYaxis()->SetRangeUser(0., 0.001);
  tight_anal3->Draw("p");
  tight_anal1->Draw("psame");
  tight_anal2->Draw("psame");


  legendH->Draw();

  TGraphAsymmErrors * g = new TGraphAsymmErrors(heff);
  g->SetLineWidth(2.0);
  g->SetMarkerSize(2.);
  //  g->Draw( "9pXsame" );
  
  
  CMS_lumi( c1, 2, 11 );
  c1->Update();
  c1->RedrawAxis();
  
  
  c1->SaveAs(("/home/jalmond/WebPlots/PreApproval/SignalPlots/Punzi_presel_lowmasscuts.pdf" ));
  return;
  
}

void Histograms::MakeL1GEfficiency(TH1F * num, TH1F * denom, const char * hname) {
  // Make an efficiency graph with the "L1G" style
  TGraphAsymmErrors * eff = new TGraphAsymmErrors( num, denom);
  eff->SetMarkerStyle(24);
  eff->SetMarkerColor(kBlue);
  eff->SetMarkerSize(1.5);
  eff->SetLineColor(kBlack);
  m_efficiencies[std::string( hname )] = eff;
  
}

TGraphAsymmErrors *fromGEpGMptoQ4GEp(TGraph *ogr) {
	if (!ogr) return 0;
	TGraphAsymmErrors *ogrE = dynamic_cast<TGraphAsymmErrors*>(ogr);
	int npts = ogr->GetN();
	Double_t x,y;
	Double_t Q2;
	Double_t Gd;
	TGraphAsymmErrors *gr = new TGraphAsymmErrors(npts);
	ogr->TAttMarker::Copy(*gr);
	ogr->TAttLine::Copy(*gr);
	ogr->TAttFill::Copy(*gr);
	gr->SetName(ogr->GetName());

	for (int i=0; i<npts; i++) {
		Double_t ex1=0,ex2=0,ey1=0,ey2=0;
		ogr->GetPoint(i,x,y);
		Q2=x;
		if (ogrE) {
			ex1 = ogrE->GetErrorXlow(i);
			ex2 = ogrE->GetErrorXhigh(i);
			ey1 = ogrE->GetErrorYlow(i);
			ey2 = ogrE->GetErrorYhigh(i);
		}
		Gd = pow( 1.0 + Q2/0.71, -2.0);
		y *= pow(Q2,2.0)*Gd;
		ey1 *= pow(Q2,2.0)*Gd;
		ey2 *= pow(Q2,2.0)*Gd;
		gr->SetPoint(i,x,y);
		gr->SetPointError(i,ex1,ex2,ey1,ey2);
	}
	return gr;
}

TGraphAsymmErrors *fromGEntoGEnGd(TGraph *ogr) {
// take a curve (with errors) of GEn and make it into (GEn/Gd)
  if (!ogr) return 0;
  TGraphAsymmErrors *ogrE = dynamic_cast<TGraphAsymmErrors*>(ogr);
  int npts = ogr->GetN();
  Double_t x,y;
  Double_t Q2;
  Double_t Gd;
  TGraphAsymmErrors *gr = new TGraphAsymmErrors(npts);
  ogr->TAttMarker::Copy(*gr);
  ogr->TAttLine::Copy(*gr);
  ogr->TAttFill::Copy(*gr);
  gr->SetName(ogr->GetName());

  for (int i=0; i<npts; i++) {
    Double_t ex1=0,ex2=0,ey1=0,ey2=0;
    ogr->GetPoint(i,x,y);
    Q2=x;
    if (ogrE) {
      ex1 = ogrE->GetErrorXlow(i);
      ex2 = ogrE->GetErrorXhigh(i);
      ey1 = ogrE->GetErrorYlow(i);
      ey2 = ogrE->GetErrorYhigh(i);
    }
    Gd = gdipole(&Q2,0);
    y /= Gd;
    ey1 /= Gd;
    ey2 /= Gd;
    gr->SetPoint(i,x,y);
    gr->SetPointError(i,ex1,ex2,ey1,ey2);
  }
  return gr;
}

//*************************************************************************************************
//
//*************************************************************************************************
TGraphAsymmErrors* MakeCurrentWPSigEffVsCutValueGraph(TH1F* signalHist, string name, Double_t myCutValue ) {
  //Make Met Plots
  const UInt_t nPoints = signalHist->GetXaxis()->GetNbins();
  double cutValue[1] = {0};
  double cutValueErr[1] = {0};
  double SigEff[1] = {0};
  double SigEffErrLow[1] = {0};
  double SigEffErrHigh[1] = {0};
  double NSigTotal = 0;
  
  Double_t effDiff = 9999;

  for (UInt_t q=0; q < nPoints+2; ++q) {
    NSigTotal += signalHist->GetBinContent(q);
  }

  for(UInt_t b=0; b < nPoints; ++b) {
    Double_t nsig = 0;
    for (UInt_t q=b; q < nPoints+2; ++q) {
      nsig += signalHist->GetBinContent(q);
    }

    Double_t ratio;
    Double_t n1 = 0;
    Double_t n2 = 0;

    n1 = TMath::Nint(nsig);
    n2 = TMath::Nint(NSigTotal);
    ratio = n1/n2;
    
      cout << myCutValue << " : " << signalHist->GetXaxis()->GetBinCenter(b) << " , " << cutValue[0] << endl;
    if (fabs(myCutValue - signalHist->GetXaxis()->GetBinCenter(b)) < fabs(myCutValue - cutValue[0])) {
      SigEff[0] = ratio;
      SigEffErrLow[0] = 0;
      SigEffErrHigh[0] = 0;
      cutValue[0] = signalHist->GetXaxis()->GetBinCenter(b);
      cutValueErr[0] = 0;
    }
  }

//   cout << "Final: " << cutValue[0] << " , " << SigEff[0] << endl;

  TGraphAsymmErrors *tmpSigEffVsCut = new TGraphAsymmErrors (1, cutValue, SigEff, cutValueErr, cutValueErr, SigEffErrLow, SigEffErrHigh  );
  tmpSigEffVsCut->SetName(name.c_str());
  tmpSigEffVsCut->SetTitle("");
  tmpSigEffVsCut->GetXaxis()->SetTitle("Cut Value");
  tmpSigEffVsCut->GetYaxis()->SetTitle("Efficiency");
  tmpSigEffVsCut->GetYaxis()->SetTitleOffset(1.1);
  tmpSigEffVsCut->GetXaxis()->SetTitleOffset(1.05);
  tmpSigEffVsCut->SetMarkerColor(kBlack);
  tmpSigEffVsCut->SetLineColor(kBlack);
  tmpSigEffVsCut->SetMarkerSize(1.5);

  return tmpSigEffVsCut;
}

TGraphAsymmErrors* fitTools::getGraphPoissonErrors( TH1D* histo, bool drawZeros, const std::string xerrType, float nSigma ) {

  
  TGraphAsymmErrors* graph = new TGraphAsymmErrors(0);

  for( unsigned iBin=1; iBin<(histo->GetXaxis()->GetNbins()+1); ++iBin ) {

    int y; // these are data histograms, so y has to be integer
    double x, xerr, yerrplus, yerrminus;
    //xerr = 0.; //no xerr for now (maybe binwidth / sqrt(12)?)

    x = histo->GetBinCenter(iBin);
    if( xerrType=="0" )
      xerr = 0.;
    else if( xerrType=="binWidth" )
      xerr = histo->GetBinWidth(iBin)/2.;
    else if( xerrType=="sqrt12" )
      xerr = histo->GetBinWidth(iBin)/sqrt(12.);
    else {
      std::cout << "Unkown xerrType '" << xerrType << "'. Setting to bin width." << std::endl;
      xerr = histo->GetBinWidth(iBin);
    }
    y = (int)histo->GetBinContent(iBin);

    if( y==0 && !drawZeros ) continue;
       
    double ym, yp;
    RooHistError::instance().getPoissonInterval(y,ym,yp,nSigma);

    yerrplus = yp - y;
    yerrminus = y - ym;

/*
    // and now poissonian errors (bayes flat prior):
    if( y==0 ) {
      yerrminus = 0.;
      yerrplus = 0.5*TMath::ChisquareQuantile(cl, 2.*(y+1.) );
    } else {  
      //float lowL = 0.5*TMath::ChisquareQuantile(1.-cl/2., 2.*y);
      //float upL = 0.5*TMath::ChisquareQuantile(cl/2., 2.*(y+1.) );
      float lowL = 0.5*TMath::ChisquareQuantile(1.-cl, 2.*y);
      float upL = 0.5*TMath::ChisquareQuantile(cl, 2.*(y+1.) );
      yerrminus = y - lowL;
      yerrplus = upL - y;
    }
*/
    int thisPoint = graph->GetN();
    graph->SetPoint( thisPoint, x, y );
    graph->SetPointError( thisPoint, xerr, xerr, yerrminus, yerrplus );

  }

  return graph;

}