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

int
plotting_output()
{

	vector< int > x, y;
	vector< double > laplace;
	string mapping_tower_file = "./numerical_solution.txt";


	/* Stream to read table from file */
	ifstream istream_mapping;

	/* Open the datafile, if it won't open return an error */
	if (!istream_mapping.is_open())
	{
		istream_mapping.open( mapping_tower_file.c_str() );
		if(!istream_mapping)
		{
			cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << mapping_tower_file << endl;
			exit(1);
		}
	}

	string line_mapping;
	int x_i, y_i;
	double laplace_i;

	while ( getline( istream_mapping, line_mapping ) )
	{
		istringstream iss(line_mapping);
		iss >> x_i >> y_i >> laplace_i;
		x.push_back( x_i );
		y.push_back( y_i );
		laplace.push_back( laplace_i );
	}

	TCanvas *c1 = new TCanvas();
	TGraph2D *gr = new TGraph2D();
	for(int i = 0; i < x.size(); i++)
	{
		gr->SetPoint( gr->GetN(), x.at(i), y.at(i), laplace.at(i) );
	}
	gr->Draw("surf1");
	cout << "graph drawn " << endl;

	TCanvas *c1 = new TCanvas();
	gr->Draw("P");





	return 0;
}

TCanvas* graph2dfit_test()
{
	gStyle->SetOptStat(0);
	gStyle->SetOptFit();

	TCanvas *c = new TCanvas("c", "Graph2D example", 0, 0, 600, 800);
	
	Double_t rnd, x, y, z;
	Double_t e = 0.3;
	Int_t nd = 400;
	Int_t np = 10000;

	TRandom r;
	Double_t fl = 6;
	TF2  *f2 = new TF2("f2", "1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200", -fl, fl, -fl, fl);
	f2->SetParameters(1, 1);
	TGraph2D *dt = new TGraph2D();

	// Fill the 2D graph
	Double_t zmax = 0;
	for (Int_t N = 0; N<nd; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		if (z>zmax) zmax = z;
		dt->SetPoint(N, x, y, z);
	}


	f2->SetParameters(0.5, 1.5);
	dt->Fit(f2);
	TF2 *fit2 = (TF2*)dt->FindObject("f2");

	f2->SetParameters(1, 1);

	for (Int_t N = 0; N<np; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		h1->Fill(f2->Eval(x, y) - z);
		z = dt->Interpolate(x, y);
		h2->Fill(f2->Eval(x, y) - z);
		z = fit2->Eval(x, y);
		h3->Fill(f2->Eval(x, y) - z);
	}

	gStyle->SetPalette(1);
	
	f2->SetTitle("Original function with Graph2D points on top");
	f2->SetMaximum(zmax);
	gStyle->SetHistTopMargin(0);
	f2->Draw("surf1");
	dt->Draw("same p0");



	return c;
}

  bool PlotBundle::drawHistsAndGraphs(const PlotStyle &plot_style) const {
    if (hasDrawables()) {
      std::vector<DrawableDataObjectDrawOptionPair<TH1*> >::const_iterator hist_it;

      for (hist_it = histograms.begin(); hist_it != histograms.end();
          hist_it++) {
        TH1* hist = hist_it->data_object;
        std::string draw_option(hist_it->draw_option);
        if (hist) {
          if (hist->GetDimension() > 1) {
            hist->GetZaxis()->SetLimits(plot_axis.z_axis_range.low,
                plot_axis.z_axis_range.high);
            hist->GetZaxis()->SetRangeUser(plot_axis.z_axis_range.low,
                plot_axis.z_axis_range.high);
          }
          draw_option.append("SAME");
          hist->Draw(draw_option.c_str());
        }
      }

      std::vector<DrawableDataObjectDrawOptionPair<TGraph*> >::const_iterator graph_it;

      for (graph_it = graphs.begin(); graph_it != graphs.end(); graph_it++) {
        TGraph* graph = graph_it->data_object;
        std::string draw_option(graph_it->draw_option);
        if (graph) {
          draw_option.append("SAME");
          graph->Draw(draw_option.c_str());
        }
      }

      std::vector<DrawableDataObjectDrawOptionPair<TGraph2D*> >::const_iterator graph2d_it;

      for (graph2d_it = graphs2d.begin(); graph2d_it != graphs2d.end();
          graph2d_it++) {
        TGraph2D* graph = graph2d_it->data_object;
        std::string draw_option(graph2d_it->draw_option);
        if (graph) {
          draw_option.append("SAME");
          graph->Draw(draw_option.c_str());
        }
      }
      return true;
    } else {
      std::cout
          << "Dude, you forgot to add a histogram or graph in the plot bundle..."
          << " Please add at least one and make sure it points to an existing object!"
          << std::endl;
      return false;
    }
  }

//g++ allDataPrint.cpp `root-config --cflags --glibs`
int main(int argc, char** argv)
#endif
{
  TCanvas c3("c3","Grafico",640,512);
  TCanvas c1("c1","Confronto",1280,512);
  c1.Divide(2,1);
  
  preparedraw myData (argv[1],
		      // preparedraw::doMax |
		      preparedraw::doFh |
		      preparedraw::doSh |
		      preparedraw::doErr);
  
  TGraph2D *g = myData.data();
  TGraph *gb = myData.firsthalf();//before
  TGraph *ga = myData.secondhalf();//after
  TGraph *gerrs = myData.errs();
  //  TGraph *maxs = myData.maximum();
  
  c3.cd();
  g->GetXaxis()->SetTitle("X");
  g->GetYaxis()->SetTitle("T");
  //g->Draw("cont1");
  g->Draw("pcol");
  //g->Draw();
  //grafo.Draw("surf1");

  cout<<"Disegno i grafici\n";
  TMultiGraph *mg = new TMultiGraph("integrali","Integrali prima e dopo la barriera");
  
  ga->SetLineColor(2);
  mg->Add(gb);
  mg->Add(ga);
  c1.cd(1);
  mg->Draw("apl");
    

  gerrs->SetTitle("Andamento degli errori");
  c1.cd(2);
  gerrs->Draw("apl");
#ifndef __CINT__
  theApp.Run(true);
  return 0;
#endif
}

//g++ allDataPrint.cpp `root-config --cflags --glibs`
int main(int argv, char** argc)
#endif
{
  /*  TCanvas c("c","Real",640,512);
  TGraph2D *grafo = new TGraph2D("out.txt");
  grafo->Draw("pcol");
  TCanvas c2("c2","Imaginary",640,512);
  TGraph2D *Cgrafo = new TGraph2D("iout.txt");
  Cgrafo->Draw("pcol");*/
  TCanvas c3("c3","Norm",640,512);
  TGraph2D *Ngrafo = new TGraph2D(argc[1]);
  Ngrafo->Draw("pcol");
  //grafo.Draw("surf1");
#ifndef __CINT__
  theApp.Run(true);
  return 0;
#endif
}

int w1() {

	TCanvas *c = new TCanvas("c1", "", 0, 0, 700, 600);
	TGraph2D *dt = new TGraph2D();

	ifstream in("/opt/workspace-cpp/simulation/result_rand3.csv");

	string buff;

	getline(in, buff);

	int meanResidenceTime;
	int products;
	float lastProductPrice;
	float shopIncome;

	int n = 0;
	while (getline(in, buff)) {
			for (unsigned int i = 0; i < buff.size(); i++) {
				if (buff[i] == ';') {
					buff[i] = ' ';
				}
			}
			stringstream ss(buff);
			ss >> meanResidenceTime;
			ss >> products;
			ss >> lastProductPrice;
			ss >> shopIncome;
			dt->SetPoint(n, products, meanResidenceTime, shopIncome);
			n++;
	}
	gStyle->SetPalette(1);
	dt->SetTitle("Zaleznosc zysku od liczby produktow oraz czasu przebywania w sklepie");
	dt->GetXaxis()->SetTitle("Liczba produktow");
	dt->GetYaxis()->SetTitle("Sredni czas przebywania w sklepie");
	dt->GetZaxis()->SetTitle("Dochod");
	dt->Draw("surf1");
}

TH2D grMu(TTree *tree, double quantileExpected){

	TGraph2D *gr = new TGraph2D();
	float X, Y, qE;
	double limit;

	std::string x = "trackedParam_proc_scaling_muV";
	std::string y = "trackedParam_proc_scaling_muF";

	tree->SetBranchAddress(Form("%s",x.c_str()),&X);
	tree->SetBranchAddress(Form("%s",y.c_str()),&Y);
	tree->SetBranchAddress("quantileExpected",&qE);
	tree->SetBranchAddress("limit",&limit);

	int pt = 0;
	for (int i=0;i<tree->GetEntries();i++){
	  tree->GetEntry(i);
	  if (TMath::Abs(qE-quantileExpected)>0.001) continue;
	  gr->SetPoint(pt,X,Y,limit); pt++;
	  //std::cout << " Lim " << X << ", " << Y << ", " << limit << std::endl;
	}
	gr->Draw("colz");
	TH2D * expected = (TH2D*) gr->GetHistogram();
	expected->SetTitle("");
	expected->GetZaxis()->SetTitleOffset(1.2);
	expected->GetYaxis()->SetTitle("#it{#mu}_{ggH}");
	expected->GetXaxis()->SetTitle("#it{#mu}_{qqH,VH}");
	expected->GetZaxis()->SetTitle("B(H #rightarrow inv.) - 95% CL upper limit");
	expected->GetXaxis()->SetRangeUser(MUVMIN,MUVMAX);
	expected->GetYaxis()->SetRangeUser(MUFMIN,MUFMAX);
	//expected->GetXaxis()->SetLimits(MUVMIN,MUVMAX);
	//expected->GetYaxis()->SetLimits(MUFMIN,MUFMAX);

	expected->SetMaximum(0.6);
	expected->SetMinimum(0.05);
	return *expected;
}

Int_t line3Dfit()
{
   gStyle->SetOptStat(0);
   gStyle->SetOptFit();


   //double e = 0.1;
   Int_t nd = 10000;


//    double xmin = 0; double ymin = 0;
//    double xmax = 10; double ymax = 10;

   TGraph2D * gr = new TGraph2D();

   // Fill the 2D graph
   double p0[4] = {10,20,1,2};

   // generate graph with the 3d points
   for (Int_t N=0; N<nd; N++) {
      double x,y,z = 0;
      // Generate a random number 
      double t = gRandom->Uniform(0,10);
      line(t,p0,x,y,z);
      double err = 1;
    // do a gaussian smearing around the points in all coordinates
      x += gRandom->Gaus(0,err);  
      y += gRandom->Gaus(0,err);  
      z += gRandom->Gaus(0,err);  
      gr->SetPoint(N,x,y,z);
      //dt->SetPointError(N,0,0,err);
   }
   // fit the graph now 
   
   ROOT::Fit::Fitter  fitter;
   
   
   // make the functor objet
   SumDistance2 sdist(gr);
#ifdef __CINT__
   ROOT::Math::Functor fcn(&sdist,4,"SumDistance2");
#else
   ROOT::Math::Functor fcn(sdist,4);
#endif
   // set the function and the initial parameter values
   double pStart[4] = {1,1,1,1};
   fitter.SetFCN(fcn,pStart);
   // set step sizes different than default ones (0.3 times parameter values)
   for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);
  
   bool ok = fitter.FitFCN();
   if (!ok) {
      Error("line3Dfit","Line3D Fit failed");
      return 1;
   }
   
   const ROOT::Fit::FitResult & result = fitter.Result();
   
   std::cout << "Total final distance square " << result.MinFcnValue() << std::endl;
   result.Print(std::cout);
   

   gr->Draw("p0");

   // get fit parameters
   const double * parFit = result.GetParams();

   
   // draw the fitted line
   int n = 1000;
   double t0 = 0;
   double dt = 10;
   TPolyLine3D *l = new TPolyLine3D(n);
   for (int i = 0; i <n;++i) {
      double t = t0+ dt*i/n;
      double x,y,z;
      line(t,parFit,x,y,z);
      l->SetPoint(i,x,y,z);
   }
   l->SetLineColor(kRed);
   l->Draw("same");
   
   // draw original line
   TPolyLine3D *l0 = new TPolyLine3D(n);
   for (int i = 0; i <n;++i) {
      double t = t0+ dt*i/n;
      double x,y,z;
      line(t,p0,x,y,z);
      l0->SetPoint(i,x,y,z);
   }
   l0->SetLineColor(kBlue);
   l0->Draw("same");
   return 0;
}

//.........这里部分代码省略.........
        laplace[1][0] = v_left / 2;
        laplace[2][0] = v_left / 2;
        laplace[3][0] = v_left / 2;
        laplace[1][40] = v_left;
        laplace[2][40] = v_left / 2;
        laplace[3][40] = v_left / 2;
        laplace[4][40] = v_left / 4;
        cout << "    (additional Constraints applied)" << endl;
    }

    //Write out the initial matrix
    if(write_matrix_init == true)
    {
        for(int i = 0; i < x_max+2; i++)
        {
            cout << " | ";
            for(int j = 0; j < y_max+2; j++)
            {
                cout << laplace[i][j] << " | ";
            }
            cout << endl;
        }
        cout << endl << endl;
    }


    //Iterate over the inner rows/columns of the matrix, averaging at each point.
    for ( int k = 0; k < iterations; k++)
    {
        //Loop over all iterations, averaging each point with those around them... skipping over boundary conditions
        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                laplace[i][j] = 0.25*( laplace[i-1][j] + laplace[i+1][j] + laplace[i][j-1] + laplace[i][j+1] );
//				laplace[i][j] = 0.125*( laplace[i-1][j] + laplace[i+1][j] + laplace[i][j-1] + laplace[i][j+1] + laplace[i-1][j-1] + laplace[i-1][j+1] + laplace[i+1][j-1] + laplace [i+1][j+1] );
            }
        }
    }
    cout << "Computations Completed. Plotting..." << endl;

    //Write out the final matrix
    if(write_matrix_final == true)
    {
        for(int i = 0; i < x_max+2; i++)
        {
            cout << " | ";
            for(int j = 0; j < y_max+2; j++)
            {
                cout << laplace[i][j] << " | ";
            }
            cout << endl;
        }
        cout << endl << endl;
    }

    //Write the output to a file
    if(output)
    {
        ofstream myfile;
        myfile.open("numerical_solution.txt");

        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                myfile << i << " " << j << " " << laplace[i][j] << endl;
//				myfile << j << " " << i << " " << laplace[i][j] << endl;
            }
        }
        myfile.close();
        cout << "Output file created." << endl;
    }


    //Plot the final output
    if(plot)
    {
        //Create the output graph
        TCanvas *c1 = new TCanvas();
        TGraph2D *gr = new TGraph2D();
        for(int i = 1; i < x_max+1; i++)
        {
            for(int j = 1; j < y_max+1; j++)
            {
                gr->SetPoint( gr->GetN(), i, j, laplace[i][j] );
            }
        }
        //gr->Draw("surf1");
        gr->Draw("COLZ");

        TCanvas *c2 = new TCanvas();
        gr->Draw("surf1");

        //	TH2 *h2 = gr->Project("xy");
        //	h2->Draw();
    }

    return 0;
}

void fillGraphsFromFilesDeltaNLL( const TString& par1name,
				  const TString& par2name,
				  const vector<TString>& fnames,
				  vector<string>&  keys,
				  map<string,TGraph2D *>& m_graphs)
{

  std::cout << "fillGraphsFromFilesDeltaNLL 1" << std::endl;

  keys.push_back("exp68");
  keys.push_back("exp95");
  keys.push_back("exp99");

  // uncommented below to plot observed!
  keys.push_back("obs95");

  TGraph2D *grobs = new TGraph2D();
  TGraph2D *grexp = new TGraph2D();

  m_graphs["obs95"] = grobs;
  m_graphs["exp95"] = grexp;

  grobs->SetName("graph2Dobs95");
  grexp->SetName("graph2Dexp95");

  Int_t nobs=0, nexp=0;

  for( size_t i=0; i<fnames.size(); i++) {
    
    TFile *f = new TFile(fnames[i]);
    TTree *t = (TTree *) f->Get("limit");

    if (!t) { 
      std::cerr<<"TFile "<<f->GetName()<<" does not contain the tree"<<std::endl;
      return;
    }
    cout << fnames[i] << " has limit tree with " << t->GetEntries() << " entries." << endl;

    Float_t deltaNLL, par1, par2;
    Int_t iToy;

    t->SetBranchAddress("iToy", &iToy);
    t->SetBranchAddress("deltaNLL", &deltaNLL);
    t->SetBranchAddress(par1name, &par1);
    t->SetBranchAddress(par2name, &par2);

    for (size_t j = 0, n = t->GetEntries(); j < n; ++j) {
      t->GetEntry(j);
      printf ("%d\r",j);
      if( !iToy){
	//	cout <<"!iToy" << endl;
	grobs->SetPoint(nobs++,par1,par2,2*deltaNLL);
	//	cout <<"grobs->SetPoint("<<nobs++<<","<<par1<<","<< par2<< ","<< 2*deltaNLL << endl;
      }
      else if (iToy == -1) {
	//	cout <<"iToy == -1" << endl;
	grexp->SetPoint(nexp++,par1,par2,2*deltaNLL);
      }
      else {
	cerr << "Unexpected value for iToy, = " << iToy << endl;
	exit(-1);
      }
    } // tree entry loop

    f->Close();
    delete f;

  } // file loop
  cout << endl;

  m_graphs["exp68"] = (TGraph2D*)grexp->Clone("graph2Dexp68");
  m_graphs["exp99"] = (TGraph2D*)grexp->Clone("graph2Dexp99");

#if 0
  TCanvas *canv = new TCanvas("tester","tester",500,500);
  cout << grexp->GetN()<<" points. " <<endl;
  grexp->Draw("TRI"); // cont 5z list");
#endif
}                                         // fillGraphsFromFilesDeltaNLL

  PadBoundaries PlotBundle::calculatePadBoundaries() const {
    PadBoundaries union_pad_boundaries;
    bool nothing_drawn(true);

    std::vector<DrawableDataObjectDrawOptionPair<TH1*> >::const_iterator hist_it;

    for (hist_it = histograms.begin(); hist_it != histograms.end(); hist_it++) {
      TH1* hist = hist_it->data_object;
      std::string draw_option(hist_it->draw_option);
      if (hist) {
        hist->Draw(draw_option.c_str());
        if (plot_axis.x_axis_range.active)
          hist->GetXaxis()->SetRangeUser(plot_axis.x_axis_range.low,
              plot_axis.x_axis_range.high);
        PadBoundaries current_pad_boundaries(*hist);

        if (nothing_drawn) {
          union_pad_boundaries = current_pad_boundaries;
          nothing_drawn = false;
        } else {
          union_pad_boundaries = union_pad_boundaries.getUnionBoundaries(
              current_pad_boundaries);
        }
      }
    }

    std::vector<DrawableDataObjectDrawOptionPair<TGraph*> >::const_iterator graph_it;

    for (graph_it = graphs.begin(); graph_it != graphs.end(); graph_it++) {
      TGraph* graph = graph_it->data_object;
      std::string draw_option("A");
      draw_option.append(graph_it->draw_option);
      if (graph) {
        graph->Draw(draw_option.c_str());
        if (plot_axis.x_axis_range.active)
          graph->GetXaxis()->SetRangeUser(plot_axis.x_axis_range.low,
              plot_axis.x_axis_range.high);
        PadBoundaries current_pad_boundaries(*graph->GetHistogram());
        if (nothing_drawn) {
          union_pad_boundaries = current_pad_boundaries;
          nothing_drawn = false;
        } else {
          union_pad_boundaries = union_pad_boundaries.getUnionBoundaries(
              current_pad_boundaries);
        }
      }
    }

    std::vector<DrawableDataObjectDrawOptionPair<TGraph2D*> >::const_iterator graph2d_it;

    for (graph2d_it = graphs2d.begin(); graph2d_it != graphs2d.end();
        graph2d_it++) {
      TGraph2D* graph = graph2d_it->data_object;
      std::string draw_option("A");
      draw_option.append(graph2d_it->draw_option);
      if (graph) {
        graph->Draw(draw_option.c_str());
        PadBoundaries current_pad_boundaries(*(TH1*) graph->GetHistogram());
        current_pad_boundaries.is_2d = true;
        current_pad_boundaries.z_min = graph->GetZmin();
        current_pad_boundaries.z_max = graph->GetZmax();
        if (nothing_drawn) {
          union_pad_boundaries = current_pad_boundaries;
          nothing_drawn = false;
        } else {
          union_pad_boundaries = union_pad_boundaries.getUnionBoundaries(
              current_pad_boundaries);
        }
      }
    }

    return union_pad_boundaries;
  }

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


			// 3D FITTING CODE (based on line3Dfit.C), draw graph and line fit
			ROOT::Fit::Fitter  fitter;
		   	SumDistance2 sdist(graph);
#ifdef __CINT__
		   	ROOT::Math::Functor fcn(&sdist,4,"SumDistance2");
#else
		   	ROOT::Math::Functor fcn(sdist,4);
#endif
			// set the function and the initial parameter values
			double pStart[4] = {param0,param1,param2,param3};
			fitter.SetFCN(fcn,pStart);
			// set step sizes different than default ones (0.3 times parameter values)
			for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);

			bool ok = fitter.FitFCN();
			if (!ok) {
			  Error("line3Dfit","Line3D Fit failed");
			  fitFailed = true;
			} else {
				const ROOT::Fit::FitResult & result = fitter.Result();
				const double * fitParams = result.GetParams();

				sumSquares = result.MinFcnValue();
				std::cout << "Sum of distance squares:  " << sumSquares << std::endl;
				std::cout << "Sum of distance squares divided by numEntries: " << sumSquares/numEntries << std::endl;
				std::cout << "Theta : " << TMath::ATan(sqrt(pow(fitParams[1], 2) + pow(fitParams[3], 2))) << std::endl;
				// result.Print(std::cout); // (un)suppress results output

				// Draw the graph
				graph->SetMarkerStyle(8);
				graph->SetMarkerSize(0.5);
				graph->Draw("pcol");

				// Draw the fitted line
				int n = 1000;
				double t0 = 0; // t is the z coordinate
				double dt = 40.96;
				TPolyLine3D *l = new TPolyLine3D(n);
				for (int i = 0; i <n;++i) {
				  double t = t0+ dt*i/n;
				  double x,y,z;
				  line(t,fitParams,x,y,z);
				  l->SetPoint(i,x,y,z);
				}
				l->SetLineColor(kRed);
				l->Draw("same");

				// Access fit params and minfcnvalue
				// cout << "FIT1: " << fitParams[1] << "\n";
				// cout << "FIT2: " << result.MinFcnValue() << "\n";
			}

			// Criteria to be a good event (if not good entry, then don't show)
			bool isGoodEvent = false;

				// the following block of code finds the mean X, Y ans Z values
				double meanX = 0;
				double meanY = 0;
				double meanZ = 0;
				reader->SetEntry(startEntryNum);
				for (int i = 0; i < endEntryNum - startEntryNum; i++) {
					meanX += graph->GetX()[i];
					meanY += graph->GetY()[i];
					meanZ += graph->GetZ()[i];

void kdTreeBinning() {

   // -----------------------------------------------------------------------------------------------
   //  C r e a t e  r a n d o m  s a m p l e  w i t h  r e g u l a r  b i n n i n g  p l o t t i n g
   // -----------------------------------------------------------------------------------------------

   const UInt_t DATASZ = 10000;
   const UInt_t DATADIM = 2;
   const UInt_t NBINS = 50;

   Double_t smp[DATASZ * DATADIM];

   double mu[2] = {0,2};
   double sig[2] = {2,3};
   TRandom3 r;
   r.SetSeed(1);
   for (UInt_t i = 0; i < DATADIM; ++i)
      for (UInt_t j = 0; j < DATASZ; ++j)
         smp[DATASZ * i + j] = r.Gaus(mu[i], sig[i]);

   UInt_t h1bins = (UInt_t) sqrt(NBINS);

   TH2D* h1 = new TH2D("h1BinTest", "Regular binning", h1bins, -5., 5., h1bins, -5., 5.);
   for (UInt_t j = 0; j < DATASZ; ++j)
      h1->Fill(smp[j], smp[DATASZ + j]);


   // ---------------------------------------------------------------------------------------------
   // C r e a t e  K D T r e e B i n n i n g  o b j e c t  w i t h  T H 2 P o l y  p l o t t i n g
   // ---------------------------------------------------------------------------------------------

   TKDTreeBinning* kdBins = new TKDTreeBinning(DATASZ, DATADIM, smp, NBINS);

   UInt_t nbins = kdBins->GetNBins();
   UInt_t dim   = kdBins->GetDim();

   const Double_t* binsMinEdges = kdBins->GetBinsMinEdges();
   const Double_t* binsMaxEdges = kdBins->GetBinsMaxEdges();

   TH2Poly* h2pol = new TH2Poly("h2PolyBinTest", "KDTree binning", kdBins->GetDataMin(0), kdBins->GetDataMax(0), kdBins->GetDataMin(1), kdBins->GetDataMax(1));

   for (UInt_t i = 0; i < nbins; ++i) {
      UInt_t edgeDim = i * dim;
      h2pol->AddBin(binsMinEdges[edgeDim], binsMinEdges[edgeDim + 1], binsMaxEdges[edgeDim], binsMaxEdges[edgeDim + 1]);
   }

   for (UInt_t i = 1; i <= kdBins->GetNBins(); ++i)
      h2pol->SetBinContent(i, kdBins->GetBinDensity(i - 1));

   std::cout << "Bin with minimum density: " << kdBins->GetBinMinDensity() << std::endl;
   std::cout << "Bin with maximum density: " << kdBins->GetBinMaxDensity() << std::endl;

   TCanvas* c1 = new TCanvas("glc1", "TH2Poly from a kdTree",0,0,600,800);
   c1->Divide(1,3);
   c1->cd(1);
   h1->Draw("lego");

   c1->cd(2);
   h2pol->Draw("COLZ L");
   c1->Update();


   /* Draw an equivalent plot showing the data points */
   /*-------------------------------------------------*/

   std::vector<Double_t> z = std::vector<Double_t>(DATASZ, 0.);
   for (UInt_t i = 0; i < DATASZ; ++i)
      z[i] = (Double_t) h2pol->GetBinContent(h2pol->FindBin(smp[i], smp[DATASZ + i]));

   TGraph2D *g = new TGraph2D(DATASZ, smp, &smp[DATASZ], &z[0]);
   gStyle->SetPalette(1);
   g->SetMarkerStyle(20);

   c1->cd(3);
   g->Draw("pcol");
   c1->Update();

   // ---------------------------------------------------------
   // make a new TH2Poly where bins are ordered by the density
   // ---------------------------------------------------------

   TH2Poly* h2polrebin = new TH2Poly("h2PolyBinTest", "KDTree binning", kdBins->GetDataMin(0), kdBins->GetDataMax(0), kdBins->GetDataMin(1), kdBins->GetDataMax(1));
   h2polrebin->SetFloat();

   /*---------------------------------*/
   /* Sort the bins by their density  */
   /*---------------------------------*/

   kdBins->SortBinsByDensity();

   for (UInt_t i = 0; i < kdBins->GetNBins(); ++i) {
      const Double_t* binMinEdges = kdBins->GetBinMinEdges(i);
      const Double_t* binMaxEdges = kdBins->GetBinMaxEdges(i);
      h2polrebin->AddBin(binMinEdges[0], binMinEdges[1], binMaxEdges[0], binMaxEdges[1]);
   }

   for (UInt_t i = 1; i <= kdBins->GetNBins(); ++i){
      h2polrebin->SetBinContent(i, kdBins->GetBinDensity(i - 1));}

   std::cout << "Bin with minimum density: " << kdBins->GetBinMinDensity() << std::endl;
//.........这里部分代码省略.........

void PlotSVMopt(){

 //  // TString FileData = "SVMoptData.dat";
//   TFile *f = new TFile("SVMoptTree.root","RECREATE");
//   TTree *T = new TTree("SVMoptTree","data from ascii file");
//   Long64_t nlines = T->ReadFile("SVMoptData.dat","Gamma:C:ROCint:SigAt1Bkg_test:SigAt1Bkg_train");
//   printf(" found %lld pointsn",nlines);
//   T->Write();

// TGraph2D *grROCint = new TGraph2D();
//  for(int i=0;i<nlines;i++){

//    grROCint->SetPoint(i,);
//  }

  TGraph2D *grROCint = new TGraph2D();
  

TGraph2D *grSigAt1Bkg_test = new TGraph2D();
TGraph2D *grOverTrain = new TGraph2D();

 TString dir = gSystem->UnixPathName(__FILE__);
 dir.ReplaceAll("PlotSVMopt.C","");
 dir.ReplaceAll("/./","/");
 ifstream in;
 // in.open(Form("%sSVMoptData.dat",dir.Data()));
 in.open(Form("%sSVMoptData_NEW2.dat",dir.Data()));

 Float_t Gamma,C,ROCint, SigAt1Bkg_test, SigAt1Bkg_train;
 Int_t nlines = 0;
 while (1) {
   in >> Gamma >> C >> ROCint >> SigAt1Bkg_test >> SigAt1Bkg_train;
   if (!in.good()) break;
   // if (in.good()){
      if (nlines < 5) printf("Gamma=%8f, C=%8f, ROCint=%8f\n",Gamma, C, ROCint);
      grROCint->SetPoint(nlines,Gamma,C,ROCint);
      grSigAt1Bkg_test->SetPoint(nlines,Gamma,C,SigAt1Bkg_test);
      grOverTrain->SetPoint(nlines,Gamma,C,fabs(SigAt1Bkg_train-SigAt1Bkg_test));
      nlines++;
      //   }
   }
   printf(" found %d points\n",nlines);
   in.close();
   TCanvas *cSVMopt = new TCanvas("cSVMopt","SVM model choice",600,600);
   cSVMopt->Divide(2,2);
   cSVMopt->cd(1);
   cSVMopt_1->SetLogx();
   cSVMopt_1->SetLogy();
   grROCint->GetXaxis()->SetLabelSize(0.04);
   grROCint->GetYaxis()->SetLabelSize(0.04);
   grROCint->GetZaxis()->SetLabelSize(0.04);
   grROCint->GetXaxis()->SetTitle("#gamma");
   grROCint->GetYaxis()->SetTitle("C");
   grROCint->GetZaxis()->SetTitle("ROC integral");
   grROCint->SetTitle("ROC integral");
   grROCint->SetMaximum(1.0);
   grROCint->SetMinimum(0.9);
   grROCint->Draw("COLZ");
   cSVMopt->cd(2);
   cSVMopt_2->SetLogx();
   cSVMopt_2->SetLogy();
   grSigAt1Bkg_test->SetTitle("Signal at 1% Bkg level (test)");
   //   grSigAt1Bkg_test->Draw("surf1");
   grSigAt1Bkg_test->Draw("COLZ");
   cSVMopt->cd(3);
   cSVMopt_3->SetLogx();
   cSVMopt_3->SetLogy();
   grOverTrain->SetTitle("Overtraining");
   grOverTrain->Draw("COLZ");
   //cSVMopt->SaveAs("SVMoptC1.root");
}

//.........这里部分代码省略.........
			std::cout << "Simpson integral: " <<s_integral <<std::endl;
			integral_hist->Fill(s_integral);
			cmp_int[i] = s_integral;
			Int_t lines = (Int_t)intb.size();
			TGraph *r_integral = new TGraph(lines, _inta, _intb);
			
			std::cout << "ROOT integral: " << r_integral->Integral() << std::endl;
			cmp_int_root[i] = r_integral->Integral();
			
				//expanding
				//expand(y, THRS_EXPAND, RATIO_EXPAND, LINES);
			
			
				//Filling TGraph2D
			
			for(Int_t j = 0; j <LINES ; j++){
				if (y[j] > max){
					max = y[j];
					maxwl = x[j];
				}
				gr->SetPoint(j+i*LINES, x[j],i,y[j]);
			}
			
			
			in.seekg(0, std::ios::beg);
			in.close();
			
				//Plotting each spectrum
			
			TGraph *_gr = new TGraph(LINES,x,y);
			_gr->GetHistogram()->GetXaxis()->SetTitle("#lambda in nm");
			_gr->GetHistogram()->GetYaxis()->SetTitle("Intensity in dB");
			c1->cd(i-NUM_ARGS);
			_gr->Draw("AP");
			_gr->GetYaxis()->SetRangeUser(-80.,-10.);
			_gr->GetXaxis()->SetRangeUser(startwl,stopwl);
			_gr->SetTitle(tmp.c_str());
			c1->Update();
			
			
				//Calculating asymmetry
			std::cout << "maximum: " << max << std::endl;
			double leftlimit, rightlimit = 1;
			leftlimit = findlower(x,y, max);
			rightlimit = findupper(x,y, max);
			if (leftlimit != 1 && rightlimit != 1){
				width_ary[i] = (leftlimit +rightlimit)/2;
			}else{
				width_ary[i] = maxwl;
			}
			double calced_asy = (maxwl-leftlimit)/(rightlimit-maxwl);
			asymmety_ary[i-NUM_ARGS] = calced_asy;
			
			std::cout << "Asymmetry: " << calced_asy << std::endl;
			
		}catch(std::exception e){
			std::cout << e.what()<< std::endl;
		}
	}
	
	
		//Setting style for 3D Plot
	TCanvas *d = new TCanvas("big","big",10,10,1500,800);
	d->Divide(2,2);
	d->cd(1);
	TGraph *the_ints = new TGraph(argc-1,argc_ary,cmp_int);