#include <cmath>	// for the function 'exp()'

#include "FCALSteppingAction.hh"
#include "FCALDetectorConstruction.hh"
#include "FCALEventAction.hh"

#include "G4StepStatus.hh"
#include "G4Track.hh"
#include "G4Step.hh"
#include "G4ios.hh"
#include "G4ParticleDefinition.hh"
#include "G4Navigator.hh" 

#include "Randomize.hh"

// parametrization of the quantum efficiency of the FEU-84-3 PMT tube
// as polynomial order 5
#define A  -2545.6
#define B  27.133
#define C -0.11521
#define D 2.45206E-04
#define E -2.60431E-07
#define F 1.09761E-10


FCALSteppingAction::FCALSteppingAction(FCALEventAction* event_Act, G4double zLG, G4int qeff)
  : event_Action(event_Act), zGlass(zLG), QE(qeff) {}

FCALSteppingAction::~FCALSteppingAction() {;}

void FCALSteppingAction::UserSteppingAction(const G4Step* aStep)
{
  //			-------- 1. Extract information of energy deposited in LeadGlass blocks
  G4StepPoint* pre_Step = aStep->GetPreStepPoint();
  G4TouchableHandle pre_TouchHandle = pre_Step->GetTouchableHandle();
  G4VPhysicalVolume* pre_Volume = pre_TouchHandle->GetVolume();
  
  G4StepPoint* post_Step = aStep->GetPostStepPoint();
  G4TouchableHandle post_TouchHandle = post_Step->GetTouchableHandle();
  G4VPhysicalVolume* post_Volume = post_TouchHandle->GetVolume();
  
  // if primary particle == the photon stops record its position
  // => this is where the shower starts
  if ((aStep->GetTrack()->GetTrackID()== 1)&&
      (aStep->GetTrack()->GetTrackStatus()!=fAlive)){
    
    if (aStep->GetTrack()->GetDefinition()->GetParticleName()!="gamma"){
      //G4cout<<"Error non gamma as primary particle found!"<< G4endl;  
      // THIS SHOULD NEVER HAPPEN!
    } else{
      //get position in global and local coordinates
      G4ThreeVector worldPos = post_Step->GetPosition();
      G4ThreeVector localPos = 
	pre_TouchHandle->GetHistory()->GetTopTransform().TransformPoint(worldPos);
      event_Action->recordShowerStart(worldPos.x(),worldPos.y(),worldPos.z(),
				      localPos.x(),localPos.y(),localPos.z());
      event_Action->recordPhotonEnergy((Float_t)aStep->GetTrack()->GetVertexKineticEnergy());
    }
  }
  
  if (pre_Volume != NULL && pre_Volume->GetName() == "LeadGlass")
    {
      G4double energy_Deposition = aStep->GetTotalEnergyDeposit();
      // G4double z_PrePosition = (pre_Step->GetPosition()).z();
      // G4double z_PostPosition = (post_Step->GetPosition()).z();
      // G4double z_AvgPosition = (z_PrePosition + z_PostPosition) / 2.;
      
      G4ThreeVector worldPos1 = pre_Step->GetPosition();
      G4ThreeVector localPos1 = 
	pre_TouchHandle->GetHistory()->GetTopTransform().TransformPoint(worldPos1);
      
      G4ThreeVector worldPos2 = post_Step->GetPosition();
      G4ThreeVector localPos2 = 
	post_TouchHandle->GetHistory()->GetTopTransform().TransformPoint(worldPos2);
      
      G4double z_AvgPos = localPos1.z()-localPos2.z(); // GEANT4 uses millimeter internally
      
      
      G4double energy_Attenuate = energy_Deposition * exp((z_AvgPos-zGlass*cm) / 1600*mm);
      
      
      if (energy_Deposition > 0)
	{
	  G4int copy_Number = pre_TouchHandle->GetCopyNumber(2);
	  event_Action->addEnergyDeposit(energy_Deposition, copy_Number);
	  event_Action->addEnergyAttenuate(energy_Attenuate, copy_Number);
	}
    }
  
  // -------- 2. Extract associated infomation of optical photons hitting the photocathode
  
  
  if (post_Volume != NULL && post_Volume->GetName() == "PMTCathode")
    {
      G4String particle_Name = aStep->GetTrack()->GetDefinition()->GetParticleName();
      if (particle_Name == "opticalphoton")	
	//////////////////////////////////////////////////////////////////////////////////////////////////
	// Actually, the above line might be meaningless. But if you generate, for example, muons, 	//
	// you can come across some occation when the associated detector counts them.			//
	// Of course, when it comes to the comparative number, this must be meaningless			//
	// since huge number of optical photons will be detected. This is just because of rigorousness.	//
	//////////////////////////////////////////////////////////////////////////////////////////////////
	{
	  G4double photon_Energy = post_Step->GetTotalEnergy();
	  G4int copy_Number = post_TouchHandle->GetCopyNumber(2);
	  
	  if (QE){
	    // convert photon energy to wave length
	    G4double wl = 1239.84/(photon_Energy*1000000.);
	    G4double qe = A+B*wl+C*wl*wl+D*wl*wl*wl+E*wl*wl*wl*wl+F*wl*wl*wl*wl*wl;
	    qe /=100.;
	    //G4cout<<"QE="<<qe<<"     Rand="<<R<<"   E="
	    //<<photon_Energy*1000000.<<"   wl="<<wl<<G4endl;

	    if (qe<0) qe=0;
	    if ((wl<300.) || (wl>620.)) qe=0;

	    G4double R = G4UniformRand();
	    if (R<=qe){
	      event_Action->countHits(copy_Number);			// count Hits
	      event_Action->addEnergy(photon_Energy, copy_Number);	// add energy
	    }

	  } else {
	      event_Action->countHits(copy_Number);			// count Hits
	      event_Action->addEnergy(photon_Energy, copy_Number);	// add energy
	  }


	}
    }
}