How to create pure MuonToMuonPairProduction events?

Physics Processes, Models and Cross Sections
1

Hi everyone,

I’m trying to implement a program to study muon (lepton) pair production from muon-nucleus scattering.

Just for context, I implemented a simple geometry with a block of Tungsten.

For testing, I have a very minimal PhysicsList:

PhysicsList::PhysicsList() {}

void PhysicsList::ConstructParticle() {
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();

  G4Gamma::GammaDefinition();
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
}

void PhysicsList::ConstructProcess() {
  AddTransportation();


  auto muonPlus = G4MuonPlus::MuonPlusDefinition();
  G4ProcessManager* pManager = muonPlus->GetProcessManager();


  // pManager->AddProcess(new G4MuIonisation());
  // pManager->AddProcess(new G4MuBremsstrahlung());
  // pManager->AddProcess(new G4MuPairProduction());
// Crashes without the three process above?
  pManager->AddProcess(new G4MuonToMuonPairProduction());
}

The strange behavior is that if I don’t add the standard muon processes such as ionisation and bremsstrahlung, the program crashes. And if I keep them on I no longer get a “pure” process anymore.

I also tried biasing like

  auto muonpair = new G4MuonToMuonPairProduction();
  muonpair->SetCrossSectionBiasingFactor(100);
  pManager->AddProcess(muonpair);

but I’m not sure whether I’m setting this up correctly and as such wanted to get some advice?

My questions are:

  1. Is it valid to have a G4MuonPlus with only my custom production process and transportation?
  2. Why would removing G4MuIonisation/G4MuBremsstrahlung cause a crash?
  3. Is there a minimum set of processes that must be registered for a muon? If so how do I find this minimum list?
  4. What should I change in the PhysicsList above to get a pure sample of mu → 3mu?

Thanks!

P.S Once I had the process I was hoping to serialize the G4Event ideally into an HepMC or atleast store the momentum/position of the incoming and outgoing muons at the vertex. Just wanted to check if there were any examples for this? (Most examples seemed to be dumping values from sensitive detectors)

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

Geant4 Version: 11.5.0
Operating System: Fedora44
Compiler/Version: 20
CMake Version:4.3.0

2

Ionisation and Brems build dE/dx tables at initialization phase, to be used later. Their registration are mandatory.

But you can use the command /process/inactivate to switch off them.

Try the example TestEm17 with the following macro pawan.mac.txt (430 Bytes)

I believe you can also bias the cross section as you suggest ( but I did not check it since a long time …)

3

I confirm that SetCrossSectionBiasingFactor(.) works well within TestEm17

4

@maire1

Thank you so much. It worked like a charm.
I have a bunch of events now showing up in the viewer and they are purely this process which is beautiful (was on this for quite a while).

Thank you so much for all the help on this.

Question on the inactivation of the processes. I’m currently doing it through the apply commands interface although would have been nice to do it through the physics constructor itself. I looked through the code and realised the gui command invoked RemoveProcess so I tried doing it in the ContructProcess. Which seemed to crash the programme.

  // I cannot do it here?
  // auto pMan = muonPlus->GetProcessManager();
  // pMan->RemoveProcess(ionisation);
  // pMan->RemoveProcess(bremsstrahlung);
  // pMan->RemoveProcess(electronPair);

Thanks again for all the help.
Attaching a plot where I was trying to extract the cross-section.

image
image762×543 34 KB

5

Removing processes is not the same as inactivating them :slight_smile: We do something similar in the CDMS simulation framework. There is an interface you can use in ConstructProcess to do this. You have to get the specific process pointer from the manager, then you can tell the manager you want that process inactive. In the notation of your commented code (but using strong typing):

    G4ProcessManager* pMan = muonPlus->GetProcessManager();
    G4VProcess* proc = pm->GetProcess("muIoni");
    if (proc) pMan->SetProcessActivation(proc, false);

and so on for the other three processes. I used the /process/list UI command just now to find the name string for the process. I think the other two for you would be “muBrems” and “muPairProd”.

6

The code suggested by Mike works only after the initialization is completed, which is usually not the case in ConstructProcess().

For example, you can put this code in your Main, after the statement runManager–>Initialize(),

or in your BeginOfRunAction() (tested in TestEm17)

7

Is that true, Michel? I’ll need to trace the RM code, but my experience up to now has been that processes don’t exist until after /run/initialize – that is, ConstructProcess() gets called during initialization, and the function I quoted does work in that state – CDMS uses it to turn off NuclearStopping, for example.