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

TH2D* GetJetCorrFunc2D_doublegaussian(int itrg, int jass)
{
                    TH2D* hcorr = (TH2D*)GetRawCorrFunc2D_ratio(itrg,jass);
                    TH2D* hcorr_clone = (TH2D*)hcorr->Clone(Form("corr_clone_itrg%d_jass%d",itrg,jass));
                    hcorr_clone->Reset();
		    for(int ietabin=1;ietabin<=hcorr->GetNbinsX();ietabin++)
		    {
                      TH1D* hcorrphi = (TH1D*)hcorr->ProjectionY(Form("corrphi_%d",ietabin),ietabin,ietabin,"e");  
                      float min = hcorrphi->GetMinimum();
                      hcorrphi->SetAxisRange(-1,1,"X");
                      float nearmax = hcorrphi->GetMaximum();
                      hcorrphi->SetAxisRange(PI-1,PI+1,"X");
                      float awaymax = hcorrphi->GetMaximum();

		      TF1* fitfunc = new TF1("fitfunc",doubleGaussian,-PI/2.,3.*PI/2.,5);		    
                      fitfunc->SetParameters(min,nearmax-min,0.3,awaymax-min,0.5);  
                      fitfunc->SetParLimits(0,0,100000);
                      fitfunc->SetParLimits(1,0,100000);
                      fitfunc->SetParLimits(2,0,100000);
                      fitfunc->SetParLimits(3,0,100000);
                      fitfunc->SetParLimits(4,0,100000);		    
		      for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%d",ietabin),"RNO");                     
		      float level = fitfunc->GetParameter(0);

		      for(int iphibin=1;iphibin<=hcorr->GetNbinsY();iphibin++) 
                        hcorr_clone->SetBinContent(ietabin,iphibin,hcorr->GetBinContent(ietabin,iphibin)-level);                      
                      delete fitfunc;
		    }

                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");
                    return hcorr_clone;
}

// ----------------------------------------------------------------------
void plotRatioRaw(string era = "2016BF") {
  gStyle->SetOptStat(0);

  if (era == "all") {
    plotRatioRaw("2011");
    plotRatioRaw("2012");
    plotRatioRaw("2016BF");
    plotRatioRaw("2016GH");
    return;
  }


  double eps(0.00001);
  TH1D *hresult = new TH1D("hresult", Form("%s", era.c_str()), 80, 0., 0.4);
  hresult->GetXaxis()->SetTitle("BDT > ");
  hresult->SetMarkerStyle(24); hresult->SetLineColor(kRed+2); hresult->SetMarkerColor(kRed+2);
  TH1D *hsresult = new TH1D("hsresult", Form("%s", era.c_str()), 80, 0.+eps, 0.4+eps);
  hsresult->SetMarkerStyle(25); hsresult->SetLineColor(kBlue);  hsresult->SetMarkerColor(kBlue);

  TH1D *hbdtcut = new TH1D("hbdtcut", Form("%s", era.c_str()), 80, 0., 0.4);
  hbdtcut->SetMarkerColor(kRed+2);  hbdtcut->SetMarkerStyle(20);
  TH1D *hsbdtcut = new TH1D("hsbdtcut", Form("%s", era.c_str()), 80, 0.+eps, 0.4+eps);
  hsbdtcut->SetMarkerColor(kBlue);  hsbdtcut->SetMarkerStyle(21);

  invertedMuonIDRaw(era, 0.15, hresult, hsresult);
  invertedMuonIDRaw(era, 0.20, hresult, hsresult);
  invertedMuonIDRaw(era, 0.25, hresult, hsresult);
  invertedMuonIDRaw(era, 0.30, hresult, hsresult);
  invertedMuonIDRaw(era, -99., hbdtcut, hsbdtcut);

  c0->Clear();
  double ymax(50.);
  double smax = hresult->GetMaximum();
  double dmax = hsresult->GetMaximum();
  double themax = (smax>dmax? smax: dmax);
  if (themax < 50.)  ymax = 50.;
  if (themax < 20.)  ymax = 30.;
  if (themax < 5.)  ymax = 10.;
  if (themax < 0.5) ymax = 1.0;
  hresult->SetMinimum(-ymax);
  hresult->SetMaximum(ymax);
  hresult->Draw();
  hsresult->Draw("same");
  hbdtcut->Draw("same");
  hsbdtcut->Draw("same");
  pl->DrawLine(0., 0., 0.4, 0.);

  TLegend *tle = new TLegend(0.25, 0.2, 0.50, 0.5);
  tle->SetFillStyle(0);
  tle->SetBorderSize(0);
  tle->SetHeader(era.c_str());
  tle->AddEntry(hresult, "yield at B0", "p");
  tle->AddEntry(hsresult, "yield at Bs", "p");
  tle->Draw();

  c0->SaveAs(Form("invertedMuonIDRaw-%s-result.pdf", era.c_str()));
}

void plotter::draw_output_pseudo(TH1* output_, TH1D* pseudotruth_, TH1D* mctruth_, bool norm, TString file_name){

  TH1* output = (TH1*) output_->Clone("output");
  TH1D* pseudotruth = (TH1D*) pseudotruth_->Clone("pseudotruth");
  TH1D* mctruth = (TH1D*) mctruth_->Clone("mctruth");

  double ymax_temp = 0;
  if(pseudotruth->GetMaximum() > ymax_temp) ymax_temp = pseudotruth->GetMaximum();
  if(mctruth->GetMaximum() > ymax_temp) ymax_temp = mctruth->GetMaximum();
  if(output->GetMaximum() > ymax_temp) ymax_temp = output->GetMaximum();
  double ymax = 1.5 * ymax_temp;


  pseudotruth->SetTitle(" ");
  pseudotruth->GetYaxis()->SetRangeUser(0., ymax);
  pseudotruth->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  if(norm) pseudotruth->GetYaxis()->SetTitle("#frac{1}{#sigma} #frac{d#sigma}{dm_{jet}} [#frac{1}{GeV}]");
  else     pseudotruth->GetYaxis()->SetTitle("events");
  pseudotruth->GetYaxis()->SetTitleOffset(1.1);
  pseudotruth->GetXaxis()->SetTitleOffset(0.9);
  pseudotruth->GetYaxis()->SetTitleSize(0.05);
  pseudotruth->GetXaxis()->SetTitleSize(0.05);
  pseudotruth->GetYaxis()->SetNdivisions(505);

  pseudotruth->SetLineWidth(4);
  pseudotruth->SetLineColor(kRed);
  mctruth->SetLineWidth(3);
  mctruth->SetLineStyle(2);
  mctruth->SetLineColor(kBlue);

  output->SetLineColor(kBlack);
  output->SetMarkerColor(kBlack);
  output->SetMarkerStyle(8);
  output->SetMarkerSize(1);

  TCanvas *c= new TCanvas("Particle Level","",600,600);
  gPad->SetLeftMargin(0.15);
  TGaxis::SetMaxDigits(3);
  pseudotruth->Draw("HIST SAME");
  mctruth->Draw("HIST SAME");
  output->Draw("E1 SAME");
  TLegend *l;
  if(mctruth->GetSize() > 20) l=new TLegend(0.2,0.6,0.4,0.85);
  else                       l=new TLegend(0.55,0.6,0.85,0.85);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  l->AddEntry(output,"pseudo data","pl");
  l->AddEntry(pseudotruth,"pseudo data truth","pl");
  l->AddEntry(mctruth,"MC truth","pl");
  l->SetTextSize(0.04);

  l->Draw();
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

// Takes as input a given process with corresponding luminosity, as well as sample CME, cuts, directory location, and output file name - saves corresponding histogram in root file.
void processtree(TString sample, Double_t weight, TLorentzVector CME, TString filename, TString directory, TString param, TString cut, std::clock_t start, Int_t all){

    //Open File
    TFile *f = new TFile(directory);
    TTree *ttree = (TTree*)f->Get(sample);

    TString recohistname, paramhistname;
    if (all==0){
        recohistname = "_recoilmassHIST";
        paramhistname = "_parameterHIST";
    }
    if (all==1){
        recohistname = "_recoilmassHIST_all";
        paramhistname = "_parameterHIST_all";
    }

    TString histname = sample + recohistname;
    TH1D *recoilmassHIST   = new TH1D(histname,histname,CME.M()/2,0,CME.M());
    //recoilmassHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, "recoilmass", cut);

    histname = sample + paramhistname;
    TH1D *parameterHIST   = new TH1D(histname,histname,CME.M()/2,1,-1);
    //parameterHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, param, cut);

    //Normalize Hist
    recoilmassHIST->Scale(1/recoilmassHIST->GetMaximum());
    parameterHIST->Scale(1/parameterHIST->GetMaximum()) ;

    //Normalize Hist
//    recoilmassHIST->Scale(weight);
//    parameterHIST->Scale(weight) ;

    // Writes yield of each histogram for each cut to text file
    ofstream myfile;
    Double_t yield = recoilmassHIST->Integral(); 
    myfile.open ("yield.txt", ios::app);
        cout << "Yield for " << sample << " for Cut " << param << " is " << yield << endl;
        myfile << param << " " << all << " " << sample << " " << yield << "\n"; //write to file
    myfile.close();

    if (recoilmassHIST->Integral() == 0){
        recoilmassHIST->Fill(0);
    }
    if (parameterHIST->Integral() == 0){
        parameterHIST->Fill(0);
    }

    TFile g(filename, "update");
    recoilmassHIST->Write();
    parameterHIST->Write();
    g.Close();
}

void plotter::draw_output_stat(TH1* output_, TH1* stat_, TH1D* truth_, bool norm, TString file_name){
  // std::vector<double> sys = get_sys_errors();
  // TH1* output_sys = add_error_bar(output, sys);

  TH1* output = (TH1*) output_->Clone("output");
  TH1* stat = (TH1*) stat_->Clone("stat");
  TH1D* truth = (TH1D*) truth_->Clone("truth");

  TCanvas *c = new TCanvas("c","",600,600);
  double ymax;
  gPad->SetLeftMargin(0.15);

  if(truth->GetMaximum() > output->GetMaximum()) ymax = 1.5 * truth->GetMaximum();
  else ymax = 1.5 * output->GetMaximum();
  TGaxis::SetMaxDigits(3);
  output->SetTitle(" ");
  output->GetYaxis()->SetRangeUser(0., ymax);
  output->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  if(norm) output->GetYaxis()->SetTitle("#frac{1}{#sigma} #frac{d#sigma}{dm_{jet}} [#frac{1}{GeV}]");
  else output->GetYaxis()->SetTitle("events");
  output->GetYaxis()->SetTitleOffset(1.1);
  output->GetXaxis()->SetTitleOffset(0.9);
  output->GetYaxis()->SetTitleSize(0.05);
  output->GetXaxis()->SetTitleSize(0.05);
  output->GetYaxis()->SetNdivisions(505);
  output->SetLineColor(kBlack);
  output->SetMarkerColor(kBlack);
  output->SetMarkerStyle(8);
  output->SetMarkerSize(1);
  output->Draw("E1");
  stat->SetLineColor(kBlack);
  stat->SetMarkerColor(kBlack);
  stat->SetMarkerStyle(8);
  stat->SetMarkerSize(1);
  gStyle->SetEndErrorSize(5);
  truth->SetLineWidth(3);
  truth->SetLineColor(kRed);
  truth->SetLineStyle(2);
  truth->Draw("HIST SAME");
  stat->Draw("E1 SAME");
  output->Draw("E1 SAME");

  TLegend *l=new TLegend(0.5,0.65,0.85,0.85);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  l->AddEntry(output,"data unfolded","pl");
  l->AddEntry(truth,"MC particle level","pl");
  l->SetTextSize(0.04);
  l->Draw();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

void plotter::draw_output_smear(std::vector<TH1*> output_, TH1D* truth_, TString file_name){
  // std::vector<double> sys = get_sys_errors();
  // TH1* output_sys = add_error_bar(output, sys);

  std::vector<TH1*> outputs;
  for(int i=0; i<output_.size(); i++){
    TH1* hist = (TH1*) output_[i]->Clone();
    outputs.push_back(hist);
  }
  TH1D* truth = (TH1D*) truth_->Clone("truth");

  TCanvas *c = new TCanvas("c","",600,600);
  double ymax;
  gPad->SetLeftMargin(0.15);

  if(truth->GetMaximum() > outputs[0]->GetMaximum()) ymax = 1.5 * truth->GetMaximum();
  else ymax = 1.5 * outputs[0]->GetMaximum();
  TGaxis::SetMaxDigits(3);
  for(auto output: outputs){
    output->SetTitle(" ");
    output->GetYaxis()->SetRangeUser(0., ymax);
    output->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
    output->GetYaxis()->SetTitle("events");
    output->GetYaxis()->SetTitleOffset(1.1);
    output->GetXaxis()->SetTitleOffset(0.9);
    output->GetYaxis()->SetTitleSize(0.05);
    output->GetXaxis()->SetTitleSize(0.05);
    output->GetYaxis()->SetNdivisions(505);
    output->SetLineColor(kBlack);
    output->SetMarkerColor(kBlack);
    output->SetMarkerStyle(8);
    output->SetMarkerSize(1);
    output->Draw("E1 SAME");
  }
  gStyle->SetEndErrorSize(5);
  truth->SetLineWidth(3);
  truth->SetLineColor(kRed);
  truth->SetLineStyle(2);
  truth->Draw("HIST SAME");
  TLegend *l;
  if(truth->GetSize() > 20) l=new TLegend(0.2,0.6,0.4,0.85);
  else                       l=new TLegend(0.5,0.6,0.85,0.85);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  l->AddEntry(outputs[0],"data unfolded","pl");
  l->AddEntry(truth,"MC particle level","pl");
  l->SetTextSize(0.04);
  l->Draw();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

//------------------------------------------------------------------------
void fitSlices(TH2* hCorr, TF1* func){

   int nBins = hCorr->GetNbinsX();

   TH1D* hMean = new TH1D(Form("%s_1",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());
   TH1D* hSigma = new TH1D(Form("%s_2",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());

   for(int i = 1; i < nBins+1; i++){
      int bin = nBins - i;
      TH1D* h = hCorr->ProjectionY(Form("%s_bin%d",hCorr->GetName(),bin),i,i);

      func->SetParameter(0,h->GetMaximum());
      func->SetParameter(1,h->GetMean());
      func->SetParameter(2,h->GetRMS());

      h->Fit(func,"Q");

      /*hMean->SetBinContent(i,func->GetParameter(1));
      hMean->SetBinError(i,func->GetParError(1));
      hSigma->SetBinContent(i,func->GetParameter(2));
      hSigma->SetBinError(i,func->GetParError(2));*/

      hMean->SetBinContent(i,h->GetMean());
      hMean->SetBinError(i,func->GetParError(1)); //errors are not use later for the actual table
      hSigma->SetBinContent(i,h->GetRMS());
      hSigma->SetBinError(i,func->GetParError(2));

   }
}

// ----------------------------------------------------------------------
void plotRatio(string era = "2016BF") {
  gStyle->SetOptStat(0);

  if (era == "all") {
    plotRatio("2011");
    plotRatio("2012");
    plotRatio("2016BF");
    plotRatio("2016GH");
    return;
  }


  TH1D *hresult = new TH1D("hresult", Form("n(gauss(5.28))/n(gauss(5.37)) %s", era.c_str()), 40, 0., 0.4);
  hresult->GetXaxis()->SetTitle("BDT > ");
  TH1D *hbdtcut = new TH1D("hbdtcut", Form("n(gauss(5.28))/n(gauss(5.37)) %s", era.c_str()), 40, 0., 0.4);
  hbdtcut->SetMarkerColor(kBlue);
  invertedMuonID(era, 0.15, hresult);
  invertedMuonID(era, 0.20, hresult);
  invertedMuonID(era, 0.25, hresult);
  invertedMuonID(era, 0.30, hresult);
  invertedMuonID(era, -99., hbdtcut);

  c0->Clear();
  double ymax(50.);
  if (hresult->GetMaximum() < 5.)  ymax = 10.;
  if (hresult->GetMaximum() < 0.5) ymax = 1.0;
  hresult->SetMinimum(-ymax);
  hresult->SetMaximum(ymax);
  hresult->Draw();
  hbdtcut->Draw("same");
  pl->DrawLine(0., 0., 0.4, 0.);
  c0->SaveAs(Form("invertedMuonID-%s-result.pdf", era.c_str()));
}

void fitSlices(TH2* hCorr, TF1* func){

   int nBins = hCorr->GetNbinsX();

   TH1D* hMean = new TH1D(Form("%s_1",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());
   TH1D* hSigma = new TH1D(Form("%s_2",hCorr->GetName()),"",nBins,hCorr->GetXaxis()->GetXmin(),hCorr->GetXaxis()->GetXmax());

   for(int i = 1; i < nBins+1; ++i){
      int bin = nBins - i;
      TH1D* h = hCorr->ProjectionY(Form("%s_bin%d",hCorr->GetName(),bin),i,i);

      func->SetParameter(0,h->GetMaximum());
      func->SetParameter(1,h->GetMean());
      func->SetParameter(2,h->GetRMS());

      if(useFits) h->Fit(func);

      hMean->SetBinContent(i,func->GetParameter(1));
      hMean->SetBinError(i,func->GetParError(1));
      hSigma->SetBinContent(i,func->GetParameter(2));
      hSigma->SetBinError(i,func->GetParError(2));
      
      if(onlySaveTable){
	 h->Delete();
      }
   }
}

void loglikdistrib(Int_t ntrials = 10000, Bool_t print = kFALSE)
{
  // compute distribution of log likelihood value
  TH1D * hmc   = gStack[gPadNr][gOrder[gPadNr][0]];
  TH1D * hdata = gStack[gPadNr][gMaxProcess-1];
  Int_t nbins = hmc->GetNbinsX();
  Double_t loglik = loglikelihood(hmc, hdata, 1, nbins);
  TH1D * htest = new TH1D(*hdata);
  TH1D * lldistrib = new TH1D("lldistrib", "log(Likelihood) distribution", 
			      1000, loglik-200, loglik+200);
  setopt(lldistrib);
  for (Int_t n = 0; n < ntrials; n++) {
    // generate poisson around theorie
    for (Int_t i = 1; i <= nbins; i++) {
      htest->SetBinContent(i, gRandom->Poisson(hmc->GetBinContent(i)));
    }
    lldistrib->Fill(loglikelihood(hmc, htest, 1, nbins));
  }
  TCanvas * llcanvas = new TCanvas("llcanvas", "Log(Likelihood) distribution", 
				   40, 40, 800, 600);
  setopt(llcanvas);
  lldistrib->SetFillColor(kYellow);
  lldistrib->Draw();
  lldistrib->GetYaxis()->SetTitle("Anzahl Ereignisse");
  lldistrib->GetXaxis()->SetTitle("-ln L");
  // autozoom
  Int_t lowbin = 1;
  while (lldistrib->GetBinContent(lowbin) == 0)
    lowbin++;
  Int_t highbin = lldistrib->GetNbinsX();
  while (lldistrib->GetBinContent(highbin) == 0)
    highbin--;
  lldistrib->SetAxisRange(lldistrib->GetBinLowEdge(lowbin), 
			  lldistrib->GetBinLowEdge(highbin));
  TH1D * hworse = (TH1D *) lldistrib->Clone();
  for (Int_t nbin = 1; nbin < 501; nbin++) {
    hworse->SetBinContent(nbin, 0);
  }
  hworse->SetFillColor(95);
  hworse->Draw("same");
  Double_t pvalue = lldistrib->Integral(501,1000) / lldistrib->Integral();
  TLatex * tex = new TLatex(0.18, 0.96, Form("-ln L_{obs} = %5.2f", loglik));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  tex = new TLatex(0.18, 0.86, Form("CL_{obs} = %.3f", pvalue));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  TLine * l = new TLine(loglik, 0, loglik, lldistrib->GetMaximum());
  l->SetLineWidth(3);
  l->SetLineColor(kBlue);
  l->Draw();
  llcanvas->Modified();
  llcanvas->Update();
  if (print)
    llcanvas->Print("lldistrib.pdf");
  cd(gPadNr+1);
}

void ReFit(int lo, int hi) {
	TF1* floo = (TF1*) gROOT->FindObject("floo");
	if (floo == 0) return;
	TH1D* hCe = (TH1D*) gROOT->FindObject("hCe");
	if (hCe == 0) return;
	hCe->GetXaxis()->SetRangeUser(lo,hi);
	double max = hCe->GetMaximum();
	int binmax = hCe->GetMaximumBin();
	double xmax = hCe->GetXaxis()->GetBinCenter(binmax);
	floo->SetParameters(max,xmax,10,max*0.1,50,20,5);
	hCe->Fit(floo,"Q","",lo,hi);
}

void plotTurnOn(TTree* inttree, TString triggerpass, TString variable, TString varname, TString varlatex)
{
  if(varname=="vtxprob")
    {
      BIN_MIN = 0;
      BIN_MAX = 1;
    }
  else if(varname=="ffls3d")
    {
      BIN_MIN = 0;
      BIN_MAX = 50;
    }
  else if(varname=="cosalpha")
    {
      BIN_MIN = 0.9;
      BIN_MAX = 1;
    }
  TH1D* hAll = new TH1D(Form("h%s_%s_All",triggerpass.Data(),varname.Data()),Form(";%s;Probability",varlatex.Data()),BIN_NUM,BIN_MIN,BIN_MAX);
  inttree->Project(Form("h%s_%s_All",triggerpass.Data(),varname.Data()),variable,prefilter);
  TH1D* hMBseed = new TH1D(Form("h%s_%s_MBseed",triggerpass.Data(),varname.Data()),Form(";%s;Probability",varlatex.Data()),BIN_NUM,BIN_MIN,BIN_MAX);
  inttree->Project(Form("h%s_%s_MBseed",triggerpass.Data(),varname.Data()),variable,Form("%s&&%s",prefilter.Data(),triggerpass.Data()));
  //cout<<hAll->Integral()<<endl;
  hAll->Scale(1./hAll->Integral());
  hMBseed->Scale(1./hMBseed->Integral());
  hAll->SetStats(0);
  hMBseed->SetStats(0);
  hAll->SetMaximum(hMBseed->GetMaximum()*1.3);
  hAll->SetLineWidth(2);
  hAll->SetLineColor(kBlue-7);
  hAll->SetFillColor(kBlue-7);
  hAll->SetFillStyle(3001);
  hMBseed->SetLineWidth(2);
  hMBseed->SetLineColor(kRed);
  hMBseed->SetFillColor(kRed);
  hMBseed->SetFillStyle(3004);
  TCanvas* c = new TCanvas(Form("c%s_%s",triggerpass.Data(),varname.Data()),"",500,500);
  hAll->Draw();
  hMBseed->Draw("same");
  TLatex* tex = new TLatex(0.18,0.96,triggerpass);
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->Draw();
  TLegend* leg = new TLegend(0.60,0.82,0.92,0.93);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->AddEntry(hAll,"all","f");
  leg->AddEntry(hMBseed,"pass trigger","f");
  leg->Draw();
  c->SaveAs(Form("triggerturnonPlots/pthat%.0f/c%s_%s.pdf",pthat,triggerpass.Data(),varname.Data()));

}

void Drawbbc(){
	TFile *f = new TFile("merged.root");
	TH1D* hbbcs;
	TH1D* hbbcn;
	TCanvas *c1;
	TCanvas *c2;
	TLegend *leg = new TLegend(0.5,0.7,0.7,0.85);
	leg->SetBorderSize(0);
	leg->SetFillColor(0);
	leg->SetTextSize(0.048);
	TH2F* bbcsbbcn = (TH2F*)f->Get(Form("bbcsbbcn_0"));
	for(int xcent=1; xcent < 6; xcent++){
		TH1F* bbcsbbcndis_t = (TH1F*)f->Get(Form("bbcsbbcn_%d",xcent));
		bbcsbbcn->Add(bbcsbbcndis_t);
	}
	hbbcs = (TH1D*)bbcsbbcn->ProjectionX("bbcs",0,-1);
	hbbcn = (TH1D*)bbcsbbcn->ProjectionY("bbcn",0,-1);
	hbbcs->Rebin(4);
	hbbcn->Rebin(4);
	c1= new TCanvas();
	c1->cd();
	c1->SetLogy();
	SetTitle(*hbbcs,"bbc charge","# of events","");
	SetRange(*hbbcs,0,1e-1,100,hbbcn->GetMaximum()*5);
	SetStyle(*hbbcs,1.2,1,20,0,0);
	SetStyle(*hbbcn,1.2,2,24,0,0);
	hbbcs->Draw("P");
	hbbcn->Draw("Psame");
	leg->AddEntry(hbbcs,"bbc south");
	leg->AddEntry(hbbcn,"bbc north");
	leg->Draw("same");
	c1->Print(Form("bbcdis.png"));
	c2= new TCanvas();
	c2->cd();
	c2->SetLogy();
	hbbcs_norm = (TH1D*)hbbcs->Clone("hbbcs_norm");
	hbbcn_norm = (TH1D*)hbbcn->Clone("hbbcn_norm");
	hbbcs_norm->Scale(1./hbbcs_norm->Integral());
	hbbcn_norm->Scale(1./hbbcn_norm->Integral());
	SetRange(*hbbcs_norm,0,1e-9,100,1);
	hbbcs_norm->Draw("P");
	hbbcn_norm->Draw("Psame");
	leg->Draw("same");
	c2->Print(Form("bbcdis_norm.png"));	
	c3= new TCanvas();
	c3->cd();
	c3->SetLogz();
	SetTitle(*bbcsbbcn,"bbc charge south","bbc charge north","# of events");
	SetRange(*bbcsbbcn,0,0,100,100);
	bbcsbbcn->Draw("colz");
	c3->Print("bbcsbbcn.png");
}

void CompareBranch(string MCfilename, string CDfilename, string MCbranchname, string CDbranchname, string xtitle, string unit, string plotname, string MCcuts, string CDcuts, string MCweight, string CDweight, double xlow, double xup, int nbins)
{
  TH1D*  MChist = MakeBranchPlot(MCfilename,MCbranchname,MCcuts,MCweight,xlow,xup,nbins);
  TH1D*  CDhist = MakeBranchPlot(CDfilename,CDbranchname,CDcuts,CDweight,xlow,xup,nbins);
  MChist->Scale(1./MChist->Integral());
  CDhist->Scale(1./CDhist->Integral());
  MChist->SetDrawOption("B");
  MChist->SetFillColor(kOrange);
  MChist->SetLineColor(kOrange);
  MChist->SetMaximum(MChist->GetMaximum()*1.3);
  MChist->SetMinimum(0);
  // Draw everything
  plotmaker plotter(MChist);
  plotter.SetTitle(xtitle, unit);
  TCanvas* plot = plotter.Draw();
  CDhist->Draw("sameE1");
  plot->SaveAs((plotname+".pdf").c_str());
}

void plotter::draw_bias(TH1* output_, TH1D* truth_, TH1* bias_, TString file_name){
  TH1D* bias = (TH1D*) bias_->Clone("bias");
  TH1D* truth = (TH1D*) truth_->Clone("truth");
  TH1D* output = (TH1D*) output_->Clone("output");

  double ymax_temp = 0;
  if(truth->GetMaximum() > ymax_temp) ymax_temp = truth->GetMaximum();
  if(bias->GetMaximum() > ymax_temp) ymax_temp = bias->GetMaximum();
  if(output->GetMaximum() > ymax_temp) ymax_temp = output->GetMaximum();
  double ymax = 1.5 * ymax_temp;

  TCanvas *c= new TCanvas("Bias Distribution","",600,600);
  TGaxis::SetMaxDigits(3);
  gPad->SetLeftMargin(0.15);
  truth->SetTitle(file_name);
  truth->GetYaxis()->SetRangeUser(0., ymax);
  truth->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  truth->GetYaxis()->SetTitle("events");
  truth->GetYaxis()->SetTitleOffset(1.1);
  truth->GetXaxis()->SetTitleOffset(0.9);
  truth->GetYaxis()->SetTitleSize(0.05);
  truth->GetXaxis()->SetTitleSize(0.05);
  truth->GetYaxis()->SetNdivisions(505);

  truth->SetLineWidth(4);
  truth->SetLineColor(kRed);
  bias->SetLineWidth(4);
  bias->SetLineColor(kAzure+7);
  output->SetLineColor(kBlack);
  output->SetMarkerColor(kBlack);
  output->SetMarkerStyle(8);
  output->SetMarkerSize(1);

  truth->Draw("HIST SAME");
  bias->Draw("HIST SAME");
  output->Draw("E1 SAME");
  TLegend *l=new TLegend(0.55,0.6,0.85,0.85);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  l->AddEntry(output,"unfolded result","pl");
  l->AddEntry(truth,"truth","pl");
  l->AddEntry(bias,"bias distribution","pl");
  l->SetTextSize(0.04);
  l->Draw();
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}