//
// global constants
//

//const Int_t kNFills = 10000000;
const Int_t kNFills = 100000;

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

   TFile *infile = TFile::Open("test.root", "read");
   TH4D *h4 = gROOT->FindObject("test"); 

   Int_t nexbins = h4->GetExtraAxis()->GetNbins();
   Double_t ulow = h4->GetExtraAxis()->GetXmin();
   Double_t uup  = h4->GetExtraAxis()->GetXmax();

   cerr << "nbins, ulow, uup = " << nexbins << " " << ulow << " " << uup << endl;

   TH4D *newhist = new TH4D("filltest", "filltest",
          h4->GetXaxis()->GetNbins(), h4->GetXaxis()->GetXmin(), h4->GetXaxis()->GetXmax(),
          h4->GetYaxis()->GetNbins(), h4->GetYaxis()->GetXmin(), h4->GetYaxis()->GetXmax(),
          h4->GetZaxis()->GetNbins(), h4->GetZaxis()->GetXmin(), h4->GetZaxis()->GetXmax(),
          h4->GetExtraAxis()->GetNbins(), h4->GetExtraAxis()->GetXmin(), h4->GetExtraAxis()->GetXmax());

   newhist->FillRandom(h4, kNFills);

   TCanvas *c1 = new TCanvas("TH4D","TH4D", 10, 10, 400, 400);
   c1->cd();
   newhist->Draw();

     
   TCanvas *c2 = new TCanvas("extra_dim","extra_dim", 10, 10, 800, 800);
   c2->Divide(2, 2);

   TH1DEX *exdim1 = newhist->GetExtraDimensionPtr(newhist->FindBin(0.90, 85, 160));
   if (exdim1) { 
      c2->cd(1); exdim1->Draw(); 
   } else { 
      cerr << "exdim1 is 0!" << endl; 
   }
   TH1DEX *exdim2 = newhist->GetExtraDimensionPtr(newhist->FindBin(0.90, 90, 180));
   if (exdim2) { 
      c2->cd(2); exdim2->Draw(); 
   } else { 
      cerr << "exdim2 is 0!" << endl; 
   }
   TH1DEX *exdim3 = newhist->GetExtraDimensionPtr(newhist->FindBin(0.90, 95, 200));
   if (exdim3) { 
      c2->cd(3); exdim3->Draw(); 
   } else { 
      cerr << "exdim3 is 0!" << endl; 
   }
   TH1DEX *exdim4 = newhist->GetExtraDimensionPtr(newhist->FindBin(0.90, 100, 220));
   if (exdim4) { 
      c2->cd(4); exdim4->Draw(); 
   } else { 
      cerr << "exdim4 is 0!" << endl; 
   }

}