C++ TH1D::SetTitleOffset方法代码示例

void ptDependence()
{
   TFile *inf = new TFile("histos/ppMC.root");
   TTree *t = (TTree*) inf->Get("nt");

   const int nBins = 4;
   double ptBin[nBins+1] = {100,120,140,160,200};
//   const int nBins = 1;
//   double ptBin[nBins+1] = {100,400};
   
   TH1D *hProb = new TH1D("hProb","",nBins,ptBin);
   TH1D *hCSV = new TH1D("hCSV","",nBins,ptBin);
   TH1D *hSVTXM = new TH1D("hSVTXM","",nBins,ptBin);
   TProfile *pGen = new TProfile("pGen","",nBins,ptBin);
   
   for (int n=0; n<nBins;n++)
   {
      RooRealVar f1 = bfractionFit("discr_prob",0,3.5,ptBin[n],ptBin[n+1]);
      RooRealVar f2 = bfractionFit("discr_csvSimple",0,1,ptBin[n],ptBin[n+1]);
      RooRealVar f3 = bfractionFit("svtxm",0,6,ptBin[n],ptBin[n+1]);
      hProb->SetBinContent(n+1,f1.getVal());    
      hProb->SetBinError(n+1,f1.getError());    
      hCSV->SetBinContent(n+1,f2.getVal());    
      hCSV->SetBinError(n+1,f2.getError());    
      hSVTXM->SetBinContent(n+1,f3.getVal());    
      hSVTXM->SetBinError(n+1,f3.getError());    
   }
   
   TCanvas *c2 = new TCanvas("c2","",600,600);
   hProb->SetAxisRange(0,0.05,"Y");
   hProb->SetXTitle("Jet p_{T} (GeV/c)");
   hProb->SetYTitle("b-jet fraction");
   hProb->SetTitleOffset(1.5,"Y");
   hProb->Draw();
   hCSV->SetLineColor(2);
   hCSV->SetMarkerColor(2);
   hCSV->SetMarkerStyle(24);
   hCSV->Draw("same");
   hSVTXM->SetLineColor(4);
   hSVTXM->SetMarkerColor(4);
   hSVTXM->SetMarkerStyle(24);
//   hSVTXM->Draw("same");
   t->Draw("abs(refparton_flavorForB)==5:jtpt","","prof same");
   
   TLegend *leg = new TLegend(0.2,0.7,0.5,0.9);
   leg->SetBorderSize(0);
   leg->SetFillStyle(0);
   leg->SetFillColor(0);
   leg->AddEntry(hProb,"Jet Probability","pl");
   leg->AddEntry(hCSV,"CSV","pl");
//   leg->AddEntry(hSVTXM,"SV mass","pl");
   leg->Draw();
}

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

  data->Draw("E same");
  data->SetMarkerStyle(20);
  
  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
   
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7, 0.6, 0.8, 0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(ttgamma , "t#bar{t}#gamma", "lf");

	if(inclTop == false){
        tleg2->AddEntry(tt_sig , "t#bar{t}#gamma (MG)", "lf");
        tleg2->AddEntry(tt_lep , "t#bar{t} l #rightarrow #gamma (MG)", "lf");
        tleg2->AddEntry(tt , "t#bar{t} fake (MG)", "lf");
        }else{
        tleg2->AddEntry(tt , "t#bar{t}", "lf");
        }
	tleg2->AddEntry(T_tW, "Single Top"      , "lf");
//	tleg2->AddEntry(Tbar_tW, "anti-single-tW"      , "lf");	
	tleg2->AddEntry(DY1 , "Z+Jets", "lf");
//	tleg2->AddEntry(DY2 , "DYJetsToLL", "lf");
	tleg2->AddEntry(ZZ, "Diboson", "lf");
//	tleg2->AddEntry(WW, "WW", "lf");
//	tleg2->AddEntry(WZ , "WZ", "lf");
	tleg2->AddEntry(wjets , "W+Jets", "lf");
	tleg2->AddEntry(QCD_all, "QCD", "lf");
	
 	tleg2->Draw("same");	

	TText* textPrelim = doPrelim(0.20, 0.96, Cut, Obj);
        textPrelim->Draw();

  c1->cd();

  if(logPlot == true){
     c1->SetLogy();
  }

   TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.3);
   pad2->SetTopMargin(0);
   pad2->SetBottomMargin(0.4);
   pad2->Draw();
   pad2->cd();

   TH1D * ratio = (TH1D*)data->Clone("ratio plot");
   //ratio->Sumw2();
   ratio->SetStats(0);
   ratio->Divide(allMC);
   ratio->SetMinimum(0);
   ratio->SetMaximum(2);
   ratio->SetFillColor(kBlack);
   ratio->SetFillStyle(3354);
   ratio->SetMarkerSize(0.);

   cout << "width: " << ratio->GetBinWidth(1) << std::endl;
   
   if(MaxX <= 0.1){
   ratio->SetAxisRange(MinX, MaxX);
   }else{ 
   ratio->SetAxisRange(MinX, MaxX-ratio->GetBinWidth(1));
   }

   ratio->SetLabelSize(0.1, "X");
   ratio->SetTitleOffset(0.5, "Y");
   ratio->GetYaxis()->SetTitle("data/MC");ratio->GetYaxis()->SetTitleSize(0.1);
   ratio->GetXaxis()->SetTitle(Xtitle);ratio->GetXaxis()->SetTitleSize(0.15);
   ratio->Draw("ep");

   TLine *line;
   if(MaxX <= 0.1){
   line = new TLine(MinX,1,MaxX,1);
   }else{
   line = new TLine(MinX,1,MaxX-ratio->GetBinWidth(1),1);
   }
   
   line->Draw();

  TString plotName("plots/Control/"+ Obj + Cut  + Type + Next );
  
  if(logPlot == true){
    plotName += Variable+"ge1b_Log.pdf";   
  }else{
    plotName += Variable+"ge1b";  
  }
 
  c1->SaveAs(plotName+".pdf");
  c1->SaveAs(plotName+".png");
  
  delete c1;
  
  }
  	
}

TCanvas *drawRatioPlot(TH1D *prediction, TH1D *sr, TH1D *data, TString xTitle, TString filename, double ecaloCut){

  gStyle -> SetPadLeftMargin(0.20);

  data->SetMarkerStyle(20);
  data->SetMarkerColor(kGreen);
  data->SetLineColor(kGreen);


  
  TCanvas *c = new TCanvas("c"+filename,"c",0,0,500,500);

  float y = 0.3;

  TPad *pad1     = new TPad("pad1", "Control Plots 1", 0.01, y, 0.99, 0.99);
  TPad *padRatio = new TPad("rp1", "Ratio1", 0.01, 0.01, 0.99, y-0.01);

  pad1->SetNumber(100);
  pad1->SetTicks(0,1);
  cout<<"number pad1     = "<<pad1->GetNumber()<<endl;
  cout<<"number padRatio = "<<padRatio->GetNumber()<<endl;

  TH1D *ratio = 0;

  //ratio = (TH1D*) sr->Clone();
  //ratio->Divide(prediction);

  ratio = (TH1D*) prediction->Clone();
  ratio->Divide(data);



  for(int i=1; i<=ratio->GetNbinsX();i++){

    if(ratio->GetBinContent(i) != 0){
      cout<<"N in CR in "<<i<<". bin ="<<prediction->GetBinContent(i)<<endl;
      cout<<"N in SR in "<<i<<". bin ="<<sr->GetBinContent(i)<<endl;
      cout<<"Rel. difference in "<<i<<". bin ="<<(1./ratio->GetBinContent(i)-1.)*100<<"%"<<endl; 
    }
    else if(sr->GetBinContent(i) == 0 && prediction->GetBinContent(i) !=0)    cout<<"Scaling Factor in "<<i<<". bin <"<<prediction->GetBinContent(i)/1.15<<" +/- "<<ratio->GetBinError(i)<<endl;
    else if(sr->GetBinContent(i) != 0 && prediction->GetBinContent(i) ==0)    cout<<"Scaling Factor in "<<i<<". bin <"<<(sr->GetEntries()/prediction->GetEntries())/sr->GetBinContent(i)<<" +/- "<<ratio->GetBinError(i)<<endl;

  }

  ratio->GetYaxis()->SetTitle("#frac{CR (MC)}{CR (data)}");
  ratio->SetTitle("");
  ratio->SetLabelSize(0.1,"X");
  ratio->SetLabelSize(0.1,"Y");
  ratio->SetTitleOffset(0.5,"Y");
  ratio->SetTitleSize(0.15,"Y");


  padRatio->cd();
  //ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->Draw("e");

  // Draw line at one!
  float xmin = ratio->GetXaxis()->GetXmin();
  float xmax = ratio->GetXaxis()->GetXmax();
  TLine *line = new TLine(xmin,1,xmax,1);
  line->SetLineWidth(2);
  line->Draw("same");
  padRatio->Modified();

  TLegend *leg = new TLegend(0.5,0.7,0.9,0.9);
  leg->AddEntry(sr,"SR (MC)","l"); 
  leg->AddEntry(prediction,"lepton CR (MC)","pel"); 
 
  pad1->cd();
  pad1->SetLogy();
  //  pad1->SetLogx();

  sr->SetLineColor(kRed);
  sr->SetMarkerColor(kRed);
  sr->SetMarkerStyle(20);
  sr->SetTitle("");
  prediction->SetMarkerStyle(20);
  prediction->SetTitle("");
  prediction->GetXaxis()->SetTitle(xTitle);
  sr->GetXaxis()->SetTitle(xTitle);

  
  

  prediction->SetTitleSize(0.07,"X");
  prediction->GetXaxis()->SetTitleOffset(0.7);
  sr->SetTitleSize(0.07,"X");
  sr->GetXaxis()->SetTitleOffset(0.7);

  double maximum = 0;
  double minimum = 1000000;
  if(sr->GetMinimum()!=0 && sr->GetMinimum()<minimum){
    minimum=sr->GetMinimum()*0.5;
  }
  if(prediction->GetMinimum()!=0 && prediction->GetMinimum()<minimum){
    minimum=prediction->GetMinimum()*0.5;
  }
  if(data->GetMinimum()!=0 && data->GetMinimum()<minimum){
    minimum=data->GetMinimum()*0.5;
  }
//.........这里部分代码省略.........

void plot(TString var, TString varlatex, TString varname, Int_t nbins, Double_t vmin, Double_t vmax)
{
  cout<<"---- Processing - "<<var<<endl;
  cout<<"  -- Fill histograms"<<endl;
  TFile* ifBkg = new TFile(infnameBkg[isChannel]);
  TTree* ntBkg = (TTree*)ifBkg->Get(texNtuple[isChannel]);
  ntBkg->AddFriend("ntHlt");
  TFile* ifSgl = new TFile(infnameSgl[isChannel]);
  TTree* ntSgl = (TTree*)ifSgl->Get(texNtuple[isChannel]);
  ntSgl->AddFriend("ntHlt");
  ntSgl->AddFriend("ntHi");

  TH1D* hBkg = new TH1D(Form("hBkg_%s",varname.Data()),"",nbins,vmin,vmax);
  TH1D* hSgl = new TH1D(Form("hSgl_%s",varname.Data()),"",nbins,vmin,vmax);

  ntBkg->Project(Form("hBkg_%s",varname.Data()),var,Form("%s&&%s",selTriggerBkg[isChannel].Data(),selBkg[isChannel].Data()));
  ntSgl->Project(Form("hSgl_%s",varname.Data()),var,TCut(weight[isChannel])*Form("%s&&%s",selTriggerSgl[isChannel].Data(),selSgl[isChannel].Data()));

  cout<<"  -- Calculate normalization"<<endl;
  Double_t normBkg=0,normSgl=0;
  //normBkg = hBkg->GetEntries();
  //normSgl = hSgl->GetEntries();
  normBkg = hBkg->Integral(vmin,vmax);
  normSgl = hSgl->Integral(vmin,vmax);
  cout<<"     normBkg: "<<normBkg<<" ;   normSgl: "<<normSgl<<endl;

  cout<<"  -- Normalize histograms"<<endl;
  hBkg->Scale(1./normBkg);
  hSgl->Scale(1./normSgl); 

  cout<<"  -- Plot"<<endl;
  hBkg->SetXTitle(varlatex);
  hBkg->SetYTitle("#Probability");
  hBkg->SetTitleOffset(1.5,"Y");
  Double_t hisMax = (hBkg->GetMaximum()>hSgl->GetMaximum())?hBkg->GetMaximum():hSgl->GetMaximum();
  hBkg->SetMaximum(hisMax*1.2);

  hBkg->SetLineColor(kBlue+1);
  hBkg->SetFillStyle(1001);
  hBkg->SetFillColor(kBlue-9);
  hBkg->SetLineWidth(3);
  hBkg->SetStats(0);

  TH1D* hSglplot = new TH1D(Form("hSglplot_%s",varname.Data()),"",nbins,vmin,vmax);
  for(int ib=0;ib<nbins;ib++) hSglplot->SetBinContent(ib+1,hSgl->GetBinContent(ib+1));
  hSglplot->SetLineColor(kRed);
  hSglplot->SetFillStyle(3004);
  hSglplot->SetFillColor(kRed);
  hSglplot->SetLineWidth(3);
  hSglplot->SetStats(0);

  TCanvas* c = new TCanvas(Form("c_%s",varname.Data()),"",600,600);
  hBkg->Draw();
  hSglplot->Draw("same");

  cout<<"  -- Plot legends"<<endl;
  TLatex* tex = new TLatex(0.18,0.935,Form("5.02TeV %s",texPP[isChannel].Data()));
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.055);
  tex->Draw();
  TLatex* texp;
  texp = new TLatex(0.68,0.935,texDecay[isChannel]);
  texp->SetNDC();
  texp->SetTextFont(42);
  texp->SetTextSize(0.055);
  texp->Draw();
  TLegend* leg = new TLegend(0.56,0.70,0.86,0.86);
  leg->AddEntry(hBkg,"Background","f");
  leg->AddEntry(hSglplot,"Signal","f");
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->Draw("same");

  cout<<"  -- Save plots"<<endl;
  c->SaveAs(Form("plots/%s_%s/c_%s.pdf",texPP[isChannel].Data(),texNtuple[isChannel].Data(),varname.Data()));
  cout<<endl;
}

void dimuonLoop(float ptmin, float ptmax, float ymin, float ymax, int i)
{
  TFile *infData = new TFile(inputdata);
  TFile *infMC = new TFile(inputmc);
  TTree *ntData = (TTree*) infData->Get("ntKp");
  TTree *ntMC = (TTree*) infMC->Get("ntKp");

  TCanvas *c = new TCanvas(Form("c%i",i),"",600,600);
  //if (logy) c->SetLogy();
  
  TH1D* hData = new TH1D("hData","",50,2.85,3.35);
  TH1D* hMC = new TH1D("hMC","",50,2.85,3.35);

  ntData->Project("hData","mumumass",Form("%s&&(pt>%f&&pt<%f)&&(y>%f&&y<%f)",cut.Data(),ptmin,ptmax,ymin,ymax));
  ntMC->Project("hMC","mumumass",Form("%s&&(pt>%f&&pt<%f)&&(y>%f&&y<%f)",cut.Data(),ptmin,ptmax,ymin,ymax));

  hData->Sumw2();

  double normData=0,normMC=0;
  normData = hData->GetEntries();
  normMC = hMC->GetEntries();
  cout<<normData<<" "<<normMC<<endl;
  
  hData->Scale(1./normData);
  hMC->Scale(1./normMC); 

  hData->Sumw2();

  hMC->SetXTitle("#mu#mu mass");
  hMC->SetYTitle("#Probability");
  hMC->SetTitleOffset(1.5,"Y");
  if(hData->GetMaximum()>hMC->GetMaximum()) hMC->SetMaximum(hData->GetMaximum()*1.1);
  else hMC->SetMaximum(hMC->GetMaximum()*1.1);

  //hData->SetLineColor(kBlue+1);
  //hData->SetFillStyle(1001);
  //hData->SetFillColor(kBlue-9);
  //hData->SetLineWidth(3);
  hData->SetMarkerStyle(8);
  hData->SetStats(0);

  hMC->SetLineColor(kRed);
  hMC->SetFillStyle(3004);
  hMC->SetFillColor(kRed);
  hMC->SetLineWidth(3);
  hMC->SetStats(0);

  hMC->Draw();
  hData->Draw("same lep");
  
  TLegend *leg = new TLegend(0.11,0.75,0.50,0.9);
  leg->AddEntry(hData,"data","lep");
  leg->AddEntry(hMC,"MC","f");
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->Draw("same");
  TLegend *leg1 = new TLegend(0.60,0.7,0.90,0.9);
  leg1->AddEntry((TObject*)0,Form("%s",particle.Data()),"");
  leg1->AddEntry((TObject*)0,Form("%.0f<p_{T}<%.0f GeV",ptmin,ptmax),"");
  leg1->AddEntry((TObject*)0,Form("%.1f<y_{CM}<%.1f",ymin,ymax),"");
  leg1->SetBorderSize(0);
  leg1->SetFillStyle(0);
  leg1->Draw("same");

  c->SaveAs(Form("MuonResult/plot_Bplus_chi2_trkPt/dimuon_yDM_%i.pdf",i));
  //c->SaveAs(Form("MuonResult/plot_Bplus_chi2_trkPt/dimuon_DM_%i.pdf",i));
  //c->SaveAs("MuonResult/plot_Bplus_chi2_trkPt/dimuon_DM_incl_1.pdf");
}

//################################################################################################################################
//################################################################################################################################
//################################################################################################################################
TCanvas *drawRatioPlot(TH1D *prediction, TH1D *sr, TString xTitle, TString filename){

  
  TCanvas *c = new TCanvas("c"+filename,"c",0,0,500,500);
  c->cd();
  float y = 0.3;

  TPad *pad1     = new TPad("pad1", "Control Plots 1", 0.01, y, 0.99, 0.99);
  TPad *padRatio = new TPad("rp1", "Ratio1", 0.01, 0.01, 0.99, y-0.01);

  TH1D *ratio = 0;

  ratio = (TH1D*) sr->Clone();
  ratio->Divide(prediction);
  ratio->GetYaxis()->SetTitle((TString) sr->GetName() + "/" + (TString) prediction->GetName());
  ratio->SetTitle("");
  ratio->SetLabelSize(0.1,"X");
  ratio->SetLabelSize(0.1,"Y");
  ratio->SetTitleOffset(0.5,"Y");
  ratio->SetTitleSize(0.15,"Y");


  padRatio->cd();
  ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->Draw();

  // Draw line at one!
  float xmin = ratio->GetXaxis()->GetXmin();
  float xmax = ratio->GetXaxis()->GetXmax();
  TLine *line = new TLine(xmin,1,xmax,1);
  line->SetLineWidth(2);
  line->Draw("same");
  padRatio->Modified();

  TLegend *leg = new TLegend(0.6,0.8,0.9,0.9);
  leg->AddEntry(sr,"bkg","l"); 
  leg->AddEntry(prediction,prediction->GetName(),"pel"); 
  pad1->cd();
  //pad1->SetLogy();

  sr->SetLineColor(kRed);
  sr->SetMarkerColor(kRed);
  sr->SetTitle("");
  prediction->SetTitle("");
  prediction->GetXaxis()->SetTitle(xTitle);
  sr->GetXaxis()->SetTitle(xTitle);

  prediction->SetTitleSize(0.07,"X");
  prediction->GetXaxis()->SetTitleOffset(0.7);
  sr->SetTitleSize(0.07,"X");
  sr->GetXaxis()->SetTitleOffset(0.7);

  prediction->Draw("e");
  sr->Draw("e same");

  if(sr->GetMaximum()>prediction->GetMaximum()){

    sr->Draw("e");
    prediction->Draw("e same");

    }

  leg->Draw("same");




  // Draw both pads to canvas
  c->cd();
  pad1->Draw();
  padRatio->SetGridy();
  padRatio->Draw();

  c->SaveAs(filename);

  return c;

}

//.........这里部分代码省略.........
    h_data->SetMarkerSize(1);
    h_data->SetMarkerStyle(24);
    h_data->SetTitle("");
    h_data->Draw("e");
    if(xmin>-9999 && xmax>-9999) {
        h_data->GetXaxis()->SetRangeUser(xmin,xmax);
    }
    cout << "Data integral = " << h_data->Integral() << endl;
    hs->Draw("histsame");
    h_data->Draw("esame");
    leg->Draw("same");


    std::string dirName = "compareDataMC_" + leptonName;
    gSystem->mkdir(dirName.data());

    std::string filename;
    std::string psname ;
    psname = dirName+ "/overlay_" + histoName;
    filename = psname + ".eps";
    c1->Print(filename.data());
    filename = psname + ".gif";
    c1->Print(filename.data());
    filename = psname + ".pdf";
    c1->Print(filename.data());


    // study the ratios

    TCanvas* c2 = new TCanvas("c2","",700,0,700,1000);
    c2->Divide(1,2,0.01,0);
    c2->cd(1);
    gPad->SetTopMargin(0.01);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.04);

    h_data->Draw("e");
    if(xmin>-9999 && xmax>-9999) {
        h_data->GetXaxis()->SetRangeUser(xmin,xmax);
    }
    cout << h_data->GetName() << " integral = " << h_data->Integral() << endl;
    hs->Draw("histsame");
    h_data->Draw("esame");
    leg->Draw("same");


    c2->cd(2);
    gStyle->SetStatW       (0.3);
    gStyle->SetStatH       (0.3);
    gStyle->SetStatX       (0.879447);
    gStyle->SetStatY       (0.939033);
    gStyle->SetStatFontSize(0.05);
    gStyle->SetStatBorderSize(0);
    gPad->SetRightMargin(0.04);
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0.2);
    gPad->SetTickx();
    gStyle->SetOptFit(1);

    TH1D* hratio = (TH1D*)h_data->Clone("hratio");
    hratio->Reset();
    hratio->Divide(h_data,h_all,1.0,1.0);
    hratio->SetTitle("");
    hratio->SetMaximum(1.5);
    hratio->SetMinimum(0.5);
    hratio->SetTitleOffset(1.2,"Y");
    hratio->GetXaxis()->SetTitle(h_data->GetXaxis()->GetTitle());
    hratio->Draw("e1");

    cout << "( " << h_data->GetBinContent(maxBin) << "+-" << h_data->GetBinError(maxBin) << " )/("
         << h_all->GetBinContent(maxBin) << "+-" << h_all->GetBinError(maxBin) <<  ")= "
         << hratio->GetBinContent(maxBin) << "+-" << hratio->GetBinError(maxBin) << endl;

    psname = dirName+ "/ratio_" + histoName;
    filename = psname + ".eps";
    c2->Print(filename.data());
    filename = psname + ".gif";
    c2->Print(filename.data());
    filename = psname + ".pdf";
    c2->Print(filename.data());


    std::string remword  =".txt";
    size_t pos  = inputFile.find(remword);

    if(pos!= std::string::npos)
        inputFile.replace(pos,remword.length(),"");

    std::string command = "recreate";

    if(update)command ="update";

    TFile* outFile = new TFile(Form("combined_%s.root",inputFile.data()),command.data());
    h_all->Write();
    hs->Write();
    h_data->Write();
    outFile->Close();


}

TF1 *fit(TTree *nt,TTree *ntMC,double ptmin,double ptmax, bool ispPb, int count){   
   //cout<<cut.Data()<<endl;
   //static int count=0;
   //count++;
   TCanvas *c= new TCanvas(Form("c%d",count),"",600,600);
   TH1D *h = new TH1D(Form("h%d",count),"",50,5,6);
   TH1D *hMC = new TH1D(Form("hMC%d",count),"",50,5,6);

   TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
   TF1 *f = new TF1(Form("f%d",count),"[0]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*("+iNP+")");
   nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   clean0(h);

   TH1D *hraw = new TH1D(Form("hraw%d",count),"",50,5,6);
   clean0(hraw);
   hraw = (TH1D*)h->Clone(Form("hraw%d",count));


   h->Draw();
   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.05);
   f->SetParLimits(8,0.01,0.05);
   f->SetParLimits(7,0,1);
   f->SetParLimits(5,0,1000);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(8,setparam3);
   f->FixParameter(1,fixparam1);
   h->GetEntries();

   hMC->Fit(Form("f%d",count),"q","",5,6);
   hMC->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L m","",5,6);

   f->FixParameter(1,f->GetParameter(1));
   f->FixParameter(2,f->GetParameter(2));
   f->FixParameter(7,f->GetParameter(7));
   f->FixParameter(8,f->GetParameter(8));
   
   h->Fit(Form("f%d",count),"q","",5,6);
   h->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L m","",5,6);
   h->SetMarkerSize(0.8);
   h->SetMarkerStyle(20);
   cout <<h->GetEntries()<<endl;

   // function for background shape plotting. take the fit result from f
   TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
   background->SetParameter(0,f->GetParameter(3));
   background->SetParameter(1,f->GetParameter(4));
   background->SetParameter(2,f->GetParameter(5));
   background->SetParameter(3,f->GetParameter(6));
   background->SetLineColor(4);
   background->SetRange(5,6);
   background->SetLineStyle(2);
   
   // function for signal shape plotting. take the fit result from f
   TF1 *Bkpi = new TF1(Form("fBkpi%d",count),"[0]*("+iNP+")");
   Bkpi->SetParameter(0,f->GetParameter(5));
   Bkpi->SetLineColor(kGreen+1);
   Bkpi->SetFillColor(kGreen+1);
//   Bkpi->SetRange(5.00,5.28);
   Bkpi->SetRange(5.00,6.00);
   Bkpi->SetLineStyle(1);
   Bkpi->SetFillStyle(3004);

   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1(Form("fmass%d",count),"[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(8));
   mass->SetParError(0,f->GetParError(0));
   mass->SetParError(1,f->GetParError(1));
   mass->SetParError(2,f->GetParError(2));
   mass->SetParError(7,f->GetParError(7));
   mass->SetParError(8,f->GetParError(8));
   mass->SetLineColor(2);
   mass->SetLineStyle(2);

//   cout <<mass->Integral(0,1.2)<<" "<<mass->IntegralError(0,1.2)<<endl;
   h->SetMarkerStyle(24);
   h->SetStats(0);
   h->Draw("e");
   h->SetXTitle("M_{B} (GeV/c^{2})");
   h->SetYTitle("Entries / (20 MeV/c^{2})");
   h->GetXaxis()->CenterTitle();
   h->GetYaxis()->CenterTitle();
   h->SetTitleOffset(1.5,"Y");
   h->SetAxisRange(0,h->GetMaximum()*1.2,"Y");
   Bkpi->Draw("same");
   background->Draw("same");   
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        hscale->SetBinContent(i,ndata/nmc);
        double err = 0;
        err=
            (ndata/nmc)*sqrt(pow(nmcerr/nmc,2)+pow(ndataerr/ndata,2));
        hscale->SetBinError(i,err);

    }


    hsherpa->GetXaxis()->SetRangeUser(xmin,xmax);
    hmadgraph->GetXaxis()->SetRangeUser(xmin,xmax);
    hscale->GetXaxis()->SetRangeUser(xmin,xmax);


    TCanvas* c1 = new TCanvas("c1","",700,1000);
    c1->Divide(1,2,0.01,0);
    c1->cd(1);
    if(logScale)
        gPad->SetLogy(1);
    gPad->SetTopMargin(0.01);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.04);


    float max_data  = hsherpa->GetBinError(hsherpa->GetMaximumBin()) + hsherpa->GetMaximum();
    float max_mc    = hmadgraph->GetBinError(hmadgraph->GetMaximumBin()) + hmadgraph->GetMaximum();

    if(max_data > max_mc)
    {
        hsherpa->Draw("e");
        hmadgraph->Draw("hesame");
    }
    else
    {   hmadgraph->Draw("he");
        hsherpa->Draw("esame");
    }


    float x1NDC = 0.691;
    float y1NDC = 0.757;
    float x2NDC = 0.894;
    float y2NDC = 0.973;

    TLegend* leg = new TLegend(x1NDC,y1NDC,x2NDC,y2NDC);

    leg->SetFillColor(0);
    leg->SetFillStyle(0);
    leg->SetTextSize(0.04);
    leg->SetBorderSize(0);
    leg->AddEntry(hsherpa, "SHERPA");
    leg->AddEntry(hmadgraph, "MADGRAPH");
    leg->Draw("same");



    c1->cd(2);
    gStyle->SetStatW       (0.3);
    gStyle->SetStatH       (0.3);
    gStyle->SetStatX       (0.879447);
    gStyle->SetStatY       (0.939033);
    gStyle->SetStatFontSize(0.05);
    gStyle->SetStatBorderSize(0);
    gPad->SetRightMargin(0.04);
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0.2);
    gPad->SetTickx();
    gStyle->SetOptFit(1);
    hscale->SetTitle("");
    hscale->SetMaximum(3.0);
    hscale->SetMinimum(-0.5);
    hscale->SetTitleOffset(1.2,"Y");
    hscale->Draw("e1");
    TF1* fline = new TF1("fline","pol1");
    TLine* l2 = new TLine(xmin,1.,xmax,1.);
    l2->SetLineColor(4);
    l2->SetLineStyle(3);
    fline->SetLineWidth(3);
    fline->SetLineColor(6);
    fline->SetNpx(2500);
    hscale->Fit("fline","","");
    l2->Draw("same");


    string dirName = "compareSherpaMadgraph";
    gSystem->mkdir(dirName.data());

    std::string filename;
    std::string psname = dirName + "/" + var;
    if(output !="test")
        psname = dirName+ "/" + output;
    else
        psname = dirName+ "/" + var;
    filename = psname + ".eps";
    c1->Print(filename.data());
    filename = psname + ".gif";
    c1->Print(filename.data());
    filename = psname + ".pdf";
    c1->Print(filename.data());
    //   c1->Close();
}

//.........这里部分代码省略.........
  cout << "Relative uncertainty is " << endl;
  for(int i=1;i<=hscale->GetNbinsX();i++)
    if(hscale->GetBinContent(i)>1e-6)
    cout << "Bin " << i << " ( " << hscale->GetBinLowEdge(i) << "~" 
	 << hscale->GetBinLowEdge(i+1) << " ): " 
	 << hscale->GetBinError(i)/hscale->GetBinContent(i) << endl;


  h1->GetXaxis()->SetRangeUser(xmin,xmax);
  h2->GetXaxis()->SetRangeUser(xmin,xmax);
  hscale->GetXaxis()->SetRangeUser(xmin,xmax);


  TCanvas* c1 = new TCanvas("c1","",700,1000);  
  c1->Divide(1,2,0.01,0);
  c1->cd(1);
  if(logScale)
    gPad->SetLogy(1);
  gPad->SetTopMargin(0.01);
  gPad->SetBottomMargin(0);
  gPad->SetRightMargin(0.04);

  
  float max_data  = h1->GetBinError(h1->GetMaximumBin()) + h1->GetMaximum();
  float max_mc    = h2->GetBinError(h2->GetMaximumBin()) + h2->GetMaximum();

  if(max_data > max_mc)
    {
      h1->Draw("e");
      h2->Draw("hesame");
    }
  else
    { h2->Draw("he");
      h1->Draw("esame");
    }


  float x1NDC = 0.798232;
  float y1NDC = 0.811949;
  float x2NDC = 0.877295;
  float y2NDC = 0.953398;

  TLegend* leg = new TLegend(x1NDC,y1NDC,x2NDC,y2NDC);
  
//   leg->SetHeader(header.data());
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->SetTextSize(0.06);
  leg->SetBorderSize(0);
  leg->AddEntry(h1, mcName1.data());
  leg->AddEntry(h2, mcName2.data());
  leg->Draw("same");



  c1->cd(2);
  gStyle->SetStatW       (0.3);
  gStyle->SetStatH       (0.3);
  gStyle->SetStatX       (0.879447);
  gStyle->SetStatY       (0.939033);
  gStyle->SetStatFontSize(0.05);
  gStyle->SetStatBorderSize(0);
  gPad->SetRightMargin(0.04);
  gPad->SetTopMargin(0);
  gPad->SetBottomMargin(0.2);
  gPad->SetTickx();
  gStyle->SetOptFit(1);
  hscale->SetTitle("");
  hscale->SetMaximum(1.2);
  hscale->SetMinimum(0.8);
  hscale->SetTitleOffset(1.2,"Y");
  hscale->Draw("e1");
  TF1* fline = new TF1("fline","pol1");
  TLine* l2 = new TLine(xmin,1.,xmax,1.);
  l2->SetLineColor(4);
  l2->SetLineStyle(3);
  fline->SetLineWidth(3);
  fline->SetLineColor(6);
  fline->SetNpx(2500);
//   hscale->Fit("fline","","");
  l2->Draw("same");


  string dirName = "unfolding";
  gSystem->mkdir(dirName.data());

  std::string filename;
  std::string psname = dirName + "/" + var1;
  if(output !="test")
    psname = dirName+ "/" + output;
  else
    psname = dirName+ "/" + var1;
  filename = psname + ".eps";
  c1->Print(filename.data());
  filename = psname + ".gif";
  c1->Print(filename.data());
  filename = psname + ".pdf";
  c1->Print(filename.data());
  //   c1->Close();
}

//.........这里部分代码省略.........
    // Print
    cout << Form("%.1f < P_{T} < %.1f GeV: ",hPt->GetBinLowEdge(iBeg),hPt->GetBinLowEdge(iEnd+1))
      << " SigSubBkg Integral - Nr: " << hNr->Integral() << " Aw: " << hAw->Integral() << endl
      << " Gen - Nr: " << hNrGen->Integral() << " Aw: " << hAwGen->Integral() << endl;
    // Styles
    hNr->SetMarkerColor(kRed);
    hNr->SetLineColor(kRed);
    hAw->SetMarkerColor(kBlue);
    hAw->SetLineColor(kBlue);
    hAwGen->SetLineStyle(2);
    // Axis Range
    if (sysMode==0) {
      hNr->SetYTitle("Background Subtracted Signal (GeV/c)");
      hNr->SetAxisRange(-5,60,"Y");
    }
    if (sysMode==2) {
      hNr->SetAxisRange(-2,6,"Y");
    }
    // Axis Label
    if (compMode==0) {
      if (sysMode==1) {
	hNr->SetYTitle("(Reco Trk)-(Sig GenP) (GeV/c)");
	//hNr->SetAxisRange(-20,20,"Y");
      }
      if (sysMode==2) {
	hNr->SetYTitle("(Reco Trk)/(Sig GenP) (GeV/c)");
      }
    }
    if (compMode==3) {
      if (sysMode==1) {
	hNr->SetYTitle("All GenP - Sig GenP (GeV/c)");
	//hNr->SetAxisRange(-3,3,"Y");
      }
      if (sysMode==2) {
	hNr->SetYTitle("(All GenP)/(Sig GenP) (GeV/c)");
      }
    }
    hNr->SetXTitle("#Delta R");
    hNr->SetAxisRange(0,0.79999,"X");
    hNr->SetTitleOffset(1.5,"X");
    hNr->GetXaxis()->CenterTitle();
    hNr->GetYaxis()->CenterTitle();
    c6->cd(i+1);
    // Fit
    if (sysMode>0) {
      TF1 * f0 = new TF1("f0","pol0");
      f0->SetLineStyle(2);
      f0->SetLineWidth(1);
      hNr->Fit("f0");
      hAw->Fit("f0");
    }
    // Draw to Inspect
    hNr->Draw();
    hAw->Draw("same");
    if (sysMode==0) {
      hNrGen->Draw("hist same");
      hAwGen->Draw("hist same");
    }
    if (sysMode==0||sysMode==1) {
      TLine *l = new TLine(0,0,0.8,0);
      l->Draw();
    }
    if (sysMode==2) {
      TLine *l = new TLine(0,1,0.8,1);
      l->Draw();
    }
    
    TLegend *leg = new TLegend(0.35,0.7,0.85,0.94);
    leg->SetBorderSize(0);
    leg->SetFillStyle(0);
    leg->AddEntry(hNr,Form("%.1f < P_{T} < %.1f GeV",hPt->GetBinLowEdge(iBeg),hPt->GetBinLowEdge(iEnd+1)),"");
    if (sysMode==0) {
      if (compMode==0) {
	leg->AddEntry(hNr,"Leading (RecoTrk)","pl");
	leg->AddEntry(hNrGen,"Leading (Sig. GenP)","l");
	leg->AddEntry(hAw,"SubLeading (RecoTrk)","pl");
	leg->AddEntry(hAwGen,"SubLeading (Sig. GenP)","l");
      }
      if (compMode==1) {
	leg->AddEntry(hNr,"Leading (RecoTrk)","pl");
	leg->AddEntry(hNrGen,"Leading (All GenP)","l");
	leg->AddEntry(hAw,"SubLeading (RecoTrk)","pl");
	leg->AddEntry(hAwGen,"SubLeading (All GenP)","l");
      }
    }
    if (sysMode>0) {
      leg->AddEntry(hNr,"Leading","pl");
      leg->AddEntry(hAw,"SubLeading","pl");
    }
    leg->SetTextSize(0.05);
    leg->Draw();
  }

  // ===============================================
  // Save
  // ===============================================
  c6->Print(Form("%s/%s.gif",outdir.Data(),tag.Data()));
  c6->Print(Form("%s/%s.eps",outdir.Data(),tag.Data()));
  c6->Print(Form("%s/%s.C",outdir.Data(),tag.Data()));
}

TCanvas *drawRatioPlot(TH1D *prediction, TH1D *sr, TString xTitle, TString filename){

  
  TCanvas *c = new TCanvas("c"+filename,"c",0,0,500,500);

  float y = 0.3;

  TPad *pad1     = new TPad("pad1", "Control Plots 1", 0.01, y, 0.99, 0.99);
  TPad *padRatio = new TPad("rp1", "Ratio1", 0.01, 0.01, 0.99, y-0.01);

  TH1D *ratio = 0;

  //ratio = (TH1D*) sr->Clone();
  //ratio->Divide(prediction);

  ratio = (TH1D*) prediction->Clone();
  ratio->Divide(sr);


  for(int i=1; i<=ratio->GetNbinsX();i++){

    cout<<"Scaling Factor in "<<i<<". bin ="<<ratio->GetBinContent(i)<<" +/- "<<ratio->GetBinError(i)<<endl; 

  }

  ratio->GetYaxis()->SetTitle("Pred/Bkg");
  ratio->SetTitle("");
  ratio->SetLabelSize(0.1,"X");
  ratio->SetLabelSize(0.1,"Y");
  ratio->SetTitleOffset(0.5,"Y");
  ratio->SetTitleSize(0.15,"Y");


  padRatio->cd();
  //ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->Draw();

  // Draw line at one!
  float xmin = ratio->GetXaxis()->GetXmin();
  float xmax = ratio->GetXaxis()->GetXmax();
  TLine *line = new TLine(xmin,1,xmax,1);
  line->SetLineWidth(2);
  line->Draw("same");
  padRatio->Modified();

  TLegend *leg = new TLegend(0.6,0.8,0.9,0.9);
  leg->AddEntry(sr,"bkg","l"); 
  leg->AddEntry(prediction,"bkg prediction","pel"); 
  pad1->cd();
  pad1->SetLogy();

  sr->SetLineColor(kRed);
  sr->SetMarkerColor(kRed);
  sr->SetMarkerStyle(20);
  sr->SetTitle("");
  prediction->SetMarkerStyle(20);
  prediction->SetTitle("");
  prediction->GetXaxis()->SetTitle(xTitle);
  sr->GetXaxis()->SetTitle(xTitle);

  prediction->SetTitleSize(0.07,"X");
  prediction->GetXaxis()->SetTitleOffset(0.7);
  sr->SetTitleSize(0.07,"X");
  sr->GetXaxis()->SetTitleOffset(0.7);

  prediction->Draw("e");
  sr->Draw("e same");

  double maximum =prediction->GetMaximum()*2.5;
  double minimum = 0.00001;
  if(sr->GetMinimum()!=0) minimum=sr->GetMinimum()*0.5;
  if(prediction->GetMinimum()!=0) minimum=prediction->GetMinimum()*0.5;
  if(sr->GetMaximum()>prediction->GetMaximum()){

    maximum=sr->GetMaximum()*2.5;
  }

  prediction->GetYaxis()->SetRangeUser(minimum,maximum);

  leg->Draw("same");




  // Draw both pads to canvas
  c->cd();
  pad1->Draw();
  padRatio->SetGridy();
  padRatio->Draw();

  c->SaveAs(filename);

  return c;

}

//.........这里部分代码省略.........
	TFile *rebinned = new TFile("~/Summer2016/root_forgraphs/figure2_0_notcorrected.root","RECREATE");
	gr_pos->Write();
	gr_neg->Write();
	gr_diff->Write();
	

/*
	gr_pos->RemovePoint(0);
	gr_pos->RemovePoint(5);

	gr_neg->RemovePoint(0);
	gr_neg->RemovePoint(5);

	gr_diff->RemovePoint(0);
	gr_diff->RemovePoint(5);
*/




	gStyle->SetLegendFont(42);
	TH1D* base = new TH1D("base","base",1,-0.20,0.20);
	//pPb
	//base->GetYaxis()->SetRangeUser(0.065, 0.075);
	base->GetYaxis()->SetRangeUser(0.06, 0.12);

	//PbPb
	//base->GetYaxis()->SetRangeUser(0.093, 0.103);
	base->GetXaxis()->SetTitle("Observed A_{ch}");
	base->GetYaxis()->SetTitle("v_{2}{2}");
	base->GetXaxis()->CenterTitle();
	base->GetYaxis()->CenterTitle();
	base->SetTitleSize  (0.040,"X");
	base->SetTitleOffset(1.4,"X");
	base->SetTitleFont  (42,"X");
	base->SetLabelOffset(0.006,"X");
	base->SetLabelSize  (0.040,"X");
	base->SetLabelFont  (42   ,"X");

	base->SetTitleSize  (0.040,"Y");
	base->SetTitleOffset(2.2,"Y");
	base->SetTitleFont  (42,"Y");
	base->SetLabelOffset(0.006,"Y");
	base->SetLabelSize  (0.040,"Y");
	base->SetLabelFont  (42   ,"Y");
	base->SetLineWidth(0);

	TH1D* base2 = new TH1D("base2","base2",1,-0.2,0.2);
	base2->GetYaxis()->SetRangeUser(-0.03, 0.03);
	base2->GetXaxis()->SetTitle("Observed A_{ch}");
	base2->GetYaxis()->SetTitle(" (v^{#minus}_{2} #minus v^{#plus}_{2})/(v^{#minus}_{2} #plus v^{#plus}_{2}) ");
	base2->GetXaxis()->CenterTitle();
	base2->GetYaxis()->CenterTitle();
	base2->SetTitleSize  (0.040,"X");
	base2->SetTitleOffset(1.4,"X");
	base2->SetTitleFont  (42,"X");
	base2->SetLabelOffset(0.006,"X");
	base2->SetLabelSize  (0.040,"X");
	base2->SetLabelFont  (42   ,"X");
	base2->SetTitleSize  (0.040,"Y");
	base2->SetTitleOffset(2.0,"Y");
	base2->SetTitleFont  (42,"Y");
	base2->SetLabelOffset(0.006,"Y");
	base2->SetLabelSize  (0.040,"Y");
	base2->SetLabelFont  (42   ,"Y");
	base2->SetLineWidth(0);

TF1 *fitSignal(TTree *ntMC,double ptmin,double ptmax){   
   //cout<<cut.Data()<<endl;
   TCanvas *cSignal= new TCanvas("cSignal","",600,600);
   TH1D *hMC = new TH1D("hMC","",50,5,6);

   //   TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
   TF1 *f = new TF1("f1","[0]*([5]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[5])*Gaus(x,[1],[6])/(sqrt(2*3.14159)*[6]))+[3]+[4]*x");
   ntMC->Project("hMC","mass",Form("%s&&pt>%f&&pt<%f",selmc.Data(),ptmin,ptmax));

   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.05);
   f->SetParLimits(6,0.01,0.05);
   f->SetParLimits(5,0,1);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(6,setparam3);
   f->FixParameter(1,fixparam1);

   hMC->Fit("f1","q","",5,6);
   hMC->Fit("f1","q","",5,6);
   f->ReleaseParameter(1);
   hMC->Fit("f1","L q","",5,6);
   hMC->Fit("f1","L q","",5,6);
   hMC->Fit("f1","L q","",5,6);
   hMC->Fit("f1","L m","",5,6);

   // function for background shape plotting. take the fit result from f
   TF1 *background = new TF1("backgroundSignal","[0]+[1]*x");
   background->SetParameter(0,f->GetParameter(3));
   background->SetParameter(1,f->GetParameter(4));
   background->SetLineColor(4);
   background->SetRange(5,6);
   background->SetLineStyle(2);
   
   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1("fmassSignal","[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(5),f->GetParameter(6));
   cout<<"[0]Signal: "<<f->GetParameter(0)<<endl;
   mass->SetLineColor(2);
   mass->SetLineStyle(2);

   hMC->SetMarkerStyle(24);
   hMC->SetStats(0);
   hMC->Draw("e");
   hMC->SetXTitle("M_{B} (GeV/c^{2})");
   hMC->SetYTitle("Entries / (20 MeV/c^{2})");
   hMC->GetXaxis()->CenterTitle();
   hMC->GetYaxis()->CenterTitle();
   hMC->SetTitleOffset(1.,"Y");
   hMC->SetAxisRange(0,hMC->GetMaximum()*1.2,"Y");

   background->Draw("same");   

   mass->SetRange(5,6);
   mass->Draw("same");
   mass->SetLineStyle(2);
   mass->SetFillStyle(3004);
   mass->SetFillColor(2);
   f->Draw("same");

   double yield = mass->Integral(5,6)/0.02;

   // Draw the legend:)   
   TLegend *leg = myLegend(0.50,0.5,0.86,0.92);
   leg->AddEntry(hMC,"CMS Preliminary","");
   leg->AddEntry(hMC,"p+Pb #sqrt{s_{NN}}= 5.02 TeV","");
   leg->AddEntry(hMC,Form("%.0f<p_{T}^{B}<%.0f GeV/c",ptmin,ptmax),"");
   leg->AddEntry(f,"Fit","l");
   leg->AddEntry(mass,"Signal","f");
   leg->AddEntry(background,"Combinatorial Background","l");
   leg->Draw();
   TLegend *leg2 = myLegend(0.44,0.33,0.89,0.50);
   leg2->AddEntry(hMC,"B meson","");
   leg2->AddEntry(hMC,Form("N_{B}=%.0f", yield),"");
   leg2->Draw();

   cSignal->SaveAs("NonpromptBplus/BMassSignal.gif");

   return mass;
}

void bfractionVsCent(char *tagger="discr_ssvHighEff", double workingPoint=2., int fixCL=0, char *taggerName="ssvHighEff", float ptlo=80, float pthi=100, float etalo=0., float etahi=2.) {

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetLabelFont(43,"xyz");
  gStyle->SetLabelSize(20,"xyz");
  gStyle->SetTitleFont(43,"xyz");
  gStyle->SetTitleSize(26,"xyz");
  gStyle->SetTitleOffset(1.0,"xy"); 
  gROOT->ForceStyle(1);
  
  int doLTJP=1;
  int doLTCSV=0;

  const int nBins = 2;
  double centBin[nBins+1] = {0,30,100};
  //const int nBins = 3;
  //double centBin[nBins+1] = {0,20,50,100};
 
  
  Double_t bPurMC, bPurData, bEffMC, bEffDataLTJP, bEffDataLTCSV, taggedFracData, bFracMC, bFracData, bFracDataLTJP, bFracDataLTCSV, bFracJPdirect;
  Double_t bPurMCError, bPurDataError, bEffMCError, bEffDataLTJPError, bEffDataLTCSVError, taggedFracDataError, bFracMCError, bFracDataError, bFracDataLTJPError, bFracDataLTCSVError, bFracJPdirectError;
  Enumerations numbers;
  
  TH1D *hBPurityData = new TH1D("hBPurityData","hBPurityData;Centrality;B-Tagging purity",nBins,centBin);
  TH1D *hBPurityMC = new TH1D("hBPurityMC","hBPurityMC;Centrality;B-Tagging purity",nBins,centBin);
  
  TH1D *hBEfficiencyMC = new TH1D("hBEfficiencyMC","hBEfficiencyMC;Centrality;B-Tagging efficiency",nBins,centBin);
  TH1D *hBEfficiencyDataLTJP = new TH1D("hBEfficiencyDataLTJP","hBEfficiencyDataLTJP;Centrality;B-Tagging efficiency",nBins,centBin);
  TH1D *hBEfficiencyDataLTCSV = new TH1D("hBEfficiencyDataLTCSV","hBEfficiencyDataLTCSV;Centrality;B-Tagging efficiency",nBins,centBin);
  
  TH1D *hBFractionMC = new TH1D("hBFractionMC","hBFractionMC;Centrality;B-jet fraction",nBins,centBin);
  TH1D *hBFractionData = new TH1D("hBFractionData","hBFractionData;Centrality;B-jet fraction",nBins,centBin);
  TH1D *hBFractionDataLTJP = new TH1D("hBFractionDataLTJP","hBFractionDataLTJP;Centrality;B-jet fraction",nBins,centBin);
  TH1D *hBFractionDataLTCSV = new TH1D("hBFractionDataLTCSV","hBFractionDataLTCSV;Centrality;B-jet fraction",nBins,centBin);
  TH1D *hBFractionJPdirect = new TH1D("hBFractionJPdirect","hBFractionJPdirect;Centrality;B-jet fraction",nBins,centBin);
  
  int ncol=1;
  int nrow=1;

  if(nBins==3||nBins==2){
    ncol=nBins;
  }
  if(nBins==4){
    ncol=nBins/2;
    nrow=nBins/2;
  }

  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  //c1->Divide(ncol,nrow,0,0);
  c1->Divide(ncol,nrow);
  TCanvas *c2=new TCanvas("c2","c2",1200,600);
  //c2->Divide(ncol,nrow,0,0);
  c2->Divide(ncol,nrow);
  TCanvas *c3=new TCanvas("c3","c3",1200,600);
  //c3->Divide(ncol,nrow,0,0);
  c3->Divide(ncol,nrow);
  TCanvas *c4=new TCanvas("c4","c4",1200,600);
  //c4->Divide(ncol,nrow,0,0);
  c4->Divide(ncol,nrow);

  TCanvas *cCount = new TCanvas("cCount","cCount",600,600);

  for (int n=0;n<nBins;n++) {

    cout<<"Processing jet centrality bin ["<<centBin[n]<<","<<centBin[n+1]<<"] ..."<<endl;
    cCount->cd();
    cout<<"centBin[n]: "<<centBin[n]<<" centBin[n+1]: "<<centBin[n+1]<<" tagger: "<<tagger<<" workingPoint: "<<workingPoint<<" ptlo: "<<ptlo<<" pthi: "<<pthi<<" etalo: "<<etalo<<" etahi: "<<etahi<<endl;
    numbers = count(centBin[n],centBin[n+1],tagger,workingPoint,ptlo,pthi,etalo,etahi);
    
    c1->cd(n+1);
    c1->GetPad(n+1)->SetLogy();
    RooRealVar fitSvtxmTag = bfractionFit(fixCL,"svtxm",0,6,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,workingPoint,6,"b-tagged sample (SSVHE > 2)",9e3);
    //RooRealVar fitSvtxmTag = bfractionFit(fixCL,"svtxm",0,6,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,workingPoint,10,Form("b-tagged sample (%s at %.1f)",taggerName,workingPoint));
    //RooRealVar fitJpDirect = bfractionFit(fixCL,"discr_prob",0,3,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,-2,10,"inclusive sample",5e4);
    c2->cd(n+1);
    c2->GetPad(n+1)->SetLogy();
    RooRealVar fitJpDirect = bfractionFit(fixCL,"discr_prob",0.0,3.,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,"discr_prob",0.,3.,"inclusive sample",4e5);

    if (doLTJP) {
      c3->cd(n+1);
      c3->GetPad(n+1)->SetLogy();
      RooRealVar fitJpBeforetag = bfractionFit(fixCL,"discr_prob",0.0,3.,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,"discr_prob",0,3.,"jets with JP info",4e5);
      c4->cd(n+1);
      c4->GetPad(n+1)->SetLogy();
      RooRealVar fitJpTag = bfractionFit(fixCL,"discr_prob",0.0,3.,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,workingPoint,6,"b-tagged sample (SSVHE > 2)",4e5);
    } 
    if (doLTCSV) {
      RooRealVar fitCsvBeforetag = bfractionFit(fixCL,"discr_csvSimple",0,1,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,-2,10,"jets with CSV info",4e5);
      RooRealVar fitCsvTag = bfractionFit(fixCL,"discr_csvSimple",0,1,centBin[n],centBin[n+1],ptlo,pthi,etalo,etahi,tagger,workingPoint,10,"b-tagged sample (SSVHE > 2)",4e5);
    } 

    taggedFracData = numbers.nTaggedJetsData / (numbers.nTaggedJetsData+numbers.nUntaggedJetsData);
    taggedFracDataError = fracError(numbers.nTaggedJetsData,numbers.nUntaggedJetsData,numbers.nTaggedJetsDataError,numbers.nUntaggedJetsDataError);
    
    //*  --- b-tagging purity --- 

    bPurMC = numbers.nTaggedBjetsMC / numbers.nTaggedJetsMC;
    cout<<" bPurMC "<<bPurMC<<" numbers.nTaggedBjetsMC "<<numbers.nTaggedBjetsMC<<" numbers.nTaggedJetsMC "<<numbers.nTaggedJetsMC<<endl;
    bPurMCError = fracError(numbers.nTaggedBjetsMC,numbers.nTaggedNonBjetsMC,numbers.nTaggedBjetsMCError,numbers.nTaggedNonBjetsMCError);
//.........这里部分代码省略.........