Failed to simulate X-ray generation from X-ray tube

[wty-bill]

Thank you in advance for your attention to this matter.I'm trying to simulate  energe spectrum of X-ray generation comprise a particle gun.Maybe i make some mistakes when i use B1 example. I'm not good at geant4 and english. Thanks again.

There’s my model.

There’s my result by use physicslist of QBBC.But my energe spectrum does not meet expectations.

There’s my goal that others’ result according to EGSnrc.

#my mac
/run/numberOfThreads 5
/run/initialize
/control/verbose 0
/run/verbose 0
/event/verbose 0
/tracking/verbose 0
#
/gun/particle e-
/gun/energy 60 keV
#
/run/beamOn 200000000

// how i get energe spectrum.This work well 
void B1SteppingAction::UserSteppingAction(const G4Step* step)
{
  if (!fScoringVolume) { 
    const B1DetectorConstruction* detectorConstruction
      = static_cast<const B1DetectorConstruction*>
        (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fScoringVolume = detectorConstruction->GetScoringVolume();   
  }

  G4StepPoint* preStepPoint = step->GetPreStepPoint();
  G4StepPoint* postStepPoint = step->GetPostStepPoint();


 if(postStepPoint->GetStepStatus() == fGeomBoundary){
    G4LogicalVolume* preStepPointVolume = preStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
    G4LogicalVolume* postStepPointVolume = postStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
    if(preStepPointVolume != fScoringVolume && postStepPointVolume == fScoringVolume){
        G4AnalysisManager *analysisManager = G4AnalysisManager::Instance();
        analysisManager->FillH1(1, preStepPoint->GetTotalEnergy() / keV);
        analysisManager->FillH1(2, preStepPoint->GetKineticEnergy() / keV);
    }
 }

[goodvin360]

Hello.
You can try to add next command to your .mac:
/run/setCut 0.001 mm
/process/em/fluo true
/process/em/auger true

[weller]

@wty-bill

There’s my result by use physicslist of QBBC

https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsListGuide/html/index.html
Would it make sense for you to try with G4EmStandardPhysics_option4, G4EmLivermorePhysics or G4PenelopePhysics? These are low energy electromagnetic things that need to be calculated…

and I would also activate auger, as @goodvin360 suggested:
/process/em/auger true

[wty-bill]

Although this does not work well after some attempts, thank you for your suggestion.

[wty-bill]

I don’t know how to use it, I will try my best to study and try again.Thank you for telling me what I need to learn.

[weller]

I assume your code is based on example B1, geant4/exampleB1.cc at master · Geant4/geant4 · GitHub
in this file, you find #include <QBBC.hh> as well as G4VModularPhysicsList* physicsList = new QBBC;.
Just replace QBBC with G4EmStandardPhysics_option4 or one of the other two (G4EmLivermorePhysics or G4EmPenelopePhysics (was missing the “Em” here in my last post)) in both lines.

What is the thickness of your Beryllium and Aluminium filters? Maybe increase the number of events to 2e9 (not more than that, because of G4int limits!) to get better statistics

[wty-bill]

Beryllium is 3 mm and aluminium is 2.364mm;
Increase the number of events is good idea;
I’m trying to replce QBBC with G4EmStandardPhysics_option4 in example B1.

[surfound]

增大粒子数；或者参考 G4VBiasingOperation 使用 biasing techniques

[idrissi_abdelghani]

it looks that you have not the characteristics x-rays of tungsten in the spectrum. You should add the de-excitation process in your physics.

[wty-bill]

Happy new year!
According to your suggestion, I tried several times, and maybe with better results. I’ll keep trying.

  // Physics list
  G4VModularPhysicsList* physicsList = new QBBC;
  physicsList->SetVerboseLevel(1);
  physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
//  physicsList->ReplacePhysics(new G4EmLivermorePhysics());
//  physicsList->ReplacePhysics(new G4EmPenelopePhysics());
  runManager->SetUserInitialization(physicsList);

some of results

[wty-bill]

好的， 我会试试G4VBiasingOperation!

[wty-bill]

That sounds like a good suggestion

[idrissi_abdelghani]

Try this physics

#include "PhysicsList.hh"
#include "G4SystemOfUnits.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

PhysicsList::PhysicsList():  G4VUserPhysicsList()
{
  defaultCutValue = 1.*micrometer;
  fCutForGamma     = defaultCutValue;
  fCutForElectron  = defaultCutValue;
  fCutForPositron  = defaultCutValue;
  fCutForProton    = defaultCutValue;
  
  SetVerboseLevel(1);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

PhysicsList::~PhysicsList()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructParticle()
{
  ConstructBosons();
  ConstructLeptons();
  ConstructBarions();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructBosons()
{ 
  // gamma
  G4Gamma::GammaDefinition();

  // optical photon
  G4OpticalPhoton::OpticalPhotonDefinition();
}
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructLeptons()
{
  // leptons
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructBarions()
{
  //  barions
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4GenericIon::GenericIonDefinition();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructProcess()
{
  AddTransportation();
  ConstructEM();
  ConstructGeneral();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

#include "G4PhotoElectricEffect.hh"
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"

#include "G4eMultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"

#include "G4MuMultipleScattering.hh"
#include "G4WentzelVIModel.hh"

#include "G4UAtomicDeexcitation.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4CoulombScattering.hh"

#include "G4hMultipleScattering.hh"
#include "G4ionIonisation.hh"
#include "G4hIonisation.hh"
#include "G4hBremsstrahlung.hh"
#include "G4hPairProduction.hh"

#include "G4StepLimiter.hh"

#include "G4LossTableManager.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructEM()
{

  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

// ****************************************************************
// Identical to G4EmStandardPhysics but added G4StepLimiter process
// ****************************************************************


  auto particleIterator=GetParticleIterator();
  particleIterator->reset();

  while( (*particleIterator)() ){

    G4ParticleDefinition* particle = particleIterator->value();

    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {

      ph->RegisterProcess(new G4PhotoElectricEffect(), particle);
      ph->RegisterProcess(new G4ComptonScattering(), particle);
      ph->RegisterProcess(new G4GammaConversion(), particle);

    } else if (particleName == "e-") {

      ph->RegisterProcess(new G4eMultipleScattering(), particle);
      ph->RegisterProcess(new G4eIonisation(), particle);
      ph->RegisterProcess(new G4eBremsstrahlung(), particle);

    } else if (particleName == "e+") {

      ph->RegisterProcess(new G4eMultipleScattering(), particle);
      ph->RegisterProcess(new G4eIonisation(), particle);
      ph->RegisterProcess(new G4eBremsstrahlung(), particle);
      ph->RegisterProcess(new G4eplusAnnihilation(), particle);
      
    } else if( particleName == "mu+" || 
               particleName == "mu-"    ) {

      G4MuMultipleScattering* msc = new G4MuMultipleScattering();
      msc->AddEmModel(0, new G4WentzelVIModel());

      ph->RegisterProcess(msc, particle);
      ph->RegisterProcess(new G4MuIonisation(), particle);
      ph->RegisterProcess(new G4MuBremsstrahlung(), particle);
      ph->RegisterProcess(new G4MuPairProduction(), particle);
      ph->RegisterProcess(new G4CoulombScattering(), particle);

    } else if (particleName == "alpha" ||
               particleName == "He3") {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4ionIonisation(), particle);

    } else if (particleName == "GenericIon") {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4ionIonisation(), particle);
     
    } else if (particleName == "proton") {
      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4hIonisation(), particle);
      ph->RegisterProcess(new G4hBremsstrahlung(), particle);
      ph->RegisterProcess(new G4hPairProduction(), particle);

      ph->RegisterProcess(new G4StepLimiter(), particle);
            
    }
  }
  //===========================================================================//
  // la deséxcitation des atomes:: pour avoir les rayonnements caractéris-
  //                tiques dans le spectre des rayons x
  //===========================================================================//
    G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
    de->SetFluo(true);
    de->SetAuger(true);
    de->SetPIXE(true);
    G4LossTableManager::Instance()->SetAtomDeexcitation(de);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructGeneral()
{ }

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::SetCuts()
{
  if (verboseLevel >0){
    G4cout << "PhysicsList::SetCuts:";
    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
  }  
  
  SetCutValue(fCutForGamma, "gamma");
  SetCutValue(fCutForElectron, "e-");
  SetCutValue(fCutForPositron, "e+");
  SetCutValue(fCutForProton, "proton");
  SetCutValue(fCutForProton, "anti_proton");
  
  if (verboseLevel>0) DumpCutValuesTable();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::SetGammaCut(G4double val)
{
 fCutForGamma=val; 
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::SetElectronCut(G4double val)
{
  fCutForElectron = val;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::SetPositronCut(G4double val)
{
  fCutForPositron = val;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::SetProtonCut(G4double val)
{
  fCutForProton = val;
}

[wty-bill]

The most successful simulation ever ！

[wty-bill]

I can’t thank you enough！！！

[weller]

there is another topic open where the author also used the approach to replace the physics list with a different one, and had unexpected results: Replace part of reference physics list - #4 by afinch

could you check if you get the same when you directly use the respective physics list? would be very interesting to see, I will probably need to simulate xray tubes as well…

 // Physics list
  G4VModularPhysicsList* physicsList = new G4EmStandardPhysics_option4();
  physicsList->SetVerboseLevel(1);
  runManager->SetUserInitialization(physicsList);

in the macro, you might then also need adjust the production cuts to have the same settings as in @idrissi_abdelghani physics list:

/run/setCutForAGivenParticle gamma 1 micrometer
/run/setCutForAGivenParticle e- 1 micrometer

[wty-bill]

Maybe can’t drectly use it in there

RunManager->SetUserInitialization()

G4EmStandardPhysics_option4

[surfound]

看你名字就不对劲。我说的这个G4VBiasingOperation 和其他人的方法都不冲突，可以结合使用。