// generate two sets of random points according to two distributions
// calculate sigma for each one and correlation pho (zero, dist. uncorrelated)
// rotate the points, calculate new correlation phocor
// calculate the angle from sigmas and phocor
//
//         jtmn, Ago, 2003


void correl()
{

 gROOT->Reset();
 gROOT->SetStyle("Plain"); 

 const Int_t nPoints = 20000;
 const Double_t xwid = 7;
 const Double_t ywid = 7;

 TRandom *r3 = new TRandom3();

 // generate x dist
 Float_t  x[nPoints];
 for (Int_t i=0;i<nPoints;i++) {
   x[i] = r3->Gaus(0,2);
  }

 // histogram x dist
 TH1F *hx = new TH1F("hx", " ", 100, -xwid,xwid); 
 for (Int_t i=0;i<nPoints;i++) {
   hx->Fill(x[i]);
  }
  Float_t sx= hx->GetRMS();
  cout << "sx = " << sx << endl; 

  // generate y dist
 Float_t y[nPoints];
 for (Int_t i=0;i<nPoints;i++) {
     y[i] = r3->Gaus(0,1);
  }
 // histogram y dist
 TH1F *hy = new TH1F("hy", " ", 100, -ywid,ywid);
  for (Int_t i=0;i<nPoints;i++) {
   hy->Fill(y[i]);
  }
 Double_t sy= hy->GetRMS();
 cout << "sy = " << sy << endl; 

 // display bidim dist (uncorrelated)
 TCanvas *c1 = new TCanvas("c1","uncorrelated data", 600,600);
 plane  = new TH2F("plane","data points", 50,-xwid,xwid,50,-ywid,ywid);
 for (Int_t i=0;i<nPoints;i++){
   plane->Fill(x[i],y[i]);
  }
 plane->SetStats(kFALSE);
 plane->GetXaxis()->SetTitle("x_{uncorr}");
 plane->GetYaxis()->SetTitle("y_{uncorr}");
 plane.Draw();
 Double_t r = plane->GetCorrelationFactor();
 cout << " r = " << r << endl;
 c1->Update();

 // theta rotation angle
 Double_t pi = TMath::Pi();
 Double_t theta = 30.0*pi/180.0;  // rad
 cout << "theta = " << theta << endl; 

 // rotate the points
 Float_t xcor[nPoints];
 Float_t ycor[nPoints];
 for(Int_t i=1;i<nPoints; i++){
   xcor[i] = x[i]*TMath::Cos(theta) - y[i]*TMath::Sin(theta);
   ycor[i] = x[i]*TMath::Sin(theta) + y[i]*TMath::Cos(theta);
 }

 // histogram xcor to get RMS
 TH1F *hxcor = new TH1F("hxcor", " ", 100, -xwid,xwid);
 for (Int_t i=0;i<nPoints;i++) {
   hxcor->Fill(xcor[i]);
  }
 Double_t sxcor= hxcor->GetRMS();
 cout << "sxcor = " << sxcor << endl; 
 
 // histogram ycor to get RMS
 TH1F *hycor = new TH1F("hycor", " ", 100, -ywid,ywid);
 for (Int_t i=0;i<nPoints;i++) {
   hycor->Fill(ycor[i]);
  }
 Double_t sycor= hycor->GetRMS();
 cout << "sycor = " << sycor << endl; 


 // display the rotated points
 TCanvas *c2 = new TCanvas("c2","correlated data", 600,600);
 planecor  = new TH2F("planecor","data points", 50,-xwid,xwid,50,-ywid,ywid);
 for (Int_t i=0;i<nPoints;i++){
   planecor->Fill(xcor[i],ycor[i]);
  } 
 planecor->SetStats(kFALSE);
 planecor->GetXaxis()->SetTitle("x_{corr}");
 planecor->GetYaxis()->SetTitle("y_{corr}");
 planecor.Draw();
 Double_t rcor = planecor->GetCorrelationFactor();
 cout << " rcor = " << rcor << endl;
 
 // check the formulas

 Double_t thetacalc = 0.0;
 if (TMath::Abs(sxcor-sycor)>=0.001){
     thetacalc = 0.5*TMath::ATan((2.0*rcor*sxcor*sycor)/(sxcor*sxcor - sycor*sycor));
     }
 cout << "thetacalc = " << thetacalc << endl;

 c2->Update();
}