#include "TRandom3.h"
#include "TColor.h"
#include "TFile.h"
#include "TH3.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TTree.h"
#include "TMath.h"
#include "TStyle.h"

void set_plot_style(){
    const Int_t NRGBs = 7;
    const Int_t NCont = 180;
    Double_t stops[NRGBs] = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00   };
    stops[1]=1./6.;
    stops[2]=2./6.;
    stops[3]=3./6.;
    stops[4]=4./6.;
    stops[5]=5./6.;
    Double_t red[NRGBs]   = { 1.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00   };
    Double_t green[NRGBs] = { 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.00   };
    Double_t blue[NRGBs]  = { 0.00, 0.00, 0.00, 1.00, 1.00, 1.00, 0.00   };
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
}

void sferiche_nt(){
  set_plot_style();
  TFile *f=new TFile("sferiche_nt.root","recreate");
  f->cd();
  TTree *t=new TTree("t","t");
  Double_t x=0; t->Branch("x",&x,"x/D");
  Double_t y=0; t->Branch("y",&y,"y/D");
  Double_t z=0; t->Branch("z",&z,"z/D");
  Double_t cost=0; t->Branch("cost",&cost,"cost/D");
  Double_t theta=0; t->Branch("theta",&theta,"theta/D");
  Double_t phi=0; t->Branch("phi",&phi,"phi/D");
  Double_t p=0; t->Branch("p",&p,"p/D");
  TRandom3 random;
//  gStyle->SetPalette(1,0);
  //t->Fill(); p=2*TMath::Pi(); t->Fill();
  for(Int_t i=0; i<20000; i++){
    //gRandom->Sphere(x,y,z,1);
    cost=random.Uniform(-1,1);
    theta=TMath::ACos(cost);
    phi=2.*TMath::Pi()*random.Uniform();
//    Double_t r=-TMath::Sqrt(3./8./TMath::Pi())*TMath::Sin(theta); Double_t angle=phi; // Y11
    Double_t r=TMath::Sqrt(3./4./TMath::Pi())*TMath::Cos(theta); Double_t angle=0; // Y10
    if(r<0){
      r=TMath::Abs(r);
      angle=angle+TMath::Pi();
    }
    Double_t pro=r*TMath::Sin(theta);
    x=pro*TMath::Cos(phi);
    y=pro*TMath::Sin(phi);
    z=r*TMath::Cos(theta);
    p=angle-2.*TMath::Pi()*TMath::Floor(angle/TMath::Pi()/2.);
    t->Fill();
  }
  t->Draw("z:y:x:p","","colz");
  TH3F *ht=(TH3F *)gROOT->FindObject("htemp");
  if(ht){
    ht->GetXaxis()->SetRangeUser(-0.5,0.5);
    ht->GetYaxis()->SetRangeUser(-0.5,0.5);
    ht->GetZaxis()->SetRangeUser(-0.5,0.5);
    ht->SetMinimum(0);
    ht->SetMaximum(2.*TMath::Pi());
    ht->Draw("colz");
  }else{
    printf("Cannot get htemp\n");
  }
}