#include <iostream>
#include <fstream>
#include <cmath>

#include "TRandom.h"
#include "TString.h"
#include "TDatime.h"
#include "TFile.h"
#include "TMath.h"
#include "TTree.h"

using namespace std;

void RandomTracks(Int_t iLOOP=10000, TString noiseTYPE="n")
{
    const Double_t PI = 3.141592654;
    
    TString nsTYPE;

    if (noiseTYPE.IsNull())
    {
        std::cout << "Include noise in layers? [y(default 0.2)/n] > ";
        std::cin >> noiseTYPE;
    }

    if      (noiseTYPE=="n") {nsTYPE="NoNoise";  }
    else if (noiseTYPE=="y") {nsTYPE="WithNoise";}
    else                     {std::cout << "Exiting script" << std::endl; exit (1);}
        
    //  Number of tracks to be created.
    std::cout << "Number of tracks to be created   = " << iLOOP << std::endl << std::endl;
    Int_t COUNTER = 1;
    
    Double_t THETA, ETA, PHI;

    TDatime *myTime = new TDatime();

    TString filename = TString("RandomTracks") + nsTYPE + TString(".root");
    TString outFILE  = TString("RandomTracks") + nsTYPE + TString(".dat");

    //filename = Form("%i_%i_" + filename,myTime->GetDate(), myTime->GetTime());
    //outFILE  = Form("%i_%i_" + outFILE,myTime->GetDate(), myTime->GetTime());

    ofstream fout(outFILE);
       
    TFile* myFile = new TFile(filename,"RECREATE");

    struct RESULTS
    {
        Int_t     ientry;

        Double_t  trkTheta;
        Double_t  trkEta;
        Double_t  trkPhi;

        Double_t  vtxX;
        Double_t  vtxY;
        Double_t  vtxZ;
        Double_t  vtxR;

    };
    
    RESULTS myResults;

    TTree *myTree = new TTree("myTree","Tree with results");

    myTree->Branch("myTree",&myResults.ientry,"ientry/I:trkTheta/D:trkEta/D:trkPhi/D:vtxX/D:vtxY/D:vtxZ/D:vtxR/D",256000);
    
    for (Int_t i=0; i<=iLOOP; i++)
    {            
        if (i%COUNTER==0) {std::cout << "Number of tracks already created = " << i << std::endl;}
        
        myResults.ientry = i;
        
        //  Create a random direction (eta - phi).  --- Flat distribution.
        THETA = PI * gRandom->Rndm(i);
        ETA   = (-1.0)*log(tan(THETA/2.0));
        PHI = 2.0 * PI * gRandom->Rndm(i);

        myResults.trkTheta = THETA;
        myResults.trkEta   = ETA;
        myResults.trkPhi   = PHI; 
        
        //  Create a random vertex position.        --- Gauss distribution
        Double_t vtx[3];
        vtx[0] = gRandom->Gaus(0.0,0.003);
        vtx[1] = gRandom->Gaus(0.0,0.003);
        vtx[2] = gRandom->Gaus(0.0,28.0);

        myResults.vtxX = vtx[0];
        myResults.vtxY = vtx[1];
        myResults.vtxZ = vtx[2];
        myResults.vtxR = sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]);
        
        fout << "ientry = "   << i      << " --- " 
             << "trkTheta = " << THETA  << " --- " 
             << "trkEta = "   << ETA    << " --- " 
             << "trkPhi = "   << PHI    << " --- " 
             << "vtxX = "     << vtx[0] << " --- " 
             << "vtxY = "     << vtx[1] << " --- " 
             << "vtxZ = "     << vtx[2] << " --- " 
             << "vtxR = "     << sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]) << endl; 
        
        myTree->Fill();
                        
    }   //  for (Int_t i=0; i<iLOOP; i++)

    std::cout << std::endl;
        
    myFile->Write();
    myFile->Close();
    
}   //  Int_t main ()