#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

G4int old=0;
G4int cnt=0;

G4ThreeVector photon_dir;
G4ThreeVector photon_pos;
G4ThreeVector photon_pos1;
G4int flag1 = 0;
G4int flag2 = 0;
G4int flag3 = 0;
G4int flag4 = 0;
G4int tid = 0;

FCALSteppingAction::FCALSteppingAction(FCALEventAction* event_Act, G4double zLG, 
				       G4int qeff, TH1F **hist, TH2F **hist2d)
  : event_Action(event_Act), zGlass(zLG), QE(qeff), Spec(hist), Spec2d(hist2d){}

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){

	//	G4cout<< "Energy Loss: "<<energy_Deposition<<"     "<<energy_Attenuate<<G4endl;

	  G4int copy_Number = pre_TouchHandle->GetCopyNumber(2);
	  event_Action->addEnergyDeposit(energy_Deposition, copy_Number);
	  event_Action->addEnergyAttenuate(energy_Attenuate, copy_Number);
      }
  }
  
  // get impact angle of photon on lead glass-PMT interace surface
  if (aStep->GetTrack()->GetDefinition()->GetParticleName()=="opticalphoton"){
    if (pre_Volume!=NULL && post_Volume!=NULL){

      photon_pos1 = post_Step->GetPosition();
      photon_pos = pre_TouchHandle->GetHistory()->GetTopTransform().
	TransformPoint(pre_Step->GetPosition());
      photon_dir = aStep->GetTrack()->GetMomentumDirection();

      if (pre_Volume->GetName() ==  "LeadGlass" &&
	  post_Volume->GetName()!= "LeadGlass" &&
	  photon_pos1.z() >= 455.){
	
	G4double x1 = photon_pos.x() - photon_dir.x()/
	  photon_dir.z()*(photon_pos.z()-225.);
	G4double y1 = photon_pos.y() - photon_dir.y()/
	  photon_dir.z()*(photon_pos.z()-225.);
	
	Spec[1]->Fill(photon_dir.cosTheta()); 
	Spec2d[0]->Fill(x1*0.1+2.,y1*0.1+2.); 

	G4ThreeVector p = photon_dir;
        G4int n = pre_TouchHandle->GetCopyNumber(2);
	if (n==36)
	  Spec[2]->Fill(p.cosTheta());  // central block in 8x8 array 
	if (n==35)
	  Spec[3]->Fill(p.cosTheta());  // block left
	if (n==37) 
	  Spec[4]->Fill(p.cosTheta());  // block right
	if (n==44)
	  Spec[5]->Fill(p.cosTheta());  // block above
	if (n==28)
	  Spec[6]->Fill(p.cosTheta());  // block below
	if (n==43)
	  Spec[7]->Fill(p.cosTheta());  // block above left
	if (n==45)
	  Spec[8]->Fill(p.cosTheta());  // block above right
	if (n==27)
	  Spec[9]->Fill(p.cosTheta());  // block below left
	if (n==29)
	  Spec[10]->Fill(p.cosTheta()); // block below right

	// factor 0.1 to convert from mm to cm

      }
      
      if (tid != aStep->GetTrack()->GetTrackID()){
	tid = aStep->GetTrack()->GetTrackID();
	flag1 = 0;
	flag2 = 0;
	flag3 = 0;
      }

      if (pre_Volume->GetName() ==  "LeadGlass" && 
	  post_Volume->GetName()== "InterfaceCookie"){
	photon_dir = aStep->GetTrack()->GetMomentumDirection();
	flag1 = 1;
	tid = aStep->GetTrack()->GetTrackID();
      }

      if (flag1==1 && pre_Volume->GetName() ==  "InterfaceCookie"){
	Spec[11]->Fill(photon_dir.cosTheta());
	Spec2d[1]->Fill(photon_pos.x()*0.1+2.,photon_pos.y()*0.1+2.); 
	flag1 = 0;
      }

      if (pre_Volume->GetName() ==  "InterfaceCookie" && 
	  post_Volume->GetName()== "PlexiGlass"){
	flag2 = 1;
      }

      if (flag2==1 && pre_Volume->GetName() == "PlexiGlass" ){
	flag2 = 0;
	Spec[12]->Fill(photon_dir.cosTheta());
      }

      if (pre_Volume->GetName() ==  "PlexiGlass" && 
	  post_Volume->GetName()== "PMTGlass"){
	flag3=1;
      }
      if (pre_Volume->GetName() ==  "PMTGlass" && 
	  post_Volume->GetName()== "PMTCathode"){
	
	G4double x1 = photon_pos.x() - photon_dir.x()/
	  photon_dir.z()*(photon_pos.z()-1.);
	G4double y1 = photon_pos.y() - photon_dir.y()/
	  photon_dir.z()*(photon_pos.z()-1.);
	
	Spec[14]->Fill(photon_dir.cosTheta());
	Spec2d[3]->Fill(x1*0.1+2.,y1*0.1+2.); 
      }
      if (flag3==1 && pre_Volume->GetName() == "PMTGlass"){
	flag3=0;
	Spec[13]->Fill(photon_dir.cosTheta());
	Spec2d[2]->Fill(photon_pos.x()*0.1+2.,photon_pos.y()*0.1+2.); 
	// GEANT4 uses mm internally hence the factor 0.1
      }
    }
  }
  // -------- 2. Extract associated infomation of optical photons hitting the photocathode
  
  //if (aStep->GetTrack()->GetDefinition()->GetParticleName()=="opticalphoton"){
  //  if (post_Volume!=NULL){
  //    if (post_Volume->GetName()!=pre_Volume->GetName()){
  //if(pre_Volume->GetName()=="LeadGlass"){
  //  G4ThreeVector w3 =pre_Step->GetPosition(); 
  //  G4ThreeVector v3 = pre_TouchHandle->GetHistory()->GetTopTransform().
  //    TransformPoint(w3);
  //  G4cout <<post_Volume->GetName()<<"/"<<pre_Volume->GetName()
  //	 <<"   track:"<<aStep->GetTrack()->GetTrackID()
  //	 <<"  x="<<v3.x()<<", y="<<v3.y()<<", z="<<v3.z()
  //	 <<",  zworld="<<w3.z()<<", zworld_poststep="<<photon_pos1.z()<<G4endl;
  //}
  //  } 
  //}
  //}
  
  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);
	  
	  // convert photon energy to wave length
	  G4double wl = 1239.84/(photon_Energy*1000000.);
	  if (QE){
	    
	    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>650) qe=0;

	    G4double R = G4UniformRand();
	    if (R<=qe){
	      Spec[0]->Fill(wl); 
	      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
	  }


	}
    }
}