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

void ratioPlots( TCanvas* c1, TH1* h_r, TH1* h_i, 
		 string xTitle, string yTitle, 	string savePath, 
		 double fitMin=-100000, double fitMax=100000, bool doubleColFit=0 ){

	double xMaximum = h_r->GetXaxis()->GetBinUpEdge(h_r->GetXaxis()->GetLast());
	double xMinimum = h_r->GetXaxis()->GetBinLowEdge(h_r->GetXaxis()->GetFirst());
	double yMaximum;
	double yMinimum;

	h_i->Sumw2();
	h_r->Sumw2();

	TLine* line1 = new TLine(xMinimum,1,xMaximum,1);
	line1->SetLineColor(1);
	line1->SetLineWidth(2);
	line1->SetLineStyle(7);
	
	TF1* fpol1  = new TF1("fpol1", "pol1", fitMin, fitMax);	
	fpol1->SetLineColor(2);
	fpol1->SetLineWidth(3);
	fpol1->SetLineStyle(7);
	
	TH1* hRatio = (TH1*)h_r->Clone("clone_record");
	hRatio->Divide(h_i);
	yMaximum = hRatio->GetMaximum();
	yMinimum = hRatio->GetMinimum(0);
	hRatio->GetYaxis()->SetRangeUser(yMinimum/2.5,yMaximum+yMaximum/5);
	hRatio->SetXTitle(xTitle.c_str());
	hRatio->SetYTitle(yTitle.c_str());
	hRatio->SetLineColor(9); 
	hRatio->SetLineWidth(2); 
	hRatio->SetMarkerStyle(8); 
	hRatio->Draw("e");
	hRatio->Fit("fpol1", "L");
	line1->Draw("SAME");
	
	if(doubleColFit){
		double p0=fpol1->GetParameter(0);
		double p1=fpol1->GetParameter(1);
		double endPoint=double(fitMax*p1)+p0;
		double p1new=(endPoint-1)/(fitMax-fitMin);
		char fun[100], text[100];
		sprintf(fun,"x*(%f)+1",p1new);	
		sprintf(text,"Tangent: %f",p1new);	
		TF1* fnew = new TF1("fnew", fun, fitMin, fitMax);
		fnew->SetLineColor(2);
		fnew->SetLineWidth(3);
		fnew->Draw("SAME");
	
			
		TText* Title = new TText( fitMax/12, yMinimum, text);
		Title->SetTextColor(2);
		Title->SetTextSize(0.035);
		Title->Draw("SAME");
	}

	c1->SaveAs(savePath.c_str());
	c1->cd();
}

void test2bMBcorr(const Char_t *datatag, const Char_t *mctag, Bool_t bayes = kFALSE)
{

  TCanvas *c = new TCanvas("c", "c");
  c->Divide(4, 3);
  for (Int_t i = 0; i < 10; i++) {
    TH1 *hnomb = UnfoldMe_TAG(datatag, mctag, "clist", i, kFALSE, kFALSE, kFALSE, 1., 4, AliUnfolding::kPowerLaw, 5000., bayes);
    TH1 *hmb = UnfoldMe_TAG(datatag, mctag, "clist", i, kTRUE, kFALSE, kFALSE, 1., 4, AliUnfolding::kPowerLaw, 5000., bayes);
    hnomb->Divide(hmb);
    c->cd(i+1)->DrawFrame(1., 0.9, 100., 1.1);
    c->cd(i+1)->SetLogx();
    hnomb->Draw("same");
    c->Update();
  }

}

  //____________________________________________________________________
  TH1* One(TDirectory* newDir, TDirectory* oldDir, Double_t c1, Double_t c2)
  {
    TString name;
    name.Form("cent%03dd%02d_%03dd%02d",
	      Int_t(c1), Int_t(c1*100)%100,
	      Int_t(c2), Int_t(c2*100)%100);
    TDirectory* newSubDir = GetD(newDir, name);
    TDirectory* oldSubDir = GetD(oldDir, name);
    if (!newSubDir || !oldSubDir) return 0;
    Int_t newDim = 0;
    if      (TString(newDir->GetName()).Contains("etaipz")) newDim = 3;
    else if (TString(newDir->GetName()).Contains("eta"))    newDim = 2;
    else if (TString(newDir->GetName()).Contains("const"))  newDim = 1;
    Int_t oldDim = 0;
    if      (TString(oldDir->GetName()).Contains("etaipz")) oldDim = 3;
    else if (TString(oldDir->GetName()).Contains("eta"))    oldDim = 2;
    else if (TString(oldDir->GetName()).Contains("const"))  oldDim = 1;

    TDirectory* newSubSubDir = GetD(newSubDir, Form("results%dd",newDim));
    TDirectory* oldSubSubDir = GetD(oldSubDir, Form("results%dd",oldDim));
    if (!newSubSubDir || !oldSubSubDir) return 0;

    TH1* newRes = GetH1(newSubSubDir, "result");
    TH1* oldRes = GetH1(oldSubSubDir, "result");
    if (!newRes || !oldRes) return 0;

    TH1* ratio = static_cast<TH1*>(newRes->Clone(name));
    ratio->SetDirectory(0);
    ratio->SetTitle(Form("%5.1f - %5.1f%%", c1, c2));
    ratio->SetYTitle("New / Old");
    ratio->Divide(oldRes);
    fMin = TMath::Min(fMin, ratio->GetMinimum());
    fMax = TMath::Max(fMax, ratio->GetMaximum());

    Printf("Calculated %s/%s", newDir->GetName(), oldDir->GetName());
    if (!fLegend) return ratio;

    
    TLegendEntry* e =
      fLegend->AddEntry("", Form("%3.0f - %3.0f%%", c1, c2), "f");
    e->SetFillStyle(1001);
    e->SetFillColor(ratio->GetMarkerColor());
  
    return ratio;
  }

void plotter::draw_delta_rel(TH1* hist_, TH1* result_, TString file_name){
  TH1* hist = (TH1*) hist_->Clone("hist");
  TH1* result = (TH1*) result_->Clone("result");
  hist->Divide(result);
  hist->Scale(100);
  TCanvas *c= new TCanvas("Particle Level","",600,600);
  gPad->SetLeftMargin(0.15);
  hist->SetTitle(file_name);
  hist->GetXaxis()->SetTitle("Leading-jet mass [GeV]");
  hist->GetYaxis()->SetTitle("relative uncertainty [%]");
  hist->GetYaxis()->SetTitleOffset(1.5);
  hist->GetYaxis()->SetNdivisions(505);
  hist->SetFillColor(810);
  hist->SetLineColor(810);
  hist->Draw("HIST");
  gPad->RedrawAxis();
  c->SaveAs(directory + file_name + ".pdf");
  delete c;
}

TH1* compRatioHistogram(const std::string& ratioHistogramName, const TH1* numerator, const TH1* denominator)
{
  TH1* histogramRatio = 0;
  
  if ( numerator->GetDimension() == denominator->GetDimension() &&
       numerator->GetNbinsX() == denominator->GetNbinsX() ) {
    histogramRatio = (TH1*)numerator->Clone(ratioHistogramName.data());
    histogramRatio->Divide(denominator);
    
    int nBins = histogramRatio->GetNbinsX();
    for ( int iBin = 1; iBin <= nBins; ++iBin ){
      double binContent = histogramRatio->GetBinContent(iBin);
      histogramRatio->SetBinContent(iBin, binContent - 1.);
    }
    
    histogramRatio->SetLineColor(numerator->GetLineColor());
    histogramRatio->SetLineWidth(numerator->GetLineWidth());
    histogramRatio->SetMarkerColor(numerator->GetMarkerColor());
    histogramRatio->SetMarkerStyle(numerator->GetMarkerStyle());
    histogramRatio->SetMarkerSize(numerator->GetMarkerSize());
  }

  return histogramRatio;
}

void playsmoothiter(const Char_t *datatag, const Char_t *mctag, const Char_t *anatag, Int_t bin, Bool_t ismc = kFALSE)
{

  TCanvas *c = new TCanvas("cs", "cs");
  c->DrawFrame(0., 0.5, 100., 1.5);

  TH1 *href = UnfoldMe_TAG(datatag, mctag, anatag, bin, kTRUE, kFALSE, ismc, 1., 4);
  TH1 *hout = NULL;
  TH1 *heff_ev = new TH1F("heff_ev", "", 12, 0, 12);
  TH1 *heff_ch = new TH1F("heff_ch", "", 12, 0, 12);
  for (Int_t j = 2; j < 10; j+=2) {
  for (Int_t i = 5; i < 16; i+=2) {
    printf("%d %d\n", j, i);
    hout = UnfoldMe_TAG(datatag, mctag, anatag, bin, kTRUE, kFALSE, ismc, 0.1 * i, j);
    heff_ev->Fill(integr_eff);
    heff_ch->Fill(dndeta_eff);

    hout->Divide(href);
    c->cd();
    hout->Draw("same");
    c->Update();
  }}

}

void DrawTwoInPad(TVirtualPad* p,
                  Int_t        sub,
                  TH1*         h1,
                  TH1*         h2,
		  Bool_t       ratio,
		  Bool_t       logy=false,
		  Bool_t       legend=false)
{
  TVirtualPad* pp = p->cd(sub);
  pp->SetRightMargin(0.02);
  pp->SetLeftMargin(0.10);
  TVirtualPad* ppp = pp;
  if (ratio) {
    pp->Divide(1,2,0,0);
    ppp = pp->cd(1);				
    ppp->SetRightMargin(0.02);
  }
  if (logy) ppp->SetLogy();
  TH1* hs[] = { h1, h2, 0 };
  if (h1->GetMaximum() < h2->GetMaximum()) {
    hs[0] = h2;
    hs[1] = h1;
  }
  TH1** ph = hs;
  Double_t size = (ratio ? 0.1 : 0.05);
  Double_t off  = (ratio ? 0.6 : 0.5);
  h1->SetFillStyle(3004);
  h2->SetFillStyle(3005);
  while (*ph) {
    TString opt("hist");
    if (ph != hs) opt.Append(" same");

    TH1* copy = (*ph)->DrawCopy(opt);
    copy->GetXaxis()->SetLabelSize(2*size);
    copy->GetYaxis()->SetLabelSize(size);
    copy->GetYaxis()->SetTitleSize(size);
    copy->GetYaxis()->SetTitleOffset(off);
    copy->SetYTitle(copy->GetTitle());
    copy->SetTitle("");
    copy->SetDirectory(0);
    ph++;
  }
  TString s1 = h1->GetYaxis()->GetTitle();
  TString s2 = h2->GetYaxis()->GetTitle();
  
  if (legend) { 
    TLegend* l = new TLegend(0.6, 0.1, 0.9, 0.9);
    l->SetBorderSize(0);
    TLegendEntry* e = l->AddEntry("dummy", s1, "lf");
    l->SetFillColor(kWhite);
    e->SetFillColor(kBlack);
    e->SetFillStyle(h1->GetFillStyle());
    e = l->AddEntry("dummy", s2, "lf");
    e->SetFillColor(kBlack);
    e->SetFillStyle(h2->GetFillStyle());
    l->Draw();
  }
  if (!ratio) return;
  ppp = pp->cd(2);
  ppp->SetRightMargin(0.02);
  TH1* r = static_cast<TH1*>(h1->Clone(Form("ratio%s", h1->GetName())));
  r->SetDirectory(0);
  r->SetTitle("");
  r->GetXaxis()->SetLabelSize(size);
  r->GetYaxis()->SetLabelSize(size);
  r->GetYaxis()->SetTitleSize(0.9*size);
  r->GetYaxis()->SetTitleOffset(0.9*off);
  r->SetMarkerStyle(20);
  r->SetMarkerColor(h1->GetFillColor()+1);
  r->SetFillStyle(3007);
  r->SetYTitle(Form("#frac{%s}{%s}", s1.Data(), s2.Data()));

  // r->Add(h2, -1);
  // r->Divide(h1);
  if (!r->IsA()->InheritsFrom(TProfile::Class())) {
    r->GetSumw2()->Set(0); // r->Sumw2(false);
    h2->GetSumw2()->Set(0); // h2->Sumw2(false);
  }
  r->Divide(h2);
  Printf("%s", r->GetName());
  for (UShort_t bin = 1; bin <= r->GetNbinsX(); bin++) {
    Printf("  bin # %2d: Diff=%g+/-%g", bin, r->GetBinContent(bin),
	   r->GetBinError(bin));
    r->SetBinError(bin, 0);
  }
  r->GetSumw2()->Set(0); //r->Sumw2(false);
  r->SetMarkerSize(4);
  r->SetMaximum(r->GetMaximum()*1.2);
  r->SetMinimum(r->GetMinimum()*0.8);
  r->Draw("hist text30");  
  p->Modified();
  p->Update();
  p->cd();
}

TH1 *
UnfoldMe_MB2(const Char_t *data, const Char_t *mc, const Char_t *anatag, Int_t bin, Bool_t useMBcorr , Bool_t usecorrfit , Bool_t ismc , Float_t smooth , Int_t iter , Int_t regul , Float_t weight , Bool_t bayesian , Int_t nloop )
{

  // MF comments:
  // usedMBcorr: changes the matrix used for unfonding, from effMatrix to bin matrix (I think this is just to use mult dependent v s mb correction_)
  // usecorrfit: if I understand correctly, fits the response matrix and uses fit to extrapolate it

  
  TFile *fdt =0;
  if (ismc)
    fdt =  TFile::Open(data);
  else
    fdt = TFile::Open(data);
  TFile *fmc = TFile::Open(mc);

  TList *ldt = (TList *)fdt->Get(Form("%s", anatag));
  TList *lmc = (TList *)fmc->Get(Form("%s", anatag));
  
  TH2 *hmatdt = (TH2 *)ldt->FindObject(Form(responseMatrix, bin));
  TH2 *hmatmc = 0;
  if (useMBcorr){
     hmatmc = (TH2 *)lmc->FindObject("effMatrix");
     std::cout << "USING MB" << std::endl;
     
  }
  else {
    hmatmc = (TH2 *)lmc->FindObject(Form(responseMatrix, bin));
  }

  TH1 *hdata = hmatdt->ProjectionY("hdata");
//  TH1 *hdata = hmatdt->ProjectionY("htrue");  // For truth Only Calculations

  hdata->Sumw2();
  hdata->SetBinContent(1, 0.);
  hdata->SetBinError(1, 0.);
  //  hdata->Scale(1. / hdata->Integral());
  hdata->SetMarkerStyle(25);
  TH1 *htrue = hmatdt->ProjectionX("htrue");
  htrue->Sumw2();
  //  htrue->Scale(1. / htrue->Integral());
  htrue->SetMarkerStyle(7);
  htrue->SetMarkerColor(2);
  htrue->SetBinContent(1, 0.);
  htrue->SetBinError(1, 0.);

  TH2 *hcorr = (TH2 *)hmatmc->Clone("hcorr");
  TH1 *hinit = (TH1 *)hdata->Clone("hinit");
  TH1 *hresu = (TH1 *)hdata->Clone("hresu");
  TH1 *hbias = (TH1 *)hdata->Clone("hbias");
  hresu->SetMarkerStyle(20);
  hresu->SetMarkerColor(4);
  hresu->Reset();

  TH1 *hnum = hcorr->ProjectionY("hnum");
  TH1 *hden = hcorr->ProjectionY("hden");
  TH1 *heff = hcorr->ProjectionY("heff");
  hnum->Reset();
  hnum->Sumw2();
  hden->Reset();
  hden->Sumw2();
  heff->Reset();
  for (Int_t i = 0; i < heff->GetNbinsX(); i++) {
    Float_t int1 = hcorr->Integral(i + 1, i + 1, 0, -1);
    if (int1 <= 0.) continue;
    Float_t int2 = hcorr->Integral(i + 1, i + 1, 2, -1);
    hnum->SetBinContent(i + 1, int2);
    hnum->SetBinError(i + 1, TMath::Sqrt(int2));
    hden->SetBinContent(i + 1, int1);
    hden->SetBinError(i + 1, TMath::Sqrt(int1));
  }
  TCanvas *cEfficiency = new TCanvas("cEfficiency", "cEfficiency");
  cEfficiency->SetLogx();
  cEfficiency->SetLogy();

  heff->Divide(hnum, hden, 1., 1., "B");
  heff->Draw();
#if 0
  for (Int_t ii = 0; ii < heff->GetNbinsX(); ii++) {
    heff->SetBinContent(ii + 1, 1.);
    heff->SetBinError(ii + 1, 0.);
  }
#endif

  for (Int_t i = 0; i < hcorr->GetNbinsX(); i++) {
    hcorr->SetBinContent(i + 1, 1, 0.);
    hcorr->SetBinError(i + 1, 1, 0.);
  }
  for (Int_t i = 0; i < hcorr->GetNbinsY(); i++) {
    hcorr->SetBinContent(1, i + 1, 0.);
    hcorr->SetBinError(1, i + 1, 0.);
  }
  TH2 *hcorrfit = ReturnCorrFromFit(hcorr);
  // Docs from AliUnfolding
  //Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
  // unfolds with unfolding method fgMethodType
  //
  // parameters:
  //  correlation: response matrix as measured vs. generated
  //  efficiency:  (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied.
//.........这里部分代码省略.........

TH1 *
UnfoldMe(Char_t *data, Char_t *mc, Char_t *anatag, Int_t bin, Bool_t useMBcorr = kTRUE, Bool_t usecorrfit = kFALSE, Bool_t ismc = kFALSE, Float_t smooth = 0.001, Int_t iter = 50, Int_t regul = AliUnfolding::kPowerLaw, Float_t weight = 100., Bool_t bayesian = kTRUE, Int_t nloop = 1)
{

  if (ismc)
    TFile *fdt = TFile::Open(data);
  else
    TFile *fdt = TFile::Open(data);
  TFile *fmc = TFile::Open(mc);
  
  TList *ldt = (TList *)fdt->Get(Form("clist_%s", anatag));
  TList *lmc = (TList *)fmc->Get(Form("clist_%s", anatag));
  
  TH2 *hmatdt = (TH2 *)ldt->FindObject(Form("b%d_corrMatrix", bin));
  if (useMBcorr)
    TH2 *hmatmc = (TH2 *)lmc->FindObject("effMatrix");
  else
    TH2 *hmatmc = (TH2 *)lmc->FindObject(Form("b%d_corrMatrix", bin));
 
  TH1 *hdata = hmatdt->ProjectionY("hdata");
  hdata->Sumw2();
  hdata->SetBinContent(1, 0.);
  hdata->SetBinError(1, 0.);
  //  hdata->Scale(1. / hdata->Integral());
  hdata->SetMarkerStyle(25);
  TH1 *htrue = hmatdt->ProjectionX("htrue");
  htrue->Sumw2();
  //  htrue->Scale(1. / htrue->Integral());
  htrue->SetMarkerStyle(7);
  htrue->SetMarkerColor(2);
  htrue->SetBinContent(1, 0.);
  htrue->SetBinError(1, 0.);
  TH2 *hcorr = (TH2 *)hmatmc->Clone("hcorr");
  TH1 *hinit = (TH1 *)hdata->Clone("hinit");
  TH1 *hresu = (TH1 *)hdata->Clone("hresu");
  TH1 *hbias = (TH1 *)hdata->Clone("hbias");
  hresu->SetMarkerStyle(20);
  hresu->SetMarkerColor(4);
  hresu->Reset();

  TH1 *hnum = hcorr->ProjectionY("hnum");
  TH1 *hden = hcorr->ProjectionY("hden");
  TH1 *heff = hcorr->ProjectionY("heff");
  hnum->Reset();
  hnum->Sumw2();
  hden->Reset();
  hden->Sumw2();
  heff->Reset();
  for (Int_t i = 0; i < heff->GetNbinsX(); i++) {
    Float_t int1 = hcorr->Integral(i + 1, i + 1, 0, -1);
    if (int1 <= 0.) continue;
    Float_t int2 = hcorr->Integral(i + 1, i + 1, 2, -1);
    hnum->SetBinContent(i + 1, int2);
    hnum->SetBinError(i + 1, TMath::Sqrt(int2));
    hden->SetBinContent(i + 1, int1);
    hden->SetBinError(i + 1, TMath::Sqrt(int1));
  }
  new TCanvas("cEfficiency");
  heff->Divide(hnum, hden, 1., 1., "B");
  heff->Draw();
#if 0
  for (Int_t ii = 0; ii < heff->GetNbinsX(); ii++) {
    heff->SetBinContent(ii + 1, 1.);
    heff->SetBinError(ii + 1, 0.);
  }
#endif
  
  for (Int_t i = 0; i < hcorr->GetNbinsX(); i++) {
    hcorr->SetBinContent(i + 1, 1, 0.);
    hcorr->SetBinError(i + 1, 1, 0.);
  }
  for (Int_t i = 0; i < hcorr->GetNbinsY(); i++) {
    hcorr->SetBinContent(1, i + 1, 0.);
    hcorr->SetBinError(1, i + 1, 0.);
  }
  TH2 *hcorrfit = ReturnCorrFromFit(hcorr);

  for (Int_t iloop = 0; iloop < nloop; iloop++) {
    if (bayesian) {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
      AliUnfolding::SetBayesianParameters(smooth, iter);
    } else {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
      AliUnfolding::SetChi2Regularization(regul, weight);
    }
    AliUnfolding::SetSkip0BinInChi2(kTRUE);
    AliUnfolding::SetSkipBinsBegin(1);
    AliUnfolding::SetNbins(150, 150);
    AliUnfolding::Unfold(usecorrfit ? hcorrfit : hcorr, heff, hdata, hinit, hresu);
    hinit = (TH1 *)hresu->Clone(Form("hinit_%d", iloop));
  }

  printf("hdata->Integral(2, -1) = %f\n", hdata->Integral(2, -1));
  printf("hresu->Integral(2, -1) = %f\n", hresu->Integral(2, -1));
  
  
  TCanvas *cUnfolded = new TCanvas ("cUnfolded", "", 400, 800);
  cUnfolded->Divide(1, 2);
  cUnfolded->cd(1)->SetLogx();
  cUnfolded->cd(1)->SetLogy();
//.........这里部分代码省略.........

void QAtracklets(const Char_t *fdata, const Char_t *fmc)
{

  style();
  
  TFile *fdtin = TFile::Open(fdata);
  TList *ldtin = (TList *)fdtin->Get("clist");
  TH2 *hdtin = (TH2 *)ldtin->FindObject("etaphiTracklets");
  TH1 *pdtin = (TH1 *)hdtin->ProjectionY("pdtin_tracklets");
  pdtin->SetMarkerStyle(20);
  pdtin->SetMarkerSize(1.);
  pdtin->SetMarkerColor(kAzure-3);
  hdtin->Scale(1. / hdtin->GetEntries());
  pdtin->Scale(1. / hdtin->GetEntries());
  
  TFile *fmcin = TFile::Open(fmc);
  TList *lmcin = (TList *)fmcin->Get("clist");

  if(!lmcin) {
      std::cout << "NOLIST" << std::endl;

  }
  lmcin->ls();
  TH2 *hmcin = (TH2 *)lmcin->FindObject("etaphiTracklets");
  if(!hmcin) {
    std::cout << "NO H!! etaphiTracklets" << std::endl;
    
  }
  TH1 *pmcin = (TH1 *)hmcin->ProjectionY("pmcin_tracklets");
  pmcin->SetLineColor(kRed+1);
  pmcin->SetFillStyle(1001);
  pmcin->SetFillColorAlpha(kRed+1, 0.1);
  hmcin->Scale(1. / hmcin->GetEntries());
  pmcin->Scale(1. / hmcin->GetEntries());
  

  /*  
  pdtin->Scale(pmcin->Integral(pmcin->FindBin(0. + 0.001),
			       pmcin->FindBin(1. - 0.001))
	       / pdtin->Integral(pdtin->FindBin(0. + 0.001),
				 pdtin->FindBin(1. - 0.001)));
  */  
  
  TCanvas *cData = new TCanvas("cTrackletData", "cTrackletData", 800, 800);
  //  cData->SetLogz();
  TH1 * hfr = cData->DrawFrame(-2.5, 0., 2.5, 2. * TMath::Pi());
  hfr->SetTitle(";#eta;#varphi");
  hdtin->Draw("same,col");
  cData->SaveAs(canvasPrefix+"trackletData.pdf");

  TCanvas *cMC = new TCanvas("cTrackletMC", "cTrackletMC", 800, 800);
  //  cMC->SetLogz();
  hfr = cMC->DrawFrame(-2.5, 0., 2.5, 2. * TMath::Pi());
  hfr->SetTitle(";#eta;#varphi");
  hmcin->Draw("same,col");
  cMC->SaveAs(canvasPrefix+"trackletMC.pdf");
  
  TCanvas *cPhi = new TCanvas("cTrackletPhi", "cTrackletPhi", 800, 800);
  //  cPhi->SetLogy();
  hfr = cPhi->DrawFrame(0., 0., 2. * TMath::Pi(), 0.01);
  hfr->SetTitle(";#varphi;");
  pdtin->DrawCopy("same");
  pmcin->DrawCopy("same,histo");
  TLegend *legend = new TLegend(0.20, 0.18+0.63, 0.50, 0.30+0.63);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(pdtin, "data", "pl");
  legend->AddEntry(pmcin, "Monte Carlo", "l");
  legend->Draw("same");
  cPhi->SaveAs(canvasPrefix+"trackletPhi.pdf");
  
  TCanvas *cPhir = new TCanvas("cTrackletPhir", "cTrackletPhir", 800, 800);
  //  cPhi->SetLogy();
  hfr = cPhir->DrawFrame(0., 0.5, 2. * TMath::Pi(), 1.5);
  hfr->SetTitle(";#varphi;data / Monte Carlo");
  pdtin->Divide(pmcin);
  pdtin->Draw("same");
  cPhir->SaveAs(canvasPrefix+"trackletPhir.pdf");

  
  
}

void makePlot(TCanvas* canvas, const std::string& outputFileName, TTree* testTree, const std::string& varName, 
	      unsigned numBinsX, double xMin, double xMax)
{
  std::cout << "creating histogramTauIdPassed..." << std::endl;
  TString histogramTauIdPassedName = TString("histogramTauIdPassed").Append("_").Append(varName.data());
  TH1* histogramTauIdPassed = fillHistogram(testTree, varName, "type==1", "",
					    histogramTauIdPassedName.Data(), numBinsX, xMin, xMax);
  std::cout << "--> histogramTauIdPassed = " << histogramTauIdPassed << ":" 
	    << " integral = " << histogramTauIdPassed->Integral() << std::endl;

  std::cout << "creating histogramTauIdFailed..." << std::endl;
  TString histogramTauIdFailedName = TString("histogramTauIdFailed").Append("_").Append(varName.data());
  TH1* histogramTauIdFailed = fillHistogram(testTree, varName, "type==0", "",
					    histogramTauIdFailedName.Data(), numBinsX, xMin, xMax);
  std::cout << "--> histogramTauIdFailed = " << histogramTauIdFailed 
	    << " integral = " << histogramTauIdFailed->Integral() << std::endl;

  std::cout << "creating histogramTauIdDenominator..." << std::endl;
  TString histogramTauIdDenominatorName = TString("histogramTauIdDenominator").Append("_").Append(varName.data());
  TH1* histogramTauIdDenominator = new TH1F(histogramTauIdDenominatorName.Data(), 
					    histogramTauIdDenominatorName.Data(), numBinsX, xMin, xMax);
  histogramTauIdDenominator->Add(histogramTauIdPassed);
  histogramTauIdDenominator->Add(histogramTauIdFailed);
  std::cout << "--> histogramTauIdDenominator = " << histogramTauIdDenominator 
	    << " integral = " << histogramTauIdDenominator->Integral() << std::endl;

  std::cout << "creating histogramFakeRate..." << std::endl;
  TString histogramFakeRateName = TString("histogramFakeRate").Append("_").Append(varName.data());
  TH1* histogramFakeRate = new TH1F(histogramFakeRateName.Data(), 
				    histogramFakeRateName.Data(), numBinsX, xMin, xMax);
  histogramFakeRate->Add(histogramTauIdPassed);
  histogramFakeRate->Divide(histogramTauIdDenominator);
  std::cout << "--> histogramFakeRate = " << histogramFakeRate 
	    << " integral = " << histogramFakeRate->Integral() << std::endl;

  std::cout << "creating histogramFakeRateWeighted..." << std::endl;
  TString histogramFakeRateWeightedName = TString("histogramFakeRateWeighted").Append("_").Append(varName.data());
  TH1* histogramFakeRateWeighted = fillHistogram(testTree, varName, "", "MVA_KNN", 
						 histogramFakeRateWeightedName.Data(), numBinsX, xMin, xMax);
  histogramFakeRateWeighted->Divide(histogramTauIdDenominator);
  std::cout << "--> histogramFakeRateWeighted = " << histogramFakeRateWeighted 
	    << " entries = " << histogramFakeRateWeighted->GetEntries() << ","
	    << " integral = " << histogramFakeRateWeighted->Integral() << std::endl;
  // Scale the weighted fake rate histogram

  histogramFakeRate->SetTitle(varName.data());
  histogramFakeRate->SetStats(false);
  histogramFakeRate->SetMinimum(1.e-4);
  histogramFakeRate->SetMaximum(1.e+1);
  histogramFakeRate->SetLineColor(2);
  histogramFakeRate->SetLineWidth(2);
  histogramFakeRate->SetMarkerStyle(20);
  histogramFakeRate->SetMarkerColor(2);
  histogramFakeRate->SetMarkerSize(1);
  histogramFakeRate->Draw("e1p");

  histogramFakeRateWeighted->SetLineColor(4);
  histogramFakeRateWeighted->SetLineWidth(2);
  histogramFakeRateWeighted->SetMarkerStyle(24);
  histogramFakeRateWeighted->SetMarkerColor(4);
  histogramFakeRateWeighted->SetMarkerSize(1);
  histogramFakeRateWeighted->Draw("e1psame");

  TLegend legend(0.11, 0.73, 0.31, 0.89);
  legend.SetBorderSize(0);
  legend.SetFillColor(0);
  legend.AddEntry(histogramFakeRate, "Tau id. discr.", "p");
  legend.AddEntry(histogramFakeRateWeighted, "Fake-Rate weight", "p");
  legend.Draw();

  canvas->Update();
  canvas->Print(outputFileName.data());
}

//.........这里部分代码省略.........
        //New #Jet
        int newNJet = HT3JetVec.size();
        if(verbose==1)printf("newNJet: %d \n ",newNJet);



      // Acceptance determination 1: Counter for all events
      // with muons at generator level
      hAccAll->Fill( binMap_mht_nj[utils2::findBin_mht_nj(evt->nJets(),evt->mht()).c_str()] );
//      hAccAll->Fill( binMap_mht_nj[utils2::findBin_mht_nj(newNJet,newMHT).c_str()] ); // this doesn't work good

      // Check if generator-level muon is in acceptance
      if( genMuPt > LeptonAcceptance::muonPtMin() && std::abs(genMuEta) < LeptonAcceptance::muonEtaMax() ) {
        if(verbose!=0)printf("Muon is in acceptance \n ");
        // Acceptance determination 2: Counter for only those events
        // with generator-level muons inside acceptance
        // hAccPass->Fill( binMap[utils2::findBin(cntNJetsPt30Eta24,nbtag,HT,template_mht).c_str()] );
        hAccPass->Fill( binMap_mht_nj[utils2::findBin_mht_nj(evt->nJets(),evt->mht()).c_str()] );
//        hAccPass->Fill( binMap_mht_nj[utils2::findBin_mht_nj(newNJet,newMHT).c_str()] );

        // Reconstruction-efficiency determination 1: Counter for all events
        // with generator-level muons inside acceptance, regardless of whether
        // the muon has also been reconstructed or not.
        // hIsoRecoAll->Fill( binMap[utils2::findBin(cntNJetsPt30Eta24,nbtag,HT,template_mht).c_str()]);
        hIsoRecoAll->Fill( binMap[utils2::findBin_NoB(evt->nJets(),evt->ht(),evt->mht()).c_str()]);

        // Check if the muon has been reconstructed: check if a reconstructed
        // muon is present in the event that matches the generator-level muon
        // Isolation-efficiency determination 1: Counter for all events with a
        // reconstructed muon that has a generator-level muon match inside the
        // the acceptance, regardless of whether the reconstructed muon is also
        // isolated or not.

//if( evt->MuPtVec_().size()>0 )printf(" RecoMu--> Pt: %g eta: %g phi: %g deltaRMax: %g ",evt->MuPtVec_()[0],evt->MuEtaVec_()[0],evt->MuPhiVec_()[0],deltaRMax); // Ahmad3
//else cout << " Muon size is 0 \n " ;
        // in R and in pt
        int matchedMuonIdx = -1;
        if(evt->MuPtVec_().size()>0 && utils->findMatchedObject(matchedMuonIdx,genMuEta,genMuPhi,evt->MuPtVec_(),evt->MuEtaVec_(),evt->MuPhiVec_(),deltaRMax,verbose) ) {
          // Muon is reconstructed
          const double relDeltaPtMu = std::abs(genMuPt - evt->MuPtVec_().at(matchedMuonIdx) ) / evt->MuPtVec_().at(matchedMuonIdx) ;
          if(verbose!=0)printf(" relDeltaPtMu: %g \n ",relDeltaPtMu);
          if( relDeltaPtMu < deltaPtMax ) {
            // and matches generated pt
            if(verbose!=0)printf("Muon is reconstructed \n ");
            // Check if the muon is also isolated: check if an isolated muon is present
            // in the event that matches the reconstructed muon in R
            if( /*muonsRelIso->at(matchedMuonIdx) <= Selection::muIso()*/ true ){
                //.................//.................//
                // Currently muons are picked if they are isolated.
                // So we don't need to put a cut here.     
                //.................//.................//
              // Muon is isolated
              if(verbose!=0)printf("Muon is isolated \n ");
              // Reconstruction-efficiency determination 2: Counter for those events
              // with generator-level muons inside acceptance where the muon has also
              // been reconstructed.
              // Isolation-efficiency determination 2: Counter for those events where
              // the muon is also isolated.
              // hIsoRecoPass->Fill( binMap[utils2::findBin(cntNJetsPt30Eta24,nbtag,HT,template_mht).c_str()] );
              hIsoRecoPass->Fill( binMap[utils2::findBin_NoB(evt->nJets(),evt->ht(),evt->mht()).c_str()] );
            } // End of muon is isolated
          } // End of pt matching
        } // End of reconstructed muon
      } // End of muon in acceptance



    } // end of loop over events


     // Compute acceptance
    TH1* hAcc = static_cast<TH1*>(hAccPass->Clone("hAcc"));
    hAcc->Divide(hAccPass,hAccAll,1,1,"B");// we use B option here because the two histograms are correlated. see TH1 page in the root manual.

    // Compute efficiencies
    TH1* hEff = static_cast<TH1*>(hIsoRecoPass->Clone("hEff"));
    hEff->Divide(hIsoRecoPass,hIsoRecoAll,1,1,"B");

    if(verbose!=0){
      for(int j=1; j<= totNbins; j++){
        printf("hAccAll: %g hAccPass: %g hAcc: %g hIsoRecoAll: %g hIsoRecoPass: %g hEff: %g \n ",hAccAll->GetBinContent(j),hAccPass->GetBinContent(j),hAcc->GetBinContent(j),hIsoRecoAll->GetBinContent(j),hIsoRecoPass->GetBinContent(j),hEff->GetBinContent(j));
      }
    }


    // --- Save the Histograms to File -----------------------------------
    sprintf(tempname,"LostLepton/LostLepton2_MuonEfficienciesFrom%s_%s.root",subSampleKey.c_str(),inputnumber.c_str());
    TFile outFile(tempname,"RECREATE");
    hAcc->Write();
    hEff->Write();
    hAccAll->Write();
    hAccPass->Write();
    hIsoRecoAll->Write();
    hIsoRecoPass->Write();
    outFile.Close();




  } // end of main

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

  THStack* histogramStack_density = new THStack("stack", "");
  histogramStack_density->Add(histogramRares_density);
  histogramStack_density->Add(histogramEWK_density);
  histogramStack_density->Add(histogramTTV_density);
  histogramStack_density->Add(histogramTT_density);
  histogramStack_density->Draw("histsame");
  
  histogramBgrUncertainty_density->SetFillColor(kBlack);
  histogramBgrUncertainty_density->SetFillStyle(3344);    
  histogramBgrUncertainty_density->Draw("e2same");
  legend->AddEntry(histogramBgrUncertainty_density, "Uncertainty", "f");

  histogramTTH_density->SetLineWidth(2);
  histogramTTH_density->SetLineStyle(1);
  histogramTTH_density->SetLineColor(kBlue);
  histogramTTH_density->Draw("histsame");
  
  histogramData_density->Draw("epsame");
  histogramData_density->Draw("axissame");
  
  legend->Draw();
  
  addLabel_CMS_luminosity(0.2050, 0.9225, 0.6850);
  
  canvas->cd();
  bottomPad->Draw();
  bottomPad->cd();
  
  TH1* histogramRatio = (TH1*)histogramData->Clone("histogramRatio");
  histogramRatio->Reset();
  if ( !histogramRatio->GetSumw2N() ) histogramRatio->Sumw2();
  checkCompatibleBinning(histogramRatio, histogramBgrSum);
  histogramRatio->Divide(histogramData, histogramBgrSum);
  int numBins_bottom = histogramRatio->GetNbinsX();
  for ( int iBin = 1; iBin <= numBins_bottom; ++iBin ) {
    double binContent = histogramRatio->GetBinContent(iBin);
    if ( histogramData && histogramData->GetBinContent(iBin) >= 0. ) histogramRatio->SetBinContent(iBin, binContent - 1.0);
    else histogramRatio->SetBinContent(iBin, -10.);
  }
  histogramRatio->SetTitle("");
  histogramRatio->SetStats(false);
  histogramRatio->SetMinimum(-0.50);
  histogramRatio->SetMaximum(+0.50);
  histogramRatio->SetMarkerStyle(histogramData_density->GetMarkerStyle());
  histogramRatio->SetMarkerSize(histogramData_density->GetMarkerSize());
  histogramRatio->SetMarkerColor(histogramData_density->GetMarkerColor());
  histogramRatio->SetLineColor(histogramData_density->GetLineColor());
  histogramRatio->Draw("ep");
  
  TAxis* xAxis_bottom = histogramRatio->GetXaxis();
  xAxis_bottom->SetTitle(xAxis_top->GetTitle());
  xAxis_bottom->SetLabelColor(1);
  xAxis_bottom->SetTitleColor(1);
  xAxis_bottom->SetTitleOffset(1.20);
  xAxis_bottom->SetTitleSize(0.13);
  xAxis_bottom->SetLabelOffset(0.02);
  xAxis_bottom->SetLabelSize(0.10);
  xAxis_bottom->SetTickLength(0.055);
  
  TAxis* yAxis_bottom = histogramRatio->GetYaxis();
  yAxis_bottom->SetTitle("#frac{Data - Simulation}{Simulation}");
  yAxis_bottom->SetTitleOffset(0.80);
  yAxis_bottom->SetNdivisions(505);
  yAxis_bottom->CenterTitle();
  yAxis_bottom->SetTitleSize(0.09);

//.........这里部分代码省略.........
  // Interface to the event content
  Event* evt = new Event(Sample::fileNameFullSample(id),nEvts);

  // Loop over the events (tree entries)
  while( evt->loadNext() ) {

    // Select only events where the W decayed into a hadronically
    // decaying tau
    if( !(evt->flgW() == 15 && evt->flgTau() == 1) ) continue;

    // Kinematic variables of generator-level tau
    const float genTauPt = evt->genLeptonPt();
    const float genTauEta = evt->genLeptonEta();
    const float genTauPhi = evt->genLeptonPhi();

    // Use only events where the tau is inside the muon acceptance
    // because lateron we will apply the response to muon+jet events
    if( genTauPt < TauResponse::ptMin() ) continue;
    if( std::abs(genTauEta) > TauResponse::etaMax() ) continue;


    // "Cross cleaning": find the jet that originates in the
    // hadronic tau-decay. Associate the jet that is closest in
    // deltaR to the tau. The distance has to be smaller than deltaRMax.

    // First, fill counters before jet-tau matching
    hNJetAll->Fill(evt->nJets());
    hHtAll->Fill(evt->ht());
    hMhtAll->Fill(evt->mht());

    // Do the matching
    int tauJetIdx = -1;		// Will store the index of the jet matched to the tau
    const float deltaRMax = genTauPt < 50. ? 0.2 : 0.1; // Increase deltaRMax at low pt to maintain high-enought matching efficiency
    if( !utils::findMatchedObject(tauJetIdx,genTauEta,genTauPhi,evt->jetsEta(),evt->jetsPhi(),evt->jetsN(),deltaRMax) ) continue;

    // Then, fill counters after matching
    hNJetPass->Fill(evt->nJets());
    hHtPass->Fill(evt->ht());
    hMhtPass->Fill(evt->mht());


    // Calculate RA2 selection-variables from "cleaned" jets,
    // i.e. jets withouth the tau-jet
    int selNJet = 0; // Number of HT jets (jets pt > 50 GeV and |eta| < 2.5)

    for(int jetIdx = 0; jetIdx < evt->jetsN(); ++jetIdx) {	// Loop over reco jets
      // Skip this jet if it is the tau
      if( jetIdx == tauJetIdx ) continue;

      // Calculate NJet 
      if( evt->jetsPt()[jetIdx] > Selection::htJetPtMin() && std::abs(evt->jetsEta()[jetIdx]) < Selection::htJetEtaMax() ) selNJet++;
    } // End of loop over reco jets
    
    
    // Select only events with at least 2 HT jets
    if( selNJet < 2 ) continue;
    

    // Fill histogram with relative visible energy of the tau
    // ("tau response template") for hadronically decaying taus
    for(int jetIdx = 0; jetIdx < evt->jetsN(); ++jetIdx) {	// Loop over reco jets
      // Select tau jet
      if( jetIdx == tauJetIdx ) {
	// Get the response pt bin for the tau
	const unsigned int ptBin = TauResponse::ptBin(genTauPt);
	// Fill the corresponding response template
	const double tauJetPt = evt->jetsPt()[jetIdx];
	hTauResp.at(ptBin)->Fill( tauJetPt / genTauPt );
	break;		// End the jet loop once the tau jet has been found
      }
    }	// End of loop over reco jets
  } // End of loop over tree entries


  // Normalize the response distributions to get the probability density
  for(unsigned int i = 0; i < hTauResp.size(); ++i) {
    if( hTauResp.at(i)->Integral("width") > 0. ) {
      hTauResp.at(i)->Scale(1./hTauResp.at(i)->Integral("width"));
    }
  }

  // Get the jet-tau matching efficiency
  TH1* hMatchEffNJet = static_cast<TH1*>(hNJetPass->Clone("hMatchEffNJet"));
  hMatchEffNJet->Divide(hNJetAll);
  TH1* hMatchEffHt = static_cast<TH1*>(hHtPass->Clone("hMatchEffHt"));
  hMatchEffHt->Divide(hHtAll);
  TH1* hMatchEffMht = static_cast<TH1*>(hMhtPass->Clone("hMatchEffMht"));
  hMatchEffMht->Divide(hMhtAll);


  // --- Save the Histograms to File -----------------------------------
  TFile outFile("../data/HadTau_TauResponseTemplates.root","RECREATE");
  for(unsigned int i = 0; i < hTauResp.size(); ++i) {
    hTauResp.at(i)->Write();
  }
  hMatchEffNJet->Write();
  hMatchEffHt->Write();
  hMatchEffMht->Write();
  outFile.Close();
}

void makeMuEff(const string filename="") 
{
	if (filename.length()<1)
	{
		cout << "Need input root file!" << endl;
		return;
	}

	vector< pair<double,double> >  vJetBins, vHtBins;
	vJetBins.push_back(make_pair(2,2));
	vJetBins.push_back(make_pair(3,5));
	vJetBins.push_back(make_pair(6,7));
	vJetBins.push_back(make_pair(8,1000));
	vHtBins.push_back(make_pair(500,8000));

	TFile f(filename.c_str());
	if (f.IsZombie())
	{
		cout << "File with name " << filename << " cannot be opened or does not exist!" << endl;
		return;
	}

	const string numname("muEff_num");
	const string denname("muEff_den");

	TCanvas *c = new TCanvas();
	size_t dotpos = filename.find(".");
	string substr = filename.substr(0,dotpos);
	stringstream epsname_o, epsname, epsname_c;
	epsname_o << substr << ".eps[";
	epsname << substr << ".eps";
	epsname_c << substr << ".eps]";
	c->Print(epsname_o.str().c_str());
	gStyle->SetOptStat(0);
	gStyle->SetPaintTextFormat("1.1f");
	for (unsigned jetbin=0; jetbin < vJetBins.size(); ++jetbin)
	{
		double minjet= vJetBins.at(jetbin).first;
		double maxjet= vJetBins.at(jetbin).second;
		for (unsigned htbin=0; htbin < vHtBins.size(); ++htbin)
		{
			double minht= vHtBins.at(htbin).first;
			double maxht= vHtBins.at(htbin).second;
			stringstream folder, numhistname, denhistname;
			folder << "muoneff/Hist/Njet" << minjet << "to" << maxjet << "Ht" << minht << "to" << maxht << "/";
			numhistname << folder.str() << numname;
			denhistname << folder.str() << denname;
			TH1* numhist = dynamic_cast<TH1*>( f.Get(numhistname.str().c_str()));
			TH1* denhist = dynamic_cast<TH1*> (f.Get(denhistname.str().c_str()));

			if (numhist == NULL) 
			{
				cout << "Numerator hist " << numhistname.str() << " is not found in file " << f.GetName() << endl;
				return;
			}

			if (denhist == NULL) 
			{
				cout << "Denominator hist " << denhistname.str() << " is not found in file " << f.GetName() << endl;
				return;
			}
			numhist->Draw("colz");
			c->Print(epsname.str().c_str());
			denhist->Draw("colz");
			c->Print(epsname.str().c_str());
			numhist->Sumw2();
			numhist->Divide(denhist);
			stringstream title;
			title << "Jets " << minjet << "-" << maxjet << ", HT>" << minht <<  ": Muon Reco+ID efficiency from " << substr << endl;
			numhist->SetTitle(title.str().c_str());
			//numhist->GetXaxis()->SetRangeUser(0,500);
			numhist->SetAxisRange(0,500,"X");
			numhist->SetAxisRange(0,3.5,"Y");
			//numhist->GetYaxis()->SetRangeUser(0,3.5);
			numhist->Draw("colzTEXT90E");
			c->Print(epsname.str().c_str());
		}
	}
	c->Print(epsname_c.str().c_str());

}