C++ TGraph2D类代码示例

void setGraphOptions(TGraph2D &g)
{
  g.SetTitle("");
  g.SetMarkerColor(1);
  g.SetMarkerStyle(24);
  g.SetMarkerSize(.5);
}

TGraph2D* preparedraw::data(){
  TGraph2D *toreturn = new TGraph2D(Z.size());
  for(int j=0;j<Nt;j++){
    for(int i=0;i<Nl;i++)
      toreturn->SetPoint(i+j*Nl,sstep*sskip*i,times.at(j),abs(Z.at(i+j*Nl)));
  }
  return toreturn;
}

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

  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
}

TGraph2D* GetInterpolatingGraph(TH2F *hold){
  float binsize = hold->GetXaxis()->GetBinWidth(1)/2.;
  TString name = hold->GetName();
  name.ReplaceAll("Org","");
  TGraph2D *g = new TGraph2D(hold);
  g->SetNpx(int(g->GetXmax()-g->GetXmin())/binsize);
  g->SetNpy(int(g->GetYmax()-g->GetYmin())/binsize);
  //cout << "name " << g->GetN() << " " << hold->Integral() << endl;
  return g;
}

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

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

TH2F* InterpolateThisHistogram(TH2F *hold/*, TH2F* hnew*/){
  float binsize = hold->GetXaxis()->GetBinWidth(1)/2.;
  TString name = hold->GetName();
  name.ReplaceAll("Org","");
  TGraph2D *g = new TGraph2D(hold);
  //cout << g->GetXmax() << " " << g->GetXmin() << " " << g->GetYmax() << " " << g->GetYmin() << " " << binsize << endl;
  g->SetNpx(int(g->GetXmax()-g->GetXmin())/binsize);
  g->SetNpy(int(g->GetYmax()-g->GetYmin())/binsize);
  TH2F *hnew = (TH2F*)g->GetHistogram();
  //TH2F *htemp = (TH2F*)hnew->Clone(name);
  //name.ReplaceAll("YXZ","");
  TH2F *h = new TH2F(name.Data(),hold->GetTitle(),hnew->GetNbinsX(),g->GetXmin()-binsize,g->GetXmax()-binsize,hnew->GetNbinsY(),g->GetYmin()-binsize,g->GetYmax()-binsize);
  for(unsigned int x = 1; x<=hnew->GetNbinsX(); ++x){
    for(unsigned int y = 1; y<=hnew->GetNbinsY(); ++y){
      h->SetBinContent(x,y,hnew->GetBinContent(x,y));
    }
  }
  delete g;
  return h;
}

pair<float, float> Analyzer::MinG(TGraph2D *g,double *min0,double*min1){
	pair<float,float> R(-99,-99);
	if(g==NULL) return R;
	if(g->GetN()==0)return R;
	double *x1,*y1,*z1;
	x1=g->GetX();
	y1=g->GetY();
	z1=g->GetZ();
	
	float a=z1[0];R=pair<float,float>(x1[0],y1[0]);
	for(int i=0;i<g->GetN();i++){if((z1[i]<a)||(a<0)){a=z1[i];    R=pair<float,float>(x1[i],y1[i]); }}
	
	if(min0!=NULL) *min0=a;
	if(min1!=NULL) *min1=z1[g2->GetN()-1]	;
	return R;
}

void Analyzer::ComputeMin(){
	g2=new TGraph2D(); //TODO
	
	alpha=1.0;beta=0;
	Loop(t_data,1);
	if(varName=="QGLMLP")
		Loop(t_mc,4);
	//scan
	alpha=1.0;beta=0;
	for(float ai=0.7; ai<=1.1; ai+=0.02)
		{
		Reset(h_mc);	
		alpha=ai;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	alpha=1.0;beta=0;
	for(float bi=-0.5; bi<=0.5; bi+=0.01)
		{
		Reset(h_mc);	
		beta=bi;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );	
		}
	
	//Find min0;min1
	float min0=1,min1=0;
	pair<float,float> R=MinG(g2);
	min0=R.first;min1=R.second;

	for(int i=-nstep;i<=nstep;i++)
	for(int j=-nstep;j<=nstep;j++)
        	{
        	alpha=min0+i*stp0;
        	beta=min1+j*stp1;
		Loop(t_mc,2);
		h_mc->Scale(h_data->Integral()/h_mc->Integral());
		g2->SetPoint(g2->GetN(),alpha,beta, h_data->Chi2Test(h_mc,opt.c_str())  );
		}
	//double min0,min1;	
	R=MinG(g2);
	printf("a=%.3f;b=%.3f;lmin=%.3f;lmax=%.3f;break;\n",R.first,R.second,lmin,lmax);
	return;
	
}

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

void view_SMEvents_3D_from_Hits() {
	/*** Displays an 3D occupancy plot for each SM Event. (stop mode event)

	Can choose which SM event to start at. (find "CHOOSE THIS" in this script)
	Input file must be a Hits file (_interpreted_Hits.root file).
	***/
	gROOT->Reset();

	// Setting up file, treereader, histogram
	TFile *f = new TFile("/home/pixel/pybar/tags/2.0.2_new/pyBAR-master/pybar/module_202_new/101_module_202_new_stop_mode_ext_trigger_scan_interpreted_Hits.root");


	if (!f) { // if we cannot open the file, print an error message and return immediately
		cout << "Error: cannot open the root file!\n";
		//return;
	}

	TTreeReader *reader = new TTreeReader("Table", f);

	TTreeReaderValue<UInt_t> h5_file_num(*reader, "h5_file_num");
	TTreeReaderValue<Long64_t> event_number(*reader, "event_number");
	TTreeReaderValue<UChar_t> tot(*reader, "tot");
	TTreeReaderValue<UChar_t> relative_BCID(*reader, "relative_BCID");
	TTreeReaderValue<Long64_t> SM_event_num(*reader, "SM_event_num");
	TTreeReaderValue<Double_t> x(*reader, "x");
	TTreeReaderValue<Double_t> y(*reader, "y");
	TTreeReaderValue<Double_t> z(*reader, "z");

	// Initialize the canvas and graph
	TCanvas *c1 = new TCanvas("c1","3D Occupancy for Specified SM Event", 1000, 10, 900, 550);
	c1->SetRightMargin(0.25);
	TGraph2D *graph = new TGraph2D();

	// Variables used to loop the main loop
	bool endOfReader = false; // if reached end of the reader
	bool quit = false; // if pressed q
	int smEventNum = 1; // the current SM-event CHOOSE THIS to start at desired SM event number
	
	// Main Loop (loops for every smEventNum)
	while (!endOfReader && !quit) {
		// Variables used in this main loop
		int startEntryNum = 0;
		int endEntryNum = 0;
		string histTitle = "3D Occupancy for SM Event ";
		string inString = "";
		bool fitFailed = false; // true if the 3D fit failed
		bool lastEvent = false;

		// Declaring some important output values for the current graph and/or line fit
		int numEntries = 0;
		double sumSquares = 0;

		// Get startEntryNum and endEntryNum
		startEntryNum = getEntryNumWithSMEventNum(reader, smEventNum);
		endEntryNum = getEntryNumWithSMEventNum(reader, smEventNum + 1);

		if (startEntryNum == -2) { // can't find the smEventNum
			cout << "Error: There should not be any SM event numbers that are missing." << "\n";
		} else if (startEntryNum == -3) { 
			endOfReader = true;
			break;
		} else if (endEntryNum == -3) { // assuming no SM event nums are skipped
			endEntryNum = reader->GetEntries(false);
			lastEvent = true;
		}

		// Fill TGraph with points and set title and axes
		graph = new TGraph2D(); // create a new TGraph to refresh

		reader->SetEntry(startEntryNum);
		for (int i = 0; i < endEntryNum - startEntryNum; i++) {
			graph->SetPoint(i, (*x - 0.001), (*y + 0.001), (*z - 0.001));
			endOfReader = !(reader->Next());
		}

		histTitle.append(to_string(smEventNum));
		graph->SetTitle(histTitle.c_str());
		graph->GetXaxis()->SetTitle("x (mm)");
		graph->GetYaxis()->SetTitle("y (mm)");
		graph->GetZaxis()->SetTitle("z (mm)");

		graph->GetXaxis()->SetLimits(0, 20); // ROOT is buggy, x and y use setlimits()
		graph->GetYaxis()->SetLimits(-16.8, 0); // but z uses setrangeuser()
		graph->GetZaxis()->SetRangeUser(0, 40.96);
		c1->SetTitle(histTitle.c_str());

		// 3D Fit, display results, draw graph and line fit, only accept "good" events, get input
		if (!endOfReader || lastEvent) {
			// Display some results
			numEntries = graph->GetN();
			cout << "Current SM Event Number: " << smEventNum << "\n";
			cout << "Number of entries:       " << numEntries << "\n";

			// Starting the fit. First, get decent starting parameters for the fit - do two 2D fits (one for x vs z, one for y vs z)
			TGraph *graphZX = new TGraph();
			TGraph *graphZY = new TGraph();
			reader->SetEntry(startEntryNum);
			for (int i = 0; i < endEntryNum - startEntryNum; i++) {
				graphZX->SetPoint(i, (*z - 0.001), (*x + 0.001));
				graphZY->SetPoint(i, (*z - 0.001), (*y + 0.001));
//.........这里部分代码省略.........

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