//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file B1SteppingAction.cc
/// \brief Implementation of the B1SteppingAction class

#include "B1SteppingAction.hh"
#include "B1EventAction.hh"
#include "B1DetectorConstruction.hh"
#include "G4DynamicParticle.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "g4analysis.hh"
#include "G4LogicalVolume.hh"
#include "G4Track.hh"
#include "G4ParticleMomentum.hh"
#include "SteppingActionMessenger.hh"

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

B1SteppingAction::B1SteppingAction(B1EventAction* eventAction)
: G4UserSteppingAction(),
  fEventAction(eventAction),
  fScoringVolume(0),
  
  //////////---------------START SETTING FOR FSCORINGVOLUME2------------//////
  
  fScoringVolume2(0)
  //////////---------------END SETTING FOR FSCORINGVOLUME2------------//////
{}

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

B1SteppingAction::~B1SteppingAction()
{}

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

void B1SteppingAction::UserSteppingAction(const G4Step* step)
{
	
	 //-----------START DECLARE ANALYSIS MANAGER-------------------------------//
   auto analysisManager = G4AnalysisManager::Instance();
   //-----------END DECLARE ANALYSIS MANAGER-------------------------------//
  if (!fScoringVolume) { 
    const B1DetectorConstruction* detectorConstruction
      = static_cast<const B1DetectorConstruction*>
        (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fScoringVolume = detectorConstruction->GetScoringVolume();   
    fScoringVolume2 = detectorConstruction->GetScoringVolume2();
}

  // get volume of the current step
  G4LogicalVolume* volume 
    = step->GetPreStepPoint()->GetTouchableHandle()
      ->GetVolume()->GetLogicalVolume();
      
      G4VPhysicalVolume* ThisVol = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
	G4VPhysicalVolume* NextVol = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume();
  // check if we are in scoring volume

  // collect energy deposited in this step
//  G4double edepStep = step->GetTotalEnergyDeposit();
//  fEventAction->AddEdep(edepStep);
  
  //////////---------------START FSCORINGVOLUME2------------//////
//   if (!fScoringVolume2) {
//    const B1DetectorConstruction* detectorConstruction2
//      = static_cast<const B1DetectorConstruction*>
//        (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
//    fScoringVolume2 = detectorConstruction2->GetScoringVolume2();
//  }

  // get volume of the current step
//  G4LogicalVolume* volume2
//    = step->GetPreStepPoint()->GetTouchableHandle()
//      ->GetVolume()->GetLogicalVolume();
      
  // check if we are in scoring volume

  
  //--------------START CHECK VOLUME WE ARE IN--------------///////
  //if (volume2 != fScoringVolume2) return;
//  if (volume != fScoringVolume2 && volume != fScoringVolume) return;
  //--------------START CHECK VOLUME WE ARE IN--------------///////

  // collect energy deposited in this step
  
	if (volume == fScoringVolume) {
  G4double edepStep = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep(edepStep);
  if(RootFlag){
  int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 analysisManager->FillNtupleDColumn(0,0, pdg);
			 double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 analysisManager->FillNtupleDColumn(0,1, kinEnergy);
			 double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 analysisManager->FillNtupleDColumn(0,2, MomDirx);
			 double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 analysisManager->FillNtupleDColumn(0,3, MomDiry);
			 double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 analysisManager->FillNtupleDColumn(0,4, MomDirz); 
			 double PosDirx = step->GetTrack()->GetPosition().x();
			 analysisManager->FillNtupleDColumn(0,5, PosDirx); 
			 double PosDiry = step->GetTrack()->GetPosition().y();
			 analysisManager->FillNtupleDColumn(0,6, PosDiry); 
			 double PosDirz = step->GetTrack()->GetPosition().z();
			 analysisManager->FillNtupleDColumn(0,7, PosDirz); 
			 double EndStepPosx = step->GetPostStepPoint()->GetPosition().x();
			 analysisManager->FillNtupleDColumn(0,8, EndStepPosx); 
			 double EndStepPosy = step->GetPostStepPoint()->GetPosition().y();
			 analysisManager->FillNtupleDColumn(0,9, EndStepPosy); 
			 double EndStepPosz = step->GetPostStepPoint()->GetPosition().z();
			 analysisManager->FillNtupleDColumn(0,10, EndStepPosz); 
			 analysisManager->FillNtupleDColumn(0,11, edepStep);
		}
			  analysisManager->AddNtupleRow();  
	}
	
	if (volume == fScoringVolume2) {
  G4double edepStep2 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep2(edepStep2);
  /*	int pdg2 = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			analysisManager->FillNtupleDColumn(1,0, pdg2);
			 double kinEnergy2 = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 analysisManager->FillNtupleDColumn(1,1, kinEnergy2);
			 double MomDirx2 = step->GetTrack()->GetMomentumDirection().x();
			 analysisManager->FillNtupleDColumn(1,2, MomDirx2);
			 double MomDiry2 = step->GetTrack()->GetMomentumDirection().y();
			 analysisManager->FillNtupleDColumn(1,3, MomDiry2);
			 double MomDirz2 = step->GetTrack()->GetMomentumDirection().z();
			 analysisManager->FillNtupleDColumn(1,4, MomDirz2);
			 double PosDirx2 = step->GetTrack()->GetPosition().x();
			 analysisManager->FillNtupleDColumn(1,5, PosDirx2); 
			 double PosDiry2 = step->GetTrack()->GetPosition().y();
			 analysisManager->FillNtupleDColumn(1,6, PosDiry2); 
			 double PosDirz2 = step->GetTrack()->GetPosition().z();
			 analysisManager->FillNtupleDColumn(1,7, PosDirz2); 
			 double EndStepPosx2 = step->GetPostStepPoint()->GetPosition().x();
			 analysisManager->FillNtupleDColumn(1,8, EndStepPosx2); 
			 double EndStepPosy2 = step->GetPostStepPoint()->GetPosition().y();
			 analysisManager->FillNtupleDColumn(1,9, EndStepPosy2); 
			 double EndStepPosz2 = step->GetPostStepPoint()->GetPosition().z();
			 analysisManager->FillNtupleDColumn(1,10, EndStepPosz2); 
			 analysisManager->FillNtupleDColumn(1,11, edepStep2);*/
			 analysisManager->AddNtupleRow(1);  
			 
	

if(NextVol && ThisVol->GetName()=="Envelope2" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && (step->GetTrack()->GetDynamicParticle()->GetPDGcode()==13 || step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-13 )) { 
//if(NextVol && ThisVol->GetName()=="Envelope2" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()==00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11) { 
			 		    int secpdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(3,0, secpdg);
			 			double seckinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(3,1, seckinEnergy);
			 			double SecMomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(3,2, SecMomDirx);
			 			double SecMomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(3,3, SecMomDiry);
			 			double SecMomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(3,4, SecMomDirz);
			 			double SecVertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(3,5, SecVertx);
						double SecVerty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(3,6, SecVerty);
						double SecVertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(3,7, SecVertz);
						double SecMom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(3,8, SecMom);
						double PosDirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(3,9, PosDirxsp); 
			 			double PosDirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(3,10, PosDirysp); 
			 			double PosDirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(3,11, PosDirzsp); 
			 			double SecMomx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(3,12, SecMomx);
						double SecMomy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(3,13, SecMomy);
						double SecMomz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(3,14, SecMomz);
						double Secxprime=acos((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetTotalMomentum()));
						analysisManager->FillNtupleDColumn(3,15, Secxprime);
						double Secyprime=acos((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetTotalMomentum()));
						analysisManager->FillNtupleDColumn(3,16, Secyprime);
						//double InvMass = step->GetInvariantMass();
						//analysisManager->FillNtupleDColumn(3,17, InvMass);
			   			analysisManager->AddNtupleRow(3);  
		}
	
}
		

  if (volume != fScoringVolume || volume != fScoringVolume2) return;
  
}

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