//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION 
//
//  ID    : LEMuSRPrimaryGeneratorAction.cc , v 1.3
//  AUTHOR: Taofiq PARAISO
//  DATE  : 2004-09-16 09:12
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                          PRIMARY GENERATOR ACTION
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//


 

#include "LEMuSRPrimaryGeneratorAction.hh"
#include "LEMuSRPgaMessenger.hh"
#include "yields.h"

//  G4 LIBRARIES

#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4UnitsTable.hh"
#include "Randomize.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh" 
#include "LEMuSRParticleGun.hh"
#include "globals.hh"
#include "Randomize.hh"
#include "G4Gamma.hh"
#include "G4ThreeVector.hh"
#include "globals.hh"
#include "G4RunManager.hh"
#include "time.h"
#include "G4ios.hh"
#include <iomanip>

#include "LEMuSRMuoniumParticle.hh"

/*!
 * In the constructor the PGA messenger is initialized 
 * the particle gun is created and the default
 * values are given to the variables for the particle definition
 * or the scanning parameters
 */
LEMuSRPrimaryGeneratorAction::LEMuSRPrimaryGeneratorAction()
{
 
   
  // Instance pointer
  thePGA=this;

  // Messenger
  messenger = new LEMuSRPgaMessenger(this);

  // get the random number engine
  theEngine = HepRandom::getTheEngine();

  m_nbxsteps = 200;
  m_nbysteps = 200;
  m_counterx = 0;
  m_countery = 0;
  m_xrange   = 2.;
  m_yrange   = 2.;
  
  lemuParticleGun = new LEMuSRParticleGun();
  lemuParticleGun->SetNumberOfParticles(1);
  
  X=0;  Y=0;  
  Z=-114.;
  momX=0.; momY=0.; momZ=1;
   
  scan=0;
  
  rnd=false;
  
  energy=5*keV;
  
  gauss=0;  mean=0;  stddev=0;

  pname="mu+";
 
  ke_offset=3.73*keV;
#if defined LEMU_TEST_ASYM  
  ke_offset=0;
#endif
  charge=+1;

  //==================================TEST
  /*  i=j=0;
  Yprint.open("test5.txt");
  if (!Yprint.is_open()) exit(8);
  while(i<1000){table[i]=0; i++;}
  */

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OO


/*!
 * The particle gun and the messenger pointers are deleted inthe destructor.
 */

LEMuSRPrimaryGeneratorAction::~LEMuSRPrimaryGeneratorAction()
{
  delete lemuParticleGun;  
  delete messenger;


  //============================0000TEST
  /*  i=0;
  do 
    {
      //      G4cout<<i*0.01<<"  " <<table[i] <<std::endl;
      Yprint<<i*0.01<<"  " <<table[i] <<std::endl;
      i++;
    }while(i<500);
    Yprint.close();*/
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OO


/*!
 * This method defines all the parameters for the particle gun and
 * calls the GeneratePrimaryVertex from the LEMuSRParticleGun class.
 * An important submethod is the GetScanningMode. which is called to 
 * get different distribution of the initial gun position.
 *
 * The decay time is preliminary assigned for muons and muonium particles.
*/

void LEMuSRPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{


  //----------------------GETTING RANDOM VALUES >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

  // the random engine
  G4RandGauss* iRndGauss = new G4RandGauss(theEngine, 20.0, 10.0);

  // the random energy
  if(rnd)
  {
      rndenergy = iRndGauss->shoot(20.0, 0.50);
      energy= rndenergy*keV;//default unit is MeV  
  }

  if(gauss==1)
  {
      rndenergy = iRndGauss->shoot(mean, stddev);
      energy= rndenergy*keV;//default unit is MeV  
  }

#ifdef DEBUG_LEMU_PGA
  G4cout<<"Energy is "<<G4BestUnit(energy,"Energy")<<G4endl;
#endif



  //----------------------GETTING POSITIVE MUONS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  G4ParticleDefinition* particle;

  particle = particleTable->FindParticle(pname);

  //  std::cout<<"Particle "<< particle->GetParticleName();
  
  if(particle->GetParticleSubType()=="mu")
    {
      // the random time of the decay for muons
  
      G4double tau=0.000002197; // the muon mean lifetime
      G4double rnddecaytime = - tau*log(1-G4UniformRand());
      decaytime = rnddecaytime*s;
      lemuParticleGun->SetDecayTime(decaytime);  
      //     G4cout<<"Decay Time is "<<decaytime/ns <<" nanoseconds."<<G4endl;
    }

  
  lemuParticleGun->SetParticleDefinition(particle);

  lemuParticleGun->SetParticleEnergy(energy+ke_offset);

  
  GetScanningMode(scan);

  //      G4cout<<"Zposition "<<Z/mm<<" [mm]"<<G4endl;
  //-------------------MOMENTUM DIRECTION AND POLARIZATION>>>>>>>>>>>>>>>>>>>
  
  lemuParticleGun->SetParticleMomentumDirection(G4ThreeVector(momX,momY,momZ));
  lemuParticleGun->SetParticlePolarization(G4ThreeVector(1.,0.,0.));
#if defined LEMU_TEST_ASYM  
  lemuParticleGun->SetParticlePosition(G4ThreeVector(0.,0.,0.));
  lemuParticleGun->SetParticlePolarization(G4ThreeVector(-1.,0.,0.)); 
  lemuParticleGun->SetParticleEnergy(0);
#endif
  
  
  lemuParticleGun->GeneratePrimaryVertex(anEvent);
  
}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OO
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OO



/*!
 * The different scanning modes are distributions of the initial position of the gun.
 * Few of them were implemented and one can modify them as he likes.
 * The user can switch the scanning mode before each run via the user interface
 * terminal thanks to the commands implemented in LEMuSRPGAMessenger.
 * 
 * So far, one can choose between the following distributions:
 * - Single position of the gun
 * - Square sweeping: the gun position will be uniformely distributed in a square. The dimensions and steps along x and y directions are defined by the user.
 * - Circular scan: the radius is scaned linearly and for all radius a circle is described according to the specified steps number. This distribution <I> is not</I> homogeneous.
 * - Gaussian scan: x and y are randomly determined according to a gaussian distribution. The mean value and the standard deviation are user parameters.
 *.
 * The following picture shows different beam profiles generated by modifying the parameters of the gaussian circular scan. The maximal radius is limited to 2.5cm around the mean value.
 */
void LEMuSRPrimaryGeneratorAction::GetScanningMode(G4int scan)
{

  //----------------------SINGLE SHOOTING POSITION >>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  if(scan==0)
    {
      lemuParticleGun->SetParticlePosition(G4ThreeVector(X*cm,Y*cm,Z*cm));
    }




  //-------------------------SQUARE SCAN >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  if(scan==1)
    {
      //scan: position of the gun
      if(m_countery<=m_nbysteps)
	{
	  if(m_counterx<=m_nbxsteps)
	    {
	      //   total scan
	      X = -m_xrange+m_counterx*(2*m_xrange)/m_nbxsteps;
	      Y = -m_yrange+m_countery*(2*m_yrange)/m_nbysteps;
	      m_counterx++;
	    }
	  
	  if(m_counterx>m_nbxsteps)
	    {
	      m_countery++;
	      m_counterx=1;
	    }
	}
      
      else if(m_countery>m_nbysteps)
	{
	  m_countery=0;
	  m_counterx=0; 
	}
      // SET POSITION
      lemuParticleGun->SetParticlePosition(G4ThreeVector(X*cm,Y*cm,-114*cm));
    }
  
  
  
  
  //----------------------CIRCULAR SCAN>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  if(scan==2)
    {
      if(m_countery<=m_nbysteps)
	{
	  if(m_counterx<=m_nbxsteps)
	    {
	      radius= m_xrange*(m_counterx/m_nbxsteps)*cm;
	      angle=2*M_PI*m_countery/(m_nbysteps)*rad;
	      X = radius*cos(angle);
	      Y = radius*sin(angle);
	      m_counterx++;
	      
	    }
	  
	  if(m_counterx>m_nbxsteps)
	    {
	      m_countery++;
	      m_counterx=1;
	    } 
	}
      
      else if (m_countery>m_nbysteps)   
	{
	  //	  X=0;
	  //	  Y=0;
	  m_countery=0;
	  m_counterx=0; 
	}
      // SET POSITION
      lemuParticleGun->SetParticlePosition(G4ThreeVector(X,Y,-114*cm));
    }  
  
 
 //---------------------- GAUSSIAN SCAN  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


  if(scan==3)
    {
      G4RandGauss* iRndGauss = new G4RandGauss(theEngine);
      
      
      //      G4cout <<"radius "<< radius/cm<<G4endl;
   do
     {   
      X = iRndGauss->shoot(sc_mean, sc_stddev)*cm;
      Y = iRndGauss->shoot(sc_mean, sc_stddev)*cm;
  
      radius= sqrt((X-sc_mean)*(X-sc_mean)+(Y-sc_mean)*(Y-sc_mean));
     }while(radius>2.5*cm);
      // SET POSITION
      lemuParticleGun->SetParticlePosition(G4ThreeVector(X,Y,-114*cm));
    }  
  
  
  
}





//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OO

LEMuSRPrimaryGeneratorAction*  LEMuSRPrimaryGeneratorAction::thePGA=0;
LEMuSRPrimaryGeneratorAction*  LEMuSRPrimaryGeneratorAction::GetPGA()
{
  return thePGA;
}