Synchrotron radiation incorrect, affected by unrelated geometry

[sudipbhattarai]

In a very simple synchrotron simulation, I have uniform magnetic field defined over 1m (-50cm to 50cm in the coordinate system). There is a planar geometry at the end with material G4_Galactic in fact every material that I have in this sim is G4_Galactic (world volume as well).

Primary generator: 11GeV electron created at (0,0,-2*m)
Magnetic field: Uniform 1.5 T over 1m.

void MagneticField::GetFieldValue(const G4double point[], G4double* field) const{
    double Bx = 0;
    double x = point[0], y = point[1], z = point[2]; 

    
    if( std::abs(z)  <= 50*cm and std::abs(y) < 5*cm and std::abs(x) < 1*cm ) {
        Bx =  1.5*tesla;
    }

    field[0] = Bx;
    field[1] = 0.0;
    field[2] = 0.0;
}

I then use: /physics_list/em/SyncRadiationAll true.

Not very clear picture but the 20cm axes are at (0,0,0).

I also added /vis/scene/add/magneticField 30 to see the field line, and the field lines are only drawn from z=0 to z=50*cm. This seems like a visual artifact.

But Instead of synchrotron radiation starting from the start of magnetic field region (-50cm), it starts only later through the magnetic field. It seemed like visual artifact because of drawing again. But when I look at origin vertices of the synchrotron photon, it confirms the visual that it starts some way later.

This is the z coordinate of vertices of the synchrotron photon (pid == 22). It says the synchrotron photons only get created around ~z=-320mm Instead of expected z=-500mm.

I was experimenting a bunch of ways. Everything else same, I add a planar G4_Galactic geometry of z=1mm thickness at z= -50*cm.

Now the synchrotron radiation correctly starts at z=-50cm to z=50cm.

I originally bumped into this problem because I am working on another more complex project where I use parallel geometry to use as virtual detectors. And it appeared, the presence and absence of some parallel planes (all G4_Galactic) seemed to affect the synchrotron radiation, completely baffling.

This is a small toy project, here is the entire source code if its helpful.
zipped source code

For more clarity here is my `ConstructSDandField` method.

void Geometry::ConstructSDandField(){
    fMagneticField = new MagneticField();
    fMagneticField->Activate();
}

and the Activate function

void MagneticField::Activate(){
    G4TransportationManager* transportationmanager = G4TransportationManager::GetTransportationManager();
    auto fFieldMgr = transportationmanager->GetFieldManager();

    fFieldMgr->SetDetectorField(this);

    auto dof = 12;
    auto equation = new G4Mag_UsualEqRhs(this);
    auto stepper = new G4ClassicalRK4(equation, dof);

    G4double minStep = 0.01 * mm;
    G4ChordFinder* chordFinder = nullptr;
    if (stepper) {
      auto intgrDriver = new G4IntegrationDriver<G4ClassicalRK4>(minStep, stepper, stepper->GetNumberOfVariables());
      if (intgrDriver) {
        chordFinder = new G4ChordFinder(intgrDriver);
      }
    }

    G4double deltaChord = 0.1 * mm;
    chordFinder->SetDeltaChord(deltaChord);

    G4double deltaIntersection = 0.001 * mm;
    fFieldMgr->SetDeltaIntersection(deltaIntersection);

    G4double deltaOneStep = 0.01 * mm;
    fFieldMgr->SetAccuraciesWithDeltaOneStep(deltaOneStep);
    //
    G4double epsMax = 1.0e-5;  // Pure number -- maximum relative integration error
    G4double epsMin = 2.5e-7;  //
    fFieldMgr->SetMinimumEpsilonStep(epsMin);
    fFieldMgr->SetMaximumEpsilonStep(epsMax);

    fFieldMgr->SetChordFinder(chordFinder);

I am not entirely sure about this function because I copied it over from one of the examples(extended/field/field06).

Also, just dropping axes shows the correct extent of magnetic field with (/vis/scene/add/magneticField 20)

Geant4 Version: 11.3.0
Operating System: Linux
Compiler/Version: 15.1.1
CMake Version: 4.0.3

[jrellin]

This feels like an issue related to the choice of stepper and/or the delta chord (see here). An 11 GeV electron is going nearly the speed of light and covering ~30 cm/ns aka 3 ns across your region. It is very possible that there is a step delay in calculating emission because the expected length scale is too long compared to the scales you are interested in i.e. the particle takes a very large step into your 1m region and then calculates emission. There are separate checks made for volume crossings (see again the link) that could force that check to happen immediately. This section is a bit more explicit about this.

There might be a quick workaround by attaching the magnetic field to a logical volume. That could force the check and be easier to visualize/debug. Create a G4_Galactic G4EllipticalTube with your target parameters and then attach your constant bx field to it by this recipe. This would also give you convenient methods to check particle parameters upon entering/exiting the volume/field.

[Josh_Gordan]

Thank you very much for this information. I am playing around now with chord finder and stepper and that seems to have the desired effect, although I haven’t fully solved it. I will come back for more help when I spend enough time with configuring the parameters. Thank you again.

[Josh_Gordan]

Just went back and looked at examples/extended/electromagnetic/TestEm16. It has ~10GeV electron in 0.5T, and the synchrotron radiation is as expected. It uses default ChordFinder and Stepper. So I am wondering if high energy electron is NOT actually the issue. The example has global field.

I also confined my magnetic field to just the logical volume of interest. It didn’t make any difference either. Beacause my entire world volume was already not very larger than the field region in this toy project.

Here is my code for manipulating the ChordFinder scale. Expand Details below for code.

  FieldManager::FieldManager(G4MagneticField* pfield):
      fSmallUnit(10*nm),
      fLargeUnit(1*mm),
      field(pfield)
  {
      SetDetectorField(field);
      fStepper = CreateStepper(1*mm);
      fLowEnergyChordFinder =  CreateChordFinder(fLargeUnit, fStepper);
      fHighEnergyChordFinder = CreateChordFinder(fSmallUnit, fStepper);
  }
  
  void FieldManager::ConfigureForTrack( const G4Track *pTrack)
  {
  
      if( pTrack->GetKineticEnergy() >= 8*GeV ) {
          HighEnergyRegime();
      } else {
          LowEnergyRegime();
      }
  }
  
  G4MagErrorStepper* FieldManager::CreateStepper(double scale){
      auto equation = new G4Mag_UsualEqRhs(field);
      auto stepper = new G4ClassicalRK4(equation,/*dof=*/ 6);
      return stepper;
  }
  
  G4ChordFinder* FieldManager::CreateChordFinder(double scale,G4MagErrorStepper* pstepper){
      auto unit = 1*scale;
      G4ChordFinder* chordFinder = nullptr;
      if (pstepper) {
          auto intgrDriver = new G4IntegrationDriver<G4MagErrorStepper>(1*unit, pstepper, pstepper->GetNumberOfVariables());
          if (intgrDriver) {
              chordFinder = new G4ChordFinder(intgrDriver);
          }
      }
      chordFinder->SetDeltaChord(1e-2*unit);
      return chordFinder;
  }
  
  
  void FieldManager::HighEnergyRegime(){
      auto unit = fSmallUnit;
      SetDeltaIntersection(1e-2*unit);
      SetAccuraciesWithDeltaOneStep(1e-2*unit);
      SetMinimumEpsilonStep(1e-3);
      SetMaximumEpsilonStep(1e-2);
      if(GetChordFinder() != fHighEnergyChordFinder) {
          SetChordFinder(fHighEnergyChordFinder);
      }
  
  }
  
  void FieldManager::LowEnergyRegime(){
      auto unit = fSmallUnit;
  
      SetDeltaIntersection(1e-2*unit);
      SetAccuraciesWithDeltaOneStep(1e-2*unit);
      SetMinimumEpsilonStep(1e-3);
      SetMaximumEpsilonStep(1e-2);
  
      if(GetChordFinder() != fLowEnergyChordFinder) {
          SetChordFinder(fLowEnergyChordFinder);
      }
  }
# Confine field manager to logical volume
void DetectorConstruction::Construct(){
            // create logical volume for magnet ...
           // ...
           /// cur_vol is LogicalVolume for the magnet.
            auto magnetic_field = new MagneticField(field_id);
            auto field_manager = new FieldManager(magnetic_field);
            cur_vol->AssignFieldManager(field_manager);
}

First of all ConfigureForTrack( ...) function tells me that it only gets called for high energy electron all 11GeV as my generator, so LowEnergyRegime() never gets called. And playing with fSmallScale < 10*nm makes the stepper stuck in loop for long time erring quite a few time about breaking the loop and in the end the same problem persists.

fSmallScale >=10*nm doesn’t complain long calculations but the Synchrotron creation problem persists.