#include <iostream>
//#include <fstream>
#include <cstdio>
#include "TFile.h"
#include "TH1.h"
#include "TRandom.h"
#include "TRandom1.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TMath.h"
#include "TClonesArray.h"
#include "TObjArray.h


 
 // Define some simple classes
using namespace std;

class Eventn : public TObject {
public:

	void EventDump();
	
	Double_t COSTHET;
        Double_t SIN2THET;
        Double_t SINTHET;
                	
	Double_t COSPSI;
	Double_t SIN2PSI;	
	Double_t PSI;
	
	Double_t COSCHI;	
        Double_t SIN2CHI;
        
	Double_t PHI;
	Double_t PHIK;
	
        Double_t SIN2THLP;
        Double_t COSTHLP;	
                
        Double_t SIN2THMU;
        Double_t COSTHMU;	
        
	Double_t dist[3];
	Double_t lnlike[3];

    	Int_t    evtn;

	Eventn() {COSTHET=0.0; SIN2THET=0.0; SINTHET=0.0; COSPSI=0.0; SIN2PSI=0.0; PSI=0.0;
	          COSCHI=0.0; SIN2CHI=0.0; PHI=0.0; PHIK=0.0; SIN2THLP=0.0; COSTHLP=0.0;
	          SIN2THMU=0.0; COSTHMU=0.0; dist[0]=0.0; dist[1]=0.0; dist[2]=0.0; 
	          lnlike[0]=0.0; lnlike[1]=0.0; lnlike[2]=0.0; evtn=0;}

	ClassDef(Eventn,1)
};

void Eventn::EventDump() // this needs to be changed for new variables 04.06.2008
{
	std::cout << "Dumping Eventn: " << std::endl
		  << "... evtn:    " << evtn    << std::endl
		  << "... COSTHET: " << COSTHET << std::endl
		  << "... COSPSI:  " << COSPSI  << std::endl
		  << "... COSCHI:  " << COSCHI  << std::endl
		  << "... PHI:     " << PHI     << std::endl
		  << "... PHIK:    " << PHIK    << std::endl
		  << "... dist:    " << dist[0] << " "
		                     << dist[1] << " "
		                     << dist[2] << std::endl
		  << "... lnlike:  " << lnlike[0] << " "
		                     << lnlike[1] << " "
		                     << lnlike[2] << std::endl
		  << "... SIN2THLP:" << SIN2THLP  << std::endl
		  << "... SIN2THET:" << SIN2THET  << std::endl
		  << "... SIN2PSI: " << SIN2PSI   << std::endl
		  << "... SIN2CHI: " << SIN2CHI   << std::endl
		  << "... SIN2THMU:" << SIN2THMU  << std::endl
		  << std::endl;
}

class Experiment : public TObject {
public:
	void ExperimentDump();
	
	Double_t SLNLIKE[3];
	Int_t    expn;

	TClonesArray* evtarray;
	
	Experiment() {SLNLIKE[0] = 0.0; SLNLIKE[1] = 0.0; SLNLIKE[2] = 0.0;
	              expn=0; evtarray=0;}

	ClassDef(Experiment,1)
};

void Experiment::ExperimentDump()
{
	std::cout << "Dumping experiment " << expn << std::endl
		  << "SLNLIKE:"
		  << " " << SLNLIKE[0]
		  << " " << SLNLIKE[1]
		  << " " << SLNLIKE[2]
		  << " expn: " << expn << std::endl;
	std::cout << " evt: size = " <<  evtarray->GetEntries() << std::endl;
	for (Int_t i=0; i< evtarray->GetEntries(); i++) {
		std::cout << "Dumping event " << i << std::endl;
		Eventn* eventn = (Eventn*) (*evtarray)[i];
		eventn->EventDump();
	}
}

//________________________Make Tree______________________________________________________
int  tree_buggy()//int argc, char **argv)
 {

//____________________Initialising constants_____________________________________ 
  Double_t  CN_EPSIL;
  Int_t  df_lux, df_rseed;
  Int_t  df_nexp, df_nevt;
  Int_t  jdist_min, jdist_max, df_jdist;
 
      CN_EPSIL = 1.0E-8;
//.. RANLUX luxury level and seed
      df_lux = 4;
      df_rseed = 0;
//.. default number of events-per-experiment, experiments, reporting quantum
      df_nevt = 5;
      df_nexp = 10;
//.. minimum, maximum, & default decay-model numbers
      jdist_min = 0;
      jdist_max = 2; 
      df_jdist = 2;
      
      Int_t jdist, evt_max, exp_max, rand_seed;   
      
      cout<<"INITIALISATION";
      
      int  yaynay;
      
      jdist       = df_jdist;
      evt_max     = df_nevt;
      exp_max     = df_nexp;
      rand_seed   = df_rseed;
      
      cout << "Do you wish to change all values y/n (1/0)";
      cin  >> yaynay;
      
      if(yaynay == 1){
      
         cout << "Dump of parameters read from datacards:\n";    
         
         
         cout << "\n Number of Experiemnts: ";
         cin  >> exp_max;
         cout << "\n Number of events per experiment: "; 
         cin  >> evt_max; 
         cout << "\n Decay Model (0-2): ";
         cin  >> jdist;
         cout << "\nRANLUX seed: ";
         cin  >> rand_seed;

      }

      cout << "Initialisation complete \n";
      
      
      if(jdist > 2 || jdist < 0){
      std::cout << "distribution number out of allowed range (" << jdist_min << "," << jdist_max << "," << jdist << ")"<< std::endl;}

//__________________Initialisation complete________________________________________________________
 
 // Create a new ROOT binary machine independent file.
 // Note that this file may contain any kind of ROOT objects, histograms,trees
 // pictures, graphics objects, detector geometries, tracks, events, etc..
 // This file is now becoming the current directory.

   TFile hfile("trial.root","RECREATE","Demo ROOT file with histograms & trees");

 // Create some histograms and a profile histogram

   TH1F *hdist0  = new TH1F("hdist0","distribution #0",100,-4,4);
   TH1F *hdist1  = new TH1F("hdist1","distribution #1",100,-4,4);
   TH1F *hdist2  = new TH1F("hdist2","distribution #2",100,-4,4);

// Create a ROOT Tree

   Experiment* experiment = new Experiment();
   TTree *tree = new TTree("T","An example of ROOT tree");
   experiment->evtarray = new TClonesArray("Eventn",1000);
   tree->Branch("experiment",&experiment);
   tree->Branch("hdist0","TH1F",&hdist0,128000,0);

//Set the seed --- this uses the TUUID function...

   TRandom1 SetSeeds(df_rseed, df_lux);
   
   for(Int_t exp_i = 0; exp_i < exp_max; exp_i++) {              

	 //experiment->evtarray = new TClonesArray("Eventn",1000);
	 Int_t expn_t = exp_i;
   
	 Int_t    event_i, evtn_t;
	 Double_t COSTHLP_t, COSTHMU_t, COSCHI_t, COSTHET_t, COSPSI_t,
		  PHI_t,PHIK_t, PSI_t,
		  SINTHET_t, SIN2THMU_t, SIN2THLP_t, SIN2CHI_t,
		  SIN2THET_t, SIN2PSI_t,    
		  dist_t[3], lnlike_t[3], slnlike_t[3];
   	
	 COSTHET_t=0.0; 
	 COSPSI_t=0.0; 
	 COSCHI_t=0.0; 
	 PHI_t=0.0; 
	 PHIK_t=0.0;
	 evtn_t=0; 
	 SIN2THLP_t=0.0; 
	 SIN2THET_t=0.0; 
	 SIN2PSI_t=0.0; 
	 SIN2CHI_t=0.0; 
	 SIN2THMU_t=0.0;

	 dist_t[0] = 0.0;
	 dist_t[1] = 0.0;
	 dist_t[2] = 0.0;   
   
	 lnlike_t[0] = 0.0;
	 lnlike_t[1] = 0.0;
	 lnlike_t[2] = 0.0;              
                                    
	 slnlike_t[0] = 0.0;
	 slnlike_t[1] = 0.0;
	 slnlike_t[2] = 0.0;
//--------------------Here we start a loop on 1000 events
	 for(event_i = 0; event_i < evt_max ; event_i++) {
   
		 Eventn* evt = new Eventn();
		 evtn_t = event_i;   
		 Int_t c = 0;

		 Double_t RAN[6];     
		 do
			 {

				 TRandom1 r(0);
				 r.RndmArray(6,RAN);
     
				 /*for(Int_t k = 0; k < 6 ; k++){
					 cout << RAN[k] << "\n";
				 }*/

				 COSTHET_t  = 2.0*RAN[0] - 1.0;
				 SIN2THET_t    = 1 - pow(COSTHET_t,2);
                    
				 if(SIN2THET_t > CN_EPSIL) {
					 SINTHET_t = sqrt(SIN2THET_t);
				 } 
				 else {
					 SINTHET_t = 0.0;
					 SIN2THET_t = 0.0;
				 }
  
				 PHI_t      = TMath::TwoPi() * RAN[1];
  
				 COSPSI_t   = 2.0*RAN[2] - 1.0;
				 SIN2PSI_t     = 1 - pow(COSPSI_t,2);
				 PSI_t         = TMath::ACos(COSPSI_t);
                                 // ..     n.b. psi, as a polar angle,
				 //        lies only in the 1st&2nd quadrants
          
				 COSCHI_t   = 2.0*RAN[3] - 1.0;
				 SIN2CHI_t  = 1 - pow(COSCHI_t,2);
  
				 PHIK_t     = TMath::TwoPi() * RAN[4];
  
				 COSTHLP_t = SINTHET_t *
					 TMath::Cos(PHI_t-PSI_t);
				 SIN2THLP_t = 1 - pow(COSTHLP_t,2);
  
				 COSTHMU_t  = SINTHET_t *
					 TMath::Sin(PHIK_t) *
					 TMath::Sin(PHI_t-PSI_t) -
					 COSTHET_t*TMath::Cos(PHIK_t);
				 SIN2THMU_t = 1 - pow(COSTHMU_t,2);

				 dist_t[0]   = SIN2THLP_t;      //J^PC = O^++
				 dist_t[1]   = SIN2THET_t * SIN2PSI_t; //J^PC = O^-+
				 dist_t[2]   = SIN2CHI_t * SIN2THMU_t; //J^PC = 1^++
     

     
				 c = c+1;
     
				 //std::cout<<"counter "<<c <<" "
				//	  <<dist_t[0] <<"\n";
     
			 }
		 while(RAN[5] > dist_t[jdist]);
          
		 lnlike_t[0] = log(TMath::Max(dist_t[0],CN_EPSIL));
		 lnlike_t[1] = log(TMath::Max(dist_t[1],CN_EPSIL));
		 lnlike_t[2] = log(TMath::Max(dist_t[2],CN_EPSIL));
     
		 c = 0;

	
	

		 evt->evtn = evtn_t;
     
		 evt->COSCHI    = COSCHI_t; 
		 evt->SIN2CHI    = SIN2CHI_t; 
		       
		 evt->COSTHET   = COSTHET_t;
		 evt->SIN2THET  = SIN2THET_t;
		 evt->SINTHET   = SINTHET_t;		 

		 evt->COSPSI    = COSPSI_t;
		 evt->SIN2PSI   = SIN2PSI_t;
		 evt->PSI       = PSI_t;
		 		 		 
		 evt->PHI       = PHI_t;
		 evt->PHIK      = PHIK_t;
     
		 evt->SIN2THMU  = SIN2THMU_t;
		 evt->COSTHMU   = COSTHMU_t;
		 
		 evt->SIN2THLP  = SIN2THLP_t;
                 evt->COSTHLP   = COSTHLP_t;
     
     
		 evt->dist[0]   = dist_t[0];
		 evt->dist[1]   = dist_t[1];
		 evt->dist[2]   = dist_t[2];
           
		 evt->lnlike[0] = lnlike_t[0];
		 evt->lnlike[1] = lnlike_t[1];
		 evt->lnlike[2] = lnlike_t[2];
     
		 slnlike_t[0]     = slnlike_t[0] + lnlike_t[0];
		 slnlike_t[1]     = slnlike_t[1] + lnlike_t[1];     
		 slnlike_t[2]     = slnlike_t[2] + lnlike_t[2];  
        
		 hdist0->Fill(evt->dist[0]);
		 hdist1->Fill(evt->dist[1]);
		 hdist2->Fill(evt->dist[2]);
     
		 // Fill the tree. For each event, save the
		 // 2 structures and 3 objects   
		 evt->EventDump();
		 //(*evtarray)[event_i] = evt;
		 
		 new ((*(experiment->evtarray))[event_i]) Eventn(*evt); 

		// Next two lines are OK
		//Eventn* evt1 = (Eventn*) (*(experiment->evtarray))[event_i];
		//evt1->EventDump();
	 }
	 //--------------End of the loop over events
	 experiment->expn = expn_t;
   
	 experiment->SLNLIKE[0] = slnlike_t[0];
	 experiment->SLNLIKE[1] = slnlike_t[1];
	 experiment->SLNLIKE[2] = slnlike_t[2];
 
	 experiment->ExperimentDump();
	 //tree->Print();
	 
	 tree->Fill(); 

         //   experiment->Clear();
   
   }
//--------------- End of loop over experiments    

   //TRandom WriteRandom("seed");
   TRandom1 GetLuxury();


   tree->Print();

 // Save all objects in this file
   hfile.Write();

 // Close the file. Note that this is automatically done when you leave
 // the application.
   hfile.Close();

//   TRandom::WriteRandom(seeds);

   return 0;
 }