/**
 * @file SolidRunAction.cc
 * @author: (modified by) Ibrahin Pinera
 * @date 2016 SoLid - University of Antwerp
 */


#include "SolidRunAction.hh"
#include "SolidRunMessenger.hh"
#include "SolidSteppingAction.hh"

//#ifdef G4MULTITHREADED
//#include "G4MTRunManager.hh"
//#else
#include "G4RunManager.hh"
//#endif
#include "G4Run.hh"
#include "G4UImanager.hh"
#ifdef G4VIS_USE
#include "G4VVisManager.hh"
#include "G4VisExecutive.hh"
#endif
#include "G4ios.hh"
#include "Randomize.hh"

#include "TFile.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TTree.h"
#include "TH2I.h"
#include "TH1D.h"
#include "TAxis.h"
#include "TError.h"
#include "TDirectory.h"
#include "TROOT.h"


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

SolidRunAction::SolidRunAction(TFile* outp_)
{
     //
     outp = outp_;

     G4RunManager::GetRunManager()->SetPrintProgress(1000);
     fSeed = 123456789;


     m_storeT2 = false;
     m_storeT6 = false;
     m_storeT7 = false;
     m_storeT8 = false;
     m_storeT9 = true;

     m_raMess = new SolidRunMessenger(this);

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

SolidRunAction::~SolidRunAction()
{;}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void SolidRunAction::BeginOfRunAction(const G4Run* aRun)
{
    //start the timer clock to calculate run times
    start = time(NULL);
    //
    if (fSeed == -1) {
       //set the random seed to the CPU clock
       fSeed = (long) start;
    }
    G4Random::setTheSeed(fSeed);
    G4cout << "Random generator seed value=" 
	   << G4Random::getTheSeed() << G4endl;

    // set random seed for rand48() in SolidIBDKinematics.cc
    srand48(fSeed);

    // Go to output file
    outp->cd();
    
    // gROOT->ProcessLine("include <vector>");  //required for the vectors in T3 Branch

    //create the ROOT tree for event data
    
    T0 = new TTree("T0","Simulation setup information");

    T0->Branch(     "Event",           &Stats.Evt,               "Evt/I");
    T0->Branch(     "NumberCubes",     &Stats.NCubes,            "NCubes/I");
    T0->Branch(     "Seed",            &fSeed,                   "Seed/L"); 


    T1 = new TTree("T1","Primary vertex");

    T1->Branch(     "pdg",         &primary_pdgid);
    T1->Branch(     "energy",      &primary_energy);
    T1->Branch(     "posx",        &primary_posx);
    T1->Branch(     "posy",        &primary_posy);
    T1->Branch(     "posz",        &primary_posz);
    T1->Branch(     "momx",        &primary_momx);
    T1->Branch(     "momy",        &primary_momy);
    T1->Branch(     "momz",        &primary_momz);
    //    T1->Branch(     "cubeID",      &primary_cube);
    T1->Branch(     "cubex",       &primary_cubex);
    T1->Branch(     "cubey",       &primary_cubey);
    T1->Branch(     "cubez",       &primary_cubez);
    T1->Branch(     "volid",       &primary_vol);


    if( m_storeT2 )
    {
        T2 = new TTree("T2","ibd information");

        T2->Branch(     "LiCapture",       &Stats.LiCapture,         "LiCapture/I");
        T2->Branch(     "PVTCapture",      &Stats.PVTCapture,        "PVTCapture/I");
        T2->Branch(     "CaptureCube",         &Stats.CapCube,           "CaptureCube/I");
        //T2->Branch(     "DepositEnergy",   &Stats.Deposit,           "DepositE/F");
        T2->Branch(     "LiCaptureEnergy", &CapEneLi,                "LiCaptureEnergy/D");
        T2->Branch(     "PvtCaptureEnergy",&CapEnePVT,               "PvtCaptureEnergy/D");
	T2->Branch(     "CapPosX",         &Stats.CapPos[0],         "InitPosX/D");
	T2->Branch(     "CapPosY",         &Stats.CapPos[1],         "InitPosY/D");
	T2->Branch(     "CapPosZ",         &Stats.CapPos[2],         "InitPosZ/D");
        //    T2->Branch(     "ThermalPosX",     &Stats.ThermPos[0],       "ThermalPosX/D");
        //    T2->Branch(     "ThermalPosY",     &Stats.ThermPos[1],       "ThermalPosY/D");
        //    T2->Branch(     "ThermalPosZ",     &Stats.ThermPos[2],       "ThermalPosZ/D");
        //  T2->Branch(     "EneCube",         &Stats.EneCube,           "EneCube[NCubes]/F");
        //    T2->Branch(     "CapE",            &Stats.CapEne,            "CapE/F");
        //    T2->Branch(     "Detected",        &Stats.Detected,          "Detected/I");
        //    T2->Branch(     "Generated1",      &Stats.Generated1,        "Generated1/I");
        //    T2->Branch(     "Generated2",      &Stats.Generated2,        "Generated2/I");
        //    T2->Branch(     "WLS",             &Stats.WLS,               "WLS/I");
            T2->Branch(     "ModerationTime",  &ModTime,                 "ModerationTime/D");
        T2->Branch(     "ThermalTime",     &ThermalTime,             "ThermalTime/D");
        T2->Branch(     "CaptureTime",     &CaptureTime,             "CaptureTime/D");
        T2->Branch(     "ThermalRadius",   &ThermalRadius,           "ThermalRadius/D");
        T2->Branch(     "ModDistance",     &ModDistance,             "ModDistance/D");
        T2->Branch(     "CapDistance",     &CaptureDistance,         "CaptureDistance/D");
        T2->Branch(     "NScatter",        &NScatter,                "NScatter/I");
        //    T2->Branch(     "LiTotCross",      &LiTotCross,              "LiTotCross/I");
        //    T2->Branch(     "LiFastCross",     &LiFastCross,             "LiFastCross/I");
        //    T2->Branch(     "LiThermalCross",  &LiThermalCross,          "LiThermalCross/I");
        //    T2->Branch(     "FurthestDistance",&Furthest,                "FurthestDistance/D");
        //    T2->Branch(     "HDPEReflection",  &NeutronReflect,          "HDPEReflection/I");
        //    T2->Branch(     "NeutronEscape",   &NeutronEscape,           "NeutronEscape/I");
        T2->Branch(     "CaptureVol",   &CaptureVol,        "CaptureVol/I");
        T2->Branch(     "CaptureCube",  &CaptureCube,       "CaptureCube/I");
        T2->Branch(     "NeutronDeltaR",&NeutronDeltaR,     "NeutronDeltaR/D");

        T2->Branch(     "PositronEndX",    &PositronEndX,           "PositronEndX/D");
        T2->Branch(     "PositronEndY",    &PositronEndY,           "PositronEndY/D");
        T2->Branch(     "PositronEndZ",    &PositronEndZ,           "PositronEndZ/D");
        T2->Branch(     "PositronDistance",&PositronDistance,       "PositronDistance/D");
        T2->Branch(     "AnnihilationTime",&AnnihilationTime,       "AnnihilationTime/D");
        T2->Branch(     "AnnihilationVol", &AnnihilationVol,        "AnnihilationVol/I");
        T2->Branch(     "AnnihilationCube",&AnnihilationCube,       "AnnihilationCube/I");
        T2->Branch(     "PositronDeltaR",  &PositronDeltaR,         "PositronDeltaR/D");

        T2->Branch(     "ibdDeltaT",&ibdDeltaT,       "ibdDeltaT/D");
        T2->Branch(     "ibdDeltaR",&ibdDeltaR,       "ibdDeltaR/D");
    }


    T5 = new TTree("T5","Cube-level information from sensitive detectors");

    //   T5->Branch(     "cubeID",          &cubenumber);
    T5->Branch(     "pdg",             &pdg);
    T5->Branch(     "trackid",         &partID);
    T5->Branch(     "parentid",        &mothID);
    T5->Branch(     "cubex",           &cubeXPosition);
    T5->Branch(     "cubey",           &cubeYPosition);
    T5->Branch(     "cubez",           &cubeZPosition);
    T5->Branch(     "edep_tot",        &edepTot);
    T5->Branch(     "edep_pvt",        &edepPvt);
    T5->Branch(     "edep_zns",        &edepLi);
 // T5->Branch(     "time",            &timeOfBirth);
    T5->Branch(     "time_last",       &timeOfDemise);
    T5->Branch(     "tracklength",     &distanceTravelled);

    if( m_storeT6 )
    {
        T6 = new TTree("T6","particles tracking info");


        T6->Branch(     "pdg",        &partPDG);
        T6->Branch(     "trackid",    &partTrackID);
        T6->Branch(     "parentid",   &partParentID);
        T6->Branch(     "energy",     &partEnergy);
        T6->Branch(     "deltaE",     &partDeposited);
        T6->Branch(     "posx",       &partX);
        T6->Branch(     "posy",       &partY);
        T6->Branch(     "posz",       &partZ);
        T6->Branch(     "time",       &partTime);
        T6->Branch(     "cubex",      &partCubeX);
        T6->Branch(     "cubey",      &partCubeY);
        T6->Branch(     "cubez",      &partCubeZ);
        T6->Branch(     "volid",      &partVolID);
    }

    if( m_storeT7 )
    {
        T7 = new TTree("T7","particles entering the detector volume");

        T7->Branch(     "pdg",        &enteringPartID);
        T7->Branch(     "trackid",    &enteringTrackID);
        T7->Branch(     "parentid",   &enteringParentID);
        T7->Branch(     "energy",     &enteringEnergy);
//      T7->Branch(     "deltaE",     &enteringDeposited);
        T7->Branch(     "posx",       &enteringX);
        T7->Branch(     "posy",       &enteringY);
        T7->Branch(     "posz",       &enteringZ);
        T7->Branch(     "momx",       &enteringXmom);
        T7->Branch(     "momy",       &enteringYmom);
        T7->Branch(     "momz",       &enteringZmom);
        T7->Branch(     "time",       &enteringTime);
        T7->Branch(     "prevolid",   &enteringPreVolID);
        T7->Branch(     "postvolid",  &enteringPostVolID);
    }


    if( m_storeT8 )
    {
        T8 = new TTree("T8","spallation neutrons from muons");

        T8->Branch(     "trackid", &nGenerationTrackID);
        T8->Branch(     "parentid", &nGenerationParentID);
        T8->Branch(     "energy",  &nGenerationEnergy);
        T8->Branch(     "posx",    &nGenerationPosX);
        T8->Branch(     "posy",    &nGenerationPosY);
        T8->Branch(     "posz",    &nGenerationPosZ);
        T8->Branch(     "time",    &nGenerationTime);
        T8->Branch(     "cubeid",  &nGenerationCube);
        T8->Branch(     "volid",   &nGenerationVol);
        T8->Branch(     "process",   &nGenerationProcess);
    }
    if( m_storeT9 )  {
        T9 = new TTree("T9","statistical process info from MC");
	T9->Branch(     "pdg",        &summaryPartPDG);
        T9->Branch(     "trackid",    &summaryPartID);
        T9->Branch(     "parentid",   &summaryParentID);
	T9->Branch(     "parentPDG",   &summaryParentPDG);	
        T9->Branch(     "initialEkin",   &summaryInitialEkin);
        T9->Branch(     "nSteps", &summaryNsteps);
        T9->Branch(     "nCubes",   &summaryNcubes);
        T9->Branch(     "firstVolID",    &summaryFirstVol);
        T9->Branch(     "lastVolID",    &summaryLastVol);

	T9->Branch(     "CreatorProcessID",   &summaryCreatorProcessID);
	T9->Branch(     "LastProcessID",   &summaryLastProcessID);
        T9->Branch(     "posx",       &summaryPosX);
        T9->Branch(     "posy",       &summaryPosY);
        T9->Branch(     "posz",       &summaryPosZ);
	/* ----- 
         the branches below are kept here commented, while doing a validation of the output. 
         They can be restored if needed. 
	*/
	//  T9->Branch(     "nIonis",  &summaryNIonisation);
	//  T9->Branch(     "nCompton",   &summaryNCompton);
        //  T9->Branch(     "nBrem",   &summaryNBrem);
	//  T9->Branch(     "CreatorProcess",   &summaryCreatorProcess);
	//  T9->Branch(     "LastProcess",   &summaryLastProcess);
	}

    
    G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
//    G4UImanager* UI = G4UImanager::GetUIpointer();
//    UI->ApplyCommand("/tracking/storeTrajectory 1");
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void SolidRunAction::EndOfRunAction(const G4Run*)
{
    // Go to output file
    outp->cd();

    Stats.Evt++;
    T0->Fill();
    T0->Write("",TObject::kOverwrite);
    T1->Write("",TObject::kOverwrite);
    if( m_storeT2 )
        T2->Write("",TObject::kOverwrite);
//    T3->Write("",TObject::kOverwrite);
//    T4->Write("",TObject::kOverwrite);
    T5->Write("",TObject::kOverwrite);
    if( m_storeT6 )
        T6->Write("",TObject::kOverwrite);
    if( m_storeT7 )
        T7->Write("",TObject::kOverwrite);
    if( m_storeT8 )
        T8->Write("",TObject::kOverwrite);
    if( m_storeT9 )
      T9->Write("",TObject::kOverwrite);

    outp->Close();
    
    //display run time and write to file Runtime.out
    time_t end = time(NULL);
    G4int elapsed = end - start;
    G4cout<<"Run Completed in "<<elapsed/3600<<":"<<(elapsed%3600)/60<<":"<<((elapsed%3600)%60)<<G4endl;

    

//#ifdef G4VIS_USE
//    if (G4VVisManager::GetConcreteInstance()) {
//        G4VisManager *visMan = (G4VisManager *) G4VVisManager::GetConcreteInstance();
//        // Apply command only when current viewer is available
//        if(visMan->GetCurrentViewer()!=0) {
//          G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
//        }
//    }
//#endif
}

void SolidRunAction::UpdateTracking(void){
//  T3->Fill();
}   

void SolidRunAction::UpdateStatistics(RunTally aRunTally){
    Stats = aRunTally;
    
    T1->Fill();
    if( m_storeT2 )
        T2->Fill();
//    T3->Fill();
//    T4->Fill();
    T5->Fill();
    if( m_storeT6 )
        T6->Fill();
    if( m_storeT7 )
        T7->Fill();
    if( m_storeT8 )
        T8->Fill();
    if( m_storeT9 )
      T9->Fill();

    // printf("Moderation Length = %e mm, \t Moderation Time = %1.2f ns , \n", ModDistance, ModTime);
    // printf("Thermal Radius = %e mm, \t Thermal Neutron Time = %1.2f ns ,\n", ThermalRadius, ThermalTime);
    // printf("Capture Distance = %e, \t Capture Time = %1.2f ns \n", CaptureDistance, CaptureTime);
    // printf("Crosses before moderated = %d, \t \t Thermal Crosses = %d, \t Total = %d \n", LiFastCross, LiThermalCross, LiTotCross);
}

void SolidRunAction::SaveFile(G4int){
    Stats.Evt++;
    T0->Fill();
    T0->Write("",TObject::kOverwrite);
    T1->Write("",TObject::kOverwrite);
    if( m_storeT2 )
        T2->Write("",TObject::kOverwrite);
//    T3->Write("",TObject::kOverwrite);
//    T4->Write("",TObject::kOverwrite);
    T5->Write("",TObject::kOverwrite);
    if( m_storeT6 )
        T6->Write("",TObject::kOverwrite);
    if( m_storeT7 )
        T7->Write("",TObject::kOverwrite);
    if( m_storeT8 )
        T8->Write("",TObject::kOverwrite);
    if( m_storeT9 )
     T9->Write("",TObject::kOverwrite);

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....