Tilt of ellipse in bigaus

ROOT
2

The “bigaus” is actually based on the “ROOT::Math::bigaussian_pdf” function (you can use the “rho”, “sigmax” and “sigmay” parameters to calculate the principal axes, which are defined by the eigenvectors of the covariance matrix, and the semidiameters, which are the square-roots of its eigenvalues):

{
  // ... create the function ...
  TF2 *f = new TF2("f", "[0] * ROOT::Math::bigaussian_pdf(x, y, [2], [4], [5], [1], [3])", -10., 10., -10., 10.);
  f->SetParNames("constant", "x0", "sigmax", "y0", "sigmay", "rho"); // note: -1 < rho < 1
  // f->Print();
  f->SetParameters(1000., 3., 2., 1., 3., -0.75); // "test" parameters' values
  // ... draw the function ...
  TCanvas *c = new TCanvas("c", "c");
  f->SetNpx(100); f->SetNpy(100);
  f->Draw("cont1z");
  gPad->SetGrid(1, 1);
  c->SetRealAspectRatio(1); // (1) or (2)
  // ... create the covariance matrix ...
  Double_t sigmax = f->GetParameter(2);
  Double_t sigmay = f->GetParameter(4);
  Double_t rho = f->GetParameter(5);
  TMatrixDSym covMatrix(2);
  covMatrix[0][0] = sigmax * sigmax; covMatrix[1][1] = sigmay * sigmay;
  covMatrix[0][1] = covMatrix[1][0] = rho * sigmax * sigmay;
  std::cout << "... covMatrix ..." << std::endl; covMatrix.Print();
  std::cout << "determinant = " << covMatrix.Determinant() << std::endl << std::endl;
  // ... calculate the principal axes ...
  TVectorD eigenValues;
  TMatrixD eigenVectors = covMatrix.EigenVectors(eigenValues);
  std::cout << "... eigenValues ..." << std::endl; eigenValues.Print();
  std::cout << "... eigenVectors ..." << std::endl; eigenVectors.Print();
  // ... draw the principal axes ...
  Double_t x0 = f->GetParameter(1);
  Double_t y0 = f->GetParameter(3);
  TArrow major(x0, y0,
               x0 + eigenVectors[0][0] * TMath::Sqrt(eigenValues[0]),
               y0 + eigenVectors[1][0] * TMath::Sqrt(eigenValues[0]),
               0.015, "|->");
  major.Draw();
  TArrow minor(x0, y0,
               x0 + eigenVectors[0][1] * TMath::Sqrt(eigenValues[1]),
               y0 + eigenVectors[1][1] * TMath::Sqrt(eigenValues[1]),
               0.015, "|->");
  minor.Draw();
}

For detailed description see:
Wikipedia → Multivariate normal distribution
Wolfram Mathworld → Bivariate Normal Distribution

1 Like