//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION 
// 
//  ID    :FocalLengthTest.cc , v 1.2
//  AUTHOR: Taofiq PARAISO
//  DATE  :  2006-01-19 15:17
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                               FOCALLENGTHTEST   
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//

#include "FocalLengthTest.hh"

#include "G4SteppingManager.hh"
#include "G4Transform3D.hh"
#include "G4DynamicParticle.hh"
#include "G4UnitsTable.hh"

#include "LEMuSRPrimaryGeneratorAction.hh"


FocalLengthTest::FocalLengthTest()
{

  BookRoot();
  muon.event=0;
  loop=0;
  pointer=this  ;
}


FocalLengthTest::~FocalLengthTest()
{
  WriteRoot();
}

FocalLengthTest*  FocalLengthTest::pointer=0;
FocalLengthTest*  FocalLengthTest::GetInstance()
{
  return pointer;
}


void FocalLengthTest::UserSteppingAction(const G4Step* aStep)
  
{  
  if( aStep->GetPreStepPoint()&& aStep->GetPreStepPoint()->GetPhysicalVolume() )
    {

      LoopKiller(aStep);
      SetParticleVolumeNames(aStep);
      
      if(CheckCondition(aStep))
	{	  
	  
	  // Get datas
	  SetPositionMomentum(aStep);
	  SetTimeEnergy(aStep);
	  
	  Update();
	  FillRoot();
#if defined G4UI_USE_ROOT 
	  myFile->Write();
#endif

	}
    }
  
  
}

void FocalLengthTest::LoopKiller(const G4Step *aStep)
{
  // loop killa
  if(aStep->GetStepLength()<0.001*mm)
    {
      loop++;
      if(loop>20)
	{
	  aStep->GetTrack()->SetTrackStatus(fStopAndKill);
	  loop=0;
	}
    }

      // kill useless particles
  if(aStep->GetTrack()->GetDefinition()->GetParticleName()=="e-"
     ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="gamma")
    {
      if(aStep->GetTrack()->GetKineticEnergy()<10.*eV)
      {
	aStep->GetTrack()->SetTrackStatus(fStopAndKill);
      }
    }
  
}




G4bool FocalLengthTest::CheckCondition(const G4Step* )
{
  G4bool condition=false;

  //  if(p_name == "geantino" && v_name=="lv_fchk") // change also in pga
  if(p_name == "mu+" && v_name=="lv_fchk")
       {	  
	  if(pv_name=="pv_f1"||pv_name=="pv_f44")// THE FIRST DUMMY PLANE 
	    {
	      condition=true;
	    }
       }

  return condition;

}  



void FocalLengthTest::SetParticleVolumeNames(const G4Step* aStep)
{
  // NAMES
  p_name = aStep->GetTrack()->GetDefinition()->GetParticleName(); // particle name 
  v_name = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName(); //lv_name
  pv_name = aStep->GetTrack()->GetVolume()->GetName(); //lv_name

}


void FocalLengthTest::SetPositionMomentum(const G4Step* aStep)
{

  // POSITION, MOMENTUM
  position = aStep->GetPreStepPoint()->GetPosition(); // position
  momentum = aStep->GetPreStepPoint()->GetMomentum(); // momentum
  momentum_direction = aStep->GetPreStepPoint()->GetMomentumDirection(); // momentum

}


void FocalLengthTest::SetTimeEnergy(const G4Step* aStep)
{
 
  // ENERGY
  kenergy= aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy(); // position
 
  tenergy= aStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy(); // position
 
  // TIME
  localtime =  aStep->GetPreStepPoint()->GetLocalTime(); // time since track creation
  globaltime =  aStep->GetPreStepPoint()->GetGlobalTime();// time since event creation
  proptime =  aStep->GetPreStepPoint()->GetProperTime(); // proper time of the particle
  time = proptime;

}


void FocalLengthTest::SetSpinDirection(const G4Step* aStep)
{
  polarization = aStep->GetTrack()->GetDynamicParticle()->GetPolarization();

}


void FocalLengthTest::Update()
{
  muon.index++;

 
  if(pv_name=="pv_f1")// THE FIRST DUMMY PLANE 
    {
      muon.index=0;
      muon.event++;
      init_kenergy=LEMuSRPrimaryGeneratorAction::GetPGA()->energy;
    }

  if(pv_name=="pv_f44")
    { 
      G4double FL, h, l, d;
      h=sqrt(position.x()*position.x()+position.y()*position.y());
      d=sqrt(momentum_direction.x()*momentum_direction.x()+momentum_direction.y()*momentum_direction.y());
      l=momentum_direction.z();
      
      FL= h*l/d+22.*cm;
      // 22 cm from the last plane pv_f44 to the center of the lense 
      //G4cout <<"Focal length for this particle : " << FL/meter <<"meter" <<G4endl;
      muon.focal = FL/centimeter;
      muon.kenergy = LEMuSRPrimaryGeneratorAction::GetPGA()->energy/keV;
      //      G4cout <<"Kinectic energy  for this particle : " << init_kenergy/keV <<"keV" <<G4endl;
      if(init_kenergy/keV<10)
	{
      muon.ratio=(-10*keV-init_kenergy)/(-init_kenergy);
	}
      else{
      muon.ratio=(10*keV-init_kenergy)/(-init_kenergy);
      }
      //      G4cout <<"Ratio for this particle : " << (15*keV-init_kenergy)/init_kenergy <<"." <<G4endl;
    }


  muon.positionx = position.x()/mm;
  muon.positiony = position.y()/mm;
  muon.positionz = position.z()/mm;
  
  muon.momdirx= momentum_direction.x();
  muon.momdiry= momentum_direction.y();
  muon.momdirz= momentum_direction.z();


  
}


//                          PRINT VALUES

void FocalLengthTest::PrintField(const G4Step* )
{

}


void FocalLengthTest::PrintDatas(const G4Step* )
{
 
}    




//----------------------------------------------------------------------------------------//
//                                            ROOT

void FocalLengthTest::BookRoot()
{
  myFile = new TFile("focal_length.root", "RECREATE");

  tree = new TTree ("tree","Muons parameters");

  tree->Branch("muon",&muon.kenergy,"kenergy/F:focal/F:ratio/F:positionx/F:positiony/F:positionz/F:momdirx/F:momdiry/F:momdirz");

}


void FocalLengthTest::FillRoot()
{
  tree->Fill();
  //  tree->Scan();

}

void FocalLengthTest::WriteRoot()
{
  myFile->Write();
  myFile->Close();

}