/*
// Version 0.1a  9/2007   Gorbunov and Butterworth
// initial version based on starlight in fortran
*/

#include <fstream>
#include <iostream>
using namespace std;
#include <math.h>
#include <time.h>

#include "main.h"
//#include "inputparameters.h"
#include "starlightconstants.h"
#include "nucleus.h"
#include "beam.h"
#include "beambeamsystem.h"
#include "twophotonluminosity.h"
#include "readinluminosity.h"
#include "gammagammaleptonpair.h"
#include "gammagammasingle.h"
#include "gammaaluminosity.h"
#include "gammaavm.h"
#include "psifamily.h"
#include <TTree.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TNtuple.h>
#include <TParticle.h>
#include "Particle.h"
#include "masses.h"

//#include "Pythia.h"
#include "PythiaStarlight.h"

//using namespace Pythia8;//Doubling up on namespaces =/

Pythia *pythia;
//Event& event; //Then it will be used in gammagammasinglemeson for pythia
		//linked via extern in that file.
//Initializing masses.h values for mac OS

const Double_t Masses::fMassMu = 0.105658;
const Double_t Masses::fMassEl = 0.00051099891;
const Double_t Masses::fMassPi = 0.13957018;



int main() {
         //Storing time.
        time_t currenttime=0;
        //Structure to read the time easily
        struct tm * readabletime;
//Reading time at the beginning of the calculation
        time(&currenttime);
        readabletime=localtime(&currenttime);
        cout<<"Beginning  at: "<<asctime(readabletime)<<endl;

	// Generator; shorthand for event.
        pythia = new Pythia("/home/butter/pythia/pythia8120/xmldoc");

	//Event& event = pythia.event;
	//event=pythia.event;        
        // Key requirement: switch off ProcessLevel, and thereby also PartonLevel.
        pythia->readString("ProcessLevel:all = off");
        //When limittau0 is on, only particles with tau0<tau0max are decayed.
        pythia->readString("ParticleDecays:limitTau0 = on");
        //Default tau0 max is 10mm/c, we are setting it to 1mm/c
        pythia->readString("ParticleDecays:tau0Max = 1");
        // Provide printout of initial information.
	pythia->settings.listChanged();//Useful for earlier stages.
	//Initialize pythia
	pythia->init();
	//Reseting the event record
	//pythia->event.reset();
	PythiaOutput=false;//If we use pythia, it will be set true elsewhere.
//	Inputparameters InputInformation; 
    InputInformation;
    Beam beam1(InputInformation.getZ1(),InputInformation.getA1(),InputInformation.getbford(),InputInformation.getincoherentorcoherent(),InputInformation);
    Beam beam2(InputInformation.getZ2(),InputInformation.getA2(),InputInformation.getbford(),InputInformation.getincoherentorcoherent(),InputInformation);
    Beambeamsystem bbs1(beam1,beam2,InputInformation); //redundant to construct two...
//	Twophotonluminosity(beam1,beam2,InputInformation.ibreakup,InputInformation);
//	Gammaainteraction(beam1,beam2,Inputparameters::ibreakup);
//	cout<<"woodsaxonskin: "<<beam1.getWoodSaxonskindepth()<<endl;
//	cout<<"A,Z: "<<beam1.getAin()<<" "<<beam1.getZin()<<endl;
//	cout<<"nofe: "<<beam1.nofe(3.)<<endl;
    streamsize precision(15);
    cout.setf(ios_base::fixed,ios_base::floatfield);
    cout.precision(15);
    bool flag;

    //Vector to store all the events
    vector<StarlightConstants::event> allevents;
    allevents.reserve(InputInformation.getnumberofevents());//This way the proper memory is allocated for the vector
	
    checkforluminositytable(flag);
	
    //Determine luminosity tables for the interaction type
    //Should put flag==true outside of switch, less code and prettier.
	
    switch(InputInformation.getinteractiontest())
    {
    case StarlightConstants::PHOTONPHOTON:
	if(flag==true){
	    cout<<"Flag is True"<<endl;
	    Twophotonluminosity(beam1,beam2,InputInformation.getbreakupmode(),InputInformation);
	}
	break;
    case StarlightConstants::PHOTONPOMERONNARROW://Both narrow and wide use the same luminosity function.
    case StarlightConstants::PHOTONPOMERONWIDE:
	if(flag==true){
	    cout<<"Flag is True"<<endl;
	    //Luminosity function
	    Gammaaluminosity(InputInformation,bbs1);
	}
	break;
    default :
	cout << "Please go back and define an appropriate interaction type. Thank you."<<endl;
    }


    //Now do it based on decay type...
    int j =0; //while counter
    StarlightConstants::event thisevent;
    Eventchannel* channel = CreateEventChannel(InputInformation,bbs1);
    while(j<InputInformation.getnumberofevents())
    {
	thisevent=channel->produceevent(j);
        //Now store this to a vector defined earlier with numberofevents. We will do the writing out later.
	allevents.push_back(thisevent);
    }//while loop

/*
  switch(InputInformation.getdecaytest())
  {
  case StarlightConstants::LEPTONPAIR:
  {
  Gammagammaleptonpair leptonpair(InputInformation,bbs1);
  cout<<"Leptonpairs being produced"<<endl;
  //leptonpair.read();
  //leptonpair.twoleptoncrosssection();
  //if(InputInformation.getparticleid()==StarlightConstants::TAUON) leptonpair.tablecalc();
	
  //StarlightConstants::event thisevent;

  while(j<InputInformation.getnumberofevents())
  {
  thisevent=leptonpair.produceevent(j);
  //Now store this to a vector defined earlier with numberofevents. We will do the writing out later.
  allevents.push_back(thisevent);
  }//while loop
	
  break;
  }//scoping the case
  case StarlightConstants::SINGLEMESON:
  {
  Gammagammasingle singlemeson(InputInformation,bbs1);
  //singlemeson.read();
  //singlemeson.singlecrosssection();
  cout<<"Single Meson In The Codeee"<<endl;	
  while(j<InputInformation.getnumberofevents())
  {
  thisevent=singlemeson.produceevent(j);
  allevents.push_back(thisevent);
  }//while loop

  break;
  }//scoping case, do this so the destruct on singlemeson object is activated
  case StarlightConstants::NARROWVMDEFAULT:
  {
  Gammaanarrowvm narrowvm(InputInformation,bbs1);
  while(j<InputInformation.getnumberofevents())
  {
  thisevent=narrowvm.produceevent(j);
  allevents.push_back(thisevent);
  }//while loop
  break;
  }

  case StarlightConstants::WIDEVMDEFAULT:
  {
  Gammaawidevm widevm(InputInformation,bbs1);
  while(j<InputInformation.getnumberofevents())
  {
  thisevent=widevm.produceevent(j);
  allevents.push_back(thisevent);
  }//while loop
  break;
  }
  case StarlightConstants::PSIFAMILY:
  {
  Psifamily psi(InputInformation,bbs1);
  while(j<InputInformation.getnumberofevents())
  {
  thisevent=psi.produceevent(j);
  allevents.push_back(thisevent);
  }//while loop
  break;
  }

  default:
  cout <<"Welcome to default mode"<<endl;
  }//mode switch
*/
    ofstream outputfile;
    outputfile.open("slight.out");
    int numberoftracks = 0;
    int eventnumber=1;
    for(vector<StarlightConstants::event>::iterator VD = allevents.begin();VD!=allevents.end();VD++){
	numberoftracks=(*VD).numberoftracks;
	outputfile <<"EVENT: "<<eventnumber<<" "<<numberoftracks<<" "<<1<<endl;
	if(PythiaOutput==false){

        	outputfile <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<endl;

		for(int i=0;i<numberoftracks;i++){
	    		outputfile << "TRACK: " <<" "<<jtog((*VD).charge[i]*(*VD).fsparticle[i])<<" "<<(*VD).px[i]<<" "<<(*VD).py[i]<<" "<<(*VD).pz[i]<<" "<<eventnumber<<" "<<i<<" "<<0<<" "<<(*VD).charge[i]*(*VD).fsparticle[i]<<endl;
		}//end of track for loop
	}//if pythiafalse loop
//In the process of making a pythia output function....checking into how to handle the verticies.
	if(PythiaOutput==true){

		for(int t=0;t<(*VD).numberofvertices;t++){
        		outputfile <<"VERTEX: "<<(*VD).vertx[t]/10.<<" "<<(*VD).verty[t]/10.<<" "<<(*VD).vertz[t]/10.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<endl; //convert from mm to cm for Geant
		}

		for(int i=0;i<numberoftracks;i++){
                        outputfile << "TRACK: " <<" "<<jtog((*VD).charge[i]*(*VD).fsparticle[i])<<" "<<(*VD).px[i]<<" "<<(*VD).py[i]<<" "<<(*VD).pz[i]<<" "<<eventnumber<<" "<<i<<" "<<(*VD).mother1[i]<<" "<<(*VD).mother2[i]<<" "<<(*VD).daughter1[i]<<" "<<(*VD).daughter2[i]<<" "<<(*VD).charge[i]*(*VD).fsparticle[i]<<endl;
                }//end of track for loop
	}//if pythiatrue loop

	eventnumber++;
    }//end of vector for loop
    outputfile.close();
	
    TFile g("slightNtuple.root", "RECREATE");
    
    TNtuple outntuple("outntuple","Data from Starlight simulation", "eventNumber:numberOfTracks:jtogFsparticle:px:py:pz:psparticle");
    /* for(vector<StarlightConstants::event>::iterator VD = allevents.begin();VD!=allevents.end();VD++){
	outdata.Fill();
	
	}*/
    
//    TTree outdata("outdata","Data from Starlight simulation");
//    TBranch *outbranch = outdata.Branch("outbranch", &);
    
    //Rootfile for the outdata
    TFile f("slight.root", "RECREATE");    
    
    //TTree for the outdata
    TTree outdata("outdata", "Data from Starlight simulation", 2);
    //TClonesArray for the tracks of each event.
    TClonesArray *particleArray = new TClonesArray("TParticle", 4);
    //Reference to particleArray
    TClonesArray &ar = *particleArray;
    //Making a branch in the tree
    outdata.Branch("branch", &particleArray);
    
    eventnumber = 1;
    
    //Looping over the events
    for(vector<StarlightConstants::event>::iterator VD = allevents.begin();VD!=allevents.end();VD++){

	numberoftracks = (*VD).numberoftracks;
	ar.Clear();
	
	//Looping over the tracks
	for(int i=0;i<numberoftracks;i++){
	  
	  //Calculating the total energy of the particle
	  Double_t mass = Masses::getMass((*VD).charge[i]*(*VD).fsparticle[i]);
	  Double_t etot = sqrt(mass*mass + (*VD).px[i]*(*VD).px[i] + (*VD).py[i]*(*VD).py[i] + (*VD).pz[i]*(*VD).pz[i]);
	  
	  //Creating a new TParticle for each track and putting it in the TClonesArray
	  new (ar[i]) TParticle((*VD).charge[i]*(*VD).fsparticle[i], -1, -1, -1, -1, -1, (*VD).px[i], (*VD).py[i], (*VD).pz[i],
				etot , 0, 0, 0, -1);   
	  
	  outntuple.Fill(eventnumber, numberoftracks,jtog((*VD).charge[i]*(*VD).fsparticle[i]),
			 (*VD).px[i],(*VD).py[i], (*VD).pz[i], (*VD).charge[i]*(*VD).fsparticle[i]);
	  
	}
	eventnumber++;
	//Filling the TTree for each event. 
	outdata.Fill();
    }
    outdata.Print();
    
    outdata.Write();
    g.Write();
    
    
        //Reading out the time at the finish
        currenttime=0;
        time(&currenttime);
        readabletime=localtime(&currenttime);
        cout<<"Time finished: "<<asctime(readabletime)<<endl;
    
    
    cout <<"Done?"<<endl;
}

void checkforluminositytable(bool &flag)
{
    cout<<"ISEED: "<<InputInformation.getseed()<<endl;
    ifstream wylumfile;
    wylumfile.precision(15);
    wylumfile.open("slight.txt");
    int Z1test=0,A1test=0,Z2test=0,A2test=0,numwtest=0,numytest=0,gg_or_gPtest=0,ibreakuptest=0,iinterferetest=0,NPTtest=0,in_or_cotest=0;
    double Gammatest=0.,Wmaxtest=0.,Wmintest=0.,Ymaxtest=0.,xinterferetest=0.,ptmaxtest=0.,bfordtest=0.;
    flag = false;
    bool b;
    wylumfile >> Z1test;
    wylumfile >> A1test;
    wylumfile >> Z2test;
    wylumfile >> A2test;
    wylumfile >> Gammatest;
    wylumfile >> Wmaxtest;
    wylumfile >> Wmintest;
    wylumfile >> numwtest;
    wylumfile >> Ymaxtest;
    wylumfile >> numytest;
    wylumfile >> gg_or_gPtest;
    wylumfile >> ibreakuptest;
    wylumfile >> iinterferetest;
    wylumfile >> xinterferetest;
    wylumfile >> in_or_cotest;
    wylumfile >> bfordtest;
    wylumfile >> ptmaxtest;
    wylumfile >> NPTtest;
    wylumfile.close();
    /* if( !(((Beam1.getZin()) == Z1test)&&((Beam1.getAin()) == A1test)&&((Beam2.getZin()) == Z2test)&&((Beam2.getAin()) == A2test)&&((inputbbs.gamma_em) == Gammatest)&&((inputbbs.numw) == numwtest)&&((inputbbs.Wmin) == Wmintest)&&((inputbbs.Ymax) == Ymaxtest)&&((inputbbs.numy) == numytest)&&((inputbbs.gg_or_gP) == gg_or_gPtest)&&((inputbbs.ibreakup) == ibreakuptest)&&((inputbbs.iinterfere) == iinterferetest)&&((inputbbs.xinterfere) == xinterferetest)&&((inputbbs.bford) == bfordtest)&&((inputbbs.in_or_co) == in_or_cotest)&&((inputbbs.ptmax) == ptmaxtest)&&((inputbbs.NPT) == NPTtest))){// == false){
     */
    if( !(((InputInformation.getZ1()) == Z1test)&&((InputInformation.getA1()) == A1test)&&((InputInformation.getZ2()) == Z2test)&&((InputInformation.getA2()) == A2test)&&((InputInformation.getgamma_em()) == Gammatest)&&((InputInformation.getnumw()) == numwtest)&&((InputInformation.getWmin()) == Wmintest)&&((InputInformation.getYmax()) == Ymaxtest)&&((InputInformation.getnumy()) == numytest)&&((InputInformation.getgg_or_gP()) == gg_or_gPtest)&&((InputInformation.getbreakupmode()) == ibreakuptest)&&((InputInformation.getinterferencemode()) == iinterferetest)&&((InputInformation.getinterferencepercent()) == xinterferetest)&&((InputInformation.getbford()) == bfordtest)&&((InputInformation.getincoherentorcoherent()) == in_or_cotest)&&((InputInformation.getmaximuminterpt()) == ptmaxtest)&&((InputInformation.getNPT()) == NPTtest))){// == false){


        cout<<"INSIDE THE LUMIN IF LOOP"<<endl;
        //okay, if we are in this loop, it means the input parameters are different than the one's used to create the last set of luminosity tables
        //Now lets create a new set
        flag=true;

    }



}

int jtog(int particleid)
{
//this function takes a jetset particle number and converts it to GEANT
    int jtogx =0;
    if(particleid==22) jtogx = 1;
    if(particleid==-11) jtogx = 2;
    if(particleid==11) jtogx = 3;
    if(particleid==12||particleid==-12||particleid==14||particleid==-14||particleid==16||particleid==-16) jtogx = 4;
    if(particleid==13) jtogx = 5;
    if(particleid==-13) jtogx = 6;
    if(particleid==111) jtogx = 7;
    if(particleid==211) jtogx = 8;
    if(particleid==-211) jtogx = 9;
    if(particleid==130) jtogx = 10;
    if(particleid==321) jtogx = 11;
    if(particleid==-321) jtogx = 12;
    if(particleid==2112) jtogx = 13;
    if(particleid==2212) jtogx = 14;
    if(particleid==-2212) jtogx = 15;
    if(particleid==310) jtogx = 16;
    if(particleid==-310) jtogx = 10;
    if(particleid==221) jtogx = 17;
    if(particleid==3122) jtogx = 18;
    if(particleid==3222) jtogx = 19;
    if(particleid==3212) jtogx = 20;
    if(particleid==3112) jtogx = 21;
    if(particleid==3322) jtogx = 22;
    if(particleid==3312) jtogx = 23;
    if(particleid==3334) jtogx = 24;
    if(particleid==-2112) jtogx = 25;
    if(particleid==-3122) jtogx = 26;
    if(particleid==-3222) jtogx = 27;
    if(particleid==-3212) jtogx = 28;
    if(particleid==-3112) jtogx = 29;
    if(particleid==-3322) jtogx = 30;
    if(particleid==-3312) jtogx = 31;
    if(particleid==-3334) jtogx = 32;
    return jtogx;
}

Eventchannel* CreateEventChannel(Inputparameters& input,Beambeamsystem& bbs)
{
    switch(input.getpidtest()){
    case StarlightConstants::ELECTRON:
    case StarlightConstants::MUON:
    case StarlightConstants::TAUON:
    {	Gammagammaleptonpair* leptonpair = new Gammagammaleptonpair(input,bbs);
	    return leptonpair;
	    break;
    }
    case StarlightConstants::A2://jetset
    case StarlightConstants::ETA://jetset
    case StarlightConstants::ETAPRIME://jetset
    case StarlightConstants::ETAC://jetset
    case StarlightConstants::F0://jetset
    {
	PythiaOutput=true;
	//This way we can output mother and daughter listings.
    }
    case StarlightConstants::F2:
    case StarlightConstants::F2PRIME:
    case StarlightConstants::ZOVERZ03:
    {
	Gammagammasingle* singlemeson = new Gammagammasingle(input,bbs);
	return singlemeson;
	break;
    }
    case StarlightConstants::RHO:
    case StarlightConstants::RHOZEUS:
    case StarlightConstants::OMEGA://Will probably be three body
    case StarlightConstants::PHI:
    case StarlightConstants::UPSILON:
    case StarlightConstants::UPSILON2S:
    case StarlightConstants::UPSILON3S:
    {
	if(input.getinteractiontest()==2){
	Gammaanarrowvm* narrowvm= new Gammaanarrowvm(input,bbs);
	return narrowvm;
	break;
	}

	if(input.getinteractiontest()==3){
        Gammaawidevm* widevm = new Gammaawidevm(input,bbs);
        return widevm;
        break;
	}
	cout<<"Please go back and adjust gg_or_gp to 2 or 3 for a VM, main.cpp"<<endl;
	return 0;
    }
    case StarlightConstants::JPSI:
    case StarlightConstants::JPSI2S:
    {
	Psifamily* Psi = new Psifamily(input,bbs);
	return Psi;
	break;
    }
    //rhoprime
    default:  cout<<"Hi and welcome to default event channel(null), main::CreateEventChannel"<<endl;
	return 0;//Maybe return empty Eventchannel object?
    }

}
