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

TH1* subtractHistograms(const std::string& newHistogramName, const TH1* histogram_central, const TH1* histogram_shift, int mode)
{
  const TH1* histogramMinuend = 0;
  const TH1* histogramSubtrahend = 0;
  if ( compIntegral(histogram_central) >= compIntegral(histogram_shift) ) {
    histogramMinuend    = histogram_central;
    histogramSubtrahend = histogram_shift;
  } else {
    histogramMinuend    = histogram_shift;
    histogramSubtrahend = histogram_central;
  }
  if ( !(histogramMinuend && histogramSubtrahend) ) return 0;
  TH1* newHistogram = (TH1*)histogramMinuend->Clone(newHistogramName.data());
  newHistogram->Reset();
  if ( !newHistogram->GetSumw2N() ) newHistogram->Sumw2();
  int numBins = newHistogram->GetNbinsX();
  for ( int iBin = 0; iBin <= (numBins + 1); ++iBin ) {
    double newBinContent = histogramMinuend->GetBinContent(iBin) - histogramSubtrahend->GetBinContent(iBin);
    double newBinError2 = square(histogramMinuend->GetBinError(iBin)) + square(histogramSubtrahend->GetBinError(iBin));
    if ( mode == kRelative ) {
      if ( histogram_central->GetBinContent(iBin) > 0. ) {
	newBinContent /= histogram_central->GetBinContent(iBin);
	newBinError2 /= square(histogram_central->GetBinContent(iBin));
      } else {
	newBinContent = 0.;
	newBinError2 = 0.;
      }
    }
    newHistogram->SetBinContent(iBin, newBinContent);
    assert(newBinError2 >= 0.);
    newHistogram->SetBinError(iBin, TMath::Sqrt(newBinError2));
  }
  return newHistogram;
}

  //____________________________________________________________________
  void Run(const char* newName,        const char* oldName,
	   const char* newTitle="New", const char* oldTitle="Old")
  {
    TFile* newFile = TFile::Open(newName,"READ");
    TFile* oldFile = TFile::Open(oldName,"READ");
    if (!newFile || !oldFile) return;

    TH1* newCent = GetH1(newFile, "realCent");
    TH1* oldCent = GetH1(oldFile, "realCent");
    if (!newCent || !oldCent) return;

    TString  t; t.Form("#it{R}=#frac{%s}{%s}", newTitle, oldTitle);
    TCanvas* c     = new TCanvas("c", t, 1200, 800);
    c->SetTopMargin(0.01);
    c->SetRightMargin(0.20);
    fLegend = new TLegend(1-c->GetRightMargin(),
			  c->GetBottomMargin(),
			  1, 1-c->GetTopMargin(),
			  t);
    fLegend->SetFillStyle(0);
    fLegend->SetBorderSize(0);
    THStack* stack = new THStack("ratios","");
			       
    fMin = +1e6;
    fMax = -1e6;
    TH1* one = 0;
    for (Int_t i = newCent->GetNbinsX(); i--;) {
      Double_t c1 = newCent->GetXaxis()->GetBinLowEdge(i+1);
      Double_t c2 = newCent->GetXaxis()->GetBinUpEdge(i+1);
      Info("", "c1=%f c2=%f", c1, c2);
      TH1*     r  = One(newFile, oldFile, c1, c2);    
      if (!r) continue;
      if (!one) {
	one = static_cast<TH1*>(r->Clone("one"));
	one->SetDirectory(0);
	one->Reset();
	for (Int_t j = 1; j <= one->GetNbinsX(); j++) {
	  one->SetBinContent(j,1);
	  one->SetBinError  (j,0);
	}
      }
      // r->Add(one, i-1);
      // r->Scale(TMath::Power(10,i));
      stack->Add(r);
    }
    stack->Draw("nostack");
    stack->SetMinimum(0.95*fMin);
    stack->SetMaximum(1.05*fMax);
    stack->GetHistogram()->SetXTitle("#eta");
    stack->GetHistogram()->SetYTitle("#it{R}");
    fLegend->Draw();
    c->Modified();
    c->Update();
    c->cd();
    c->SaveAs(Form("%sover%s.png", newTitle, oldTitle));
  }  

TH1* makeCDF(TH1* h) {
	TString sName(TString(h->GetName())+TString("_CDF"));
	TString sTitle(TString(h->GetTitle())+TString(" CDF"));
	TH1* hOut = (TH1*) h->Clone(sName);
	hOut->SetTitle(sTitle);
	hOut->Reset();
	double cdf = 0;
	for (int ibin=0; ibin < h->GetNbinsX()+2; ++ibin) {
		cdf += h->GetBinContent(ibin);
		hOut->SetBinContent(ibin,cdf);
	}
	hOut->Scale(1.0/(h->Integral(0,h->GetNbinsX()+1)));
	return hOut;
}

TH1 *
YieldMean_ReturnRandom(TH1 *hin)
{
  TH1 *hout = (TH1 *)hin->Clone("hout");
  hout->Reset();
  Double_t cont, err;
  for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
    if (hin->GetBinError(ibin + 1) <= 0.) continue;
    cont = hin->GetBinContent(ibin + 1);
    err = hin->GetBinError(ibin + 1);
    hout->SetBinContent(ibin + 1, gRandom->Gaus(cont, err));
    hout->SetBinError(ibin + 1, err);
  }
  return hout;
}

void drawTimeDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstTime");
  TH1I* hLast = (TH1I*)directory->Get("lastTime");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaT");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaT");
  h->Reset();
  h->SetTitle("DeltaT");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaT");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    std::time_t t1 = hFirst->GetAt(i);
    std::time_t t2 = hLast->GetAt(i);
    TTimeStamp ts1(hFirst->GetAt(i));
    std::cout << "Fit started at " << ts1.AsString() << std::endl;
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,difftime(t2,t1));
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

void drawEventDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstEvent");
  TH1I* hLast = (TH1I*)directory->Get("lastEvent");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaE");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaE");
  h->Reset();
  h->SetTitle("DeltaE");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaE");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    int e1 = hFirst->GetAt(i);
    int e2 = hLast->GetAt(i);
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,e2-e1);
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

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 CalcAcceptanceSystematicsFromScale(){

  char tempname[200];
  TFile * infile = new TFile("LostLepton2_MuonEfficienciesFromAll_Elog408_.root","R");
  TH1* hScaleAccPass, * hScaleAccAll, * hScaleAccPass_lowDphi, * hScaleAccAll_lowDphi;
  vector<TH1*> hScaleAccVec, hScaleAcc_lowDphiVec, hScaleAcc_DeviationFromNomVec, hScaleAcc_DeviationFromNom_lowDphiVec;

  int Scalesize = 9; // 101 pdf
  for(int iacc=0; iacc < Scalesize; iacc++){

    // get the pass
    sprintf(tempname,"ScaleAcc/hScaleAccPass_%d",iacc);
    hScaleAccPass = (TH1*) infile->Get(tempname)->Clone();

    // get the all
    sprintf(tempname,"ScaleAcc/hScaleAccAll_%d",iacc);
    hScaleAccAll = (TH1*) infile->Get(tempname)->Clone();    
    
    // divide the pass and all to get acc
    sprintf(tempname,"hScaleAccVec_%d",iacc);
    hScaleAccVec.push_back(static_cast<TH1*>(hScaleAccPass->Clone(tempname)));
    hScaleAccVec[iacc]->Divide(hScaleAccPass,hScaleAccAll,1,1,"B");

    cout << " bin5=> acc: " << hScaleAccVec[iacc]->GetBinContent(5) << endl;

    // lowDphi
    sprintf(tempname,"ScaleAcc/hScaleAccPass_lowDphi_%d",iacc);
    hScaleAccPass_lowDphi = (TH1*) infile->Get(tempname)->Clone();

    sprintf(tempname,"ScaleAcc/hScaleAccAll_lowDphi_%d",iacc);
    hScaleAccAll_lowDphi = (TH1*) infile->Get(tempname)->Clone();

    sprintf(tempname,"hScaleAcc_lowDphiVec_%d",iacc);
    hScaleAcc_lowDphiVec.push_back(static_cast<TH1*>(hScaleAccPass_lowDphi->Clone(tempname)));
    hScaleAcc_lowDphiVec[iacc]->Divide(hScaleAccPass_lowDphi,hScaleAccAll_lowDphi,1,1,"B");


    // calculate the deviation from nominal acceptance
    hScaleAcc_DeviationFromNomVec.push_back(static_cast<TH1*>(hScaleAccVec[iacc]->Clone("hScaleAcc_DeviationFromNomVec")));// copy
    hScaleAcc_DeviationFromNomVec[iacc]->Add(hScaleAccVec[0],-1.0); // subtract the nominal from each acceptance: Acc - Acc_nom

    // do the same for lowDphi
    hScaleAcc_DeviationFromNom_lowDphiVec.push_back(static_cast<TH1*>(hScaleAcc_lowDphiVec[iacc]->Clone("hScaleAcc_DeviationFromNom_lowDphiVec")));// copy
    hScaleAcc_DeviationFromNom_lowDphiVec[iacc]->Add(hScaleAcc_lowDphiVec[0],-1.0); // subtract the nominal from each acceptance: Acc - Acc_nom

  }// end of loop over pdf



  TH1 * hScaleAccSysMax = static_cast<TH1*>(hScaleAccAll->Clone("hScaleAccSysMax"));
  hScaleAccSysMax->Reset();
  TH1 * hScaleAccSysMax_lowDphi = static_cast<TH1*>(hScaleAccAll_lowDphi->Clone("hScaleAccSysMax_lowDphi"));
  hScaleAccSysMax_lowDphi->Reset();
  TH1 * hScaleAccSysMin = static_cast<TH1*>(hScaleAccAll->Clone("hScaleAccSysMin"));
  hScaleAccSysMin->Reset();
  TH1 * hScaleAccSysMin_lowDphi = static_cast<TH1*>(hScaleAccAll_lowDphi->Clone("hScaleAccSysMin_lowDphi"));
  hScaleAccSysMin_lowDphi->Reset();

  for(int ibin=0; ibin < hScaleAccSysMax->GetNbinsX()+2; ibin++){
    //cout << "bincontent: " <<  hScaleAccSysMax->GetBinContent(ibin) << endl;
    for(int iacc=0; iacc < Scalesize; iacc++){
      if(hScaleAccSysMax->GetBinContent(ibin)<hScaleAcc_DeviationFromNomVec[iacc]->GetBinContent(ibin)){
        hScaleAccSysMax->SetBinContent(ibin,hScaleAcc_DeviationFromNomVec[iacc]->GetBinContent(ibin)); // get the max of all acc
      }
      if(hScaleAccSysMax_lowDphi->GetBinContent(ibin)<hScaleAcc_DeviationFromNom_lowDphiVec[iacc]->GetBinContent(ibin)){
        hScaleAccSysMax_lowDphi->SetBinContent(ibin,hScaleAcc_DeviationFromNom_lowDphiVec[iacc]->GetBinContent(ibin)); 
      }
      if(hScaleAccSysMin->GetBinContent(ibin)>hScaleAcc_DeviationFromNomVec[iacc]->GetBinContent(ibin)){
        hScaleAccSysMin->SetBinContent(ibin,hScaleAcc_DeviationFromNomVec[iacc]->GetBinContent(ibin)); // get the max of all acc
      }
      if(hScaleAccSysMin_lowDphi->GetBinContent(ibin)>hScaleAcc_DeviationFromNom_lowDphiVec[iacc]->GetBinContent(ibin)){
        hScaleAccSysMin_lowDphi->SetBinContent(ibin,hScaleAcc_DeviationFromNom_lowDphiVec[iacc]->GetBinContent(ibin));
      }
    }
  }


  TFile * outFile = new TFile("AcceptanceSystematicsFromScale_AllSamples.root","RECREATE");
  hScaleAccSysMax->Write();
  hScaleAccSysMax_lowDphi->Write();
  hScaleAccSysMin->Write();
  hScaleAccSysMin_lowDphi->Write();

  for(int iacc=0; iacc < Scalesize ; iacc++){
    hScaleAccVec[iacc]->Write();
    hScaleAcc_lowDphiVec[iacc]->Write();
  }

  outFile->Close();

}

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

  histogramRares_density->SetFillColor(kBlue - 8); 
  legend->AddEntry(histogramRares_density, "Rares", "f");

  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);

void drawBeamSpotGraph (TDirectory* directory, TH1* refHisto, const char* name,
			const char* fname, float* runSummary = 0)
{
  TGraphErrors* graph = directory->Get(name);
  if ( graph==0 )  return;

  gStyle->SetOptFit(1);
  gStyle->SetOptStat(0);

  char newName[64];
  newName[0] = 'c';
  strcpy(&newName[1],name);
  std::string fullName(name);
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  newName[0] = 'h';
  TH1* h = refHisto->Clone(newName);
  h->Reset();
  h->SetTitle(name);
  std::string xTitle("Luminosity block / ");
  xTitle += directory->GetName();
  h->GetXaxis()->SetTitle(xTitle->c_str());
//   int nb = h->GetNbinsX();
//   if ( nb>50 ) {
//     int iscale = nb/50+1;
//     TAxis* xaxis = h->GetXaxis();
//     for ( int i=1; i<=h->GetNbinsX(); ++i ) {
//       if ( (i-1)%iscale )  xaxis->SetBinLabel(i,"");
//     }
//   }
  TAxis* yaxis = h->GetYaxis();
  float scale(1);
  if ( strcmp(name,"x")==0 ) {
    yaxis->SetTitle("PV x position [mm]");
    scale = 10.;
  }
  else if ( strcmp(name,"y")==0 ) {
    yaxis->SetTitle("PV y position [mm]");
    scale = 10.;
  }
  else if ( strcmp(name,"z")==0 ) {
    yaxis->SetTitle("PV z position [cm]");
  }
  else if ( strcmp(name,"ex")==0 ) {
    yaxis->SetTitle("PV x width [#mum]");
    scale = 10000.;
  }
  else if ( strcmp(name,"ey")==0 ) {
    yaxis->SetTitle("PV y width [#mum]");
    scale = 10000.;
  }
  else if ( strcmp(name,"ez")==0 ) {
    yaxis->SetTitle("PV z width [cm]");
  }
  else if ( strcmp(name,"corrxy")==0 ) {
    yaxis->SetTitle("PV x-y correlation");
  }
  else if ( strcmp(name,"dxdz")==0 ) {
    yaxis->SetTitle("PV slope dx/dz [10^{-3}]");
    scale = 1000.;
  }
  else if ( strcmp(name,"dydz")==0 ) {
    yaxis->SetTitle("PV slope dy/dz [10^{-3}]");
    scale = 1000.;
  }
  scaleGraph(graph,scale);

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum((ymax+ymin)/2.-2.*(ymax-ymin)/2.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
  graph->Fit("pol1","same");
  graph->GetFunction("pol1")->SetLineStyle(2);
  graph->GetFunction("pol1")->SetLineWidth(2);

  TF1* fit = graph->GetFunction("pol1");
  cout << "Pol1 fit chi2 = " << fit->GetChisquare() 
       << " " << fit->GetNDF() << endl;

//   string epsName = fullName + ".eps";
//   c->SaveAs(epsName.c_str());
//   string pngName = fullName + ".png";
//   c->SaveAs(pngName.c_str());

  if ( runSummary ) {
    TF1* fit = graph->GetFunction("pol1");
    runSummary[0] = fit->GetChisquare();
    runSummary[1] = fit->GetNDF();
    runSummary[2] = fit->GetParameter(0);
    runSummary[3] = fit->GetParError(0);
    runSummary[4] = fit->GetParameter(1);
    runSummary[5] = fit->GetParError(1);
  }
}

void Fit::nll(int& ndim, double* gout, double& result, double* par, int flags) {
  size_t n = Fit::signals.size();
  result = 0;

  // make summed fit histogram
  TH1* hfit = (TH1*) signals[0].histogram->Clone("hfit");
  hfit->Reset();
  for (size_t i=0; i<n; i++) {
    TH1* h = Fit::signals.at(i).histogram;
    h->Scale(1.0/h->Integral());
    hfit->Add(h, Fit::norms[i] * par[i]);
  }

  // loop over bins
  if (hfit->IsA() == TH2F::Class()) {
    TH2F* h2 = dynamic_cast<TH2F*>(hfit);
    for (int i=1; i<h2->GetNbinsX(); i++) {
      if (i < h2->GetXaxis()->FindBin(Fit::r_range.min) || i > h2->GetXaxis()->FindBin(Fit::r_range.max)) {
        continue;
      }
      for (int j=1; j<h2->GetNbinsY(); j++) {
        if (j < h2->GetYaxis()->FindBin(Fit::e_range.min) || j > h2->GetYaxis()->FindBin(Fit::e_range.max)) {
          continue;
        }
        double nexp = h2->GetBinContent(i, j);
        double nobs = dynamic_cast<TH2F*>(Fit::data)->GetBinContent(i, j);
        result += (nexp - nobs * TMath::Log(TMath::Max(1e-12, nexp)));
        //std::cout << "- " << nexp << " " << nobs << " " << TMath::Log(nexp) << " // " << result << std::endl;
      }
    }
  }
  else {
    for (int i=1; i<hfit->GetNbinsX(); i++) {
      if (i < hfit->FindBin(Fit::e_range.min) || i > hfit->FindBin(Fit::e_range.max)) {
        continue;
      }
      double nexp = hfit->GetBinContent(i);
      double nobs = Fit::data->GetBinContent(i);
      result += (nexp - nobs * TMath::Log(TMath::Max(1e-12, nexp)));
      //std::cout << "- " << nexp << " " << nobs << " " << TMath::Log(nexp) << " // " << result << std::endl;
    }
  }

  // constraints
  for (size_t i=0; i<n; i++) {
    if (Fit::signals.at(i).constraint > 0) {
      result += 0.5 * TMath::Power((par[i]-1.0)
                / Fit::signals.at(i).constraint, 2);
    }
  }
  delete hfit;

#ifdef DEBUG
  // print parameters at each iteration
  std::cout << "+ ";
  for (size_t i=0; i<n; i++) {
    std::cout << par[i] * Fit::signals.at(i).rate << " (" << par[i] << ") \t";
  }
  std::cout << result << std::endl;
#endif
}

void treeClient(Bool_t evol=kFALSE) 
{
   // Client program which creates and fills 2 histograms and a TTree. 
   // Every 1000000 fills the histograms and TTree is send to the server which displays the histogram.
   //
   // To run this demo do the following:
   //   - Open at least 2 windows
   //   - Start ROOT in the first windows
   //   - Execute in the first window: .x fastMergeServer.C
   //   - Execute in the other windows: root.exe -b -l -q .x treeClient.C
   //     (You can put it in the background if wanted).
   // If you want to run the hserv.C on a different host, just change
   // "localhost" in the TSocket ctor below to the desired hostname.
   //
   //Author: Fons Rademakers, Philippe Canal
   
   gBenchmark->Start("treeClient");

   // Open connection to server
   TSocket *sock = new TSocket("localhost", 9090);
   if (!sock->IsValid()) {
      Error("treeClient","Could not establish a connection with the server %s:%d.","localhost",9090);
      return;
   }

   // Wait till we get the start message
   // server tells us who we are
   Int_t status, version, kind;
   sock->Recv(status, kind);
   if (kind != 0 /* kStartConnection */) 
   {
      Error("treeClient","Unexpected server message: kind=%d status=%d\n",kind,status);
      delete sock;
      return;
   }
   sock->Recv(version, kind);
   if (kind != 1 /* kStartConnection */) 
   {
      Fatal("treeClient","Unexpected server message: kind=%d status=%d\n",kind,status);
   } else {
      Info("treeClient","Connected to fastMergeServer version %d\n",version);
   }
   
   int idx = status;
   
   Float_t messlen  = 0;
   Float_t cmesslen = 0;

   TMemFile *file = new TMemFile("mergedClient.root","RECREATE");
   TH1 *hpx;
   if (idx == 0) {
      // Create the histogram
      hpx = new TH1F("hpx","This is the px distribution",100,-4,4);
      hpx->SetFillColor(48);  // set nice fillcolor
   } else {
      hpx = new TH2F("hpxpy","py vs px",40,-4,4,40,-4,4);
   }
   Float_t px, py;
   TTree *tree = new TTree("tree","tree");
   tree->SetAutoFlush(4000000);
   tree->Branch("px",&px);
   tree->Branch("py",&py);
 
   TMessage::EnableSchemaEvolutionForAll(evol);
   TMessage mess(kMESS_OBJECT);

   // Fill histogram randomly
   gRandom->SetSeed();
   const int kUPDATE = 1000000;
   for (int i = 0; i < 25000000; ) {
      gRandom->Rannor(px,py);
      if (idx%2 == 0)
         hpx->Fill(px);
      else
         hpx->Fill(px,py);
      tree->Fill();
      ++i;
      if (i && (i%kUPDATE) == 0) {
         file->Write();
         mess.Reset(kMESS_ANY);              // re-use TMessage object
         mess.WriteInt(idx);
         mess.WriteTString(file->GetName());
         mess.WriteLong64(file->GetEND());   // 'mess << file->GetEND();' is broken in CINT for Long64_t
         file->CopyTo(mess);
         sock->Send(mess);          // send message
         messlen  += mess.Length();
         cmesslen += mess.CompLength();
         
         file->ResetAfterMerge(0);  // This resets only the TTree objects.
         hpx->Reset();
      }
   }
   sock->Send("Finished");          // tell server we are finished

   if (cmesslen > 0)
      printf("Average compression ratio: %g\n", messlen/cmesslen);

   gBenchmark->Show("hclient");

   // Close the socket
//.........这里部分代码省略.........