musrsim/geant4/spin_rot/src/sr1MuFormation.cc
2008-12-22 17:53:30 +00:00

107 lines
3.6 KiB
C++

//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium Formation according to yield.cc function (through GetYields method).
// Id : sr1MuFormation.cc, v 1.4
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#include "sr1MuFormation.hh"
using namespace std;
sr1MuFormation::sr1MuFormation(const G4String& name, G4ProcessType aType)
: G4VDiscreteProcess(name, aType){}
sr1MuFormation::~sr1MuFormation(){}
G4VParticleChange* sr1MuFormation::PostStepDoIt(const G4Track& trackData,
const G4Step& aStep)
{ // Initialize ParticleChange (by setting all its members equal to
// the corresponding members in G4Track)
fParticleChange.Initialize(trackData);
G4Track theNewTrack;
if(CheckCondition(aStep))
{
GetDatas(&aStep);
G4Step theStep;
PrepareSecondary( trackData);
fParticleChange.AddSecondary(aSecondary);
fParticleChange.ProposeTrackStatus(fStopAndKill) ;
}
else
{
fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ;
}
return &fParticleChange;
}
G4bool sr1MuFormation::CheckCondition(const G4Step& aStep)
{ // Decide when to call the MuFormation process - i.e. for muons going through the C foil.
G4bool condition=false;
p_name = aStep.GetTrack()->GetDefinition()->GetParticleName(); // particle name
//if(p_name == "mu+"&&aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName()=="log_CFoil")
std::string logVolName = aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
if(p_name == "mu+" && ((logVolName=="log_coulombCFoil")||(logVolName=="log_CFoil")))
{
condition=true;
}
return condition;
}
G4double sr1MuFormation::GetMeanFreePath(const G4Track&,
G4double,
G4ForceCondition* condition)
{
*condition = Forced;
return DBL_MAX;
}
void sr1MuFormation::GetDatas(const G4Step* aStep)
{ // Particle generation according to yield table
particleTable=G4ParticleTable::GetParticleTable();
rnd=G4UniformRand();
G4double E = aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy()/keV;
Gonin.GetYields(E,105.658369*1000,yvector); // Energy [keV], muon mass [keV/c2], yield table
G4String p_new = "Mu";
// Positive muon
if(p_name=="mu+")
{
if(rnd<yvector[0])
{
particle = particleTable->FindParticle(p_name) ;
}
else
{
particle = particleTable->FindParticle(p_new);
}
// Set the new dynamic particle DP
DP = new G4DynamicParticle(particle,
aStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection(),
aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy());
// IMPORTANT : COPY THOSE DATA TO GET THE SAME PARTICLE PROPERTIES!!!
// SHOULD BE KEPT WHEN BUILDING A PARTICLE CHANGE
DP->SetProperTime( aStep->GetTrack()->GetDynamicParticle()->GetProperTime());
DP->SetPolarization(aStep->GetTrack()->GetDynamicParticle()->GetPolarization().x(),
aStep->GetTrack()->GetDynamicParticle()->GetPolarization().y(),
aStep->GetTrack()->GetDynamicParticle()->GetPolarization().z());
DP->SetPreAssignedDecayProperTime(aStep->GetTrack()->GetDynamicParticle()->GetPreAssignedDecayProperTime());
}
}
void sr1MuFormation::PrepareSecondary(const G4Track& track)
{
if(p_name=="mu+")
{
aSecondary = new G4Track(DP,track.GetGlobalTime(),track.GetPosition());
}
}