C++ THStack::GetHists方法代码示例

void DrawEmpirical(const char* filename="Empirical.root", 
		   Bool_t fmd=true)
{
  gStyle->SetOptTitle(0);

  TFile* file = TFile::Open(filename, "READ");
  if (!file) return;

  Double_t yr = 0.3;
  TCanvas* c  = new TCanvas("c","c", 1000,1000);
  TPad*    p1 = new TPad("p1","p1",0,0,1,yr);
  TPad*    p2 = new TPad("p2","p2",0,yr,1,1);
  c->cd(); p1->Draw();
  c->cd(); p2->Draw();
  
  gDirectory->cd("Forward");
  THStack* r = DrawOne(p1, yr, false, gDirectory, "ratios");  
  THStack* e = DrawOne(p2, yr, true, gDirectory, "empirical");

  r->SetMinimum(0.945);
  r->SetMaximum(1.055);
  r->GetXaxis()->SetTitle("#it{#eta}");
  r->GetYaxis()->SetTitle("Ratio to mean");
  e->SetMinimum(0.005);
  e->GetYaxis()->SetTitle("#it{E_{c}}(#it{#eta})");
  TIter nextE(e->GetHists());
  TIter nextR(r->GetHists());
  TH1*  hist = 0;
  Color_t cols[]  = { kRed+2, kGreen+2, kBlue+2, kMagenta+2, 0 };
  Color_t *ptr    = cols;
  Style_t stys[]  = { 20, 21, 22, 23 };
  Style_t* sty    = stys;
  while (*ptr) { 
    hist = static_cast<TH1*>(nextE()); 
    hist->SetMarkerColor(*ptr);
    hist->SetMarkerSize(2);
    hist->SetMarkerStyle(*sty);
    hist = static_cast<TH1*>(nextR()); 
    hist->SetMarkerColor(*ptr);
    hist->SetMarkerSize(2);
    hist->SetMarkerStyle(*sty);
    ptr++;
    sty++;
  }


  TLegend* l = p2->BuildLegend(0.35, .2, .65, .8);
  l->SetFillColor(0);
  l->SetFillStyle(0);
  l->SetBorderSize(0);

  c->Modified();
  c->Update();
  c->cd();
  c->Print("empirical.png");
  
}

void CalcQCDNormFactor() {
  //TFile *f = TFile::Open("results/Plotter_out_2016_05_29_22h19m32.root"); // 76X Silver JSON
  TFile *f = TFile::Open("results/Plotter_out_2016_06_21_15h27m59.root"); // 76X Golden JSON
  TCanvas *c = (TCanvas*)f->Get("NJet/PlotSamples_Pass5Cuts_PassHLT");
  THStack *s = (THStack*)c->GetListOfPrimitives()->At(1);
  TH1D *data = (TH1D*)c->GetListOfPrimitives()->At(3);
  double MC_integral = 0;
  double QCD_count = 0;
  for (int i=s->GetHists()->GetEntries()-1; i>=0; --i) {
    TH1D* h = (TH1D*)s->GetHists()->At(i);
    if (i==s->GetHists()->GetEntries()-1) QCD_count = h->Integral();
    std::cout<<h->GetName()<<" "<<h->Integral()<<std::endl;
    MC_integral += h->Integral();
  }
  double data_QCD_estimate = data->Integral()- (MC_integral-QCD_count);
  double QCD_scale = data_QCD_estimate/QCD_count;
  std::cout<<"MC:                  "<<MC_integral<<std::endl;
  std::cout<<"Data:                "<<data->Integral()<<std::endl;
  std::cout<<"MC   (QCD only):     "<<QCD_count<<std::endl;
  std::cout<<"Data (QCD only est): "<<data_QCD_estimate<<std::endl;
  std::cout<<"QCD Scale: "<<QCD_scale<<std::endl;

  TH1D* qcd = (TH1D*)s->GetHists()->At(s->GetHists()->GetEntries()-1);
  qcd->Scale(QCD_scale);
  c->Draw();
  
  //TCanvas *c = (TCanvas*)f->Get("NJet/PlotSamples_Pass5Cuts_PassHLT_Ratio");
  //
  //TH1D* ratio = (TH1D*)((TVirtualPad*)c->cd(2))->GetListOfPrimitives()->At(0);
  //TF1* fit = new TF1("fit","pol1", 400, 2000);
  //ratio->Fit("fit","RBQ0");
  //fit->SetLineColor(2);
  //fit->SetLineWidth(1);
  //TF1* fit_up   = (TF1*)fit->Clone("fit_up");
  //TF1* fit_down = (TF1*)fit->Clone("fit_down");
  //fit_up  ->SetParameter(0,fit->GetParameter(0)+fit->GetParError(0)*1);
  //fit_down->SetParameter(0,fit->GetParameter(0)-fit->GetParError(0)*1);
  //fit_up  ->SetParameter(1,fit->GetParameter(1)+fit->GetParError(1)*1);
  //fit_down->SetParameter(1,fit->GetParameter(1)-fit->GetParError(1)*1);
  //fit_up  ->SetLineColor(4); fit_up  ->Draw("SAME");
  //fit_down->SetLineColor(4); fit_down->Draw("SAME");
  //fit->Draw("SAME");
  //
  //std::cout<<"Fit result:"<<std::endl;
  //std::cout<<"p0: "<<fit->GetParameter(0)<<" +- "<<fit->GetParError(0)*1<<std::endl;
  //std::cout<<"p1: "<<fit->GetParameter(1)<<" +- "<<fit->GetParError(1)*1<<std::endl;
  //f->Close();  
}

TH1* GetOne(UShort_t sNN, const TString& trigger)
{
  Long_t   p = gROOT->ProcessLine(Form("Drawer::GetStack(0, \"pp\", %d, "
				       "\"%s\", false, true)", 
				       sNN, trigger.Data()));
  THStack* s = (THStack*)p;
  TList*   l = s->GetHists();
  TH1*     h = 0;
  TIter    n(l);
  l->ls();
  while ((h = static_cast<TH1*>(n()))) {
    TString m(h->GetName());
    if (m.EqualTo("dndetaForward_all")) break;
  }

  if (h) {
    switch (sNN) { 
    case  900: h->SetTitle("900GeV");  h->SetMarkerColor(kRed+2);   break;
    case 2760: h->SetTitle("2.76TeV"); h->SetMarkerColor(kGreen+2); break;
    case 7000: h->SetTitle("7TeV");    h->SetMarkerColor(kBlue+2);  break;
    case 8000: h->SetTitle("8TeV");    h->SetMarkerColor(kBlack);   break;
    }
  }

  return h;
}

void fit() {
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      fitmmp(h);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

TH1F* GetHisto(TFile* fin, string region, string process, string varname, float& norm, bool do_norm, float input_norm)
{
    string cname = CHANNEL_NAME+string("/")+region+"/"+varname;
    TCanvas* c = (TCanvas*) fin->Get(cname.c_str());
    string hname = "v:"+varname+"|p:"+process+"|r:"+region+string("|c:")+CHANNEL_NAME+string("|t:1DEntries");
    TH1F* h = 0;
    if(VERBOSE>0){
 	 cerr<<"cname :"<<cname<<endl;
   	 cerr<<"histo name: "<<hname<<endl;
   	 cerr<<"pointer: "<<c<<endl;
    } 
    TList* l = c->GetListOfPrimitives();
    TPad* pad = (TPad*) l->At(0);
    THStack* stack = (THStack*) pad->GetPrimitive("");
    h = (TH1F*) stack->GetHists()->FindObject(hname.c_str());
    if(do_norm) h->Scale(input_norm/h->Integral());
    norm = h->Integral();
    return (TH1F*) h->Clone();
}

//.........这里部分代码省略.........
  for (int iSubDet=1; iSubDet<=6; ++iSubDet) {

    // TEC requires special care since rings 1-4 and 5-7 are plotted separately
    bool isTEC = (iSubDet==6);

    // if subdet is specified, skip other subdets
    if (plotSubDetN!=0 && iSubDet!=plotSubDetN)
      continue;

    // Skips plotting too high layers
    if (plotLayerN > numberOfLayers[iSubDet-1]) {
      continue;
    }

    int minlayer = plotLayers ? 1 : plotLayerN;
    int maxlayer = plotLayers ? numberOfLayers[iSubDet-1] : plotLayerN;
    
    for (int layer = minlayer; layer <= maxlayer; layer++) {

      // two plots for TEC, skip first 
      for (int iTEC = 0; iTEC<2; iTEC++) {
	if (!isTEC && iTEC==0) continue;
	
	char  selection[1000];
	if (!isTEC){
	  if (layer==0)
	    sprintf(selection,"subDetId==%d",iSubDet); 
	  else
	    sprintf(selection,"subDetId==%d && layer == %d",iSubDet,layer); 
	}
	else{	  // TEC
	  if (iTEC==0)  // rings 
	    sprintf(selection,"subDetId==%d && ring <= 4",iSubDet); 
	  else
	    sprintf(selection,"subDetId==%d && ring > 4",iSubDet); 
	}


	// Title for plot and name for the file

	TString subDetName;
	switch (iSubDet) {
	case 1: subDetName = "BPIX"; break;
	case 2: subDetName = "FPIX"; break;
	case 3: subDetName = "TIB"; break;
	case 4: subDetName = "TID"; break;
	case 5: subDetName = "TOB"; break;
	case 6: subDetName = "TEC"; break;
	}

	TString myTitle = "Surface Shape, ";
	myTitle += subDetName;
	if (layer!=0) {
	  myTitle += TString(", layer ");
	  myTitle += Form("%d",layer); 
	}
	if (isTEC && iTEC==0)
	  myTitle += TString(" R1-4");
	if (isTEC && iTEC>0)
	  myTitle += TString(" R5-7");

	// Save plot to file
	std::ostringstream plotName;
	plotName << outputDir << "/SurfaceShape_" << subDetName << "_";
	plotName << residType; 
	if (layer!=0)
	  plotName << "_" << "layer" << layer;
	if (isTEC && iTEC==0)
	  plotName << "_" << "R1-4";
	if (isTEC && iTEC>0)
	  plotName << "_" << "R5-7";
	plotName << ".eps";

	// Generate histograms with selection
	THStack *hs = addHists(selection, residType);
	if (!hs || hs->GetHists()==0 || hs->GetHists()->GetSize()==0) {
	  std::cout << "No histogram for " << subDetName << ", perhaps not enough data?" << std::endl; 
	  continue; 
	}
	hs->SetTitle( myTitle ); 
	hs->Draw("nostack PE");  

	// Adjust Labels
	TH1* firstHisto = (TH1*) hs->GetHists()->First();
	TString xName = firstHisto->GetXaxis()->GetTitle();
	TString yName = firstHisto->GetYaxis()->GetTitle();
	hs->GetHistogram()->GetXaxis()->SetTitleColor( kBlack ); 
	hs->GetHistogram()->GetXaxis()->SetTitle( xName ); 
	hs->GetHistogram()->GetYaxis()->SetTitleColor( kBlack ); 
	hs->GetHistogram()->GetYaxis()->SetTitle( yName ); 

	// Save to file
	c.Update(); 
	c.Print(plotName.str().c_str());
      }
    }
  }

  return;
}

bool check(int n = 2) {
   TFile *file = TFile::Open(TString::Format("merged%d.root",n));

   bool result = true;
   TH1F *h; file->GetObject("h1",h);
   if (!h) {
      Error("execFileMerger","h1 is missing\n");
      result = false;
   }
   if (h->GetBinContent(2) != n || h->GetBinContent(3) != n) {
      Error("execFileMerger","h1 not added properly");
      result = false;
   }
   
   THnSparseF *sparse; file->GetObject("sparse",sparse);
   if (!sparse) {
      Error("execFileMerger","sparse is missing\n");
      result = false;
   } else {
      Int_t coordIdx[5] = {1, 2, 3, 4, 5};
      Double_t cont = sparse->GetBinContent(coordIdx);
      if (cont > n + 0.4 || cont < n - 0.4) {
         Error("execFileMerger","sparse merge failed: expected bin content %g, read %g\n",
               (Double_t)n, cont);
         result = false;
      }
      Double_t entries = sparse->GetEntries();
      if (entries > n + 0.4 || entries < n - 0.4) {
         Error("execFileMerger","sparse merge failed: expected %g entries, read %g\n",
               (Double_t)n, entries);
         result = false;
      }
   }
   
   THStack *stack; file->GetObject("stack",stack);
   if (!stack) {
      Error("execFileMerger","stack is missing\n");
      result = false;
   }
   h = (TH1F*)stack->GetHists()->FindObject("hs_1");
   if (!h) {
      Error("execFileMerger","hs_1 is missing\n");
      result = false;
   }
   if (h->GetBinContent(2) != n || h->GetBinContent(3) != n) {
      Error("execFileMerger","hs_1 not added properly");
      result = false;
   }
   h = (TH1F*)stack->GetHists()->FindObject("hs_2");
   if (!h) {
      Error("execFileMerger","hs_2 is missing\n");
      result = false;
   }
   if (h->GetBinContent(4) != n || h->GetBinContent(5) != n) {
      Error("execFileMerger","hs_2 not added properly");
      result = false;
   }

   TGraph *gr; file->GetObject("exgraph",gr);
   if (!gr) {
      Error("execFileMerger","exgraph is missing\n");
      result = false;
   }
   if (gr->GetN() != ( n * 3)) {
      Error("execFileMerger","exgraph not added properly n=%d rather than %d",gr->GetN(),n*3);
      result = false;            
   } else {
      for(Int_t k = 0; k < gr->GetN(); ++k) {
         double x,y;
         gr->GetPoint(k,x,y);
         if ( x != ( (k%3)+1 ) ||  y != ( (k%3)+1 ) ) {
            Error("execFileMerger","exgraph not added properly");
            result = false;            
         }
      }
   }
   
   TTree *tree; file->GetObject("tree",tree);
   if (!tree) {
      Error("execFileMerger","tree is missing\n");
      result = false;
   }
   if (tree->GetEntries() != n*2) {
      Error("execFileMerger","tree does not have the expected number of entries: %lld rather than %d",tree->GetEntries(),n*2);
      result = false;            
   } else {
      if ( tree->GetEntries("data==1") != n ) {
         Error("execFileMerger","tree does not have the expected data. We got %lld entries with 'data==1' rather than %d",tree->GetEntries("data==1"),n);
         tree->Scan();
         result = false;
      }
   }   
   return result;
}

//.........这里部分代码省略.........
	  if(density) area = Xrange*Yrange; //normalize each bin by its area

	  if(bintotal[binN+1]>0)
	    histUnrolledPercentage[i]->SetBinContent(binN+1, (value/area)/bintotal[binN+1]);

	  binN++;
	}
		
    if ( histUnrolled[i]->Integral() > 0) {
	stackUnrolledPercentage->Add(histUnrolledPercentage[i]);
    }

    cout << "Unrolling Percentage for Process : " << processLabels[i] << "\n";	  
  }

  /// Unrolled plots in bins of R&MR
  TLatex t1(0.1,0.92, "CMS Preliminary");
  TLatex t2(0.6,0.92, "#sqrt{s}=13 TeV, L = 2 fb^{-1}");
  TLatex t3(0.4,0.92, Form("%s",latexlabel.c_str()) );
  t1.SetNDC();
  t2.SetNDC();
  t3.SetNDC();
  t1.SetTextSize(0.05);
  t2.SetTextSize(0.05);
  t3.SetTextSize(0.02);
  t1.SetTextFont(42);
  t2.SetTextFont(42);
  t3.SetTextFont(42);

  stackUnrolled->Draw();
  stackUnrolled->SetMinimum(0.0001);
  // stackUnrolled->SetMaximum(1000);
  cv->SetLogy();
  stackUnrolled->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolled->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolled->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolled->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackUnrolled->Draw();
  if(hasSignal) histUnrolled[0]->Draw("same hist");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled%s.pdf",Label.c_str()));

  // Unrolled plots in percentages
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.85,0.20,0.95,0.80);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }
  stackUnrolledPercentage->Draw();
  stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetXaxis()->GetTitle());
  // stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetRangeUser(0, 35);
  stackUnrolledPercentage->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetYaxis()->GetTitle());
  if(hasSignal) histUnrolledPercentage[0]->Draw("same hist"); 
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("UnrolledPercentage%s.pdf",Label.c_str()));

void fit(float bgpar2smudge=1.0) {
  TF1 *fsig = new TF1("fsig",&d_sig,0.4,2,5);
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  //int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\tstepx0\tstepx1\txsectFn\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      delete [] wq;
      int wbin = h3->GetYaxis()->FindBin(wval);
      float wlow = h3->GetYaxis()->GetBinLowEdge(wbin);
      float step_x0 = sqrt(wlow*wlow+MASS_P*MASS_P-2*wlow*MASS_P);
      float whigh = h3->GetYaxis()->GetBinLowEdge(wbin+1);
      float step_x1 = sqrt(whigh*whigh+MASS_P*MASS_P-2*whigh*MASS_P);
      fitmmp(h,step_x0,step_x1,wval);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      //int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);

      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      fsig->SetParameter(0,ftmp->GetParameter(5));
      fsig->SetParameter(1,ftmp->GetParameter(6));
      fsig->SetParameter(2,ftmp->GetParameter(7));
      fsig->SetParameter(3,step_x0);
      fsig->SetParameter(4,step_x1);
      fsig->SetLineWidth(2);
      fsig->SetLineColor(kBlue+1);
      h->GetListOfFunctions()->Add((TF1*)fsig->Clone("fsig"));
      //fsig->Print();
      double Nfn = 0;
      for (int b = 1; b < h->GetNbinsX(); b++) {
        double x = h->GetXaxis()->GetBinCenter(b);
        Nfn += fsig->Eval(x);
      }
      //printf("**** %.3e\t\%.3e\n",Nfn,N);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double xsectFn = Nfn/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile,"\t%.3e",xsectFn/(1e6));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

  /** 
   * Process a single type - i.e., one of <i>symmetric</i>,
   * <i>positive</i>, <i>negative</i>, or <i>other</i> - by looping
   * over all contained objects and call ProcessBin for each found
   * bin.
   * 
   * @param measured     Input collection of measured data
   * @param corrections  Input collection of correction data
   * @param method       Unfolding method to use 
   * @param regParam     Regularisation parameter
   * @param out          Output directory. 
   * @param sys          Collision system
   * @param sNN          Collision energy 
   */
  void ProcessType(TCollection* measured, 
		   TCollection* corrections,
		   UInt_t       method,
		   Double_t     regParam,
		   TDirectory*  out,
		   UShort_t     sys, 
		   UShort_t     sNN)
  {
    Printf("  Processing %s ...", measured->GetName());
    TDirectory* dir = out->mkdir(measured->GetName());
    
    // Make some summary stacks 
    THStack*  allMeasured  = new THStack("measured",      
					 "Measured P(#it{N}_{ch})");
    THStack*  allTruth     = new THStack("truth",        
					 "MC 'truth' P(#it{N}_{ch})");
    THStack*  allTruthA    = new THStack("truthAccepted",
					 "MC 'truth' accepted P(#it{N}_{ch})");
    THStack*  allUnfolded  = new THStack("unfolded",
					 "Unfolded P(#it{N}_{ch})");
    THStack*  allCorrected = new THStack("corrected",
					 "Corrected P(#it{N}_{ch})");
    THStack*  allRatio     = (sys != 1 ? 0 : 
			      new THStack("ratios", "Ratios to other"));
    TMultiGraph* allALICE  = (sys != 1 ? 0 : 
			      new TMultiGraph("alice", "ALICE Published"));
    TMultiGraph* allCMS    = (sys != 1 ? 0 : 
			      new TMultiGraph("cms", "CMS Published"));

    // Loop over the list of objects. 
    static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
    TIter          next(measured);
    TObject*       o = 0;
    Int_t          i = 0;
    Double_t       r = regParam;
    while ((o = next())) {
      // Go back to where we where 
      dir->cd();
      
      // if not a collection, don't bother 
      if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
    
      // If it doesn't match our regular expression, don't bother 
      TString n(o->GetName());
      if (n.Index(regex) == kNPOS) { 
	// Warning("ScanType", "%s in %s doesn't match eta range regexp", 
	//         n.Data(), real->GetName());
	continue;
      }
      TCollection* mBin = static_cast<TCollection*>(o);
      TCollection* cBin = GetCollection(corrections, n.Data());
      if (!cBin) continue;

      THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
      if (!binS) continue;

      TH1* result = 0;
      Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
		allUnfolded, allCorrected, result);

      TGraph* alice = 0;
      TGraph* cms   = 0;
      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
      Ratio2Stack(i, result, alice, cms, allRatio);
      i++;
    }
    dir->Add(allMeasured);
    dir->Add(allTruth);
    dir->Add(allTruthA);
    dir->Add(allUnfolded);
    dir->Add(allCorrected);
    if (allALICE && allALICE->GetListOfGraphs()) {
      if (allALICE->GetListOfGraphs()->GetEntries() > 0)
	dir->Add(allALICE);
      else 
	delete allALICE;
    }
    if (allCMS && allCMS->GetListOfGraphs()) {
      if (allCMS->GetListOfGraphs()->GetEntries() > 0) 
	dir->Add(allCMS);
      else 
	delete allCMS;
    }
    if (allRatio && allRatio->GetHists()) { 
      if (allRatio->GetHists()->GetEntries() > 0) 
	dir->Add(allRatio);
//.........这里部分代码省略.........

void browseStacks( bool makePictures=false, bool wait=true , bool addHistName = false, Double_t maxYScaleF = 1., 
                  bool logScale = false, bool setMinZero = true) {


  gStyle->SetOptTitle(0);

  bool keep2D=false;

  //fix the hNJet histos
  TList *list = gDirectory->GetList();
  TIterator *iter = list->MakeIterator();
  TObject *obj = 0;
  while(obj = iter->Next()) {
  
    if(TString(obj->GetName()).Contains("hnJet") && obj->InheritsFrom(TH1::Class())) {
      
      int nbins = ((TH1F*)obj)->GetNbinsX();
      float overflow = ((TH1F*)obj)->GetBinContent(nbins+1);
      float lastbinval = ((TH1F*)obj)->GetBinContent(nbins);
      ((TH1F*)obj)->SetBinContent(nbins, overflow+lastbinval);
      ((TH1F*)obj)->GetXaxis()->SetBinLabel(nbins, "#geq4");
    }
  }
  
    
    
    
  // Find out what the names of the existing histograms are
  // The histogram names are XX_YY_ZZ, where XX is the sample,
  // eg, "tt", YY is the actual name, ZZ is the final state, eg, "ee"
  TObjArray* myNames = getMyHistosNames("ttdil","ee",keep2D);
    

  // Now loop over histograms, and make stacks
  TCanvas *c = new TCanvas();
  c->Divide(2,2);
  char* suffix[4];
  suffix[0] = "ee";
  suffix[1] = "mm";
  suffix[2] = "em";
  suffix[3] = "all";
  if (makePictures) c->Print("out/stacks.ps[");
  for (int i=0; i<myNames->GetEntries(); i++) {
     
    for (int sample=0; sample<4; sample++) {
       
       
      hist::stack(Form("st_%s_%s",myNames->At(i)->GetName(),suffix[sample]),
		  Form("%s_%s$",myNames->At(i)->GetName(), suffix[sample]));
      THStack* thisStack = (THStack*) gROOT->FindObjectAny(
							   Form("st_%s_%s", myNames->At(i)->GetName(), suffix[sample]));
       
      thisStack->SetMaximum(thisStack->GetMaximum()*maxYScaleF);
      if(TString(myNames->At(i)->GetName()).Contains("hnJet")) {
	TList* histolist = thisStack->GetHists();
	int hatchcount = 0;
	// 	for(int j = 0; j<histolist->GetSize();j++) {
	// 	  if(TString(histolist->At(j)->GetName()).Contains("tt") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("tautau") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("ww") ) continue;
	// 	  hatch(histolist->At(j)->GetName(), FavoriteHatches[hatchcount]);
	// 	  hatchcount++;
	// 	}
      }
	 
	 
      TLegend* thisLeg = hist::legend(thisStack, "lpf", 0, 0, 0.75, 0.65, 0.99, 0.99);
      c->cd(sample+1);
      if (logScale) gPad->SetLogy(); else gPad->SetLogy(0);
      double stackMax = ((TH1*)thisStack->GetHists()->At(0))->GetMaximum();
      double stackMin = ((TH1*)thisStack->GetHists()->At(0))->GetMinimum();
      thisStack->SetMinimum(stackMin);
      if (setMinZero) thisStack->SetMinimum(0);
      if (logScale && stackMin <=0) thisStack->SetMinimum(1e-2*stackMax);
      if (logScale && stackMax == 0) thisStack->SetMinimum(1e-12); 
      thisStack->Draw("hist");
      string xtitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetXaxis()->GetTitle());
      string ytitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetYaxis()->GetTitle());
      thisStack->GetXaxis()->SetTitle(xtitle.c_str());
      thisStack->GetYaxis()->SetTitle(ytitle.c_str());
      TString hname = thisStack->GetName();
      if(hname.Contains("hnJet")) {
	thisStack->GetXaxis()->SetLabelSize(0.075);
	thisStack->GetYaxis()->SetLabelSize(0.05);
	thisStack->GetXaxis()->SetTitle("N_{jets}");
      }
      thisLeg->Draw();
	
      TPaveText *pt1 = new TPaveText(0.1, 0.95, 0.4, 0.999, "brNDC");
      pt1->SetName("pt1name");
      pt1->SetBorderSize(0);
      pt1->SetFillStyle(0);
	
      TText *blah;
      if (addHistName) blah = pt1->AddText(hname);
      else blah = pt1->AddText("CMS Preliminary");
      blah->SetTextSize(0.05);
      pt1->Draw();
      c->Modified();
	
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TLegend *legend = 0;
  TLatex *tex = 0;
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.7,0.53,0.90,0.88);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }

  /// Unrolled plots in bins of R&MR
  TLatex t1(0.1,0.92, "CMS Preliminary");
  TLatex t2(0.6,0.92, "#sqrt{s}=13 TeV, L = 2.1 fb^{-1}");
  TLatex t3(0.4,0.92, Form("%s",latexlabel.c_str()) );
  t1.SetNDC();
  t2.SetNDC();
  t3.SetNDC();
  t1.SetTextSize(0.05);
  t2.SetTextSize(0.05);
  t3.SetTextSize(0.02);
  t1.SetTextFont(42);
  t2.SetTextFont(42);
  t3.SetTextFont(42);
  stackUnrolled->Draw();
  stackUnrolled->SetMinimum(0.01);
  stackUnrolled->SetMaximum(1000);
  cv->SetLogy();
  stackUnrolled->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolled->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolled->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolled->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackUnrolled->Draw();
  if(hasSignal) histUnrolled[0]->Draw("same PE");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled_QCD%s.root",Label.c_str()));
  ////

  //*******************************************************************************************
  //MR
  //*******************************************************************************************
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.50,0.54,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);

  tex = new TLatex();
  tex->SetNDC();
  tex->SetTextSize(0.030);
  tex->SetTextFont(42);
  tex->SetTextColor(kBlack);
  tex->DrawLatex(0.2, 0.92, Form("CMS Simulation #sqrt{s} = 13 TeV, #int L = %d fb^{-1}, %s",int(intLumi/1000), latexlabel.c_str()));

  THStack *stackMR = new THStack("stackMR", "");
  THStack *stackRsq = new THStack();

  //*******************************************************************************************
  //MR Before and After DPhi Cut
  //*******************************************************************************************

void view()
{
  TFile* f = TFile::Open("result.root");

  int signalColorTable[20], backgroundColorTable[20];

  for ( int i=0; i<20; ++i )
  {
    signalColorTable[i] = kAzure+10-i;
    backgroundColorTable[i] = kOrange+10-i;
  }

  TList* signalPlots = makePlots((TDirectory*)f->Get("MC_Signal_EMEM"), signalColorTable);
  TList* backgroundPlots = makePlots((TDirectory*)f->Get("MC_Background_EMEM"), backgroundColorTable, true);

  if ( signalPlots == 0 || backgroundPlots == 0 ) return;

  const int nPlots = signalPlots->GetSize();
  for ( int i=0; i<nPlots; ++i )
  {
    THStack* hSignal = (THStack*)signalPlots->At(i);
    THStack* hBackground = (THStack*)backgroundPlots->At(i);

    TString histName = hSignal->GetName();
    bool doLog = histName.Contains("Pt");// || histName.Contains("RelIso");

    TCanvas* c = new TCanvas(TString("c")+hSignal->GetName(), hSignal->GetTitle(), 1200, 600);
    TPad* pad;

    c->Divide(2,1);

    TString xTitle, yTitle;

    pad = (TPad*)c->cd(1);
    if ( doLog ) pad->SetLogy();
    pad->SetBorderSize(0);
    pad->SetBorderMode(0);
    hBackground->Draw();

    xTitle = ((TH1*)hBackground->GetHists()->At(0))->GetXaxis()->GetTitle();
    yTitle = ((TH1*)hBackground->GetHists()->At(0))->GetYaxis()->GetTitle();
    hBackground->GetXaxis()->SetTitle(xTitle);
    hBackground->GetYaxis()->SetTitle(yTitle);

    pad->BuildLegend(0.6, 0.6, 0.98, 0.98);

    pad = (TPad*)c->cd(2);
    if ( doLog ) pad->SetLogy();
    pad->SetBorderSize(0);
    pad->SetBorderMode(0);
    hSignal->Draw("nostack");

    xTitle = ((TH1*)hSignal->GetHists()->At(0))->GetXaxis()->GetTitle();
    yTitle = ((TH1*)hSignal->GetHists()->At(0))->GetYaxis()->GetTitle();
    hSignal->GetXaxis()->SetTitle(xTitle);
    hSignal->GetYaxis()->SetTitle(yTitle);

    pad->BuildLegend(0.6, 0.7, 0.98, 0.98);

    c->Print(TString(c->GetName())+".png");
  }
}

void BiasDiff () {
  vector<TH1F*> AIC;
  
  for (UInt_t i = 1 ; i <= 7; i++){
    AIC.push_back( new TH1F( Form("AIC_%d",i), " ", 9, 0,10 ));
    AIC[i-1]->GetXaxis()->SetBinLabel(1,"category 0");
    AIC[i-1]->GetXaxis()->SetBinLabel(2,"");
    AIC[i-1]->GetXaxis()->SetBinLabel(3,"category 1");
    AIC[i-1]->GetXaxis()->SetBinLabel(4,"");
    AIC[i-1]->GetXaxis()->SetBinLabel(5,"category 2");
    AIC[i-1]->GetXaxis()->SetBinLabel(6,"");
    AIC[i-1]->GetXaxis()->SetBinLabel(7,"category 3");
    AIC[i-1]->GetXaxis()->SetBinLabel(8,"");
    AIC[i-1]->GetXaxis()->SetBinLabel(9,"Inclusive");
    AIC[i-1]->SetFillColor(bkgColors[i-1]);
    AIC[i-1]->SetLineColor(bkgColors[i-1]);
  }
  
  AIC[0]->Fill(1.0,0.02);
  AIC[0]->Fill(3.0,0.3089);
  AIC[0]->Fill(5.0,0.62);
  AIC[0]->Fill(7.0,0.36);
  AIC[0]->Fill(9.0,0.03);
  
  AIC[1]->Fill(1.0,0.08);
  AIC[1]->Fill(3.0,0.21);
  AIC[1]->Fill(5.0,0.08);
  AIC[1]->Fill(7.0,0.16);
  AIC[1]->Fill(9.0,0.33);
  
  AIC[2]->Fill(1.0,0.01);
  AIC[2]->Fill(3.0,0.02);
  AIC[2]->Fill(5.0,0.01);
  AIC[2]->Fill(7.0,0.02);
  AIC[2]->Fill(9.0,0.04);
  
  AIC[3]->Fill(1.0,0.20);
  AIC[3]->Fill(3.0,0.38);
  AIC[3]->Fill(5.0,0.24);
  AIC[3]->Fill(7.0,0.451);
  AIC[3]->Fill(9.0,0.60);
  
  AIC[4]->Fill(1.0,0.03);
  AIC[4]->Fill(3.0,0.08);
  AIC[4]->Fill(5.0,0.017);
  AIC[4]->Fill(7.0,0.001);
  AIC[4]->Fill(9.0,0.0001);
  
  AIC[5]->Fill(1.0,0.58);
  AIC[5]->Fill(3.0,0.001);
  AIC[5]->Fill(5.0,0.03);
  AIC[5]->Fill(7.0,0.007);
  AIC[5]->Fill(9.0,0.000003);
  
  AIC[6]->Fill(1.0,0.08);
  AIC[6]->Fill(3.0,0.0001);
  AIC[6]->Fill(5.0,0.003);
  AIC[6]->Fill(7.0,0.001);
  AIC[6]->Fill(9.0,0.0000004);
  
  TCanvas *cv = 0;
  TLegend *legend = 0;
  bool firstdrawn = false;
  
  // ===================================
  //           AIC
  // ===================================
  
  cv = new TCanvas("cv","cv",800,600);
  legend = new TLegend(0.64,0.64,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);
  
  THStack *stackAIC = new THStack();
  for (Int_t i = AIC.size()-1; i>=0; i--) {
    if (AIC[i]->Integral()>0) {
      stackAIC->Add(AIC[i]);
      legend->AddEntry(AIC[i],modelLegendLabels[i].c_str(),"F");
    }
  }
  stackAIC->Draw();
  stackAIC->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackAIC->GetHists()->At(0)))->GetXaxis()->GetTitle());
  legend->Draw();
  cv->SaveAs("AICvalues.pdf");
  
  // =======================================================

  //double x[7] = {1.0,2.0,3.0,4.0,5.0,6.0,7.0};
  double x[4] = {1.0,2.0,3.0,4.0};
  double comps[4],singEs[4],doubEs[4],tripEs[4],modEs[4],polys[4],pows[4],dpows[4];

 //Composite  
  //comps[0]  = ( (0.0042-0.0033)/0.0033+(0.0007-0.0036)/0.0036+(0.0180-0.012)/0.012+(0.0090-0.0023)/0.0023+ (0.0033-0.0031)/0.0031 ) /5;
  comps[0] = ( (0.0045-0.0043)/0.0043 + (0.0035-0.0018)/0.0018 + (0.0074-0.0084)/0.0084 + (0.0053-0.0065)/0.0065 + (0.0016-0.0029)/0.0029 ) /5;

  //comps[2] = ( (0.0036-0.0049)/0.0049 + (0.0075-0.0029)/0.0029 + (0.0125-0.0126)/0.0126 + (0.0041-0.0031)/0.0031 + (0.0053-0.0027)/0.0027 )/ 5;

  comps[1] = ( (0.0088-0.0122)/0.0122  + (0.0072-0.0064)/0.0064 + (0.0064-0.0057)/0.0057 ) / 3;

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

  //____________________________________________________________________
  void DrawResCollection(TCollection* top, const TString& name)
  {
    TCollection* c = GetCollection(top, name, false);
    if (!c) return;

    THStack* s = GetStack(c, "all");
    s->SetTitle("");
    DrawInPad(fBody, 0, s, "nostack", kLogy);
    TLegend* l = new TLegend(.5, .75, .98, .98, "P(#it{N}_{ch})");
    l->SetBorderSize(0);
    // l->SetBorderMode(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);
    TIter next(s->GetHists());
    TH1*  h = 0;
    Bool_t hasTrue = false;
    while ((h = static_cast<TH1*>(next()))) { 
      TString n(h->GetTitle());
      if (n.BeginsWith("True")) { hasTrue = true; continue; }
      n.ReplaceAll("Raw P(#it{N}_{ch}) in ", "");
      TLegendEntry* e = l->AddEntry("dummy", n, "p");
      e->SetMarkerStyle(h->GetMarkerStyle());
    }
    if (hasTrue) {
      TLegendEntry* e = l->AddEntry("dummy", "Raw", "p");
      e->SetMarkerStyle(20);
      e->SetMarkerColor(kRed+1);
      e = l->AddEntry("dummy", "MC truth", "p");
      e->SetMarkerStyle(24);
      e->SetMarkerColor(kBlue+1);
      e = l->AddEntry("dummy", "MC truth selected", "p");
      e->SetMarkerStyle(24);
      e->SetMarkerColor(kOrange+1);
    }
    fBody->cd();
    l->Draw();
    
    PrintCanvas(Form("%s results", name.Data()));

    // return;

    TIter       nextO(c);
    TObject*    o = 0;
    while ((o = nextO())) { 
      Double_t etaMin = 999;
      Double_t etaMax = 999;
      TCollection* bin = GetEtaBin(o, etaMin, etaMax);
      if (!bin) continue;

      fBody->Divide(2,3);
      DrawInPad(fBody, 1, GetH1(bin, "rawDist"),      "",     kLogy);
      DrawInPad(fBody, 1, GetH1(bin, "truthAccepted", false),
		"same", kSilent);
      DrawInPad(fBody, 1, GetH1(bin, "truth", false),"same", kSilent|kLegend);
      DrawInPad(fBody, 2, GetH1(bin, "coverage"));
      DrawInPad(fBody, 3, GetH2(bin, "corr"),            "colz");
      DrawInPad(fBody, 4, GetH2(bin, "response", false), "colz", kLogz|kSilent);
      DrawInPad(fBody, 5, GetH1(bin, "triggerVertex", false), "", kSilent);
      
      PrintCanvas(Form("%+5.1f < #eta < %+5.1f", etaMin, etaMax));
    }
  }