C++ TGraph::GetMean方法代码示例

int RampCorrellation_single(std::string file_name, TGraphErrors &graph, TGraphErrors &graph2, vector<double> &vec_errors)
{

	TTree *t = new TTree();
	const std::string file_path = "/home/jlab/github/MagCloak_Analysis/calibration/data-calib/DATA_MegaVIEW/";
//	std::string file_name = "DataFile_2016-12-08_06-59-11.csv";
	std::string file = file_path+file_name;

	t->ReadFile(file.c_str());
//	t->Print();

	TCanvas *c = new TCanvas();
	t->Draw("TMath::Abs(B1/B2):time");
	TGraph *gh = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double ratio_mean = gh->GetMean(2);
	double ratio_std = gh->GetRMS(2);

	t->Draw("Bnom:time");
	TGraph *gh1 = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double nom_mean = gh1->GetMean(2);
	double nom_std = 0;

//	cout << "At " << nom_mean << " mT, B1/B2 is: " << ratio_mean << " +/- " << ratio_std << endl;
	
	c->Close();

	int n = graph.GetN();
	graph.SetPoint(n,nom_mean,ratio_mean);
	graph.SetPointError(n,nom_std,ratio_std);

//-------------------------------------------------------------------------------

	TCanvas *c2 = new TCanvas();
	t->Draw("B3:time");
	TGraph *gh2 = (TGraph*)c2->GetListOfPrimitives()->FindObject("Graph");
	double B3_mean = -1*(gh2->GetMean(2));
	double B3_std = gh2->GetRMS(2);
	
	int n2 = graph2.GetN();
	graph2.SetPoint(n2, nom_mean, B3_mean);
	graph2.SetPointError(n2, nom_std, B3_std);
	
	c2->Close();

	if(nom_mean < 500) vec_errors.push_back(B3_std);

	return 0;
}

/// used in display_beamprofile() : not working !
void getEllipseParameters(const float * x_data, const float * y_data, const unsigned int N, float& x_width, float& y_width, float& angle) {
	// In order to fit a good ellipse on the scattered plot :
	// 1) The TH2 is copied into a TGraph, to fit it with y(x) = ax => to retrieve the angle
	// 2) Rotation of the Graph to get the RMS in X and Y
	// 3) Creation of the final ellipse, with the good widths and angle

	TCanvas * ca0 = new TCanvas;
	ca0->Divide(2,1);
	ca0->cd(1);

	TGraph * draft = new TGraph(N,x_data,y_data);
	draft->Draw("AP");
	draft->Fit("pol1","Q");
	TF1 * pol1 = draft->GetFunction("pol1");
	pol1->Draw("same");

	// gets the angle [rad]
	angle = asin(1.) - atan(pol1->GetParameter(1));
	
	double x_datarot[N], y_datarot[N];
	for (unsigned int i=0; i<N; i++) {
		x_datarot[i]=  x_data[i]*cos(angle) - y_data[i]*sin(angle);
		y_datarot[i]=  x_data[i]*sin(angle) + y_data[i]*cos(angle);
	}

	ca0->cd(2);
	TGraph * draft2 = new TGraph(N,x_datarot,y_datarot);
	draft2->Draw("AP");
	x_width =  draft2->GetRMS(1); 
	y_width =  draft2->GetRMS(2);
	angle = 180-90*angle/asin(1.);
//	draft->Draw("AP");
	ca0->cd(1);
	TEllipse * ell = new TEllipse(draft->GetMean(1),draft2->GetMean(2),x_width*3,y_width*3);
	ell->SetTheta(angle);
	ell->Draw("same");

	//cout << "x = " << x_width << "\t y = " << y_width << "\t angle = " << angle << endl;

//	delete draft2;
//	delete draft;
//	delete ca0;
	return;
}

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

  }

  sumMean = 0;
  sumMean2 = 0;
  sumGrads = 0;
  for(int ants = 0; ants < 32; ants++){
    count8 = 0;
    for(int events = 1; events < phiAngleArray2[ants].size(); events++){

      if( deltaTArrayLoop[count8]<1){      
	deltaTArrayLoop[count8] = deltaTArray2[ants][events];
	phiAngleArrayLoop[count8] = phiAngleArray2[ants][events];
	count8++;
      }
      
    }
    
    if(count8==0){
      tempAntGraph  = new TGraph(1, dummyArray[ants], dummyArray[ants]);
    }else{
      
      tempAntGraph  = new TGraph(count8-1,  phiAngleArrayLoop, deltaTArrayLoop);
      
      if(ants == 7){  
	
	canSurf->cd(1);
	
	tempAntGraph->SetMinimum(-0.5);
	tempAntGraph->SetMaximum(0.5);
	tempAntGraph->Draw("ap");    

	tempAntGraph->SetMarkerStyle(1);
	tempAntGraph->GetXaxis()->SetLimits(0,360);
    
	sumMean = sumMean + tempAntGraph->GetMean(2)*tempAntGraph->GetMean(2);

	if(ants == 8 || ants == 16 || ants == 12 || ants == 24){
	  tempAntGraph->SetMarkerColor(8);
	}

	if(ants == 3 || ants == 7 || ants == 23 || ants == 31){
	  tempAntGraph->SetMarkerColor(1);
	}
    
	if(ants == 9 || ants == 13 || ants == 18 || ants == 26){
	  tempAntGraph->SetMarkerColor(2);
	}
    
    
	if(ants == 10 || ants == 14 || ants == 20 || ants == 28){
	  tempAntGraph->SetMarkerColor(3);
	}
    
	if(ants == 11 || ants == 15 || ants == 22 || ants == 30){
	  tempAntGraph->SetMarkerColor(4);
	}
    
	if(ants == 0 || ants == 4 || ants == 17 || ants == 25){
	  tempAntGraph->SetMarkerColor(5);
	}


	if(ants == 2 || ants == 6 || ants == 21 || ants == 29){
	  tempAntGraph->SetMarkerColor(6);
	}

	if(ants == 1 || ants == 5 || ants == 19 || ants == 27){
	  tempAntGraph->SetMarkerColor(7);
	}
    
	tempAntGraph->GetXaxis()->SetTitle("phi (degrees)");
	tempAntGraph->GetYaxis()->SetTitle("actual - expected time");
	myMG2->Add(tempAntGraph);
	myMG2->Draw("p");
	
 	vector<double> myFit = leastSquares(phiAngleArrayLoop, deltaTArrayLoop, count8-1);
	double slope = myFit[0];
	double intercept = myFit[1];
	
	double tempX[2] = {slope*(meanPhi[ants]-20)+intercept,slope*(meanPhi[ants]+20)+intercept};
	double tempY[2] = {(meanPhi[ants]-20),(meanPhi[ants]+20)};
	sumGrads = sumGrads + myFit[0]*myFit[0]*10000;
      
      }
    } 
    

  }
  
  cout << " " << endl;
  cout << sumMean << "  " << sumGrads <<endl;
  canSurf->Update();
  cout << "  " << endl;  

  theReturn = sumMean+sumGrads;

  return theReturn;
   
}

void fitcorr(){

   setTDRStyle();

   TFile *f = new TFile("rootfiles/fitresults_2D_data.root");
   TTree *t = (TTree*)f->Get("FitResults");

   double mt=0, jsf=0, mcmass=0;
   t->SetBranchAddress("mt", &mt);
   t->SetBranchAddress("jesfactor", &jsf);
   t->SetBranchAddress("mcmass", &mcmass);

   TH2D *h2D = new TH2D("h2D",";M_{t} [GeV];JSF", 20, 170, 173, 20, 0.99, 1.04 );
   TGraph *g = new TGraph();
   for( int i=0; i < t->GetEntries(); i++ ){
      t->GetEntry(i);
      //if( mcmass != 172.5 ) continue;

      h2D->Fill(mt, jsf);
      g->SetPoint(g->GetN(), mt, jsf);

   }

   // uncertainty contour
   double mean_mt = g->GetMean(1);
   double mean_jsf = g->GetMean(2);
   double sigma_mt = g->GetRMS(1);
   double sigma_jsf = g->GetRMS(2);
   double rho = g->GetCorrelationFactor();
   std::cout << "Mean Mt = " << mean_mt << " +- " << sigma_mt << std::endl;
   std::cout << "Mean JSF = " << mean_jsf << " +- " << sigma_jsf << std::endl;
   std::cout << "Correlation = " << rho << std::endl;

   TMatrixDSym m(2);
   m(0,0) = sigma_mt*sigma_mt;
   m(1,1) = sigma_jsf*sigma_jsf;
   m(0,1) = rho*sigma_mt*sigma_jsf;
   m(1,0) = rho*sigma_mt*sigma_jsf;
   TMatrixDSymEigen eigen(m);
   TVectorD eigenval = eigen.GetEigenValues();
   std::cout << eigenval[0] << " " << eigenval[1] << endl;

   // convert to ellipse tilt and radius
   //double phi = 0.5*TMath::ATan( (3/0.03)*(2*rho*sigma_mt*sigma_jsf)/(sigma_mt*sigma_mt-sigma_jsf*sigma_jsf) );
   double phi = (180.0/TMath::Pi())*0.5*TMath::ATan( (2*rho*sigma_mt*sigma_jsf)/(sigma_mt*sigma_mt-sigma_jsf*sigma_jsf) );
   std::cout << "phi = " << phi << std::endl;

   TEllipse *el = new TEllipse(mean_mt, mean_jsf, sqrt(eigenval[0]), sqrt(eigenval[1]), 0, 360, phi);
   TEllipse *el2 = new TEllipse(mean_mt, mean_jsf, 2*sqrt(eigenval[0]), 2*sqrt(eigenval[1]), 0, 360, phi);
   
   /*
   TF2* fell = new TF2("fell","pow((x-[1])*cos([0])+(y-[2])*sin([0]),2)/[3] + pow((x-[1])*sin([0])-(y-[2])*cos([0]),2)/[4]",171,174,0.985,1.015);
   fell->SetParameter(0,phi);
   fell->SetParameter(1,mean_mt);
   fell->SetParameter(2,mean_jsf);
   fell->SetParameter(3,eigenval[0]);
   fell->SetParameter(4,eigenval[1]);

   double cont[1] = {1};
   fell->SetContour(1, cont);
   */
   
   const Int_t Number = 2;
   Double_t Red[Number]    = { 0.90, 0.10};
   Double_t Green[Number]  = { 0.90, 0.10};
   Double_t Blue[Number]   = { 0.90, 0.10};
   Double_t Stops [Number] = {0, 1};
   Double_t Length [Number] = {0, 1};
   Int_t nb=50;
   TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);

   gStyle->SetPadRightMargin(0.16);
   gStyle->SetPadTopMargin(0.08);

   TCanvas *c = new TCanvas("c","c",800,600);
   h2D->Draw("colzC");
   el->SetFillStyle(0);
   el->SetLineColor(2);
   el->SetLineWidth(3);
   el->Draw();
   el2->SetFillStyle(0);
   el2->SetLineColor(2);
   el2->SetLineWidth(3);
   el2->Draw();

   int iPeriod = 2;
   int iPos = 33;
   lumi_sqrtS = "8 TeV";
   writeExtraText = false;
   CMS_lumi( c, iPeriod, iPos );

   c->Update();
   c->RedrawAxis();
   c->GetFrame()->Draw();

   TLatex latex;
   latex.SetNDC();
   latex.SetTextSize(0.06);
   latex.SetTextFont(42);
   latex.DrawLatex(0.2, 0.84, "2D fit");
//.........这里部分代码省略.........