| Message: Help setting up an Electric Field | Not Logged In (login) |
|
I have a xenon cylinder surrounded by air. I want to make an electric field in only the xenon. I tried borrowing the class from the field02 extended example. My field setup class is this:
#include "ExN02ElectricFieldSetup.hh" #include "ExN02FieldMessenger.hh"
#include "G4UniformElectricField.hh" #include "G4UniformMagField.hh" #include "G4MagneticField.hh" #include "G4FieldManager.hh" #include "G4TransportationManager.hh" #include "G4EquationOfMotion.hh" #include "G4EqMagElectricField.hh" #include "G4Mag_UsualEqRhs.hh" #include "G4MagIntegratorStepper.hh" #include "G4MagIntegratorDriver.hh" #include "G4ChordFinder.hh"
#include "G4ExplicitEuler.hh" #include "G4ImplicitEuler.hh" #include "G4SimpleRunge.hh" #include "G4SimpleHeum.hh" #include "G4ClassicalRK4.hh" #include "G4HelixExplicitEuler.hh" #include "G4HelixImplicitEuler.hh" #include "G4HelixSimpleRunge.hh" #include "G4CashKarpRKF45.hh" #include "G4RKG3_Stepper.hh" ////////////////////////////////////////////////////////////////////////// // // Constructors:
ExN02ElectricFieldSetup::ExN02ElectricFieldSetup()
: fChordFinder(0), fStepper(0), fIntgrDriver(0)
{
fEMfield = new G4UniformElectricField(
G4ThreeVector(0.0,1000.0*kilovolt/cm,0.0));
fFieldMessenger = new ExN02FieldMessenger(this) ;
fEquation = new G4EqMagElectricField(fEMfield);
fMinStep = 0.010*mm ; // minimal step of 10 microns
fStepperType = 4 ; // ClassicalRK4 -- the default stepper
fFieldManager = GetGlobalFieldManager(); UpdateField(); }
/////////////////////////////////////////////////////////////////////////////////
ExN02ElectricFieldSetup::ExN02ElectricFieldSetup(G4ThreeVector pFV)
: fChordFinder(0),
fStepper(0),
fIntgrDriver(0)
{
fEMfield = new G4UniformElectricField(pFV);
// GetGlobalFieldManager()->CreateChordFinder(this);
fFieldMessenger = new ExN02FieldMessenger(this) ; fEquation = new G4EqMagElectricField(fEMfield); fMinStep = 0.010*mm ; // minimal step of 10 microns fStepperType = 4 ; // ClassicalRK4 -- the default stepper
fFieldManager = GetGlobalFieldManager(); UpdateField(); } ////////////////////////////////////////////////////////////////////////////////
ExN02ElectricFieldSetup::~ExN02ElectricFieldSetup()
{
if(fChordFinder) delete fChordFinder;
if(fStepper) delete fStepper;
if(fEquation) delete fEquation;
if(fEMfield) delete fEMfield;
}
///////////////////////////////////////////////////////////////////////////// // // Register this field to 'global' Field Manager and // Create Stepper and Chord Finder with predefined type, minstep (resp.) //
void ExN02ElectricFieldSetup::UpdateField()
{
SetStepper();
G4cout<<"The minimal step is equal to "<<fMinStep/mm<<" mm"<<G4endl ;
fFieldManager->SetDetectorField(fEMfield );
if(fChordFinder) delete fChordFinder; // fChordFinder = new G4ChordFinder( fEMfield, fMinStep, fStepper);
fIntgrDriver = new G4MagInt_Driver(fMinStep,
fStepper,
fStepper->GetNumberOfVariables() );
fChordFinder = new G4ChordFinder(fIntgrDriver); fFieldManager->SetChordFinder( fChordFinder ); } ///////////////////////////////////////////////////////////////////////////// // // Set stepper according to the stepper type //
void ExN02ElectricFieldSetup::SetStepper()
{
G4int nvar = 8;
if(fStepper) delete fStepper;
switch ( fStepperType )
{
case 0:
fStepper = new G4ExplicitEuler( fEquation, nvar );
G4cout<<"G4ExplicitEuler is calledS"<<G4endl;
break;
case 1:
fStepper = new G4ImplicitEuler( fEquation, nvar );
G4cout<<"G4ImplicitEuler is called"<<G4endl;
break;
case 2:
fStepper = new G4SimpleRunge( fEquation, nvar );
G4cout<<"G4SimpleRunge is called"<<G4endl;
break;
case 3:
fStepper = new G4SimpleHeum( fEquation, nvar );
G4cout<<"G4SimpleHeum is called"<<G4endl;
break;
case 4:
fStepper = new G4ClassicalRK4( fEquation, nvar );
G4cout<<"G4ClassicalRK4 (default) is called"<<G4endl;
break;
case 5:
fStepper = new G4CashKarpRKF45( fEquation, nvar );
G4cout<<"G4CashKarpRKF45 is called"<<G4endl;
break;
case 6:
fStepper = 0; // new G4RKG3_Stepper( fEquation, nvar );
G4cout<<"G4RKG3_Stepper is not currently working for Electric Field"<<G4endl;
break;
case 7:
fStepper = 0; // new G4HelixExplicitEuler( fEquation );
G4cout<<"G4HelixExplicitEuler is not valid for Electric Field"<<G4endl;
break;
case 8:
fStepper = 0; // new G4HelixImplicitEuler( fEquation );
G4cout<<"G4HelixImplicitEuler is not valid for Electric Field"<<G4endl;
break;
case 9:
fStepper = 0; // new G4HelixSimpleRunge( fEquation );
G4cout<<"G4HelixSimpleRunge is not valid for Electric Field"<<G4endl;
break;
default: fStepper = 0;
}
}
///////////////////////////////////////////////////////////////////////////// // // Set the value of the Global Field to fieldValue along Z //
void ExN02ElectricFieldSetup::SetFieldValue(G4double fieldValue)
{
G4ThreeVector fieldVector( 0.0, 0.0, fieldValue );
SetFieldValue( fieldVector ); } /////////////////////////////////////////////////////////////////////////////// // // Set the value of the Global Field value to fieldVector //
void ExN02ElectricFieldSetup::SetFieldValue(G4ThreeVector fieldVector)
{
// Find the Field Manager for the global field
G4FieldManager* fieldMgr= GetGlobalFieldManager();
if(fieldVector != G4ThreeVector(0.,0.,0.))
{
if(fEMfield) delete fEMfield;
fEMfield = new G4UniformElectricField(fieldVector);
fEquation->SetFieldObj(fEMfield); // must now point to the new field
// UpdateField();
fieldMgr->SetDetectorField(fEMfield);
}
else
{
// If the new field's value is Zero, then it is best to
// insure that it is not used for propagation.
if(fEMfield) delete fEMfield;
fEMfield = 0;
fEquation->SetFieldObj(fEMfield); // As a double check ...
G4MagneticField* fEMfield = 0;
fieldMgr->SetDetectorField(fEMfield);
}
}
////////////////////////////////////////////////////////////////////////////////
//
// Utility method
G4FieldManager* ExN02ElectricFieldSetup::GetGlobalFieldManager()
{
return G4TransportationManager::GetTransportationManager()
->GetFieldManager();
}
I added the line: ExN02ElectricFieldSetup* field = new ExN02ElectricFieldSetup(); to my main. When the simulation starts, it seems as though particles are leaving the xenon and keep going on forever until they exit the world volume, which gives me a segmentation fault. I know that to isolate the field to only the xenon cylinder, I need to use the line: xenon_log->SetFieldManager(field->GetGlobalFieldManager(), false); Is this right? If it is right, in what class do I put this line? The documentation isn't very specific about this. |
| Inline Depth: | Outline Depth: | Add message: |
|
to: |