C++ TPad::DrawFrame方法代码示例

void LHCHiggsExample() 
{ 

#ifdef __CINT__
  gROOT->LoadMacro("LHCHiggsUtils.C");
#endif

  SetLHCHiggsStyle();

  TCanvas* c1 = new TCanvas("c1","Higgs Cross Section",50,50,600,600);
  TPad* thePad = (TPad*)c1->cd();
  thePad->SetLogy();

  Double_t ymin=1.e-2;  Double_t ymax=1.e2;
  Double_t xmin=90.00;  Double_t xmax=600.;
  TH1F *h1 = thePad->DrawFrame(xmin,ymin,xmax,ymax);
  h1->SetYTitle("#sigma_{pp #rightarrow H} [pb]");
  h1->SetXTitle("M_{H}  [GeV]");
  h1->GetYaxis()->SetTitleOffset(1.4);
  h1->GetXaxis()->SetTitleOffset(1.4);
  //h1->GetXaxis()->SetNdivisions(5);
  h1->Draw();

  myText(0.2,0.88,1,"#sqrt{s}= 13 TeV");
  myBoxText(0.55,0.67,0.05,5,"NNLO QCD");

  LHCHIGGS_LABEL(0.98,0.725);
  myText(0.2,0.2,1,"Preliminary");

  c1->Print("LHCHiggsExample.eps");
  c1->Print("LHCHiggsExample.png");
  c1->Print("LHCHiggsExample.pdf");

}

// -----------------------------------------------------------------------------
//
void test() {

  time_t start = TTimeStamp().GetSec();

  set_plot_style();

  bool draw = true;
  //bool debug = true;

  // Define analysis configuration
  PSet ps;
  defaultPSet(ps);

  // Response plots
  if (false) {
    //xSectDistr(ps);
    responseProfile();
    return;
  }
  
  // Print configuration
  std::stringstream ss;
  printPSet(ps,ss);
  std::cout << ss.str() << std::endl;
  
  // Params to store
  DoubleVV ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err;
  IntV length;
  clear( ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err, length );
  init( ps, ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err, length );
  
  // Loop through Meff bins 
  int loop = 0;
  int nloops = ps.nmeff;
  for ( int imeff = 0; imeff < ps.nmeff; ++imeff ) {

    // Generate numbers in (x1,x2) plane
    DoubleVV dalitz;
    generateTruth( ps, imeff, dalitz, true );

    // Integrate across dalitz plane
    integrate( ps, imeff, dalitz, ratio, ratio_errh, ratio_errl, 
	       pass, pass_err, fail, fail_err, length );
    
    // Labeling
    std::stringstream ss;
    ss << "Meff" << int( ps.meff_bins[imeff] );
    
    // New canvas for plots
    TCanvas* c1 = 0;
    if (draw) c1 = new TCanvas( TString("Canvas"+ss.str()), "" );
    
    // Pad for cross-section plot
    TPad* pad = 0;
    if (draw) pad = new TPad(TString("Pad"+ss.str()),"",0.,0.,1.,1.);
    if (pad) {
      pad->SetGrid();
      pad->Draw();
      pad->cd();
      pad->SetLogz();
    }
    TH1F* hr = 0;
    if (draw) hr = pad->DrawFrame(ps.min,ps.min,ps.max,ps.max);
    
    // Histo title
    if (hr) {
      std::stringstream sss;
      sss << "M_{eff}=" << ps.meff_bins[imeff] << " GeV"
	  << ", p_{T1}=" << dr(ps.pt1_bins[imeff],1) << " GeV"
	  << ", p_{T2}=" << dr(ps.pt2_bins[imeff],1) << " GeV"
	  << ", p_{T3}=" << dr(ps.pt3_bins[imeff],1) << " GeV";
      hr->SetTitle( sss.str().c_str() );
      hr->GetXaxis()->SetTitle( "x_{2}" );
      hr->GetYaxis()->SetTitle( "x_{1}" );
    }
    
    // Create 2D cross-section plot
    TH2D* his = 0;
    if (draw) his = new TH2D(TString("Histo"+ss.str()),"",
			     ps.nbins,ps.min,ps.max,
			     ps.nbins,ps.min,ps.max);
    
    //double x3 = ( 2. * ps.pt3_bins[imeff] ) / ( ps.meff_bins[imeff] + ps.pt3_bins[imeff] );
    
    // Fill 2D cross-section plot
    for ( int x2_bin = 0; x2_bin < ps.nbins; ++x2_bin ) { 
      for ( int x1_bin = 0; x1_bin < ps.nbins; ++x1_bin ) { 
// 	std::cout << " Fill:"
// 		  << " x2_bin: " << x2_bin 
// 		  << " x2: " << val(x2_bin,nbins) 
// 		  << " x1_bin: " << x1_bin 
// 		  << " x1: " << val(x1_bin,nbins) 
// 		  << " val: " << dalitz[x2_bin][x1_bin]
// 		  << std::endl;
	if (his) his->Fill( val(x2_bin,ps)+ps.width/2.,
			    val(x1_bin,ps)+ps.width/2.,
			    dalitz[x2_bin][x1_bin] ); 
      }
//.........这里部分代码省略.........

// -----------------------------------------------------------------------------
//
void talk() {
  
  time_t start = TTimeStamp().GetSec();

  bool draw = true;
  bool debug = true;
  
  // Binning
  int xbins = 100;
  int ybins = 100;
  double xmin = 0.0;
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.0;
  double xbin_centre = ( ( xmax - xmin ) / xbins ) / 2.;
  double ybin_centre = ( ( ymax - ymin ) / ybins ) / 2.;
  
  // AlphaT values
  const int nat = 1;
  double at[nat];
  for ( int ii = 0; ii < nat; ++ii ) { at[ii] = 0.55 + ii * 0.001; } 
  
  // HT regions
  const int nht = 3;
  double ht_min[nht] = { 250., 300., 350. };

  // Jet pT thresholds
  double pt1_min[nht] = { 71.4, 85.7, 100. };
  double pt2_min[nht] = { 71.4, 85.7, 100. };
  double pt3_min[nht] = { 35.7, 42.9., 50. };

  // x fractions
  double x1_min[nht];
  double x2_min[nht];
  double x3_max[nht];
  for ( int ii = 0; ii < nht; ++ii ) { x1_min[ii] = ( 2. * pt1_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  for ( int ii = 0; ii < nht; ++ii ) { x2_min[ii] = ( 2. * pt2_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  for ( int ii = 0; ii < nht; ++ii ) { x3_max[ii] = ( 2. * pt3_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  
  // Loop through bins 
  int loop = 0;
  int nloops = nht;
  for ( int iht = 0; iht < nht; ++iht ) {
    std::cout << "Completed " 
	      << 100.*float(loop)/float(nloops) 
	      << "%..." 
	      << std::endl; 
    loop++;
    
    // Labeling
    std::stringstream ss;
    ss << "HT" << int(ht_min[iht]);
	  
    // New canvas for plots
    TCanvas* c1 = 0;
    if (draw) c1 = new TCanvas( TString("Canvas"+ss.str()), "" );
	  
    // Pad for cross-section plot
    TPad* pad = 0;
    if (draw) pad = new TPad(TString("Pad"+ss.str()),"",0.,0.,1.,1.);
    if (pad) {
      pad->SetGrid();
      pad->Draw();
      pad->cd();
      pad->SetLogz();
    }
    TH1F* hr = 0;
    if (draw) hr = pad->DrawFrame(0.,0.,1.,1.);

    // Histo title
    if (hr) {
      std::stringstream sss;
      sss << "H_{T}=" << ht_min[iht]
	  << "(p_{T1},p_{T2},p_{T3})=" 
	  << pt1_min[iht] << ","
	  << pt2_min[iht] << ","
	  << pt3_min[iht] << ")"
	  << ", (x_{1},x_{2},x_{3})=" 
	  << x1_min[iht] << ","
	  << x2_min[iht] << ","
	  << x3_max[iht] << ")";
      hr->SetTitle( sss.str().c_str() );
      hr->GetXaxis()->SetTitle( "x_{2}" );
      hr->GetYaxis()->SetTitle( "x_{1}" );
    }
	  
    // Create 2D cross-section plot
    TH2D* his = 0;
    if (draw) his = new TH2D(TString("Histo"+ss.str()),"",
			     xbins,xmin,xmax,
			     ybins,ymin,ymax);
	  
    // Fill 2D cross-section plot
    for ( int xbin = 0; xbin < xbins; ++xbin ) { 
      for ( int ybin = 0; ybin < ybins; ++ybin ) { 
	double x2 = ( ( xmax - xmin ) / xbins ) * xbin + xmin;
	double x1 = ( ( ymax - ymin ) / ybins ) * ybin + ymin;
	double val = ( x1*x1 + x2*x2 ) / ( ( 1 - x1 ) * ( 1 - x2 ) ); 
//.........这里部分代码省略.........

void GE11sEfficiencyScan(int RunNumber, string RunName, string path)
{
   
   ifstream InGE11_IV_GIF, InGE11_IV, InGE11_V;

   //string path = "/home/ramkrishna/TEMP/LogFiles_TB/LogFiles306To407";	    

   string gif	= path+"/Efficiency_LC1_"+std::to_string(RunNumber)+".log";
   string IV	= path+"/Efficiency_LC2_"+std::to_string(RunNumber)+".log";
   string V	= path+"/Efficiency_LC3_"+std::to_string(RunNumber)+".log";

   cout<<"gif = "<<gif<<endl;
   InGE11_IV_GIF.open(gif);
   InGE11_IV.open(IV);
   InGE11_V.open(V);

   string rootFile = "Efficiency_Run"+std::to_string(RunNumber)+".root";
   const char *CharrootFile = rootFile.c_str();
   TFile *f = new TFile(CharrootFile,"RECREATE");
   //TTree *tree = new TTree("Run306", "Detector info for Run 306");
   TNtuple *GE11_IV_GIF = new TNtuple("GE11_IV_GIF","data from text file LC1","MeanPosOfSector:Efficiency:EfficiencyError:Nevents");
   TNtuple *GE11_IV = new TNtuple("GE11_IV","data from text file LC2","MeanPosOfSector:Efficiency:EfficiencyError:Nevents");
   TNtuple *GE11_V = new TNtuple("GE11_V","data from text file LC3","MeanPosOfSector:Efficiency:EfficiencyError:Nevents");

   Int_t nlines = 0;

   vector<double> GIF_MeanPosOfSector, GIF_Efficiency, GIF_EfficiencyError;
   vector<unsigned int> GIF_Nevents;
   unsigned int temp_Nevents;
   double temp_MeanPosOfSector, temp_Efficiency, temp_EfficiencyError;

   vector<double> IV_MeanPosOfSector, IV_Efficiency, IV_EfficiencyError;
   vector<unsigned int> IV_Nevents;

   vector<double> V_MeanPosOfSector, V_Efficiency, V_EfficiencyError;
   vector<unsigned int> V_Nevents;

   string NameOfDet, xRange;

   while (1) 
   {
     InGE11_IV_GIF >> NameOfDet >> xRange >> temp_MeanPosOfSector >> temp_Efficiency >> temp_EfficiencyError >> temp_Nevents;
	if (!InGE11_IV_GIF.good()) break;

     GIF_MeanPosOfSector.push_back(temp_MeanPosOfSector+(nlines*5));
     GIF_Efficiency.push_back(temp_Efficiency);
     GIF_EfficiencyError.push_back(temp_EfficiencyError);
     GIF_Nevents.push_back(temp_Nevents);
    	GE11_IV_GIF->Fill(temp_MeanPosOfSector+(nlines*5),temp_Efficiency,temp_EfficiencyError,temp_Nevents);
	nlines++;
	if (nlines > 20) 
	{
	    cout<<"Check the input text file for run number "<< RunNumber << endl;
	    exit(EXIT_SUCCESS);
	}
   }

   InGE11_IV_GIF.close();
   nlines=0;
   while (1) 
   {
     InGE11_IV >> NameOfDet >> xRange >> temp_MeanPosOfSector >> temp_Efficiency >> temp_EfficiencyError >> temp_Nevents;
	if (!InGE11_IV.good()) break;

     IV_MeanPosOfSector.push_back(temp_MeanPosOfSector+(nlines*5));
     IV_Efficiency.push_back(temp_Efficiency);
     IV_EfficiencyError.push_back(temp_EfficiencyError);
     IV_Nevents.push_back(temp_Nevents);
    	GE11_IV->Fill(temp_MeanPosOfSector+(nlines*5),temp_Efficiency,temp_EfficiencyError,temp_Nevents);
	nlines++;
   }

   InGE11_IV.close();
   nlines=0;
   while (1) 
   {
     InGE11_V >> NameOfDet >> xRange >> temp_MeanPosOfSector >> temp_Efficiency >> temp_EfficiencyError >> temp_Nevents;
	if (!InGE11_V.good()) break;

     V_MeanPosOfSector.push_back(temp_MeanPosOfSector+(nlines*5));
     V_Efficiency.push_back(temp_Efficiency);
     V_EfficiencyError.push_back(temp_EfficiencyError);
     V_Nevents.push_back(temp_Nevents);
    	GE11_V->Fill(temp_MeanPosOfSector+(nlines*5),temp_Efficiency,temp_EfficiencyError,temp_Nevents);
	nlines++;
   }

   InGE11_V.close();
   string CanvasName = "RunNumber"+std::to_string(RunNumber);
   const char * CharCanvasName = CanvasName.c_str();
   TCanvas* c1 = new TCanvas(CharCanvasName,"Efficiency Scan Plot",200,10,700,500);
   TPad *pad = new TPad("pad","",0,0,1,1);
   //pad->SetFillColor(42);
   pad->SetGrid();
   pad->Draw();
   pad->cd();

      // draw a frame to define the range
   TH1F *hr = pad->DrawFrame(0,-0.5,100,1.1);
   hr->SetXTitle("Detector Position (mm)");
//.........这里部分代码省略.........

void AtlasExample() 
{ 

#ifdef __CINT__
  gROOT->LoadMacro("AtlasUtils.C");
#endif

  SetAtlasStyle();

  Int_t icol1=5;
  Int_t icol2=5;

  TCanvas* c1 = new TCanvas("c1","single inclusive jets",50,50,600,600);
  // TCanvas* c1 = new TCanvas("c1","single inclusive jets");
  TPad* thePad = (TPad*)c1->cd();
  thePad->SetLogy();

  Double_t ymin=1.e-3;  Double_t ymax=2e7;
  Double_t xmin=60.00;  Double_t xmax=3500.;
  TH1F *h1 = thePad->DrawFrame(xmin,ymin,xmax,ymax);
  h1->SetYTitle("d#sigma_{jet}/dE_{T,jet} [fb/GeV]");
  h1->SetXTitle("E_{T,jet}  [GeV]");
  h1->GetYaxis()->SetTitleOffset(1.4);
  h1->GetXaxis()->SetTitleOffset(1.4);
  //h1->GetXaxis()->SetNdivisions(5);
  h1->Draw();

  const Int_t ncut=1;
  TGraphErrors *data[ncut];

  for (Int_t icut=0; icut<ncut; icut++) { // loop over cuts

    TGraphErrors *g1[nren][ncut];
    for (Int_t ir=0; ir<nren; ir++) { // loop over ren scale
      g1[ir][icut]= GetGraph(ir,ir,icut,0);
      if (g1[ir][icut]) 
	cout << g1[ir][icut]->GetTitle() << " found "  << g1[ir][icut]->GetName()  << endl;
      else { 
	cout << " g1 not  found " << endl; 
	return; 
      } 
      g1[ir][icut]->SetLineColor(1);
      g1[ir][icut]->SetMarkerStyle(0);
      //g1[ir][icut]->Draw("C");
    }

    char daname[100];
    sprintf(daname,"data_%d",icut); 
    data[icut]=(TGraphErrors*)g1[0][icut]->Clone(daname); 
    data[icut]->SetMarkerStyle(20);
    data[icut]->SetMarkerColor(1);

    // just invent some data
    for (Int_t i=0; i< data[icut]->GetN(); i++) {
      Double_t x1,y1,e,dx1=0.;
      data[icut]->GetPoint(i,x1,y1);
      Double_t r1 = 0.4*(gRandom->Rndm(1)+2);
      Double_t r2 = 0.4*(gRandom->Rndm(1)+2);
      //cout << " i= " << i << " x1= " << x1 << " y1= " << y1 << " r= " << r << endl;
      Double_t y;
      if (icut==0) y=r1*y1+r1*r2*r2*x1/50000.;
      else         y=r1*y1;
      e=sqrt(y*1000)/200;
      data[icut]->SetPoint(i, x1,y);
      data[icut]->SetPointError(i,dx1,e);
    }
    //data[icut]->Print();
  
    TGraphAsymmErrors* scale[ncut];
    TGraphAsymmErrors* scalepdf[ncut];

    scale[icut]=  myMakeBand(g1[0][icut],g1[1][icut],g1[2][icut]);
    //printf(" band1: \n");
    //scale->Print();

    scalepdf[icut]=(TGraphAsymmErrors* ) scale[icut]->Clone("scalepdf");

    TGraphErrors *gpdf[NUMPDF][ncut];
    for (Int_t ipdf=0; ipdf<NUMPDF; ipdf++) {
      gpdf[ipdf][icut]= GetGraph(0,0,icut,ipdf);
      if (gpdf[ipdf][icut]) 
	cout << gpdf[ipdf][icut]->GetTitle() << " found "  << gpdf[ipdf][icut]->GetName() << endl;
      else { 
	cout << " gpdf not  found " << endl; 
	return; 
      } 
      gpdf[ipdf][icut]->SetLineColor(2);
      gpdf[ipdf][icut]->SetLineStyle(1);
      gpdf[ipdf][icut]->SetMarkerStyle(0);
      myAddtoBand(gpdf[ipdf][icut],scalepdf[icut]); 
    }

    scalepdf[icut]->SetFillColor(icol2);
    scalepdf[icut]->Draw("zE2"); 
    scale[icut]->SetFillColor(icol1);
    scale[icut]->Draw("zE2");
    g1[0][icut]->SetLineWidth(3);
    g1[0][icut]->Draw("z");
    data[icut]->Draw("P");
    
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ratio->GetYaxis()->SetTitle("S/(#sqrt{S+B})");
    ratio->SetMarkerSize(1.1);

    ratioW = new TH1F(("ratioW_"+string(denominator.at(0)->GetName())).c_str(),"",numTotal->GetNbinsX(),numTotal->GetBinLowEdge(1),numTotal->GetBinLowEdge(numTotal->GetNbinsX()+1));
    ratioW->GetYaxis()->SetTitle("weighted S/(#sqrt{S+B})");
    ratioW->SetMarkerSize(1.1);

    TString name = "norm_" ;
    name += denTotal->GetName () ;
    TH1F * norm_denTotal = (TH1F *) denTotal->Clone (name) ;
    norm_denTotal->Scale (1. / norm_denTotal->GetMaximum ()) ; 
    // weight the S/sqrt (B) by the shape of the total, 
    // so that only bins with a lot of stats become visibly significant
    for(int iBin = 0; iBin < ratio->GetNbinsX()+1; iBin++){
      if(denTotal->GetBinContent(iBin) !=0){ 
        ratioW->SetBinContent(iBin,
			      norm_denTotal->GetBinContent (iBin) * numTotal->GetBinContent(iBin) / 
			      sqrt(denTotal->GetBinContent(iBin)));
        ratio->SetBinContent(iBin,
			     numTotal->GetBinContent(iBin) / sqrt(denTotal->GetBinContent(iBin)));
      }
      else 
	ratio->SetBinContent(iBin,0.);
    }
 
    ratio->GetXaxis()->SetTitle("");
    ratio->SetLineColor(kBlue);
    ratio->SetLineStyle(2);
    ratio->SetLineWidth(2);
    ratio->GetXaxis()->SetLabelOffset(999);
    ratio->GetXaxis()->SetLabelSize(0);
    ratio->GetYaxis()->SetLabelSize(0.15);
    ratio->GetYaxis()->SetTitleSize(0.15);
    ratio->GetYaxis()->SetTitleOffset(0.30);
    ratio->GetYaxis()->SetNdivisions(504);

    ratioW->GetXaxis()->SetTitle("");
    ratioW->SetLineColor(kBlack);
    ratioW->SetLineWidth(2);
    ratioW->GetXaxis()->SetLabelOffset(999);
    ratioW->GetXaxis()->SetLabelSize(0);
    ratioW->GetYaxis()->SetLabelSize(0.15);
    ratioW->GetYaxis()->SetTitleSize(0.15);
    ratioW->GetYaxis()->SetTitleOffset(0.30);
    ratioW->GetYaxis()->SetNdivisions(504);

    ratio->GetYaxis()->SetRange(min(ratio->GetMinimum(),ratioW->GetMinimum())*0.9,max(ratio->GetMaximum(),ratioW->GetMaximum())*1.1);    


    TH1F * frame = lowerPad->DrawFrame (ratio->GetXaxis ()->GetXmin (), 0., 
                                        ratio->GetXaxis ()->GetXmax (), 2.) ;
    frame->GetXaxis()->SetTitle (ratio->GetXaxis ()->GetTitle ()) ;
    frame->GetYaxis()->SetTitle (ratio->GetYaxis ()->GetTitle ()) ;
    frame->GetXaxis()->SetLabelOffset(999);
    frame->GetXaxis()->SetLabelSize(0);
    frame->GetYaxis()->SetLabelSize(0.15);
    frame->GetYaxis()->SetTitleSize(0.15);
    frame->GetYaxis()->SetTitleOffset(0.30);
    frame->GetYaxis()->SetNdivisions(504);

    ratio->Draw("P");    
    ratioW->Draw("Lsame");    
    
    upperPad->cd();
    
    tex->Draw("same");
    tex2->Draw("same");
    tex3->Draw("same");
    legend->Draw("same");

    cCanvas->SaveAs(string("output/"+outputPlotDirectory+"/xs/"+variableList.at(iVar).variableName+".pdf").c_str(),"pdf");
    cCanvas->SaveAs(string("output/"+outputPlotDirectory+"/xs/"+variableList.at(iVar).variableName+".png").c_str(),"png");
    cCanvas->SaveAs(string("output/"+outputPlotDirectory+"/xs/"+variableList.at(iVar).variableName+".root").c_str(),"root");

    cCanvasNorm->cd();

    tex->Draw("same");
    tex2->Draw("same");
    tex3->Draw("same");
    legend->Draw("same");

    cCanvasNorm->SaveAs(string("output/"+outputPlotDirectory+"/norm/"+variableList.at(iVar).variableName+".pdf").c_str(),"pdf");
    cCanvasNorm->SaveAs(string("output/"+outputPlotDirectory+"/norm/"+variableList.at(iVar).variableName+".png").c_str(),"png");
    cCanvasNorm->SaveAs(string("output/"+outputPlotDirectory+"/norm/"+variableList.at(iVar).variableName+".root").c_str(),"root");

    
    legend->Clear();
    
  } // loop on var

  cout<<"LHE filter efficiency : "<<passingLHEFilter<<" totEvent "<<totEvent<<" efficiency "<<float(passingLHEFilter)/float(totEvent)*100<<" % "<<endl;
  
  //Normalize histograms
  for(size_t ihisto = 0; ihisto < plotVector.size(); ihisto++){
    if(plotVector.at(ihisto).varName == "DeltaPhi_LL")
      cout<<"Events Histo "<<plotVector.at(ihisto).histogram->GetName()<<" unweighted "<<plotVector.at(ihisto).histogram->GetEntries()<<" weighted "<<plotVector.at(ihisto).histogram->Integral(0,plotVector.at(ihisto).histogram->GetNbinsX()+1)<<endl;
  }	          

  return 0 ;
}  

// -----------------------------------------------------------------------------
//
void rob() {

  // Binning
  
  int xbins = 100;
  int ybins = 100;
  
  double xmax = 1.0;
  double xmin = 0.0;
//   double xrange = xmax - xmin;
//   xmax += xrange / xbins / 2;
//   xmin -= xrange / xbins / 2;
//   //xbins++;
  
  double ymax = 1.0;
  double ymin = 0.0;
//   double yrange = ymax - ymin;
//   ymax += yrange / ybins / 2;
//   ymin -= yrange / ybins / 2;
//   //ybins++;

  std::cout << " Binning: "
	    << " xbins: " << xbins
	    << " ybins: " << xbins
	    << " xmin: " << xmin
	    << " xmax: " << xmax
	    << " ymin: " << ymin
	    << " ymax: " << ymax
	    << std::endl;
  
  // Examples values of pt1, pt2, mht, x1, x2, x3, sigma and alpha_t
//   double pt1 = 50.;
//   double pt2 = 50.;
//   double mht = 50.;
//   double x1 = ( 2. * pt1 ) / ( pt1 + pt2 + mht );
//   double x2 = ( 2. * pt2 ) / ( pt1 + pt2 + mht );
//   double x3 = 2 - x1 - x2;
//   double sigma =  ( x1*x1 + x2*x2 ) / ( (1-x1) * (1-x2) );
//   double alpha_t =  x2 / ( 2 * sqrt(x1+x2-1) );
//   std::cout << " pt1: " << pt1
// 	    << " pt2: " << pt2
// 	    << " mht: " << mht 
// 	    << " x1: " << x1
// 	    << " x2: " << x2
// 	    << " x3: " << x3
// 	    << " sigma: " << sigma
// 	    << " alpha_t: " << alpha_t
// 	    << std::endl;
  
  // Cross section
  TCanvas* c1 = new TCanvas( "Contours", "" );
  //c1->SetGridx(1);
  //c1->SetGridy(1);

  TPad* pad = new TPad("pad","",0.,0.,1.,1.);
  pad->SetGrid();
  pad->Draw();
  pad->cd();
  pad->SetLogz(1);

  TH1F* hr = pad->DrawFrame(0.,0.,1.,1.);
  
  const int nx = 3;
  const int ny = 10;
  double pt[nx] = { 30., 50., 100. };
  double ht[ny];
  for ( int ii = 0; ii < ny; ++ii ) { ht[ii] = 150. + ii * 50.; }
  double ratio[nx][ny];
  
  double xbin_centre = ( ( xmax - xmin ) / xbins ) / 2.;
  double ybin_centre = ( ( ymax - ymin ) / ybins ) / 2.;
  TH2D* his = new TH2D("Contours","",xbins,xmin,xmax,ybins,ymin,ymax);
  for ( int ii = 0; ii < nx; ++ii ) {
    for ( int jj = 0; jj < ny; ++jj ) {
      double x3 = ( 2. * pt[ii] ) / ( ht[jj] + pt[ii] );
      double n = 0.;
      double d = 0.;
      for ( int xbin = 0; xbin < xbins; ++xbin ) { 
	for ( int ybin = 0; ybin < ybins; ++ybin ) { 
	  double x2 = ( ( xmax - xmin ) / xbins ) * xbin + xmin;
	  double x1 = ( ( ymax - ymin ) / ybins ) * ybin + ymin;
	  double val = ( x1*x1 + x2*x2 ) / ( ( 1 - x1 ) * ( 1 - x2 ) ); 
	  if ( x1 < x2 ||      // jet ordering by Pt
	       x1 + x2 > 2. || // from relation "x1 + x2 + x3 = 2"
	       x1 > 1.0 ||     // from "lost jet" and relation "xmiss = -x1 -x2"
	       x1 + x2 < 1.    // from "lost jet" and relation "xmiss = -x1 -x2"
	       ) { continue; }
	  if ( ( x1 + x2 ) < ( 2 - x3 ) ) { continue; }
	  d += val;
	  double alpha_t =  x2 / ( 2 * sqrt(x1+x2-1) );
	  if ( alpha_t > 0.5 ) n+= val;
	  if ( ii == 2 && jj == 9 ) { his->Fill( x2+xbin_centre, x1+ybin_centre, val ); }
	}
      }
      double r = 0.;
      if ( d > 0. ) { r = n/d; }
      ratio[ii][jj] = r;
      std::cout << " Pt: " << pt[ii]
//.........这里部分代码省略.........

void QA_Draw_Jet_Summary(const char *jet_family = "AntiKt_Tower",
                         const char *qa_file_name_new =
                             "data/G4sPHENIXCells_2000jets25GeV.root_qa.root",
                         const char *qa_file_name_ref =
                             "data/G4sPHENIXCells_250jets25GeV.root_qa.root")
{
  //! drawing energy range
  const double min_Et = 10;
  const double max_Et = 80;

  SetsPhenixStyle();
  TVirtualFitter::SetDefaultFitter("Minuit2");

  // file IO
  TFile *qa_file_new = new TFile(qa_file_name_new);
  assert(qa_file_new->IsOpen());

  // buffer for results
  vector<float> vec_radius;
  vector<TGraphErrors *> vec_phi_res;
  vector<TGraphErrors *> vec_eta_res;
  vector<TGraphErrors *> vec_e_res;
  vector<TGraphErrors *> vec_et_res;
  vector<TGraphErrors *> vec_reco_eff;
  vector<TGraphErrors *> vec_purity;

  // list and process all jets
  TList *hist_key_list = qa_file_new->GetListOfKeys();
  for (int i = 0; i < hist_key_list->GetSize(); ++i)
  {
    TString key_name = hist_key_list->At(i)->GetName();

    TString s_re_fullname = Form(
        "h_QAG4SimJet_.*_r[0-9]*_%s_r[0-9]*_Matching_Count_Truth_Et",
        jet_family);  // regular expression for search
    TRegexp re_fullname(s_re_fullname, false);
    if (key_name.Index(re_fullname) == kNPOS)
      continue;

    //      cout << " key_name = " << key_name << endl;
    TString jet_pair_name = key_name(0,
                                     key_name.Length() - TString("_Matching_Count_Truth_Et").Length());  // remove suffix

    //      cout << " jet_pair_name = " << jet_pair_name << endl;

    //get jet radius
    TRegexp re_jetradius("_r[0-9]*", false);
    Ssiz_t index_radius = key_name.Index(re_jetradius);             // first radius
    index_radius = key_name.Index(re_jetradius, index_radius + 1);  // second radius
    assert(index_radius != kNPOS);
    float radius = 0;
    sscanf(key_name(index_radius, 100).Data(), "_r%f", &radius);
    //      cout << " index_radius = " << index_radius << endl;
    assert(radius != 0);
    radius /= 10;  // jet radius convention in DST names

    cout << "QA_Draw_Jet_Summary - process jet pair " << jet_pair_name
         << " with radius = " << radius << endl;

    vector<TGraphErrors *> resolution_efficiency_summary(
        QA_Draw_Jet_TruthMatching(jet_pair_name, qa_file_name_new,
                                  qa_file_name_ref));

    //save results
    vec_radius.push_back(radius);
    vec_phi_res.push_back(resolution_efficiency_summary[0]);
    vec_eta_res.push_back(resolution_efficiency_summary[1]);
    vec_e_res.push_back(resolution_efficiency_summary[2]);
    vec_et_res.push_back(resolution_efficiency_summary[3]);
    vec_reco_eff.push_back(resolution_efficiency_summary[4]);
    vec_purity.push_back(resolution_efficiency_summary[5]);

    //      break;
  }

  // plot
  TCanvas *c1 = new TCanvas(
      TString("QA_Draw_Jet_Summary_") + TString(jet_family),
      TString("QA_Draw_Jet_Summary_") + TString(jet_family), 1800, 900);
  c1->Divide(3, 2);
  int idx = 1;
  TPad *p;

  // ------------------------------------
  p = (TPad *) c1->cd(idx++);
  c1->Update();
  //  p->SetLogz();

  TH1 *h_frame =
      p->DrawFrame(min_Et, -.1, max_Et, .1,
                   TString(jet_family) + " #phi Reconstruction;E_{T, Truth} (GeV);#phi_{Reco} - #phi_{Truth} (rad)");
  //  h_frame->GetYaxis()->SetTitleOffset(1.01);
  TLine *l = new TLine(min_Et, 0, max_Et, 0);
  l->Draw();
  p->SetGridx(0);
  p->SetGridy(0);
  TLegend *legend = new TLegend(0.7, 0.2, .95, 0.5);
  legend->SetFillColor(kWhite);
  legend->SetFillStyle(1001);
  legend->SetLineWidth(2);
//.........这里部分代码省略.........