283 lines
8.2 KiB
C++
283 lines
8.2 KiB
C++
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
// LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION
|
|
//
|
|
// ID : LEMuSRAtRestSpinRotation.cc , v 1.3
|
|
// AUTHOR: Taofiq PARAISO
|
|
// DATE : 2006-01-19 16:15
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
//
|
|
// & &&&&&&&&&& &&&&&&& &&&&&&&&
|
|
// & & && && & &&
|
|
// & & & & & & &&
|
|
// & &&&&&&& & & &&&&&& &&&&&&&&
|
|
// & & & && & & &&
|
|
// & & && & & && && & &
|
|
// &&&&&&&&&& &&&&&&&&&& & &&&&& && &&&&&&& & &&
|
|
// &
|
|
// &
|
|
// &
|
|
// &
|
|
// AT REST SPIN ROTATION
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
|
|
#include "LEMuSRAtRestSpinRotation.hh"
|
|
|
|
#include "G4ThreeVector.hh"
|
|
#include "G4MagneticField.hh"
|
|
#include "G4Track.hh"
|
|
#include "G4FieldManager.hh"
|
|
#include <iomanip>
|
|
#include "globals.hh"
|
|
#include "G4ParticleChange.hh"
|
|
#include "G4ios.hh"
|
|
#include "G4Transform3D.hh"
|
|
#include "G4UnitsTable.hh"
|
|
|
|
#include "G4ParticleTable.hh"
|
|
#include "G4ParticleDefinition.hh"
|
|
|
|
#include "G4ProcessVector.hh"
|
|
#include "G4ProcessManager.hh"
|
|
|
|
|
|
LEMuSRAtRestSpinRotation::LEMuSRAtRestSpinRotation(const G4String& processName)
|
|
: G4VRestProcess (processName)
|
|
{
|
|
G4cout << GetProcessName() << " is created "<< G4endl;
|
|
pointer = this;
|
|
}
|
|
|
|
|
|
LEMuSRAtRestSpinRotation::~LEMuSRAtRestSpinRotation()
|
|
{;}
|
|
|
|
|
|
LEMuSRAtRestSpinRotation* LEMuSRAtRestSpinRotation::pointer=0;
|
|
LEMuSRAtRestSpinRotation* LEMuSRAtRestSpinRotation::GetInstance()
|
|
{
|
|
return pointer;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$ AT REST DO IT $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
/*!
|
|
* This method contains the list of modifications to do when the particle come to rest state. It is the implementation of a virtual method of G4VProcess. Therefore, it's name should not be changed.
|
|
*
|
|
* The procedure is the following:
|
|
* -# Initialize the particle change object according to the particle's properties taken from the track object.
|
|
* -# Propose a status for the track, a position, momentum, energy, and time change.
|
|
* -# Check is the volume where the particle is at rest contains a magnetic field
|
|
* -# If there is a magnetic field, get the B vector and call the RotateSpin method to rotate the spin.
|
|
* -# Propose a spin change
|
|
* -# Return the particle change object.
|
|
* .
|
|
* Then the AtRestUpdateStep method is called cf LEMuSRParticleChangeForSR.
|
|
*/
|
|
G4VParticleChange* LEMuSRAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack, const G4Step& aStep)
|
|
{
|
|
|
|
theParticleChange.Initialize(theTrack);
|
|
theParticleChange.ProposeTrackStatus(fStopButAlive);
|
|
theParticleChange.ProposePosition(theTrack.GetPosition());
|
|
theParticleChange.ProposeMomentumDirection(theTrack.GetMomentumDirection());
|
|
theParticleChange.ProposeEnergy(theTrack.GetKineticEnergy());
|
|
|
|
// Get Time
|
|
itime = theTrack.GetProperTime();
|
|
G4double gtime = theTrack.GetGlobalTime();
|
|
ftime = theTrack.GetDynamicParticle()->GetPreAssignedDecayProperTime();
|
|
|
|
|
|
|
|
deltatime = ftime - itime;
|
|
theParticleChange.ProposeGlobalTime(deltatime + itime -gtime);
|
|
theParticleChange.ProposeProperTime(deltatime);
|
|
|
|
polar = aStep.GetTrack()->GetPolarization();
|
|
|
|
|
|
if(aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetFieldManager())
|
|
{
|
|
|
|
|
|
G4FieldManager *fMgr = theTrack.GetVolume()->GetLogicalVolume()->GetFieldManager();
|
|
|
|
|
|
// if(!fMgr->DoesFieldChangeEnergy())//then we have a magnetic field
|
|
if(fMgr->FieldHasMagComponent())//then we have a magnetic field
|
|
{
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout<<"AT REST::: MAGNETIC FIELD HERE" <<G4endl;
|
|
#endif
|
|
|
|
// Get B Field vector
|
|
point[0]=theTrack.GetPosition().x();
|
|
point[1]=theTrack.GetPosition().y();
|
|
point[2]=theTrack.GetPosition().z();
|
|
|
|
mfield = fMgr->GetDetectorField();
|
|
|
|
mfield->GetFieldValue(point,B);
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout <<"AT REST::: MAGNETIC FIELD B="<< B[0] <<" " << B[1] <<" " << B[2] <<G4endl;
|
|
|
|
#endif
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout <<"AT REST::: TIME= proper: "<< itime <<" ; global: " << gtime <<" decay: " << ftime <<G4endl;
|
|
#endif
|
|
|
|
|
|
G4ThreeVector magField(B[0],B[1],B[2]);
|
|
if(magField!=G4ThreeVector(0.0,0.0,0.0))
|
|
{
|
|
RotateSpin(aStep,magField,deltatime);
|
|
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout<<"AT REST::: spin rotated";
|
|
#endif
|
|
}
|
|
|
|
}
|
|
|
|
theParticleChange.ProposePolarization(polar);
|
|
|
|
|
|
|
|
#ifdef G4SRVERBOSE
|
|
theParticleChange.DumpInfo();
|
|
#endif
|
|
|
|
|
|
return &theParticleChange;
|
|
|
|
// then the AtRestUpdateStep method is called cf LEMuSRParticleChangeForSpinRotation (SR)
|
|
|
|
}
|
|
|
|
else
|
|
{
|
|
// G4cout <<"No MAGNETIC FIELD do not call AtRestPrecession :: RotateSpin!!! ";
|
|
return &theParticleChange;
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//$$$$$$$$$$$$$$$$$$$$$$$$$$ ROTATE SPIN $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
/*!
|
|
* This method rotates the spin of the particle according to the field value.It also takes the gyromagnetic parameters from the particle definition. It will make the spin precess until the particle decays. The gyromagnetic factor is taken from the particle definition.
|
|
*/
|
|
void LEMuSRAtRestSpinRotation::RotateSpin( const G4Step& aStep, G4ThreeVector B, G4double deltatime )
|
|
{
|
|
|
|
|
|
G4Transform3D Spin_rotation;
|
|
|
|
G4double Bnorm = sqrt(sqr(B[0]) + sqr(B[1]) +sqr(B[2]) );
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout<< "AT REST::: PARAMETERS\n"
|
|
<< "Magnetic Field Norm : " << G4BestUnit(Bnorm,"Magnetic flux density") <<"\n";
|
|
#endif
|
|
|
|
//! Gyromagnetic Ratio
|
|
G4double gamma = GetGyromagneticRatio();
|
|
//! Compute precession frequency
|
|
G4double omega= - (gamma) * Bnorm;
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout<< "AT REST::: PARAMETERS\n"
|
|
<< "Frequency: " << G4BestUnit(fqz*gauss/(2*M_PI*rad),"Frequency") <<"\n";
|
|
G4cout<< "AT REST::: PARAMETERS\n"
|
|
<< "FrequencyG: " << G4BestUnit(gamma*gauss/(2*M_PI*rad),"Frequency") <<"\n";
|
|
#endif
|
|
|
|
|
|
rotation_angle = deltatime*omega;
|
|
|
|
|
|
Spin_rotation= G4Rotate3D(rotation_angle,B/Bnorm);
|
|
|
|
HepVector3D spin = aStep.GetTrack()->GetPolarization();
|
|
HepVector3D newspin;
|
|
newspin = Spin_rotation*spin;
|
|
|
|
G4double x,y,z,alpha;
|
|
x = sqrt(spin*spin);
|
|
y = sqrt(newspin*newspin);
|
|
z = spin*newspin/x/y;
|
|
alpha = acos(z);
|
|
|
|
#ifdef G4SRVERBOSE
|
|
G4cout<< "AT REST::: PARAMETERS\n"
|
|
<< "Initial spin : " << spin <<"\n"
|
|
<< "Delta time : " << deltatime <<"\n"
|
|
<< "Rotation angle: " << rotation_angle/(M_PI*rad) <<"\n"
|
|
<< "New spin : " << newspin <<"\n"
|
|
<< "Checked norms : " << x <<" " <<y <<" \n"
|
|
<< G4endl;
|
|
#endif
|
|
|
|
polar = newspin;
|
|
|
|
}
|
|
|
|
|
|
|
|
// Return the gyromagnetic ratio
|
|
/*!
|
|
* This method call the event manager to get the tracking manager and then the current track in order to get the particle name. According to the name of the particle it returns the gyromagnetic ratio.
|
|
*/
|
|
G4double LEMuSRAtRestSpinRotation::GetGyromagneticRatio()
|
|
{
|
|
//! Get the event manager
|
|
G4EventManager* evtMgr = G4EventManager::GetEventManager();
|
|
|
|
//! Get the track from the tracking manager
|
|
G4Track* theTrack = evtMgr->GetTrackingManager()->GetTrack();
|
|
|
|
//! Get the particle name
|
|
G4String particle = theTrack->GetDefinition()->GetParticleName();
|
|
|
|
/*! Arbitrary initialisation of \f$ \gamma \f$ as muons plus
|
|
* gyromagnetic ratio
|
|
*/
|
|
G4double gamma = 8.5062e+7*rad/(s*kilogauss);
|
|
|
|
//! Set gamma according to the particle. One can add other particles at will.
|
|
|
|
if(particle== "Mu")
|
|
{
|
|
gamma= gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light)));
|
|
}
|
|
|
|
else if (particle == "mu+")
|
|
{
|
|
gamma= 8.5062e+7*rad/(s*kilogauss);
|
|
}
|
|
|
|
else
|
|
{
|
|
gamma= 8.5062e+7*rad/(s*kilogauss);
|
|
}
|
|
|
|
return gamma;
|
|
}
|