LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add muonium spin rotation in flight.

Browse Source
This commit is contained in:
paraiso
2005-11-23 10:00:01 +00:00
parent 462d5ce591
commit 1a17b0b10e
1 changed files with 147 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

155
geant4/LEMuSR/src/LEMuSRMUONIUMScatt.cc
View File

@ -4,6 +4,10 @@
#include "G4TransportationManager.hh"
#include "Randomize.hh"
#include "G4ProductionCutsTable.hh"
#include "G4FieldManager.hh"
#include "G4Field.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -26,11 +30,11 @@ G4VParticleChange* LEMuSRMUONIUMScatt::PostStepDoIt(
fParticleChange.Initialize(trackData);
// tao :: Get Time
G4double itime = trackData.GetProperTime();
G4double gtime = trackData.GetGlobalTime();
G4double ftime = trackData.GetDynamicParticle()->GetPreAssignedDecayProperTime();
itime = trackData.GetProperTime();
gtime = trackData.GetGlobalTime();
ftime = trackData.GetDynamicParticle()->GetPreAssignedDecayProperTime();
G4double deltatime = ftime - itime;
deltatime = ftime - itime;
fParticleChange.ProposeGlobalTime(deltatime + itime -gtime);
// set position momentum energy
@ -46,15 +50,18 @@ G4VParticleChange* LEMuSRMUONIUMScatt::PostStepDoIt(
// fParticleChange.ProposeEnergy(0);
fParticleChange.ProposePosition(trackData.GetStep()->GetPreStepPoint()->GetPosition());
fParticleChange.ProposePolarization(trackData.GetPolarization());
// fParticleChange.ProposePolarization(trackData.GetPolarization());
// G4cout<<" at rest "<< trackData.GetStep()->GetPreStepPoint()->GetPosition() <<G4endl;
fParticleChange.ProposeTrackStatus(fStopButAlive) ;
}
}
else
{
fParticleChange.ProposePolarization(trackData.GetPolarization());
// fParticleChange.ProposePolarization(trackData.GetPolarization());
}
polar = RotateSpinIfMag(trackData, aStep);
fParticleChange.ProposePolarization(polar);
// fParticleChange.DumpInfo();
return &fParticleChange;
@ -78,7 +85,8 @@ G4bool LEMuSRMUONIUMScatt::CheckCondition( const G4Step& aStep)
return condition;
}
G4double LEMuSRMUONIUMScatt:: GetMeanFreePath(const G4Track& ,
G4double ,
@ -123,3 +131,134 @@ void LEMuSRMUONIUMScatt::PrepareSecondary(const G4Track& track)
;}
G4ThreeVector LEMuSRMUONIUMScatt::RotateSpinIfMag( const G4Track& theTrack, const G4Step& aStep)
{
G4ThreeVector theSpin = theTrack.GetPolarization() ;
if(theTrack.GetKineticEnergy()!=0&&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<<"IN FLIGHT::: MAGNETIC FIELD HERE" <<G4endl;// getchar();
#endif
// tao :: Get Field
point[0]=theTrack.GetPosition().x();
point[1]=theTrack.GetPosition().y();
point[2]=theTrack.GetPosition().z();
const G4Field* mfield;
mfield = fMgr->GetDetectorField();
mfield->GetFieldValue(point,B);
#ifdef G4SRVERBOSE
G4cout <<"IN FLIGHT::: MAGNETIC FIELD B="<< B[0]/gauss <<" G, " << B[1]/gauss <<" G, " << B[2]/gauss<<" G" <<G4endl;
#endif
#ifdef G4SRVERBOSE
G4cout <<"IN FLIGHT::: TIME= proper: "<< itime <<" ; global: " << gtime <<" decay: " << ftime <<G4endl;
#endif
G4ThreeVector magField(B[0],B[1],B[2]);
theSpin= RotateSpin(aStep,magField,deltatime);
#ifdef G4SRVERBOSE
G4cout<<"IN FLIGHT::: spin rotated";
#endif
}
}
return theSpin;
}
G4ThreeVector LEMuSRMUONIUMScatt::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<< "IN FLIGHT::: PARAMETERS\n"
<< "Magnetic Field Norm : " << G4BestUnit(Bnorm,"Magnetic flux density") <<"\n";
#endif
G4double omega,q,a,fqz;
G4double gamma;
q= aStep.GetTrack()->GetDefinition()->GetPDGCharge();
a= 1.165922e-3;
fqz = 8.5062e+7*rad/(s*kilogauss);
gamma = aStep.GetTrack()->GetDefinition()->GetGyromagneticRatio()*rad;
// G4cout<< fqz*(s*tesla)<<G4endl;
// G4cout<< gamma*(s*tesla)<<G4endl;
#ifdef G4SRVERBOSE
G4cout<< "IN FLIGHT::: PARAMETERS\n"
<< "Charge : " << q <<"\n";
#endif
// omega= - (fqz)*(1.+a) * Bnorm;
omega= - (gamma) * Bnorm;
#ifdef G4SRVERBOSE
G4cout<< "IN FLIGHT::: PARAMETERS\n"
<< "Frequency: " << G4BestUnit(fqz*gauss/(2*M_PI*rad),"Frequency") <<"\n";
G4cout<< "IN FLIGHT::: PARAMETERS\n"
<< "FrequencyG: " << G4BestUnit(gamma*gauss/(2*M_PI*rad),"Frequency") <<"\n";
#endif
G4double 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<< "IN FLIGHT::: 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
return newspin;
}