/**
 * @file SolidEventAction.cc
 * @author: (modified by) Ibrahin Pinera
 * @date 2016 SoLid - University of Antwerp
 * @author: (modified by) Mariangela Settimo
 * @date 2017 SoLid - Subatech
 */


#include "SolidEventAction.hh"
#include "SolidRunAction.hh"
#include "SolidEventMessenger.hh"
#include "SolidTrajectory.hh"
#include "SolidSteppingAction.hh"
#include "SolidPrimaryGeneratorAction.hh"
#include "SolidUtils.hh"

//#ifdef G4MULTITHREADED
//#include "G4MTRunManager.hh"
//#else
#include "G4RunManager.hh"
//#endif
#include "G4SteppingManager.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4HCofThisEvent.hh"
#include "G4SystemOfUnits.hh"
#include "G4VHitsCollection.hh"
#include "G4TrajectoryContainer.hh"
#include "G4Trajectory.hh"
#include "G4VVisManager.hh"
#include "G4SDManager.hh"
#include "G4UImanager.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"
#include "G4PrimaryVertex.hh"


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
SolidEventAction::SolidEventAction(SolidRunAction *ra, SolidPrimaryGeneratorAction* GePGA):
run_action(ra), genAction(GePGA), fVerbose(true)
{

    pvtCollID=-1;
    /*  const unsigned int nParticleTypes=12;
	std::string particleTypes[nParticleTypes] = {
	  "Neutron",		"Gamma",
	  "Electron", 	        "Positron",
	  "Alpha", 		"Triton",
	  "MuPlus",		"MuMinus",
	  "Proton",             "C12",
          "Deuteron",
	  "Tot"
	};
	const unsigned int nSDTypes=4;
	std::string SDTypes[nSDTypes] = {
	  "Edep", "Time", "Track", "TrackID"
	};
    */
	const unsigned int nSDVolumes=1;
	std::string SDVolumes[nSDVolumes] = {"ScintSD"};

	std::string SDName;
	
	G4cout << "Sensitive detectors:" << G4endl;
	
    for (unsigned int idSDVolume=0; idSDVolume < nSDVolumes; idSDVolume++)
    {
	//	for (unsigned int idSDType=0; idSDType < nSDTypes; idSDType++){
	//	for (unsigned int idParticle=0; idParticle < nParticleTypes; idParticle++){
	  SDName = "SolidDet/" + SDVolumes[idSDVolume];
	  //	  SDName = SDVolumes[idSDVolume]+"/"+SDTypes[idSDType]+particleTypes[idParticle];
	   G4cout <<"   "<< SDName << G4endl;
	   SDNames.push_back(SDName);
		  //	}
	       //	}
	}
    G4cout << "Number of sensitive detectors: " << SDNames.size() << G4endl;
      
    fgeEvMess = new SolidEventMessenger(this);
    
}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
SolidEventAction::~SolidEventAction()
{;}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void SolidEventAction::BeginOfEventAction(const G4Event* evt)
{
    
    event_id = evt->GetEventID();
//    G4cout << "######################################" << G4endl;
//    G4cout << "# Begin of event: " << event_id << G4endl;
//    G4cout << "######################################" << G4endl;

    Statistics.Evt = event_id;

    run_action->triggerTime = 0.;
    
    run_action->ThermalRadius = 0.;
    run_action->ModDistance = 0.;
    run_action->CaptureDistance = 0.;
//  run_action->LiCheck = 0;
    run_action->istherm = 0;
    run_action->NScatter = 0;
    run_action->NeutronStepCount = 0;
    run_action->PositronStepCount = 0;
    run_action->CapEneLi = 0.;
    run_action->CapEnePVT = 0.;
    run_action->CaptureVol = 0;
    run_action->CaptureCube = -10000;
    run_action->NeutronDeltaR = -10.;

    run_action->InitPosX = 0.;
    run_action->InitPosY = 0.;
    run_action->InitPosZ = 0.;
    run_action->PositronEndX = 0.;
    run_action->PositronEndY = 0.;
    run_action->PositronEndZ = 0.;
    run_action->PositronDistance = 0.;
    run_action->AnnihilationTime = 0.;
    run_action->AnnihilationCube = -10000;
    run_action->PositronDeltaR = -10.;
    run_action->AnnihilationVol = 0;

    run_action->ibdDeltaT = -1.e4;
    run_action->ibdDeltaR = -10.;


    // Vectors need to be cleared!

    // T1
    run_action->primary_pdgid.clear();
    run_action->primary_energy.clear();
    run_action->primary_posx.clear();
    run_action->primary_posy.clear();
    run_action->primary_posz.clear();
    run_action->primary_momx.clear();
    run_action->primary_momy.clear();
    run_action->primary_momz.clear();
    run_action->primary_vol.clear();
    run_action->primary_cube.clear();
    run_action->primary_cubex.clear();
    run_action->primary_cubey.clear();
    run_action->primary_cubez.clear();

    // T5
    run_action->cubenumber.clear();
    run_action->cubeXPosition.clear();
    run_action->cubeYPosition.clear();
    run_action->cubeZPosition.clear();
    run_action->pdg.clear();
    run_action->partID.clear();
    run_action->mothID.clear();
    run_action->edepTot.clear();
    run_action->edepPvt.clear();
    run_action->edepLi.clear();
    run_action->timeOfDemise.clear();
    run_action->distanceTravelled.clear();

    // T6
    run_action->partPDG.clear();
    run_action->partTrackID.clear();
    run_action->partParentID.clear();
    run_action->partEnergy.clear();
    run_action->partDeposited.clear();
    run_action->partX.clear();
    run_action->partY.clear();
    run_action->partZ.clear();
    run_action->partTime.clear();
    run_action->partCubeID.clear();
    run_action->partCubeX.clear();
    run_action->partCubeY.clear();
    run_action->partCubeZ.clear();    
    run_action->partVolID.clear();

    // T7
    run_action->checkEnteredDetector = false;
    run_action->checkTrackID = 1;
    run_action->enteringPartID.clear();
    run_action->enteringTrackID.clear();
    run_action->enteringParentID.clear();
    run_action->enteringEnergy.clear();
//    run_action->enteringDeposited.clear();
    run_action->enteringX.clear();
    run_action->enteringY.clear();
    run_action->enteringZ.clear();
    run_action->enteringXmom.clear();
    run_action->enteringYmom.clear();
    run_action->enteringZmom.clear();
    run_action->enteringTime.clear();
    run_action->enteringPreVolID.clear();
    run_action->enteringPostVolID.clear();

    // T8
    run_action->nGeneratedTrackID = 0;
    run_action->nGenerationTrackID.clear();
    run_action->nGenerationParentID.clear();
    run_action->nGenerationEnergy.clear();
    run_action->nGenerationPosX.clear();
    run_action->nGenerationPosY.clear();
    run_action->nGenerationPosZ.clear();
    run_action->nGenerationTime.clear();
    run_action->nGenerationCube.clear();
    run_action->nGenerationVol.clear();
    run_action->nGenerationProcess.clear();
    run_action->particleStepCount = 0;

        // T9
    run_action->summaryAllPartID.clear();
    run_action->summaryAllPDG.clear();
    run_action->summaryPartPDG.clear(); 
    run_action->summaryPartID.clear(); 
    run_action->summaryParentID.clear();
    run_action->summaryParentPDG.clear();
    run_action->summaryNsteps.clear(); 
    run_action->summaryNcubes.clear(); 
    //  run_action->summaryNCompton.clear(); 
    //  run_action->summaryNIonisation.clear(); 
    //  run_action->summaryNBrem.clear(); 
    run_action->summaryInitialEkin.clear(); 
    run_action->summaryFirstVol.clear();  
    run_action->summaryLastVol.clear(); 
    //    run_action->summaryCreatorProcess.clear();
    //  run_action->summaryLastProcess.clear(); 
    run_action->summaryCreatorProcessID.clear(); 
    run_action->summaryLastProcessID.clear();
    run_action->summaryGoodEntry.clear();
     run_action->summaryPosX.clear();
     run_action->summaryPosY.clear();
     run_action->summaryPosZ.clear();
    
    run_action->listCubes.clear(); 


    IsEdep0Set = false;

    Statistics.Deposit = 0;
    Statistics.Generated = 0;
    Statistics.Generated1 = 0;
    Statistics.Generated2 = 0;
    Statistics.WLS = 0;
    Statistics.BulkAbs = 0;
    Statistics.Detected = 0;
    Statistics.LiCapture = 0;
    Statistics.PVTCapture = 0;
    Statistics.CapCube = 0;
    Statistics.CapEne = 0;
   
    //Position Reset
    for(int i=0; i<3; i++)
    {
        Statistics.CapPos[i] = -650.; //so neutrons that are not captured are represented separately away from those which are captured
        Statistics.ThermPos[i] = 0.;
    }
    for(int i=0; i<Statistics.NCubes; i++)
    {
        Statistics.EneCube[i] = 0.;
    }

      
  if ( pvtCollID == -1 )
  {
    G4SDManager *SDman = G4SDManager::GetSDMpointer();
    pvtCollID = SDman->GetCollectionID("pvtCollection");
  }
  
}

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

G4THitsMap<G4double>*
SolidEventAction::GetHitsCollection(G4int hcID,
                                    const G4Event* event) const
{
  G4THitsMap<G4double>* hitsCollection
    = static_cast<G4THitsMap<G4double>*>(
        event->GetHCofThisEvent()->GetHC(hcID));

  if ( ! hitsCollection )
  {
    G4ExceptionDescription msg;
    msg << "Cannot access hitsCollection ID " << hcID;
    G4Exception("B4dEventAction::GetHitsCollection()",
      "MyCode0003", FatalException, msg);
  }

  return hitsCollection;
}

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

G4double SolidEventAction::GetSum(G4THitsMap<G4double>* hitsMap) const
{
  hitsMap->PrintAllHits();
  G4double sumValue = 0.;
  std::map<G4int, G4double*>::iterator it;
  G4int counter = 0;
  for ( it = hitsMap->GetMap()->begin(); it != hitsMap->GetMap()->end(); it++)
  {
	  counter++;
      //	  G4cout << "Counter: " << counter << " Id: " << it->first << "Value: " << *(it->second) << G4endl;
      sumValue += *(it->second);
  }
  return sumValue;
}

void SolidEventAction::WriteHit( G4int cubeId, G4String layer, G4int pdg, G4double energy, G4double track, G4double time, G4int pID, G4int mID)
{

      G4int cubeZPosition = cubeId/10000;
      G4int cubeYPosition = cubeId%100;
      G4int cubeXPosition = (cubeId%10000)/100;

	  run_action->cubenumber.push_back(cubeId);
	  run_action->cubeXPosition.push_back(cubeXPosition);
	  run_action->cubeYPosition.push_back(cubeYPosition);
	  run_action->cubeZPosition.push_back(cubeZPosition);
      run_action->pdg.push_back(pdg);
	  run_action->partID.push_back(pID);
	  run_action->mothID.push_back(mID);
	  run_action->timeOfDemise.push_back(time);
	  run_action->distanceTravelled.push_back(track);
	  
	  
      if ( layer == "PVT" )
      {
		  run_action->edepPvt.push_back(energy);
		  run_action->edepLi.push_back(0.0);
	  }
      if ( layer == "Li" )
      {
		  run_action->edepPvt.push_back(0.0);
		  run_action->edepLi.push_back(energy);
	  }
	  run_action->edepTot.push_back(energy);


}


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

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void SolidEventAction::EndOfEventAction(const G4Event* evt)
{



  ////////////////////////////////////////////////////
  // Get data from primary vertex (T1)              //
  ////////////////////////////////////////////////////
  G4int nPrimaryVertexes = evt->GetNumberOfPrimaryVertex ();
  
  for (G4int i=0; i<nPrimaryVertexes; i++)
  {
    G4int pdgid = evt->GetPrimaryVertex(i)->GetPrimary()->GetPDGcode();
    G4ThreeVector position = evt->GetPrimaryVertex(i)->GetPosition();
    G4ThreeVector momentum = evt->GetPrimaryVertex(i)->GetPrimary()->GetMomentum();
        
    G4double mass = evt->GetPrimaryVertex(i)->GetPrimary()->GetMass();
    G4double totMomentum = pow( pow( momentum[0],2) + pow(momentum[1],2) + pow(momentum[2],2), 0.5 );
    G4double kineticEnergy = pow( pow(totMomentum,2) + pow(mass,2), 0.5 ) - mass;

    G4int cube = run_action->InitCube;

    G4int cubex=-10, cubey=-10, cubez=-10;
    if( cube >= 0 )
    {
        cubez = cube/10000;
        cubex = (cube%10000)/100;
        cubey = cube%100;
    }

    run_action->primary_pdgid.push_back(pdgid);
    run_action->primary_energy.push_back(kineticEnergy);
    run_action->primary_posx.push_back(position[0]);
    run_action->primary_posy.push_back(position[1]);
    run_action->primary_posz.push_back(position[2]);
    run_action->primary_momx.push_back(momentum[0]/totMomentum);
    run_action->primary_momy.push_back(momentum[1]/totMomentum);
    run_action->primary_momz.push_back(momentum[2]/totMomentum);
    run_action->primary_vol.push_back(run_action->InitVol);
    //    run_action->primary_cube.push_back(run_action->InitCube);
    //    run_action->primary_cube.push_back(cube);
    run_action->primary_cubex.push_back(cubex);
    run_action->primary_cubey.push_back(cubey);
    run_action->primary_cubez.push_back(cubez);

  }
  
  if( run_action->AnnihilationCube >= 0 && run_action->CaptureCube >= 0 )
  {
      run_action->ibdDeltaT = run_action->CaptureTime - run_action->AnnihilationTime;
      run_action->ibdDeltaR = CubesDistance(run_action->AnnihilationCube, run_action->CaptureCube);
  }
  
  if( pvtCollID < 0 )
  {
    G4cout << " pvtCollId < 0" << G4endl;
    return; 
  }


  //=============================
   // address hits collections
  SolidPVTHitsCollection* SHC = NULL;
  //SolidPVTHitsCollection* SHC2 = NULL;
  
  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
  if(HCE)
  {
    SHC = (SolidPVTHitsCollection*)(HCE->GetHC(pvtCollID));

 
    // SHC2 = (SolidPVTHitsCollection*)(HCE->GetHC(liCollID));
    // G4cout << "Get collection liCollID " << liCollID << G4endl;  
  }

  /*
  for(unsigned int nameId = 0; nameId < SDNames.size(); nameId++) {
    std::string SDName = SDNames.at(nameId);
    int HCId= G4SDManager::GetSDMpointer()->GetCollectionID(SDName);
    hitmap_type* hitmap=GetHitsCollection(HCId,evt)->GetMap();
    hitcollection[SDName]=hitmap;
  }
  */
  
  
  //  std::map<G4int, G4double*>::iterator it;
  //////////////////////////////////////////////////////////////////////////
  // Read PVT sensitive detector information and
  // write it into the vectors in ROOT tree
  //////////////////////////////////////////////////////////////////////////
  if(SHC) 
    S_hits = SHC->entries();

  //  G4cout << "S_hits = " << S_hits << G4endl;
  G4double totalCountedEnergy=0.0;
  G4double totalCountedTrack=0.0;
  G4double totalCountedTime=0.0;

 
  G4double energy = 0;
  G4double track = 0;
  G4double time = 0;
  int klast=0;

//  int cubeId_last = -100;
//  int pId_last = -100;
  int ncount = 0;
  
  for (int i=0; i<S_hits; ++i)
  {
    G4int cubeId = (*SHC)[i]->GetCubeID();
    G4int pID = (*SHC)[i]->GetTrackID();
    G4int pdg = (*SHC)[i]->GetParticlePDG();
    G4int mID = (*SHC)[i]->GetMotherID();

    G4String layer = (*SHC)[i]->GetLayerName();
    //if (i>0 && (cubeId == cubeId_last && pID == pId_last)) continue;

    for (int k=0; k<S_hits; ++k)
    {
      
      if ((*SHC)[k]->GetCubeID() == cubeId && (*SHC)[k]->GetTrackID() == pID)
      {
          energy += (*SHC)[k]->GetEdep();
          time += (*SHC)[k]->GetTime();
          track += (*SHC)[k]->GetTrackLength();
//          cubeId_last = cubeId;
//          pId_last = pID;
          klast = k;
          ncount++;
      }
      
    }
    totalCountedEnergy += energy;
    totalCountedTrack += track;
    totalCountedTime += time;
    time /= (double)ncount;
    
    //      G4cout <<  " WriteHit("<<cubeId << " " << layer << " " << pdg << " " << energy << " " << track << " " << time << " " << pID << " " << mID << ") " << G4endl;
    WriteHit(cubeId, layer, pdg, energy, track, time, pID, mID);
    i = klast;
    energy = time = track = ncount = 0;
    
  }
   
    
    //Open run manager and update statistics
//      if(run_action->edepTot.size()>0)
//      if(run_action->nGenerationTrackID.size()>0)
      {
//	  G4cout <<"######################################"<< G4endl;
//	  G4cout <<"# END of event: "<< event_id << G4endl;
//	  G4cout <<"######################################"<< G4endl;

          SolidRunAction *runac = (SolidRunAction*)(G4RunManager::GetRunManager()->GetUserRunAction());
          runac->UpdateStatistics(Statistics);
	  
	  // Save trees to output root file each 10000 events
//	  if( event_id%1000 == 0 && event_id > 0 )
//	  {
//	      runac->SaveFile(event_id+1);
//	  }
      }

}

void SolidEventAction::CountDetected(RunTally)
{
    Statistics.Detected++;
    //BranchFiller = PhotonTally;
    //EventTree->Fill();
}