Forum: Fields: Magnetic and Otherwise Not Logged In (login)
Forums by Category
Forums by Time Order
Request a New Forum
Recent Postings
Search in Forums
Subscribe to Forums
Member Info
Members List
New Member
Overview
Contact Admin
Show subscribers

This is a discussion of uniform and non-uniform fields, magnetic, electric and otherwise, in which particles are propagated, and their relation to particle steppers.

The email gateway for this forum is: emfields-g4hn@slac.stanford.edu

Inline Depth:
 0 0
 1 1
 Add message: (add)

Question How turn on the magnetic field  by Paul <Paul>,   Nov 20, 06:18 
Hi Users,
Two weeks I try to add magnetic field to my simulation. I used to N03 and in this case it was verry simple. I wrote line /N03/det/setField 3 tesla in file.mac and gave control/execute file.mac in Idle> and it worked. Next I tried to do this same method in TestEm3. When I wrote /testem/det/setField 5 tesla in file.mac and started in Idle> I got the command: COMMAND NOT FOUND.
On the other hand iI tried to modify the codes in DetectorConstruction.cc file.
void DetectorConstruction::SetMagField(G4double fieldValue)
{
  //apply a global uniform magnetic field along Z axis
  //

fieldValue = 3. * tesla;
  G4FieldManager* fieldMgr
   = G4TransportationManager::GetTransportationManager()->GetFieldManager();

  if(magField) delete magField;		//delete the existing magn field

  if(fieldValue!=0.)			// create a new one if non nul
  { 
	  magField = new G4UniformMagField(G4ThreeVector(0.,fieldValue,0.));
    fieldMgr->SetDetectorField(magField);
    fieldMgr->CreateChordFinder(magField);
  } else {
    magField = 0;
    fieldMgr->SetDetectorField(magField);
  }
}

But the results from simulations without magnetic field are this same like 3 tesla and 20 tesla. Maybe I should to switch on magnetic filed but I don't know how. Why command from .mac file is don't working. Help me please Pawel

None Re: How turn on the magnetic field  by michel maire <michel maire>,   Nov 20, 07:01
Re: Question How turn on the magnetic field (Paul)
On Fri, 20 Nov 2009 14:18:18 GMT, Paul wrote: 
> Hi Users,
> Two weeks I try to add magnetic field to my simulation. I used to N03 and in this case it was verry simple. I wrote line /N03/det/setField 3 tesla in file.mac and gave control/execute file.mac in Idle> and it worked. Next I tried to do this same method in TestEm3. When I wrote /testem/det/setField 5 tesla in file.mac and started in Idle> I got the command: COMMAND NOT FOUND.


 I tried TestEm3 with the attached macro.
 The command setField works both in PreInit and Idle state

 Are you sure of the syntax of your commands ?

    Michel 

   Attachment:
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/11/20/06.59-82246-paul.txt

None =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc  by Hiromi Iinuma <Hiromi Iinuma>,   Nov 04, 23:32 
Hello,

I am checking spin direction as a function of time in electromagnetic 
field.
G4EqEMFieldWithSpin.cc includes Thomas-BMT equation as follows:
http://www-geant4.kek.jp/lxr/source//geometry/magneticfield/src/
G4EqEMFieldWithSpin.cc#L27
------snip---------
128    G4ThreeVector BField(Field[0],Field[1],Field[2]);
129    G4ThreeVector EField(Field[3],Field[4],Field[5]);
130 
131    G4ThreeVector u(y[3], y[4], y[5]);
132    u *= pModuleInverse;
133 
134    G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
135    G4double ucb = (anomaly+1./gamma)/beta;
136    G4double uce = anomaly + 1./(gamma+1.);
137 
138    G4ThreeVector Spin(y[9],y[10],y[11]);
139 
140    G4ThreeVector dSpin
141      = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.
cross(u))
142                                // from Jackson
143                                // -uce*Spin.cross(u.cross(EField)) );
144                                // but this form has one less 
operation
145                                - uce*(u*(Spin*EField) - EField*(Spin
*u)) );
146 
----snip--------

I think "EField" or "usb" should be divided by "c_light".
Relativistic charged particle feels electric field as magnetic field via
"\vec{B} =\vec{beta} x \vec{E} /light_c".
In the source code (G4EqEMFieldWithSpin.cc), dimension does not match.
"EField [volt/m]=[kg*m/(sec^3*A)]"
"BField[tesla]=[kg/(sec^2*A)]" 

Without applying 
"EField *=1/light_c " or "uce = uce/light_c",
simulated spin rotation frequency do not agree with the correct 
frequency from Thomas-BMT equation.

If you look at a different part of "G4EqEMFieldWithSpin.cc",
this problem may be clear.
-----snip-------
101    G4double cof2     = Energy/c_light ;
-----snip--------
117    dydx[3] = cof1*(cof2*Field[3] + (y[4]*Field[2] - y[5]*Field[1])) ;
118    
119    dydx[4] = cof1*(cof2*Field[4] + (y[5]*Field[0] - y[3]*Field[2])) ;
 
120  
121    dydx[5] = cof1*(cof2*Field[5] + (y[3]*Field[1] - y[4]*Field[0])) ;
------snip--------
In this part, Electric field; Field[3-5] is multiplied by 1/c_light.
Here,
128    G4ThreeVector BField(Field[0],Field[1],Field[2]);
129    G4ThreeVector EField(Field[3],Field[4],Field[5]);

Does anybody has such a problem?
Thanks in advance,
Hiromi
More Re: =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc  by Kevin Lynch <Kevin Lynch>,   Nov 05, 06:42
Re: None =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc (Hiromi Iinuma)
It's possible ... I asked the Geant collaboration to add this code, so I'm responsible for any bug is this piece of code. Before 4.9.2.p02, the whole beta x E terms was missing. I could have screwed it up. It's been a while since I wrote this, so let me trace back through what I did and see how the argument went....

When you say the "simulated spin rotation frequency do not agree with the correct frequency from Thomas-BMT equation", did you check that? Are you sure the C is not hiding somewhere else? Are you sure you're using the Maxwell equations in the right units? I cribbed most of the equations from Jackson, 2nd edition, and matched them onto the existing code to try and get the units right.

I'm not sure I got it all right, but it would be nice to do so before I need to start relying on my own simulation results :-)

None Re: =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc  by Hiromi Iinuma <Hiromi Iinuma>,   Nov 05, 08:06
Re: More Re: =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc (Kevin Lynch)
Dear Kevin,

Hello. Thank you for a reply.
Let me explain what I have done:
1) I checked so-called "spin-frozen" condition.
There is a good example article hep-ph/0012087v1.
This article discusses about how to cancel the muon g-2 precession by 
applying an electric field.

Here, I am going to use these parameters:
1) beta is muon velocity,
2) B is an uniform magnetic field and \vec{beta}.\vec{B}=0,
  "." means scalar product,
3) Radial electric field (E) in the lab frame and \vec{beta}.\vec{E}=0,
4) a=(g-2)/2 is muon anomalous magnetic moment.

The required electric field to cancel the g-2 precession is,
E=a*B*light_c*gamma^2*beta.

In case of gamma=5 and B=0.24 Tesla, the required electric field is E=2 
MV/m.
"Spin-frozen" happens when spin rotation cycle and muon rotation cycle 
are same.
In this case, both cycles should be 149.5 nsec. 
This is my expectation from the article hep-ph/0012087v1 and Jackson, 
3rd edition.
(Actually, my Jackson is Japanese version.)

Then I tried to apply proper electromagnetic field and muon momentum in 
GEANT4 simulation.
Although I got correct muon rotation cycle,I did not get correct spin 
rotation cycle.
Therefore I modified "G4EqEMFieldWithSpin.cc" as I wrote in the previous 
e-mail,
and finally I got expected results.
I can show you more detail analysis note as needed.

Best regards,
Hiromi



----- Original Message -----
> 
> *** Discussion title: Fields: Magnetic and Otherwise
> 
> It's possible ... I asked the Geant collaboration to add this code, so
> I'm responsible for any bug is this piece of code. Before 4.9.2.p02, 
the
> whole beta x E terms was missing. I could have screwed it up. It's 
been
> a while since I wrote this, so let me trace back through what I did 
and
> see how the argument went....
> 
> When you say the "simulated spin rotation frequency do not agree with
> the correct frequency from Thomas-BMT equation", did you check that? 
Are
> you sure the C is not hiding somewhere else? Are you sure you're using
> the Maxwell equations in the right units? I cribbed most of the
> equations from Jackson, 2nd edition, and matched them onto the 
existing
> code to try and get the units right.
> 
> I'm not sure I got it all right, but it would be nice to do so before 
I
> need to start relying on my own simulation results :-)
> 
> -------------------------------------------------------------
> Visit this GEANT4 at hypernews.slac.stanford.edu message (to reply or 
unsubscribe) at: 
> http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161/1.
html 
>
None Re: =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc  by Kevin Lynch <Kevin Lynch>,   Nov 06, 06:28
Re: None Re: =?UTF-8?B?QnVnPz8=?= A question of G4EqEMFieldWithSpin.cc (Hiromi Iinuma)
Just a followup for anyone that comes across this in the future: Hiromi, Peter Gumplinger, and I have talked about this off-forum, and we all agree with Hiromi that my code was faulty. A patch has been submitted to Peter that will hopefully appear in 9.3 to fix this. Just in case it doesn't for some reason, I attach the patch to this post. 

   Attachment:
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/11/06/06.27-11221-EMFieldWithSpin.patch.txt

Sad the proper time goes back  by Olivier Dadoun <Olivier Dadoun>,   Sep 10, 07:04
Hello,

I am trying to put an Electromagnetic field (in a localized volume) after an Tungsten target. I have followed the recommendation gave by Peter Gumplinger http://tinyurl.com/nfpk7t So inside my DetectorConstruction I have 

    G4ThreeVector positionTarget= G4ThreeVector(0,0,samplerLength+targetSizeZ/2);
    G4Box* solidTarget = new G4Box("Target",targetSizeXY/2.,targetSizeXY/2.,targetSizeZ/2.);
    G4LogicalVolume* logicTarget = new G4LogicalVolume(solidTarget, aMaterial,"Target");
    pTarget = new G4PVPlacement(0,positionTarget,logicTarget,"Target",lWorld,false,0);
    G4VisAttributes* LogVisAttTarget= new G4VisAttributes(G4Colour(.0,1.,0.));
    logicTarget ->SetVisAttributes(LogVisAttTarget);
    // register logical Volume in PolarizationManager with zero polarization
    G4PolarizationManager * polMgr = G4PolarizationManager::GetInstance();
    polMgr->SetVolumePolarization(logicTarget,G4ThreeVector(0.,0.,0.));

    //Target Calo
    PPSTargetCaloSD* TargetCaloSD= new PPSTargetCaloSD("TargetCaloSD",this);
    SDman->AddNewDetector(TargetCaloSD);
    logicTarget->SetSensitiveDetector(TargetCaloSD);

 static G4bool fieldIsInitialized = false;

    if(!fieldIsInitialized)

    {
        Field = new PPSField();
        pEquation = new G4EqMagElectricField(Field);
        pStepper = new G4ClassicalRK4 (pEquation);
        pFieldMgr = new G4FieldManager();
        //pFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
        pIntgrDriver = new G4MagInt_Driver(0.000001*mm,pStepper,pStepper->GetNumberOfVariables() );
        pChordFinder = new G4ChordFinder(pIntgrDriver);
        pFieldMgr->SetChordFinder( pChordFinder );
        pFieldMgr->SetDetectorField(Field);
        fieldIsInitialized = true;
    }

    G4ThreeVector positionCapture= G4ThreeVector(0,0,targetSizeZ+CaptureLength/2.);
    G4Box *solidCapture = new G4Box("CaptureBox",AMD_XY_size,AMD_XY_size,CaptureLength/2);
    G4LogicalVolume* logicCapture = new G4LogicalVolume(solidCapture,Vacuum,"theCaptureBox");
    pCapture=new G4PVPlacement(0,positionCapture,logicCapture,"CaptureBox",lWorld,false,0);
    G4VisAttributes* LogVisAttCapture= new G4VisAttributes(G4Colour(1.0,0.5,0.5));
    logicCapture ->SetVisAttributes(LogVisAttCapture);
    logicCapture->SetFieldManager(pFieldMgr,true);

And my PPSField class is only: 

PPSField::PPSField() 
{     }
void PPSField::GetFieldValue(const double point[4], double *Bfield ) const
{ 
	// Magnetic field
	Bfield[0] = 0;
	Bfield[1] = 0;
	Bfield[2] = 0.5*tesla;

	// Electric field
	Bfield[3] = 0;
	Bfield[4] = 0;
	Bfield[5] = 17.*megavolt/m;

}

When I run a electron beam of Several MeV impinging on my target I have those following error: 

G4ParticleChange::CheckIt    : the global time goes back  !!
  Difference:  0.074773488261472[ns] 
  G4ParticleChange::CheckIt    : the proper time goes back  !!
  Difference:  0.0026325540767444[ns] 
 G4ParticleChange::CheckIt 
      -----------------------------------------------
        G4ParticleChange Information  
      -----------------------------------------------
        # of 2ndaries       :                    0
      -----------------------------------------------
        Energy Deposit (MeV):                    0
        Non-ionizing Energy Deposit (MeV):                    0
        Track Status        :                Alive
        True Path Length (mm) :                 11.1
        Stepping Control      :                    0
        Mass (GeV)   :                    0
        Charge (eplus)   :                    0
        MagneticMoment   :                    0
                :  =                    0*[e hbar]/[2 m]
        Position - x (mm)   :                   -9
        Position - y (mm)   :                -2.96
        Position - z (mm)   :                 14.9
        Time (ns)           :                    0
        Proper Time (ns)    :             -0.00249
        Momentum Direct - x :               -0.795
        Momentum Direct - y :                -0.55
        Momentum Direct - z :                0.255
        Kinetic Energy (MeV):                   14
        Polarization - x    :                    0
        Polarization - y    :                    0
        Polarization - z    :                    0
        Touchable (pointer) :            0x4dbd0a0

Do you have an idea of what I am doing wrong ? by advance thank you Olivier

Ok Re: the proper time goes back  by Olivier Dadoun <Olivier Dadoun>,   Sep 10, 08:44
Re: Sad the proper time goes back (Olivier Dadoun)
pIntgrDriver = new G4MagInt_Driver(0.000001*mm,pStepper,pStepper->GetNumberOfVariables())

pStepper->GetNumberOfVariables() is equal to 6

If I change pIntgrDriver to

pIntgrDriver = new G4MagInt_Driver(0.000001*mm,pStepper,8)

it seems to fixe my problem (don't know why ).... olivier

Warning Re: the proper time goes back  by Gumplinger Peter <Gumplinger Peter>,   Sep 17, 19:45
Re: Ok Re: the proper time goes back (Olivier Dadoun)
Dear Oliver,

> I have followed the recommendation gave by

> Peter Gumplinger http://tinyurl.com/nfpk7t

Thanks for bringing this to my/users' attention. The code in the above forum link is still not entierly correct; it would be for a localized pure magnetic field but for EM-fields we need: 

   G4int nvar = 8;  // To integrate time & energy 
                    //    in addition to position, momentum

and the correct code for EM-fields should read: 

    G4int nvar = 8;
    fStepper = new G4ClassicalRK4( fEquation, nvar );

for the default in the G4ClassicalRK4 constructor is nvar=6 and hence tailored to pure mag. fields.

you can then subsequently use: 

    fIntgrDriver = new G4MagInt_Driver(fMinStep, 
                                     fStepper, 
                                     fStepper->GetNumberOfVariables() );

This should be the right approach rather than just: 

    fIntgrDriver = new G4MagInt_Driver(fMinStep, 
                                     fStepper, 
                                     8);

Please, also see the code fragment (Example 4.12) in the Application Developers Guide:

http://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch04s03.html#sect.EMField.Pract

Regards, Peter

Question A problem on the creation of electric field  by <lesward>,   23 May, 2009
In the examples provided in Geant4 which I've read, all the electric fields are created in the whole detector geometry. Therefore, I find it hard to describe a field for only a part of the detector. Is there anybody who can show me an example which sets electric fields for specific parts of the volume hierarchy? Thanks very much!

Feedback Re: A problem on the creation of electric field  by Gumplinger Peter <Gumplinger Peter>,   25 May, 2009
Re: Question A problem on the creation of electric field
Please, see in this forum. There are several postings dealing with 'local electric fields' and how to specify them. Among them are: 

http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/124.html
http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/114.html

Question Problems with Radial Electric Field  by Clarisse <Clarisse>,   27 Apr, 2009
Hi,

I designed a class derived from G4ElectricField to setup a non-uniform, radial electric field inside a cylindrical He3 Tube (E = Voltage / ( radius * ln(b/a)), where Voltage, b, and a are 3 parameters and r is the radius at which the field is evaluated and contains the position dependence).

I am attaching the relevant files (4 of them). In addition, I have the following few lines in myDetDetectorConstruction:

#include "myDetFieldSetup.hh"

logicGasTube = new G4LogicalVolume(solidGasTube , gasmixture, "gastube_log"); 

  // Set local field manager and local Electric Field in he3 tube
  G4bool allLocal = true ;

  logicGasTube->SetFieldManager( fEmFieldSetup->GetLocalFieldManager(), 
                                  allLocal ) ;

The code compiles fine, but I am getting the following error as I throw a thermal neutron at my tube: 

Start Run processing.
  G4ParticleChange::CheckIt    : the global time goes back  !!
  Difference:  11577.237025344[ns] 
  G4ParticleChange::CheckIt    : the proper time goes back  !!
  Difference:  11576.449246183[ns] 
 G4ParticleChange::CheckIt 
      -----------------------------------------------
        G4ParticleChange Information  
      -----------------------------------------------
        # of 2ndaries       :                    0
      -----------------------------------------------
        Energy Deposit (MeV):                    0
        Non-ionizing Energy Deposit (MeV):                    0
        Track Status        :                Alive
        True Path Length (mm) :                  0.8
        Stepping Control      :                    0
        Mass (GeV)   :                    0
        Charge (eplus)   :                    0
        MagneticMoment   :                    0
                :  =                    0*[e hbar]/[2 m]
        Position - x (mm)   :                 14.2
        Position - y (mm)   :                 20.8
        Position - z (mm)   :                -4.77
        Time (ns)           :                    0
        Proper Time (ns)    :            -1.16e+04
        Momentum Direct - x :                0.493
        Momentum Direct - y :               -0.863
        Momentum Direct - z :               -0.115
        Kinetic Energy (MeV):                0.191
        Polarization - x    :                    0
        Polarization - y    :                    0
        Polarization - z    :                    0
        Touchable (pointer) :            0xab1b008

*** G4Exception : 200
      issued by : G4ParticleChange::CheckIt
momentum, energy, and/or time was illegal
*** Event Must Be Aborted 
Run terminated.

I would appreciate your input in what I may be doing wrong here.

Thank you,

Clarisse 

   Attachment:
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/04/27/07.35-23670-nclude_myDetFieldSetup.hh
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/04/27/07.35-32084-DetRadialElectricField.hh
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/04/27/07.35-88341-3_src_myDetFieldSetup.cxx
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/04/27/07.35-98902-etRadialElectricField.cxx

Question The Strangest Problem I've ever met  by Zukai Wang <Zukai Wang>,   19 Feb, 2009 
Dear All:
    Actually, I really have no idea which category my problem belongs to.
    If my magnetic field definition is in a simple way like this:
------------------------------------------------------------------------- 
void ExN02MagneticField::GetFieldValue(const double Point[3],double *Bfield) const
{
     Bfield[0]=0.;
     Bfield[1]=0.;
     Bfield[2]=0.;
if(fabs(Point[0])<75.*mm && fabs(Point[1])<405.*mm && fabs(Point[2])<705.*mm)
     Bfield[0]=-6450.*gauss;}
--------------------------------------------------------------------------
    The whole project works well.

    But when I put on a real magnetic field like this:

------------------------------------------------------------------
   void ExN02MagneticField::GetFieldValue(const double Point[3],double *Bfield) const
{
     Bfield[0]=0.;
     Bfield[1]=0.;
     Bfield[2]=0.;

 fstream in;

 G4double x=0.,y=0.,z=0.,Bx=0.,By=0.,Bz=0.;

 if (fabs(fabs(Point[0])-70.*mm)<5.*mm)
   in.open("Bmap36.dat");

  else if (fabs(fabs(Point[0])-60.*mm)<5.*mm)
   in.open("Bmap37.dat");

 else if (fabs(fabs(Point[0])-50.*mm)<5.*mm)
   in.open("Bmap38.dat");  

 else  if (fabs(fabs(Point[0])-40.*mm)<5.*mm)
   in.open("Bmap39.dat");

 else  if (fabs(fabs(Point[0])-30.*mm)<5.*mm)

 {  if (fabs(Point[1])<75.*mm)
    in.open("Bmap401.dat");

 else if (fabs(Point[1])<145.*mm)
   in.open("Bmap4022.dat");

 else in.open("Bmap4021.dat");
  }

  else  if (fabs(fabs(Point[0])-20.*mm)<5.*mm)

 {if (fabs(Point[1])<65.*mm)
    in.open("Bmap411.dat");

  else in.open("Bmap412.dat");}

 else  if (fabs(fabs(Point[0])-10.*mm)<5.*mm)
   in.open("Bmap42.dat");

 else  if (fabs(Point[0])<5.*mm)
   in.open("Bmap43.dat");

  while ( fabs(Point[1])<405.*mm && fabs(Point[2])<705.*mm) 

{
     in>>x>>y>>z>>Bx>>By>>Bz;

 if (!in.good()) break;

 if(Point[0]>0 && fabs(Point[1]-y*mm)<5.*mm && fabs(Point[2]-z*mm)<5.*mm)

{    Bfield[0]=Bx*gauss;
     Bfield[1]=By*gauss;
     Bfield[2]=Bz*gauss;

     break;
    }

else if (Point[0]<0 && fabs(Point[1]-y*mm)<5.*mm && fabs(Point[2]-z*mm)<5.*mm) 

{     Bfield[0]=Bx*gauss;
     Bfield[1]=-By*gauss;
     Bfield[2]=-Bz*gauss;

     break;
       }

else if (Point[0]>0 && fabs(-Point[1]-y*mm)<5.*mm && fabs(Point[2]-z*mm)<5.*mm) 

{     Bfield[0]=Bx*gauss;
     Bfield[1]=-By*gauss;
     Bfield[2]=Bz*gauss;

     break;
       }

else if (Point[0]<0 && fabs(-Point[1]-y*mm)<5.*mm && fabs(Point[2]-z*mm)<5.*mm) 

{     Bfield[0]=Bx*gauss;
     Bfield[1]=By*gauss;
     Bfield[2]=-Bz*gauss;

     break;
       }

else if (Point[0]>0 && fabs(Point[1]-y*mm)<5.*mm && fabs(-Point[2]-z*mm)<5.*mm) 

{    Bfield[0]=Bx*gauss;
     Bfield[1]=By*gauss;
     Bfield[2]=-Bz*gauss;

     break;
       }

   else if (Point[0]<0 && fabs(Point[1]-y*mm)<5.*mm && fabs(-Point[2]-z*mm)<5.*mm)

{     Bfield[0]=Bx*gauss;
     Bfield[1]=-By*gauss;
     Bfield[2]=Bz*gauss;

     break;
       }

  else if (Point[0]>0 && fabs(-Point[1]-y*mm)<5.*mm && fabs(-Point[2]-z*mm)<5.*mm)

{     Bfield[0]=Bx*gauss;
     Bfield[1]=-By*gauss;
     Bfield[2]=-Bz*gauss;

     break;
       }

  else if (Point[0]<0*mm && fabs(-Point[1]-y*mm)<5.*mm && fabs(-Point[2]-z*mm)<5.*mm)

{     Bfield[0]=Bx*gauss;
     Bfield[1]=By*gauss;
     Bfield[2]=Bz*gauss;

     break;
       }
 }

if (in.good())
   in.close();

} -------------------------------------------------------------------------- 

    The magnetic field values must be read from several files. The magnetic field values of more than 17000*8 points should be defined. I have edited the data files in a certain way to make the running of the program as fast as possible. The codes compiled and nothing wrong has been found.
    I shoot Krypton 78 in the particle gun and it will hit a detector which can record its kinetic energy and other information. The information has been showed on the screen. And at the end of an event, the event id should show up on the screen according to my codes. However, the event id has never appeared on the screen. Never one event has been finished!
    I wrote something in SteppingAction to track the Krypton 78 and found it pass through the magnetic field quickly (only 3 or 4 min) and hit the detector.
    And then, this ion kept on going. It reached the end of World and since then, the program ceased running. I mean, it was not aborted, it just stopped as if the poor ion did not know what to do next!
    I then increased the size of World. The ion (Kr78) lost all its energy in the last step, and the program stopped again! 
    No errors and no segment faults have been reported. If I didn't track the ion from the particle gun, I would think that the whole program might last for hours!
    So, I doubt that the pointers (or something else) has exceed the maximum of Geant4.8.3(my version). If you have come across similar problems, I'd like to thank you in advance for your direction.
                                                         Zukai Wang

Note Re: The Strangest Problem I've ever met  by Tom Roberts <Tom Roberts>,   19 Feb, 2009
Re: Question The Strangest Problem I've ever met (Zukai Wang)
You most definitely do not want to open and read a file during EVERY call to GetFieldValue() -- this will slow your simulation down by a truly enormous factor! You consider "only 3 or 4 min" to be "quickly" -- I routinely track hundreds of particles per second for many meters through magnetic field maps. Even a single sprintf() or sscanf() during tracking slows it down unacceptably.

Arrange the code so you open and read the files ONCE during initialization and keep the results around to use in GetFieldValue(). It is generally better (more general, easier to modify) to use a single file with a field map, rather than multiple files. A good design has the field map file describe its own geometrical structure (boundaries of the map, grid within the map, units, etc.) and the code uses this rather than hard-coded values for the geometry. I strongly recommend an ASCII format, so you can edit it easily, and can write conversion programs from the output of an EM modeling program (superfish, tosca, ...) to the format your simulation reads. Field maps of several tens of megabytes are not a problem on modern computers.

Whenever a track reaches the edge of the world, the track is killed. Similarly, when a track stops (kinetic energy goes to 0 or below a limit set in the program) it is killed unless there is an applicable physics process that makes it get captured or decay (which also kills the track). Just think about it -- how could a track possibly be tracked outside the known world? Or tracked forever while stopped?

Normally, when all tracks have been processed (i.e. tracked until they are killed), the run is ended and the program exits. But you must do this in your main program, as in the examples.

Question Re: The Strangest Problem I've ever met  by Zukai Wang <Zukai Wang>,   20 Feb, 2009
Re: Note Re: The Strangest Problem I've ever met (Tom Roberts)
On Thu, 19 Feb 2009 15:24:38 GMT, Tom Roberts wrote:
> You most definitely do not want to open and read a file during EVERY
> call to GetFieldValue() -- this will slow your simulation down by a
> truly enormous factor! You consider "only 3 or 4 min" to be "quickly" --
> I routinely track hundreds of particles per second for many meters
> through magnetic field maps. Even a single sprintf() or sscanf() during
> tracking slows it down unacceptably.
> 
> Arrange the code so you open and read the files ONCE during
> initialization and keep the results around to use in GetFieldValue(). It
> is generally better (more general, easier to modify) to use a single
> file with a field map, rather than multiple files. A good design has the
> field map file describe its own geometrical structure (boundaries of the
> map, grid within the map, units, etc.) and the code uses this rather
> than hard-coded values for the geometry. I strongly recommend an ASCII
> format, so you can edit it easily, and can write conversion programs
> from the output of an EM modeling program (superfish, tosca, ...) to the
> format your simulation reads. Field maps of several tens of megabytes
> are not a problem on modern computers.
> 
> Whenever a track reaches the edge of the world, the track is killed.
> Similarly, when a track stops (kinetic energy goes to 0 or below a limit
> set in the program) it is killed unless there is an applicable physics
> process that makes it get captured or decay (which also kills the
> track). Just think about it -- how could a track possibly be tracked
> outside the known world? Or tracked forever while stopped?
> 
> Normally, when all tracks have been processed (i.e. tracked until they
> are killed), the run is ended and the program exits. But you must do
> this in your main program, as in the examples.
> 
Thank you so much for the direction. I have modified my code this way:
------------------------------------------------------------------------
ExN02MagneticField::ExN02MagneticField()
{
  fstream in;
  G4int n=0;
 in.open("Bmap.dat");

  while (in.good()&&!in.eof()) {
  in>>  x[n]>>y[n]>>z[n]>>Bx[n]>>By[n]>>Bz[n];

 G4cout<<  x[n]<<  y[n]<<  z[n]<<  Bx[n]<<  By[n]<<  Bz[n]<< G4endl;
 n++;}
  in.close();
  GetGlobalFieldManager()->SetDetectorField(this);
  GetGlobalFieldManager()->CreateChordFinder(this);

 }
ExN02MagneticField::~ExN02MagneticField()
{;}

void ExN02MagneticField::GetFieldValue(const double Point[3],double *Bfield) const
{
     Bfield[0]=0.;
     Bfield[1]=0.;
     Bfield[2]=0.;
    //G4int number=0;
// if (fabs(Point[0])<75.*mm && fabs(Point[1])<405.*mm && fabs(Point[2])<705.*mm)
  // Bfield[0]=-6450.*gauss;
    for (G4int i=0;i<23288;i++) 
{
   if (fabs(Point[0])>75.*mm || fabs(Point[1])>405.*mm || fabs(Point[2])>705.*mm) break;
   if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
{   

    Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
    Bfield[2]=Bz[i]*gauss;
    break;}
   else if (-fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
    Bfield[2]=-Bz[i]*gauss;
    break;}
    else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
    Bfield[2]=Bz[i]*gauss;
    break;}
    else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}
      else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}

     else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
     Bfield[2]=Bz[i]*gauss;
    break;}
    else if (fabs(Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=-By[i]*gauss;
     Bfield[2]=-Bz[i]*gauss;
    break;}
    else if (fabs(-Point[0]-x[i]*mm)<5.*mm && fabs(-Point[1]-y[i]*mm)<5.*mm && fabs(-Point[2]-z[i]*mm)<5.*mm)
  { Bfield[0]=Bx[i]*gauss;
    Bfield[1]=By[i]*gauss;
     Bfield[2]=Bz[i]*gauss;
    break;}
 }
}
-------------------------------------------------------------------------
   Now, I have arranged the code to open and read the files once during initialization and keep the results around to use in GetFieldValue(). The program runs much faster.
   The funny thing is: IT GOT STUCK AGAIN!
   Just like the situation before I made this improvement. I'd like to mention that the main program is fine. Because if I change the magnetic field in a simple way, the first event ends quickly and the second comes up.
   But now, there never has been an end for my first event!
   I really want to know what the matter is.

Question Is This the Reason?  by Zukai Wang <Zukai Wang>,   20 Feb, 2009
Re: Question Re: The Strangest Problem I've ever met (Zukai Wang)
    After modifying the code of magnetic field, the simulation got a lot faster. The reason why it still takes a lot of time after the ion hit the detector might due to the great many secondary particles generated. The magnetic field is also applied to them. That is perhaps why my simulation got a lot faster when a simple magnetic field was applied.
    So far, I only need the kinetic energy, the position of the Kr78 and the global time (about 47ns) when it hits the detector. Can I just kill all the secondary particles without affecting the information I need? Or I may end one event as soon as Kr 78 reaches the end of the world.
    For further work, I will need to know the total energy deposition of every chamber. Do I have any approaches to make it faster besides improving the method of getting the magnetic field value?
Question Spin tracking of muons at rest  by Kevin Lynch <Kevin Lynch>,   18 Feb, 2009
As a follow-up to my previous question on muon spin tracking of moving muons (post 155, http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/155.html), I have a new question about hybrid systems where moving muons come to rest in a target. To help myself and others in future, I'm going to start with a few observations that aren't handled in much detail in the documentation, and then get to my question:

1) To decay a muon, you need to instantiate a "Decay" process, and Geant4 provides 2: G4Decay, and G4DecayWithSpin. The Decay processes can be attached to both moving and resting particles by appropriate calls to G4ProcessManager::SetProcessOrdering, with the ProcessManager of the given particle type.

2) These processes permit decay, but somewhat confusingly (at least to me!) have nothing to do with whether or not the daughter products are produced in correlation with the spin; that determination is made by the selection of DecayChannels that are attached to the G4Muon{Plus,Minus} class. Geant4 provides a number of choices, depending on your needs: G4MuonDecayChannel, G4MuonDecayChannelWithSpin, and G4MuonRadiativeDecayChannelWithSpin.

3) This (well, and the code) leads me to conclude that the only real difference between the two Decay processes is whether the muon spin is evolved. The G4Decay class provides no evolution, while the G4DecayWithSpin class provides the static BMT \vec{s} \cross \vec{B} component. My experiments confirm that this does something to the spin of at rest muons, although I haven't tested the precession in detail.

4) For muons in motion, the spin is evolved differently, by the G4Transportation processes. In particular, the partial (full) BMT spin evolution equations given in the G4Mag_SpinEqRhs (G4EqEMFieldWithSpin) class.

I hope my understanding is correct. Assuming it is, here's my question:

For muons at rest, I clearly need to attach the G4DecayWithSpin class, via a call equivalent to G4MuonPlus::MuonPlus()->SetProcessOrdering(new G4DecayWithSpin(), idxAtRest), to get at rest spin evolution. But if I attach G4DecayWithSpin to idxPostStep, do I simultaneously get BOTH G4Mag_EqRhs BMT evolution, AND the G4DecayWithSpin evolution? In other words, do I double count the spin rotation for muons in motion? Should I attach G4Decay to idxPostStep instead of G4DecayWithSpin? The only example in the 4.9.2 distribution that uses G4DecayWithSpin is extended/field04, which attaches G4DecayWithSpin in both cases. Assuming I'm wrong, can you point out how the code in the distribution ensures the G4DecayWithSpin terms don't come in to play for moving muons?

Setting up a test environment to figure this out in detail is somewhat onerous, hence my hope that someone can quickly set me straight. Thanks in advance for you comments!

Agree Re: Spin tracking of muons at rest  by Gumplinger Peter <Gumplinger Peter>,   18 Feb, 2009
Re: Question Spin tracking of muons at rest (Kevin Lynch)
Thanks Kevin, for both of your two insightful postings.

1) yes

2) The G4Decay class itself (G4DecayWithSpin inherits from G4Decay) does not define decay modes of the particle. Geant4 provides two ways of doing this: using G4DecayChannel in G4DecayTable, and using thePreAssignedDecayProducts of G4DynamicParticle. The G4Decay class calculates the PhysicalInteractionLength and boosts decay products created by G4VDecayChannel or event generators.

3) Yes, G4DecayWithSpin retrieves or sets a particles polarization and evolves it through the remainder of a particle's lifetime. (this is crucial to your ultimate question).

4) yes

> But if I attach G4DecayWithSpin to idxPostStep, do I simultaneously get BOTH

> G4Mag_EqRhs BMT evolution, AND the G4DecayWithSpin evolution?

The method G4DecayWithSpin::Spin_Precession takes the argument fRemainderLifeTime. For decays in flight, this will be zero, and consequently, the calculated remaining spin precession angle will also be zero.

> Should I attach G4Decay to idxPostStep instead of G4DecayWithSpin?

This will not work because only G4DecayWithSpin has the code:

parent_polarization = aParticle->GetPolarization(); decaychannel->SetPolarization(parent_polarization);

necessary to set the spin direction in the DecayChannelWithSpins.

> can you point out how the code in the distribution ensures

> the G4DecayWithSpin terms don't come in to play for moving muons?

(see above)

rotationangle = fRemainderLifeTime * omega = 0. * omega

You can verify that fRemainderLifeTime has already been updated (reduced to zero) by the time the G4DecayWithSpin::DecayIt is called for a decay in flight.

Peter

None Re: Spin tracking of muons at rest  by Kevin Lynch <Kevin Lynch>,   18 Feb, 2009
Re: Agree Re: Spin tracking of muons at rest (Gumplinger Peter)
On Thu, 19 Feb 2009 00:30:27 GMT, Gumplinger Peter wrote: 
> > can you point out how the code in the distribution ensures
> 
> > the G4DecayWithSpin terms don't come in to play for moving muons?
> 
> (see above)
> 
> rotationangle = fRemainderLifeTime * omega = 0. * omega
> 
> You can verify that fRemainderLifeTime has already been updated (reduced
> to zero) by the time the G4DecayWithSpin::DecayIt is called for a decay
> in flight.
> 
> Peter
>

This is the key point I was missing; now that I understand it, the rest of the code makes perfect sense. Thanks again for your help.

Question Spin tracking of moving muons  by Kevin Lynch <Kevin Lynch>,   12 Feb, 2009
I have a number of questions about spin tracking of moving particles (muons particularly).

To set up the questions, a few words about our experimental environment: We have a polarized (of course:-) muon beam travelling in vacuum through a set of magnetic fields, and coming to rest in a target, also immersed in a magnetic field. We need to track the muon spin, and apply spin dependent muon decay. During transport, we use the G4Mag_SpinEqRhs equations, and see the spin precessing correctly (we think). We use G4DecayWithSpin, and see the proper distribution of muon decays around the muon spin direction. I'm currently using both 4.9.1 and 4.9.2 on Linux machines. So far so good.

Here are some of the issues that I don't quite understand:

1) If I understand correctly, the units of the evolution equation have d\vec{s}/dx = (1/\beta)(d\vec{s}/dt). The dominant term on the RHS in the BMT spin evolution for low energy muons is the \vec{s} \cross \vec{B} term, which has no problem as \beta->0. But, because we multiply by a number that diverges (1/\beta) in the stopping limit, we in principle could have a big problem with inaccurate evolution. Practically, given the way Geant4 handles the energy loss, should I be worried about this? Do we ever end up in a condition where 1/\beta gets big enough to be a problem? If so, is there any good way to deal with this?

2) In some of the EqRhs equation classes, the EvaluatRhsGivenB term includes the statements: 

   dydx[6] = dydx[8] = 0.;//not used
   // Lab Time of flight
   dydx[7] = inverse_velocity;

for example, G4EqEMFieldWithSpin. But others, particularly the G4Mag_SpinEqRhs that I need, has instead 

   // Initialise the values of dydx that we do not update.
   dydx[6] = dydx[7] = dydx[8] = 0.0;

Is the difference meaningful? I definitely need the lab TOF to have the right value.

3) In another experiment, we have higher momentum muons and also electric fields. All else being correct, I should be able to use G4EqEMFieldWithSpin to evolve the spin in these fields. But I think there's a problem with the implementation in G4EqEMFieldWithSpin::EvaluateRhsGivenB: the BMT equation has a term of the form \vec{s} \cross (\vec{\beta} \cross \vec{E}), but the implementation of the spin evolution in this class is identical to that of G4Mag_SpinEqRhs. Is this intentional (am I missing something), or an oversight (ie. bug)?

Feedback Re: Spin tracking of moving muons  by Gumplinger Peter <Gumplinger Peter>,   18 Feb, 2009
Re: Question Spin tracking of moving muons (Kevin Lynch)
Hi Kevin,

Thanks for your detailed and constructive forum posting.

3) As for your last question, the missing third term is an oversight. I noticed it a little while ago when I went back to Jackson. I made mental note to investigate and fix it. Actually, I have no feeling as to when this term (size of E field and beta) becomes sizable relative to the others. Until now, I only thought of E fields inside drift chambers. According to Jackson: "If the particle is relativistic (beta->1), even the presence of an electric field causes the longitudinal polarization to change only very slowly, at a rate proportional to gamma**(-2) times the electric field component perpendicular to v_vector." (I assume you are dealing with 'surface muons' - longitudinal polarized muons?)

2) G4EqEMFieldWithSpin copies from G4EqMagElectricField while G4Mag_SpinEqRhs is a copy of G4Mag_UsualEqRhs. Only if the field can change energy, then the time must be integrated, else:

http://www-geant4.kek.jp/lxr/source/processes/transportation/src/G4Transportation.cc#L480

the time is incremented from the particle's velocity and step length.

1) > Practically, given the way Geant4 handles the energy loss, should I be

> worried about this?

> Do we ever end up in a condition where 1/\beta gets big enough to be a

> problem? If so, is there any good way to deal with this?

You have a very good point. I don't know the answers.... I suggest you flag large final step spin variations - should they occur - in your application. My experience has only ever dealt with muons stopping in a rather homogeneous, albeit strong magnetic field. We never observed a weired spin. We were very sensitive to any false prediction (depolarization) at the 10-3 level and consequently plotted the final spin predicition for each event history for many millions of such histories.

A group who extensively uses spin tracking in G4 is from PSI: 

A general paper:
              T. Shiroka et al.,
              "Geant4 as a simulation framework in muSR"
              proceedings of the muSR 2008 conference,  Tsukuba, Japan
              (21-25th July 2008), to be published in Physica B.

None Re: Spin tracking of moving muons  by Kevin Lynch <Kevin Lynch>,   18 Feb, 2009
Re: Feedback Re: Spin tracking of moving muons (Gumplinger Peter)
On Thu, 19 Feb 2009 01:29:47 GMT, Gumplinger Peter wrote: 
> Hi Kevin,
> 
> Thanks for your detailed and constructive forum posting.
>

And thanks for your detailed answers. I've now got everything I need to make forward progress. 

> According to Jackson: "If the particle is relativistic (beta->1), even
> the presence of an electric field causes the longitudinal polarization
> to change only very slowly, at a rate proportional to gamma**(-2) times
> the electric field component perpendicular to v_vector." (I assume you
> are dealing with 'surface muons' - longitudinal polarized muons?)
>

Yes and no. I'm working with a student on a project with surface muons, with no E fields, but where stopped muon spin evolution is important. But I'm also working on a next generation g-2 experiment, with large quadrupole E fields that rotate spins, requiring a non-negligible correction to the final experimental results. And, I've just been asked to consider doing some spin tracking for a potential muon EDM experiment, with very large E fields. So, all these issues are very interesting to me.

I have been able to implement the missing terms in the G4EqEMFieldWithSpin class, and will send them to you and post them here once I've tested them. I also plan to implement the EDM field equations; although I doubt anyone else will have a use for them, I plan on offering the implementation for inclusion in the Geant4 distribution if you would be interested.

Ok Re: Spin tracking of moving muons  by Gumplinger Peter <Gumplinger Peter>,   19 Feb, 2009
Re: None Re: Spin tracking of moving muons (Kevin Lynch)
Hi Kevin,

>I have been able to implement the missing terms in the G4EqEMFieldWithSpin

>class, and will send them to you and post them here once I've tested them. I

>also plan to implement the EDM field equations; although I doubt anyone else

>will have a use for them, I plan on offering the implementation for inclusion in

>the Geant4 distribution if you would be interested.

Thanks for offering your assistance in improving Geant4. - Looks like there is great interest in the G4 community and this once obscure part of G4 is getting some attention. I am happy to see that!

Much appreciated (and anticipated) Peter

None Implementation of missing terms and EDM spin evolution  by Kevin Lynch <Kevin Lynch>,   25 Feb, 2009
Re: Ok Re: Spin tracking of moving muons (Gumplinger Peter)
As promised, I attach:

1) G4EqEMFieldWithSpin.cc: a revised version of the Geant 4.9.2 distribution's BMT spin evolution equation in the presence of both B and E fields. This includes the previously missing \vec{\beta} \cross \vec{E} term, as well as some extensive comments on the various terms.

2) G4EqEMFieldWithEDM.{hh,cc}: My own little contribution to the Geant4 distribution. Based on the G4EqEMFieldWithSpin class, this class implements the full BMT spin evolution equation as well as the associated EDM spin evolution equation. Some references are listed, as well as the implementation. The default value of the EDM is set to 0, and needs to be set by the user. These files need correct CVS ID Tags inserted at the noted places if they are accepted. 

   Attachment:
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/02/25/11.38-58609-G4EqEMFieldWithSpin.cc
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/02/25/11.38-81748-G4EqEMFieldWithEDM.hh
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2009/02/25/11.38-16941-G4EqEMFieldWithEDM.cc

None Re: Spin tracking of moving muons  by Tom Roberts <Tom Roberts>,   18 Feb, 2009
Re: Feedback Re: Spin tracking of moving muons (Gumplinger Peter)
> 3) As for your last question, the missing third term is an oversight. I
> noticed it a little while ago when I went back to Jackson. I made mental
> note to investigate and fix it. Actually, I have no feeling as to when
> this term (size of E field and beta) becomes sizable relative to the
> others. Until now, I only thought of E fields inside drift chambers.
> According to Jackson: "If the particle is relativistic (beta->1), even
> the presence of an electric field causes the longitudinal polarization
> to change only very slowly, at a rate proportional to gamma**(-2) times
> the electric field component perpendicular to v_vector." (I assume you
> are dealing with 'surface muons' - longitudinal polarized muons?)

In a few months I will be interested in simulating polarized muons in a muon cooling channel. This involves RF cavities with electric fields of 10-30 Mega-Volts/meter, perhaps more. There are also magnetic fields up to 30 Tesla, perhaps more. The muon momentum starts 200-600 MeV/c or so, and depending on application remains above ~100 MeV/c or decreases all the way down to zero. So I hope you can implement this without approximation, including terms which might be negligible in some physical situations.

The highly-polarized muons come from forward or backward decays of pions, so they are essentially longitudinally polarized.

Question recommendation for atmosphere mag field setup   by Scott Nutter <Scott Nutter>,   01 Feb, 2009 
I am working on a simulation of the production of synchrotron radiation by very high energy electrons traversing the Earth's magnetic field.  I have to have times of synchrotron photons arriving at a certain height in the atmosphere good to better than a nanosecond. I have been having trouble getting this kind of accuracy.  I base this conclusion on an algorithm to reconstruct the zenith angle looking at the time difference of two spatially nearby photons at a certain height, which arrive on essentially parallel trajectories.  The resulting calculated zenith angles indicate the photons are all arriving from nearly vertically above. This algorithm works fine in a simple monte carlo using geant3 RK tracking and simple straight projection of photon tracks from production positions returned by the g3synge generator.

In reading about magnetic field setup in the application users manual, I noticed a small reference to "integrating over time" in the section on stepper choice. (The code fragment has "nvar=8" and subsequent comments in the version for G4.9.1, the section is 4.3.2.4.)  Since the atmosphere is a completely different scale than the typical HEP experiment, I think it best if I get some recommendations on setting up a field properly, so that times of synchrotron photon production are very carefully computed.  

The field characteristics are:
o The purely magnetic field B~0.5 gauss strength, typical Bperp ~0.1 gauss.  The field varies very little (<25% over 360 km). 
o Atmosphere geometry is G4Orbs setup like "Russian dolls" with radii the size of the Earth + different heights.  The resulting layer thicknesses range from several meters to 50 km.  The thickness of each atmosphere layer is determined by requiring it to have a column density of 20 mg/cm2. Thus the G4Orbs are filled with air of different densities. The lowest density of the air is 1e-22 g/cm3, the lowest limit I could get G4.9.2 to accept. I track from 400 km height above the Earth's surface to 38 km height.
o Electron energies are from 2-50 TeV. 
o My test case is a straightdown electron. With all physics off, the deflection for a constant perpendicular field is in agreement with calculations, and is ~19.5 m for Bperp=0.1 gauss over 362 km.
o My current code for defining the field:
        myField = new CrestMagField(latitude, longitude, time);
        G4FieldManager* fieldMgr
            = G4TransportationManager::GetTransportationManager()
              ->GetFieldManager();
        fieldMgr->SetDetectorField(myField);
        fieldMgr->CreateChordFinder(myField);

Any recommendations? There is much discussion of special cases and fine tuning of the field.  I believe the default tracker is tracking correctly now, but I do not believe the times yet.
Thanks.
Sad electron stuck on volume boundary with (local) uniform magnetic field  by Brian Fisher <Brian Fisher>,   21 Jul, 2008 
I'm experiencing a situation similar to posting #96 in this forum. That thread seems to have gone dead, so I'm posting a new thread.

I'm using Geant4 9.1p02

In a very simple setup, electrons sometimes are "bouncing" between a volume (with a local magnetic field) and its mother (which does not have a magnetic field).  Both volumes are filled with "vacuum."

I get many (tens of thousands) of messages like these:

WARNING - G4Navigator::ComputeStep()
          Track stuck, not moving for 10 steps
          in volume -expHall- at point (-8.257550732,-9.52505,18.49195399)
          direction: (-0.01057780492,7.285534641e-05,0.9999440508).
          Potential geometry or navigation problem !
          Trying pushing it of 9e-10 mm ...

I've been modeling charged particle detector telescopes, and in the process of hunting this problem down, I have reduced the geometry of my problem to:

1) a 1x1x1 m^3 box filled with vacuum. This is the "expHall" volume
2) a "BFieldZone" box also filled with vacuum with the following parameters:

 Solid type: G4Box
 Parameters:  half length X: 1.39526 cm 
              half length Y: 0.952505 cm 
              half length Z: 0.952505 cm 
-----------------------------------------------------------
 Center of Box: X = 0 cm Y = 0 cm Z = 1.98525 cm 

3) the "BFieldZone" box has a uniform 2 kilogauss magnetic field in the +X direction.

This is how I've defined the (local) magnetic field for the "BFieldZone" volume:

	G4double = fFieldValue = 2.0*kilogauss;
	fMagField = new G4UniformMagField(G4ThreeVector(fFieldValue, 0.0, 0.0));
	fEquation = new G4Mag_UsualEqRhs(fMagField);
	fStepper = new G4ExactHelixStepper(fEquation);
	fIntgrDriver = new G4MagInt_Driver(1.0*micrometer, fStepper, fStepper->GetNumberOfVariables(), 2);
	fChordFinder = new G4ChordFinder(fIntgrDriver);
	fFieldMgr = new G4FieldManager(fMagField, fChordFinder);

	logicBFieldZone->SetFieldManager(fFieldMgr, true);

4) The "expHall" box has no magnetic field.


When I turn the tracking verbosity up to 2, the 'offending' event looks like this:


*********************************************************************************************************
* G4Track Information:   Particle = e-,   Track ID = 1,   Parent ID = 0
*********************************************************************************************************

Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0  -8.06 mm  -9.69 mm 7.87e-55 fm   1.97 MeV     0 eV      0 fm      0 fm      expHall    initStep
    1  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV4.82e-19 eV   1.81 cm   1.81 cm      expHall  Transportation
    2  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm   physiBFieldZone  Transportation
    3  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm      expHall  Transportation
    4  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm   physiBFieldZone  Transportation
    5  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm      expHall  Transportation
    6  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm   physiBFieldZone  Transportation
    7  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm      expHall  Transportation
    8  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm   physiBFieldZone  Transportation
    9  -8.25 mm  -9.53 mm   1.81 cm   1.97 MeV     0 eV      0 fm   1.81 cm      expHall  Transportation
   10  -8.25 mm  -9.52 mm   1.82 cm   1.97 MeV2.39e-21 eV   89.9 um   1.82 cm   physiBFieldZone  Transportation
   11  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV1.91e-22 eV   7.16 um   1.82 cm   physiBFieldZone  Transportation
   12  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm      expHall  Transportation
   13  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm   physiBFieldZone  Transportation
   14  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm      expHall  Transportation
   15  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm   physiBFieldZone  Transportation
   16  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm      expHall  Transportation
   17  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm   physiBFieldZone  Transportation
   18  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm      expHall  Transportation
   19  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm   physiBFieldZone  Transportation
   20  -8.25 mm  -9.53 mm   1.82 cm   1.97 MeV     0 eV      0 fm   1.82 cm      expHall  Transportation
   21  -8.26 mm  -9.52 mm   1.83 cm   1.97 MeV2.39e-21 eV   89.9 um   1.83 cm   physiBFieldZone  Transportation
   22  -8.26 mm  -9.53 mm   1.83 cm   1.97 MeV1.32e-22 eV   4.96 um   1.83 cm   physiBFieldZone  Transportation
   23  -8.26 mm  -9.53 mm   1.83 cm   1.97 MeV     0 eV      0 fm   1.83 cm      expHall  Transportation
   24  -8.26 mm  -9.53 mm   1.83 cm   1.97 MeV     0 eV      0 fm   1.83 cm   physiBFieldZone  Transportation
   25  -8.26 mm  -9.53 mm   1.83 cm   1.97 MeV     0 eV      0 fm   1.83 cm      expHall  Transportation
   26  -8.26 mm  -9.53 mm   1.83 cm   1.97 MeV     0 eV      0 fm   1.83 cm   physiBFieldZone  Transportation

[and so on, seemingly forever]


If the verbosity is turned up to 3, I get this:

[snip]

#Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
30132  -8.26 mm  -9.53 mm   1.85 cm   1.97 MeV     0 eV      0 fm   1.85 cm   physiBFieldZone  Transportation

 >>AlongStepDoIt (after all invocations):
    ++List of invoked processes 
      1) Transportation
      2) msc
      3) eIoni
      4) eBrem

    ++G4Step Information 
      Address of G4Track    : 0x3046970
      Step Length (mm)      : 0
      Energy Deposit (MeV)  : 0
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -8.257551022101929  -8.257551022101929
         Position - y (mm)   :             -9.52505            -9.52505
         Position - z (mm)   :    18.49198141201282   18.49198141201282
         Global Time (ns)    :  0.06304591361302052 0.06304591361302052
         Local Time (ns)     :  0.06304591361302052 0.06304591361302052
         Proper Time (ns)    :  0.01300913309755812 0.01300913309755812
         Momentum Direct - x : -0.01057780492354854-0.01057780492354854
         Momentum Direct - y : 7.285445476253925e-057.285445476253925e-05
         Momentum Direct - z :   0.9999440508024575  0.9999440508024575
         Momentum - x (MeV/c): -0.02563161644372119-0.02563161644372119
         Momentum - y (MeV/c): 0.00017653733021041490.0001765373302104149
         Momentum - z (MeV/c):    2.423015224859268   2.423015224859268
         Total Energy (MeV)  :    2.476444998409542   2.476444998409542
         Kinetic Energy (MeV):    1.965445938409542   1.965445938409542
         Velocity (mm/ns)    :    293.3407906754122   293.3407906754122
         Volume Name         :              expHall             expHall
         Safety (mm)         :                5e-10               5e-10
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :           Geom Limit          Geom Limit
         Process defined Step:       Transportation      Transportation
      -----------------------------------------------------------------------

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0

 **PostStepDoIt (after all invocations):
    ++List of invoked processes 
      1) Transportation
      2) msc
      3) eIoni (Forced)

    ++G4Step Information 
      Address of G4Track    : 0x3046970
      Step Length (mm)      : 0
      Energy Deposit (MeV)  : 0
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -8.257551022101929  -8.257551022101929
         Position - y (mm)   :             -9.52505            -9.52505
         Position - z (mm)   :    18.49198141201282   18.49198141201282
         Global Time (ns)    :  0.06304591361302052 0.06304591361302052
         Local Time (ns)     :  0.06304591361302052 0.06304591361302052
         Proper Time (ns)    :  0.01300913309755812 0.01300913309755812
         Momentum Direct - x : -0.01057780492354854-0.01057780492354854
         Momentum Direct - y : 7.285445476253925e-057.285445476253925e-05
         Momentum Direct - z :   0.9999440508024575  0.9999440508024575
         Momentum - x (MeV/c): -0.02563161644372119-0.02563161644372119
         Momentum - y (MeV/c): 0.00017653733021041490.0001765373302104149
         Momentum - z (MeV/c):    2.423015224859268   2.423015224859268
         Total Energy (MeV)  :    2.476444998409542   2.476444998409542
         Kinetic Energy (MeV):    1.965445938409542   1.965445938409542
         Velocity (mm/ns)    :    293.3407906754122   293.3407906754122
         Volume Name         :              expHall     physiBFieldZone
         Safety (mm)         :                5e-10               5e-10
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :           Geom Limit          Geom Limit
         Process defined Step:       Transportation      Transportation
      -----------------------------------------------------------------------

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0
      [Note]Secondaries from AlongStepDoIt included.

#Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
30133  -8.26 mm  -9.53 mm   1.85 cm   1.97 MeV     0 eV      0 fm   1.85 cm      expHall  Transportation

 >>AlongStepDoIt (after all invocations):
    ++List of invoked processes 
      1) Transportation
      2) msc
      3) eIoni
      4) eBrem

    ++G4Step Information 
      Address of G4Track    : 0x3046970
      Step Length (mm)      : 0
      Energy Deposit (MeV)  : 0
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -8.257551022101929  -8.257551022101929
         Position - y (mm)   :             -9.52505            -9.52505
         Position - z (mm)   :    18.49198141201282   18.49198141201282
         Global Time (ns)    :  0.06304591361302052 0.06304591361302052
         Local Time (ns)     :  0.06304591361302052 0.06304591361302052
         Proper Time (ns)    :  0.01300913309755812 0.01300913309755812
         Momentum Direct - x : -0.01057780492354854-0.01057780492354854
         Momentum Direct - y : 7.285445476253925e-057.285445476253925e-05
         Momentum Direct - z :   0.9999440508024575  0.9999440508024575
         Momentum - x (MeV/c): -0.02563161644372119-0.02563161644372119
         Momentum - y (MeV/c): 0.00017653733021041490.0001765373302104149
         Momentum - z (MeV/c):    2.423015224859268   2.423015224859268
         Total Energy (MeV)  :    2.476444998409542   2.476444998409542
         Kinetic Energy (MeV):    1.965445938409542   1.965445938409542
         Velocity (mm/ns)    :    293.3407906754122   293.3407906754122
         Volume Name         :      physiBFieldZone     physiBFieldZone
         Safety (mm)         :                5e-10               5e-10
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :           Geom Limit          Geom Limit
         Process defined Step:       Transportation      Transportation
      -----------------------------------------------------------------------

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0

 **PostStepDoIt (after all invocations):
    ++List of invoked processes 
      1) Transportation
      2) msc
      3) eIoni (Forced)

    ++G4Step Information 
      Address of G4Track    : 0x3046970
      Step Length (mm)      : 0
      Energy Deposit (MeV)  : 0
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -8.257551022101929  -8.257551022101929
         Position - y (mm)   :             -9.52505            -9.52505
         Position - z (mm)   :    18.49198141201282   18.49198141201282
         Global Time (ns)    :  0.06304591361302052 0.06304591361302052
         Local Time (ns)     :  0.06304591361302052 0.06304591361302052
         Proper Time (ns)    :  0.01300913309755812 0.01300913309755812
         Momentum Direct - x : -0.01057780492354854-0.01057780492354854
         Momentum Direct - y : 7.285445476253925e-057.285445476253925e-05
         Momentum Direct - z :   0.9999440508024575  0.9999440508024575
         Momentum - x (MeV/c): -0.02563161644372119-0.02563161644372119
         Momentum - y (MeV/c): 0.00017653733021041490.0001765373302104149
         Momentum - z (MeV/c):    2.423015224859268   2.423015224859268
         Total Energy (MeV)  :    2.476444998409542   2.476444998409542
         Kinetic Energy (MeV):    1.965445938409542   1.965445938409542
         Velocity (mm/ns)    :    293.3407906754122   293.3407906754122
         Volume Name         :      physiBFieldZone             expHall
         Safety (mm)         :                5e-10               5e-10
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :           Geom Limit          Geom Limit
         Process defined Step:       Transportation      Transportation
      -----------------------------------------------------------------------

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0
      [Note]Secondaries from AlongStepDoIt included.

#Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
30134  -8.26 mm  -9.53 mm   1.85 cm   1.97 MeV     0 eV      0 fm   1.85 cm   physiBFieldZone  Transportation


[and so on... ]

So I'm wondering why is my step length at the boundary being set to zero?  
Do I have a bad set of physics lists?  
Is the fact that this particle is nearly perpendicular to the field the problem? 

Any help anyone can provide is GREATLY appreciated.

Thanks,

Brian
Feedback Re: electron stuck on volume boundary with (local) uniform magnetic field  by Gumplinger Peter <Gumplinger Peter>,   21 Jul, 2008
Re: Sad electron stuck on volume boundary with (local) uniform magnetic field (Brian Fisher)
Hi Brian,

I can see nothing wrong with your code snippet. Unfortunately, the G4 developer who can best assist you with your problem is going 'off line' for some weeks summer vacation, starting today. One suggestion, what happens when you change:

fIntgrDriver = new G4MagInt_Driver(1.0*micrometer, fStepper, fStepper->GetNumberOfVariables(), 2);

the 1.0*micrometer to some larger value?

I have also noticed that none of our 'official examples' use G4ExactHelixStepper.

What about if you used the method:

fieldMgr->CreateChordFinder(fMagField);

instead?

This seems to be using G4ClassicalRK4 even for constant magn. fields.

More Re: electron stuck on volume boundary with (local) uniform magnetic field  by Brian Fisher <Brian Fisher>,   21 Jul, 2008
Re: Feedback Re: electron stuck on volume boundary with (local) uniform magnetic field (Gumplinger Peter)
Thanks for your help. Sadly, neither solution worked.

Try 1: I adjusted the minimum steplength in 

 IntgrDriver = new G4MagInt_Driver(1.0*micrometer, fStepper, fStepper->GetNumberOfVariables(), 2);

from 1 m to 1 cm to 1 mm to 1 micrometer, and the electron seems to get stuck at the same point each time.

I then tried explicitly changing the stepper used. I tried G4ExactHelixStepper (as before), G4ClassicalRK4, and G4HelixImplicitEuler. The electron gets stuck at the same place.

Try 2: I then commented out the stuff that manually sets the stepper and driver and was left with just: 

	fMagField = new G4UniformMagField(G4ThreeVector(fFieldValue, 0.0, 0.0));
	fEquation = new G4Mag_UsualEqRhs(fMagField);
	fFieldMgr = new G4FieldManager(fMagField);
	fFieldMgr->CreateChordFinder(fMagField); 
	logicBFieldZone->SetFieldManager(fFieldMgr, true);

and got the same results.

I can't help but notice that the momentum vector of this particular electron is nearly entirely in the +z direction: 

WARNING - G4Navigator::ComputeStep()
          Track stuck, not moving for 10 steps
          in volume -expHall- at point (-8.257530138,-9.52505,18.49000716)
          direction: (-0.01057780492,7.290631171e-05,0.9999440508).
          Potential geometry or navigation problem !
          Trying pushing it of 9e-10 mm ...

The Lorentz force on it wants to push it back into the BFieldZone volume, though the electron's momentum is carrying it across the boundary into a zero-field region. Could this be why this is happening?

Thanks for any more help you can give.

More Re: electron stuck on volume boundary with (local) uniform magnetic field  by Brian Fisher <Brian Fisher>,   21 Jul, 2008
Re: More Re: electron stuck on volume boundary with (local) uniform magnetic field (Brian Fisher)
err... whoops. I got my directions (and charges) confused. The Lorentz force is is pushing the electron in the -Y direction in this case, out of the BFieldZone volume.

Question non-uniform electric field in chamber  by <janitor0405@126.com>,   28 Jun, 2008 
Dear sir,
   I am a beginner of Geant4,I want to setup a non-uniform electric field in my chamber filled with Xenon,I follow the example from/advanced/purging_magnet,I just made a conversion from the magnetic field to electric field,compiling is completed but when beamOn I receive some strange information as follows:

---> Begin Of Event: 0
G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper
  Step's start x=7.24397 and end x= 7.24397 are equal !!
  Due to step-size= 5.62767e-20 . Note that input step was 0.562767
G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper
  Step's start x=0.410347 and end x= 0.410347 are equal !!
  Due to step-size= 1.08637e-20 . Note that input step was 10.8637
G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper
                           ....
I really don't know what is the problem,I need some help.can you help me?thank you very much.
                                                                   janitor
None Re: non-uniform electric field in chamber  by <janitor0405@126.com>,   01 Jul, 2008
Re: Question non-uniform electric field in chamber 
Hello,everyone
     I modify some of my codes as follows: 
********my detectorconstruction.cc:

     G4FieldManager *pFieldMgrEF3 = new G4FieldManager();
    G4ElectricField* PurgEl3Field= new TabulatedElField2D3D("myelectric.txt"); 

      pFieldMgrEF3=G4TransportationManager::GetTransportationManager()->GetFieldManager();

     G4EqMagElectricField *myEquationEF3 = new G4EqMagElectricField(PurgEl3Field);

      G4MagIntegratorStepper* myStepperEF3 = new G4ClassicalRK4(myEquationEF3,8);
      G4MagInt_Driver* myIntgrDriverEF3 = new G4MagInt_Driver(1.0*mm,    myStepperEF3,
                                           myStepperEF3->GetNumberOfVariables());
      G4ChordFinder* myChordFinderEF3 = new G4ChordFinder(myIntgrDriverEF3);

      pFieldMgrEF3->SetChordFinder(myChordFinderEF3);
      pFieldMgrEF3->SetDetectorField(PurgEl3Field);
      pFieldMgrEF3->SetFieldChangesEnergy(true);

      bool fieldIsInitialized = true;
*******myelectric.txt
     6 6 10
    1 X
    2 Y
    3 Z
    4 EX
    5 EY
    6 EZ
    0 [MILLIMETRE]
-35.93 -35.93 47.62 -18566.7849 -21733.4192 -0.6929
-35.93 -35.93 58.91 -26282.7828 -30768.8520 -10.2787
  .........
   but when I run it,I receive Segmentation fault,I really don't know why,can you give me a hand? Thank you very much.
                                                               jane

Feedback Re: non-uniform electric field in chamber  by Gumplinger Peter <Gumplinger Peter>,   02 Jul, 2008
Re: None Re: non-uniform electric field in chamber
Hi Jane,

You should change the line:

G4FieldManager *pFieldMgrEF3 = new G4FieldManager();

to simply read:

G4FieldManager *pFieldMgrEF3;

The former is only required if you need a new G4FieldManager (as is the case for fields attached to logical volumes). When you get a pointer to the global field manager, with:

G4TransportationManager::GetTransportationManager()->GetFieldManager();

the pointer already exists and you don't have to get it from the heap.

I don't know if this is the cause of your seg-fault in the absence of knowing the message you get when the program crashes. It could be in your private code 'TabulatedElField2D3D'.

Please, also consult the /examples/extended/field02 and field04

None Re: non-uniform electric field in chamber  by <janitor0405@126.com>,   04 Jul, 2008
Re: Feedback Re: non-uniform electric field in chamber (Gumplinger Peter)
  Thank you for your reply.
  I modify my codes as you said and everything goes well without seg-fault.Thank you very much,but when I run my programme with beamOn,the old problem comes,It says:

  ---> Begin Of Event: 1
  G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper
  Step's start x=9.10312 and end x= 9.10312 are equal !!
  Due to step-size= 8.85047e-20 . Note that input step was 0.885047

  and it stops without any messages.I try to change my input electric values,but that problem still exits.what is the possible reason of it?Can you give me some suggestion?Thank you for your time!

                                                 jane

None Gradient in G4QuadrupoleMagField  Keywords: quadrupole magnetic field
by Mariusz <Mariusz>,   27 Feb, 2008 

Hello,
  Im just asking because I have not found any example.
  The gradient in G4QuadrupoleMagField(gradient) can be given in tesla/m?
  G4double gradient=223*tesla/m; 

   Regards, 

        Mariusz Sapinski

PS actually is this code correct for setting up the Quadrupole mag field in
only a part of geometry (logicSpace_Pipe):
   G4double fGradient = -223.*tesla/m;
   G4QuadrupoleMagField* pipeField = new G4QuadrupoleMagField(fGradient);
   G4Mag_UsualEqRhs* myEquation = new G4Mag_UsualEqRhs(pipeField);
   G4MagIntegratorStepper* myStepper = new G4ClassicalRK4(myEquation);
   G4ChordFinder* myChordFinder = new G4ChordFinder(pipeField,0.1*mm,myStepper);

   G4FieldManager* fieldMgr = new G4FieldManager(pipeField,myChordFinder,false);

   logicSpace_Pipe->SetFieldManager(fieldMgr,true);

I'm asking because it gives warnings in the parts of the geometry which are
quite far from Space_Pipe:
WARNING - G4Navigator::ComputeStep()
          Track stuck, not moving for 10 steps
          in volume -Austenitic_Steel_Collar- at point (89.47743818,103.7805885,1461.432287)
          direction: (0.1819718324,-0.1776272256,0.9671270966).
          Potential geometry or navigation problem !
          Trying pushing it of 9e-10 mm ...
None Re: Gradient in G4QuadrupoleMagField  Keywords: quadrupole magnetic field
by Tatiana Nikitina <Tatiana Nikitina>,   28 Feb, 2008
Re: None Gradient in G4QuadrupoleMagField (Mariusz)

Dear Mariusz

Yes, you can use tesla*m for gradient value. Important thing is to use the units, this make you sure that used value is correct. Using only local Field Manager is also possible and from your code it is good defined.

From your message is difficult to say why the warning message is appears, but it can be a overlapping or problem with a solid.

Can you,please, give more details about your geometry and the volume -Austenitic_Steel_Collar-?

Best Regards, Tatiana

Question Local Electric Field: GetFieldValue() problem  Keywords: coordinate systems nonuniform electric field
by Mario <Mario>,   15 Feb, 2008 
Hello Geant4 Gurus!!

I am using geant4.8.2.p01 currently and am implementing an electric field,
and later a magnetic field, to a specified local volume.  I have defined
my own electric field class, all of the code is below for that class.
The main concern is point[4] which comes in at the GetFieldValue() function.
I have displayed point[0], point[1], point[2] and calculated the magnitude
of this vector and in some instances it is larger than my entire world volume, 
which obviously makes no sense.  I guess I am curious about the coordinate 
system that point[4] is using... does it use the logicalVolume in which 
the fieldManager of the field belongs to as the coordinate system, or is it 
with respect to the global volume??  

Also, when I place this field_vol into the global and it overlaps with other
volumes already there or new ones that are placed on, it does not give me a
'overlap' warning, and the protons (which is what I am using) go through 
all volumes, overlapping or not, and it seems okay... is there something
special about a volume that is made of Vacuum and holds a field manager??

The main issue is what this point[4] is, because the EField is calculated
based on that number, so if I am using it incorrectly then there is a problem.
Thanks in advance for your help.

Cheers,
Mario 






HEADER:
#ifndef ChargeShell_hh
#define ChargeShell_hh 1

#include "G4Types.hh"
#include "G4ThreeVector.hh"
#include "G4ElectricField.hh"

class ChargeShell : public G4ElectricField
{
public: 

ChargeShell();
virtual ~ChargeShell() ;
virtual void GetFieldValue(const G4double point[4],
 G4double *field) const;

G4ThreeVector	CalcVolumeCharge(const G4double[4]) const;
G4ThreeVector	CalcSurfaceCharge(const G4double[4]) const;

G4double	GetLength()		{return 10.0*ro;};
G4ThreeVector	GetPos()		{return srcToShellPos;};

private:

	void ReadInputFile();

	G4double		ri;
	G4double		ro;
	G4double		Q;
	G4ThreeVector	srcToShellPos;

	G4double		fourPiEpsNot;
};
#endif


CC FILE:
#include "ChargeShell.hh"

extern int readFileVerbose;

ChargeShell::ChargeShell()
{
	ReadInputFile();
	fourPiEpsNot = 4.0*3.14159265359*8.85418781762e-12*farad/meter;
}

// ------------------------------------------------------------------------

ChargeShell::~ChargeShell()
{
}

void ChargeShell::GetFieldValue (const G4double point[4],G4double *fieldOut) const
{
	if(ri > ro){
		G4cout << "ri cannot be larger than ro" << G4endl;
		return;
	}
	
	G4ThreeVector fieldVal = G4ThreeVector(0.0,0.0,0.0);
	if(ro == ri){fieldVal = CalcSurfaceCharge(point);}
	else{fieldVal = CalcVolumeCharge(point);}

	G4cout << "fieldVal.x() = " << fieldVal.x()/(volt/m) << G4endl;
	G4cout << "fieldVal.y() = " << fieldVal.y()/(volt/m) << G4endl;
	G4cout << "fieldVal.z() = " << fieldVal.z()/(volt/m) << G4endl << G4endl;

	fieldOut[0]= 0.0;
	fieldOut[1]= 0.0;
	fieldOut[2]= 0.0;
	fieldOut[3]= fieldVal.x();
	fieldOut[4]= fieldVal.y();
	fieldOut[5]= fieldVal.z();
	return;
}

G4ThreeVector ChargeShell::CalcSurfaceCharge(const G4double p[4]) const 
{
	G4double x = p[0];
	G4double y = p[1];
	G4double z = p[2];
	G4double r = std::sqrt(x*x + y*y + z*z);

	G4double Er;
	G4cout << "r = " << r/mm << G4endl;
	if(r <= ro){Er = 0.0*volt/m;G4cout << "----->>>>>Er = 0<<<<<-----" << G4endl;}
	else{Er = Q/(fourPiEpsNot*r*r);}
	G4cout << "Er = " << Er/(volt/m) << G4endl;

	G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
	return out;
}

G4ThreeVector ChargeShell::CalcVolumeCharge(const G4double p[4]) const
{
	G4double x = p[0];
	G4double y = p[1];
	G4double z = p[2];
	G4double r = std::sqrt(x*x + y*y + z*z);

	double riOverro = ri/ro;
	double riOverr = ri/r;

	G4double Er;
	if(r >= 0.0 && r <= ri){Er = 0.0*volt/m;}
	else if(r > ri && r < ro){
		Er = (Q/(fourPiEpsNot*ro*ro))*(r/ro)*(1-riOverr*riOverr*riOverr)/(1-riOverro*riOverro*riOverro);
	}
	else{Er = Q/(fourPiEpsNot*r*r);} 
	
	G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
	return out;
}

void ChargeShell::ReadInputFile()
{
	//read in parameters from text file
	FILE *chargeShellInput;
	char* inputFile = "C:\\g4work\\PRSim\\src\\Param_Files\\Param_ChargeShellField.txt";
	chargeShellInput = fopen(inputFile,"rt");

	if(chargeShellInput == NULL){
		G4cout << "Param_ChargeShellField.txt did not come in correctly, check path" << G4endl;
		G4cout << "chargeShellInput = " << chargeShellInput << G4endl;
		G4cout << "errno = " << errno << G4endl;
		return;
	}	

	char readWholeLineBuff[500];
	float tmpFloat;
	char tmpString[50];

	//read line of text for ri
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		ri = tmpFloat*um;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << ri/um << G4endl;
		}
	}
	else{G4cout << "Did not read in ri" << G4endl; return;}

	//read line of text for ro
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		ro = tmpFloat*um;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << ro/um << G4endl;
		}
	}
	else{G4cout << "Did not read in ro" << G4endl; return;}

	//read line of text for Q
	fgets(readWholeLineBuff,500,chargeShellInput);
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		Q = tmpFloat*coulomb;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << Q/coulomb << G4endl;
		}
	}
	else{G4cout << "Did not read in Q" << G4endl; return;}

	//read line of text for srcToShellPos
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f",&tmpFloat);
		srcToShellPos.setX(tmpFloat*cm);
		fscanf(chargeShellInput,"%f",&tmpFloat);
		srcToShellPos.setY(tmpFloat*cm);
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		srcToShellPos.setZ(tmpFloat*cm);
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << srcToShellPos << G4endl;
		}
	}
	else{G4cout << "Did not read in srcToShellPos" << G4endl; return;}
	fclose(chargeShellInput);	
	return;
}
Feedback Re: Local Electric Field: GetFieldValue() problem  Keywords: coordinate systems nonuniform electric field
by Gumplinger Peter <Gumplinger Peter>,   15 Feb, 2008
Re: Question Local Electric Field: GetFieldValue() problem (Mario)

Hi Mario,

My guess is that none of us will speculate what happens when you have overlapping volumes because it is a NO NO in Geant4, precisely because what the G4Navigator will do is unpredictable. You have to use special (G4) tools to test your geometry for overlaps. Please, see the documentation about that. 

> The main concern is point[4] which comes in at the GetFieldValue()
> function.

point[3] is 'time' - ergo, you can have a time dependent field. 

> I have displayed point[0], point[1], point[2] and calculated the
> magnitude of this vector and in some instances it is larger than my
> entire world volume

You can also turn on step-wise debugging (/tracking/verbose 1) to follow where your tracks go and since you seem to have a debugging print statement in your SetFieldValue method you can find out when it's called - should only be called when in your 'logicalVolume' - and at what global coordinate points. If you print Point[0,1,2] you should be able to verify whether Point is local or global. I think it is global - but I am not 100% sure. So, if you want to return a value easier calculated from the local system, you need to use a GlobalToLocal transform:

G4StepPoint* p = aStep->GetPreStepPoint();

G4ThreeVector globalPoint(Point[0],Point[1],Point[2]);

const G4AffineTransform global2local = p->GetTouchable()->GetHistory()->GetTopTransform();

G4ThreeVector localPoint= global2local.TransformPoint(globalPoint);

Except I am at a loss how you get a handle to G4Step aStep in your SetFieldValue method. One lengthy way might be:

G4Step aStep = G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager()->GetStep(); 

> I guess I am curious about the coordinate system that point[4] is
> using... does it use the logicalVolume in which the fieldManager of
> the field belongs to as the coordinate system, or is it with respect
> to the global volume??

I am afraid it's the global coordinate system. 

> is there something special about a volume that is made of Vacuum and
> holds a field manager??

Not that I know off.

None Re: Local Electric Field: GetFieldValue() problem  Keywords: coordinate systems nonuniform electric field
by Mario <Mario>,   15 Feb, 2008
Re: Feedback Re: Local Electric Field: GetFieldValue() problem (Gumplinger Peter)
Thanks Peter, that helps...knowing about the global coordinate system.  I 
also just found a .pdf presentation from John Apostolakis where it says
that the point[] and field are in global coordinates.  

But to clarify something, I was printing point[0,1,2] and magnitude whenever
the GetFieldValue() function was being called.  So I saw a print out every
time the function was called, but after looking at the /tracking/verbose 1
output, I have seen that GetFielValue() is called a large number of times
during one track step... am I right? ... Therefore, when I was seeing these
values I was thinking that it was calculating incorrect field values... but
if the program calculates many values of the E-Field at radically different 
places, perhaps based on curvature or something of this nature, then it makes
sense that I am seeing all of these values, even if they are not directly 
used as the E-Field value... but looking at the tracking/verbose 1 I see that
the protons are still only spending one step through the volume which holds
the field... see implementation:

void EMFieldConstruction::SetFieldParameters(G4ElectricField* field)
{
    // Create an equation of motion for this field
    equation = new G4EqMagElectricField(field); 
    
    //set up a 4th Order Runge Kutta stepper for this equation
    G4int nvar = 8;
    stepper = new G4ClassicalRK4(equation, nvar);       
    
    //set up the driver to make sure that errors are under control by talking with stepper
    SetMinStep(1.0*um);
    intgrDriver = new G4MagInt_Driver(minFieldStep, stepper, stepper->GetNumberOfVariables());
    
    //define the chord finder
    chordFinder = new G4ChordFinder(intgrDriver);
    
    // Set this field and chord to the global field manager 
    fieldManager = new G4FieldManager(field,chordFinder);
    
}
void    EMFieldConstruction::SetMinStep(G4double s){minFieldStep = s;}

I have changed the minFieldStep to various values, but there is still 
only one step within the volume... is there another variable I should be 
setting?? The volume is Vacuum, but has a radially varying E-Field, so 
I would like there to be more than one step inside of the volume for accuracy.

Thank you Peter for a very timely response to the previous post,it is 
very appreciated.

Cheers,
Mario
Feedback Re: Local Electric Field: GetFieldValue() problem  Keywords: coordinate systems nonuniform electric field
by Gumplinger Peter <Gumplinger Peter>,   18 Feb, 2008
Re: None Re: Local Electric Field: GetFieldValue() problem (Mario)

> but after looking at the /tracking/verbose 1 output, I have seen

> that GetFielValue() is called a large number of times during one track

> step... am I right?

Yes, this is a feature of the Runge Kutta stepper. 

> I have changed the minFieldStep to various values, but there is still

> only one step within the volume... is there another variable I should

> be setting??

Looks to me like you want to set the max. step (not the min. step). I would have thought that in a highly varying field the RK method will reduce the step size automatically to adjust for this. I don't know why that's not the case for your field.

Anyway, you can force a max. step with:

G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();

tmanager->GetPropagatorInField()->SetLargestAcceptableStep(MaxTrackingStep);

None Problem with magnetic field from a table  Keywords: field map
by Mariusz <Mariusz>,   25 Jan, 2008 

 Hi,

  I have a 2d table with values of Bx and By (Bz is 0). It is quite large: 4900 points
with spacing of 7.7mm. I use Geant 4.9.0.p01.
  I have a class:

---------------------VVV ------------------------
class WSMagneticField: public G4MagneticField
{
  public:

   WSMagneticField(G4ThreeVector);  //  The value of the field
   WSMagneticField();               //  A zero field
  ~WSMagneticField();  

   void SetFieldValue(G4double fieldValue);
   void SetFieldValue(G4ThreeVector fieldVector);

   void GetFieldValue(const double Point[3],
                double *Bfield) const;           

  protected:
      // Find the global Field Manager
      G4FieldManager* GetGlobalFieldManager();   
  private:
      static const int imesh=4900;
      float x[imesh];
      float y[imesh];
      float Bx[imesh];
      float By[imesh];   
};
----------------------AAA---------------------
In the constructor I fill the tables x,y,Bx,By:

---------------------VVV----------------------
WSMagneticField::WSMagneticField() 
{
        // here open and read file      
        ifstream bmap;
        string line;
        const char *charline;
        int ni1,ni2,idx=0;
        float modB;
        char ni3;

        bmap.open("Bfield.data");
        while(bmap) {
                getline(bmap,line);
                charline=line.c_str();
                sscanf(charline,"  %d  %d  %f  %f  %f  %f  %f  %c\n",
                   &ni1,&ni2,&x[idx],&y[idx],&Bx[idx],&By[idx],&modB,&ni3);
                G4cout << "mag field check i = " << idx << " x = " << x[idx] << " Bx = " << Bx[idx] <<G4endl;
                idx++;                   
        }       

  GetGlobalFieldManager()->SetDetectorField(this);
  GetGlobalFieldManager()->CreateChordFinder(this);
}
--------------------AAA-------------------------

and I have GetFieldValue function: --------------------VVV------------------------- 

// argument is position in mm
void WSMagneticField::GetFieldValue(const double Point[3], double *Bfield) const
{
        Bfield[0] = 0.;
        Bfield[1] = 0.;
        Bfield[2] = 0.;
        G4float dmesh=7.93; // [mm]
        G4double cPoint[3];
        cPoint[0]=Point[0];
        cPoint[1]=Point[1]-77.; // [mm] shift back the Quad to World ref.    
        cPoint[2]=Point[2];     
        G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);
        if(cPoint[0]>0 && radius < 240. && fabs(cPoint[2]) < 1400.) {
         for(G4int i=0; i<3; i++) if(cPoint[i]<0) cPoint[i]*=(-1); // because of symmetry
         for(G4int i=0; i<4900; i++) {

                if(fabs(cPoint[0]-x[i])<dmesh/2 && fabs(cPoint[1]-y[i])<dmesh/2) {
                                Bfield[0]=Bx[i];
                                Bfield[1]=By[i];
                                break;
                }
         }
        }       
}
-------------------AAA------------------------- 

And then, in DetectorConstruction I have:
-------------------VVV-------------------------
        WSMagneticField* myField = new WSMagneticField();
        G4FieldManager* fieldMgr
                                = G4TransportationManager::GetTransportationManager()
                                ->GetFieldManager();

// I tried this and it did not work either // fieldMgr->SetDetectorField(myField); // fieldMgr->CreateChordFinder(myField);

/* ... geometry declaration ... */ 

   G4FieldManager *localFieldMgr = new G4FieldManager(myField);
   logicQ5->SetFieldManager(localFieldMgr,true);
--------------------AAA-------------------------

The program compiles, the file with table of values is read but never ever any event has finished being calculated. I've left it even for the whole night even if in normal conditions one event takes about 15 minutes.

I would appreciate any suggestion leading to solution;-) 

 Mariusz Sapinski

None Re: Problem with magnetic field from a table  by John Apostolakis <John Apostolakis>,   25 Jan, 2008
Re: None Problem with magnetic field from a table (Mariusz)
Dear Mariusz,

Geant4 uses the GetFieldValue method that you create many times during 
every step.  Looking at the method that you have written, I can see that 
it is making a very large number of operations, in order to calculate 
the value.  So it is not surprising that it your program is very slow.

I have two suggestion to speed up your key GetFieldValue method:

1) Avoid large loops:

        for(G4int i=0; i<4900; i++) {
                if(fabs(cPoint[0]-x[i])<dmesh/2 ...

     Use the functions for converting a floating point number to a 
integer, in order to calculate the nearest 'i' integer to your point.  
Some code along the lines of:
                 double ratio = (cPoint[0]/dmesh);
                 int i = round( ratio );   // Double check this - it 
should round to nearest integer maybe it is std::round
                 if( fabs( cpoint[0]- x[i]) < dmesh/2 ) ...
etc

2)  to replace the line

G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);

with

G4double radius = std::sqrt(Point[0]*cPoint[0]+cPoint[1]*cPoint[1]);

since
  - the power function is very costly compared to a multiplication, and
  - square root is faster than pow( .. , 0.5)
The power method generally requires a logarithm and an exponential.

I expect that this will speed up your field method, and your 
application, by a factor of 100 or more.  I expect that your particles 
will actually move.

Best regards,
John

Mariusz wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/146"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
>  Hi,
>
>   I have a 2d table with values of Bx and By (Bz is 0). It is quite large: 4900 points
> with spacing of 7.7mm. I use Geant 4.9.0.p01.
>   I have a class:
>
> ---------------------VVV ------------------------
> class WSMagneticField: public G4MagneticField
> {
>   public:
>
>    WSMagneticField(G4ThreeVector);  //  The value of the field
>    WSMagneticField();               //  A zero field
>   ~WSMagneticField();  
>
>    void SetFieldValue(G4double fieldValue);
>    void SetFieldValue(G4ThreeVector fieldVector);
>
>    void GetFieldValue(const double Point[3],
>                 double *Bfield) const;           
>
>   protected:
>       // Find the global Field Manager
>       G4FieldManager* GetGlobalFieldManager();   
>   private:
>       static const int imesh=4900;
>       float x[imesh];
>       float y[imesh];
>       float Bx[imesh];
>       float By[imesh];   
> };
> ----------------------AAA---------------------
> In the constructor I fill the tables x,y,Bx,By:
>
> ---------------------VVV----------------------
> WSMagneticField::WSMagneticField() 
> {
>         // here open and read file      
>         ifstream bmap;
>         string line;
>         const char *charline;
>         int ni1,ni2,idx=0;
>         float modB;
>         char ni3;
>
>         bmap.open("Bfield.data");
>         while(bmap) {
>                 getline(bmap,line);
>                 charline=line.c_str();
>                 sscanf(charline,"  %d  %d  %f  %f  %f  %f  %f  %c\n",
>                    &ni1,&ni2,&x[idx],&y[idx],&Bx[idx],&By[idx],&modB,&ni3);
>                 G4cout << "mag field check i = " << idx << " x = " << x[idx] << " Bx = " << Bx[idx] <<G4endl;
>                 idx++;                   
>         }       
>
>   GetGlobalFieldManager()->SetDetectorField(this);
>   GetGlobalFieldManager()->CreateChordFinder(this);
> }
> --------------------AAA-------------------------
>
> and I have GetFieldValue function:
> --------------------VVV-------------------------
>
> // argument is position in mm
> void WSMagneticField::GetFieldValue(const double Point[3], double *Bfield) const
> {
>         Bfield[0] = 0.;
>         Bfield[1] = 0.;
>         Bfield[2] = 0.;
>         G4float dmesh=7.93; // [mm]
>         G4double cPoint[3];
>         cPoint[0]=Point[0];
>         cPoint[1]=Point[1]-77.; // [mm] shift back the Quad to World ref.    
>         cPoint[2]=Point[2];     
>         G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);
>         if(cPoint[0]>0 && radius < 240. && fabs(cPoint[2]) < 1400.) {
>          for(G4int i=0; i<3; i++) if(cPoint[i]<0) cPoint[i]*=(-1); // because of symmetry
>          for(G4int i=0; i<4900; i++) {
>
>                 if(fabs(cPoint[0]-x[i])<dmesh/2 && fabs(cPoint[1]-y[i])<dmesh/2) {
>                                 Bfield[0]=Bx[i];
>                                 Bfield[1]=By[i];
>                                 break;
>                 }
>          }
>         }       
> }
> -------------------AAA------------------------- 
>
> And then, in DetectorConstruction I have:
> -------------------VVV-------------------------
>         WSMagneticField* myField = new WSMagneticField();
>         G4FieldManager* fieldMgr
>                                 = G4TransportationManager::GetTransportationManager()
>                                 ->GetFieldManager();
>
> // I tried this and it did not work either //
> fieldMgr->SetDetectorField(myField); //
> fieldMgr->CreateChordFinder(myField);
>
> /* ... geometry declaration ... */
>
>    G4FieldManager *localFieldMgr = new G4FieldManager(myField);
>    logicQ5->SetFieldManager(localFieldMgr,true);
> --------------------AAA-------------------------
>
> The program compiles, the file with table of values is read but never
> ever any event has finished being calculated. I've left it even for the
> whole night even if in normal conditions one event takes about 15
> minutes.
>
> I would appreciate any suggestion leading to solution;-)
>
>  Mariusz Sapinski  
>
>
None Re: Problem with magnetic field from a table  by Mariusz <Mariusz>,   27 Jan, 2008
Re: None Re: Problem with magnetic field from a table (John Apostolakis)

Actually I met a problem when running simulation for higher energy particles. My program gave core dump. It looks like it happend at the very end of an event.

Here is the output: 

(...)
MRoot Event::Clear() Event had 5517 particles and 36687 CCells
WARNING - G4PropagatorInField::LocateIntersectionPoint()
          Convergence is requiring too many substeps: 10010
          Abandoning effort to intersect.
          Information on start & current step follows in cout.
WARNING - G4PropagatorInField::LocateIntersectionPoint()
          Convergence is requiring too many substeps: 10010
          Found intersection = 0
          Intersection exists = 1
          Start and Endpoint of Requested Step:
       Current Position  and  Direction
Step#       s        X(mm)      Y(mm)      Z(mm)    N_x     N_y     N_z   Delta|N|   StepLen  StartSafety   PhsStep
Start         0  118.79256  77.074619      -1600 -0.1061 -6.902e-05  0.9944 8.33e-16         0         15.3 Init/NotKnown
   0  361985.68 -38277.146  52.089425  358343.59 -0.1061 -6.902e-05  0.9944 8.33e-16  3.62e+05         15.3 Init/NotKnown

'Bracketing' starting and endpoint of current Sub-Step 10009 3108.2408 97.600574 76.871888 1369.5735 0.1818 -6.081e-05 0.9833 1.55e-15 0 15.3 Init/NotKnown 10010 361985.68 -38277.146 52.089425 358343.59 -0.1061 -6.902e-05 0.9944 -3.23 3.59e+05 15.3 Init/NotKnown 

ERROR - G4PropagatorInField::LocateIntersectionPoint()
        Undertaken only length: 3108.240753 out of 361985.6821 required.
        Remaining length = 358877.4414

*** G4Exception : UnableToLocateIntersection
      issued by : G4PropagatorInField::LocateIntersectionPoint()
Too many substeps while trying to locate intersection.
*** Fatal Exception *** core dump ***

*** G4Exception: Aborting execution *** 

zsh: abort      bin/Linux-g++/mqy
None Re: Problem with magnetic field from a table  by John Apostolakis <John Apostolakis>,   28 Jan, 2008
Re: None Re: Problem with magnetic field from a table (Mariusz)
The problem that you report seems strongly related to your geometry:  
The particle seems not to encounter a volume boundary even at X= -38277 
(mm) and Z= 358343 (mm).   Is your geometry as large as this ?

Else likely there is a problem in navigating in your geometry.  Likely 
this is due to an error in constructing it.  To help identify it, I 
think that the volume involved is one of the volumes to which the 
initial point belongs:

X (mm)          Y(mm)          Z(mm)
118.79256 77.074619    -1600

I expect that it must lie on a boundary between volumes.  Can you check 
the volumes on either side, starting with the volume on the side of the 
initial direction ?

Could you explain the complexity of the geometry model used?

Regards,
John

Mariusz wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/146/1/3"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
> Actually I met a problem when running simulation for higher energy
> particles. My program gave core dump. It looks like it happend at the
> very end of an event.
>
> Here is the output:
>
> (...)
> MRoot Event::Clear() Event had 5517 particles and 36687 CCells
> WARNING - G4PropagatorInField::LocateIntersectionPoint()
>           Convergence is requiring too many substeps: 10010
>           Abandoning effort to intersect.
>           Information on start & current step follows in cout.
> WARNING - G4PropagatorInField::LocateIntersectionPoint()
>           Convergence is requiring too many substeps: 10010
>           Found intersection = 0
>           Intersection exists = 1
>           Start and Endpoint of Requested Step:
>        Current Position  and  Direction
> Step#       s        X(mm)      Y(mm)      Z(mm)    N_x     N_y     N_z   Delta|N|   StepLen  StartSafety   PhsStep
> Start         0  118.79256  77.074619      -1600 -0.1061 -6.902e-05  0.9944 8.33e-16         0         15.3 Init/NotKnown
>    0  361985.68 -38277.146  52.089425  358343.59 -0.1061 -6.902e-05  0.9944 8.33e-16  3.62e+05         15.3 Init/NotKnown
>
> 'Bracketing' starting and endpoint of current Sub-Step 10009 3108.2408
> 97.600574 76.871888 1369.5735 0.1818 -6.081e-05 0.9833 1.55e-15 0 15.3
> Init/NotKnown 10010 361985.68 -38277.146 52.089425 358343.59 -0.1061
> -6.902e-05 0.9944 -3.23 3.59e+05 15.3 Init/NotKnown
>
> ERROR - G4PropagatorInField::LocateIntersectionPoint()
>         Undertaken only length: 3108.240753 out of 361985.6821 required.
>         Remaining length = 358877.4414
>
> *** G4Exception : UnableToLocateIntersection
>       issued by : G4PropagatorInField::LocateIntersectionPoint()
> Too many substeps while trying to locate intersection.
> *** Fatal Exception *** core dump ***
>
> *** G4Exception: Aborting execution ***
>
> zsh: abort      bin/Linux-g++/mqy
>
>
None Re: Problem with magnetic field from a table  by Mariusz <Mariusz>,   26 Jan, 2008
Re: None Re: Problem with magnetic field from a table (John Apostolakis)

Thank you John, in fact it helped. The reduction of the computation time is by factor about 50, so you still owe me factor 2 ;-) 

I get quite a lot of warnings:
WARNING - G4MagIntegratorDriver::AccurateAdvance()
          Proposed step is zero; hstep = 0 !
WARNING - G4MagIntegratorDriver::AccurateAdvance()
          Proposed step is zero; hstep = 0 !
 G4Transportation is killing track that is looping or stuck
   This track has 333.65357 MeV energy.
   Number of trials = 10   No of calls to AlongStepDoIt = 13894548

Should I be worried? Tracks with 300 MeV can be important for me... 

  My best regards,

       Mariusz
None Re: Problem with magnetic field from a table  by John Apostolakis <John Apostolakis>,   28 Jan, 2008
Re: None Re: Problem with magnetic field from a table (Mariusz)
It is hard to know whether you are seeing something harmless or 
problematic.  It depends on what is causing the two problems - which may 
be independent or related.

It could related to the problem you reported today - and something like 
an overlap or volume protruding from a daughter volume past the 
boundaries of a mother.

To understand better, I/we would need some more information.  Can you 
try to run with the command /tracking/verbose 1, and put the output into 
a local file (eg in /tmp ) for a short/medium size run?  Then send us 
the lines preceding each of the three warning messages that you see, and 
the type of particle involved?

Could you give some more information about the particle type, and 
whether the location has zero or finite magnetic field?  I note that in 
the Geant4 tracking, if a field exists the code assumes it has a 
non-zero value.  To get the fastest possible conditions, you should give 
each volume a field manager (at least for each volume whose field you 
need to set to zero).  Then you can create a field manager for the null 
field, and set its field pointer to zero.  This ensures that the 
propagation code recognised these volumes as not having any field.

-- In case it is a real, occasional problem --
In case it is a simple difficulty it may be possible for you to overcome 
it if you change a parameter in G4Transportation which controls the 
number of trials (which is 10 by default, as printed).   You can try to 
double it, or more, by calling the following method of G4Transportation:
   void SetThresholdTrials(G4int <ident?v=9.0.p1;i=G4int> newMaxTrials );

If you are running locally on a machine with a debug version of Geant4, 
probably the fastest way to do this is using a debugger. Else get a 
pointer to the Transportation, you need to find the one created in your 
physics lists.  There is probably a better way to do it, but the 
following should work:
  - Choose a particle (eg electron), and get its particle definition eg 
G4Electron::Electron()
  - Get its physics manager
  - Ask it for its process vector  (using the method G4ProcessVector 
<ident?i=G4ProcessVector>* GetAlongStepProcessVector)
  - Go through this vector checking each process to see its name is 
G4Transportation.
Do this only once, eg in your main, before the event loop, and only 
after creating your physics list.

Best regards,
John Apostolakis

Mariusz wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/146/1/1"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
> I get quite a lot of warnings:
> WARNING - G4MagIntegratorDriver::AccurateAdvance()
>           Proposed step is zero; hstep = 0 !
> WARNING - G4MagIntegratorDriver::AccurateAdvance()
>           Proposed step is zero; hstep = 0 !
>  G4Transportation is killing track that is looping or stuck
>    This track has 333.65357 MeV energy.
>    Number of trials = 10   No of calls to AlongStepDoIt = 13894548
>
> Should I be worried? Tracks with 300 MeV can be important for me...
>
>   My best regards,
>
>        Mariusz
>
>
Question non-uniform magnetic field  by <adammachona@yahoo.co.uk>,   15 Nov, 2007

I'm a new Geant4 user and so you will have to forgive me if I ask a stupid question. I want to create a non-uniform magnetic field in my detector. What is the simplest way to do this? Thank you for your time.

None Re: non-uniform magnetic field  by John Apostolakis <John Apostolakis>,   15 Nov, 2007
Re: Question non-uniform magnetic field 
Dear user,

You will need to create a class that describes the field. The simplest 
example for doing this is in the N04 novice example.

You can find it on our LXR source code browser at
  http://www-geant4.kek.jp/lxr/source/examples/novice/N04/
or in the examples directory in the release.

You can copy this code, and change it for your purposes.

For more information please consult the documentation (Application 
Developers Guide).

Best regards,
John Apostolakis

adammachona@yahoo.co.uk wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/144"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
> 
> I'm a new Geant4 user and so you will have to forgive me if I ask a
> stupid question. I want to create a non-uniform magnetic field in my
> detector. What is the simplest way to do this? Thank you for your time.
None non-uniform electric field  by <janitor0405@126.com>,   12 Nov, 2007 

Dear sir,
   I am a beginner of Geant4,I want to setup a non-uniform electric field in my chamber filled with Xenon,
I follow the example from/advanced/purging_magnet,I just made a conversion from the magnetic field to electric field,
when gmakeing I receive some strange information about syntax error such as <<  and so on,
I already include the head file G4ios.hh and globals.hh,I really don't know what is the problem,I need some help.
thank you very much. 

                                                          Jane
Feedback Re: non-uniform electric field  by Gumplinger Peter <Gumplinger Peter>,   13 Nov, 2007
Re: None non-uniform electric field

Hello Jane,

Looks like you are having issues with C++ and less so with G4. - We are not really in the business of solving people's compilation errors but I'll take a look (in private) if you send me your relevant code and the syntax error messages; that is, if you can't find anyone physically near you to help you with this.

Peter

None Electric potential / non-uniform electric field  Keywords: non-uniform electric field, electric potential
by Hannele <hatueron@jyu.fi>,   17 Sep, 2007

Is it possible to add an electric potential to components? And if not, how can I create a non-uniform electric field?

None Re: Electric potential / non-uniform electric field  by John Apostolakis <John Apostolakis>,   17 Sep, 2007
Re: None Electric potential / non-uniform electric field (Hannele)
Dear Hanelle,

You cannot add an electrical potential to a component.

Regarding how to create a non-uniform electric field, did you look at 
the previous answers in HyperNews for answers to this question? 

As I recall, it has been answered a few times.  If you do not find this, 
let us know.

Regards,
John A.

Hannele wrote:
> *** Discussion tiwtle: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/140"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
> Is it possible to add an electric potential to components? And if not,
> how can I create a non-uniform electric field?
>
None Re: Electric potential / non-uniform electric field  by Hannele <hatueron@jyu.fi>,   18 Sep, 2007
Re: None Re: Electric potential / non-uniform electric field (John Apostolakis)

Dear John,

I searched for the answers in the previous messages to the forum, but unfortunately in vain. I maybe expressed myself not entirely exactly. What I meant by "adding a field to component" should have been "defining the field inside a volume". I have used an example from /advanced/purging_magnet to define a non-uniform magnetic field which is read from file containing information about the field. In order to do the same with the non-uniform electric field I have made a "conversion" of the magnetic field to electric (as shown below) feeding to the program a file with information about the electric field: 

G4ThreeVector positionEl3Field = G4ThreeVector(0,0,1*m);
   G4double fieldValueEF3 = 0.2*kilovolt/cm;
      G4FieldManager   *pFieldMgrEF3;

      //Field grid in A9.TABLE. File must be in accessible from run urn directory. 
      G4ElectricField* PurgEl3Field= new TabulatedElField2D3D("MMPS_in_edit.txt", 0,0,0);

      pFieldMgrEF3=G4TransportationManager::GetTransportationManager()->GetFieldManager();

      G4cout<< "DeltaStep "<<pFieldMgrEF3->GetDeltaOneStep()/mm <<"mm" <<endl;

      G4EqMagElectricField *myEquationEF3 = new G4EqMagElectricField(PurgEl3Field);

      G4int nvar3 = 8; // Need to integrate 8 variables: x,y,z,p[xyz],E,t 

      G4MagIntegratorStepper* myStepperEF3 = new G4ClassicalRK4(myEquationEF3,nvar3);
      G4MagInt_Driver* myIntgrDriverEF3 = new G4MagInt_Driver(1.0e-3*mm, myStepperEF3,
                                           myStepperEF3->GetNumberOfVariables());
      G4ChordFinder* myChordFinderEF3 = new G4ChordFinder(myIntgrDriverEF3);

      pFieldMgrEF3->SetDetectorField(PurgEl3Field);

      bool fieldIsInitialized = true;

Here, MMPS_in_edit.txt is the file with the electric field components (x,y,z,E). The program compiles, however, when I try to run it, the program returns "Segmentation fault" and quits without giving a hint of what could go wrong.

I would really appreaciate if you help me with this problem, Thank you.

Feedback Re: Electric potential / non-uniform electric field  by Gumplinger Peter <Gumplinger Peter>,   26 Sep, 2007
Re: None Re: Electric potential / non-uniform electric field (Hannele)

Try this:

G4FieldManager* pFieldMgrEF3 = new G4FieldManager();

G4ElectricField* PurgEl3Field= new TabulatedElField2D3D("MMPS_in_edit.txt",0,0,0);

G4EqMagElectricField* myEquationEF3 = new G4EqMagElectricField(PurgEl3Field);

G4MagIntegratorStepper* myStepperEF3 = new G4ClassicalRK4(myEquationEF3,8) 

G4MagInt_Driver* myIntgrDriverEF3 = new G4MagInt_Driver(1.0e-3*mm,
                                               myStepperEF3,
                                               myStepperEF3->GetNumberOfVariables());

G4ChordFinder* myChordFinderEF3 = new G4ChordFinder(myIntgrDriverEF3);

pFieldMgrEF3->SetDetectorField(PurgEl3Field);

pFieldMgrEF3->SetChordFinder(myChordFinderEF3);

pFieldMgrEF3->SetFieldChangesEnergy(true);

logicVolume_in_question->SetFieldManager(pFieldMgrEF3,allLocal);

see also /examples/extended/field/field02 for how to define electric fields (albeit G4UniformElectricField) and /field03 of how to attach it to a logical volume (instead of using the GlobalFieldManager)

None Overlapping Fields  by Maurizio <ungaro@jlab.org>,   31 Jul, 2007

I read in #74 that one can have electric and magnetic field at the same time by using the return values of bfield (components 0-2 and 3-5)

What about overlapping magnetic fields? If I multiple magnetic fields at one point in space, can I add them together somehow?

This is important when working with strong fringe fields.

thanks, mauri

Feedback Re: Overlapping Fields  by Gumplinger Peter <Gumplinger Peter>,   31 Jul, 2007
Re: None Overlapping Fields (Maurizio)

I understand why you are asking. I don't have the answer except to say that you, yourself, have to calculate the total field by adding the contributing fields at every point in space. A lot of work and coding is required and it's a shame since G4 already "knows" where all the field components are and their distance to point x/y. Now, in G4 you must not have overlapping volumes (field elements)!

For a better start, please take a look at posting #116 below by Tom Roberts.

Tom Roberts is the author of URL http://g4beamline.muonsinc.com . Its code is available, source and binary. I understand there is a workable solution to the problem of overlapping fringe fields from different beam elements.

More Re: Overlapping Fields  by Gumplinger Peter <Gumplinger Peter>,   07 Nov, 2007
Re: Feedback Re: Overlapping Fields (Gumplinger Peter)

G4 9.1 will have an extended example (/examples/extended/field/field04) which will show how to define/use OVERLAPPING field elements in Geant4. Fields might be either magnetic, electric or both. 

     Credit goes to Tom Roberts and Muons Inc. since much of the code
     and ideas were taken at liberty from the (GNU GPL) source of
     G4BEAMLINE release 1.12.

     http://g4beamline.muonsinc.com

None Questions about SetLargestAcceptableStep()   Keywords: SetLargestAcceptableStep
by <throwpen@163.com>,   21 Jul, 2007 
When do i need to set the maximum acceptable step using the 
SetLargestAcceptableStep() method ?
I know that if i want to get a smooth trajectory for a charged particle 
in the magnetic field, i need to set the maximum acceptable step to small
one.However, the tracking consumes too much time.

whether the maximum acceptable step setting will affect the final 
result or not? The dimension of detector in my simulation is <1 meter.
Feedback Re: Questions about SetLargestAcceptableStep()   Keywords: SetLargestAcceptableStep
by Gumplinger Peter <Gumplinger Peter>,   21 Jul, 2007
Re: None Questions about SetLargestAcceptableStep() 

> When do i need to set the maximum acceptable step using the 

> SetLargestAcceptableStep() method ?

for the visualization of smooth trajectories. 

> whether the maximum acceptable step setting will affect the final 

> result or not?

should not affect your results.

What can affect your results and why is illustrated in this presentation by John Apostolakis:

http://geant4.in2p3.fr/2007/prog/JohnApostolakis/EMField.pdf

None RE: Questions about SetLargestAcceptableStep()   by John Apostolakis <John Apostolakis>,   22 Jul, 2007
Re: Feedback Re: Questions about SetLargestAcceptableStep() (Gumplinger Peter)
The original use of  SetLargestAcceptableStep() was to make sure that a
low energy track that is spiraling in a field does not try step size
that would be impossible for the integration methods (eg in lab vacuum a
particle would likely be given a very very large step before the next
interaction, and this is a protection to ensure that work is not wasted
trying to integrate really far - whereas the setup/aparatus/detector is
of much smaller distances.)

It used to be that such a 'lever' could be needed to create smooth
tracks (for visualisation purposes only -- the tracks really are already
simulated as being curved!) - but now there are better ways to ensure
that the visualisation of tracks shows their curvature.   If you need
that, I suggest that you check the documentation for smooth trajectories
-- and ask in the forum for visualisation. 

In summary this parameter is not relevant for the purposes of accuracy.
It is meant to make sure that the field module is more robust -- and in
cases where the setup is very much larger than typical lab setups (eg
planet size) or much smaller, it is necessary to change its value.  Only
then.

Best regards,
John Apostolakis
------------------------------------------------------------------------
John Apostolakis                   Email:  John.Apostolakis@cern.ch
  

> -----Original Message-----
> From: Peter Gumplinger [mailto:gum@triumf.ca] 
> Sent: Saturday, July 21, 2007 7:10 PM
> To: PublicHyperNews@slac.stanford.edu
> Subject: Re: Questions about SetLargestAcceptableStep() 
> 
> *** Discussion title: Fields: Magnetic and Otherwise Email 
> replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/136/1"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
> 
> > When do i need to set the maximum acceptable step using the
> 
> > SetLargestAcceptableStep() method ?
> 
> for the visualization of smooth trajectories.
> 
> > whether the maximum acceptable step setting will affect the final
> 
> > result or not?
> 
> should not affect your results.
> 
> What can affect your results and why is illustrated in this 
> presentation by John Apostolakis:
> 
> http://geant4.in2p3.fr/2007/prog/JohnApostolakis/EMField.pdf
>
None Re: Questions about SetLargestAcceptableStep()   Keywords: SetLargestAcceptableStep
by <throwpen@163.com>,   22 Jul, 2007
Re: Feedback Re: Questions about SetLargestAcceptableStep() (Gumplinger Peter)
Today i tried it again. 

if i didn't set the maximum acceptable step in the program. 
I encountered the same problem as Ilya reported in 
http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/emfields/96.html
which was called "Electrons bouncing on the volume boundary in pure magnetic field"
When it happened, electrons didn't lost their energies in the volume(not vacumm).
After more than 1000000 steps, then suddenly, electrons stopped bounding and lost their energies by
MSC process. But a great deal of time was wasted in the boucing

*********************************************************************************************************
* G4Track Information:   Particle = e-,   Track ID = 11,   Parent ID = 1
*********************************************************************************************************

Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
    0       18     22.6    -27.8     0.001        0        0         0    collim21 initStep
    1       18     22.6    -27.8   0.00081 0.000195   5e-006    5e-006    collim21 msc
    2       18     22.6    -27.8  0.000758 5.14e-005   5e-006    1e-005    collim21 msc
    3       18     22.6    -27.8  0.000722 3.63e-005   5e-006  1.5e-005    collim21 msc
    4       18     22.6    -27.8  0.000714 8.37e-006 4.85e-007 1.55e-005     expHall Transportation
    5       18     22.6    -27.8  0.000714        0        0 1.55e-005    collim21 Transportation
    6       18     22.6    -27.8  0.000714        0        0 1.55e-005     expHall Transportation
    7       18     22.6    -27.8  0.000714        0        0 1.55e-005    collim21 Transportation
    8       18     22.6    -27.8  0.000714        0        0 1.55e-005     expHall Transportation
    9       18     22.6    -27.8  0.000714        0 4.19e-007 1.59e-005     expHall Transportation
   10       18     22.6    -27.8  0.000714        0 2.79e-006 1.87e-005    collim21 Transportation
   11       18     22.6    -27.8  0.000714        0        0 1.87e-005     expHall Transportation
   12       18     22.6    -27.8  0.000714        0        0 1.87e-005    collim21 Transportation
   13       18     22.6    -27.8  0.000714        0        0 1.87e-005     expHall Transportation
   14       18     22.6    -27.8  0.000714        0 3.51e-007  1.9e-005     expHall Transportation
   15       18     22.6    -27.8  0.000714        0 2.44e-006 2.15e-005    collim21 Transportation
   ....
   ....
   ....
3605172       18     22.8    -27.5  0.000491        0        0     0.491    collim21 Transportation
3605173       18     22.8    -27.5  0.000491        0        0     0.491     expHall Transportation
3605174       18     22.8    -27.5  0.000491        0 4.08e-007     0.491    collim21 Transportation
3605175       18     22.8    -27.5  0.000491        0        0     0.491     expHall Transportation
3605176       18     22.8    -27.5  0.000491        0        0     0.491    collim21 Transportation
3605177       18     22.8    -27.5  0.000491        0        0     0.491     expHall Transportation
3605178       18     22.8    -27.5  0.000491        0 3.11e-007     0.491    collim21 Transportation
3605179       18     22.8    -27.5  0.000491        0        0     0.491     expHall Transportation
3605180       18     22.8    -27.5  0.000491        0        0     0.491    collim21 Transportation
3605181       18     22.8    -27.5  0.000491        0        0     0.491     expHall Transportation
3605182       18     22.8    -27.5  0.000491        0 1.02e-007     0.491    collim21 Transportation
3605183       18     22.8    -27.5  0.000322 0.000169   5e-006     0.491    collim21 msc
3605184       18     22.8    -27.5 5.71e-005 0.000265   5e-006     0.491    collim21 msc
3605185       18     22.8    -27.5 3.97e-005 1.74e-005 7.72e-007     0.491     expHall Transportation
3605186       18     22.8    -27.5 3.97e-005        0        0     0.491    collim21 Transportation
3605187       18     22.8    -27.5         0 3.97e-005 2.32e-006     0.491    collim21 LowEnergyIoni
Track (trackID 11, parentID 1) is processed with stopping code 2
### pop requested out of 9 stacked tracks.


But if i set the maximum acceptable step as the following:

G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();
tmanager->GetPropagatorInField()->SetLargestAcceptableStep(0.1*mm);

the problem didn't show up.


Something about my program:

Primary particle: electrons(0.1 to 0.5MeV)
the direction of momentum: +x
Magnetic field: uniform field(in +Z direction,about 4000 Gauss)

void ExN01Field::GetFieldValue(const double Point[3],double *Bfield) const
{
  
  
  G4double px=(Point[0]/mm);
  G4double py=(Point[1]/mm);
  G4double pz=(Point[2]/mm);

	  Bfield[0]=0.001*tesla;
	  Bfield[1]=0.001*tesla;
	  Bfield[2]=0.001*tesla;

  if(px>=0&&px<=45)
	  if(py>=-30&&py<=30)
	  if(pz>=-30&&pz<=30)
	{
	  Bfield[0]=0*tesla;
	  Bfield[1]=0.4*tesla;
	  Bfield[2]=0*tesla;
	}
  

}
Question simultaneous B and E fields  Keywords: magnetic field, electric field
by Abigail Bickley <bickley@nscl.msu.edu>,   16 Jul, 2007

Hello,

I am trying to simulate a detector that has a local electric field applied to a subset of the detector volume and a global magnetic field. I understand from message #74 that this is possible using the G4ElectroMagneticField class, but it is unclear to me from the response to this message how to apply this class. If I setup global and local fields as in example field03 then do I simply convert 

  fField = new G4UniformElectricField(G4ThreeVector(0.0,0.0*volt,0.0 ));
  fLocalField = new G4UniformElectricField(G4ThreeVector(0.0,-5000.0*volt/m,0.0 ));

to 

  fField = new G4ElectroMagneticField(gbl_Bfield_x, gbl_Bfield_y, gbl_Bfieldz, gbl_Efield_x, gbl_Efield_y, gbl_Efield_z);
  fLocalField = new G4ElectroMagneticField(lcl_Bfield_x, lcl_Bfield_y, lcl_Bfieldz, lcl_Efield_x, lcl_Efield_y, lcl_Efield_z ));

Thank you for your advice.

ps. In reality I would prefer to apply both a local electric field and a local magnetic field to different subsets of my geometry. Is this really not possible in Geant4?

Feedback Re: simultaneous B and E fields  Keywords: magnetic field, electric field
by Gumplinger Peter <Gumplinger Peter>,   16 Jul, 2007
Re: Question simultaneous B and E fields (Abigail Bickley)

Let's assume you have a global field:

gbl_Bfield_x = 3.3*tesla;

gbl_Bfield_y = 0.0*tesla;

gbl_Bfield_z = 0.0*tesla;

fieldVector = G4ThreeVector(gbl_Bfield_x,gbl_Bfield_y,gbl_Bfield_z);

fMagneticField = new G4UniformMagField(fieldVector);

GetGlobalFieldManager()->SetDetectorField(fMagneticField);

Now you want in addition an electric field inside a subvolume:

lcl_Efield_x = 0.0*volt/m;

lcl_Efield_y = -5000.0*volt/m

lcl_EField_z = 0.0*volt/m; 

fLocalField = 
new G4ElectroMagneticField(gbl_Bfield_x, gbl_Bfield_y, gbl_Bfield_z,
                           lcl_Efield_x, lcl_Efield_y, lcl_Efield_z);

fLocalFieldManager = new G4FieldManager();

fLocalFieldManager->SetDetectorField(fLocalField );

and assign this to the local LogicalVolume as usual. 

> ps. In reality I would prefer to apply both a local electric field

> and a local magnetic field to different subsets of my geometry.

> Is this really not possible in Geant4?

There must be a misunderstanding. You can define as many local fields as you need and assign each to its respective logical volume. If there is no field in the space between them, then don't assign a field to the global field manager. You have to take care that in each subset of your geometry you define a field that is already the superposition of all fields in that space. What you must not do is overlap positioned logical volumes each assigned their respective field then hope that G4 somehow adds the fields in the overlap region.

Question Re: simultaneous B and E fields  Keywords: magnetic field, electric field
by Abigail Bickley <bickley@nscl.msu.edu>,   16 Jul, 2007
Re: Feedback Re: simultaneous B and E fields (Gumplinger Peter)

Thank you for your suggestion. If I understand correctly I still need to have a global field defined as a G4UniformMagField and then the local field is a superposition of the global B and local E fields.

Following your code I added the following lines to my DetectorConstruction.cc. 

  // Global Magnetic Field
  G4double gbl_Bfield_x = 0.0*tesla;
  G4double gbl_Bfield_y = 2.0*tesla;
  G4double gbl_Bfield_z = 0.0*tesla;
  G4ThreeVector globalfieldVector = G4ThreeVector(gbl_Bfield_x,gbl_Bfield_y,gbl_Bfield_z);
  G4UniformMagField* fGlobalField = new G4UniformMagField(globalfieldVector);
  G4FieldManager* fGlobalFieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager();
  fGlobalFieldManager->SetDetectorField(fGlobalField);
  fGlobalFieldManager->CreateChordFinder(fGlobalField);
  fGlobalFieldManager->GetChordFinder()->SetDeltaChord(0.01*mm);

  // Local Electric Field
  G4double lcl_Efield_x = 0.0*volt/cm;
  G4double lcl_Efield_y = -5000.0*volt/cm;
  G4double lcl_Efield_z = 0.0*volt/cm;
  G4ElectroMagneticField* fLocalField = new G4ElectroMagneticField(gbl_Bfield_x,gbl_Bfield_y,gbl_Bfield_z,lcl_Efield_x,lcl_Efield_y,lcl_Efield_z);
  G4FieldManager* fLocalFieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager();
  fLocalFieldManager->SetDetectorField(fLocalField);

My observation is that the global magnetic field is applied properly, but the local electric field causes the code to fail to compile. The error messages I receive seem to indicate that G4ElectroMagneticField is expecting an initialization with two doubles (not 6)? Do you have any suggestions for a way around this problem? 

Compiling TPCDetectorConstruction.cc ...
src/TPCDetectorConstruction.cc: In member function `G4VPhysicalVolume* 
   TPCDetectorConstruction::ConstructTPC()':
src/TPCDetectorConstruction.cc:385: error: cannot allocate an object of type `
   G4ElectroMagneticField'
src/TPCDetectorConstruction.cc:385: error:   because the following virtual 
   functions are abstract:
/opt/geant4.8.2/include/G4ElectroMagneticField.hh:75: error:    virtual void 
   G4ElectroMagneticField::GetFieldValue(const G4double*, G4double*) const
/opt/geant4.8.2/include/G4ElectroMagneticField.hh:79: error:    virtual G4bool 
   G4ElectroMagneticField::DoesFieldChangeEnergy() const
make: *** [/projects/proj2/EOS_sim/TPC/tmp/Linux-g++-3.3/tpc/TPCDetectorConstruction.o] Error 1

Thank you

Warning Re: simultaneous B and E fields  Keywords: magnetic field, electric field
by Gumplinger Peter <Gumplinger Peter>,   17 Jul, 2007
Re: Question Re: simultaneous B and E fields (Abigail Bickley)

Sorry, I gave the wrong advise! This constructor does NOT exist: 

> new G4ElectroMagneticField(gbl_Bfield_x, gbl_Bfield_y, gbl_Bfield_z,
                           lcl_Efield_x, lcl_Efield_y, lcl_Efield_z);

Maybe, it should, to make life easy for situations with both a constant B and E field; or a class G4UniformElectroMagneticField.

Unfortunately, there is no example in /examples/extended/field that would show you how to go about this. The only example I found (but I didn't check/follow the code in detail) that might give you hints of what to do is:

/examples/advanced/microbeam/

where the MicrobeamEMField inherits from G4ElectroMagneticField. You now can write a constructor for 'your field class' with the above signature to initialize both constant B and E fields.

Also note, you MUST NOT get a new local field manager from: 

  G4FieldManager* fLocalFieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager();

and then set: 

  fLocalFieldManager->SetDetectorField(fLocalField);

This will just assign the fLocalField to the World (GlobalFieldManager which is what the G4TransporationManager returns) and hence replace the fGlobalField.

What you need to do is actually make a new FieldManager with:

fLocalFieldManager = new G4FieldManager();

Note, that the examples/field/field03 has this correct only since the G4.9 release.

Question Re: simultaneous B and E fields  Keywords: magnetic field, electric field
by Abigail Bickley <bickley@nscl.msu.edu>,   19 Jul, 2007
Re: Warning Re: simultaneous B and E fields (Gumplinger Peter)

Thank you for your suggestion. Following the example provided by advanced/microbeam/ I was able to create a new class that inherits from G4ElectroMagneticField and have confirmed that it does simultaneously apply an E and B field to the same geometrical volume. However, I have encountered one problem for which I have a followup question. If I apply the E&B field to the whole detector in my DetectorConstruction as follows everything works fine. I see the curvature of charged particle tracks and e- drift in the field gradient. 

      Field = new TPCEMField();
      pEquation = new G4EqMagElectricField(Field);
      pStepper = new G4ClassicalRK4 (pEquation);
      pFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
      pIntgrDriver = new G4MagInt_Driver(0.000001*mm,pStepper,pStepper->GetNumberOfVariables() );
      pChordFinder = new G4ChordFinder(pIntgrDriver);
      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->GetChordFinder()->SetDeltaChord(1e-5*m);
      pFieldMgr->SetDetectorField(Field);

But, I would rather apply it only to a localized volume. You mentioned that this should not be done with a global field manager so instead I changed the code to the following as shown in field03. However, now I find that the field is not being applied anywhere (ie no curved tracks and no e- drift). What am I doing wrong? 

      Field = new TPCEMField();
      pEquation = new G4EqMagElectricField(Field);
      pStepper = new G4ClassicalRK4 (pEquation);
      pFieldMgr = new G4FieldManager();
      pIntgrDriver = new G4MagInt_Driver(0.000001*mm,pStepper,pStepper->GetNumberOfVariables() );
      pChordFinder = new G4ChordFinder(pIntgrDriver);
      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->GetChordFinder()->SetDeltaChord(1e-5*m);
      detector_log->SetFieldManager(pFieldMgr, true);

Thank you for your help.

Feedback Re: simultaneous B and E fields  Keywords: magnetic field, electric field
by Gumplinger Peter <Gumplinger Peter>,   20 Jul, 2007
Re: Question Re: simultaneous B and E fields (Abigail Bickley)

Abigail,

One line is missing in the second code block. You still have to: 

 pFieldMgr->SetDetectorField(Field);

then 

 detector_log->SetFieldManager(pFieldMgr, true);

Peter

Warning Re: simultaneous B and E fields  by Gumplinger Peter <Gumplinger Peter>,   Sep 17, 19:50
Re: Feedback Re: simultaneous B and E fields (Gumplinger Peter)
We also need for EM-fields: 

   G4int nvar = 8;  // To integrate time & energy 
                    //    in addition to position, momentum

   fStepper = new G4ClassicalRK4( fEquation, nvar );

Question Magnetic Field from field01 Example  Keywords: magnetic field ield01
by Matthew Middione <middionematt@csufresno.edu>,   13 Jul, 2007 
It is probabbly apparent by my multiple magnetic field postings on the 
forum that I am having difficulties getting my magnetic field to work.

The basic run down of my situation is as follows.  I am using code from 
TestEm1 (Geant/examples/extended/electromagnetic/TestEm1) and exGPS 
(Geant/examples/extended/eventgenerator/exGPS).  I am using the 
GeneralParticleSource to create positron emitting point sources and 
examining their the distance that these positrons travel prior to 
annihilation with electrons in a tissue like material.   I also want to 
utilize a uniform transverse magnetic field to see how the field alters
the distance the positrons travel.  The distance should be reduced, 
because the magnetic field will cause the charged particles to spiral 
inwards due to the transverse nature of the field. 

For the magnetic field code I have used the F01FieldSetup.cc/hh 
and F01FieldMessenger.cc/hh found in example field01 
(Geant/examples/extended/field/field01).

Everything compiles, links, and runs beautifully.  The problem is that 
I do not believe that the magnetic field is taking any effect.  I know 
I have introduced a magnetic field because the simulation is hard coded 
to print out the magnetic field paramaters.  These are the paramaters 
that were printed:

PreInit> /control/execute Mono.in

F01FieldSetup: magnetic field set to Uniform( 0, 0, 9.4 ) 

EpsilonStep: set min= 1e-05 max= 0.0001

G4ClassicalRK4 (default) is called

The minimal step is equal to 0.01 mm


The Box is 2 m   of Tissue

PhysicsList::SetCuts:CutLength : 1 mm 

To simplyify a test of the magnetic field's effect, I use the following 
macro file to shoot monoentergetic positrons:

/pr/det/setMat Tissue
/pr/det/setSize 2 m
/run/initialize
#
/process/eLoss/verbose 1
/tracking/verbose 1
/pr/stepMax 6 mm
/gps/particle e+
/gps/ene/type Mono
/gps/ene/mono 4 MeV
/gps/direction 1 0 0
/run/beamOn 1 

I do this with both the field set to 0 Tesla and then 9.4 Tesla and I 
change nothing else.  The outputs are identical.  The track length, 
ending position in x-y-z, and range of the particles are the same.  
This indicated that the effect of the magnetic field is not being 
realized.  

I intend on shooting ions using the gps commands, but it is easier to 
test the magnetic field via positrons for now.

Do you have any suggestions or assistance that can help in this matter.
I have been working on getting this magnetic field to operate properly
for weeks.  I sincerely appreciate any input you may provide.

Thanks,
Matthew Middione
Feedback Re: Magnetic Field from field01 Example  Keywords: magnetic field ield01
by Tatiana Nikitina <Tatiana Nikitina>,   27 Sep, 2007
Re: Question Magnetic Field from field01 Example (Matthew Middione)
 Hello,

 To check if the magnetic field is On, 
you can use visualisation macro 'vis.mac' from TestEm1 example and
particle called 'chargedgeantino'. This particle has no physics, 
it is used to check transportation in magnetic field. 
By using this particle with magnetic field equal to 0 and 
with magnetic field not zero, you can easily see
 if trajectory of particle is changed.

If you don't see any difference,
 you can look in examples/novice/N04/ExN04 to see how to create magnetic field
or look in Geant4 tutorials about magnetic field :
  
http://geant4.web.cern.ch/geant4/support/training.shtml

Tatiana Nikitina
None Placing a magnetic field inside of a volume  Keywords: magnetic field, logical volume, testem1
by Matthew Middione <middionematt@csufresno.edu>,   11 Jul, 2007 
I need to place a magnetic field in a volume.  I am a Geant novice so a step by step
approach would be much appreciated.  

I want to apply a local uniform transverse magnetic field in my logical volume, but
I do not know how to go about doing so.  I would like to control its strength via a
detector messenger and I already have that taken care of.  I just need to know how to 
apply a field to a logical volume.  I have tried to follow some of the other postings
regarding this topic on the HyperNews forum, but when I do so the results are dismal 
and it doesn't work.

Thanks,
Matthew Middione
Feedback Re: Placing a magnetic field inside of a volume  Keywords: magnetic field, logical volume, testem1
by Gumplinger Peter <Gumplinger Peter>,   27 Sep, 2007
Re: None Placing a magnetic field inside of a volume (Matthew Middione)

Please, follow examples/extended/field/field03 of version G4.8.2 or G4.9.0.

None Help setting up an Electric Field  by Adam <ax_blais@laurentian.ca>,   10 Jul, 2007

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.

Feedback Re: Help setting up an Electric Field  by Gumplinger Peter <Gumplinger Peter>,   10 Jul, 2007
Re: None Help setting up an Electric Field (Adam)

> 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.

My guess is the segmentation fault has nothing to do with your field setup. You probably have code like this somewhere: 

  G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
  G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
  ...... thePostPV->GetName();

You must check that the pointer thePostPV is not NULL: 

  if(thePostPV)thePostPV->GetName();

It is NULL for the last step when the particle leaves the world volume. 

> I know that to isolate the field to only the xenon cylinder, I need to

> use the line:

> xenon_log->SetFieldManager(field->GetGlobalFieldManager(), false);

No, that's not correct! While yes, you have to set a field manager to the logical volume, if you set the 'global field manager' to it, you achieved nothing. Moreover, you assigned your field to the global field manager, thus your field is everywhere. The 'false/true' means: 

       // Sets FieldManager and propagates it:

       //  i) only to daughters with G4FieldManager = 0

       //     if forceToAllDaughters=false

       // ii) to all daughters

       //     if forceToAllDaughters=true

I know there was an error in example:

examples/extended/field/field03/src/F03DetectorConstruction.cc

until G4.9. It was supposed to have exhibited tracking in a magnetic field associated to selected logical volumes (besides a global field). Now corrected, you want to follow this example, except you don't want to set any field to the global field manager. The code: 

  fLocalFieldManager = new G4FieldManager();
  fLocalFieldManager->SetDetectorField(....
  fLocalFieldManager->SetChordFinder(

  // Set local field manager and local field

  xenon_log->SetFieldManager(fLocalFieldManager,true);

> in what class do I put this line? The documentation isn't very

> specific about this.

In DetectorConstruction.cc where you have the pointer xenon_log.

None Re: Help setting up an Electric Field  by Adam <ax_blais@laurentian.ca>,   10 Jul, 2007
Re: Feedback Re: Help setting up an Electric Field (Gumplinger Peter)

The segmentation fault is gone. How do I confine the field to inside the cylinder without having it outside?

None Re: Help setting up an Electric Field  by John Apostolakis <John Apostolakis>,   10 Jul, 2007
Re: None Re: Help setting up an Electric Field (Adam)
Dear Adam,

In an application in which you need a local electric field, you should 
have two field managers:
  - a global field manager, which holds no field
  - a local field manager, which holds the electric field and the 
classes which define the equation of motion for charged particles.

If you have copied examples/extended/field/field02 you will need to 
modify it so that it does not use the GlobalFieldManager, but instead 
creates only a localFieldManager.

You can see this in the field03 example too - but it suffices simply to 
the calls to GetGlobalFieldManager() in the Setup method and instead add 
a localFieldManager object in the setup method and create such an object.

After that it does not matter in what class you set the logical volume's 
field manager - so long as it has access to both of  xenon_log and the 
field manager.  But a logical place (no pun intended) is to use your 
detector Construction.

Best regards,
John.

> 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);
>
>   

Although I would change
> Is this right? If it is right, in what class do I put this line? The
> documentation isn't very specific about this.


Adam wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/124/1/1"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
> The segmentation fault is gone. How do I confine the field to inside the
> cylinder without having it outside?
>
None Re: Help setting up an Electric Field  by Adam <ax_blais@laurentian.ca>,   11 Jul, 2007
Re: None Re: Help setting up an Electric Field (John Apostolakis)

Alright, thank you both for you assistance. I believe that I have set the field up correctly. When I run the simulation, an electron will remain in the xenon, but keep going and going and going... even when it's energy is down to 10^-8 keV. I set the E field to 1kV/cm. The cylinder is about 2.4m high and 2.4m in diameter. It won't stop tracking it! I seem to be stuck because of this. Any ideas?

None Re: Help setting up an Electric Field  by John Apostolakis <John Apostolakis>,   13 Jul, 2007
Re: None Re: Help setting up an Electric Field (Adam)
Could you please explain a bit more what you are seeing ?  It could be 
an unforeseen side-effect from revisions.

If you can produce some verbose output for that event (eg using 
/tracking/verbose 1), it will be very helpful.

Best regards,
John Apostolakis.

Adam wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/124/1/1/1/1"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
> Alright, thank you both for you assistance. I believe that I have set
> the field up correctly. When I run the simulation, an electron will
> remain in the xenon, but keep going and going and going... even when
> it's energy is down to 10^-8 keV. I set the E field to 1kV/cm. The
> cylinder is about 2.4m high and 2.4m in diameter. It won't stop tracking
> it! I seem to be stuck because of this. Any ideas?
>
None Re: Help setting up an Electric Field  by Adam <ax_blais@laurentian.ca>,   13 Jul, 2007
Re: None Re: Help setting up an Electric Field (John Apostolakis)

Thank you for your help. Peter Gumplinger was right. I imposed a minimum kinetic energy of 1.0e-5keV and the behavior has returned to normal.

Warning Re: Help setting up an Electric Field  by Gumplinger Peter <Gumplinger Peter>,   11 Jul, 2007
Re: None Re: Help setting up an Electric Field (Adam)

Hi Adam,

Now, this is interesting .... my guess is ......

This may be a feature (you could call it a 'bug') of G4. You see that in G4 every track, once on the G4 tracking stack, will be tracked down to the end of its range; i.e. as far as it can go. Now, if the track (charge) remains inside an electric field, it will never come to rest as it is always accelerated.

What you need is some "user limits"; please see:

http://geant4.in2p3.fr/2007/prog/MakotoAsai/Kernel2.pdf

and examples where you find this used:

fMinEkine/G4UserSpecialCuts

Peter

Question Magnetic Field doesn't work when point source positioined at origin  Keywords: GPS, Magnetic Field, Point Source at Origin, TestEm1, exgps
by Matthew Middione <middionematt@csufresno.edu>,   02 Jul, 2007

I am using the following two examples in Geant4 to model my simulation: TestEm1 and exgps

I have effectively created a hybrid model of the two to allow TestEm1 to use the GPS as opposed ot the ParticleGun. I did this so I could model a positron emitting point source.

I have a transverse magnetic field applied in the simulation as well and its strength can be changed. When I run the simulation with a 0.0 T magnetic field the calculated range of a 4 MeV positron is 2.0205 g/cm2. Then I apply a 9.4 T magnetic field expecting to see a significant reduction in the range of the same positron and the resulting range is 2.0279 g/cm2.

If I run the normal TestEm1 which does not utilize the GPS and instead uses the Particle Gun the results are what I expect. However, if I use the same example and change the location of the positron being emitted to the origin as opposed to the face of the detector, which was previously used, the results are consistent with what I first described above in the hybrid example. Basically there is the magnetic field takes no effect. What is causing this problem?

Any help would be much appreciated...

Thanks in anticipation, Matthew Middione

Warning endless loop in magnetic field  by Anton <Anton>,   16 Jun, 2007

Hello,

I found problem of charged particles propagation in magnetic field. In test program I define only world volume with magnetic field and electron/positron as primary particle. This simple setup leads to endless loop. (I use geant4 8.3)

I try different electron energies (200 eV, 200 keV), different world volume materials (vacuum, aluminum), different directions of magnetic field - in all this cases I obtain endless loop.

You can find short program, which reproduce it here: http://cern.ch/begemot/testb.zip

None Re: endless loop in magnetic field  Keywords: magnetic field ield01
by Tatiana Nikitina <Tatiana Nikitina>,   27 Sep, 2007
Re: Warning endless loop in magnetic field (Anton)
 Hello,
 
I tried your test program, it is working correctly (no endless loop).

I have tried different energies and different particles. 
In all cases the program gives expected results and no endless loop.

You can try examples/novice/N01/ExN01 without changes and
Compare the output with  original: examples/novice/N01/exampleN01.out.
Then add the changes step by step to see were it is going wrong. 

Tatiana Nikitina
None Re[2]: endless loop in magnetic field  by Anton <Anton>,   31 Oct, 2007
Re: None Re: endless loop in magnetic field (Tatiana Nikitina)
Hello Tatiana,

  I find that this problem raise on cygwin + g++ (G4SYSTEM=WIN32-g++)
  platform. In the same time all does fine on linux + g++.
  Actually, endless loop not appear in any program with magnetic field
  and charged particles, examples\novice\N02\ works fine with magnetic
  field turned on, but for some cases (combinations of geometry/field)
  I can see the problem.

  Although WIN32-g++ not officially supported now, may be you can help
  me?

> I tried your test program, it is working correctly (no endless loop).
> I have tried different energies and different particles.
> In all cases the program gives expected results and no endless loop.
> You can try examples/novice/N01/ExN01 without changes and
> Compare the output with  original: examples/novice/N01/exampleN01.out.
> Then add the changes step by step to see were it is going wrong. 
> Tatiana Nikitina

-- 
Best regards,
 Anton                            mailto:tosha_korneev@tut.by
None Problem with local magnetic field  by alex <caraboy@gmail.com>,   11 Mar, 2007

Hello,

I am trying to add a local magnetic field into a box, and I encounter several compiler errors that I can`t fix. I managed to add a global uniform magnetic field but it is not what I need, so after searching the forum and studying some .ppt I found at some workshops I managed to add the code, but I get the following errors:

--------------- root@ubuntu:~/GeantProjects/teste/N01# make

Making dependency for file src/ExN01DetectorConstruction.cc ...

Compiling ExN01DetectorConstruction.cc ...

src/ExN01DetectorConstruction.cc: In member function virtual G4VPhysicalVolume* ExN01DetectorConstruction::Construct():

src/ExN01DetectorConstruction.cc:72: error: expected type-specifier before G4ClassicalRK4

src/ExN01DetectorConstruction.cc:72: error: cannot convert int* to G4MagIntegratorStepper* in initialization

src/ExN01DetectorConstruction.cc:72: error: expected , or ; before G4ClassicalRK4

src/ExN01DetectorConstruction.cc:70: warning: unused variable myEquation 

make: *** [/home/oem/GeantProjects/tmp/Linux-g++/testN01/ExN01DetectorConstruction.o] Error 1
------------------

Can you help me figure out what am I doing wrong? I have the following detectorconstruction source:

-------------------- #include "ExN01DetectorConstruction.hh" 

#include "G4Material.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "globals.hh"
#include "G4NistManager.hh"


#include <G4UniformMagField.hh>
#include <G4FieldManager.hh>
#include <G4TransportationManager.hh>


#include <G4Mag_UsualEqRhs.hh>
#include <G4MagIntegratorStepper.hh>
#include <G4ChordFinder.hh>


ExN01DetectorConstruction::ExN01DetectorConstruction()
 :  experimentalHall_log(0), tracker_log(0),
    calorimeterBlock_log(0), calorimeterLayer_log(0),
    experimentalHall_phys(0), calorimeterLayer_phys(0),
    calorimeterBlock_phys(0), tracker_phys(0)
{;}

ExN01DetectorConstruction::~ExN01DetectorConstruction() { }

G4VPhysicalVolume* ExN01DetectorConstruction::Construct() { 

  //------------------------------------------------------ materiale - definire materiale

  G4double a;  
  G4double z; 
  G4double density;

  G4Material* Pb = new G4Material("Lead", z= 82., a= 207.19*g/mole, density= 11.35*g/cm3); // Plumb sau Fier

  G4NistManager* man = G4NistManager::Instance();
  G4Material* Air  = man->FindOrBuildMaterial("G4_AIR"); // Aer
  G4cout <<Air<<"\r";                               
 G4cout << *(G4Material::GetMaterialTable()); // lista materiale

  //------------------------------------------------------ volumes

  //------------------------------ experimental hall (world volume)
  //------------------------------ beam line along x axis

  G4double expHall_x = 6.0*m;
  G4double expHall_y = 3.0*m;
  G4double expHall_z = 3.0*m;
  G4Box* experimentalHall_box = new G4Box("World_Box_Shape",expHall_x,expHall_y,expHall_z); // dedinesc shape pt World Volume
  experimentalHall_log = new G4LogicalVolume(experimentalHall_box, Air,"World_Box__Logical",0,0,0); // definesc proprietati pentru World Volume
  experimentalHall_phys = new G4PVPlacement(0,G4ThreeVector(), experimentalHall_log,"World_Box",0,false,0); //pozitionez in scena World Volume

  G4double block_x = 1.0*m;
  G4double block_y = 50.0*cm;
  G4double block_z = 50.0*cm;

  G4double fieldValue = 10.0*tesla;
  G4UniformMagField* myField =new G4UniformMagField(G4ThreeVector(0.0,fieldValue,0.0));

  G4Mag_UsualEqRhs* myEquation = new G4Mag_UsualEqRhs(myField);

  G4MagIntegratorStepper* myStepper = new G4ClassicalRK4(myEquation);

  G4ChordFinder* myChordFinder = new G4ChordFinder(myField,0.01*mm,myStepper);

  G4FieldManager* fieldMgr = new G4FieldManager(myField,myChordFinder);

  G4Box* calorimeterBlock_box = new G4Box("calorimetru_box", block_x, block_y, block_z);
  calorimeterBlock_log = new G4LogicalVolume(calorimeterBlock_box,Pb,"calorimetru_log", fieldMgr,0,0); //camp mg

  G4double blockPos_x = 1.0*m;
  G4double blockPos_y = 0.0*m;
  G4double blockPos_z = 0.0*m;

  calorimeterBlock_phys = new G4PVPlacement(0,G4ThreeVector(blockPos_x,blockPos_y,blockPos_z),calorimeterBlock_log,"calorimetru",experimentalHall_log,false,0);

return experimentalHall_phys; }

------------------------

Thank you for you time.

None Re: Problem with local magnetic field  by Vladimir IVANTCHENKO <vnivanch@mail.cern.ch>,   11 Mar, 2007
Re: None Problem with local magnetic field (alex)
On Sun, 11 Mar 2007, alex wrote:

> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/118"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
> 
> Hello,
> 
> I am trying to add a local magnetic field into a box, and I encounter
> several compiler errors that I can`t fix. I managed to add a global
> uniform magnetic field but it is not what I need, so after searching the
> forum and studying some .ppt I found at some workshops I managed to add
> the code, but I get the following errors:
> 
> --------------- root@ubuntu :~/GeantProjects/teste/N01# make
> 
> Making dependency for file src/ExN01DetectorConstruction.cc ...
> 
> Compiling ExN01DetectorConstruction.cc ...
> 
> src/ExN01DetectorConstruction.cc: In member function virtual
> G4VPhysicalVolume* ExN01DetectorConstruction::Construct():
> 
> src/ExN01DetectorConstruction.cc:72: error: expected type-specifier
> before G4ClassicalRK4
> 
> src/ExN01DetectorConstruction.cc:72: error: cannot convert int* to
> G4MagIntegratorStepper* in initialization
> 
> src/ExN01DetectorConstruction.cc:72: error: expected , or ; before
> G4ClassicalRK4
> 
> src/ExN01DetectorConstruction.cc:70: warning: unused variable myEquation
> 
> make: *** [/home/oem/GeantProjects/tmp/Linux-g++/testN01/ExN01DetectorConstruction.o] Error 1
> 
> ------------------
> 
> Can you help me figure out what am I doing wrong? I have the following
> detectorconstruction source:
> 
> -------------------- #include "ExN01DetectorConstruction.hh"
> 
> #include "G4Material.hh"
> #include "G4Box.hh"
> #include "G4Tubs.hh"
> #include "G4LogicalVolume.hh"
> #include "G4ThreeVector.hh"
> #include "G4PVPlacement.hh"
> #include "globals.hh"
> #include "G4NistManager.hh"
> 
> #include <G4UniformMagField.hh>
> #include <G4FieldManager.hh>
> #include <G4TransportationManager.hh>
> 
> #include <G4Mag_UsualEqRhs.hh>
> #include <G4MagIntegratorStepper.hh>
> #include <G4ChordFinder.hh>
> 
> ExN01DetectorConstruction::ExN01DetectorConstruction()
>  :  experimentalHall_log(0), tracker_log(0),
>     calorimeterBlock_log(0), calorimeterLayer_log(0),
>     experimentalHall_phys(0), calorimeterLayer_phys(0),
>     calorimeterBlock_phys(0), tracker_phys(0)
> {;}
> 
> ExN01DetectorConstruction::~ExN01DetectorConstruction() { }
> 
> G4VPhysicalVolume* ExN01DetectorConstruction::Construct() {
> 
>   //------------------------------------------------------ materiale - definire materiale
> 
>   G4double a;  
>   G4double z; 
>   G4double density;
> 
>   G4Material* Pb = new G4Material("Lead", z= 82., a= 207.19*g/mole, density= 11.35*g/cm3); // Plumb sau Fier
> 
>   G4NistManager* man = G4NistManager::Instance();
>   G4Material* Air  = man->FindOrBuildMaterial("G4_AIR"); // Aer
>   G4cout <<Air<<"\r";                               
>  G4cout << *(G4Material::GetMaterialTable()); // lista materiale
> 
>   //------------------------------------------------------ volumes
> 
>   //------------------------------ experimental hall (world volume)
>   //------------------------------ beam line along x axis
> 
>   G4double expHall_x = 6.0*m;
>   G4double expHall_y = 3.0*m;
>   G4double expHall_z = 3.0*m;
>   G4Box* experimentalHall_box = new G4Box("World_Box_Shape",expHall_x,expHall_y,expHall_z); // dedinesc shape pt World Volume
>   experimentalHall_log = new G4LogicalVolume(experimentalHall_box, Air,"World_Box__Logical",0,0,0); // definesc proprietati pentru World Volume
>   experimentalHall_phys = new G4PVPlacement(0,G4ThreeVector(), experimentalHall_log,"World_Box",0,false,0); //pozitionez in scena World Volume
> 
>   G4double block_x = 1.0*m;
>   G4double block_y = 50.0*cm;
>   G4double block_z = 50.0*cm;
> 
>   G4double fieldValue = 10.0*tesla;
>   G4UniformMagField* myField =new G4UniformMagField(G4ThreeVector(0.0,fieldValue,0.0));
> 
>   G4Mag_UsualEqRhs* myEquation = new G4Mag_UsualEqRhs(myField);
> 
>   G4MagIntegratorStepper* myStepper = new G4ClassicalRK4(myEquation);
> 
>   G4ChordFinder* myChordFinder = new G4ChordFinder(myField,0.01*mm,myStepper);
> 
>   G4FieldManager* fieldMgr = new G4FieldManager(myField,myChordFinder);
> 
>   G4Box* calorimeterBlock_box = new G4Box("calorimetru_box", block_x, block_y, block_z);
>   calorimeterBlock_log = new G4LogicalVolume(calorimeterBlock_box,Pb,"calorimetru_log", fieldMgr,0,0); //camp mg
> 
>   G4double blockPos_x = 1.0*m;
>   G4double blockPos_y = 0.0*m;
>   G4double blockPos_z = 0.0*m;
> 
>   calorimeterBlock_phys = new G4PVPlacement(0,G4ThreeVector(blockPos_x,blockPos_y,blockPos_z),calorimeterBlock_log,"calorimetru",experimentalHall_log,false,0);
> 
> return experimentalHall_phys; }
> 
> ------------------------
> 
> Thank you for you time.
> 

Hello,

There are extended examples showing how to implement magnetic field.  
Please, have a look:

$G4INSTALL/examples/extended/field

VI
None Re: Problem with local magnetic field  by alex <caraboy@gmail.com>,   17 Mar, 2007
Re: None Re: Problem with local magnetic field (Vladimir IVANTCHENKO )

Thanks, I managed to solve this problem after studying those examples.

Question multiple field maps  Keywords: multiple field maps magnetic field
by Kazutaka Nakahara <nakahara@post.kek.jp>,   27 Feb, 2007 
Hi,

I'd like to use multiple field maps in my simulation.  Does the G4FieldManager
allow this?

I basically have:

----------------------------------------------
     G4FieldManager   *pFieldMgr;
//	 G4MagIntegratorStepper *pStepper;

      //Field grid in A9.TABLE. File must be in accessible from run urn directory
      G4MagneticField* MuonMagField;
      G4MagneticField* MuonMagField1;
      G4MagneticField* MuonMagField2;

      MuonMagField= new MuonField3D("front_field_original.Table", zOffset);
      MuonMagField1= new MuonField3D("transport_file1.Table", zOffset); 
      MuonMagField2= new MuonField3D("transport_file2.Table", zOffset); 

      pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      pFieldMgr1=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      pFieldMgr2=G4TransportationManager::GetTransportationManager()->GetFieldManager();    
  
      G4ChordFinder *pChordFinder = new G4ChordFinder(MuonMagField);
      G4ChordFinder *pChordFinder1 = new G4ChordFinder(MuonMagField1);
      G4ChordFinder *pChordFinder2 = new G4ChordFinder(MuonMagField2);

      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->SetChordFinder( pChordFinder1 );
      pFieldMgr->SetChordFinder( pChordFinder2 );

      pFieldMgr->SetDetectorField(MuonMagField);
      pFieldMgr->SetDetectorField(MuonMagField1);
      pFieldMgr->SetDetectorField(MuonMagField2);
-------------------------------------------------------------

But this doesn't seem to work at all (none of the fields register in the 
simulation).  However, this does work when I use only one field map.

Is it possible to make G4FieldManager accept multiple maps, or do I need to 
combine them into one file?

Thanks in advance,
Kaz
Warning Re: multiple field maps  Keywords: multiple field maps magnetic field
by Gumplinger Peter <Gumplinger Peter>,   30 Apr, 2007
Re: Question multiple field maps (Kazutaka Nakahara)

Hi Kaz,

I am sorry that nobody replied to your posting until now. I am not the most expert on this subject, but I do know that:

G4TransportationManager::GetTransportationManager()->GetFieldManager();

returns a pointer to a G4FieldManager:

G4TransportationManager::G4TransportationManager(){ . . fFieldManager = new G4FieldManager(); . .} 

G4FieldManager* G4TransportationManager::GetFieldManager() const
{
   return fFieldManager;
}

So, the G4 kernel only ever has one such G4FieldManager. If you request
it:
      pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      pFieldMgr1=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      pFieldMgr2=G4TransportationManager::GetTransportationManager()->GetFieldManager();

you get the same pointer every time. And yes, subsequent calls to: 

      pFieldMgr->SetDetectorField(MuonMagField);
      pFieldMgr->SetDetectorField(MuonMagField1);
      pFieldMgr->SetDetectorField(MuonMagField2);

simply replace the previous; i.e. the G4FieldManager doesn't keep a list of G4Fields. It only has one pointer:

G4Field* fDetectorField;

Maybe, if there is need we should extend this to allow the G4FieldManager to keep a list of fields and each subsequent call to SetDetectorField will add automatically to the list.

For now, as I said and unless I am completely wrong with my C++, user code should never have multiple lines:

G4TransportationManager::GetTransportationManager()->GetFieldManager();

in the same method.

None Re: multiple field maps  Keywords: multiple field maps magnetic field
by Tom Roberts <Tom Roberts>,   01 May, 2007
Re: Warning Re: multiple field maps (Gumplinger Peter)

 > Maybe, if there is need we should extend this to allow the 
 > G4FieldManager to keep a list of fields and each subsequent call
 > to SetDetectorField will add automatically to the list.

It takes more than a mere list -- EM fields overlap in general, and they must be summed. But in most cases the evaluation of the field dominates the CPU time to track particles, and summing a bunch of fields only makes that worse.

What I did is implement a GlobalField class (derived from G4ElectroMagneticField) that permits other objects to register themselves as defining a field and a bounding box. Now the GlobalField::GetFieldValue() can efficiently loop over all registered objects, checking their bounding boxes before calling their AddFieldValue() to add its result into the sum. It is important that the bounding box be in global coordinates, because the global->local transform is often the most expensive part of the field computation.

For one of my simulations with 38 individual field objects, implementing these bounding boxes sped it up by a factor of 8 (!).

None Re: multiple field maps  Keywords: multiple field maps magnetic field
by maurizio <ungaro@jlab.org>,   01 Aug, 2007
Re: None Re: multiple field maps (Tom Roberts)

Dear Tom,

could you please provide a working example for what you described? I believe this should be one of GEANT4 official examples.

thanks, mauri

Question insert Electric fileld to logical volume  Keywords: insert Electric fileld to logical volume
by Lee, Hyun-Seung <sherufa@naver.com>,   20 Feb, 2007

With reference guide, I tried to make an uniform electric field and attach that to a logical volume(detector) However, when it is building, error is occured.

("singleGap_log" is the logical volume of detector)

--------***code***---------------------------

G4UniformElectricField* fUniformElectricField = new G4UniformElectricField(G4ThreeVector(0.,0.,100.0*kilovolt/cm));

singleGap_log->SetFieldManager(fUniformElectricField->GetLocalFieldManager(),allLocal);

--------***error message***-----------------------------

src/MyDetectorConstruction.cc: In member function `virtual G4VPhysicalVolume* MyDetectorConstruction::Construct()':

src/MyDetectorConstruction.cc:366: error: 'class G4UniformElectricField' has no member named 'GetLocalFieldManager' ---------------------------------------------------------

below, there is the text of "Geant4 User's Guide For Application Developers Detector Definition and Response" (4.3.2 Practical Aspects)

---------***text of user`s guide***------------------------------------

Creating a Field for a Part of the Volume Hierarchy

It is possible to create a field for a part of the detector. In particular it can describe the field (with pointer fEmField, for example) inside a logical volume and all its daughters. This can be done by simply creating a G4FieldManager and attaching it to a logical volume (with pointer, logicVolumeWithField, for example) or set of logical volumes. 

  G4bool allLocal = true;       
  logicVolumeWithField->SetFieldManager(fEmField->GetLocalFieldManager(),allLocal);

---------------------------------------------

what class is "fEmField" is instance of? I can`t find "GetLocalFieldManager()"method in reference guide..

None Re: insert Electric fileld to logical volume  Keywords: insert Electric fileld to logical volume
by Dan Shapira <Dan Shapira>,   22 Feb, 2008
Re: Question insert Electric fileld to logical volume (Lee, Hyun-Seung)
That was a nice and comprehensive response. Should add that the different
classes such as fLocalEquation, fLocalChord, etc must be defined either in 
header or but declaring them in line - for example:
fLocalEquation =  new G4EqMagElectricField(MyEField);
should be replaced by:
G4EqMagElectricField* fLocalEquation = 
                                 new G4EqMagElectricField(MyEField);
etc....
Warning Re: insert Electric fileld to logical volume  Keywords: insert Electric fileld to logical volume
by Gumplinger Peter <Gumplinger Peter>,   30 Apr, 2007
Re: Question insert Electric fileld to logical volume (Lee, Hyun-Seung)

Again, sorry nobody addressed your problem until now.

The application manual is wrong in stating that you could use: 

  G4bool allLocal = true;       
  logicVolumeWithField->SetFieldManager(fEmField->GetLocalFieldManager(),allLocal);

This line is taken from an example code:

http://www-geant4.kek.jp/lxr/source//examples/extended/field/field03/include/F03FieldSetup.hh#L69

Here, fEmField is actually not a field but of class F03FieldSetup and the code in examples/extended/field/field03 is: 

   logicRadiator->SetFieldManager( fEmFieldSetup->GetLocalFieldManager(), 
                                   allLocal ) ;

Instead, what you need to do is: 

  fLocalEquation = new G4EqMagElectricField(fEMfield);

  G4int nvar = 8;
  fLocalStepper = new G4ClassicalRK4( fLocalEquation, nvar );

  G4double fMinStep = 0.010*mm;
  fLocalChordFinder = new G4ChordFinder(fEmField,
                                        fMinStep, fLocalStepper);

  fLocalFieldManager = new G4FieldManager();
  fLocalFieldManager->SetDetectorField(fEmField);
  fLocalFieldManager->SetChordFinder(fLocalChordFinder);

  G4bool allLocal = true ;
  logicVolumeWithField->SetFieldManager(fLocalFieldManager,allLocal);

None how to add non-uniform M-field and rf EM field  Keywords: non-uniform field, rf field,
by long <long>,   08 Jan, 2007 

Dear Friends, I want to tranport eletron beam through a non-uniform magnetic 
field, then pass through a cavity filled with rf field. How to use these two
fiels type? I have not found them in the manual. Could anybody give some details?
   Thanks a lot!
Feedback Re: how to add non-uniform M-field and rf EM field  Keywords: non-uniform field, rf field,
by John Apostolakis <John Apostolakis>,   09 Jan, 2007
Re: None how to add non-uniform M-field and rf EM field (long)

I would suggest an alternative approach to "transport an electron beam through a non-uniform magnetic field, then pass through a cavity filled with [an] rf field." 

There are a few conditions that would need to be met, for this 
solution to work:
 - the two fields must cover different volumes in space, and these must be well
separated,
 - there must be no conflicts/overlaps between these 'field' volumes and 
'real' volumes which contain parts of the setup

If these conditions are fulfilled it should be possible to create two volumes and put attach a different fields to each one. There would be no problem if one volume contained the other, or if they were apart.

Please remember that each field must use global coordinates only - both for the position and for the field it returns.

For the volume with the magnetic field, the standard way to create the field manager and chord finder can be used.

As Peter mentions, for the RF field, a time-dependent field will be needed. Depending on the field frequency and particle momentum/velocity, such a field could be challenging for the built-in integration. It could be necessary, or just better, to integrate the RF field yourself and give a kick to the particle.

Best regards, John Apostolakis

None Re: how to add non-uniform M-field and rf EM field  Keywords: non-uniform field, rf field,
by Tom Roberts <Tom Roberts>,   09 Jan, 2007
Re: Feedback Re: how to add non-uniform M-field and rf EM field (John Apostolakis)

I have successfully implemented overlaping time-varying electric and magnetic fields, such as an RF cavity inside a solenoid. Email me if you want the code I use. I just use G4EqMagElectricField and G4ClassicalRK4. 

Key point: if an electric field is present anywhere, you must call
    fieldMgr->SetFieldChangesEnergy(true);

For multiple fields generated by multiple components (e.g. magnets and RF cavities), the GlobalField object I implemented permits the individual components to register their field object and a bounding box; to compute the field at a given point it sums all of the individual fields for registered objects with bounding boxes containing the point. In most cases the conversion from global to loacl coordinates is more expensive than computing the local field, so the bounding boxes are in global coordinates.

None Re: how to add non-uniform M-field and rf EM field  Keywords: non-uniform field, rf field,
by long <long>,   09 Jan, 2007
Re: None Re: how to add non-uniform M-field and rf EM field (Tom Roberts)

Thank you Peter,John and Tom! Now I got direction to do it! 

 With best regards!

 Long
Feedback Re: how to add non-uniform M-field and rf EM field  Keywords: non-uniform field, rf field,
by Gumplinger Peter <Gumplinger Peter>,   08 Jan, 2007
Re: None how to add non-uniform M-field and rf EM field (long)

Please, look at postings 39, 41, 55, 63, 72 and 74 on this forum. Basically, you have to supply the field equation. If it depends on time, you need to get it from G4Track.

None How to view helical trajectories in a solenoid  Keywords: helical trajectories, visualization, magnetic field, solenoid
by Ken Teh <teh@anl.gov>,   05 Dec, 2006 
I am unable to "see" my helical trajectories in a solenoid with a
constant B field unless I put in cylindrical slices of vacuum in my
solenoid.  With these slices in place, the trajectory crosses several
"volumes" of vacuum and I can see the helical trajectory.

This is rather hokey.  Is there a way of drawing the trajectories 
without resorting to these slices of vacuum?

Thanks! 

Ken
None Re: How to view helical trajectories in a solenoid  by Anton <Anton>,   09 Dec, 2006
Re: None How to view helical trajectories in a solenoid (Ken Teh)
Hello, Ken,
use this code in DetectorConstruction class to make more smooth visulisation:

G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();
tmanager->GetPropagatorInField()->SetLargestAcceptableStep(1*mm);

-- 
Best regards,
 Anton                            mailto:tosha_korneev@tut.by
None Rotating a field & obtaining the EM field in a certain point  by Giovanni Marchiori <giovanni.marchiori@pi.infn.it>,   26 Sep, 2006 
I have two unrelated questions on electromagnetic fields:
1) when I associate a field to a logical volume with SetFieldManager, when
I then rotate the logical volume and place it inside the world, the field
is not rotated at all: is this behaviour expected? Do I have to rotate the
field by hand?
What I do is:

  G4double fieldValue = 1.5*tesla;
  G4UniformMagField* intField = new G4UniformMagField(G4ThreeVector(0.,0.,fieldValue));
  G4Mag_UsualEqRhs* iEquation = new G4Mag_UsualEqRhs(intField);
  G4MagIntegratorStepper* iStepper = new G4ClassicalRK4(iEquation);
  G4ChordFinder* iChordFinder = new G4ChordFinder(intField,1.0e-6*mm,iStepper);
  G4FieldManager* fieldMgr = new G4FieldManager(intField,iChordFinder);  

  G4RotationMatrix* RotPos = new G4RotationMatrix(
                                                  G4ThreeVector( 0., 0. , 1 ),
                                                  G4ThreeVector( 0., 1. , 0.),
                                                  G4ThreeVector(-1., 0. , 0. )
						  );

  G4Tubs* solidInnerSol = new G4Tubs("sInnerSolenoid",
				     0.,
				     rmin,
				     length/2.0,
				     startphi,
				     deltaphi);
 
  G4LogicalVolume* logicInnerSol = new G4LogicalVolume(solidInnerSol,
						       CentralVacuum,
						       "lInnerSolenoid",
						       fieldMgr);
                    
  G4VPhysicalVolume* physiInnerSol = new G4PVPlacement(RotPos,
						       G4ThreeVector(),
						       logicInnerSol,
						       "pInnerSolenoid",
						       logicWorld,
						       false,
						       0);

  What I obtain is a cilinder with axis along x, but the field inside it
  is still parallel to z. Isn't there any connection between the
  orientation of a volume and that of its associated field?

2) What is the easiest way to know, once I have setup the
geometry of the fields, what is the total electromagnetic field - in global
coordinates - in a generic point in space? Does Geant4 provide such a
functionality?

I am using Geant4 8.1.patch01 on a Linux machine.

Thanks
   Giovanni
Feedback Re: Rotating a field & obtaining the EM field in a certain point  by John Apostolakis <John Apostolakis>,   27 Sep, 2006
Re: None Rotating a field & obtaining the EM field in a certain point (Giovanni Marchiori)

Dear Giovanni,

Thank for your good questions. 

> 1) ... when I [] rotate the logical volume and place it inside the world, the field
> is not rotated at all: is this behaviour expected?

It is indeed the expected (and I think documented) behaviour that the field is expressed in global coordinates, despite being hooked into the logical volume. 

I emphasize that both position and field are in global coordinates:
  - the position you give to the field must be in global coordinates, and
  - the value of the field returned will be in global coordinates.

I will try to inspect your provided code to provide feedback soon. 

> 2) What is the easiest way to know, once I have setup the
> geometry of the fields, what is the total electromagnetic field - in 
> global coordinates - in a generic point in space? Does Geant4 
> provide such a functionality?

This is a very good question - and likely should be in our FAQ.

First I note that only a single field exists for any point of space in Geant4 currently. It can be either a global field (the one chosen for all the setup) or a local field (that is set for a logical volume and can be propagated to ones below, and which overrides any global field.) [ Extensions which allow different field descriptions have been created, including BeamTools from D. Elvira at Fermilab. ]

A single method to get the global field value is not currently available. You will need to call several objects to do this currently and tailor it to your current field needs (eg B field only, EM field, custom field.) I will communicate privately a trial version of this, and once it is ironed out post it.

Best regards, John Apostolakis.

None Re: Rotating a field & obtaining the EM field in a certain point  by Giovanni Marchiori <giovanni.marchiori@pi.infn.it>,   27 Sep, 2006
Re: Feedback Re: Rotating a field & obtaining the EM field in a certain point (John Apostolakis)
Dear John,
thanks for answering my questions and clarifying my doubts.

>> 2) What is the easiest way to know, once I have setup the
>> geometry of the fields, what is the total electromagnetic field - in 
>> global coordinates - in a generic point in space? Does Geant4 
>> provide such a functionality?
>
>This is a very good question - and likely should be in our FAQ.
>
>First I note that only a single field exists for any point of space in 
>Geant4 currently. It can be either a global field (the one chosen for 
>all the setup) or a local field (that is set for a logical volume and
>can be propagated to ones below, and which overrides any global field.) 
>[ Extensions which allow different field descriptions have been created, 
>including BeamTools from D. Elvira at Fermilab. ]

Thanks for emphasizing the "override" issue: at first glance I had the 
wrong idea that the local field aws added to the global one (which would
have made my life easier...). So is there a way to add two fields in Geant4?

>A single method to get the global field value is not currently available.
>You will need to call several objects to do this currently and tailor it
>to your current field needs (eg B field only, EM field, custom field.) 
>I will communicate privately a trial version of this, and once it is
>ironed out post it.

Yesterday I was trying with the following code

  G4StepPoint* pt = new G4StepPoint;
  x[0]=0.;
  x[1]=1.45*m;
  x[2]=0.;
  x[3]=0.;
  pt->SetPosition(G4ThreeVector(x[0],x[1],x[2]));
  G4FieldManager* fMgr = pt->GetPhysicalVolume()->GetLogicalVolume()->GetFieldManager();
  if (fMgr) {
    const G4Field* field = fMgr->GetDetectorField();
    if (field)
      field->GetFieldValue(x,B);
    else
      B[0]=B[1]=B[2]=0.0;
  }

but the application crashed. I will now try your recipe and let you know.

Best regards,

   Giovanni
Question Electrons tracking in high magnetic fields with LowEn. el/m physics  Keywords: electron, magnetic field, low-energy el/m physics.
by Ilya <ilya.kraev@fys.kuleuven.ac.be>,   16 May, 2006

Hi everybody,

I am busy with a simulation of 60Co beta-decay (energies up to 320 keV) in a set-up with magnetic field B = 13 T. In my program I used low-energy electromagnetic processes for this purpose. The final spectrum in a detector that I obtained in sumulations has events that are above end-point energy of 60Co (above 320 keV) and number of these events can not be explained by physics (the detector is too thin to stop two gamma-quanta emitted together with a beta-particle and also the previous simulations in a test set-up did not show such a behaviour).

To check whether it is not related to some bugs in my program I switched to the Standard el/m processes. The spectrum that I see is OK! Then I tried to simulate monoenergetic electron (300 keV). With LOW ENERGY el/m processes I again get events with energy above initaial energy of the particle while for the STANDARD el/m physics the spectrum is fine.

In previous simulations in absence of the magnetic field I do not see such effect.

Did somebody encounter such a problem and is there a bug in LOW ENERGY el/m physics when it has to deal with magnetic fields?

Thanks a lot for any help, Ilya.

None Re: Electrons tracking in high magnetic fields with LowEn. el/m physics  Keywords: electron, magnetic field, low-energy el/m physics.
by michel maire <michel maire>,   17 May, 2006
Re: Question Electrons tracking in high magnetic fields with LowEn. el/m physics (Ilya)

 See forum Physics List, item 199

     Michel
Question Problems with Tracking in Vacuum  Keywords: tracking, magnetic field
by Jens Hasper <hasper@ikp.tu-darmstadt.de>,   27 Mar, 2006

Hi, I'm using a magnetic field inside a vacuum box and having problems to get proper particle tracks. I've already tried to set the UserLimits with the methods below, which has already been suggested in several threads:

G4UserLimits* userLimits=new G4UserLimits(1.*mm); expHall_log->SetUserLimits(userLimits);

However, this didn't improve the simulation at all. The tracks look better when using Air instead of Vacuum (I guess because the TrackingCuts are different). Am I doing something wrong with the assignment of the UserLimits or do you have any suggestion what else to do?

Thanks a lot,

Jens

None Re: Problems with Tracking in Vacuum  by John Apostolakis <John Apostolakis>,   28 Mar, 2006
Re: Question Problems with Tracking in Vacuum (Jens Hasper)
Subject: Re: Problems with Tracking in Vacuum

Dear Jens,

I am unclear whether you are interested to visualise the tracks or you 
are having a problem with a physical observable, boundary crossing or 
something else. You also do not specify what typical value of field you 
have and the average momentum.  Could you let us know these and the 
direction of the momentum with respect to the field ? 

I note that you need to set the user limits to all volume in which you 
want to see (visualise) nice tracks.  One possibility is to set the user 
limits to the Region of the detector (potentially the global region) if 
you want it to apply to many volumes at the same time. Also a smaller 
value for the  could be needed depending on the radius of curvature of 
the tracks. In case the pictures are a concern, I note that the 
visualisation is currently approximate.

We can try to provide further assistance, if you could let us know a bit 
more about the difficulty you are experiencing.

Regards,
John Apostolakis

Jens Hasper wrote:

>*** Discussion title: Fields: Magnetic and Otherwise
>Email replies to PublicHyperNews@slac.stanford.edu must include:
>  In-Reply-To: <"/emfields/98"@geant4-hn.slac.stanford.edu>
>  Subject: ...change this to be about your reply.
>
>Hi, I'm using a magnetic field inside a vacuum box and having problems
>to get proper particle tracks. I've already tried to set the UserLimits
>with the methods below, which has already been suggested in several
>threads:
>
>G4UserLimits* userLimits=new G4UserLimits(1.*mm);
>expHall_log->SetUserLimits(userLimits);
>
>However, this didn't improve the simulation at all. The tracks look
>better when using Air instead of Vacuum (I guess because the
>TrackingCuts are different). Am I doing something wrong with the
>assignment of the UserLimits or do you have any suggestion what else to
>do?
>
>Thanks a lot,
>
>Jens
>  
>
None Re: Problems with Tracking in Vacuum  by Jens Hasper <hasper@ikp.tu-darmstadt.de>,   28 Mar, 2006
Re: None Re: Problems with Tracking in Vacuum (John Apostolakis)

Dear John,

sorry for not having specified my needs properly. I want to track electrons in a constant magnetic field. The electrons have an energy of 30 MeV (momentum in z-direction), the field is 0.4 tesla in y-direction.

I wrote the tracks to a file and they look correct. So I guess it's mainly a problem with the visualization. I'm using HepRepFile and WIRED. As I have already mentioned, the visualization looks much better, when using Air instead of Vacuum as material. Setting the UserLimits to smaller values gave no improvement. For testing I am only using one volume in the whole simulation and assigned the UserLimits to its LogicalVolume.

Is there anything I can do to get a proper visualization?

Regards,

Jens

None Re: Problems with Tracking in Vacuum  Keywords: Step limit
by Makoto Asai <Makoto Asai>,   28 Mar, 2006
Re: None Re: Problems with Tracking in Vacuum (Jens Hasper)
Have you defined G4StepLimiter process class in your physics list?
It should be set to all particle types you want to limit the step
length. G4UserLimit is not effective unless G4StepLimiter is set.
None Re: Problems with Tracking in Vacuum  Keywords: Step limit
by Jens Hasper <hasper@ikp.tu-darmstadt.de>,   28 Mar, 2006
Re: None Re: Problems with Tracking in Vacuum (Makoto Asai)

Hello Makoto,

thanks for your response. Actually I hadn't defined the G4StepLimiter. But even after having defined it now, it still does not have any effect on the UserLimits. I used the following commands for every particle manager: 

  G4StepLimiter* StepLimiter=new G4StepLimiter("StepLimiter");

  pmanager->AddProcess(StepLimiter);

Did I do any mistake or could there be anything else I might have forgotten when using the UserLimits?

Best regards,

Jens

None Re: Problems with Tracking in Vacuum  Keywords: Step limit
by Anton <Anton>,   25 May, 2006
Re: None Re: Problems with Tracking in Vacuum (Jens Hasper)
Hello, Jens,

you can use this in DetectorConstruction class to make more smooth visulisation:

G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();
tmanager->GetPropagatorInField()->SetLargestAcceptableStep(1*mm);
Question Electrons bouncing on the volume boundary in pure magnetic field  Keywords: Magnetic field, electrons, ions, tracking
by Ilya <ilya.kraev@fys.kuleuven.ac.be>,   20 Feb, 2006

(GEANT 4.7.0.p1)

Dear GEANTers,

I have encountered a problem with tracking of the electrons in magnetic field in a simple set-up. I have a source of electrons (radioactive 60Co) located in 20 mum thick foil which is made of VACUUM. Mother volume - WORLD (also vacuum) Electron emitted in a random direction hits the boundary of the foil and gets reflected into the opposite direction. It can even leave the volume and then it is reflected again into oposite direction in vacuum at some point (that does not have any specific properties) in the WORLD volume. It continues to loop like this for a huge number of steps and then it finaly leaves the volume and reaches the detector The same happens to the recoil ion.

The energy of electron is typical beta-spectrum with end-point energy of 320 keV. Recoil ion: 60Ni (energy is about several eV). I'm using the low-energy physics processes.

The processes defined for electron: G4MultipleScattering (with FacRange = 0.01) G4LowEnergyIonisation G4LowEnergyBremsstrahlung

For ion: G4MultipleScattering G4ionIonisation

The magnetic field is user-defined: I have a field profile with 13T central field and decreasing to approx. 1T over 20 cm. The field does not have any sharp points in profile. In the field definition I customize following parameters: 

////////////////////////////////////////////////
fEquation = new G4Mag_UsualEqRhs(this);
fMinStep     = 0.01*mm; // minimal step of 1 mm is default
fMaxStep = 500.*mm;
fStepperType = 4 ;      // ClassicalRK4 is default stepper
fFieldManager = G4TransportationManager::GetTransportationManager()
                                      ->GetFieldManager();
G4TransportationManager::GetTransportationManager()->GetPropagatorInField()->SetMaxLoopCount(500);
G4TransportationManager::GetTransportationManager()->GetPropagatorInField()->SetLargestAcceptableStep(fMaxStep);
//////////////////////////////////////////////////

Please can somebody point out the problem (why it bounces like this) and possible solution (it takes ages to simulate few events or sometimes it stucks forever)!!!

Thanks a lot, Ilya.

Below is a part of output of tracking for one event: 

*********************************************************************************************************
* G4Track Information:   Particle = Co60[0.0],   Track ID = 1,   Parent ID = 0
*********************************************************************************************************


Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0   -620 mum -1.18 mm   7.08 mum     0 eV      0 eV      0 fm      0 fm         Foil    initStep
    1   -620 mum -1.18 mm   7.08 mum     0 eV      0 eV      0 fm      0 fm         Foil  RadioactiveDecay
    :----- List of 2ndaries - #SpawnInStep=  5(Rest= 5,Along= 0,Post= 0), #SpawnTotal=  5 ---------------
    :   -620 mum -1.18 mm   7.08 mum 0.959 eV  Ni60[0.0]
    :   -620 mum -1.18 mm   7.08 mum  1.33 MeV     gamma
    :   -620 mum -1.18 mm   7.08 mum  1.17 MeV     gamma
    :   -620 mum -1.18 mm   7.08 mum   263 keV anti_nu_e
    :   -620 mum -1.18 mm   7.08 mum  55.4 keV        e-
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------
Track (trackID 1, parentID 0) is processed with stopping code 2
A new track 0x928f760 (trackID 2, parentID 1) is passed to G4StackManager.
A new track 0x928f808 (trackID 3, parentID 1) is passed to G4StackManager.
A new track 0x928f8b0 (trackID 4, parentID 1) is passed to G4StackManager.
A new track 0x928f610 (trackID 5, parentID 1) is passed to G4StackManager.
A new track 0x9291030 (trackID 6, parentID 1) is passed to G4StackManager.
### pop requested out of 5 stacked tracks.
Selected G4StackedTrack : 0x928fa68 with G4Track 0x9291030 (trackID 6, parentID 1)
Track 0x9291030 (trackID 6, parentID 1) is passed to G4TrackingManager.

*********************************************************************************************************
* G4Track Information:   Particle = e-,   Track ID = 6,   Parent ID = 1
*********************************************************************************************************


Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0   -620 mum -1.18 mm   7.08 mum  55.4 keV     0 eV      0 fm      0 fm         Foil    initStep
    1   -611 mum -1.18 mm     10 mum  55.4 keV     0 eV   9.28 mum  9.28 mum        Foil  Transportation
    2   -671 mum  -1.1 mm     10 mum  55.4 keV     0 eV   49.3 cm   49.3 cm        World  Transportation
    3   -668 mum -1.15 mm    -10 mum  55.4 keV     0 eV   61.2 mum  49.3 cm         Foil  Transportation
    4   -561 mum -1.11 mm    -10 mum  55.4 keV     0 eV   48.9 cm   98.2 cm        World  Transportation
    5   -599 mum -1.07 mm     10 mum  55.4 keV     0 eV   59.4 mum  98.2 cm         Foil  Transportation
    6   -609 mum -1.18 mm     10 mum  55.4 keV     0 eV   48.1 cm   1.46 m         World  Transportation
    7   -566 mum -1.15 mm    -10 mum  55.4 keV     0 eV   57.7 mum  1.46 m          Foil  Transportation
    8   -639 mum -1.07 mm    -10 mum  55.4 keV     0 eV   47.7 cm   1.94 m         World  Transportation
    9   -674 mum -1.11 mm     10 mum  55.4 keV     0 eV   57.4 mum  1.94 m          Foil  Transportation
   10   -559 mum -1.12 mm     10 mum  55.4 keV     0 eV   47.4 cm   2.41 m         World  Transportation
   11   -583 mum -1.08 mm    -10 mum  55.4 keV     0 eV   55.7 mum  2.41 m          Foil  Transportation
   12   -619 mum -1.07 mm    -10 mum  55.4 keV     0 eV   46.8 cm   2.88 m         World  Transportation
   13   -662 mum -1.09 mm     10 mum  55.4 keV     0 eV   54.4 mum  2.88 m          Foil  Transportation
   14   -622 mum -1.18 mm     10 mum  55.4 keV     0 eV   46.8 cm   3.35 m         World  Transportation
   15   -576 mum -1.17 mm    -10 mum  55.4 keV     0 eV   53.6 mum  3.35 m          Foil  Transportation
   16   -667 mum  -1.1 mm    -10 mum  55.4 keV     0 eV   46.4 cm   3.81 m         World  Transportation
   17   -668 mum -1.15 mm     10 mum  55.4 keV     0 eV     53 mum  3.81 m          Foil  Transportation
   18   -578 mum -1.08 mm     10 mum  55.4 keV     0 eV     46 cm   4.27 m         World  Transportation
   19   -624 mum -1.07 mm    -10 mum  55.4 keV     0 eV   52.5 mum  4.27 m          Foil  Transportation
   20   -572 mum -1.16 mm    -10 mum  55.4 keV     0 eV   46.1 cm   4.74 m         World  Transportation
..................
  524    367 mum -1.32 mm     10 mum  55.4 keV     0 eV   46.6 cm    113 m         World  Transportation
  525    372 mum -1.32 mm    -10 mum  55.4 keV     0 eV   20.7 mum   113 m          Foil  Transportation
  526    300 mum -1.38 mm    -10 mum  55.4 keV     0 eV   37.3 cm    113 m         World  Transportation
  527    301 mum -1.38 mm     10 mum  55.4 keV     0 eV   20.4 mum   113 m          Foil  Transportation
  528   21.8 mum   660 mum  19.3 cm   55.4 keV     0 eV   20.7 cm    114 m         World  Transportation
  529   21.6 mum   660 mum  19.3 cm   54.9 keV   434 eV    292 nm    114 m    FrontDeadLayer         msc

HERE IT REACHED DETECTOR............................ (then it goes fine) next part is for ion 

*********************************************************************************************************
* G4Track Information:   Particle = Ni60[0.0],   Track ID = 2,   Parent ID = 1
*********************************************************************************************************


Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0    590 mum   257 mum  5.63 mum  3.47 eV      0 eV      0 fm      0 fm         Foil    initStep
    1    593 mum   260 mum   -10 mum  3.47 eV      0 eV   29.9 mum  29.9 mum        Foil  Transportation
    2    593 mum   260 mum   -10 mum  3.47 eV 2.74e-16 eV   49.9 cm   49.9 cm        World  Transportation
    3    589 mum   256 mum    10 mum  3.47 eV      0 eV   42.4 mum  49.9 cm         Foil  Transportation
    4    588 mum   266 mum    10 mum  3.47 eV 2.73e-16 eV   49.7 cm   99.6 cm        World  Transportation
    5    592 mum   264 mum   -10 mum  3.47 eV      0 eV   41.8 mum  99.6 cm         Foil  Transportation
    6    586 mum   265 mum   -10 mum  3.47 eV 2.72e-16 eV   49.6 cm   1.49 m         World  Transportation
    7    592 mum   264 mum    10 mum  3.47 eV      0 eV   43.5 mum  1.49 m          Foil  Transportation
    8    588 mum   266 mum   797 mum  3.47 eV 2.74e-16 eV     50 cm   1.99 m         World  Transportation
    9    593 mum   261 mum    10 mum  3.47 eV 1.03e-18 eV   1.76 mm   1.99 m         World  Transportation
   10    589 mum   256 mum   -10 mum  3.47 eV      0 eV   44.9 mum  1.99 m          Foil  Transportation
   11    593 mum   261 mum -6.38 mm   3.47 eV 2.74e-16 eV     50 cm   2.49 m         World  Transportation
   12    585 mum   257 mum   -10 mum  3.47 eV 7.86e-18 eV   1.42 cm   2.51 m         World  Transportation
   13    584 mum   264 mum    10 mum  3.47 eV      0 eV   44.6 mum  2.51 m          Foil  Transportation
   14    593 mum   260 mum    10 mum  3.47 eV 2.71e-16 eV   49.4 cm      3 m         World  Transportation
   15    589 mum   256 mum   -10 mum  3.47 eV      0 eV   42.4 mum     3 m          Foil  Transportation
   16    589 mum   266 mum   -10 mum  3.47 eV 2.74e-16 eV   49.9 cm    3.5 m         World  Transportation
   17    594 mum   261 mum    10 mum  3.47 eV      0 eV   44.1 mum   3.5 m          Foil  Transportation
   18    588 mum   256 mum  4.26 mm   3.47 eV 2.74e-16 eV     50 cm      4 m         World  Transportation
   19    583 mum   261 mum    10 mum  3.47 eV 5.38e-18 eV   9.68 mm   4.01 m         World  Transportation
   20    589 mum   267 mum   -10 mum  3.47 eV      0 eV   45.6 mum  4.01 m          Foil  Transportation
   21    584 mum   259 mum -1.79 mm   3.47 eV 2.74e-16 eV     50 cm   4.51 m         World  Transportation
   22    588 mum   267 mum   -10 mum  3.47 eV 2.27e-18 eV   4.03 mm   4.51 m         World  Transportation
   23    593 mum   261 mum    10 mum  3.47 eV      0 eV   45.2 mum  4.51 m          Foil  Transportation
   24    584 mum   265 mum  1.59 mm   3.47 eV 2.74e-16 eV     50 cm   5.01 m         World  Transportation
   25    585 mum   266 mum    10 mum  3.47 eV 1.86e-18 eV   3.41 mm   5.02 m         World  Transportation
   26    591 mum   265 mum   -10 mum  3.47 eV      0 eV   43.1 mum  5.02 m          Foil  Transportation
   27    583 mum   261 mum -1.46 mm   3.47 eV 2.74e-16 eV     50 cm   5.52 m         World  Transportation
..........................
and so on

Note Re: Electrons bouncing on the volume boundary in pure magnetic field  Keywords: Magnetic field, electrons, ions, tracking
by John Apostolakis <John Apostolakis>,   07 Mar, 2006
Re: Question Electrons bouncing on the volume boundary in pure magnetic field (Ilya)

Please find below a few suggestions for your reported issue/problem. I have also a few questions regarding your report.

1. You do not mention what direction the field is in. From the printout it appears to be in the z direction -- and the particles appearing to bounce or loop have no momentum in the z direction. In this case for a strong field it is physical for tracks to loop ( and in a constant field in a closed circular trajectory when neglecting interactions), and only energy loss or multiple scattering causing a substantial change in the track.

The track potentially feels the effect of scattering (the track length in the world, vacuum volume varies). For clearer information about the 'orbit' of the particle, I suggest limiting the maximum step size in the world volume to 10 cm or, and to see other points on the particle's trajectory. A related point: I am intrigued at the closeness between your value for the maximum step ( 50 cm ) and the step size in vacuum ( ~ 49 cm ). Did you choose your simulation run (or its initial values) out of several runs or are you reproducing a problem you discovered ?

If further investigation is needed of the field tracking, it will be good to revise the SteppingVerbose to print also the momentum direction and the position with higher precisions. But first I suggest investigating the physics and materials:

2. I notice that your electron does not lose energy even in the foil - in which it seems to step several microns: 

 Step#  X        Y       Z      KineE   dEStep  StepLeng  TrakLeng  Volume Process
...
 1  -611 mum -1.18 mm   10 mum  55.4 keV  0 eV   9.28 mum  9.28 mum  Foil  Transportation
...
 3  -668 mum -1.15 mm  -10 mum  55.4 keV  0 eV  61.2 mum   49.3 cm   Foil  Transportation

For me this means that there is a problem with either the physics (the list or a process), with the material of the foil in your setup. A problem in an energy-loss process implementation is possible, but not certain.

a. Physics lists: Could you please let me know whether your physics list is a copy (exact or revised) of one from an example -- or if it is your own creation ? In any case, especially if it is revised or created by you, I think that sharing it would help check whether there are any issues related to this causing your difficulties.

b. Foil material: Could you check it and let us know ?

c. If neither a) nor b) is the cause, then probably there is a problem in a process. I will be glad to forward it to our responsibles -- but I will request your assistance to document it in our problem reporting system to help us track its progress and more easily inform you.

3. Other In a first reading your field parameters look ok. One potential exception is your line 

 > fStepperType = 4 ;      // ClassicalRK4 is default stepper

which requires will change the type of stepper only when used together with other code as in our field examples.

I hope these suggestions help. Best regards, John Apostolakis

Question Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1  Keywords: magnetic field, low energy ion, GEANT4.7.0.p1, looping
by Ilya <ilya.kraev@fys.kuleuven.ac.be>,   20 Jan, 2006

Dear GEANTers,

I have a problem with tracking of the low energy ions in the magnetic field. In my program I simulate the beta-decay of 60Co which is at rest. As a secondary particles (decay products) I obtain one electron (energies < 318 keV), two gammas (1173keV, 1332keV) anti-neutrino and recoil ion (60Ni). When the tracking of these decay products starts electron, gammas and anti-neutrino are tracked without problems but with the recoil ion which has an energy typically about 1 eV the program starts looping. With "/tracking/verbose 2" I see that after few (2-3) steps the ion stops to change it's z-coordinate. The process responsible for the tracking is transportation. (part of the output is in the end of the message.) The field value is 13 tesla but it does not work for 0.1 tesla as well. In zero magnetic field everything works fine. In my definition of the magnetic field I apply constraint for the maximum of the loop count = 1000 but it does not help. The 60Ni ion propagates in copper. The stepper that is used for the tracking is G4ClassicalRK4.

Could somebody help how to solve this problem or point out a reason why it happens?

Thanks a lot, Ilya.

Below I give the part of the output: 

*********************************************************************************************************
* G4Track Information:   Particle = Ni60[0.0],   Track ID = 2,   Parent ID = 1
*********************************************************************************************************


Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0  -12.6 mum -1.16 mm      0 fm  0.912 eV      0 eV      0 fm      0 fm         Foil    initStep
    1  -9.97 mum -1.16 mm     -5 mum 0.912 eV      0 eV   11.8 mum  11.8 mum        Foil  Transportation
particle being tracked is...Ni60[0.0]
    2  -9.99 mum -1.16 mm   -5.5 mum 0.912 eV      0 eV   1.13 mum  12.9 mum       World  Transportation
particle being tracked is...Ni60[0.0]
    3    -12 mum -1.16 mm   -7.5 mum 0.912 eV      0 eV   3.48 mum  16.4 mum      Solder  Transportation
particle being tracked is...Ni60[0.0]
    4  -12.1 mum -1.16 mm     -8 mum 0.912 eV      0 eV   1.14 mum  17.5 mum       World  Transportation
particle being tracked is...Ni60[0.0]
    5  -62.9 mum -1.16 mm  -3.01 mm  0.912 eV      0 eV   6.41 cm   6.41 cm   SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
    6   23.4 mum -1.21 mm  -3.01 mm  0.912 eV      0 eV   99.8 cm   1.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
    7   23.4 mum -1.21 mm  -3.01 mm  0.912 eV      0 eV      0 fm   1.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
    8  -2.86 mum -1.13 mm  -3.01 mm  0.912 eV      0 eV   99.8 cm   2.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
    9  -23.8 mum -1.21 mm  -3.01 mm  0.912 eV      0 eV   2.01 mm   2.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
   10  -43.5 mum -1.25 mm  -3.01 mm  0.912 eV      0 eV   99.9 cm   3.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
   11  -35.7 mum  -1.2 mm  -3.01 mm  0.912 eV      0 eV   2.46 mm   3.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
   12   -102 mum -1.15 mm  -3.01 mm  0.912 eV      0 eV   99.6 cm   4.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
   13   -134 mum -1.22 mm  -3.01 mm  0.912 eV      0 eV   2.03 mm   4.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
   14  -57.9 mum -1.19 mm  -3.01 mm  0.912 eV      0 eV   99.9 cm   5.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
   15   -140 mum -1.19 mm  -3.01 mm  0.912 eV      0 eV   1.99 mm   5.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
   16   -163 mum -1.26 mm  -3.01 mm  0.912 eV      0 eV   99.9 cm   6.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
   17   -106 mum  -1.2 mm  -3.01 mm  0.912 eV      0 eV      2 mm   6.06 m    SourceHolderBottom  Transportation
particle being tracked is...Ni60[0.0]
   18  -27.7 mum -1.16 mm  -3.01 mm  0.912 eV      0 eV   99.7 cm   7.06 m         World  Transportation
particle being tracked is...Ni60[0.0]
Feedback Re: Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1  Keywords: magnetic field, low energy ion, GEANT4.7.0.p1, looping
by John Apostolakis <John Apostolakis>,   22 Jan, 2006
Re: Question Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1 (Ilya)

Dear Ilya,

I think that I have an idea what is causing your track to be seen always at z= 3.01 mm. I am confused, though, why the track shows no energy loss, and why it survives to go several meters in these volumes.

My guess is that z= -3.01 mm is the surface between your World and 'Surface Holder' modules. As the tracks spiraling in a helix, it encounters this surface repeatedly.

I cannot guess the direction of the field that you are applying from this track. Could you please let me know ? In case there is a problem could you also show how you assign it to the world or a particular volume ?

If you would like to see that the particle is progress, despite what the print-out appears to say, I suggest putting a step limit that is between 5 mm and 5 cm, globally or in the world volume. I think that you will see better the helix of the particle.

My main concern, though, is that this particle is not losing energy. If either of the volumes that it encounters has material other than vacuum I expect some energy loss, but none is seen.

I expect that this is because the ion does not have a good physics list. What physics list are you using for your program?

Best regards, John Apostolakis

Feedback Re: Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1  Keywords: magnetic field, low energy ion, GEANT4.7.0.p1, looping
by Ilya <ilya.kraev@fys.kuleuven.ac.be>,   25 Jan, 2006
Re: Feedback Re: Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1 (John Apostolakis)

Dear John,

I checked the physics processes defined for the ion and indeed there was a problem in it: I did not define processes for ion. So I've added processes as it's done in example /extended/radioactivedecay/exrdm. There I found an error: while adding losses due to ionisation for ions instead of defining G4ionIonisation there is G4hIonisation defined (file ExRdmPhysicsListEmStandard.cc, line 99). After assigning correct processes for the ion tracking the problem disappeared.

Thanks a lot for help,

Ilya.

None Re: Seg. fault running a simple app with geant4-07-01  Keywords: GEANT4.7.1 Linux SL4 segmentation fault
by Vladimir Ivanchenko <Vladimir Ivanchenko>,   27 Jan, 2006
Re: Feedback Re: Problem with tracking low energy ions in magnetic field in GEANT 4.7.0.p.1 (Ilya)

Hello,

To clear the situation:

For ions with charge>1: He3, alpha, GenericIons in the PhysicsList one can use G4ionIonisation or G4hLowEnergyIonisation.

For other hadrons G4hIonisation or G4hLowEnergyIonisation.

In the case if G4hIonisation is used for GenericIons no crash will expected but ion ranges will be incorrect.

VI

Question Local Electric Field Segmentation Fault  Keywords: local electric field
by <mcgeac00@uvic.ca>,   12 Jan, 2006

Hello, I have been trying to set up a local electric field inside an enclosed volume.

I followed the example for creating a local magnetic field from the geant4 application developer's manual.

It compiles fine but the simulation exhibits strange behaviour.

The gun can shoot particles parallel to the field with no problems but when I try and shoot a beam in any other direction the simulation crashes with a segmantation fault.

Changing the strength of the field changes the behaviour of the program. If it's 250volts/cm it runs fine (although I don't see it affecting the system at all). If the strength is 250 kilovolts/cm the program segfaults again. If the strength is 100250 kilovolts/cm the program just hangs.

Can anyone suggest a solution to this problem?

Below is how I'm creating the local EField

Thanks,

Jason McGeachie

// Electric Field for the mother volume 

  G4ElectricField*	 	motherEField;
  G4FieldManager* 		motherFieldManager;
  G4EqMagElectricField* 	EFieldEq;
  G4MagIntegratorStepper* 	fStepper;
  G4ChordFinder*          	fChordFinder = 0;
  G4MagInt_Driver*		fIntgrDriver;

  motherEField = new G4UniformElectricField(G4ThreeVector(0.,0.,0.));

  EFieldEq = new G4EqMagElectricField(motherEField);

  //motherFieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager();

  motherFieldManager = new G4FieldManager();

  fStepper = new G4ClassicalRK4(EFieldEq,8);

  motherFieldManager->SetDetectorField(motherEField);

  if(fChordFinder) delete fChordFinder;

  fIntgrDriver = new G4MagInt_Driver(0.010*mm, 	// minimum step of 0.01mm
                                     fStepper, 
                                     fStepper->GetNumberOfVariables() );

  fChordFinder = new G4ChordFinder(fIntgrDriver);

  motherFieldManager->SetChordFinder(fChordFinder);

  TPCmotherLogical->SetFieldManager(motherFieldManager,
  				    true);	// do daughter volumes inherit this field?

None Crazy Point in G4Field::GetFieldValue()  by Tom Roberts <Tom Roberts>,   17 Dec, 2005

In Geant4 7.1, my program occasionally crashes with a segment fault in my GetFieldValue(). This is caused by its argument Point[] being an array of 4 NaN-s. That's crazy -- why should the Geant4 tracking ever generate a Point that is four not-a-number?

This is quite rare, occuring a few times out of 400 million events (800 500k-event jobs on a cluster). I can trace it to tracking a pi+ inside an Al beam pipe (much of the rarity is because most particles don't hit the beam pipes). As I am not interested in particles that hit beam pipes, I have a workaround by killing any particle inside a beam pipe. But soon I will be looking at backgrounds from that....

I am using the physics use case LHEP_BIC.

When I compile and run with G4DEBUG=1, the problem does not occur. This implies that Point[] is probably not initialized when this fault happens; that is deep inside Geant4 code. My program seeds the random number generator with the event# and resets RanGauss(), so I can reliably re-run specific events. I put check code in GetFieldValue() to print "NaN" and force a coredump whenever Point[0] != Point[0] (i.e. Point[0] is NaN); the backtrace below did that.

With G4DEBUG=0 all I can get is the stack traceback from the core file (levels 0,14,15 are mine):

(gdb) backtrace 

#0  0x080ad57e in BLGlobalField::GetFieldValue ()
#1  0x08a29768 in G4ClassicalRK4::DumbStepper ()
#2  0x08a2a65f in G4MagErrorStepper::Stepper ()
#3  0x08a2c764 in G4MagInt_Driver::QuickAdvance ()
#4  0x08a28881 in G4ChordFinder::FindNextChord ()
#5  0x08a2826b in G4ChordFinder::AdvanceChordLimited ()
#6  0x08a3a9c4 in G4PropagatorInField::ComputeStep ()
#7  0x08327e17 in G4Transportation::AlongStepGetPhysicalInteractionLength ()
#8  0x082fda88 in G4SteppingManager::DefinePhysicalStepLength ()
#9  0x082fbded in G4SteppingManager::Stepping ()
#10 0x082f027d in G4TrackingManager::ProcessOneTrack ()
#11 0x082e3ee6 in G4EventManager::DoProcessing ()
#12 0x082cd9ee in G4RunManager::DoEventLoop ()
#13 0x082cd4c1 in G4RunManager::BeamOn ()
#14 0x080c239e in BLManager::trackBeam ()
#15 0x0805ee07 in main ()

None Re: Crazy Point in G4Field::GetFieldValue()  by <vnivanch@mail.cern.ch>,   17 Dec, 2005
Re: None Crazy Point in G4Field::GetFieldValue() (Tom Roberts)
On Sat, 17 Dec 2005, Tom Roberts wrote:

> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/90"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
> 
> In Geant4 7.1, my program occasionally crashes with a segment fault in
> my GetFieldValue(). This is caused by its argument Point[] being an
> array of 4 NaN-s. That's crazy -- why should the Geant4 tracking ever
> generate a Point that is four not-a-number?
> 
> This is quite rare, occuring a few times out of 400 million events (800
> 500k-event jobs on a cluster). I can trace it to tracking a pi+ inside
> an Al beam pipe (much of the rarity is because most particles don't hit
> the beam pipes). As I am not interested in particles that hit beam
> pipes, I have a workaround by killing any particle inside a beam pipe.
> But soon I will be looking at backgrounds from that....
> 
> I am using the physics use case LHEP_BIC.
> 
> When I compile and run with G4DEBUG=1, the problem does not occur. This
> implies that Point[] is probably not initialized when this fault
> happens; that is deep inside Geant4 code. My program seeds the random
> number generator with the event# and resets RanGauss(), so I can
> reliably re-run specific events. I put check code in GetFieldValue() to
> print "NaN" and force a coredump whenever Point[0] != Point[0] (i.e.
> Point[0] is NaN); the backtrace below did that.
> 
> With G4DEBUG=0 all I can get is the stack traceback from the core file
> (levels 0,14,15 are mine):
> 
> (gdb) backtrace
> 
> #0  0x080ad57e in BLGlobalField::GetFieldValue ()
> #1  0x08a29768 in G4ClassicalRK4::DumbStepper ()
> #2  0x08a2a65f in G4MagErrorStepper::Stepper ()
> #3  0x08a2c764 in G4MagInt_Driver::QuickAdvance ()
> #4  0x08a28881 in G4ChordFinder::FindNextChord ()
> #5  0x08a2826b in G4ChordFinder::AdvanceChordLimited ()
> #6  0x08a3a9c4 in G4PropagatorInField::ComputeStep ()
> #7  0x08327e17 in G4Transportation::AlongStepGetPhysicalInteractionLength ()
> #8  0x082fda88 in G4SteppingManager::DefinePhysicalStepLength ()
> #9  0x082fbded in G4SteppingManager::Stepping ()
> #10 0x082f027d in G4TrackingManager::ProcessOneTrack ()
> #11 0x082e3ee6 in G4EventManager::DoProcessing ()
> #12 0x082cd9ee in G4RunManager::DoEventLoop ()
> #13 0x082cd4c1 in G4RunManager::BeamOn ()
> #14 0x080c239e in BLManager::trackBeam ()
> #15 0x0805ee07 in main ()
> 
> 

Please, try to use G4 7.1 patch 01 - in this patch the bug was fixed, 
which provided nans with some probability.

VI
None Re: Crazy Point in G4Field::GetFieldValue()  by <test@gmail.com>,   20 Jun, 2007
Re: None Re: Crazy Point in G4Field::GetFieldValue()

test

None Re: Crazy Point in G4Field::GetFieldValue()  by Tom Roberts <Tom Roberts>,   17 Dec, 2005
Re: None Re: Crazy Point in G4Field::GetFieldValue()
Thank you. I downloaded the Geant4 7.1p01 source, built it and my program (no errors), and re-ran my test. Now instead of a segment fault it loops forever on the same event. Here are the last few lines of my customized steppingVerbose print (in UserSteppingAction()): 
=================== Event 167408188 ==================
=========== Event 167408188 Track 1 pi+  Parent 0 =========
 Step   X(mm)   Y(mm)   Z(mm)    T(ns)      Px,Py,Pz (MeV/c)       StepLen This Volume      Process
... 
   46   -53.2    83.0  3856.5    13.55      8.8     20.2    424.0   100.00 vacuum2          StepLimiter
   47   -52.2    85.3  3905.7    13.72      8.7     19.9    424.1    49.32 vacuum2          Transportation
   48   -51.8    86.7  3928.2    13.80     -2.3     -6.0   -114.4    22.53 pipe2            PionPlusInelastic
   49   -51.9    86.1  3922.1    13.83      nan      nan      nan     6.18 pipe2            hIoni
Following this is an endless number of diagnostic prints stating there is a nan in point[] to G4Field::GetFieldValue() -- that is within a step so I could not even kill the track if I tried.

In the prints above, the nan-s came from step->GetTrack()->GetMomentum(), "This Volume" comes from step->GetPreStepPoint(), and "Process" comes from step->GetPostStepPoint(). vacuum2 has material=Vacuum, and pipe2 has material=Al. I'm using the physics use case LHEP_BIC. There is a fringe B field ~0.03 Tesla here.

None Re: Crazy Point in G4Field::GetFieldValue()  by <vnivanch@mail.cern.ch>,   18 Dec, 2005
Re: None Re: Crazy Point in G4Field::GetFieldValue() (Tom Roberts)
On Sat, 17 Dec 2005, Tom Roberts wrote:

> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/90/1/1"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
> 
> Thank you. I downloaded the Geant4 7.1p01 source, built it and my
> program (no errors), and re-ran my test. Now instead of a segment fault
> it loops forever on the same event. Here are the last few lines of my
> customized steppingVerbose print (in UserSteppingAction()):
> 
> =================== Event 167408188 ==================
> =========== Event 167408188 Track 1 pi+  Parent 0 =========
>  Step   X(mm)   Y(mm)   Z(mm)    T(ns)      Px,Py,Pz (MeV/c)       StepLen This Volume      Process
> ... 
>    46   -53.2    83.0  3856.5    13.55      8.8     20.2    424.0   100.00 vacuum2          StepLimiter
>    47   -52.2    85.3  3905.7    13.72      8.7     19.9    424.1    49.32 vacuum2          Transportation
>    48   -51.8    86.7  3928.2    13.80     -2.3     -6.0   -114.4    22.53 pipe2            PionPlusInelastic
>    49   -51.9    86.1  3922.1    13.83      nan      nan      nan     6.18 pipe2            hIoni
> 
> Following this is an endless number of diagnostic prints stating there
> is a nan in point[] to G4Field::GetFieldValue() -- that is within a step
> so I could not even kill the track if I tried.
> 
> In the prints above, the nan-s came from
> step->GetTrack()->GetMomentum(), "This Volume" comes from
> step->GetPreStepPoint(), and "Process" comes from
> step->GetPostStepPoint(). vacuum2 has material=Vacuum, and pipe2 has
> material=Al. I'm using the physics use case LHEP_BIC. There is a fringe
> B field ~0.03 Tesla here.
> 

Hello Tom,

It is likely problem of a some process in vacuum. To localize the problem 
more information is needed. As John explain, a detailed log file may be 
huge. One of the methods to have limited printout is following:

in the User action (Tracking, Stepping, or even Event actions) check the 
event ID and when it will be ID of the probematic event it is 
possible to perform:

G4UImanager::GetUIpointer()->ApplyCommand("/tracking/verbose 6");

Another method is to identify random seed for the event with the infinite
loop and start from this seed next job, of course, with high verbosity
level.
 
VI
Feedback Re: Crazy Point in G4Field::GetFieldValue()  Keywords: process NaN
by John Apostolakis <John Apostolakis>,   18 Dec, 2005
Re: None Re: Crazy Point in G4Field::GetFieldValue() (Tom Roberts)

Dear Tom,

In both cases (problem with point and direction) the cause can be outside the field module, and in particular with a physics process.

I will suggest a 'brute-force' method that almost certain to succeed - if the problem comes from a process and not the field module. [ I think that there is a more reliable way to track it down, but it will take a bit more time for me/us to find it and let you know. ] 

For the event in which this problem occurs, please use 
  /tracking/verbose 6

I suggest to save the output in a file in a file-system that has lots of space and not to copy it directly to the screen (with tee) in order to speed execution. It will work so long as the program continues to print (as you do for the field value) in order to flush the output buffers.

Then please find the lines in which "NaN" or "nan" first appear to identify the culprit process. If you can copy the relevant section of the printout, we could assist in finding the process.

A more advanced way of handling the output would keep only the relevant 1000 lines or so around the NaN, but I leave it to you to experiment on this.

I will also check for our more straight-forward way of identifying such problems, and try to have it posted here too.

Best regards, John Apostolakis

None Re: Crazy Point in G4Field::GetFieldValue()  Keywords: process NaN
by Tom Roberts <Tom Roberts>,   21 Dec, 2005
Re: Feedback Re: Crazy Point in G4Field::GetFieldValue() (John Apostolakis)
Thank you. This is clearly in the hIoni process, and not in the field 
at all (but that's where the initial segfault happened - NaNs don't 
cause faults, but using them to compute an index does).

It is not at all clear to me how the hIoni process is put together, 
and I have no idea what is wrong. Below are relevant parts of the 
logfile (I can send you the whole thing if necessary). I did 
/tracking/verbose 6 and seeded the RNG for the event with the problem.

Some questions:
 A) Why is the ProposedStep(Post Step) larger than for AlongStep?
    Step 48: (6.18 vs 6.17)
 B) Why do the PostStep momentum values not equal the PreStep
    values for the following step? (here within a step these
    seem to be printed for successive processes)

Note the problem occurs for a pi+ going BACKWARDS. PionPlusInelastic 
turned it around in the previous step (#47, from Pz=454 to Pz=-114).
I don't think that is a problem.

This is Geant4 7.1p01. Physics use case LHEP-BIC, with the Decay 
process disabled. The volume pipe2 has material=Al.

Any help would be appreciated.


--- START ---
G4LEDATA=/home/tjrob/g4mice/g4beamline/data/G4EMLOW2.3
G4LEVELGAMMADATA=/home/tjrob/g4mice/g4beamline/data/PhotonEvaporation
g4beamline ver 0.96beta,  built December 18, 2005
               G4INSTALL=geant4.7.1-gcc3.2.1-Coin
               OIVLIBS=-L Coin-2.4.1-gcc3.2.1/lib -lSoXt -lCoin

*************************************************************
 Geant4 version Name: geant4-07-01    (30-June-2005)
                      Copyright : Geant4 Collaboration
                      Reference : NIM A 506 (2003), 250-303
                            WWW : http://cern.ch/geant4
*************************************************************

@(#)Histo-Scope API     Version 5.0   May 22 2001 - Initialized -


--- BIG SKIP ---

hIoni:   tables are built for  proton
      dE/dx and range tables from 100 eV  to 100 TeV in 120 bins.
      Lambda tables from threshold to 100 TeV in 120 bins.
      Scaling relation is used to proton dE/dx and range
      Bether-Bloch model for Escaled > 2 MeV, ICRU49 parametrisation for protons below.
      Step function: finalRange(mm)= 1, dRoverRange= 0.2, integral: 1


--- BIG SKIP TO START OF STEP 48 ---

#Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
   48    -51.8     86.7 3.93e+03        41     8.59     22.5  3.93e+03       pipe2 PionPlusInelastic

 >>DefinePhysicalStepLength (List of proposed StepLengths):
    ++ProposedStep(PostStep ) = 221.8337849509732 : ProcName = PionPlusInelastic (No ForceCondition)
    ++ProposedStep(PostStep ) = 244.2746510767323 : ProcName = LElastic (No ForceCondition)
    ++ProposedStep(PostStep ) =       100 : ProcName = StepLimiter (No ForceCondition)
    ++ProposedStep(PostStep ) = 6.180367282009078 : ProcName = hIoni (No ForceCondition)
    ++ProposedStep(PostStep ) = 1.797693134862316e+308 : ProcName = msc (Forced)
    ++ProposedStep(PostStep ) = 1.797693134862316e+308 : ProcName = Transportation (Forced)
    ++ProposedStep(AlongStep) = 31.68938164595265 : ProcName = hIoni (CandidateForSelection)
    ++ProposedStep(AlongStep) = 6.170823351947385 : ProcName = msc (NotCandidateForSelection)
    ++ProposedStep(AlongStep) = 6.170823351947385 : ProcName = Transportation (CandidateForSelection)


--- SMALL SKIP To PostStepProcess before the problem ---

>>PostStepDoIt (process by process):    Process Name = msc

    ++G4Step Information
      Address of G4Track    : 0xb24e640
      Step Length (mm)      : 6.180367282009078
      Energy Deposit (MeV)  : 4.764824693109617
      -----------------------------------------------------------------------
       StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -51.75781903625004  -51.86570541540852
         Position - y (mm)   :    86.68243848862836   86.07330254835335
         Position - z (mm)   :    3928.199851403338    3922.05321498811
         Global Time (ns)    :    13.80049770665614   13.83292980324902
         Local Time (ns)     :    13.80049770665614   13.83292980324902
         Proper Time (ns)    :    4.309067393714032   4.334130398274898
         Momentum Direct - x : -0.02008590734084448-0.03755507509820165
         Momentum Direct - y : -0.05210216468989706-0.009777664297793955
         Momentum Direct - z :  -0.9984397431797885  -0.999246723094576
         Momentum - x (MeV/c):   -2.302357774540915  -4.016918062662423
         Momentum - y (MeV/c):   -5.972238241897225  -1.045826062809252
         Momentum - z (MeV/c):   -114.4466847383142  -106.8801540286907
         Total Energy (MeV)  :    180.6067469075589   175.8419222144493
         Kinetic Energy (MeV):    41.03674690755896   36.27192221444931
         Velocity (mm/ns)    :    190.2690235975647   182.3570754640252
         Volume Name         :                pipe2               pipe2
         Safety (mm)         :   0.3529215700473277  0.7476569960578692
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :        PostStep Proc       PostStep Proc
         Process defined Step:    PionPlusInelastic               hIoni
      -----------------------------------------------------------------------


    ++G4ParticleChange Information
      -----------------------------------------------
        G4ParticleChange Information
      -----------------------------------------------
        # of 2ndaries       :                    0
      -----------------------------------------------
        Energy Deposit (MeV):                    0
        Track Status        :                Alive
        True Path Length (mm) :                 6.18
        Stepping Control     :                    0

        Position - x (mm)   :                -51.9
        Position - y (mm)   :                 86.1
        Position - z (mm)   :             3.92e+03
        Momentum Direct - x :              -0.0376
        Momentum Direct - y :             -0.00978
        Momentum Direct - z :               -0.999

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0

 >>PostStepDoIt (process by process):    Process Name = hIoni

    ++G4Step Information
      Address of G4Track    : 0xb24e640
      Step Length (mm)      : 6.180367282009078
      Energy Deposit (MeV)  : 4.764824693109617
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -51.75781903625004  -51.86570541540852
         Position - y (mm)   :    86.68243848862836   86.07330254835335
         Position - z (mm)   :    3928.199851403338    3922.05321498811
         Global Time (ns)    :    13.80049770665614   13.83292980324902
         Local Time (ns)     :    13.80049770665614   13.83292980324902
         Proper Time (ns)    :    4.309067393714032   4.334130398274898
         Momentum Direct - x : -0.02008590734084448                 nan
         Momentum Direct - y : -0.05210216468989706                 nan
         Momentum Direct - z :  -0.9984397431797885                 nan
         Momentum - x (MeV/c):   -2.302357774540915                 nan
         Momentum - y (MeV/c):   -5.972238241897225                 nan
         Momentum - z (MeV/c):   -114.4466847383142                 nan
         Total Energy (MeV)  :    180.6067469075589   175.2471911177023
         Kinetic Energy (MeV):    41.03674690755896   35.67719111770232
         Velocity (mm/ns)    :    190.2690235975647   182.3570754640252
         Volume Name         :                pipe2               pipe2
         Safety (mm)         :   0.3529215700473277  0.4578441892248865
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :        PostStep Proc       PostStep Proc
         Process defined Step:    PionPlusInelastic               hIoni
     -----------------------------------------------------------------------


    ++G4ParticleChange Information
      -----------------------------------------------
        G4ParticleChange Information
      -----------------------------------------------
        # of 2ndaries       :                    1
        Pointer to 2ndaries :            0xb24e988
        (Showed only 1st one)
      -----------------------------------------------
        Energy Deposit (MeV):                    0
        Track Status        :                Alive
        True Path Length (mm) :                    0
        Stepping Control     :                    0

        Charge (eplus)   :                    1
        Kinetic Energy (MeV):                 35.7
        Momentum Direct - x :                  nan
        Momentum Direct - y :                  nan
        Momentum Direct - z :                  nan

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 0

 **PostStepDoIt (after all invocations):
    ++List of invoked processes
      1) Transportation
      2) msc
      3) hIoni (Forced)

    ++G4Step Information
      Address of G4Track    : 0xb24e640
      Step Length (mm)      : 6.180367282009078
      Energy Deposit (MeV)  : 4.764824693109617
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :   -51.75781903625004  -51.86570541540852
         Position - y (mm)   :    86.68243848862836   86.07330254835335
         Position - z (mm)   :    3928.199851403338    3922.05321498811
         Global Time (ns)    :    13.80049770665614   13.83292980324902
         Local Time (ns)     :    13.80049770665614   13.83292980324902
         Proper Time (ns)    :    4.309067393714032   4.334130398274898
         Momentum Direct - x : -0.02008590734084448                 nan
         Momentum Direct - y : -0.05210216468989706                 nan
         Momentum Direct - z :  -0.9984397431797885                 nan
         Momentum - x (MeV/c):   -2.302357774540915                 nan
         Momentum - y (MeV/c):   -5.972238241897225                 nan
         Momentum - z (MeV/c):   -114.4466847383142                 nan
         Total Energy (MeV)  :    180.6067469075589   175.2471911177023
         Kinetic Energy (MeV):    41.03674690755896   35.67719111770232
         Velocity (mm/ns)    :    190.2690235975647   181.2984889452028
         Volume Name         :                pipe2               pipe2
         Safety (mm)         :   0.3529215700473277  0.4578441892248865
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :        PostStep Proc       PostStep Proc
         Process defined Step:    PionPlusInelastic               hIoni
      -----------------------------------------------------------------------

    ++List of secondaries generated (x,y,z,kE,t,PID):  No. of secodaries = 5
      [Note]Secondaries from AlongStepDoIt included.
      -51.75781903625004 86.68243848862836 3928.199851403338 43.98470745932611 13.80049770665614            neutron
      -51.75781903625004 86.68243848862836 3928.199851403338 11.67567080606411 13.80049770665614             proton
      -51.75781903625004 86.68243848862836 3928.199851403338 45.85114518000712 13.80049770665614            neutron
      -51.75781903625004 86.68243848862836 3928.199851403338 11.60944364114667 13.80049770665614            neutron
      -51.86570541540852 86.07330254835335 3922.05321498811 0.5947310967469912 13.83292980324902                 e-
None Re: Crazy Point in G4Field::GetFieldValue()  Keywords: process NaN
by Tom Roberts <Tom Roberts>,   03 Jan, 2006
Re: None Re: Crazy Point in G4Field::GetFieldValue() (Tom Roberts)

OK. I discovered an error on my part: I carefully patched Geant4 7.1 and recompiled it, and then linked my program with the old (unpatched) version (:-(). The change in behavior (segment fault => infinite loop) was due to a change I made to protect GetFieldValue() from NaNs in its input parameters.

Using Geant4 7.1p01 changes the events (different multiple scatter angles and ionization losses), and the problem event no longer even hits the beam pipe where the error occured. But I have since run 30 times more than the previous average number of events between errors, with no errors.

So I believe the patch fixed the problem for me, too.

Thank you very much for your assistance. Sorry for the confusion.

Warning Re: Crazy Point in G4Field::GetFieldValue()  Keywords: process NaN
by Gumplinger Peter <Gumplinger Peter>,   04 Jan, 2006
Re: None Re: Crazy Point in G4Field::GetFieldValue() (Tom Roberts)

I don't know what suddenly fixed your problem (or not), but just so everybody knows, I have an outstanding bug report about NaN momentum directions returned from a process:

http://pcitapiww.cern.ch/asdcgi/geant4/problemreport/show_bug.cgi?id=778

In the bug report I suggest code that will flag this.

In our case we found ProtonInelastic and NeutronInelastic to sometimes cause NaNs but it well also be true for PionPlusInelastic.

None Re: Crazy Point in G4Field::GetFieldValue()  Keywords: process NaN
by Tom Roberts <Tom Roberts>,   04 Jan, 2006
Re: Warning Re: Crazy Point in G4Field::GetFieldValue() (Gumplinger Peter)

What fixed my problem was patch01 to Geant4 7.1.

The design changes in Geant4 8.0 make it difficult for me to migrate to it (:-(), so I'll be on 7.1p01 for a while.

Question cold neutron spin tracking  Keywords: spin, neutrons, magnetic field
by Emil Frlez <Emil Frlez>,   22 Nov, 2005 
 Dear GEANT-ers,

 I would appreciate if people who have a better insight
 into the GEANT4 tracking could help me increase my productivity ;8-)

 My problem is the cold neutron spin tracking in the non-homogeneous
 magnetic field. I am using geant4.6.2.p01. Reading the exchange
 between Bertalan and Peter it appeared to me (at the first sight)
 that this problem is (basically) solved in G4.

 So I used the Peter's recipe with "G4Mag_SpinEqRhs.cc". First thing
 I notice is that the anomaly is hard-wired for muons. Second issue,
 the code will only change the spin of the charged particles
 because dSpin is multiplied by particle charge. I checked 
 that code compiles and works for electrons, it does run, and
 does the spin procession, though I have not checked that 
 the results are OK.

 So, I thought that for cold neutrons only thing necessary is
 to copy and rename "G4Mag_SpinEqRhs.cc" into, say,
 "G4Mag_Neutron_SpinEqRhs.cc", change "anomaly = 1.165923e-3;"
 into  "anomaly = 2.913;" and set "ParticleCharge = -1;"
 to reflect the negative neutron magnetic moment and get non-zero
 dSpin ...

 But then I noticed that my "GetFieldValue" is not called at 
 all for neutron, so there will be no tracking for (presumably)
 any neutral particle. Does that means that the G4 FieldManager
 checks the charge of the tracked particle and quits if it is dealing
 with neutral particle? Can this behavior be changed in an elegant
 way? And are the simple changes I was thinking of, outlined above, 
 only ones necessary for the neutron spin tracking?
 (With a rush of cold neutron experiments I find google-ing
 that a number of people are working on the spin tracking
 simulations. If the full machinery already exists in G4
 and an experienced C++ programmer can make the fixes in 10 mins,
 what are those people working on ...;8-?) Or is the whole problem
 more complicated and I am missing something?

                                        Best Regards, Emil Frlez
Question Re: cold neutron spin tracking  Keywords: spin, neutrons, magnetic field
by Emil Frlez <Emil Frlez>,   11 Jan, 2006
Re: Question cold neutron spin tracking (Emil Frlez)
 Thanks Peter,
 
 your comment was helpful, but has enabled me to only
 partially solve the problem. Maybe you or
 others can think about this: 

 I did not want to modify and recompile
 the Geant4 source so that it would track the neutral
 particle in magnetic field, but I fudged it:

 (1) I define the new type of neutron that has
     a charge "1.0*eplus" and call this new partice "bdneutron".
     I am interested only how its spin changes as
     a function of time and position, and I want it
     to move straight, field or no field.

 (2) Therefore, I copy G4Mag_SpinEqRhs.cc into
     my new user file G4Mag_Neutron_SpinEqRhs.cc where I do
     necessary modifications:

     (i)    // Reset the anomaly for Neutron
            anomaly = -2.913;

     (ii)   // "turn off" the force on the "bdneutron" ;8-)
 
           dydx[3] = 0.; //cof*(y[4]*B[2] - y[5]*B[1]) ;   // Ax = a*(Vy*Bz - Vz*By)
           dydx[4] = 0.; //cof*(y[5]*B[0] - y[3]*B[2]) ;   // Ay = a*(Vz*Bx - Vx*Bz)
           dydx[5] = 0.; //cof*(y[3]*B[1] - y[4]*B[0]) ;   // Az = a*(Vx*By - Vy*Bx)

 Now, I thought I would have a bdneutron tracked straight in a magnetic field,
 with spin equations integrated. Well, it actually works, I tried it
 for a polarized neutron with Ekine=1.0001*eV. It steps
 thru my "World" from (0,1799.9*mm,0) to (0,-1999.5*mm,0) and
 spin is tracked properly, I tried several different spin states, eg.:

sx   sy       sz        x y/mm   z particle  Ekine/MeV  time/ns

0.99 0.140759 0.0093242 0 1799.9 0 bdneutron 1.0001e-06 2167.34
0.99 0.139835 0.0186078 0 1799.8 0 bdneutron 1.0001e-06 4334.68
0.99 0.138299 0.0278099 0 1799.7 0 bdneutron 1.0001e-06 6502.02
0.99 0.136158 0.0368904 0 1799.6 0 bdneutron 1.0001e-06 8669.36
0.99 0.133422 0.0458095 0 1799.5 0 bdneutron 1.0001e-06 10836.7
0.99 0.130103 0.0545283 0 1799.4 0 bdneutron 1.0001e-06 13004
0.99 0.126214 0.0630085 0 1799.3 0 bdneutron 1.0001e-06 15171.4
     ...
0.99 0.140802 -0.00865308 0 -1999.4 0 bdneutron 1.0001e-06 8.23459e+07
0.99 0.140802 -0.00865308 0 -1999.5 0 bdneutron 1.0001e-06 8.2348e+07

 But then, in my application, I need to track a cold neutron
 with Ekine=1E-6*eV. And my code does not track this kind of
 "charged" bdneutron. I try Ekine=0.9999*eV and
 the Geant4 then does just one step and then it terminates, 
 and it seems there is also a memory problem, because it thinks
 that the kinetic energy of the particle is 2.11758e-22/MeV:

 0.99 0.140759 0.00932521 0 1799.9 0 bdneutron 2.11758e-22 2167.56
 0.99 0.140759 0.00932521 0 1799.9 0 bdneutron 2.11758e-22 2167.56

 So there is some cut somewhere, I assume, in the integration of equations 
 of motion, there is some "if statement", I am not able to locate 
 with a value Ekine_cut=1*eV. 

 I have checked that the default, neutral neutron with Ekine=9.09495E-7*eV, 
 is transported thru my detector correctly:

 0.99 0.141067 0 0 1799.9 0 neutron 9.09495e-13 7581.02
 0.99 0.141067 0 0 1799.8 0 neutron 9.09495e-13 15162
 0.99 0.141067 0 0 1799.7 0 neutron 9.09495e-13 22743.1
 0.99 0.141067 0 0 1799.6 0 neutron 9.09495e-13 30324.1
 0.99 0.141067 0 0 1799.5 0 neutron 9.09495e-13 37905.1
 0.99 0.141067 0 0 1799.4 0 neutron 9.09495e-13 45486.1
 0.99 0.141067 0 0 1799.3 0 neutron 9.09495e-13 53067.2
 0.99 0.141067 0 0 1799.2 0 neutron 9.09495e-13 60648.2
 0.99 0.141067 0 0 1799.1 0 neutron 9.09495e-13 68229.2
....
0.99 0.141067 0 0 -1999.5 0 neutron 9.09495e-13 2.88041e+08

 But now, the equations of motion are not integrated and
 spin components do not change ;8-( A bummer!

 So, it appears to me that even if I edited out and recompiled
 "if( (particleCharge != 0.0) )" in the code quoted by Peter above in
 G4double G4Transportation::AlongStepGetPhysicalInteractionLength()
 I would still have non-working code for Ekine<1*eV ...

					Best Regards, Emil
Note Re: cold neutron spin tracking  Keywords: spin, neutrons, magnetic field
by Emil Frlez <Emil Frlez>,   03 Feb, 2006
Re: Question Re: cold neutron spin tracking (Emil Frlez)

 Hi all,

 my thanks to Peter and John for helping me out with
 an advice. Here is the summary:

 Peter was right: in order to do a neutral particle
 spin tracking G4Transportation.cc needs to be modified.

 The simplest way I find to accomplish that was:

 (i) copy the geant4.6.2.p01 sources into geant4.6.2.p02
 (ii) in  geant4.6.2.p02 modify the G4Transportation,
      change 

   if( (particleCharge != 0.0) )
   { ... }

      into 

   fNeutron = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
   if( (particleCharge != 0.0 || pParticle->GetDefinition() == fNeutron ) )
   { ... }

 (iii) Modify G4Mag_SpinEqRhs.cc, change 

   dSpin = ParticleCharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));

   into

   dSpin = omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));

   or, for nonrelativistic case:     
   dSpin = (3.8260837/2.) * omegac*(Spin.cross(BField));

   Set the anomaly for neutron, set the straight propagation:

   anomaly = -2.913;
   if ( particleCharge == 0. ) {dydx[3] = 0.;   dydx[4] = 0.;   dydx[5] = 0.;}

 (iv) recompile the sources and set 
      setenv G4INSTALL $HOME/geant4.6.2.p1

 With these changes I get the correct neutron spin precession.

                                               Cheers, Emil

Feedback Re: cold neutron spin tracking  Keywords: spin, neutrons, magnetic field
by Gumplinger Peter <Gumplinger Peter>,   22 Nov, 2005
Re: Question cold neutron spin tracking (Emil Frlez)

Hi Emil

.... and I still haven't written the documentation for spin tracking.

So great, you found the thread on this forum which advertises this functionality in Geant4!

http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/emfields/53/1.html 

> First thing I notice is that the anomaly is hard-wired for muons.

Yes, while the class G4Mag_SpinEqRhs constructor 

G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
   : G4Mag_EqRhs( MagField ) 
 {
    anomaly = 1.165923e-3;
 }

sets the muon anomaly, but no, the class also comes with the public method:

inline void SetAnomaly(G4double a) { anomaly = a; }

with which you can set the anomaly to any value you like. 

> But then I noticed that my "GetFieldValue" is not called at 

> all for neutron, so there will be no tracking for (presumably)

> any neutral particle. Does that means that the G4 FieldManager

> checks the charge of the tracked particle and quits if it is dealing

> with neutral particle?

John Apostolakis is the expert here, but I suspect you are correct in that the code is written such that field tracking is done only for particles with non-zero charge. I believe it is in G4Transportation: 

   G4double G4Transportation::
      AlongStepGetPhysicalInteractionLength(

   G4ParticleDefinition* pParticleDef   = pParticle->GetDefinition() ;
   G4double              particleCharge = pParticleDef->GetPDGCharge() ; 

   G4FieldManager* fieldMgr=0;
   G4bool          fieldExertsForce = false ;
   if( (particleCharge != 0.0) )
   {
      fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); 
      if (fieldMgr != 0) {
         fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
      }
   }

   if( !fieldExertsForce ) 
   {
      G4double linearStepLength ;
      linearStepLength = fLinearNavigator->ComputeStep( .....

etc. 

> Can this behavior be changed in an elegant way?

Good point! We should try hard to extend the whole idea of field tracking to include magnetic moments of neutral particles. John ...??? 

> And are the simple changes I was thinking of, outlined above, only
> ones necessary for the neutron spin tracking?

I think so, and so long as we manage to propagate the neutron straight (zero electric charge) and yet precess its magnetic moment in a magn. field. 

> With a rush of cold neutron experiments I find google-ing

> that a number of people are working on the spin tracking

> simulations. If the full machinery already exists in G4

> and an experienced C++ programmer can make the fixes in 10 mins,

> what are those people working on ...;8-?) Or is the whole problem

> more complicated and I am missing something?

I don't know about 10 mins ... but I agree with you, the problem shouldn't be complicated with what's already in place in G4. If is is, then only because magn. moment tracking wasn't considered in the class design from the outset.

Question Problem with Quadrupole Magnetic Field!!  by S.Fonseca UERJ/Brazil <sfonseca@uerj.br>,   24 Oct, 2005

Hello people,

I'm trying to simulate 4 low beta quadrupole magnets that are being used in the Tevatron beam line. For this I used G4QuadrupoleMagField method for simulate the Quadrupole Magnetic Field. However particles that enter in contact with this field disappear and I get the following message during Event: 

"G4Transportation is killing track that is looping or stuck.
 This track has XXX MeV energy"

best regards,

Sandro Fonseca

Graduate student of the State University of Rio de Janeiro (UERJ/Brazil)

None Re: Problem with Quadrupole Magnetic Field!!  by Tom Roberts <Tom Roberts>,   17 Dec, 2005
Re: Question Problem with Quadrupole Magnetic Field!! (S.Fonseca UERJ/Brazil)

Check your units!

Whenever this has happened to me, I found I had >1000 T fields.

None Re: Problem with Quadrupole Magnetic Field!!  by John Apostolakis <John Apostolakis>,   25 Oct, 2005
Re: Question Problem with Quadrupole Magnetic Field!! (S.Fonseca UERJ/Brazil)

I suggest using the following methods of PropagatorInField to 

  void    SetMaxLoopCount( G4int new_max );

to enable it to go further in case it is simply taking too many substeps in order to integrate in the much-varying field. The current default is 1000. 

  G4int  SetVerboseLevel( G4int verbose );

if you give it a value of 2 or 3 it print some more information which could be useful in identifying the problem. 

To find the G4PropagatorInField, please use the G4TransportationManager. Its method
   static G4TransportationManager* GetTransportationManager()
gives the singleton's pointer and the method
   G4PropagatorInField*  GetPropagatorInField() 
provides that.

If you can provide more information from these I will be glad to try to assist.

Best regards, John Apostolakis

Warning Trajectory in dipole field independent of energy !  by Tim Chambers <Tim Chambers>,   13 Jul, 2005 
So I have my Earth with its dipole working, and running pretty speedily.  
Now comes the big problem: the paths of charged particles in the field 
are independent of energy!  I have gone through the following steps, 
changing no code except the GetFieldValue() function...

First, uniform B-field along z axis.  As expected, the radius of curvature
of the particle's path increases with energy.  As a side note, I found 
that above a certain energy, about 8.5MeV for an electron in a 3 nT field,
the particle would go absolutely straight, zero curvature.
Then, a B-field along the z-axis where field strength is a function of the
particle's z coordinate.  Again, radius of curvature increases with energy.
Now, a more complicated function, looking like
Bfield[0] = y
Bfield[1] = z
Bfield[2] = x

And still the path depends on energy.  Adding another layer of complexity:
Bfield[0] = y*z
Bfield[1] = z*x
Bfield[2] = x*y

And as always, more energy makes the particle take larger loops with less
curvature.  

So then I implemented my dipole field:
Bfield[0] = 3*coeff*z*x
Bfield[1] = 3*coeff*z*y
Bfield[2] = coeff*(3*z*z - r*r)

and suddenly, the trajectory becomes independent of energy.  For the same 
starting position and direction, an electron ranging anywhere from 10 
eV to 10 TeV will take exactly the same path.  The initial direction has 
a small effect; it determines whether the particle travels forward or 
backward along the field line trapping it but otherwise does not matter.
For example, these two run macros:

/gps/position -8000. -8000. -8000. km
/gps/energy 100 eV
/gps/direction 1. 0. 0.

and

/gps/position -8000. -8000. -8000. km
/gps/energy 10 MeV
/gps/direction 1. 0. 1.

produce identical output.  What is the cause of this problem?
None Addendum to above...  by Tim Chambers <Tim Chambers>,   14 Jul, 2005
Re: Warning Trajectory in dipole field independent of energy ! (Tim Chambers)
Since posting the last message I made an obvious check that I overlooked 
earlier: I ran the simulation three times, with an electron, positron, 
and gamma.  The gamma followed a straight line (as expected), but the 
electron and positron trajectories were identical.  So I've managed to 
create a simulation where the path in a dipole field is independent of 
energy and now charge as well.  Advice on how to un-break this would be 
greatly appreciated...

-Tim Chambers
None Re: Trajectory in dipole field ... what is happening ?  by John Apostolakis <John Apostolakis>,   15 Jul, 2005
Re: None Addendum to above... (Tim Chambers)

Your results are very surprising. It is hard to see how the simpler fields work in your setup but the more complex one does not. The field integration sub-module treats all magnetic fields in the same way (using 'standard' Runge-Kutta integration) unless someone chooses otherwise. So the results that you see are very unusual. Likely there is an issue with the code creating the field.

I would suggest a few more checks:

Can you please let us know the name of the Equation of motion that you are using ? ( I would expect that you do not need to set this (that you are using that provided automatically. )

Also could you post (or give a web link) to the code of your field ?

In addition, it would be helpful to see both the 'physics' step reported by 'tracking/verbose 1' and the exact 'sub-steps' taken in the field -- I will try to provide a way to do this.

Best regards, John Apostolakis

None program hang on in magnetic field  by Zhongdong Zhang <zh.zhang@fz-juelich.de>,   11 Jul, 2005

Dear all,

In my magnetic field, the world has global magnetic field zero and a cylinder in world has local magnetic field 3 Tesla along z-axis. When I run my programme, the verbose output as follows: 

21   -1.9 mm   29.9 cm  -67.6 cm   17.1 MeV     0 eV   2.69 cm     74 cm    physiBore    LElastic    :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  1 ---------------
    :   -1.9 mm   29.9 cm  -67.6 cm    966 keV  N14[0.0]
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------
   22  -1.53 cm   28.6 cm    -70 cm   17.1 MeV     0 eV   3.04 cm     77 cm    physiBore  Transportation   23  -15.6 cm   14.8 cm    -95 cm   17.1 MeV     0 eV   31.8 cm   1.09 m    OutOfWorld  Transportation
*********************************************************************************************************
* G4Track Information:   Particle = N14[0.0],   Track ID = 190,   Parent ID = 132
*********************************************************************************************************

    1  -1.73 mm   29.9 cm  -67.6 cm    958 keV  7.53 keV   275 mum   275 mum   physiBore       hIoni    :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  1 ---------------
    :  -1.73 mm   29.9 cm  -67.6 cm    150 eV         e-
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------
    2  -1.31 mm   29.9 cm  -67.6 cm    939 keV  18.5 keV   664 mum   939 mum   physiBore       hIoni    :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  2 ---------------
    :  -1.31 mm   29.9 cm  -67.6 cm    147 eV         e-
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------
    3  -1.21 mm   29.9 cm  -67.6 cm    935 keV  3.85 keV   151 mum  1.09 mm    physiBore       hIoni    :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  3 ---------------
    :  -1.21 mm   29.9 cm  -67.6 cm    147 eV         e-
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------
    4     -1 mm     30 cm  -67.6 cm    927 keV  8.68 keV   323 mum  1.41 mm    physiBore       hIoni    :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  4 ---------------
    :     -1 mm     30 cm  -67.6 cm    145 eV         e-
    :----------------------------------------------------------------- EndOf2ndaries Info ---------------

Then the program seems enter a infinite loop and did nothing (no exit, no killing loop).

Is anything wrong or should I setup some parameter to control the field?

Thank you.

Zhongdong

Idea Re: program hang on in magnetic field  by John Apostolakis <John Apostolakis>,   12 Jul, 2005
Re: None program hang on in magnetic field (Zhongdong Zhang)

Output from /tracking/verbose shows that a process makes all the momentum direction cartesian components 'NaN'.

Looking at further verbose output (from /tracking/verbose 6) that Zhongdong has provided, we see that it is not the transportation (and thus the field module) that is responsible for this NaN in momentum, but an Ionisation process (see verbosity below).

In order to identify the process and the underlying reason, we are looking into the physics lists used. 

Best regards,
John.
------------------------------------------------------------------------------------
 >>PostStepDoIt (process by process):    Process Name = hIoni

    ++G4Step Information
      Address of G4Track    : 0x94509a8
      Step Length (mm)      : 0.1831550063464724
      Energy Deposit (MeV)  : 0.007035918703070614
      -----------------------------------------------------------------------
        StepPoint Information               PreStep            PostStep
      -----------------------------------------------------------------------
         Position - x (mm)   :     845.595933610975   845.6682652455552
         Position - y (mm)   :   -143.0060193537246  -142.8497442884144
         Position - z (mm)   :    150.7255513994749   150.7874230557541
         Global Time (ns)    :    19.26187279975728   19.32340311636678
         Local Time (ns)     :                    0 0.06153031660949893
         Proper Time (ns)    :                    0 0.06152728927959584
         Momentum Direct - x :   0.3952969563502984                 nan
         Momentum Direct - y :   0.8540531128474157                 nan
         Momentum Direct - z :   0.3381325135739386                 nan
         Momentum - x (MeV/c):    51.13571432925922                 nan
         Momentum - y (MeV/c):    110.4805268520178                 nan
         Momentum - z (MeV/c):    43.74090754250312                 nan
         Total Energy (MeV)  :    13040.85002131427   13040.84288593296
         Kinetic Energy (MeV):   0.6416179440887584  0.6344825627746238
         Velocity (mm/ns)    :    2.973826720314346   2.957477625205532
         Volume Name         :                World               World
         Safety (mm)         :                    0   71.33887903644329
         Polarization - x    :                    0                   0
         Polarization - y    :                    0                   0
         Polarization - Z    :                    0                   0
         Weight              :                    1                   1
         Step Status         :            Undefined       PostStep Proc
         Process defined Step:            Undefined               hIoni
      -----------------------------------------------------------------------
None Large system, very slow field tracking... how to improve?  by Tim Chambers <Tim Chambers>,   27 Jun, 2005 
I have built a full-size Earth (~6400km radius) with its ~10^22 J/T 
dipole, and I'm having a bit of trouble with run speed.  I know 
the field is the cause of the slowness because when I set the field =0, 
the simulation runs quickly and with about the results I expect.  
I also tested a small dipole with a 1m sphere and got the usual 
behavior one expects from a dipole magnet, with a short run time.

I currently have the field initialized in the DetectorConstruction class:

	SolarField* theField = new SolarField();
	G4FieldManager* fieldMgr =
		G4TransportationManager::GetTransportationManager()
		->GetFieldManager();
	
	fieldMgr->CreateChordFinder(theField);
	fieldMgr->SetDetectorField(theField);

like so.  I figured there were two things I could do to speed up the 
performance: use a lower-order integrator, and increase step size.  
When I tried to use the following code, the program compiles and runs...
somewhat.  It works fine until I command /run/beamOn at which point it 
says "Run Start" and then locks up.  It doesn't crash, it simply doesn't
do anything after that.  The code:

	SolarField* theField = new SolarField();
	G4FieldManager* fieldMgr =
		G4TransportationManager::GetTransportationManager()
		->GetFieldManager();
	
	G4EqMagElectricField* theEqn = new G4EqMagElectricField(theField);
	G4MagIntegratorStepper* theStepper = new G4SimpleHeum(theEqn);	
	
	G4double MinStep = 10*m;
	G4MagInt_Driver* theDriver = new G4MagInt_Driver(MinStep,theStepper, 
					theStepper->GetNumberOfVariables() );
	G4ChordFinder* theFinder = new G4ChordFinder(theDriver);
	
	fieldMgr->SetChordFinder(theFinder);
	fieldMgr->SetDetectorField(theField);

Can anybody help me figure out why this refuses to run, or give me any 
other advice for speeding up the program?  In the time it's taken me to
write this, one electron has managed to travel a whole centimeter in 
vacuum...

Thanks,
Tim Chambers
None Re: Large system, very slow field tracking... how to improve?  by Tom Roberts <Tom Roberts>,   27 Jun, 2005
Re: None Large system, very slow field tracking... how to improve? (Tim Chambers)
You also need to set the tracking parameters much larger -- the defaults
are appropriate for a typical HEP detector, so the tracking will be 
integrating the equations of motion with ~millimeter parameters. 
See the Geant4 user's guide and physics manual for details.

The default values are (units are mm):
        deltaChord=3.0
 deltaIntersection=0.1
      deltaOneStep=0.01
            epsMax=0.05
            epsMin=2.5e-7
           maxStep=100.0
           minStep=0.01

There may well be others....

Tom Roberts    tjrob@fnal.gov
Ok Re: Large system, very slow field tracking... how to improve?  by Tim Chambers <Tim Chambers>,   29 Jun, 2005
Re: None Re: Large system, very slow field tracking... how to improve? (Tom Roberts)
Just in case anyone else has the same problem I thought I should post 
my solution... I implemented this code yesterday and it works great.  
It would not be an exaggeration to say the program runs more than 100x 
faster.
in DetectorConstruction.cc:


	SolarField* theField = new SolarField();
	G4FieldManager* fieldMgr =
		G4TransportationManager::GetTransportationManager()
		->GetFieldManager();
	G4double minStep = 1*mm;
	
	G4Mag_UsualEqRhs* theEqn = new G4Mag_UsualEqRhs(theField);
	G4MagIntegratorStepper* theStepper = new G4HelixExplicitEuler(theEqn);	
	G4MagInt_Driver* theDriver = new G4MagInt_Driver(minStep, theStepper, 
					theStepper->GetNumberOfVariables());
	G4ChordFinder* theFinder = new G4ChordFinder(theDriver);
	
	fieldMgr->SetChordFinder(theFinder);
	fieldMgr->SetDetectorField(theField);


In addition to using the Helix integrator, scaling the tracking 
parameters to match the simulation size gave me another noticeable 
improvement in the running time.  Thanks all!

-Tim Chambers
None G4MagIntegratorDriver::OneGoodStep error  by John Carter <John Carter>,   13 Jun, 2005 
Hi all,

I'm having problems with a user-defined solenoid field. I have set up the
field and tested it with a few dry runs and it seems to work ok. However, 
some events produce the following 'error' output:

G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper 
  Step's start x=16376 and end x= 16376 are equal !! 
  Due to step-size= 2.5e-13 . Note that input step was 2.5e-06

This is something that randomly happens and I call it an 'error' because it
doesn't cause the program to crash, merely locks it up. I have tried changing 
my MaxAllowedStep for my volume and I don't seem to get the error if I reduce
the step length to something on the order of 10mm. However, this vastly reduces
the speed of the tracking (18 hours to track 1000 events compared to a few seconds...)

Has anyone seen this error and know the cause or a workaround?

Cheers,

John
None Re: G4MagIntegratorDriver::OneGoodStep error  by John Carter <John Carter>,   14 Jun, 2005
Re: None G4MagIntegratorDriver::OneGoodStep error (John Carter)
A reply to myself with the solution for the benefit of the Forum:

After speaking to Geant4 developer we realised that the problem came from 
having a discontinuity in my user-defined field map. So Geant4 was having
significant problems in getting accuracy at certain locations.

The solution was to provide an interpolation of the field between the set
field map points - essentially smoothing out discontinuities.

Regards,

John
None Re: G4MagIntegratorDriver::OneGoodStep error  by Tom Roberts <Tom Roberts>,   21 Jun, 2005
Re: None Re: G4MagIntegratorDriver::OneGoodStep error (John Carter)

For the record, there is another way to cause an enormous number of microscopic steps and huge tracking times (minutes per track): get the units wrong. Tracking in a 2,500 Tesla field takes almost forever!

Question magnetic field in a vacuum box  Keywords: magnetic field in a vacuum box
by <maryam@jlab.org>,   08 Jun, 2005

Dear all, 

 I have a 3 tesla uniform magnetic field in a box of vacuum which is rotated 30 degrees with respect to my beamline. 
The field properly bends the beam but sharply not like a smooth curve.
 But the e- beam also gets bent another 30 degrees sharply, when it exits the volume. 
Could anyone help me figure out why this is happening?

Many thanks,

Maryam

None Re: magnetic field in a vacuum box  Keywords: magnetic field in a vacuum box
by John Carter <John Carter>,   13 Jun, 2005
Re: Question magnetic field in a vacuum box 
Someone may correct me if I'm wrong, but I think you are seeing sharp
turns because of the MaxAllowedStep set for your volume. You'll need to 
set some user limits and apply them to your logical volume for the box.

e.g.

 G4UserLimits* BoxUserLimits = new G4UserLimits();
 BoxUserLimits->SetMaxAllowedStep(10*mm);
 BoxLogicalVolume->SetUserLimits(BoxUserLimits);

Where the length set should be something small enough to be able to see
the curvature as the curve will be made up of a number of straight lines
whose maximum length is set by your MaxAllowedStep.

Hope that helps.
Question Simultaneous E and B Fields???  Keywords: electric magnetic field simultaneous
by Daniel McDevitt <Daniel McDevitt>,   07 Jun, 2005 
I’m trying to create an electric AND magnetic field, simultaneously 
inside a portion of my geometry into which I will be shooting positrons.  
I’m using a global magnetic field (cannibalized from /novice/exampleN02)
and a local electric field.  

  
// For global mag field
    G4double fieldValue = 0.5687*tesla; 
    fpMagField = new ExN02MagneticField(G4ThreeVector (0.,0.,fieldValue));

//  Local Electric Field
     G4double pfieldValue = 0.001*volt/cm;
     G4UniformElectricField* myElectricField = new G4UniformElectricField(G4ThreeVector(pfieldValue,0,0));
     G4EqMagElectricField *myeEquation= new G4EqMagElectricField(myElectricField);
     G4int nvar = 8;
     G4MagIntegratorStepper* myStepper = new G4ClassicalRK4(myeEquation,nvar);
     G4MagInt_Driver* myIntgrDriver=new G4MagInt_Driver(1.0e-3*mm, myStepper, myStepper->GetNumberOfVariables());
     G4ChordFinder* myChordFinder = new G4ChordFinder(myIntgrDriver);
     G4bool fieldChangesEnergy = true;
     G4FieldManager* pFieldMgr = new G4FieldManager(myElectricField,myChordFinder,fieldChangesEnergy);


When the electric field is turned on inside the geometry (ie pFieldMgr 
is inserted into the Logical volume as shown below), the magnetic field
no longer works.  However, the electric field functions.

// Aerogel Hi
 solidAerogelhi1 = new G4Box("Aerogel_hi1Box", aerogelx,aerogely,aerogelz);
 logicAerogelhi1 = new G4LogicalVolume(solidAerogelhi1 , AerogelhiMater, "Aerogel_hi1",0,0,0);
 physiAerogelhi1 = new G4PVPlacement(0,              // no rotation
                                  positionAerogelhi1, // at (x,y,z)
                                  logicAerogelhi1,    // its logical volume
                                  "Aerogel_hi1",  // its name
                                  logicWorld,      // its mother  volume
                                  false,           // no boolean operations
                                  0);              // copy number
// Aerogel hi replicated
   solidAeroghi1Column
     = new G4Box("Aeroghi1ColumnBox",1.0*cm,0.25*cm,0.025*cm);
   logicAeroghi1Column
     = new G4LogicalVolume(solidAeroghi1Column,AerogelhiMater,"logicAeroghi1Column",pFieldMgr,0,0);
  physiAeroghi1Column = new G4PVReplica("physiAeroghi1Column",logicAeroghi1Column,
                  logicAerogelhi1, kZAxis,100,0.05*cm);



It appears that the electric field pre-empts the magnetic field and I 
am left with the options of either one or the other, depending on which
one I do not implement.

How does one create both a magnetic and electric field, simultaneously 
in a geometry?

Thanks,

Daniel
Feedback Re: Simultaneous E and B Fields???  Keywords: electric magnetic field simultaneous
by Gumplinger Peter <Gumplinger Peter>,   07 Jun, 2005
Re: Question Simultaneous E and B Fields??? (Daniel McDevitt)

> When the electric field is turned on inside the geometry (ie pFieldMgr
> is inserted into the Logical volume as shown below), the magnetic field
> no longer works.


> It appears that the electric field pre-empts the magnetic field

Yes, this is so by design. The solution is for you to specify a local electric AND magnetic field by inheriting from G4ElectroMagneticField

Unfortunately, there is not a single example in the G4 distribution showing how this is done but I think it should be easy from the comments in G4ElectroMagneticField.hh:

// class G4ElectroMagneticField

//

// Class description:

//

// A full Electromagnetic field, containing both electric and magnetic fields.

// It is an abstract class, and a derived type of this field must be

// created by the user to describe his/her field configuration.

// We have established a convention for the electromagnetic field components:

// In the GetValue method, the return values of Bfield will have

// the following meaning

// - Components 0, 1 and 2 are the Magnetic Field (x, y, z respectively);

// - Components 3, 4 and 5 are the Electric field (x, y, z respectively).

//

// Note 1: one or the other field could optional, depending on the Equation

// Note 2: such a convention is required between any field and its

// corresponding equation of motion.

On the other hand, I support your idea that maybe there should be the option of having the local field superimposed on the global field; or for that matter, allowing local fields attached to various logical volumes to overlap; i.e. in effect to purposely create overlapping volumes for the purpose of, for example, simulating how the fringe field of one mag. element overlaps with the fringe field of another, where each field exists in its respective volume.

None How to make custom magnetic field class?  by Tim Chambers <Tim Chambers>,   25 May, 2005 
I am building a simulation that has a static magnetic 
field that is a function of the x,y,z coordinates;

B[0] = cos(x);  //just examples, the real functions are uglier
B[1] = y*y;

that sort of thing.

I have gathered that for a non-uniform magnetic field it is necessary to 
create your own magnetic field class.  By looking at the source code 
for G4UniformMagField I get the general idea of how this works, but 
I can't figure out how to make the magnetic field values depend on 
position.  Any advice will be greatly appreciated!

-Tim Chambers
None Trajectories in a Magnetic field.  by Ryan Romero <ryan3523@neo.tamu.edu>,   31 Jul, 2006
Re: None How to make custom magnetic field class? (Tim Chambers)

Hi. I was hoping someone could help me with a few questions about magenetic fields. I think I have one defined, I've included that code below.

Have I propely defined a field that exists within a particular rectangular region and drops off with first power distance from the z axis? Just to keep things simple I only care about distance along a single axis(y).

Supposing I've defined the magnetic field as I like it, how do I incorporate it into my code so that it will affect trajectories? I saw a method AddTrajectory() but I didn't see it implemented in the N02 code. I assume its defined in the source code for GEant4 but I want to avoid altering that if possible.

Thanks for any help.

#include "ExN01Field.hh" 

ExN01Field::ExN01Field()
{
  BzMax = 6.87*tesla;
  rmax_sq = sqr(50.*cm);
  zmax = 450.*cm;
  xmax= 50.*cm;
  ymax= 275.*cm;
  ymin=250.*cm;
}

ExN01Field::~ExN01Field() {;} 

void ExN01Field::GetFieldValue(const double Point[3],double *Bfield) const
{
  Bfield[0] = 0.;
  Bfield[1] = 0.;
  if(std::abs(Point[2])<zmax && std::abs(Point[0])<xmax &&
  std::abs(Point[1]-525.*cm)<ymax )
  { Bfield[2] = BzMax*ymin/Point[1]; }
  else
  { Bfield[2] = 0.; }
}
Feedback Re: How to make custom magnetic field class?  by Gumplinger Peter <Gumplinger Peter>,   25 May, 2005
Re: None How to make custom magnetic field class? (Tim Chambers)

The simplest example is probably:

/examples/novice/N04/src/ExN04Field.cc, line 36

I agree, your simple question should be answered in:

http://wwwasd.web.cern.ch/wwwasd/geant4/G4UsersDocuments/UsersGuides/ForApplicationDeveloper/html/Detector/electroMagneticField.html

None Different radii for the 'same' electron in B-field  Keywords: GEATN4 magnetic field electrons curve
by Valentin <Valentin.Kozlov@fys.kuleuven.ac.be>,   10 May, 2005

Dear GEANT4ers,

I have the following situation: there are 3MeV electrons in a homogenious magnetic field (in vacuum). The field is defined via GetFieldValue() function setting B[0]=0.; B[1]=0.; B[2]=9*tesla; (in 'real' case there is a fieldmap, but for the test the field is described as above). Electrons start from the (0,0,0) point in direction (0,1.,0), i.e. perpendicular to the magnetic field. For this configuration I calculate the radius of e- motion as 1.3mm. To store the "trajectory" of electrons in GEANT4, (x,y,z) position is saved in GetFieldValue() function. Now, when these positions are plotted, one can see different circles of different radii. The major part of events is on the right radius of ~1.3mm but there are more circules from ~1.6mm to 1.1mm radius and some individual events far from the center. Going from fMinStep of 1mm to 0.1mm changed the situation so that the structure of circles became more 'fine', i.e. more visually different circles appeared near the right radius. As I do not see any physics (it has to be only one curve, isn't it?!), it should be some settings of GEANT4 which are used or something else from how it is calculated. Could you tell me, where to look for the solution?

Thank you, Valentin

P.S. GEANT4 version: 6.1

None Re: Different radii for the 'same' electron in B-field  Keywords: GEATN4 magnetic field electrons curve
by John Apostolakis <John Apostolakis>,   10 May, 2005
Re: None Different radii for the 'same' electron in B-field (Valentin)

I think that you are not extracting the points on the trajectory correctly. The points at which GetFieldValue is called are not on the trajectory of the particle. Following the RungeKutta method, the field is sampled at different trial points, coming from linear and higher order approximations of the final and intermediate points of a step. Then these are used to estimate the endpoint of a step, which is accepted if its estimate is deemed accurate enough (see below).

Of course if you used a method that assumed that the field was uniform, then it would only evaluate it at the step's starting point. But since your real problem has a non-uniform field it makes sense to use the general Runge-Kutta based methods.

To get the actual points of the trajectory, you should create and use a G4Trajectory of the endpoints of the steps or another way to store or print only the endpoints of each step. To find out how to create a trajectory please see section 5.1.6 of the Users Guide for Application Developers (UGAD), which you can find at http://cern.ch/geant4/G4UsersDocuments/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/tracking.html

On a related topic, for more information about obtaining the accuracy you require for your application, please see the UGAD, section 4.3 and in particular the part of 4.3.2 on "How to Adjust the Accuracy of Propagation" http://cern.ch/geant4/G4UsersDocuments/UsersGuides/ForApplicationDeveloper/html/Detector/electroMagneticField.html

In particular if you require high accuracy for the endpoints, you should choose small values for the epsilon parameters. In general these should be smaller than the ratio of the error in the radius (or endpoint) that you are willing to tolerate over the total distance travelled by the track (along the helix). I would expect good results with values the order 10^-5, and do not know fields which are measured better than about one part in 10^4 or so anyway. I note that values below 10^-9 or so will be very challenging to achieve and likely not meaningful due to precision issues.

Best regards, John Apostolakis

None