// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler

#include "Analysis.hh"
#include "Parameters.hh"
#include <ctime>
#include <fstream>


Analysis* Analysis::singleton = 0;

Analysis::Analysis()
{
    // Set file name for simulation output
    G4int threadNumber;
    std::ifstream inFile;
    inFile.open("threadNumber.txt");
    inFile >> threadNumber;
    inFile.close();

    sprintf(fileName,"outFile_%i.root", threadNumber);

    detectedPhotons = new TTree("DetectedPhotons","Photons detected at the detector strip.");

    detectedPhotons->Branch("runID", &runIDBuffer, "runID/F"); //modifed by me
    detectedPhotons->Branch("eventID", &eventIDBuffer, "eventID/F"); //modifed by me
    detectedPhotons->Branch("detNumb", &detNumbBuffer, "detNumb/F"); //modifed by me

    detectedPhotons->Branch("xPixel", &xPixelBuffer, "xPixel/F"); //modifed by me
    detectedPhotons->Branch("yPixel", &yPixelBuffer, "yPixel/F"); //modifed by me
    detectedPhotons->Branch("energy",&energyBuffer,"energy/F");
    detectedPhotons->Branch("wavelength",&wavelengthBuffer,"wavelength/F");
    detectedPhotons->Branch("localtime",&timeBuffer,"localtime/F");
    detectedPhotons->Branch("abstime",&absTimeBuffer,"abstime/F");
    detectedPhotons->Branch("length",&lengthBuffer,"length/F");
    detectedPhotons->Branch("x",&xBuffer,"x/F");
    detectedPhotons->Branch("y",&yBuffer,"y/F");
    detectedPhotons->Branch("z",&zBuffer,"z/F");
    detectedPhotons->Branch("px",&pxBuffer,"px/F");
    detectedPhotons->Branch("py",&pyBuffer,"py/F");
    detectedPhotons->Branch("pz",&pzBuffer,"pz/F");
    detectedPhotons->Branch("vertexX",&vertexXBuffer,"vertexX/F");
    detectedPhotons->Branch("vertexY",&vertexYBuffer,"vertexY/F");
    detectedPhotons->Branch("vertexZ",&vertexZBuffer,"vertexZ/F");
    detectedPhotons->Branch("vertexPx",&vertexPxBuffer,"vertexPx/F");
    detectedPhotons->Branch("vertexPy",&vertexPyBuffer,"vertexPy/F");
    detectedPhotons->Branch("vertexPz",&vertexPzBuffer,"vertexPz/F");
    detectedPhotons->Branch("gpsPosX",&gpsPositionX,"gpsPosX/F");
    detectedPhotons->Branch("gpsPosY",&gpsPositionY,"gpsPosY/F");
    detectedPhotons->Branch("gpsPosZ",&gpsPositionZ,"gpsPosZ/F");
    detectedPhotons->Branch("gpsDirX",&gpsDirectionX,"gpsPosX/F");
    detectedPhotons->Branch("gpsDirY",&gpsDirectionY,"gpsPosY/F");
    detectedPhotons->Branch("gpsDirZ",&gpsDirectionZ,"gpsPosZ/F");
    detectedPhotons->Branch("runId",&runIdBuffer,"runId/I");
    detectedPhotons->Branch("eventId",&eventIdBuffer,"eventId/I");
    detectedPhotons->Branch("trackId",&trackIdBuffer,"trackId/I");
    detectedPhotons->Branch("creatorProcess",&creatorProcessBuffer,"creatorProcess/I");
    detectedPhotons->Branch("parentId",&parentIdBuffer,"parentId/I");
    detectedPhotons->Branch("reflMirr",&reflMirrBuffer,"reflMirr/I");
    detectedPhotons->Branch("reflSurf",&reflSurfBuffer,"reflSurf/I");
    detectedPhotons->Branch("reflTotalCladClad",&reflTotalCladCladBuffer,"reflTotalCladClad/I");
    detectedPhotons->Branch("reflTotalCoreClad",&reflTotalCoreCladBuffer,"reflTotalCoreClad/I");
    detectedPhotons->Branch("reflFresnelCladClad",&reflFresnelCladCladBuffer,"reflFresnelCladClad/I");
    detectedPhotons->Branch("reflFresnelCoreClad",&reflFresnelCoreCladBuffer,"reflFresnelCoreClad/I");
    detectedPhotons->Branch("refracCladClad",&refracCladCladBuffer,"refracCladClad/I");
    detectedPhotons->Branch("refracCoreClad",&refracCoreCladBuffer,"refracCoreClad/I");
    detectedPhotons->Branch("rayleighScatterings",&rayleighScatteringsBuffer,"rayleighScatterings/I");
    detectedPhotons->Branch("lengthInCore",&lengthInCoreBuffer,"lengthInCore/F");
    detectedPhotons->Branch("lengthInInnerCladding",&lengthInInnerCladdingBuffer,"lengthInInnerCladding/F");
    detectedPhotons->Branch("lengthInOuterCladding",&lengthInOuterCladdingBuffer,"lengthInOuterCladding/F");

    trigger = new TTree("Trigger","Hits in the trigger.");
  
    trigger->Branch("runID", &runIDBuffer, "runID/F");
    trigger->Branch("eventID", &eventIDBuffer, "eventID/F");
    trigger->Branch("edep", &edepBuffer, "edep/F");
    trigger->Branch("xPos", &xPosBuffer, "xPos/F");
    trigger->Branch("yPos", &yPosBuffer, "yPos/F");
    trigger->Branch("zPos", &zPosBuffer, "zPos/F");
	
    energyTrack = new TTree("EnergyTrack", "Deposited energy and tracklength in core of fibre.");

    energyTrack->Branch("runID", &runIDBuffer, "runID/F");
    energyTrack->Branch("eventID", &eventIDBuffer, "eventID/F");
    energyTrack->Branch("edep", &energyBuffer, "edep/F");
    energyTrack->Branch("trackL", &lengthBuffer, "trackL/F");

    h_energy = new TH1F("h_energy","h_energy",200,0.0,1.0);
    h_trackL = new TH1F("h_trackL","h_trackL",400,0.0,2.0);

    primaryParticleTrack = new TTree("PrimaryParticleTrack", "Track of the primary particles");

    primaryParticleTrack->Branch("runID", &runIDBuffer, "runID/F");
    primaryParticleTrack->Branch("eventID", &eventIDBuffer, "eventID/F");
    primaryParticleTrack->Branch("xPos", &xPosBuffer, "xPos/F");
    primaryParticleTrack->Branch("yPos", &yPosBuffer, "yPos/F");
    primaryParticleTrack->Branch("zPos", &zPosBuffer, "zPos/F");

    initialParticle = new TTree("InitialParticle", "Initial Particle");
  
    initialParticle->Branch("runID", &runIDBuffer, "runID/F");
    initialParticle->Branch("eventID", &eventIDBuffer, "eventID/F");
    initialParticle->Branch("energy", &energyBuffer, "energy/F");
    initialParticle->Branch("xMom", &xMomBuffer, "xMom/F");
    initialParticle->Branch("yMom", &yMomBuffer, "yMom/F");
    initialParticle->Branch("zMom", &zMomBuffer, "zMom/F");
 


}


void Analysis::PrepareNewRun(const G4Run* aRun)
{
    runIdBuffer = aRun->GetRunID();
}


void Analysis::EndOfRun()
{
    dataFile = new TFile(fileName,"update");

    //detectedPhotons->Write("",TObject::kOverwrite);
    //trigger->Write("", TObject::kOverwrite);
    // primaryParticleTrack->Write("", TObject::kOverwrite);
    //initialParticle->Write("", TObject::kOverwrite);
    //   energyTrack->Write("", TObject::kOverwrite);
    h_energy->Write("", TObject::kOverwrite);
    h_trackL->Write("", TObject::kOverwrite);
    G4cout << "Data written to file " << fileName << G4endl;

    dataFile->Close();

    G4cout << "Data file closed: " << fileName << G4endl;

    delete dataFile;
}


void Analysis::Close()
{
    delete detectedPhotons;
    delete trigger;
    delete initialParticle;
    delete energyTrack;
}


void Analysis::PrepareNewEvent(const G4Event* anEvent)
{
    eventIdBuffer = anEvent->GetEventID();

    reflectionsAtMirror.resize(0);
    reflectionsAtFibreSurface.resize(0);
    reflectionsTotalAtCladCladInterface.resize(0);
    reflectionsTotalAtCoreCladInterface.resize(0);
    reflectionsFresnelAtCladCladInterface.resize(0);
    reflectionsFresnelAtCoreCladInterface.resize(0);

    refractionsAtCladCladInterface.resize(0);
    refractionsAtCoreCladInterface.resize(0);

    rayleighScatterings.resize(0);

    lengthInCore.resize(0);
    lengthInInnerCladding.resize(0);
    lengthInOuterCladding.resize(0);
}

void Analysis::FillInitialParticle(Double_t runID, Double_t eventID, Double_t energy, Double_t xMom, Double_t yMom, Double_t zMom)
{
    if(initialParticle != NULL)
    {
        runIDBuffer = runID;
        eventIDBuffer = eventID;
        energyBuffer = energy;
        xMomBuffer = xMom;
        yMomBuffer = yMom;
        zMomBuffer = zMom;
	//  initialParticle->Fill();
    }
    else
        G4cout << "Nullpointer to initialParticle tree." << G4endl;
}

void Analysis::FillTrigger(Double_t runID, Double_t eventID, Double_t edep, Double_t xPos, Double_t yPos, Double_t zPos)
{
    if(trigger != NULL)
    {
        runIDBuffer = runID;
        eventIDBuffer = eventID;
        edepBuffer = edep;
        xPosBuffer = xPos;
        yPosBuffer = yPos;
        zPosBuffer = zPos;

	//   trigger->Fill();
    }
    else
    {
        G4cout << "Nullpointer to trigger tree." << G4endl;
    }
}

void Analysis::FillEnergyTrack(Double_t runID, Double_t eventID, Double_t energy, Double_t trackL)
{
    if(energyTrack != NULL)
    {
        runIDBuffer = runID;
        eventIDBuffer = eventID;
        energyBuffer = energy;
        lengthBuffer = trackL;
	//        energyTrack->Fill();
	h_energy->Fill(energy);
	h_trackL->Fill(trackL);
    }
    else
    {
        G4cout << "Nullpointer to energyTrack tree" << G4endl;
    }
}

void Analysis::FillPrimaryParticleTrack(Double_t runID, Double_t eventID,
                                        Double_t xPos, Double_t yPos, Double_t zPos)
{
    if(primaryParticleTrack != NULL)
    {
        runIDBuffer = runID;
        eventIDBuffer = eventID;
        xPosBuffer = xPos;
        yPosBuffer = yPos;
        zPosBuffer = zPos;
	//  primaryParticleTrack->Fill();
    }
    else
    {
        G4cout << "Nullpointer to primaryParticleTrack tree" << G4endl;
    }
}

void Analysis::FillDetectedPhotons(Double_t runID, Double_t eventID, Double_t detNumb, Double_t xPixel,
                                   Double_t yPixel, Double_t energy,
                                   Double_t time, Double_t length, Double_t absTime,
                                   Double_t x, Double_t y, Double_t z,
                                   Double_t px, Double_t py, Double_t pz,
                                   Double_t vertexX, Double_t vertexY, Double_t vertexZ,
                                   Double_t vertexPx, Double_t vertexPy, Double_t vertexPz,
                                   Int_t trackId,
                                   Int_t creatorProcess, Int_t parentId)
{
    if(detectedPhotons!=NULL)
    {
        runIDBuffer   = runID;
        eventIDBuffer = eventID;
        detNumbBuffer = detNumb;
        xPixelBuffer = xPixel;
        yPixelBuffer = yPixel;

        energyBuffer = energy;
        timeBuffer = time;
        lengthBuffer = length;

        wavelengthBuffer = 1239.842/energy;

        absTimeBuffer = absTime;
        xBuffer = x;
        yBuffer = y;
        zBuffer = z;

        pxBuffer = px;
        pyBuffer = py;
        pzBuffer = pz;

        vertexXBuffer = vertexX;
        vertexYBuffer = vertexY;
        vertexZBuffer = vertexZ;

        vertexPxBuffer = vertexPx;
        vertexPyBuffer = vertexPy;
        vertexPzBuffer = vertexPz;

        trackIdBuffer = trackId;

        creatorProcessBuffer = creatorProcess;
        parentIdBuffer = parentId;

        reflMirrBuffer = reflectionsAtMirror[trackId];
        reflSurfBuffer = reflectionsAtFibreSurface[trackId];
        reflTotalCladCladBuffer = reflectionsTotalAtCladCladInterface[trackId];
        reflTotalCoreCladBuffer = reflectionsTotalAtCoreCladInterface[trackId];
        reflFresnelCladCladBuffer = reflectionsFresnelAtCladCladInterface[trackId];
        reflFresnelCoreCladBuffer = reflectionsFresnelAtCoreCladInterface[trackId];

        refracCladCladBuffer = refractionsAtCladCladInterface[trackId];
        refracCoreCladBuffer = refractionsAtCoreCladInterface[trackId];

        rayleighScatteringsBuffer = rayleighScatterings[trackId];

        lengthInCoreBuffer = lengthInCore[trackId];
        lengthInInnerCladdingBuffer = lengthInInnerCladding[trackId];
        lengthInOuterCladdingBuffer = lengthInOuterCladding[trackId];

	//   detectedPhotons->Fill();
    }
    else
    {
        G4cout << G4endl << "Pointer to detectedPhotons tree is NULL!" << G4endl;
    }
}

char* Analysis::FileName()
{
    return fileName;
}

void Analysis::SetGpsPosition(G4ThreeVector position)
{
    gpsPositionX = position[0];
    gpsPositionY = position[1];
    gpsPositionZ = position[2];
}

void Analysis::IncreaseReflectionsAtMirror(Int_t trackId)
{
    reflectionsAtMirror[trackId]++;
}


void Analysis::SetGpsDirection(G4ThreeVector direction)
{
    gpsDirectionX = direction[0];
    gpsDirectionY = direction[1];
    gpsDirectionZ = direction[2];
}

void Analysis::IncreaseReflectionsAtFibreSurface(Int_t trackId)
{
    reflectionsAtFibreSurface[trackId]++;
}


void Analysis::IncreaseTotalReflectionsAtCladCladInterface(Int_t trackId)
{
    reflectionsTotalAtCladCladInterface[trackId]++;
}


void Analysis::IncreaseTotalReflectionsAtCoreCladInterface(Int_t trackId)
{
    reflectionsTotalAtCoreCladInterface[trackId]++;
}


void Analysis::IncreaseFresnelReflectionsAtCladCladInterface(Int_t trackId)
{
    reflectionsFresnelAtCladCladInterface[trackId]++;
}


void Analysis::IncreaseFresnelReflectionsAtCoreCladInterface(Int_t trackId)
{
    reflectionsFresnelAtCoreCladInterface[trackId]++;
}


void Analysis::IncreaseRefractionsAtCladCladInterface(Int_t trackId)
{
    refractionsAtCladCladInterface[trackId]++;
}


void Analysis::IncreaseRefractionsAtCoreCladInterface(Int_t trackId)
{
    refractionsAtCoreCladInterface[trackId]++;
}


void Analysis::IncreaseRayleighScatterings(Int_t trackId)
{
    rayleighScatterings[trackId]++;
}


void Analysis::IncreaseReflectionRefractionAndScatteringVectors(Int_t trackId)
{
    for(int i=reflectionsAtMirror.size(); i<trackId+1; i++)
        reflectionsAtMirror.push_back(0);

    for(int i=reflectionsAtFibreSurface.size(); i<trackId+1; i++)
        reflectionsAtFibreSurface.push_back(0);

    for(int i=reflectionsTotalAtCladCladInterface.size(); i<trackId+1; i++)
        reflectionsTotalAtCladCladInterface.push_back(0);

    for(int i=reflectionsTotalAtCoreCladInterface.size(); i<trackId+1; i++)
        reflectionsTotalAtCoreCladInterface.push_back(0);

    for(int i=reflectionsFresnelAtCladCladInterface.size(); i<trackId+1; i++)
        reflectionsFresnelAtCladCladInterface.push_back(0);

    for(int i=reflectionsFresnelAtCoreCladInterface.size(); i<trackId+1; i++)
        reflectionsFresnelAtCoreCladInterface.push_back(0);

    for(int i=refractionsAtCladCladInterface.size(); i<trackId+1; i++)
        refractionsAtCladCladInterface.push_back(0);

    for(int i=refractionsAtCoreCladInterface.size(); i<trackId+1; i++)
        refractionsAtCoreCladInterface.push_back(0);

    for(int i=rayleighScatterings.size(); i<trackId+1; i++)
        rayleighScatterings.push_back(0);
}


void Analysis::IncreaseLengthInCore(Int_t trackId, Float_t lengthValue)
{
    lengthInCore[trackId] += lengthValue;
}


void Analysis::IncreaseLengthInInnerCladding(Int_t trackId, Float_t lengthValue)
{
    lengthInInnerCladding[trackId] += lengthValue;
}


void Analysis::IncreaseLengthInOuterCladding(Int_t trackId, Float_t lengthValue)
{
    lengthInOuterCladding[trackId] += lengthValue;
}


void Analysis::IncreaseLengthVectors(Int_t trackId)
{
    for(int i=lengthInCore.size(); i<trackId+1; i++)
        lengthInCore.push_back(0);

    for(int i=lengthInInnerCladding.size(); i<trackId+1; i++)
        lengthInInnerCladding.push_back(0);

    for(int i=lengthInOuterCladding.size(); i<trackId+1; i++)
        lengthInOuterCladding.push_back(0);
}