#include <iostream>
#include "TRandom3.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TTree.h"
#include "TF1.h"

int isotropic_4d(){
  TRandom3 *ran = new TRandom3();
  const double pi = 4*atan(1.0);
  
  TTree *t = new TTree("t","tree for isotropic 4d");
  TF1   *f = new TF1("f","-(1./4.)*sin(2*x)+(1./2.)*x",0.0,pi);
  
  Double_t theta, phi, psi, x1, x2, x3, x4;
  t->Branch("theta",&theta,"theta/D");
  t->Branch("phi",&phi,"phi/D");
  t->Branch("psi",&psi,"psi/D");
  t->Branch("x1",&x1,"x1/D");
  t->Branch("x2",&x2,"x2/D");
  t->Branch("x3",&x3,"x3/D");
  t->Branch("x4",&x4,"x4/D");
  
  for (int i=0; i<30000; i++){
    theta = f->GetX(ran->Uniform(0.0,pi/2.));
    phi   = acos(-ran->Uniform(-1.0,1.0));
    psi   = ran->Uniform(0.0,2*pi);
    x1    = sin(theta)*sin(phi)*cos(psi);
    x2    = sin(theta)*sin(phi)*sin(psi);
    x3    = sin(theta)*cos(phi);
    x4    = cos(theta);
    t->Fill();
  }

  TCanvas *c1 = new TCanvas("c1", "c1");
  gStyle->SetOptLogz(0);
  gStyle->SetPadTopMargin(0.15);
  gStyle->SetLabelSize(0.05,"XYZ");
  gStyle->SetTitleSize(0.10,"");
  gStyle->SetTitleXSize(0.05);
  gStyle->SetTitleXSize(0.05);
  gStyle->SetTitleYSize(0.05);
  
  c1->Divide(7,3);
  c1->cd(1);
  t->Draw("x4:x3:x2:x1","x1>0&&x2>0&&x3>0&&x4>0","colz");
  c1->cd(2);
  t->Draw("x3:x2:x1","x1>0&&x2>0&&x3>0");
  c1->cd(3);
  t->Draw("x4:x2:x1","x1>0&&x2>0&&x4>0");
  c1->cd(4);
  t->Draw("x4:x3:x1","x1>0&&x3>0&&x4>0");
  c1->cd(5);
  t->Draw("x4:x3:x2","x2>0&&x3>0&&x4>0");
  c1->cd(8);
  t->Draw("theta");
  c1->cd(9);
  t->Draw("phi");
  c1->cd(10);
  t->Draw("psi");
  c1->cd(11);
  t->Draw("x1");
  c1->cd(12);
  t->Draw("x2");
  c1->cd(13);
  t->Draw("x3");
  c1->cd(14);
  t->Draw("x4");
  c1->cd(15);
  t->Draw("x2:x1");
  c1->cd(16);
  t->Draw("x3:x1");
  c1->cd(17);
  t->Draw("x4:x1");
  c1->cd(18);
  t->Draw("x3:x2");
  c1->cd(19);
  t->Draw("x4:x2");
  c1->cd(20);
  t->Draw("x4:x3");
  c1->cd(21);

  return 0;
}