#include "ELYSE/EventAction.hh"

//Geant4
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4EventManager.hh"
#include "G4UImanager.hh"
#include "G4TrajectoryContainer.hh"
#include "G4VVisManager.hh"
#include "G4ios.hh"
#include "globals.hh"
#include "G4ThreeVector.hh"
#include "G4TransportationManager.hh" 
#include "G4Navigator.hh" 
#include "G4SDManager.hh"
#include "G4DigiManager.hh"
#include "G4UnitsTable.hh"
#include "G4UIcmdWith3VectorAndUnit.hh"

//std
#include <set>
#include <iomanip>
#include <string>
#include <vector>
#include <stdlib.h> //abs

// AIDA :
#include <AIDA/ITree.h>
#include <AIDA/IManagedObject.h>
#include <AIDA/ITuple.h>

//ELYSE
#include "ELYSE/Cast.hh"
#include "ELYSE/Analysis.hh"
#include "ELYSE/Trajectory.hh"
#include "ELYSE/RunAction.hh"
#include "ELYSE/PrimaryGeneratorAction.hh"
//#include "ELYSE/TrappVolHit.hh"
#include "ELYSE/DetectorConstruction.hh"


ELYSE::EventAction::EventAction(ELYSE::Analysis& aAnalysis,
				ELYSE::RunAction* myRun, 
				ELYSE::DetectorConstruction* myDetector, 
				ELYSE::PrimaryGeneratorAction* myGenerator)
:fAnalysis(aAnalysis),
 runAction(myRun), 
 generatorAction(myGenerator), 
 detectorConstructor(myDetector),
 eventTuple(0) {
  
  //Get User Histo pointers
  AIDA::ITree* usrTree = aAnalysis.tree(); //get the tree
  if (!usrTree) {
    G4cout << "ELYSE::EventAction: cannot get Analysis Tree" << G4endl;
    return;
  }
  //see AIDA::IManagedObject a comment on dynamic_cast
  AIDA::IManagedObject* obj; //generic Object managed by a Tree

  obj = usrTree->find("Event");
  if(!obj) {
    G4cout << "EventAction: WARNING: no tuple Event" << G4endl;
    usrTree->ls();
    //exit(0);
  } else {
    eventTuple =  CAST(obj,AIDA::ITuple);
    if (!eventTuple) {
      G4cout << "EventAction: FATAL: eventTuple not a Tuple" << G4endl;
      usrTree->ls();
      exit(0);
    }
  }

}//Ctor

//-------------------------------------------------------------------------------------------

ELYSE::EventAction::~EventAction(){
}//Dtor

//-------------------------------------------------------------------------------------------

void ELYSE::EventAction::BeginOfEventAction(const G4Event* evt){

  G4cout << " (JEC) EventAction::Begin EventAction" << G4endl;
  if (!evt) evt->Print();

}//BeginOfEventAction

//-------------------------------------------------------------------------------------------

void ELYSE::EventAction::EndOfEventAction(const G4Event* evt) {
  G4cout << " (JEC) EventAction::End EventAction" << G4endl;
    
  //see Analysis.cxx to get the Tuple variables

  // --------------------
  //  Get Event Information
  // --------------------
  //JEC FIXME: difference between event_id &  vecRecNumber
  G4int         event_id = evt->GetEventID();                      

  G4ThreeVector vtx = generatorAction->GetVtx();                   
  G4int         vtxvol   = EventFindStartingVolume(vtx);

  G4cout << ">>>>>>>>>>>>> New Event ";
  G4cout << " evt Id " << event_id 
	 << " start in volume " << vtxvol
	 << G4endl;
  G4cout << "Vertex (" 
	 << vtx.x()/cm << " , "
	 << vtx.y()/cm << " , "
	 << vtx.z()/cm << " , "
	 << ")"
	 << G4endl;

  //JEC FIXME: introduce enumeration for the column shared by Analysis/EventAction &  namespace protected

  if(!eventTuple) return;

  eventTuple->fill(0, event_id);                                    //eventId
  eventTuple->fill(1, vtxvol);                                      //vtxVol

  AIDA::ITuple* vtxPos = eventTuple->getTuple( 2 );
  vtxPos->fill(0, vtx.x()/cm);                                       //vtxPos
  vtxPos->fill(1, vtx.y()/cm);
  vtxPos->fill(2, vtx.z()/cm);
  vtxPos->addRow();
  

  // --------------------
  //  Get Track Information
  // --------------------

  //variables used to fill the tuple (track part)
  G4int pId, parent; 
  G4float timeStart;
  G4double dx, dy, dz;
  G4double ETot, px, py, pz;
  G4int startVol, stopVol;

  //----------------
  // The beam: 1 particle for the moment
  //----------------
  AIDA::ITuple* track = eventTuple->getTuple( 4 );

  G4int  beampdg = generatorAction->GetBeamPDG();
  pId =  beampdg;                                                  //pId
  track->fill(0, pId);

  parent = 0;                                                      //parent (none)
  track->fill(1, parent);  
  
  timeStart = 0;                                                   //creation time
  track->fill(2, timeStart);

  G4ThreeVector beamdir    = generatorAction->GetBeamDir();        //direction
  AIDA::ITuple* direction = track->getTuple( 3 );
  dx = beamdir.x();
  dy = beamdir.y();
  dz = beamdir.z();

  direction->fill(0, dx);
  direction->fill(1, dy);
  direction->fill(2, dz);
  direction->addRow();

  
  G4double beamenergy = generatorAction->GetBeamEnergy();
  G4double mass =    G4ParticleTable::GetParticleTable()->FindParticle(beampdg)->GetPDGMass();
  
  track->fill(4, mass);

  ETot = beamenergy;                         
  G4double pTot = ETot;
  if ( (beampdg != 22)       && //gamma & Optical Photon (To be verified)
       (std::abs(beampdg) != 12) && //(anti)nu-e
       (std::abs(beampdg) != 14) && //(anti)nu-mu
       (std::abs(beampdg) != 16) ) { //(anti)nu-tau  
    pTot = ETot*ETot - mass*mass;
    if (pTot < 1e-6){ 
      pTot  = 0.;
    } else {
      pTot = std::sqrt(pTot);
    }
  }

  track->fill(5, pTot);  
  track->fill(6, ETot);
  

  AIDA::ITuple* momentum = track->getTuple( 7 );                           //momentum
  px = pTot * dx;
  py = pTot * dy;
  pz = pTot * dz;
  momentum->fill(0, px);
  momentum->fill(1, py);
  momentum->fill(2, pz);
  momentum->addRow();

  AIDA::ITuple* startPos = track->getTuple( 8 );                          //start position
  startPos->fill(0, vtx.x()/cm);
  startPos->fill(1, vtx.y()/cm);
  startPos->fill(2, vtx.z()/cm);
  startPos->addRow(); 

  AIDA::ITuple* stopPos = track->getTuple( 9 );                          //stop position
  stopPos->fill(0, vtx.x()/cm);
  stopPos->fill(1, vtx.y()/cm);
  stopPos->fill(2, vtx.z()/cm);
  stopPos->addRow(); 

  startVol = stopVol = -1;
  track->fill(10, startVol);                                   //Starting Volume
  track->fill(11, stopVol);                                   //Stopping Volume

  //add the Beam track to tuple
  track->addRow();

  G4cout << "----> Tk{Beam}: " 
	 << " pId " << pId
	 << " parent " << parent
	 << " creation time " << timeStart 
	 << " Volumes " << startVol << " " << stopVol << "\n"
	 << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n"
	 << " ETot " << ETot
	 << " Mass " << mass
	 << " pTot " << pTot
	 << " px,py,pz " << px << " " << py << " " << pz << "\n"
	 << G4endl;
  
  // --------------------------
  //  Loop over Trajectories
  // --------------------------

  G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer();

  G4int n_trajectories = 0;
  if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();

  G4int trjIndex;    //this is the trajectory Index
  G4int parentIndex; //this the parent trajectory index

  G4ThreeVector mom;        //px,py,pz
  G4double      mommag;     //p = sqrt(px*px + ...)
  G4ThreeVector stop;       //stopping point
  G4ThreeVector start;      //starting point
  G4String stopVolumeName;  //Volume name of the stopping point for retreive

  eventTuple->fill(5, n_trajectories);                                           //nPart
  G4cout << "Final # of tracks : " <<  n_trajectories << G4endl;


  for (G4int i=0; i < n_trajectories; i++) {
    Trajectory* trj = (Trajectory*)((*(evt->GetTrajectoryContainer()))[i]);
        
    pId      =  trj->GetPDGEncoding();                                   //pId
    trjIndex =  trj->GetTrackID();
      
    if ( ! trj->GetSaveFlag() ) continue;  //only consider trajectories which are flaged (see TrackingAction, SteppingAction)

    // initial point of the trajectory
    G4TrajectoryPoint* aa =   (G4TrajectoryPoint*)trj->GetPoint(0) ;   
    
    track->fill(0, pId);                                         //pId
    
    parentIndex = trj->GetParentID();
    track->fill(1,parentIndex);                                  //parent
    
    timeStart =  trj->GetGlobalTime()/ns;
    track->fill(2, timeStart);                                   //creation time
    
    mom    = trj->GetInitialMomentum();
    mommag = mom.mag();
    direction = track->getTuple( 3 );                            //direction
    dx = mom.x()/mommag;
    dy = mom.y()/mommag;
    dz = mom.z()/mommag;
    direction->fill(0, dx);
    direction->fill(1, dy);
    direction->fill(2, dz);
    direction->addRow();
    
    mass   = trj->GetParticleDefinition()->GetPDGMass();
    track->fill(4, mass);                                        //mass
    
    pTot = mommag;
    track->fill(5, pTot);                                        //p Total
    
    ETot = sqrt(mom.mag2() + mass*mass);
    track->fill(6, ETot);                                        //Total energy
    
    momentum = track->getTuple( 7 );                             //3-momentum
    px = mom.x();
    py = mom.y();
    pz = mom.z();
    momentum->fill(0, px);
    momentum->fill(1, py);
    momentum->fill(2, pz);
    momentum->addRow();
    
    start  = aa->GetPosition();
    startPos = track->getTuple( 8 );                             //start position
    startPos->fill(0, start.x()/cm);
    startPos->fill(1, start.y()/cm);
    startPos->fill(2, start.z()/cm);
    startPos->addRow(); 
    
    stop   = trj->GetStoppingPoint();
    stopPos = track->getTuple( 9 );                              //stop position
    stopPos->fill(0, stop.x()/cm);
    stopPos->fill(1, stop.y()/cm);
    stopPos->fill(2, stop.z()/cm);
    stopPos->addRow(); 
    
    startVol = EventFindStartingVolume(start);
    track->fill(10, startVol);                                   //Starting Volume
    
    
    stopVolumeName = trj->GetStoppingVolume()->GetName();
    stopVol  = EventFindStoppingVolume(stopVolumeName);
    track->fill(11, stopVol);                                    //Stopping Volume
    
    //add the new track to tuple
    track->addRow(); 

    G4cout << "----> Tk{"<<i<<"}: " 
	   << " pId " << pId
	   << " parent " << parent
	   << " creation time " << timeStart 
	   << " Volumes " << startVol << " " << stopVol << "\n"
	   << " Start Pos (" << start.x()/cm << "," << start.y() << "," << start.z() << ")\n"
	   << " Stop Pos (" << stop.x()/cm << "," << stop.y() << "," << stop.z() << ")\n"
	   << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n"
	   << " m " << mass
	   << " ETot " << ETot
	   << " pTot " << pTot
	   << " px,py,pz " << px << " " << py << " " << pz << "\n"
	   << G4endl;
       
  }//end of loop on trajectories
  

  
  // --------------------
  //  Get Hit Collection
  // --------------------
  //JEC 22/2/07 not yet available
    
//   G4SDManager* SDman = G4SDManager::GetSDMpointer(); //JEC FIXME: use data member
    
//   // Get Hit collection of this event
//   G4HCofThisEvent* HCE         = evt->GetHCofThisEvent();
//   TrappVolHitsCollection* TrappVolHC = 0;

//   AIDA::ITuple* hit = eventTuple->getTuple(6);                         //hit
//   G4int nHits=0;

//   if (HCE) { 
//     G4String name = "TrappingVolume"; //see TrappingVolume.cxx Ctor
//     G4int collectionID = SDman->GetCollectionID(name);
//     TrappVolHC = (TrappVolHitsCollection*)HCE->GetHC(collectionID);
//   }
  
//   if (TrappVolHC) {

//     nHits = TrappVolHC->entries();
//     eventTuple->fill(5, nHits);                                          //nHits
    
//     TrappVolHit* aHit;

//     for (G4int i=0; i<nHits  ;i++) {
//       aHit = (*TrappVolHC)[i];

//       //dump
//       aHit->Print();

//       hit->fill(0,aHit->GetTime());                  //time of hit creation
//       hit->fill(1,aHit->GetEtot());                  //Total energy of the particle creating the hit 
//       hit->fill(2,aHit->GetCoordinates().x());       //position of the hit
//       hit->fill(3,aHit->GetCoordinates().y());
//       hit->fill(4,aHit->GetCoordinates().z());
//       hit->fill(5,aHit->GetTrkId());                 //Id of the particle
//       hit->fill(6,aHit->GetParentId());              //Id of the parent of the particle
//       hit->fill(7,aHit->GetPDGEncoding());           //PDG code of the particle
//       hit->fill(8,aHit->GetEdep());                  //Energy deposit

//       hit->addRow();
//     }//loop on hits
//   }//Hit container

  //Save the Event
  eventTuple->addRow();

}//EndOfEventAction

//--------------------------------------------------------------------------------------------------

G4int ELYSE::EventAction::EventFindStartingVolume(G4ThreeVector vtx) {
  // Get volume of starting point (see GEANT4 FAQ)

  G4Navigator* tmpNavigator = 
    G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();

  G4VPhysicalVolume* tmpVolume = tmpNavigator->LocateGlobalPointAndSetup(vtx);
  G4String       vtxVolumeName = tmpVolume->GetName(); 
   
  return FindVolume(vtxVolumeName);

}//EventFindStartingVolume

//--------------------------------------------------------------------------------------------------

G4int ELYSE::EventAction::EventFindStoppingVolume(G4String stopVolumeName) {
  return FindVolume(stopVolumeName);
}//EventFindStoppingVolume


//--------------------------------------------------------------------------------------------------
G4int ELYSE::EventAction::FindVolume(G4String volumeName) {
  G4int volumeCode = -1;

  G4cout << "(FindVolume): VolumeName: <"
	 << volumeName
	 << ">" 
	 << G4endl;

  if ( volumeName == "expHall" ) {
    volumeCode = 0;
  } else if ( volumeName == "SensiVol" ) {
    volumeCode = 1;
  } else if (volumeName == "TyvekSheet") {
    volumeCode = 2;
  } else if (volumeName == "waterTank") {
    volumeCode = 3;
  }

  return volumeCode;
}//FindVolume