C++ TGraphAsymmErrors::BayesDivide方法代码示例

void printHisto( TChain *data , TChain *mc , TCut num , TCut denom , char* var , int nbins , float xmin , float xmax , char* xtitle , char* ytitle){

  TH1F* hpass_data = new TH1F(Form("hpass_data_%i",iplot),Form("hpass_data_%i",iplot),nbins,xmin,xmax);
  TH1F* hall_data  = new TH1F(Form("hall_data_%i" ,iplot),Form("hall_data_%i" ,iplot),nbins,xmin,xmax);
  TH1F* hpass_mc   = new TH1F(Form("hpass_mc_%i"  ,iplot),Form("hpass_mc_%i"  ,iplot),nbins,xmin,xmax);
  TH1F* hall_mc    = new TH1F(Form("hall_mc_%i"   ,iplot),Form("hall_mc_%i"   ,iplot),nbins,xmin,xmax);

  TCanvas *can = new TCanvas(Form("can_%i",iplot),Form("can_%i",iplot),600,600);
  can->cd();
  
  data->Draw(Form("min(%s,%f)>>hpass_data_%i"  , var,xmax-0.0001,iplot),denom+num);
  data->Draw(Form("min(%s,%f)>>hall_data_%i"   , var,xmax-0.0001,iplot),denom);
  mc->Draw  (Form("min(%s,%f)>>hpass_mc_%i"    , var,xmax-0.0001,iplot),denom+num);
  mc->Draw  (Form("min(%s,%f)>>hall_mc_%i"     , var,xmax-0.0001,iplot),denom);

  TGraphAsymmErrors *grdata = new TGraphAsymmErrors();
  grdata->BayesDivide(hpass_data,hall_data);

  TGraphAsymmErrors *grmc = new TGraphAsymmErrors();
  grmc->BayesDivide(hpass_mc,hall_mc);

  cout << "data all  " << hall_data->GetBinContent(8) << endl;
  cout << "data pass " << hpass_data->GetBinContent(8) << endl;
  cout << "data eff  " << hpass_data->GetBinContent(8) / hall_data->GetBinContent(8) << endl;

  Double_t x;
  Double_t y;
  grdata->GetPoint(7,x,y);
  cout << "data eff2 " << y << endl;

  gPad->SetGridx();
  gPad->SetGridy();

  grdata->SetMarkerColor(2);
  grdata->SetLineColor(2);
  grmc->SetMarkerColor(4);
  grmc->SetLineColor(4);
  grmc->SetMarkerStyle(25);

  grdata->GetXaxis()->SetTitle(xtitle);
  grdata->GetYaxis()->SetTitle(ytitle);
  grdata->Draw("AP");
  grmc->Draw("sameP");

  TLegend *leg = new TLegend(0.5,0.2,0.7,0.4);
  leg->AddEntry(grdata ,"data","lp");
  leg->AddEntry(grmc   ,"mc"  ,"lp");
  leg->SetBorderSize(1);
  leg->SetFillColor(0);
  leg->Draw();

  iplot ++;
}

void plotEvtSelEff() {

  gStyle->SetOptStat(0);

  TFile *f = new TFile("../test/ROOTupleMC_HighPurity.root");
  TH1F *hNSD = (TH1F*) f->Get("preTrgAna/hGenMultNSD");
  TH1F *hNSDtrg = (TH1F*) f->Get("postTrgAna/hGenMultNSD");
  TH1F *hNSDevt = (TH1F*) f->Get("postEvtSelAna/hGenMultNSD");
  TH1F *hNSDvtx = (TH1F*) f->Get("postVtxAna/hGenMultNSD");
  TH1F *hNSDtrkvtx = (TH1F*) f->Get("postTrkVtxAna/hGenMultNSD");

  TGraphAsymmErrors *gNSDtrg = new TGraphAsymmErrors();
  TGraphAsymmErrors *gNSDevt = new TGraphAsymmErrors();
  TGraphAsymmErrors *gNSDvtx = new TGraphAsymmErrors();
  TGraphAsymmErrors *gNSDtrkvtx = new TGraphAsymmErrors();

  gNSDtrg->BayesDivide(hNSDtrg,hNSD);
  gNSDevt->BayesDivide(hNSDevt,hNSD);
  gNSDvtx->BayesDivide(hNSDvtx,hNSD);
  gNSDtrkvtx->BayesDivide(hNSDtrkvtx,hNSD);

  TH1F *dum1 = new TH1F("dum1",";Charged-particle multiplicity;Fraction of events",100,0,100);
  dum1->SetMaximum(0.05);
  dum1->GetYaxis()->SetTitleOffset(1.8);
  TH1F *dum2 = new TH1F("dum2",";Charged-particle multiplicity;Selection efficiency",100,0,60);

  TCanvas *c1 = new TCanvas("c1","Event Selection",900,500);
  c1->Divide(2,1);

  c1->cd(1);
  dum1->Draw();
  hNSDtrkvtx->Sumw2();
  hNSDtrkvtx->Scale(1./hNSDtrkvtx->GetEntries());
  hNSDtrkvtx->SetMarkerStyle(24);
  hNSDtrkvtx->Draw("pzsame");
  hNSDvtx->Sumw2();
  hNSDvtx->Scale(1./hNSDvtx->GetEntries());
  hNSDvtx->SetMarkerStyle(20);
  hNSDvtx->Draw("pzsame");

  c1->cd(2);
  dum2->Draw();
  gNSDtrkvtx->SetMarkerStyle(24);
  gNSDtrkvtx->Draw("pzsame");
  gNSDvtx->SetMarkerStyle(20);
  gNSDvtx->Draw("pzsame");
}

// Take a tree, a variable and calculate the efficiency
TGraphAsymmErrors* getEfficiency(TTree *t, char *variable, TCut preselection, TCut cut, int nBin, Float_t *bins)
{
   static int count = 0;
   count++;
   TH1D *hPass = new TH1D (Form("hPass%d",count),"",nBin,bins);
   TH1D *hAll = new TH1D (Form("hAll%d",count),"",nBin,bins);
   t->Draw(Form("%s>>hAll%d",variable,count),preselection);
   t->Draw(Form("%s>>hPass%d",variable,count),preselection&&cut);

   TGraphAsymmErrors *g = new TGraphAsymmErrors;
   g->BayesDivide(hPass,hAll);
   return g;
}

//____________________________________________________________________________________
TGraphAsymmErrors *calcEfficiency(TH1* h1, TH1* h2)
{

  TH1 *phUp   = (TH1 *)h2->Clone("phUp");
  TH1 *phDown = (TH1 *)h1->Clone("phDown");
  
  phUp->SetDirectory(0);
  phDown->SetDirectory(0);
  
  TGraphAsymmErrors *pgEfficiency = new TGraphAsymmErrors();
  pgEfficiency->BayesDivide(phUp,phDown,"");
  return pgEfficiency;
}

TGraphAsymmErrors *calcEffpythia(TH1* h1, TH1* hCut,double *npart)
{
   TGraphAsymmErrors *gEfficiency = new TGraphAsymmErrors();
   gEfficiency->BayesDivide(hCut,h1);
   cout <<gEfficiency->GetN()<<endl;
   for (int i=0;i<gEfficiency->GetN();i++)
   {
      double x,y;
      gEfficiency->GetPoint(i,x,y);
      double errYL = gEfficiency->GetErrorYlow(i);
      double errYH = gEfficiency->GetErrorYhigh(i);
      gEfficiency->SetPointError(i,0,0,errYL,errYH);
      gEfficiency->SetPoint(i,2,y);
      cout<<" Setting PYTHIA point to npart=2 by hand, Yen-Jie please check -Matt "<<endl;
   }
   return gEfficiency;
}

TGraphAsymmErrors *calcEff(TH1* h1, TH1* h2,double *npart)
{
   TGraphAsymmErrors *gEfficiency = new TGraphAsymmErrors();
   gEfficiency->BayesDivide(h2,h1);
   cout <<gEfficiency->GetN()<<endl;
   for (int i=0;i<gEfficiency->GetN();i++)
   {
      double x,y;
      gEfficiency->GetPoint(i,x,y);
      double errYL = gEfficiency->GetErrorYlow(i);
      double errYH = gEfficiency->GetErrorYhigh(i);
      gEfficiency->SetPointError(i,0,0,errYL,errYH);
      gEfficiency->SetPoint(i,npart[h1->FindBin(x)-1],y);
      cout <<x<<" "<<h1->FindBin(x)<<" "<<npart[h1->FindBin(x)-1]<<endl;
   }
   return gEfficiency;
}

void LLGAnalysis::MakeEfficiencyPlot( TH1D hpass, TH1D htotal, TCanvas *c, string triggerName ) {
    
    TGraphAsymmErrors geff;
    geff.BayesDivide( &hpass, &htotal );
    geff.GetXaxis()->SetTitle( hpass.GetXaxis()->GetTitle() );
    string ytitle = "#varepsilon (" + triggerName + ")";
    geff.GetYaxis()->SetTitle( ytitle.c_str() );
    string efftitle = "efficiency_" + triggerName;
    geff.SetNameTitle(efftitle.c_str(), efftitle.c_str());
    geff.SetMarkerColor(kBlue);
    geff.Draw("APZ"); 
    for( vector<string>::iterator itr_f = _plotFormats.begin(); itr_f != _plotFormats.end(); ++itr_f ) {
        string thisPlotName = efftitle + (*itr_f);
        c->Print( thisPlotName.c_str() );
    }

}

float getBinomialError( float num , float den ){

  TGraphAsymmErrors* grtemp = new TGraphAsymmErrors();

  TH1F* hnum = new TH1F("hnum","",1,0,1);
  TH1F* hden = new TH1F("hden","",1,0,1);

  hnum->SetBinContent(1,num);
  hden->SetBinContent(1,den);

  grtemp->BayesDivide(hnum,hden);

  float err = 0.5 * ( grtemp->GetErrorYlow(0) + grtemp->GetErrorYhigh(0) );

  delete hnum;
  delete hden;

  return err;
}

TGraphAsymmErrors* getEfficiencyGraph( TChain* ch , TCut num , TCut denom , char* var , int nbins , float xmin , float xmax , char* xtitle , char* ytitle){


  //float ptbin[] = {10., 15., 20., 30., 40., 50., 7000.};
  float ptbin[] = { 20.0 , 25.0 , 30.0 , 35.0 , 40.0 , 45.0 , 50.0 , 55.0 , 60.0 , 70.0 , 80.0 , 90.0 , 100.0 , 150.0 , 200.0 , 250.0 , 300.0};
  int   nptbin  = 16;

  TH1F* hpass   = new TH1F(Form("hpass_%i",iplot),Form("hpass_%i",iplot),nptbin,ptbin);
  TH1F* hall    = new TH1F(Form("hall_%i" ,iplot),Form("hall_%i" ,iplot),nptbin,ptbin);

  //TH1F* hpass   = new TH1F(Form("hpass_%i",iplot),Form("hpass_%i",iplot),nbins,xmin,xmax);
  //TH1F* hall    = new TH1F(Form("hall_%i" ,iplot),Form("hall_%i" ,iplot),nbins,xmin,xmax);

  TCanvas *ctemp = new TCanvas();
  ctemp->cd();  
  ch->Draw(Form("min(%s,%f)>>hpass_%i"  , var,xmax-0.0001,iplot),denom+num);
  ch->Draw(Form("min(%s,%f)>>hall_%i"   , var,xmax-0.0001,iplot),denom);
  delete ctemp;

  /*
  cout << "PASS " << hpass->Integral() << endl;
  cout << "ALL  " << hall->Integral() << endl;
  cout << "EFF  " << (float) hpass->Integral() / (float) hall->Integral() << endl;

  TCanvas *c1 = new TCanvas();
  c1->cd();
  hall->Draw();
  hpass->SetLineColor(2);
  hpass->Draw("samehist");
  */

  TGraphAsymmErrors *gr = new TGraphAsymmErrors();
  gr->BayesDivide(hpass,hall);

  gr->GetXaxis()->SetTitle(xtitle);
  gr->GetYaxis()->SetTitle(ytitle);

  iplot++;

  return gr;
}

void jetEfficiencies_V21( bool printplots = false ){

  gStyle->SetOptFit(0);

  TFile* f     = TFile::Open("../output/V00-02-21/wz_summer11_madgraph.root");
  TFile* ftt   = TFile::Open("../output/V00-02-21/ttbar.root");

  //TChain *ch = new TChain("T1");
  //ch->Add("../output/V00-02-20/wzsms_gen_baby.root");
  //TFile* f = TFile::Open("../output/V00-02-20/wzsms.root");
  //ch->Add("../output/V00-02-21/wz_summer11_madgraph_gen_baby.root");
  //ch->Add("../output/V00-02-20/ttbar_gen_baby.root");
  //TFile* f = TFile::Open("../output/V00-02-20/ttbar.root");

  //---------------------------------------
  // jets
  //---------------------------------------

  TH1F* hjetpt_all      = (TH1F*) f->Get("hjetpt_all");
  TH1F* hjetpt_q_all    = (TH1F*) f->Get("hjetpt_q_all");
  TH1F* hjetpt_c_all    = (TH1F*) f->Get("hjetpt_c_all");
  TH1F* hjetpt_b_all    = (TH1F*) ftt->Get("hjetpt_b_all");
  TH1F* hjetpt_g_all    = (TH1F*) f->Get("hjetpt_g_all");

  TH1F* hjetpt_pass     = (TH1F*) f->Get("hjetpt_pass");
  TH1F* hjetpt_q_pass   = (TH1F*) f->Get("hjetpt_q_pass");
  TH1F* hjetpt_c_pass   = (TH1F*) f->Get("hjetpt_c_pass");
  TH1F* hjetpt_b_pass   = (TH1F*) ftt->Get("hjetpt_b_pass");
  TH1F* hjetpt_g_pass   = (TH1F*) f->Get("hjetpt_g_pass");

  TH1F* hbtag_q_all     = (TH1F*) f->Get("hbtag_q_all");
  TH1F* hbtag_c_all     = (TH1F*) f->Get("hbtag_c_all");
  TH1F* hbtag_b_all     = (TH1F*) ftt->Get("hbtag_b_all");
  TH1F* hbtag_g_all     = (TH1F*) f->Get("hbtag_g_all");

  TH1F* hbtag_q_passL   = (TH1F*) f->Get("hbtag_q_passL");
  TH1F* hbtag_c_passL   = (TH1F*) f->Get("hbtag_c_passL");
  TH1F* hbtag_b_passL   = (TH1F*) ftt->Get("hbtag_b_passL");
  TH1F* hbtag_g_passL   = (TH1F*) f->Get("hbtag_g_passL");

  TH1F* hbtag_q_passM   = (TH1F*) f->Get("hbtag_q_passM");
  TH1F* hbtag_c_passM   = (TH1F*) f->Get("hbtag_c_passM");
  TH1F* hbtag_b_passM   = (TH1F*) ftt->Get("hbtag_b_passM");
  TH1F* hbtag_g_passM   = (TH1F*) f->Get("hbtag_g_passM");

  int rebin = 5;

  hjetpt_all->Rebin(rebin);
  hjetpt_q_all->Rebin(rebin);
  hjetpt_b_all->Rebin(rebin);
  hjetpt_c_all->Rebin(rebin);
  hjetpt_g_all->Rebin(rebin);

  hjetpt_pass->Rebin(rebin);
  hjetpt_q_pass->Rebin(rebin);
  hjetpt_b_pass->Rebin(rebin);
  hjetpt_c_pass->Rebin(rebin);
  hjetpt_g_pass->Rebin(rebin);

  int rebinb = 5;

  hbtag_q_all->Rebin(rebinb);
  hbtag_c_all->Rebin(rebinb);
  hbtag_b_all->Rebin(rebinb);
  hbtag_g_all->Rebin(rebinb);

  hbtag_q_passL->Rebin(rebinb);
  hbtag_c_passL->Rebin(rebinb);
  hbtag_b_passL->Rebin(rebinb);
  hbtag_g_passL->Rebin(rebinb);

  hbtag_q_passM->Rebin(rebinb);
  hbtag_c_passM->Rebin(rebinb);
  hbtag_b_passM->Rebin(rebinb);
  hbtag_g_passM->Rebin(rebinb);

  // cout << "Total entries   " << hjetpt_all->GetEntries()   << endl;
  // cout << "Total q entries " << hjetpt_q_all->GetEntries() << endl;
  // cout << "Total c entries " << hjetpt_b_all->GetEntries() << endl;
  // cout << "Total b entries " << hjetpt_c_all->GetEntries() << endl;
  // cout << "Total g entries " << hjetpt_g_all->GetEntries() << endl;

  // cout << "Total jet entries   " << hjetpt_pass04->GetEntries()   << endl;
  // cout << "Total jet q entries " << hbtag_q_all->GetEntries() << endl;
  // cout << "Total jet c entries " << hbtag_b_all->GetEntries() << endl;
  // cout << "Total jet b entries " << hbtag_c_all->GetEntries() << endl;
  // cout << "Total jet g entries " << hbtag_g_all->GetEntries() << endl;

  //---------------------------------------
  // jet reco
  //---------------------------------------

  TGraphAsymmErrors* grjet = new TGraphAsymmErrors();
  grjet->BayesDivide(hjetpt_pass,hjetpt_all);

  TGraphAsymmErrors* grjet_q = new TGraphAsymmErrors();
  grjet_q->BayesDivide(hjetpt_q_pass,hjetpt_q_all);

  TGraphAsymmErrors* grjet_c = new TGraphAsymmErrors();
  grjet_c->BayesDivide(hjetpt_c_pass,hjetpt_c_all);
//.........这里部分代码省略.........

//.........这里部分代码省略.........

  for (unsigned int i=1; i<=Nbins/nRebin;i++){
    // All Muon
    if (denMB->GetBinContent(i)>0.) {
      fakeMB->SetBinContent(i,double(numMB->GetBinContent(i)/denMB->GetBinContent(i)));

      xMB[i-1]=fakeMB->GetBinCenter(i);
      yMB[i-1]=fakeMB->GetBinContent(i);
      exlMB[i-1]=0.;
      exhMB[i-1]=0.;
      //eyhMB[i-1]=+0.5 + sqrt(fakeMB->GetBinContent(i)+0.25);
      //eylMB[i-1]=-0.5 + sqrt(fakeMB->GetBinContent(i)+0.25); 
      
      if (double(numMB->GetBinContent(i)/denMB->GetBinContent(i))<1.){
	fakeMB->SetBinError(i, sqrt( double(numMB->GetBinContent(i)/denMB->GetBinContent(i))*(1-double(numMB->GetBinContent(i)/denMB->GetBinContent(i)))/double(denMB->GetBinContent(i)) ) );}
      else fakeMB->SetBinError(i, sqrt( double(numMB->GetBinContent(i)/pow(denMB->GetBinContent(i),2)+1./pow(denMB->GetBinContent(i),3))));
      eyhMB[i-1]=fakeMB->GetBinError(i)/2.;
      eylMB[i-1]=fakeMB->GetBinError(i)/2.;      
    }
    else {
      fakeMB->SetBinContent(i,0.);
      fakeMB->SetBinError(i,0.);
      xMB[i-1] = 0.;
      eylMB[i-1] = 0.;
      eyhMB[i-1] = 0.;
      exlMB[i-1]=0.;
      exhMB[i-1]=0.;
    }

  }
 
  //TGraphAsymmErrors *grMB = new TGraphAsymmErrors(Nbins/nRebin,xMB,yMB,exlMB,exhMB,eylMB,eyhMB);
  TGraphAsymmErrors *grMB = new TGraphAsymmErrors();
  grMB->BayesDivide(numMB,denMB);

  //TGraphAsymmErrors *grMB = new TGraphAsymmErrors(numMB,denMB);
  grMB->SetMarkerColor(2);
  grMB->SetMarkerStyle(20);
  grMB->SetMarkerSize(0.95);

  //TF1* fn1 = new TF1("fn1","pol1",10,100);
  //fn1->SetLineColor(1);
  //grMB->Fit("fn1","R,smaes");

  //TF1* fn2 = new TF1("fn2",FittingFunc1,100,2000.0,3);
  //fn2->SetLineColor(1);
  //grMB->Fit("fn2","R");


  // ME
  Float_t xME[arraysize[0]],yME[arraysize[0]],exlME[arraysize[0]],exhME[arraysize[0]],eylME[arraysize[0]],eyhME[arraysize[0]];

  for (unsigned int i=1; i<=Nbins/nRebin;i++){
    // All Muon                                                                                                                                                                            
    if (denME->GetBinContent(i)>0.) {
      fakeME->SetBinContent(i,double(numME->GetBinContent(i)/denME->GetBinContent(i)));

      xME[i-1]=fakeME->GetBinCenter(i);
      yME[i-1]=fakeME->GetBinContent(i);
      exlME[i-1]=0.;
      exhME[i-1]=0.;
      //eyhME[i-1]=+0.5 + sqrt(fakeME->GetBinContent(i)+0.25);                                                                                                                               
      //eylME[i-1]=-0.5 + sqrt(fakeME->GetBinContent(i)+0.25);                                                                                

      if (double(numME->GetBinContent(i)/denME->GetBinContent(i))<1.){
        fakeME->SetBinError(i, sqrt( double(numME->GetBinContent(i)/denME->GetBinContent(i))*(1-double(numME->GetBinContent(i)/denME->GetBinContent(i)))/double(denME->GetBinContent(i)) ) );

//__________________________________________________________________________
TGraphAsymmErrors *calcEff(TH1* h1, TH1* h2)
{
  TGraphAsymmErrors *gEfficiency = new TGraphAsymmErrors();
  gEfficiency->BayesDivide(h1,h2);
  return gEfficiency;
}

void Draw3paneleff() {
    gStyle->SetOptStat(kFALSE);
    TString outG="G1.root";
    TCanvas *c1  = new TCanvas("","",1000,400);
    c1->Divide(3);
    gStyle->SetErrorX(0);
    double hfbin[]= {0,1,2,3,4,6,8,10,13,16,20,25,30,40,55,70,90};
    int nhfbin = 16;
    TString name, dirname;
    TFile *fDSeff = TFile::Open("/scratch/xuq7/Centrality/pPbHijing_EffCorr_forNBD.root");
    TFile *ftreff = TFile::Open("/scratch/xuq7/Centrality/pPbHist_Hijing_TrandEs.root");
    TH1D* hbef = (TH1D*)ftreff->Get("hHFEnergy4");
    TH1D* rehbef = (TH1D*)hbef->Rebin(nhfbin,"rehHFEnergy4",hfbin);
    TH1D* haft = (TH1D*)ftreff->Get("hHFEnergy4_tr");
    TH1D* rehaft = (TH1D*)haft->Rebin(nhfbin,"rehHFEnergy4_tr",hfbin);
    TGraphAsymmErrors *gtreff = new TGraphAsymmErrors();
    gtreff->BayesDivide(rehaft,rehbef);

    TGraphAsymmErrors *geff = (TGraphAsymmErrors*)fDSeff->Get("regEffHF4");
    for(int i=0; i<geff->GetN(); i++) {
        geff->SetPointEXlow(i,0);
        geff->SetPointEXhigh(i,0);
        gtreff->SetPointEXlow(i,0);
        gtreff->SetPointEXhigh(i,0);
    }

    geff->SetTitle("");
    geff->SetMarkerStyle(33);
    geff->SetMarkerColor(5);
    geff->SetMarkerSize(1.3);

    gtreff->SetTitle("");
    gtreff->SetMarkerStyle(21);
    gtreff->SetMarkerColor(6);
    gtreff->SetMarkerSize(1.3);

    for(int cenvar=0; cenvar<3; cenvar++) {
        TFile *f = TFile::Open(Form("%s/%s",cenvardir[cenvar].Data(),outG.Data()));
        c1->cd(cenvar+1);
        TObjString* dataname = (TObjString*)f->Get(Form("dataname"));
        TObjString* histoname = (TObjString*)f->Get(Form("histoname"));
        TFile *fdata = TFile::Open(dataname->GetString());
        TH1D *histo_obs = (TH1D*)fdata->Get(histoname->GetString());
        histo_obs->Sumw2();
        TH1D* hFrame = histo_obs->Clone();
        hFrame->Reset();
        hFrame->SetTitle("");
        hFrame->SetMaximum(1.2);
        hFrame->SetMinimum(0);
        if(cenvar==0)hFrame->GetYaxis()->SetTitle("ratio");
        else hFrame->GetYaxis()->SetTitle("");
        hFrame->GetXaxis()->SetRangeUser(0,100);
        hFrame->GetXaxis()->CenterTitle(0);
        hFrame->GetYaxis()->CenterTitle(1);
        hFrame->GetYaxis()->SetTitleOffset(1.1);
        hFrame->GetXaxis()->SetTitleOffset(1.1);
        hFrame->GetXaxis()->SetTitleSize(0.056);
        hFrame->GetYaxis()->SetTitleSize(0.056);
        hFrame->GetXaxis()->SetLabelSize(0.05);
        hFrame->GetYaxis()->SetLabelSize(0.05);

        int binnum = histo_obs->GetNbinsX();
        double Minx = histo_obs->GetXaxis()->GetXmin();
        double Maxx = histo_obs->GetXaxis()->GetXmax();
        double binsize = (Double_t)(Maxx-Minx)/binnum;
        hFrame->Draw();
        TLegend *leg = new TLegend(0.40, 0.2, 0.65, 0.4);

        for(int sth=0; sth<3; sth++) {
            int Gth=0;
            if(sth==0) {
                dirname = "std";
            }
            else if(sth==1) {
                dirname ="Gri055";
            }
            else {
                dirname ="Gri101";
            }
            hFrame->GetXaxis()->SetTitle(cenvariable[cenvar]);

            if(Gth==0)
                name = "G0";
            else if(Gth<nGlau)
                name = Form("Glau_%d",Gth);
            else
                name = Form("bin_%d",Gth-nGlau+1);
            TObjString* Glaubername = (TObjString*)f->Get(Form("%s/%s/Glaubername",dirname.Data(),name.Data()));
            TVectorD* xmin = (TVectorD*)f->Get(Form("%s/%s/xmin",dirname.Data(),name.Data()));
            TVectorD* xmax = (TVectorD*)f->Get(Form("%s/%s/xmax",dirname.Data(),name.Data()));
            TVectorD* mubest = (TVectorD*)f->Get(Form("%s/%s/mubest",dirname.Data(),name.Data()));
            TVectorD* kbest = (TVectorD*)f->Get(Form("%s/%s/kbest",dirname.Data(),name.Data()));
            TVectorD* Ndf = (TVectorD*)f->Get(Form("%s/%s/Ndf",dirname.Data(),name.Data()));
            TVectorD* chis = (TVectorD*)f->Get(Form("%s/%s/chis",dirname.Data(),name.Data()));
            TVectorD *kpoint = (TVectorD*)f->Get(Form("%s/%s/kpoint",dirname.Data(),name.Data()));
            int xbinmin=(int)(((*xmin)[0]-Minx)/binsize);
            int xbinmax=(int)(((*xmax)[0]-Minx)/binsize);
            TFile *fGlauber = TFile::Open(Glaubername->GetString());
            TH1D *histo_exp = new TH1D("histo_exp","Simulated distribution;Multiplicity;Event Fraction",binnum,Minx,Maxx);
            histo_exp->Sumw2();
//.........这里部分代码省略.........

void plotHist3D() {

  RiceStyle();

  gStyle->SetOptStat(0);

  TFile *f = new TFile("../rootfile/PbPb_eff_MC_v1.root");

  char ndir[256] = "HITrackCorrections";
  double ptmax = 300.;

  // sim-to-reco hists
  TH3F *hSim = (TH3F*) f->Get(Form("%s/hsim3D",ndir)); hSim->GetYaxis()->SetRangeUser(0.2,ptmax);
  TH3F *hEff = (TH3F*) f->Get(Form("%s/heff3D",ndir)); hEff->GetYaxis()->SetRangeUser(0.2,ptmax);
  TH3F *hMul = (TH3F*) f->Get(Form("%s/hmul3D",ndir)); hMul->GetYaxis()->SetRangeUser(0.2,ptmax);

  // reco-to-sim hists
  TH3F *hRec = (TH3F*) f->Get(Form("%s/hrec3D",ndir)); hRec->GetYaxis()->SetRangeUser(0.2,ptmax);
  TH3F *hFak = (TH3F*) f->Get(Form("%s/hfak3D",ndir)); hFak->GetYaxis()->SetRangeUser(0.2,ptmax);
  TH3F *hSec = (TH3F*) f->Get(Form("%s/hsec3D",ndir)); hSec->GetYaxis()->SetRangeUser(0.2,ptmax);

  // ratio histograms
  TH3F *rEff = (TH3F*) hEff->Clone("rEff");
  TH3F *rMul = (TH3F*) hMul->Clone("rMul");
  TH3F *rFak = (TH3F*) hFak->Clone("rFak");
  TH3F *rSec = (TH3F*) hSec->Clone("rSec");

  //---------------------------------------------
  //---------------------------------------------

  // find bins corresponding to projections for below
  int ptbin04=hSim->GetYaxis()->FindBin(0.41);
  int ptbin10=hSim->GetYaxis()->FindBin(1.01);
  int ptbins=hSim->GetYaxis()->GetNbins();

  int etabin24m=hSim->GetXaxis()->FindBin(-2.39);
  int etabin24p=hSim->GetXaxis()->FindBin(2.39);
  int etabin10m=hSim->GetXaxis()->FindBin(-0.99);
  int etabin10p=hSim->GetXaxis()->FindBin(0.99);

  int occbin0 = hSim->GetZaxis()->FindBin(0);
  int occbin50 = hSim->GetZaxis()->FindBin(99);
  int occbin51 = hSim->GetZaxis()->FindBin(100);
  int occbin100 = hSim->GetZaxis()->FindBin(199);

  cout << "etabin10m: " << etabin10m << " etabin10p: " << etabin10p << endl;
  cout << "etabin10m: " << etabin24m << " etabin10p: " << etabin24p << endl;
  cout << "occbin0: " << occbin0 << "occbin50: " << occbin50 << endl;

  // projected hists: pt > 1.0 GeV/c, cBin (0,50%)
  TH1D* hSimEta = (TH1D*) hSim->ProjectionX("hSimEta",ptbin10,ptbins,occbin0,occbin50,"e");
  TH1D* hEffEta = (TH1D*) hEff->ProjectionX("hEffEta",ptbin10,ptbins,occbin0,occbin50,"e");
  TH1D* hMulEta = (TH1D*) hMul->ProjectionX("hMulEta",ptbin10,ptbins,occbin0,occbin50,"e");
  TH1D* hRecEta = (TH1D*) hRec->ProjectionX("hRecEta",ptbin10,ptbins,occbin0,occbin50,"e");
  TH1D* hFakEta = (TH1D*) hFak->ProjectionX("hFakEta",ptbin10,ptbins,occbin0,occbin50,"e");
  TH1D* hSecEta = (TH1D*) hSec->ProjectionX("hSecEta",ptbin10,ptbins,occbin0,occbin50,"e");

  // projected hists: pt > 1.0 GeV/c. cBin (50-100%)
  TH1D* hSimEta2 = (TH1D*) hSim->ProjectionX("hSimEta2",ptbin10,ptbins,occbin50,occbin100,"e");
  TH1D* hEffEta2 = (TH1D*) hEff->ProjectionX("hEffEta2",ptbin10,ptbins,occbin50,occbin100,"e");
  TH1D* hMulEta2 = (TH1D*) hMul->ProjectionX("hMulEta2",ptbin10,ptbins,occbin50,occbin100,"e");
  TH1D* hRecEta2 = (TH1D*) hRec->ProjectionX("hRecEta2",ptbin10,ptbins,occbin50,occbin100,"e");
  TH1D* hFakEta2 = (TH1D*) hFak->ProjectionX("hFakEta2",ptbin10,ptbins,occbin50,occbin100,"e");
  TH1D* hSecEta2 = (TH1D*) hSec->ProjectionX("hSecEta2",ptbin10,ptbins,occbin50,occbin100,"e");

  TH1D* hDumEta = new TH1D("hDumEta",";#eta",60,-2.4,2.4); hDumEta->SetMaximum(1.0); hDumEta->SetTitle("p_{T} > 1.0");
  hDumEta->GetXaxis()->CenterTitle(); hDumEta->GetYaxis()->SetTitleOffset(1.8);
  TH1D* hDumEta2 = (TH1D*) hDumEta->Clone("hDumEta2"); hDumEta2->SetMaximum(0.1); 
  TH1D* hDumEta3 = (TH1D*) hDumEta->Clone("hDumEta3"); hDumEta3->SetMaximum(0.00049); 

  // projected hists: abs(eta) < 1.0, cBin(0-50%)
  TH1D* hSimPt  = (TH1D*) hSim->ProjectionY("hSimPt",etabin10m,etabin10p,occbin0,occbin50,"e");
  TH1D* hEffPt  = (TH1D*) hEff->ProjectionY("hEffPt",etabin10m,etabin10p,occbin0,occbin50,"e");
  TH1D* hMulPt  = (TH1D*) hMul->ProjectionY("hMulPt",etabin10m,etabin10p,occbin0,occbin50,"e");
  TH1D* hRecPt  = (TH1D*) hRec->ProjectionY("hRecPt",etabin10m,etabin10p,occbin0,occbin50,"e");
  TH1D* hFakPt  = (TH1D*) hFak->ProjectionY("hFakPt",etabin10m,etabin10p,occbin0,occbin50,"e");
  TH1D* hSecPt  = (TH1D*) hSec->ProjectionY("hSecPt",etabin10m,etabin10p,occbin0,occbin50,"e");

  // projected hists: abs(eta) < 1.0, cBin(50-100%)
  TH1D* hSimPt2  = (TH1D*) hSim->ProjectionY("hSimPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  TH1D* hEffPt2  = (TH1D*) hEff->ProjectionY("hEffPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  TH1D* hMulPt2  = (TH1D*) hMul->ProjectionY("hMulPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  TH1D* hRecPt2  = (TH1D*) hRec->ProjectionY("hRecPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  TH1D* hFakPt2  = (TH1D*) hFak->ProjectionY("hFakPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  TH1D* hSecPt2  = (TH1D*) hSec->ProjectionY("hSecPt2",etabin10m,etabin10p,occbin50,occbin100,"e");
  
  TH1D* hDumPt = new TH1D("hDumPt",";p_{T} [GeV/c]",80,0.1,300.0); hDumPt->SetMaximum(1.0);
  hDumPt->GetXaxis()->CenterTitle(); hDumPt->GetYaxis()->SetTitleOffset(1.8); hDumPt->SetTitle("|#eta| < 1");
  TH1D* hDumPt2 = (TH1D*) hDumPt->Clone("hDumPt2"); hDumPt2->SetMaximum(0.1); 
  TH1D* hDumPt3 = (TH1D*) hDumPt->Clone("hDumPt3"); hDumPt3->SetMaximum(0.00049); 

  // Efficiency
  TGraphAsymmErrors *gEffEta = new TGraphAsymmErrors(); gEffEta->SetName("gEffEta");
  gEffEta->BayesDivide(hEffEta,hSimEta);
  gEffEta->SetMarkerStyle(25);
  gEffEta->SetLineStyle(2);
  gEffEta->SetLineColor(2);
  gEffEta->SetMarkerColor(2);

  TGraphAsymmErrors *gEffPt = new TGraphAsymmErrors(); gEffPt->SetName("gEffPt");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TH1D* hSecEta2 = (TH1D*) hSec->ProjectionX("hSecEta2",ptbin20,ptbins,"e");

  TH1D* hDumEta = new TH1D("hDumEta",";#eta",60,-2.4,2.4); hDumEta->SetMaximum(1.0);
  hDumEta->GetXaxis()->CenterTitle(); hDumEta->GetYaxis()->SetTitleOffset(1.8);
  TH1D* hDumEta2 = (TH1D*) hDumEta->Clone("hDumEta2"); hDumEta2->SetMaximum(0.1); 
  TH1D* hDumEta3 = (TH1D*) hDumEta->Clone("hDumEta3"); hDumEta3->SetMaximum(0.00049); 


  // projected hists: abs(eta) < 1.0
  TH1D* hSimPt  = (TH1D*) hSim->ProjectionY("hSimPt",etabin10m,etabin10p,"e");
  TH1D* hAccPt  = (TH1D*) hAcc->ProjectionY("hAccPt",etabin10m,etabin10p,"e");
  TH1D* hEffPt  = (TH1D*) hEff->ProjectionY("hEffPt",etabin10m,etabin10p,"e");
  TH1D* hMulPt  = (TH1D*) hMul->ProjectionY("hMulPt",etabin10m,etabin10p,"e");
  TH1D* hRecPt  = (TH1D*) hRec->ProjectionY("hRecPt",etabin10m,etabin10p,"e");
  TH1D* hFakPt  = (TH1D*) hFak->ProjectionY("hFakPt",etabin10m,etabin10p,"e");
  TH1D* hSecPt  = (TH1D*) hSec->ProjectionY("hSecPt",etabin10m,etabin10p,"e");

  // projected hists: abs(eta) < 2.4
  TH1D* hSimPt2  = (TH1D*) hSim->ProjectionY("hSimPt2",etabin24m,etabin24p,"e");
  TH1D* hAccPt2  = (TH1D*) hAcc->ProjectionY("hAccPt2",etabin24m,etabin24p,"e");
  TH1D* hEffPt2  = (TH1D*) hEff->ProjectionY("hEffPt2",etabin24m,etabin24p,"e");
  TH1D* hMulPt2  = (TH1D*) hMul->ProjectionY("hMulPt2",etabin24m,etabin24p,"e");
  TH1D* hRecPt2  = (TH1D*) hRec->ProjectionY("hRecPt2",etabin24m,etabin24p,"e");
  TH1D* hFakPt2  = (TH1D*) hFak->ProjectionY("hFakPt2",etabin24m,etabin24p,"e");
  TH1D* hSecPt2  = (TH1D*) hSec->ProjectionY("hSecPt2",etabin24m,etabin24p,"e");
  
  TH1D* hDumPt = new TH1D("hDumPt",";p_{T} [GeV/c]",80,0.0,10.0); hDumPt->SetMaximum(1.0);
  hDumPt->GetXaxis()->CenterTitle(); hDumPt->GetYaxis()->SetTitleOffset(1.8);
  TH1D* hDumPt2 = (TH1D*) hDumPt->Clone("hDumPt2"); hDumPt2->SetMaximum(0.1); 
  TH1D* hDumPt3 = (TH1D*) hDumPt->Clone("hDumPt3"); hDumPt3->SetMaximum(0.00049); 

  // Acceptance
  TGraphAsymmErrors *gAccEta = new TGraphAsymmErrors(); gAccEta->SetName("gAccEta");
  gAccEta->BayesDivide(hAccEta,hSimEta);
  gAccEta->SetMarkerStyle(25);
  gAccEta->SetLineStyle(2);
  gAccEta->SetLineColor(2);
  gAccEta->SetMarkerColor(2);

  TGraphAsymmErrors *gAccPt = new TGraphAsymmErrors(); gAccPt->SetName("gAccPt");
  gAccPt->BayesDivide(hAccPt,hSimPt);
  gAccPt->SetMarkerStyle(24);
  gAccPt->SetLineColor(4);
  gAccPt->SetMarkerColor(4);

  TGraphAsymmErrors *gAccEta2 = new TGraphAsymmErrors(); gAccEta2->SetName("gAccEta2");
  gAccEta2->BayesDivide(hAccEta2,hSimEta2);
  gAccEta2->SetMarkerStyle(24);
  gAccEta2->SetLineColor(4);
  gAccEta2->SetMarkerColor(4);

  TGraphAsymmErrors *gAccPt2 = new TGraphAsymmErrors(); gAccPt2->SetName("gAccPt2");
  gAccPt2->BayesDivide(hAccPt2,hSimPt2);
  gAccPt2->SetMarkerStyle(25);
  gAccPt2->SetLineStyle(2);
  gAccPt2->SetLineColor(2);
  gAccPt2->SetMarkerColor(2);

  TLegend *legEta = new TLegend(0.35,0.3,0.65,0.5);
  legEta->SetFillColor(0); legEta->SetBorderSize(0);
  legEta->AddEntry(gAccEta,"p_{T} > 0.9 GeV/c","lp");
  legEta->AddEntry(gAccEta2,"p_{T} > 2.0 GeV/c","lp");

  TLegend *legPt = new TLegend(0.4,0.3,0.6,0.5);
  legPt->SetFillColor(0); legPt->SetBorderSize(0);
  legPt->AddEntry(gAccPt2,"|#eta| < 2.4","lp");