//
// global constants
//

const Long64_t    kNFills = 1000000;

const Int_t    kNTimeBins       = 500;
const Int_t    kNZenithAngBins  = 18;
const Int_t    kNAzimuthAngBins = 36;
const Int_t    kNCosThetaBins = 200; 

//const Int_t    kNTimeBins       = 50;
//const Int_t    kNZenithAngBins  = 18;
//const Int_t    kNAzimuthAngBins = 36;
//const Int_t    kNCosThetaBins = 20; 

//const Int_t    kNTimeBins       = 500;
//const Int_t    kNZenithAngBins  = 4;
//const Int_t    kNAzimuthAngBins = 4;
//const Int_t    kNCosThetaBins = 4; 

const Double_t kCosThetaBinSize = (2.0 / kNCosThetaBins); 
const Double_t kCosThetaBinEdge = (1.0 - kCosThetaBinSize);
const Double_t kPi = TMath::Pi();
const Double_t kRadToDeg = (180. / kPi);

// Test program for TH4D class
Int_t TH4DTest()
{
   gSystem->Load("libPhysics.so");
   gSystem->Load("libPhysics.so");
   gSystem->Load("../../lib/libutils.so");
   gStyle->SetOptStat("ne");
   //gDebug = 100;

   TFile ofile("test.root", "recreate");

   TH4D *h4 = new TH4D("test","test",
                        kNCosThetaBins, -1., 1.,
                        kNZenithAngBins, 0., 180.,
                        kNAzimuthAngBins, 0., 360.,
                        kNTimeBins, 0., 1500.);

   // fill gaussian at specific bin

   Double_t costheta = 0.90;
   Double_t meanzenith   = 90;
   Double_t sigmazenith  = 20;
   Double_t meanazimuth  = 180;
   Double_t sigmaazimuth = 20;
   Double_t sigmatime    = 5.;  
   Double_t meantime, time, zenith, azimuth;

   for (Long64_t i=0; i<kNFills; i++) {
      zenith   = gRandom->Gaus(meanzenith, sigmazenith);
      azimuth  = gRandom->Gaus(meanazimuth, sigmaazimuth);
      meantime = (azimuth / 3.0);
      time = gRandom->Gaus(meantime, sigmatime);
      h4->Fill(costheta, zenith, azimuth, time);
      if (i % 100000 == 0) cerr << i << " event finished! " << endl;
   }      

   //Int_t bin = h4->FindBin(costheta, zenith, azimuth);
   //TH1DEX* exdim = h4->GetExtraDimensionPtr(bin);
   //exdim->Draw();

   ofile.Write(); 
   ofile.Close();

}