musrsim/geant4/LEMuSR/src/FocalLengthTest.cc
2006-02-16 17:20:45 +00:00

267 lines
6.3 KiB
C++

//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
// 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();
}