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

Fixed Bugs

Browse Source
This commit is contained in:
paraiso
2006-03-15 02:37:11 +00:00
parent da93a5a456
commit 299699b906
25 changed files with 272 additions and 341 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
geant4/LEMuSR/src/LEMuSRAtRestSpinRotation.cc
View File

@@ -82,7 +82,7 @@ G4VParticleChange* LEMuSRAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack,
{
theParticleChange.Initialize(theTrack);
theParticleChange.ProposeTrackStatus(fAlive);
theParticleChange.ProposeTrackStatus(fStopButAlive);
theParticleChange.ProposePosition(theTrack.GetPosition());
theParticleChange.ProposeMomentumDirection(theTrack.GetMomentumDirection());
theParticleChange.ProposeEnergy(theTrack.GetKineticEnergy());
@@ -128,24 +128,27 @@ G4VParticleChange* LEMuSRAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack,
#ifdef G4SRVERBOSE
G4cout <<"AT REST::: MAGNETIC FIELD B="<< B[0] <<" " << B[1] <<" " << B[2] <<G4endl;
#endif
#endif
#ifdef G4SRVERBOSE
G4cout <<"AT REST::: TIME= proper: "<< itime <<" ; global: " << gtime <<" decay: " << ftime <<G4endl;
#endif
G4ThreeVector magField(B[0],B[1],B[2]);
RotateSpin(aStep,magField,deltatime);
if(magField!=G4ThreeVector(0.0,0.0,0.0))
{
RotateSpin(aStep,magField,deltatime);
#ifdef G4SRVERBOSE
G4cout<<"AT REST::: spin rotated";
#endif
}
}
theParticleChange.ProposePolarization(polar);
theParticleChange.ProposePolarization(polar);

9
geant4/LEMuSR/src/LEMuSRCryoSD.cc
View File

@@ -197,6 +197,13 @@ void LEMuSRCryoSD::getHit()
theHit.positron = e;
theHit.gamma = g;
theHit.runid = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
theParticle.muon = mu;
theParticle.positron = e;
theParticle.gamma = g;
theParticle.runid = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
}
// ROOT METHODS
@@ -211,6 +218,8 @@ void LEMuSRCryoSD::BookRoot()
tree->Branch("cryoHit",&theHit.kenergy,"kenergy/F:tenergy/F:edeposit/F:localtime/F:globaltime:propertime/F:positionx/F:positiony:positionz:momdirx:momdiry:momdirz/F:foil/F:muon/I:muonium/I:gamma/I:runID/I");
tree->Branch("particle",&theParticle.muon,"muon/I:muonium/I:gamma/I:runID/I");
}
void LEMuSRCryoSD::FillRoot()

42
geant4/LEMuSR/src/LEMuSRDetectorConstruction.cc
View File

@@ -42,6 +42,7 @@
#include "G4PropagatorInField.hh"
#include "G4ChordFinder.hh"
#include "G4El_UsualEqRhs.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ClassicalRK4.hh"
#include "G4SimpleHeum.hh"
#include "G4SimpleRunge.hh"
@@ -169,13 +170,13 @@ G4VPhysicalVolume* LEMuSRDetectorConstruction::lemuDetector() // !mind the V in
//! Definition of the top mothe volume: the laboratory. Refered to as Wolrd.
// +++++++++++++++++++++++++++++++DEFINE THE MOTHER VOLUME: THE LABORATOY++++++++++++++++++++++
// solid
G4double LABO_x = 1*m;
G4double LABO_y = 1*m;
G4double LABO_z = 2*m;
G4double LABO_x = 2*m;
G4double LABO_y = 2*m;
G4double LABO_z = 4*m;
LABO_box = new G4Box("World_box",LABO_x,LABO_y,LABO_z);
// logical volume
LABO_material = G4Material::GetMaterial("Air");
LABO_material = G4Material::GetMaterial("vacuum");
lv_LABO = new G4LogicalVolume(LABO_box,LABO_material,"lv_World",0,0,0);
// physical volume
@@ -354,7 +355,7 @@ void LEMuSRDetectorConstruction::lemuMCP2()
// gate valve chamber
GATV_tube= new G4Tubs( "sGATV",7.65*cm, 10.325*cm, 9.25*cm, dSPhi, dEPhi);
GATV_tube= new G4Tubs( "sGATV",7.655*cm, 10.32*cm, 9.25*cm, dSPhi, dEPhi);
GATS_tube= new G4Tubs( "sGATS", 10.325*cm, 12.65*cm, 9.25*cm, dSPhi, dEPhi);
G4Transform3D tr_gatv;
@@ -391,13 +392,12 @@ void LEMuSRDetectorConstruction::lemuMCP2()
if(AsymCheck!=1){
pv_MCPS = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, 0*cm),lv_MCPS ,"pv_MCPS",lv_LABO, false, 0 );
pv_F160 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, -15.10*cm),lv_F160 ,"pv_F160",lv_LABO, false, 0 );
pv_F100 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, +16.2*cm),lv_F100 ,"pv_F100",lv_LABO, false, 0 );
pv_F100 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, +16.2*cm),lv_F100 ,"pv_F100",lv_MCPV, false, 0 );
G4VPhysicalVolume* pv_F101;
pv_F101 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, -124.*cm),lv_F100 ,"pv_F101",lv_LABO, false, 0 );
GATVz= -25.45;
pv_F200 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, 8.05*cm),lv_F200 ,"pv_F200",lv_GATV, false, 0 );
@@ -416,7 +416,7 @@ void LEMuSRDetectorConstruction::lemuMCP2()
// Apply attributes and user limit
lv_MCPV->SetVisAttributes(G4VisAttributes::Invisible);
lv_MCPV->SetVisAttributes(VTBB_style);//G4VisAttributes::Invisible);
lv_MCPS->SetVisAttributes(Blue_style);
lv_F160->SetVisAttributes(fBlue_style);
lv_F100->SetVisAttributes(fBlue_style);
@@ -428,7 +428,6 @@ void LEMuSRDetectorConstruction::lemuMCP2()
// user limits
G4UserLimits* RAV_lim = new G4UserLimits();
RAV_lim -> SetMaxAllowedStep(FieldStepLim );
//RAV_lim->SetUserMaxTrackLength(1.2*m);
RAV_lim-> SetUserMinEkine(1.*eV);
lv_MCPV->SetUserLimits(RAV_lim);
//
@@ -453,13 +452,13 @@ void LEMuSRDetectorConstruction::lemuANODE()
RA_Ebox = new G4Box("RA_Ebox", 1.*cm,12.5*cm,2.25*cm );
RA_Mbox = new G4Box("RA_Mbox", 1.*cm,12.5*cm,2.25*cm );
RA_G_tube= new G4Tubs( "sRA_G", 5.8*cm,6.25*cm,5.55*cm, dSPhi, dEPhi);
RA_G_tube = new G4Tubs( "sRA_G", 5.8*cm,6.25*cm,5.55*cm, dSPhi, dEPhi);
// logical volumes
lv_RA_E = new G4LogicalVolume( RA_E_cone , SSteel ,"lv_RA_E",0,0,0);
lv_RA_M = new G4LogicalVolume( RA_M_cone , SSteel ,"lv_RA_M",0,0,0);
lv_RA_G = new G4LogicalVolume( RA_G_tube , SSteel ,"lv_RA_G",0,0,0);
lv_RA_E = new G4LogicalVolume( RA_E_cone ,G4Material::GetMaterial("stainless_steel") ,"lv_RA_E",0,0,0);
lv_RA_M = new G4LogicalVolume( RA_M_cone ,G4Material::GetMaterial("stainless_steel") ,"lv_RA_M",0,0,0);
lv_RA_G = new G4LogicalVolume( RA_G_tube , G4Material::GetMaterial("copper") ,"lv_RA_G",0,0,0);
// physical volumes
@@ -469,8 +468,12 @@ void LEMuSRDetectorConstruction::lemuANODE()
// Mother Volume is lv_MCPV!!!!
pv_RA_E = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Ez*cm-mcpv_z),lv_RA_E ,"pv_RA_E",lv_MCPV, false, 0 );
pv_RA_M = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Mz*cm-mcpv_z),lv_RA_M ,"pv_RA_M",lv_MCPV, false, 0 );
if(halfview==0)
{
pv_RA_E = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Ez*cm-mcpv_z),lv_RA_E ,"pv_RA_E",lv_MCPV, false, 0 );
pv_RA_M = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Mz*cm-mcpv_z),lv_RA_M ,"pv_RA_M",lv_MCPV, false, 0 );
}
pv_RA_G = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Gz*cm-mcpv_z),lv_RA_G ,"pv_RA_G",lv_MCPV, false, 0 );
rotation1=G4RotateZ3D(180*deg) ;
@@ -493,7 +496,6 @@ void LEMuSRDetectorConstruction::lemuANODE()
lv_RA_E->SetVisAttributes(lBlue_style);
lv_RA_M->SetVisAttributes(dBlue_style);
lv_RA_G->SetVisAttributes(lBlue_style);
}
@@ -534,7 +536,7 @@ void LEMuSRDetectorConstruction:: lemuCRYO()
Guard_Rings = new G4Tubs ("sGuard",3.0*cm,3.5*cm,0.25*cm, dSPhi, dEPhi);
// logicals volumes
// logicals volumes
// sample MAGNETIC field
@@ -1164,7 +1166,7 @@ void LEMuSRDetectorConstruction::buildAnodeField()
// 2 case of electro magnetic field
G4UniformMagField* mcField;
G4MagneticField* mcField;
if(magfield==1)
{
@@ -1235,7 +1237,7 @@ void LEMuSRDetectorConstruction::buildAnodeField()
if(magfield==1)
{
G4UniformMagField* mcField;
G4MagneticField* mcField;
#ifdef LEMU_TRANSVERSE_FIELD

12
geant4/LEMuSR/src/LEMuSRDetectorMessenger.cc
View File

@@ -403,8 +403,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu
{
if(newvalue=="total")
{
theDetector->dSPhi=0*deg;
theDetector->dEPhi=360*deg;
theDetector->dSPhi=0.*deg;
theDetector->dEPhi=360.*deg;
theDetector->halfview=0;
newDetector = theDetector->Construct();
@@ -413,8 +413,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu
}
else if(newvalue=="half")
{
theDetector->dSPhi=+90*deg;
theDetector->dEPhi=180*deg;
theDetector->dSPhi=90.*deg;
theDetector->dEPhi=180.*deg;
theDetector->halfview=1;
newDetector = theDetector->Construct();
@@ -423,8 +423,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu
}
else if(newvalue=="quarter")
{
theDetector->dSPhi=0*deg;
theDetector->dEPhi=270*deg;
theDetector->dSPhi=0.*deg;
theDetector->dEPhi=270.*deg;
theDetector->halfview=1;
newDetector = theDetector->Construct();

103
geant4/LEMuSR/src/LEMuSRIonPhysics.cc
View File

@@ -1,103 +0,0 @@
// $Id$
// GEANT4 tag $Name: $
//
//
#include "LEMuSRIonPhysics.hh"
#include "globals.hh"
#include "G4ios.hh"
#include <iomanip>
LEMuSRIonPhysics::LEMuSRIonPhysics(const G4String& name)
: G4VPhysicsConstructor(name)
{
}
LEMuSRIonPhysics::~LEMuSRIonPhysics()
{
}
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
// Nuclei
#include "G4IonConstructor.hh"
void LEMuSRIonPhysics::ConstructParticle()
{
// Construct light ions
G4IonConstructor pConstructor;
pConstructor.ConstructParticle();
}
#include "G4ProcessManager.hh"
void LEMuSRIonPhysics::ConstructProcess()
{
G4ProcessManager * pManager = 0;
// Elastic Process
theElasticModel = new G4LElastic();
theElasticProcess.RegisterMe(theElasticModel);
// Generic Ion
pManager = G4GenericIon::GenericIon()->GetProcessManager();
// add process
pManager->AddDiscreteProcess(&theElasticProcess);
pManager->AddProcess(&fIonMultipleScattering, -1, 1, 1);
pManager->AddProcess(&fIonIonisation, -1, 2, 2);
// Deuteron
pManager = G4Deuteron::Deuteron()->GetProcessManager();
// add process
pManager->AddDiscreteProcess(&theElasticProcess);
fDeuteronModel = new G4LEDeuteronInelastic();
fDeuteronProcess.RegisterMe(fDeuteronModel);
pManager->AddDiscreteProcess(&fDeuteronProcess);
pManager->AddProcess(&fDeuteronMultipleScattering, -1, 1, 1);
pManager->AddProcess(&fDeuteronIonisation, -1, 2, 2);
// Triton
pManager = G4Triton::Triton()->GetProcessManager();
// add process
pManager->AddDiscreteProcess(&theElasticProcess);
fTritonModel = new G4LETritonInelastic();
fTritonProcess.RegisterMe(fTritonModel);
pManager->AddDiscreteProcess(&fTritonProcess);
pManager->AddProcess(&fTritonMultipleScattering, -1, 1, 1);
pManager->AddProcess(&fTritonIonisation, -1, 2, 2);
// Alpha
pManager = G4Alpha::Alpha()->GetProcessManager();
// add process
pManager->AddDiscreteProcess(&theElasticProcess);
fAlphaModel = new G4LEAlphaInelastic();
fAlphaProcess.RegisterMe(fAlphaModel);
pManager->AddDiscreteProcess(&fAlphaProcess);
pManager->AddProcess(&fAlphaMultipleScattering, -1, 1, 1);
pManager->AddProcess(&fAlphaIonisation, -1, 2, 2);
// He3
pManager = G4He3::He3()->GetProcessManager();
// add process
pManager->AddDiscreteProcess(&theElasticProcess);
pManager->AddProcess(&fHe3MultipleScattering, -1, 1, 1);
pManager->AddProcess(&fHe3Ionisation, -1, 2, 2);
}

2
geant4/LEMuSR/src/LEMuSRMUONIUM.cc
View File

@@ -67,7 +67,7 @@ G4VParticleChange* LEMuSRMUONIUM::PostStepDoIt(
{
fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ;
}
// G4cout<<"MU Formation"<<G4endl;getchar();
return &fParticleChange;
}

14
geant4/LEMuSR/src/LEMuSRMag_SpinEqRhs.cc
View File

@@ -63,13 +63,10 @@ LEMuSRMag_SpinEqRhs::SetChargeMomentumMass(G4double particleCharge, // in e+ uni
G4double gratio=GetGyromagneticRatio();
//G4cout <<"g ratio [MHz]/[T]"<<gratio<<G4endl;
omegac =gratio;// 0.105658387*GeV/mass * 2.837374841e-3*(rad/cm/kilogauss);
oldomegac = 0.105658387*GeV/mass * 2.837374841e-3*(rad/cm/kilogauss);
anomaly = 1.165923e-3;
omegac =gratio;
#ifdef DEBUG
anomaly = 1.165923e-3;
oldomegac= 0.105658387*GeV/mass * 2.837374841e-3*(rad/cm/kilogauss);
//G4cout<< "Old FrequencyG: " << G4BestUnit(oldomegac*gauss/(2*M_PI*rad)*cm,"Frequency") <<G4endl;
#endif
@@ -92,6 +89,13 @@ LEMuSRMag_SpinEqRhs::EvaluateRhsGivenB( const G4double y[],
G4double dydx[] ) const
{
G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]);
if (momentum_mag_square==0)
{
G4cout<<"LEMuSRMag_SpinEqRhs says: SEVERE ERROR:: MOMENTUM SQUARE ==0 !!!"<<G4endl;
return ;
}
G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
G4double cof = FCof()*inv_momentum_magnitude;
// G4cout << "\n------------ LEMuSRMAg_SpinEqRhs :: mommagnitude : "<< sqrt(momentum_mag_square) <<"\n";

37
geant4/LEMuSR/src/LEMuSRMagneticField.cc
View File

@@ -26,7 +26,6 @@
LEMuSRMagneticField::LEMuSRMagneticField(const G4ThreeVector FieldVector)
:G4UniformMagField(FieldVector )
{
BField=FieldVector;
// G4cout<<"MAGNETIC FIELD DEFINED AS "<<BField.x()/gauss<<"gauss"<<BField.y()/gauss<< "gauss"<<BField.z()/gauss<< "gauss"<<G4endl;
@@ -42,27 +41,25 @@ void LEMuSRMagneticField::GetFieldValue (const G4double pos[4],
G4double *field ) const
{
field[0]= 0.0;
field[1]= 0.0;
field[2]= 0.0;
field[0]= 0.0*gauss;
field[1]= 0.0*gauss;
field[2]= 0.0*gauss;
G4double X,Y,Z,factor;
G4double X,Y,Z,factor;
X= pos[0]*mm;Y=pos[1]*mm;Z=pos[2]*mm;
// G4cout<<"\n"<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<G4endl;
if(Z<20*cm&&Z>-20*cm)
{ //G4cout<<"true!";
factor=exp((-Z*Z)/(10*cm*10*cm));
field[0]= BField.x()*factor ;//TAO
field[1]= BField.y()*factor ;
field[2]= BField.z()*factor ;
//G4cout<<"true!"<<field[0]/gauss<<"gauss"<<field[1]/gauss<< "gauss"<<field[2]/gauss<< "gauss"<<G4endl;getchar();
// getchar();
}
// G4cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<G4endl;
/* if(Z<20*cm&&Z>-20*cm)
{ //G4cout<<"true!";*/
factor=exp((-Z*Z)/(10*cm*10*cm));
field[0]= BField.x()*factor*0 ;//TAO
field[1]= BField.y()*factor*0 ;
field[2]= BField.z()*factor*0 ;
// G4cout<<"true: "<<field[0]/gauss<<" gauss, "<<field[1]/gauss<< " gauss, "<<field[2]/gauss<< " gauss."<<G4endl;// getchar();
// getchar();
/*
}
*/
}

2
geant4/LEMuSR/src/LEMuSRPgaMessenger.cc
View File

@@ -56,7 +56,7 @@ LEMuSRPgaMessenger::LEMuSRPgaMessenger(LEMuSRPrimaryGeneratorAction *thePGA)
posCmd->SetDefaultValue(Hep3Vector(0., 0., -114));
momCmd = new G4UIcmdWith3Vector("/lemuGun/MomentumDirection",this);
momCmd->SetGuidance("Set momentum direction >> norm must be one.");
momCmd->SetGuidance("Set momentum direction >> sum must equal one.");
momCmd->SetParameterName("momX ","momY ","momZ ",true,true);
momCmd->SetDefaultValue(Hep3Vector(0., 0., 1.));

92
geant4/LEMuSR/src/LEMuSRPhysicsList.cc
View File

@@ -51,8 +51,8 @@
LEMuSRPhysicsList::LEMuSRPhysicsList(): G4VUserPhysicsList()
{
defaultCutValue = 1.0*mm;
SetVerboseLevel(1);
defaultCutValue = 5.0*mm;
SetVerboseLevel(2);
}
LEMuSRPhysicsList::~LEMuSRPhysicsList()
@@ -78,11 +78,13 @@ void LEMuSRPhysicsList::ConstructParticle()
// for all particles which you want to use.
// This ensures that objects of these particle types will be
// created in the program.
SetCuts();
ConstructBosons();
ConstructLeptons();
ConstructBaryons();
// ConstructMesons();
}
@@ -169,12 +171,14 @@ void LEMuSRPhysicsList::ConstructBaryons()
#include "G4PhotoElectricEffect.hh"
#include "G4MultipleScattering.hh"
#include "G4MultipleScattering52.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4MuIonisation.hh"
#include "G4MuIonisation52.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
@@ -210,7 +214,9 @@ void LEMuSRPhysicsList::ConstructProcess()
*/
ConstructEM();
/*! Construction of the Decay process when applicable. */
/*! Construction of the Decay process when applicable.
*
* Construction of User Cuts and Step Limitation. */
ConstructGeneral();
}
@@ -291,7 +297,7 @@ void LEMuSRPhysicsList::ConstructEM()
}
/*! \anchor muplusphysics*/
//! Muon Plus Physics is implemented here.
else if( particleName == "mu+" ) {
@@ -299,7 +305,7 @@ void LEMuSRPhysicsList::ConstructEM()
// PROCESSES FOR MUON PLUS
//! Multiple Scattering including Meyer'Algorithm: LEMuSRMSC
G4VProcess* aMultipleScattering = new LEMuSRMSC();
G4VProcess* aMultipleScattering = new LEMuSRMSC();//new G4MultipleScattering52();
//! Muonium formation: LEMuSRMUONIUM
G4VProcess* aMuoniumFormation = new LEMuSRMUONIUM();
@@ -313,7 +319,7 @@ void LEMuSRPhysicsList::ConstructEM()
G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung();
G4VProcess* aPairProduction = new G4MuPairProduction();
G4VProcess* anIonisation = new G4MuIonisation();
// ADD WANTED PROCESSES TO PROCESS MANAGER
@@ -332,20 +338,20 @@ void LEMuSRPhysicsList::ConstructEM()
// SET PROCESSES ORDERING (ALONG STEP. POST STEP. AT REST).
// set ordering for AlongStepDoIt
pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,1);
pmanager->SetProcessOrdering(anIonisation, idxAlongStep,2);
pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,3);
pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,4);
pmanager->SetProcessOrdering(anIonisation, idxAlongStep,3);
pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,4);
pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,5);
// set ordering for PostStepDoIt
pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep,1);
pmanager->SetProcessOrdering(aMuoniumFormation, idxPostStep,2);
pmanager->SetProcessOrdering(anIonisation, idxPostStep,3);
pmanager->SetProcessOrdering(anIonisation, idxPostStep,3);
pmanager->SetProcessOrdering(aBremsstrahlung, idxPostStep,4);
pmanager->SetProcessOrdering(aPairProduction, idxPostStep,5);
// set ordering for AtRestDoIt
pmanager->SetProcessOrdering(aRestSpinRotation, idxAtRest,1);
pmanager->SetProcessOrderingToFirst(aRestSpinRotation, idxAtRest);
}
@@ -360,42 +366,19 @@ void LEMuSRPhysicsList::ConstructEM()
//! Spin Precession Process At Rest: LEMuSRAtRestSpinRotation
G4VProcess* aRestSpinRotation = new LEMuSRAtRestSpinRotation();
G4VProcess* aMuScatt = new LEMuSRMUONIUMScatt();
/*
//! Unchanged \gf processes
G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung();
G4VProcess* aPairProduction = new G4MuPairProduction();
G4VProcess* anIonisation = new G4MuIonisation();
*/
// ADD WANTED PROCESSES TO PROCESS MANAGER
// add lemu processes
pmanager->AddProcess(aMuScatt);
pmanager->AddProcess(aRestSpinRotation);
// add \gf processes
// pmanager->AddProcess(aBremsstrahlung);
// pmanager->AddProcess(anIonisation);
// pmanager->AddProcess(aPairProduction);
// SET PROCESSES ORDERING (ALONG STEP. POST STEP. AT REST).
// set ordering for AlongStepDoIt
// pmanager->SetProcessOrdering(anIonisation, idxAlongStep,1);
// pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,2);
// pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,3);
// set ordering for PostStepDoIt
pmanager->SetProcessOrdering(aMuScatt, idxPostStep,1);
// pmanager->SetProcessOrdering(anIonisation, idxPostStep,2);
// pmanager->SetProcessOrdering(aBremsstrahlung, idxPostStep,3);
// pmanager->SetProcessOrdering(aPairProduction, idxPostStep,4);
// set ordering for AtRestDoIt
pmanager->SetProcessOrdering(aRestSpinRotation, idxAtRest,1);
pmanager->SetProcessOrderingToFirst(aRestSpinRotation, idxAtRest);
}
@@ -405,7 +388,7 @@ void LEMuSRPhysicsList::ConstructEM()
//! Muon Minus Physics is implemented here.
else if( particleName == "mu-" ) {
//muon minus
G4VProcess* aMultipleScattering = new G4MultipleScattering();
G4VProcess* aMultipleScattering = new G4MultipleScattering52();
G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung();
G4VProcess* aPairProduction = new G4MuPairProduction();
G4VProcess* anIonisation = new G4MuIonisation();
@@ -468,7 +451,8 @@ void LEMuSRPhysicsList::ConstructEM()
#include "G4Decay.hh"
#include "G4UserSpecialCuts.hh"
#include "G4DecayTable.hh"
// Step Limiter Process
#include "G4StepLimiter.hh"
/*! \anchor constructdecayprocess*/
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
@@ -481,9 +465,12 @@ void LEMuSRPhysicsList::ConstructEM()
void LEMuSRPhysicsList::ConstructGeneral()
{
/*! Declare Decay Processes (from \lemu and \gf).*/
G4Decay* theDecayProcess = new G4Decay();
// G4Decay* theDecayProcess = new G4Decay();
LEMuSRDecay* theLEMuSRDecayProcess = new LEMuSRDecay();
/*! Declare User Cut.*/
G4StepLimiter* aStepLimiter = new G4StepLimiter();
G4UserSpecialCuts* aUserCut = new G4UserSpecialCuts();
/*! Initialization of the particle iterator.*/
theParticleIterator->reset();
@@ -512,10 +499,14 @@ void LEMuSRPhysicsList::ConstructGeneral()
if (theLEMuSRDecayProcess->IsApplicable(*particle))
{
pmanager ->AddProcess(theLEMuSRDecayProcess);
// set ordering for PostStepDoIt and AtRestDoIt
pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxPostStep);
pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxAtRest);
pmanager ->AddProcess(aStepLimiter);
// set ordering
pmanager ->SetProcessOrderingToLast(aStepLimiter, idxPostStep);
pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep);
pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest);
}
pmanager->DumpInfo();
}
else if (particle->GetParticleName()=="Mu")
{
@@ -534,11 +525,15 @@ void LEMuSRPhysicsList::ConstructGeneral()
if (theLEMuSRDecayProcess->IsApplicable(*particle))
{
pmanager ->AddProcess(theLEMuSRDecayProcess);
// set ordering for PostStepDoIt and AtRestDoIt
pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxPostStep);
pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxAtRest);
pmanager ->AddProcess(aStepLimiter);
// set ordering
// pmanager ->SetProcessOrdering(aStepLimiter, idxPostStep, 1); // not needed for the muonium
pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep);
pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest);
}
pmanager->DumpInfo();
}
@@ -553,6 +548,8 @@ void LEMuSRPhysicsList::ConstructGeneral()
else
{
G4ProcessManager* pmanager = particle->GetProcessManager();
pmanager ->AddProcess(aUserCut);
pmanager ->SetProcessOrdering(aUserCut, idxPostStep);
if (theLEMuSRDecayProcess->IsApplicable(*particle))
{
pmanager ->AddProcess(theLEMuSRDecayProcess);
@@ -580,6 +577,9 @@ void LEMuSRPhysicsList::SetCuts()
// These values are used as the default production thresholds
// for the world volume.
SetCutsWithDefault();
SetCutValue(1.*mm, "e+");
SetCutValue(1.*mm, "e-");
// SetCutValue(1*mm, "proton");
}

2
geant4/LEMuSR/src/LEMuSRPrimaryGeneratorAction.cc
View File

@@ -78,7 +78,7 @@ LEMuSRPrimaryGeneratorAction::LEMuSRPrimaryGeneratorAction()
lemuParticleGun->SetNumberOfParticles(1);
X=0; Y=0;
Z=0.;
Z=-114.;
momX=0.; momY=0.; momZ=1;
scan=0;

44
geant4/LEMuSR/src/LEMuSRStackingAction.cc
View File

@@ -38,6 +38,8 @@ LEMuSRStackingAction::LEMuSRStackingAction()
:ScintHits(0), stage(0)
{
theMessenger = new LEMuSRStackingActionMessenger(this);
kill_e = true; kill_gamma = true; kill_nu_e = kill_nu_mu = true;
}
LEMuSRStackingAction::~LEMuSRStackingAction()
@@ -48,10 +50,50 @@ LEMuSRStackingAction::~LEMuSRStackingAction()
G4ClassificationOfNewTrack
LEMuSRStackingAction::ClassifyNewTrack(const G4Track *)
LEMuSRStackingAction::ClassifyNewTrack(const G4Track *track)
{
G4ClassificationOfNewTrack classification;
classification = fUrgent;
G4String p_name;
p_name = track->GetDefinition()->GetParticleName(); // particle name
if (p_name=="e-")
{
if (kill_e) classification=fKill;
}
else if (p_name=="gamma")
{
if (kill_gamma) classification=fKill;
}
else if (p_name=="nu_e")
{
if (kill_nu_e) classification=fKill;
}
else if (p_name=="anti_nu_e")
{
if (kill_nu_e) classification=fKill;
}
else if (p_name=="nu_mu")
{
if (kill_nu_mu) classification=fKill;
}
else if (p_name=="anti_nu_mu")
{
if (kill_nu_mu) classification=fKill;
}
return classification;
}

29
geant4/LEMuSR/src/LEMuSRSteppingAction.cc
View File

@@ -66,6 +66,7 @@ void LEMuSRSteppingAction::UserSteppingAction(const G4Step* aStep)
{
LoopKiller(aStep);
// G4cout<<"LEMuSR STEPPING ACTION INVOKED!"<<G4endl ;
G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
G4String p_name;
@@ -87,6 +88,9 @@ void LEMuSRSteppingAction::UserSteppingAction(const G4Step* aStep)
pVVisManager -> Draw(polyline);
}
// ParticleInfo(aStep);
}
@@ -105,23 +109,14 @@ void LEMuSRSteppingAction::UserSteppingAction(const G4Step* aStep)
void LEMuSRSteppingAction::LoopKiller(const G4Step *aStep)
{
// loop killa
if(aStep->GetTrack()->GetCurrentStepNumber()>2500)
if(aStep->GetTrack()->GetCurrentStepNumber()>1500)
{
aStep->GetTrack()->SetTrackStatus(fStopAndKill);
G4cout<<"killed "<<G4endl;
G4cout<<"killed:: step number > 1500 "<<G4endl;
ParticleInfo(aStep);
// getchar();
}
/* if(aStep->GetTrack()->GetDefinition()->GetParticleName()=="e-"
||aStep->GetTrack()->GetDefinition()->GetParticleName()=="gamma"
||aStep->GetTrack()->GetDefinition()->GetParticleName()=="nu_e"
||aStep->GetTrack()->GetDefinition()->GetParticleName()=="anti_nu_e"
||aStep->GetTrack()->GetDefinition()->GetParticleName()=="nu_mu"
||aStep->GetTrack()->GetDefinition()->GetParticleName()=="anti_nu_mu")
{
aStep->GetTrack()->SetTrackStatus(fStopAndKill);
}
*/
}
@@ -246,8 +241,8 @@ void LEMuSRSteppingAction:: ParticleInfo(const G4Step *aStep)
G4double spin= aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetPDGSpin(); // spin in units of 1
//b
G4ThreeVector position = aStep->GetPostStepPoint()->GetPosition(); // position
G4ThreeVector momentum = aStep->GetPostStepPoint()->GetMomentum(); // momentum
G4ThreeVector position = aStep->GetTrack()->GetPosition(); // position
G4ThreeVector momentum = aStep->GetTrack()->GetMomentumDirection(); // momentum
//c
G4double tof = aStep->GetTrack()->GetLocalTime(); // time since track creation
G4double globaltime = aStep->GetTrack()->GetGlobalTime();// time since the event in which the track belongs is created.
@@ -263,8 +258,8 @@ void LEMuSRSteppingAction:: ParticleInfo(const G4Step *aStep)
// G4cout << "Parent particles are : " << aStep->GetTrack()->GetCreatorProcess()->GetProcessName() <<" ;\n ";
}
G4cout << "particle name : " << p_name <<" ;\n "
G4cout.precision(6);
G4cout << "particle name : " << p_name <<" ;\n "
<< "volume name : " << v_name <<" ;\n "
<< "next volume name : " << nv_name <<" ;\n "
<< "spin : " << spin <<" ;\n "

9
geant4/LEMuSR/src/LEMuSRTrack.cc
View File

@@ -1,9 +0,0 @@
#include "LEMuSRTrack.hh"
void LEMuSRTrack::SetDynamicParticle(G4DynamicParticle particle)
{
G4DynamicParticle *theParticle;
// theParticle = this->GetDynamicParticle();
*theParticle = particle;
}

8
geant4/LEMuSR/src/LEMuSRdummydets.cc
View File

@@ -96,10 +96,12 @@ void LEMuSRDetectorConstruction::lemuAsym()
G4Sphere *Asym_sphr = new G4Sphere("sphAsymR",25*cm,25.000001*cm,0.*deg,360*deg,0*deg,90*deg);
// logical volumes
lv_Asym = new G4LogicalVolume(Asym_tube,G4Material::GetMaterial("vacuum"),"lv_Asym",0,0,0);
G4String mat = "vacuum";
lv_AsymL = new G4LogicalVolume(Asym_sphl,G4Material::GetMaterial("vacuum"),"lv_AsymL",0,0,0);
lv_AsymR = new G4LogicalVolume(Asym_sphr,G4Material::GetMaterial("vacuum"),"lv_AsymR",0,0,0);
lv_Asym = new G4LogicalVolume(Asym_tube,G4Material::GetMaterial(mat),"lv_Asym",0,0,0);
lv_AsymL = new G4LogicalVolume(Asym_sphl,G4Material::GetMaterial(mat),"lv_AsymL",0,0,0);
lv_AsymR = new G4LogicalVolume(Asym_sphr,G4Material::GetMaterial(mat),"lv_AsymR",0,0,0);
// G4Translate3D trl =G4ThreeVector(+50.*cm,0.*cm,0.*cm); //uncomment
Asym_rotation = G4RotateY3D(90*deg);