void example_CovarianceEstimator() {

   int nevts = 200;

   TRobustEstimator rob(nevts, 2);

   auto f1 = new TF2("f2","bigaus",-10,10,-10,10);
   f1->SetParameters(1,0,2,0,3,0.7);

   std::vector<double> data(2);
   std::vector<double> vx(nevts), vy(nevts);
   for (int i =0; i < nevts; i++) {
      double x,y;
      f1->GetRandom2(x,y);
      data[0] = x; data[1] = y;
      rob.AddRow(data.data());
      vx[i] = x; vy[i] = y;
   }

   rob.Evaluate();

   rob.GetMean()->Print();

   rob.GetCovariance()->Print();

   auto s = new TScatter(nevts, vx.data(), vy.data());

   s->SetMarkerStyle(20);
   s->Draw("A");


   auto covMatrix = *rob.GetCovariance();
   //TMatrixD eigenVectors = covMatrixRoot->EigenVectors();
   TMatrixDSymEigen e(covMatrix);
   TVectorD eigenValues = e.GetEigenValues();
   auto eigenVectors = e.GetEigenVectors();
   eigenVectors.Print();
   eigenValues.Print();

   // this is half of ellipse major axis
   double majorAxis = sqrt(eigenValues(0));
   double minorAxis = sqrt(eigenValues(1));

   double alpha = atan2(eigenVectors(1,0),eigenVectors(0,0));

   std::cout << "alpha " << alpha * 180./TMath::Pi() << " major axis (sqrt eigenvector) " << majorAxis << std::endl;

   
   // ellipse with probability have axis = c * eigenvalues where  c= sqrt(chisquared_quantile( prob, 2)
   //ellipse at 50
   double c50 = sqrt(ROOT::Math::chisquared_quantile(0.5,2));
   TEllipse *ellipse1 = new TEllipse(0, 0, majorAxis * c50, minorAxis * c50, 0., 360., alpha * 180./TMath::Pi());
   //ellipse->SetFillColor(0);
   ellipse1->SetFillColorAlpha(9, 0.1);
   ellipse1->SetLineColor(kRed);

    // Draw ellipse
    ellipse1->Draw();

   double c90 = sqrt(ROOT::Math::chisquared_quantile(0.9,2));
   TEllipse *ellipse2 = new TEllipse(0, 0, c90 * majorAxis, c90 * minorAxis, 0., 360., alpha * 180./TMath::Pi());

   ellipse2->SetFillColorAlpha(3, 0.1);
   ellipse2->SetLineColor(kRed+3);
   ellipse2->Draw();

   double c99 = sqrt(ROOT::Math::chisquared_quantile(0.99,2));
   TEllipse *ellipse3 = new TEllipse(0, 0, c99 * majorAxis , c99 * minorAxis, 0., 360., alpha * 180./TMath::Pi());
   ellipse3->SetFillColorAlpha(4, 0.1);
   ellipse3->SetLineColor(kRed+4);
   ellipse3->Draw();

   
}