C++ TH2::Draw方法代码示例

void plot_corr(TH2* phin, float cmin, float cmax) {
  ostringstream ssnam;
  ssnam << phin->GetName() << "_" << cmin << "_" << cmax;
  new TCanvas;
  TH2* ph = phin->Clone(ssnam.str().c_str());
  ph->GetXaxis()->SetRangeUser(cmin, cmax);
  ph->GetYaxis()->SetRangeUser(cmin, cmax);
  ph->Draw("colz");
  addaxis();
}

void
polar()
{
   TCanvas *c1 = new TCanvas("c1","c1",500,500);
   TGraphPolar * grP1 = new TGraphPolar();
   grP1->SetTitle("TGraphPolar example");

   for (Int_t i = 0; i < 2; i++) {
     grP1->SetPoint(i, (i+1) * TMath::Pi() / 4, (i+1) * 0.05);
     grP1->SetPointError(i, 9*TMath::Pi()/180, 0.0);
     Double_t rr = grP1->GetY()[i];
     Double_t tt = grP1->GetX()[i];
     Double_t x  = rr * TMath::Cos(tt);
     Double_t y  = rr * TMath::Sin(tt);
     Printf("(x,y)=(%f,%f)", x, y);
   }
   Double_t r = 1;
   TH2* frame = new TH2F("frame", "Frame", 100, -r,r, 100, -r, r);
   frame->Draw();

   grP1->SetMarkerStyle(1);
   grP1->SetMarkerSize(1.);
   grP1->SetMarkerColor(4);
   grP1->SetLineColor(4);
   grP1->SetLineWidth(3);
   grP1->SetFillColor(kRed+1);
   grP1->SetFillStyle(3001);
   grP1->Draw("PNEF same");
   // grP1->Draw("APNEF");

   // Update, otherwise GetPolargram returns 0
   c1->Update();
   TGraphPolargram* gram = grP1->GetPolargram();
   gram->SetLineWidth(0);
   // gram->SetLineColor(kWhite);
   gram->SetNdivPolar(20);
   gram->SetNdivRadial(10);   
   gram->SetTextSize(0);
   gram->SetRadialLabelSize(0);
   gram->SetPolarLabelSize(0);
   gram->SetAxisAngle(9*TMath::Pi()/180);
   gram->SetTwoPi();
   gram->SetToRadian();
   c1->SetGridx();
   c1->SetGridy();
   c1->Update();
}

void root_hist_test()
{
  // Project3D uses overflows?
  TH3* h3D = new TH3D("h3D", "h3D", 2,0,2, 2,0,2, 2,0,2);
  for(int i=0; i<=h3D->GetNbinsX()+1; i++) {
    for(int j=0; j<=h3D->GetNbinsY()+1; j++) {
      for(int k=0; k<=h3D->GetNbinsZ()+1; k++) {
        h3D->SetBinContent(i,j,k, 1);
      }
    }
  }

  printf("h3D integral = %d\tentries = %d\n", h3D->Integral(), h3D->GetEntries());

  h3D->GetZaxis()->SetRange(3,3);
  TH2* hYX = (TH2*) h3D->Project3D("yx");

  printf("hYX integral = %d\tentries = %d\n", hYX->Integral(), hYX->GetEntries());
  hYX->Draw("colz");
}

  void Draw() 
  {
    float min = 100000;
    float max = 0;
    std::vector<TGraph*> gs;
    float_array bg(fTKine.size());
    for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass();
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      for (size_t i = 0; i < fTKine.size(); i++) {
	if ((*j)->fValues[i] == 0) continue;
	min = std::min(min, (*j)->fValues[i]);
	max = std::max(max, (*j)->fValues[i]);
      }
    }

    TCanvas* c = new TCanvas("c", "C", 700, 700);
    c->SetFillColor(0);
    c->SetLogy();
    c->SetLogx();
    c->SetGridy();
    c->SetGridx();
    float_array y(101);
    float ymin = log10(min);
    float dy   = (log10(max)+.5 - log10(min)) / y.size();
    for (size_t i = 1; i < y.size(); i++)  y[i] = pow(10, ymin + i * dy);
    TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0]));
    f->SetXTitle("#beta#gamma");
    f->SetDirectory(0);
    f->SetStats(kFALSE);
    f->Draw();
    TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
    l->SetFillColor(0);
    // l->SetFillStyle(0);
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      TGraph& g = (*j)->Draw();
      l->AddEntry(&g, g.GetTitle(), "l");
    }
    l->Draw("same");
  }

// Make 2D efficiency plots
void makeplots2D( TH2& eff, TString name) {
  gROOT->ProcessLine(".L ~/tdrstyle.C");
  setTDRStyle();

  const Int_t NRGBs = 5;
  const Int_t NCont = 200;
  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  tdrStyle->SetNumberContours(NCont);

  
  if(name.Contains("_Pt")) {
    eff.GetXaxis()->SetTitle("p_{T} (GeV/c)");
    eff.GetYaxis()->SetTitle("#eta     ");
  }

  if(name.Contains("_Phi")) {
    eff.GetXaxis()->SetTitle("#phi     ");
    eff.GetYaxis()->SetTitle("#eta     ");
  }

  eff.GetYaxis()->SetTitleOffset(1);


  tdrStyle->SetPadLeftMargin(0.08);
  tdrStyle->SetPadRightMargin(0.1);
 
  TCanvas canvas("canvas",name,600,600);
  eff.Draw("colz");
  gPad->Update();
  TPaletteAxis* palette = 
    (TPaletteAxis*)eff.GetListOfFunctions()->FindObject("palette");
  palette->SetLabelSize(0.02);
  canvas.SaveAs(name + TString(".eps"));
  canvas.SaveAs(name + TString(".gif"));
  canvas.Close();

}

void DrawGrid()
{

Int_t ncol = 100;
  Int_t nrow = 200;
  
  
   TPad *grid = new TPad("grid","",0,0,1,1);
   grid->Draw();
   grid->cd();
   grid->SetGrid();
   grid->SetFillStyle(4000);
   grid->SetFrameFillStyle(0);


   TH2 *hgrid = new TH2C("hgrid","", ncol+1, -5, ncol+5, nrow, -5, nrow-1+5);
   hgrid->Draw();
   hgrid->GetXaxis()->SetNdivisions(6,100);
   hgrid->GetYaxis()->SetNdivisions(200);
   hgrid->GetYaxis()->SetLabelOffset(999.);
   hgrid->GetXaxis()->SetLabelOffset(999.);
}

void IslandLoss()
{
  /*****************************************************************/
  // Prepare the canvas
  gStyle->SetOptStat(""); 
  TCanvas *canvas = (TCanvas *) gROOT->GetListOfCanvases()->At(0);
  canvas->Clear();

  /*****************************************************************/
  TH2 *hExpected = (TH2 *)gDirectory->Get("hExpectedNumberOfIslands");
  TH2 *hActual = (TH2 *)gDirectory->Get("hNumberOfIslands");
  if (hExpected && hActual)
    {
      TVirtualPad *pad = canvas->cd(1);

      hExpected->Add(hActual, -1);
      hExpected->SetTitle("Difference between Expected Number of TPulseIslands and Actual Number of TPulseIslands");
      hExpected->GetXaxis()->SetRange(9,9);
      hExpected->GetYaxis()->SetRangeUser(0, 1200);
      hExpected->Draw("COLZ");
    }
}

void DrawHistos(TString infile = "histos.root"){
    set_plot_style();
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);
    
    TFile *_file0 = TFile::Open(infile);
    TH1* hHardPartonspT = (TH1*)gROOT->FindObject("hHardPartonspT");
    TH1* hHardPartonsName = (TH1*)gROOT->FindObject("hHardPartonsName");
    TH1* hParticlepT = (TH1*)gROOT->FindObject("hParticlepT");
    TH2* hParticleEtaPhi = (TH2*)gROOT->FindObject("hParticleEtaPhi");

    TCanvas* can1 = new TCanvas("can1","",600,600);
    gPad->SetLeftMargin(0.15);
    gPad->SetLogy();
    SetHist(hParticlepT,"p_{T} (GeV/c)","Counts",2);
    SetHist(hHardPartonspT,"p_{T} (GeV/c)","Counts",1);
    hParticlepT->Draw();
    hHardPartonspT->Draw("same");
    TLegend* leg1 = new TLegend(0.6,0.7,0.85,0.85);
    SetLeg(leg1);
    leg1->AddEntry(hHardPartonspT,"Hard Parton","L");
    leg1->AddEntry(hParticlepT,"Final Particles","L");
    leg1->Draw();
    can1->SaveAs("pTgraph.png");
    
    TCanvas* can2 = new TCanvas("can2","",600,600);
    gPad->SetLeftMargin(0.15);
    SetHist(hHardPartonsName,"parton ID","Counts",1);
    hHardPartonsName->Draw();
    
    TCanvas* can3 = new TCanvas("can3","",600,600);
    gPad->SetLeftMargin(0.15);
    SetHist(hParticleEtaPhi,"p_{T} (GeV/c)","Counts",1);
    hParticleEtaPhi->Draw("COLZ");
   
    
  }

TList* getContours(const TH2* hist,double contourLevel,const TString filename)
{
  TH2* h = (TH2*)hist->Clone("_clone");
  double limitValue[1] = {contourLevel};
  h->SetContour(1,limitValue);
  TCanvas* c = new TCanvas("contour_canvas","Contour List",0,0,600,600);
  h->Draw("CONT LIST");
  c->Update();
  TList* contours = (TList*)((TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours"))->At(0);
  TGraph* contour = (TGraph*)contours->First();
  if(filename!="") 
    {
      for(unsigned int j = 0; j < contours->GetSize(); j++)
	{
	  TString newFilename = filename+"_";
	  newFilename+=j;
	  contour->SaveAs(newFilename+".C");
	  contour = (TGraph*)contours->After(contour); // Get Next graph
	}
    }
  delete h;
  delete c;
  return contours;
}

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

  
  
}

int drawPerformanceTPCQAMatch(const char* inFile = "perf.root") {
  //
  // Draw control histograms
  // and generate output pictures
  //

  gSystem->Load("libSTAT");
  gSystem->Load("libANALYSIS");
  gSystem->Load("libANALYSISalice");
  gSystem->Load("libANALYSIScalib");
  gSystem->Load("libCORRFW");
  gSystem->Load("libTPCcalib");
  gSystem->Load("libTRDcalib");
  gSystem->Load("libT0calib");
  gSystem->Load("libTOFcalib");
  gSystem->Load("libANALYSISalice");
  gSystem->Load("libANALYSIScalib");
  gSystem->Load("libTender");
  gSystem->Load("libPWGPP");

  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  gStyle->SetTitleSize(0.025);
  TH1::AddDirectory(kFALSE);

  //
  // set criteria
  //
  Float_t mineta = -0.8;
  Float_t maxeta = 0.799;
  Float_t minNclust = 70.0;
  Double_t ptMax = 10;
  //
  // open input file
  //
  TFile *file = TFile::Open(inFile);
  if(!file)
    return -9;
  cout<<"QA file opened"<<endl;
  file->cd();
  // get the TPC list
  if(gROOT->FindObject("TPC_PerformanceQA")) TPC_PerformanceQA->cd();
  cout<<"TPC_PerformanceQA opened"<<endl;
  TList *TPC = (TList*)gROOT->FindObject("TPCQA");
  if(TPC==NULL) TPC = (TList*)gROOT->FindObject("TPCQA_v0_c0");
  if(TPC==NULL) TPC = (TList*)gROOT->FindObject("TPCQA_v0_c30");
  if(TPC==NULL) TPC = (TList*)gROOT->FindObject("TPCQA_v0_c70");
  if(TPC==NULL) TPC = (TList*)gROOT->FindObject("TPC");
  if(TPC==NULL) return(0);
  cout<<"TPCQA list found"<<endl;
  // get TPC performance object
  AliPerformanceTPC *obj = (AliPerformanceTPC*)TPC->FindObject("AliPerformanceTPC");
  if(obj==NULL) return(0);
  cout<<"what about here after obj  "<<endl;
  // get folder with histograms
  TFolder *fold = obj->GetAnalysisFolder();
  if(!fold) return(0);
  cout<<"what about here after folder  "<<endl;
  //
  // get the HLT list
  //	file->cd();
  //	if(gROOT->FindObject("HLT_PerformanceQA")) HLT_PerformanceQA->cd();
  //	TList *HLT = (TList*)gROOT->FindObject("HLTQA");


  //
  // Draw histograms
  //

  //
  // event level
  //



  TH1 *h1D = 0;
  TH2 *h2D = 0;
  TH3 *h3D = 0;

  h1D = (TH1*)fold->FindObject("h_tpc_event_1");
  Double_t NEvents = h1D->GetEntries();

  cout<<"number of events    "<<NEvents<<endl;

  TCanvas *can1 = new TCanvas("can1","TPC event information",1200,800);
  can1->Divide(3,2);

  can1->cd(1);
  fold->FindObject("h_tpc_event_6")->Draw("histe");

  can1->cd(2);
  gPad->SetLogy();
  h1D = (TH1*)fold->FindObject("h_tpc_event_recvertex_0");
  h1D->GetXaxis()->SetRangeUser(-1.,1.);
  h1D->Draw("histe");

  can1->cd(3);
  gPad->SetLogy();
//.........这里部分代码省略.........

int Execute(int argc, char *argv[]){
	if(argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0)){
		help(argv[0]);
		return 1;
	}
	else if(argc < 3){
		std::cout << " Error: Invalid number of arguments to " << argv[0] << ". Expected 2, received " << argc-1 << ".\n";
		help(argv[0]);
		return 1;
	}

	std::string draw_string = "";
	std::string save_filename = "";
	std::string save_histname = "";
	bool create_hist = false;
	int index = 3;
	while(index < argc){
		if(strcmp(argv[index], "--draw") == 0){
			if(index + 1 >= argc){
				std::cout << " Error! Missing required argument to '--draw'!\n";
				help(argv[0]);
				return 1;
			}
			draw_string = std::string(argv[++index]);
			create_hist = true;
		}
		else if(strcmp(argv[index], "--save") == 0){
			if(index + 2 >= argc){
				std::cout << " Error! Missing required argument to '--save'!\n";
				help(argv[0]);
				return 1;
			}
			save_filename = std::string(argv[++index]);
			save_histname = std::string(argv[++index]);
		}
		else{ 
			std::cout << " Error! Unrecognized option '" << argv[index] << "'!\n";
			help(argv[0]);
			return 1;
		}
		index++;
	}

	TApplication* rootapp = new TApplication("rootapp", 0, NULL);
	
	TCanvas *can1 = new TCanvas("can1", "Canvas");
	can1->cd()->SetLogz();

	TH2 *hist = NULL;
	TFile *file = new TFile(argv[1], "READ");
	if(!file->IsOpen()){
		std::cout << " Error! Failed to open input file '" << argv[1] << "'.\n";
		return 1;
	}

	if(!create_hist){
		hist = (TH2*)file->Get(argv[2]);
		if(!hist){
			std::cout << " Error! Failed to load input histogram '" << argv[2] << "'.\n";
			file->Close();
			return 1;
		}
	}
	else{	
		TTree *tree = (TTree*)file->Get(argv[2]);
		if(!tree){
			std::cout << " Error! Failed to load input histogram '" << argv[2] << "'.\n";
			file->Close();
			return 1;
		}
	
		std::cout << " " << argv[2] << "->Draw(\"" << draw_string << "\", \"\", \"COLZ\");\n";
		std::cout << " Filling TH2... " << std::flush;	
		//tree->Draw("vandle.loc:vandle.ctof>>(250,-80,-20,42,16,58)"); // VANDLE
		//tree->Draw("liquidbar.loc:liquidbar.ctof>>(250,-80,-20,10,64,74)"); // NeuDLES
		//tree->Draw("generic.loc:generic.tof>>(250,-80,-20,10,80,90)"); // HAGRiD
		tree->Draw(draw_string.c_str(), "", "COLZ");
		std::cout << "DONE\n";
	
		hist = (TH2*)tree->GetHistogram();	
	}

	if(!save_filename.empty()){
		TFile *ofile = new TFile(save_filename.c_str(), "RECREATE");
		ofile->cd();
		hist->Write(save_histname.c_str());
		ofile->Close();
		std::cout << " Saved histogram '" << save_histname << "' to file '" << save_filename << "'.\n";
	}

	hist->Draw("COLZ");
	can1->WaitPrimitive();

	TCanvas *can2 = new TCanvas("can2", "Canvas");

	Process(hist, can2);

	file->Close();
	can1->Close();
	can2->Close();
//.........这里部分代码省略.........

//______________________________________________________________________________
void DrawEvolution(const char* file, bool normalized=true)
{

  TFile* f = TFile::Open(gSystem->ExpandPathName(file));

  if (!f) return;

  TCanvas* c = new TCanvas("mch-status-evolution","mch-status-evolution");

  c->SetGridy();
  c->SetTicky();

  c->Draw();

  TLegend* l = new TLegend(0.1,0.7,0.3,0.95,"ch evolution");

  TGraph* g = static_cast<TGraph*>(f->Get("nbad"));
  if (!g) return;

  int runmin = TMath::Nint(g->GetX()[0]);
  int runmax = TMath::Nint(g->GetX()[g->GetN()-1]);

  cout << Form("Run range found in file %s = %d - %d",file,runmin,runmax) << endl;

  double ymax(0.4);

  TH2* h = new TH2F("hframe","hframe;Run number;Fraction of dead channels",100,runmin-200,runmax+200,100,0,ymax);

  gStyle->SetOptStat(kFALSE);
  h->Draw();
  h->GetXaxis()->SetNoExponent();
  h->GetXaxis()->SetNdivisions(505);

  gStyle->SetOptTitle(kFALSE);

  DrawPeriod(runmin,runmax,115881,117222,0,ymax,"10b");
  DrawPeriod(runmin,runmax,119159,120824,0,ymax,"10c");
  DrawPeriod(runmin,runmax,122374,126424,0,ymax,"10d");
  DrawPeriod(runmin,runmax,127724,130850,0,ymax,"10e");
  DrawPeriod(runmin,runmax,133005,134929,0,ymax,"10f");
  DrawPeriod(runmin,runmax,135658,136376,0,ymax,"10g");
  DrawPeriod(runmin,runmax,137133,139513,0,ymax,"10h");
  DrawPeriod(runmin,runmax,143856,146860,0,ymax,"11a");
  DrawPeriod(runmin,runmax,148370,150702,0,ymax,"11b");
  DrawPeriod(runmin,runmax,151566,154583,0,ymax,"11c");
  DrawPeriod(runmin,runmax,158084,159606,0,ymax,"11d");
  DrawPeriod(runmin,runmax,160677,162717,0,ymax,"11e");
  DrawPeriod(runmin,runmax,162933,165744,0,ymax,"11f");
  DrawPeriod(runmin,runmax,167703,170593,0,ymax,"11h");

  // 2012

  DrawPeriod(runmin,runmax,176661,177295,0,ymax,"12a");
  DrawPeriod(runmin,runmax,177384,178053,0,ymax,"12b");
  DrawPeriod(runmin,runmax,179603,180569,0,ymax,"12c");
  DrawPeriod(runmin,runmax,183913,186320,0,ymax,"12d");
  DrawPeriod(runmin,runmax,186365,186602,0,ymax,"12e");
  DrawPeriod(runmin,runmax,186668,188123,0,ymax,"12f");
  DrawPeriod(runmin,runmax,188362,188503,0,ymax,"12g");
  DrawPeriod(runmin,runmax,189122,190110,0,ymax,"12h");

  // 2013

  DrawPeriod(runmin,runmax,195344,195483,0,ymax,"13b");
  DrawPeriod(runmin,runmax,195529,195677,0,ymax,"13c");
  DrawPeriod(runmin,runmax,195681,195873,0,ymax,"13d");
  DrawPeriod(runmin,runmax,195949,196311,0,ymax,"13e");
  DrawPeriod(runmin,runmax,196433,197388,0,ymax,"13f");

  // 2015
  // periods are from the logbook, taking only PHYSICS% partitions
  // into account
  DrawPeriod(runmin,runmax,213329,215151,0,ymax,"15a");
  DrawPeriod(runmin,runmax,215160,215579,0,ymax,"15b");
  DrawPeriod(runmin,runmax,215580,219968,0,ymax,"15c");
  DrawPeriod(runmin,runmax,220050,223195,0,ymax,"15d");
  DrawPeriod(runmin,runmax,223227,224779,0,ymax,"15e");
  DrawPeriod(runmin,runmax,224826,226606,0,ymax,"15f");
  DrawPeriod(runmin,runmax,228936,231321,0,ymax,"15g");
  DrawPeriod(runmin,runmax,232914,234050,0,ymax,"15h");
  DrawPeriod(runmin,runmax,235196,236866,0,ymax,"15i");
  DrawPeriod(runmin,runmax,236965,238621,0,ymax,"15j");
  DrawPeriod(runmin,runmax,238890,239025,0,ymax,"15k");
  DrawPeriod(runmin,runmax,239319,241544,0,ymax,"15l");
  DrawPeriod(runmin,runmax,244340,244628,0,ymax,"15n");
  DrawPeriod(runmin,runmax,244918,246994,0,ymax,"15o");

  // 2016
  DrawPeriod(runmin,runmax,247189,247193,0,ymax,"16a");
  DrawPeriod(runmin,runmax,247723,250585,0,ymax,"16b");
  DrawPeriod(runmin,runmax,250647,252091,0,ymax,"16c");
  
  DrawPeriod(runmin,runmax,252191,252518,0,ymax,"16d");
  DrawPeriod(runmin,runmax,252598,253609,0,ymax,"16e");
  DrawPeriod(runmin,runmax,253613,253979,0,ymax,"16f");
  DrawPeriod(runmin,runmax,253991,254332,0,ymax,"16g");
  DrawPeriod(runmin,runmax,254373,255469,0,ymax,"16h");
  DrawPeriod(runmin,runmax,255515,255650,0,ymax,"16i");
  DrawPeriod(runmin,runmax,256083,256420,0,ymax,"16j");
//.........这里部分代码省略.........

void simple2Dmodel_ex(Int_t statexp, Int_t statmod)
{
   //------------Suppress RooFit info messages-----------
   RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);

   //If RooWorkspace already generated
//    TFile *file = new TFile("./workspace_PARA_par1[0.0,10.0,10]_par2[-20.0,20.0,20]_OBS_obs1[-45.0,45.0,45]_obs2[-30.0,30.0,60].root");
//    RooWorkspace *ww = (RooWorkspace *) file->Get("init_workspace");


   //----------Init weights histogram--------
//    TFile *fileh = new TFile("./histo_par1_par2.root");
//    TH2F *histo = (TH2F*) fileh->Get("fittedParameters");
   //TH2F *histo = (TH2F*) fileh->Get("dummy");

   model = new BackTrack::Simple2DModel_Binned();
   model->InitWorkspace();
   model->InitParObs();

   //===========================
   //==     "EXP" DATASET     ==
   //===========================

   //Generate data with correlated gaussian-distribution parameters
   model->SetNumGen(statexp);
   RooDataSet data("data", "data to be fitted", model->GetObservables());
   RooDataSet params("params", "input parameters", model->GetParameters());
   RooRealVar& PAR1 = *(model->GetParameter("par1"));
   RooRealVar& PAR2 = *(model->GetParameter("par2"));

   for (int i = 0; i < model->GetNumGen(); i++) {
      //Double_t par1 = gRandom->Gaus(3.5,1.);
      Int_t par1     = (int) gRandom->Gaus(3.5, 1.);
      Double_t par2m = 15. - 0.5 * pow(par1, 2.);
      Double_t par2 = gRandom->Gaus(par2m, 1.);
      PAR1.setVal(par1);
      PAR2.setVal(par2);
      model->generateEvent(RooArgList(PAR1, PAR2), data);
      params.add(model->GetParameters());
   }

   //Create the "experimental" RooDataHist used for the fit  with the RooDataSet
   //generated just before
   RooDataHist datahist("datahist_exp", "input parameters", model->GetObservables(), data);

//    Modify the ranges for the fit
//    model->GetObservable("obs1")->setRange("RANGE",-10,32);
//    model->GetObservable("obs2")->setRange("RANGE",-20,15);

   //===========================
   //==    BACKTRACKING       ==
   //===========================
   //Initialise the BackTracking
   model->SetExperimentalDataHist(datahist);
   model->SaveInitWorkspace();
   model->SetNumInt(kTRUE);
   model->SetExtended(kFALSE);
   model->SetUniformInitWeight(kTRUE);
   model->SetNumGen(statmod);
   model->ImportAllModelData();
   model->ImportParamInitWeight();
   model->ConstructPseudoPDF(kFALSE);
   //model->fitTo(datahist, Save(), NumCPU(6), SumW2Error(kTRUE), PrintLevel(1), Minimizer("TMinuit","migrad"), Extended(kFALSE), Offset(kTRUE));
   model->fitTo(datahist, Save(), NumCPU(6), SumW2Error(kTRUE), PrintLevel(1), Minimizer("TMinuit", "migrad"), Extended(kFALSE), Offset(kTRUE), SetMaxIter(500000), SetMaxCalls(500000));
   //model->fitTo(datahist, Save(), NumCPU(6), SumW2Error(kTRUE), PrintLevel(1), Minimizer("GSLMultiMin","conjugatefr"), Extended(kFALSE), Offset(kTRUE), SetMaxCalls(500000), SetMaxIter(500000));
   //model->fitTo(datahist, Save(), NumCPU(4), SumW2Error(kTRUE), PrintLevel(1), Minimizer("GSLSimAn"), Extended(kFALSE), Offset(kTRUE)); //, SetEpsilon(0.01), SetMaxCalls(500000));


   //==================
   //==    PLOTS     ==
   //==================
   //---------------Results compared to experimental data----------------
   RooPlot* plot = model->GetObservable("obs1")->frame();
   data.plotOn(plot);
   model->GetPseudoPDF()->plotOn(plot);
   TCanvas* c = new TCanvas("Observables_results", "Observables_results");
   c->Divide(2, 1);
   c->cd(1);
   plot->Draw();

   plot = model->GetObservable("obs2")->frame();
   data.plotOn(plot);
   model->GetPseudoPDF()->plotOn(plot);
   c->cd(2);
   plot->Draw();

   //For binning
   RooRealVar* o1 = dynamic_cast<RooRealVar*>(model->GetObservables().at(0));
   RooAbsBinning& bins1 = o1->getBinning();
   Int_t N1 = bins1.numBins();

   RooRealVar* o2 = dynamic_cast<RooRealVar*>(model->GetObservables().at(1));
   RooAbsBinning& bins2 = o2->getBinning();
   Int_t N2 = bins2.numBins();


   TH2* hh = (TH2*)(model->GetPseudoPDF())->createHistogram("obs1,obs2", N1, N2);
   TH2* hhdata = (TH2*) data.createHistogram("obs1,obs2", N1, N2);

   if (hh && hhdata) {
//.........这里部分代码省略.........

void drawSupport(){

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


    char name[200];

    //sprintf(name,"Mydecorr_0_Calib1.root");
    sprintf(name,"Mydecorr_0_Calib1_1.root"); 
    //sprintf(name,"support.root");
    fin = TFile::Open(name);

    sprintf(name, "hcentrality");
    TH1* hcent = (TH1D*)fin->Get(name);

    sprintf(name, "hntrk_Fcal_%d",0);
    hntrk_Fcal[0] = (TH2*)fin->Get(name);
    hntrk_Fcal[0]->Rebin2D(5,5);

    sprintf(name, "hntrk_Fcal_%d",1);
    hntrk_Fcal[1] = (TH2*)fin->Get(name);
    sprintf(name, "hFCalAC");
    TH2* hFCalAC = (TH2*)fin->Get(name);

    TF1 *fa = new TF1("fa","[0] + [1]*x",0,7);
    double par2[] = {-1300,1125};
    fa->SetLineColor(kRed);
    fa->SetParameters(par2);
    fa->SetLineWidth(2);
    fa->SetLineStyle(2);

    TF1 *fa2 = new TF1("fa2","[0] + [1]*x",0,7);
    double par[] = {4.1,-1};
    fa2->SetLineColor(kRed);
    fa2->SetParameters(par);
    fa2->SetLineWidth(2);
    fa2->SetLineStyle(2);

    for (int ic=0; ic<NCENT; ic++) {
        sprintf(name,"hntrkloose_ic%d", ic); hntrkloose[ic] =(TH1D*)fin->Get(name);
        cout<<name<< " "<<hntrkloose[ic]->GetBinContent(3325)<< endl;
        sprintf(name,"hntrktight_ic%d", ic); hntrktight[ic] =(TH1D*)fin->Get(name);
    }

    for (int ic=0; ic<20; ic++) {
        sprintf(name,"hntrkcut_ic%d", ic); hntrkcut[ic] =(TH1D*)fin->Get(name);
    }



       sprintf(name,"sfcalAC");
       can = SetCan2(name, 2,2);
       hFCalAC->Draw("colz");
       hFCalAC->GetXaxis()->SetTitle("#Sigma E_{T}^{FCal A side} [TeV]");
       hFCalAC->GetYaxis()->SetTitle("#Sigma E_{T}^{FCal C side} [TeV]");
       hFCalAC->Rebin2D(10,10);
       hFCalAC->GetXaxis()->SetRangeUser(0,3.5);
       hFCalAC->GetYaxis()->SetRangeUser(0,3.5);
       gPad->SetLogz();
       fa2->Draw("same");
       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);


       sprintf(name,"sfcal");
       can = SetCan2(name, 2,2);
       hntrk_Fcal[0]->Draw("colz");
       hntrk_Fcal[0]->GetXaxis()->SetTitle("#Sigma E_{T}^{FCal} [TeV]");
       hntrk_Fcal[0]->GetYaxis()->SetTitle("N_{track}^{Loose}  [p_{T}>0.5GeV]");

       gPad->SetLogz();
       fa->Draw("same");
       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);



       sprintf(name,"sCentrality");
       can = SetCan2(name, 2,2);
       hcent->Draw();

       hcent->GetXaxis()->SetTitle("Centrality %");
       hcent->GetYaxis()->SetTitle("Nevents");

       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);



       int Nx=5, Ny=4;
       sprintf(name,"sntrkLoose");
       can = SetCan2D2(name,Nx, Ny);
       for (int ix=0; ix<Nx; ix++) {
       for (int iy=0; iy<Ny; iy++) {
       int ic = ix + iy*Nx;
       can->cd(ic+1);
       if(ic>=NCENT) continue;
       hntrkloose[ic] ->Rebin(10);
       hntrkloose[ic] ->DrawClone();
       gPad->SetLogy();
//.........这里部分代码省略.........