How to modify B1 example

faca87 · November 23, 2019, 7:02pm

Hi everyone, I should simulate a 45GeV positron beam hitting a fixed (Carbon or beryllium) target having a fixed size (2cm or 6 cm) and I want to produce 22GeV electron/positron and muons couples. I was looking the B1 example, I opened the B1.in file
# Macro file for example B1 test

/run/initialize

# gamma 6 MeV 
/gun/particle gamma
/gun/energy 6 MeV
#
/run/printProgress 100
/run/beamOn 1000
#
# proton 210 MeV 
/gun/particle proton
/gun/energy 210 MeV
#
/run/beamOn 1000

and I modified it as following

# Macro file for example B1 test

/run/initialize

# positron 22 GeV 
/gun/particle positron
/gun/energy 22 GeV
#
/run/printProgress 100
/run/beamOn 1000
#
# electron 0 GeV 
/gun/particle electron
/gun/energy 0 GeV
#
/run/beamOn 1000

But

in this way I don’t have the target but free 0 GeV electrons
I don’t think that i really changed the B1 example particles.

What should I do to modify the example for my goal?

Moreover how can I write the Geant4 results in a Root file?
Thank you

dnevill · November 24, 2019, 5:36am

The geometry of the problem is defined in B1DetectorConstruction.cc. This is the file where you would define your target’s material and shape and place it in your world volume. The B1 example comes with some geometry already in place that you may not want for what you’re trying to do, so be sure to remove the structures in there that you do not want.

For your positron beam, you can define it using macro commands with the /gun/ family of commands, although if it is a 45 GeV beam you would want to do /gun/energy 45 GeV not 22 GeV. You can also define your beam in the B1PrimaryGeneratorAction.cc class. Note that the B1 example’s gun by default fires a beam that is as wide as the “Envelope” geometry object, which may not be what you want, so you will probably want to edit this file as well.

faca87 · November 24, 2019, 9:50am

Hi @dnevill , thank you for your reply

I think this is the code for the target setting

// Envelope parameters
  //
  G4double env_sizeXY = 20*cm, env_sizeZ = 30*cm;
  G4Material* env_mat = nist->FindOrBuildMaterial("G4_WATER");

then, for example, if I want to simulate a graphite target 6 cm thikness, I should modify it as following

// Envelope parameters
  //
  G4double env_sizeXY = 20*cm, env_sizeZ = 6*cm;
  G4Material* env_mat = nist->FindOrBuildMaterial("G4_GRAPHITE");

In this case, I should have a 3D box 20*20*6 in the XYZ space. True?

You are right! It was a careless mistake, anyway, I think that I should modify the B1.in file as following

    # Macro file for example B1 test

/run/initialize

# positron 45 GeV 
/gun/particle positron
/gun/energy 45 GeV
#
/run/printProgress 100
/run/beamOn 1000
#
# proton 210 MeV 
#/gun/particle proton
#/gun/energy 210 MeV
#
#/run/beamOn 1000

to have just a 45GeV positron beam.

Unfortunately, I compiled the executable and I runned it but it didn’t work (in attachment the report)
Report.txt (36.3 KB)

dnevill · November 24, 2019, 5:19pm

You should read the entire B1DetectorConstruction file. The example normally consists of a world volume (a big box of air), containing an envelope volume (a big box of water), which in turn contains a couple shapes made of different kinds of tissue. If you changed the size of the envelope but didn’t change anything else in the file, it would be too small to contain the shapes that the example normally wants there to be. For what you’ve described that you want, you probably don’t want these shapes at all.

I’ll also point out that the example has a usersteppingaction that scores energy deposition in one of the volumes. When you remove Shape1 and Shape2, you will also need to change the scoring volume to point to something that still exists (like the envelope, if that is what you are going to use as your target), or else you will need to remove all the code that depends on having the scoring volume set.

faca87 · November 28, 2019, 1:00pm

Hi @dnevill thank you for your reply.

Yes, I deleted this code (see attachment).

Where can I find this code?

B1DetectorConstruction.cc (7.6 KB)

dnevill · November 28, 2019, 4:54pm

If you mean where you set the scoring volume, it is the line that you now have commented out that sets fScoringVolume. If you mean, where is the code that is or depends upon the stepping action, that’s B1SteppingAction for the actual stepping action, although it is referenced in several places: in the user action initialization, in the user event action, and the user run action, as well as the already referenced setting of the fScoringVolume in the detector construction.

faca87 · November 29, 2019, 11:21am

Hi @dnevill I meant the code depending on the fscoring! anyway, By searching “fScoring”, I just found it B1SteppingAction.cc and I commented all the file (see attachment)

B1SteppingAction.cc (3.1 KB)

Then I modified the example B1.in file to have just a 45GeV positron beam hitting on the target (see attachment) , then I tried to compile and run the exe but I got this error

dnevill · November 29, 2019, 5:25pm

If you’re going to strip out the user stepping action in this example, you also need to remove the areas that depend on it. As I mentioned, there are code dependencies surrounding that that in the action initialization, run action and event action. You need to stop the action initialization from trying to set B1SteppingAction as a new user stepping action, because now you’ve commented out the constructor for it, this is likely the cause of your error. You’re also going to have a bunch of now unused basically unused code in the run action and event action. You should read and understand this code and change it to fit your needs, since you clearly don’t want all the dose calculation it was normally doing in the example.

faca87 · December 5, 2019, 10:04am

Hi @dnevill, I’m reading the files B1EventAction.cc, B1RunAction.cc and B1ActionInitialization.cc , but I can’t find a call to B1SteppingAction.cc or to its variable as for example edepStep. I just see the

#include "B1SteppingAction.hh"

in B1ActionInitalizzation.cc file, but I don’t think that it can generate the error.

faca87 · December 26, 2019, 1:26pm

@dnevill or anyone please?

dnevill · January 26, 2020, 6:53am

You probably still have SetUserAction(new B1SteppingAction(eventAction)); in your code, that is now invalid because you’re deleting things before checking the implications. To make progress on this you will really need to go through the whole example, and learn what all of the code present in this first example does and why it is there.

faca87 · January 26, 2020, 10:36am

Hi @dnevill thank for your reply. I looked for the code

SetUserAction(new B1SteppingAction(eventAction));

and I found it in the B1ActionInitialization.cc file in the last function

void B1ActionInitialization::Build() const
{
  SetUserAction(new B1PrimaryGeneratorAction);

  B1RunAction* runAction = new B1RunAction;
  SetUserAction(runAction);
  
  B1EventAction* eventAction = new B1EventAction(runAction);
  SetUserAction(eventAction);
  
  SetUserAction(new B1SteppingAction(eventAction));
}

then I commented that line so currently the last function is:

void B1ActionInitialization::Build() const
{
  SetUserAction(new B1PrimaryGeneratorAction);

  B1RunAction* runAction = new B1RunAction;
  SetUserAction(runAction);
  
  B1EventAction* eventAction = new B1EventAction(runAction);
  SetUserAction(eventAction);
  
  //SetUserAction(new B1SteppingAction(eventAction));
}

I compiled the code and it worked, but I wasn’t sure to get what I want! I have to study the e-e+ —> mu-mu+ interaction. Then I want to simulate a 45GeV positron beam hitting a carbon fixed target having thickness 6cm. What I did is:

To simulate just a 45 positron beam I modified the exampleB1.in file as following:

#Macro file for example B1 test

/run/initialize

#positron 45 GeV
/gun/particle positron
/gun/energy 45 GeV

/run/printProgress 100
/run/beamOn 1000

#proton 210 MeV
#/gun/particle proton
#/gun/energy 210 MeV

#/run/beamOn 1000
To have a Carbon target 6cm thickness I modified the B1DetectorConstruction.cc file as following:

// Envelope parameters
//
G4double env_sizeXY = 20cm, env_sizeZ = 6cm;
G4Material* env_mat = nist->FindOrBuildMaterial(“G4_GRAPHITE”);

But during the run (see attachment first log.txt) I read, about gamma particles for example:

G4Track Information: Particle = gamma, Track ID = 1, Parent ID = 0

then it looked like to simulate gamma beam instead of positron beam! Then I also modified the B1PrimaryGeneratorAction.cc file as following:
{
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);

  // default particle kinematic
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  G4ParticleDefinition* particle
    = particleTable->FindParticle(particleName="positron");
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
  fParticleGun->SetParticleEnergy(45.*GeV);
}

I compiled the the code but when I runned the exe file, I got this error
*** G4Exception : Event0101
issued by : G4ParticleGun::SetParticleDefinition()
Null pointer is given.
*** Fatal Exception *** core dump ***
**** Track information is not available at this moment
**** Step information is not available at this moment

See attachment “second log.txt” for the full error.

first log.txt (37.8 KB) second log.txt (1.2 KB)

anna · January 28, 2020, 8:56am

The error message indicates the particle particle is not found in the database:

G4ParticleDefinition* particle
    = particleTable->FindParticle(particleName="positron");

The proper name for positron is "e+" (see G4Positron.cc).

faca87 · January 28, 2020, 9:33am

Hi @anna, thank you for your reply, I wrote

  // default particle kinematic
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  G4ParticleDefinition* particle
    = **particleTable->FindParticle(particleName="e+");**
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
  **fParticleGun->SetParticleEnergy(45.*GeV);**

and I also modified the run1.mac in this way

  #Macro file for example B1
# 
# Can be run in batch, without graphic
# or interactively: Idle> /control/execute run1.mac
#
# Change the default number of workers (in multi-threading mode) 
#/run/numberOfWorkers 4
#
# Initialize kernel
/run/initialize
#
/control/verbose 2
/run/verbose 2
/event/verbose 0
/tracking/verbose 1
# 
# positron 45 GeV to the direction (0.,0.,1.)
#
/gun/particle e+
/gun/energy 45 GeV
#
/run/beamOn 5
# 
# proton 210 MeV to the direction (0.,0.,1.)
#
#/gun/particle proton
#/gun/energy 210 MeV
#/tracking/verbose 2
#
#/run/beamOn 1

I compiled the source and I runned the executable file. In attachment the log. Do you think that it is simulating the correct process (45GeV e+ hitting a fixed 6cm Carbon target)?

log.txt (65.9 KB)