Peculiar Beta energy spectrum

[sanchitsharma]

Hello,

I am trying to simulate the energy spectrum attained in NaI due to I-128 decay. For this, I have modeled a volumetric beta source in NaI Volume. The source has an energy distribution that I attained from ICRP… shown below (source.JPG)

.

However, after running my simulations, I get a weird looking Beta energy spectrum (betaspectrum.JPG) … shown below. To my understanding and literature, the spectrum has to be smooth. The hard cutoffs at certain energies seem peculiar to me.

As this is the first time I am simulating a source with a user-defined distribution I believe I did some mistake in my macro file, I am attaching my macro file and source file from ICRP below.

ICRPsource.txt (3.7 KB)

macro.txt (5.2 KB)

I will be grateful for your suggestions for resolving this issue.

Thanks very much in advance!

-Sanchit

[maire]

Why not to use directly an I-128 source ?
/gps/particle ion
/gps/ion 53 128
/gps/energie 0. eV

I128.mac.txt (259 Bytes)

[sanchitsharma]

Yes I agree that will do a better job.
Which physics list did you use? I am using QGSP_BIC_HP. and I dont see any positrons (e+) in /tracking/verbose 1. I know that the branching ration of e+ is only 6.9 %. I did /run/beamOn 1000 still did not see any positrons. However I do see Te-128? Do I need to add any other physics?

• G4Track Information: Particle = Xe128, Track ID = 2, Parent ID = 1
• G4Track Information: Particle = I128, Track ID = 1, Parent ID = 0
• G4Track Information: Particle = e-, Track ID = 4, Parent ID = 1
• G4Track Information: Particle = nu_e, Track ID = 14, Parent ID = 1
• G4Track Information: Particle = Te128, Track ID = 13, Parent ID = 1
• G4Track Information: Particle = gamma, Track ID = 12, Parent ID = 1

Source macro:
/gps/particle ion
/gps/ion 53 128
/gps/energy 0. eV
/gps/pos/type Volume
/gps/pos/shape Cylinder
/gps/pos/centre -0.450 0.5334 -0.161 m
/gps/pos/rot1 0 1 0
/gps/pos/rot2 0 0 90
/gps/pos/radius 0.0381 m
/gps/pos/halfz 0.0375 m
/gps/ang/type iso
/gps/number 1
/run/beamOn 1000000000

[maire]

The macro I have sent is for rdecay01 : G4radioctiveDecay alone.
on the output there are only 3 e+ over 100000 events. The decay I128 -> Te128 is by electronic capture

[sanchitsharma]

Dr. Maire, Thanks very much for the clarification.

[sanchitsharma]

Respected Dr. Maire,

I modeled the I-128 source in NaI and attained energy spectrum, It seems that except at two bins it looks perfectly fine, at two low energy bins it shows overshoots. My physics list is QGSP_BIC_HP which has G4RadioactiveDecay. I do not understand the reason for these two peaks, I would be grateful to know the resolution for this.

I also tried using e- filter to get rid of these overshoots. But it doesn’t help.

Thanks very much in advance!

[maire]

I send again the macro I128.mac.txt
I guess the answer is here.
I128.mac.txt (303 Bytes)

[sanchitsharma]

Hello Dr. Maire,

I tried the above attached macro and compared it with my G4 spectrum. The peak in the beta spectrum in NaI is at 0.0408 MeV while the peak in gamma spectracfrom your macro is at 0.444 MeV. I ave attached both images.

Even though I tried using particle filters in my sensitive detector definition as:

  G4MultiFunctionalDetector* BM = new G4MultiFunctionalDetector("BM");
  G4SDManager::GetSDMpointer()->AddNewDetector(BM);
  G4VPrimitiveScorer* primitiv1 = new G4PSEnergyDeposit("eDep");
  G4SDParticleFilter* PFilter = new G4SDParticleFilter("PFilter");
  PFilter->add("e-");
  PFilter->add("e+");
  primitiv1->SetFilter(PFilter);
  BM->RegisterPrimitive(primitiv1);
  SetSensitiveDetector("NaI_FARLV",BM);

This huge overshoot in my spectrum is still there. MY physics definition in main.cc is:

G4VModularPhysicsList* phys = new QGSP_BIC_HP(1);
phys->RegisterPhysics(new G4ThermalNeutrons(1));
phys ->RegisterPhysics(new G4RadioactiveDecayPhysics);
runManager->SetUserInitialization(phys);

Could this be a physics issue? Or is this a genuine artifact of NaI detectors? The K-edge for Na is approx 1KeV and K-edge for Iodine is approx. 33 KeV. So I do not think this is due to that. Please let me know what could be done to resolve this. For your reference I am sending you my source code. I would be grateful if you can provide me some insight.

Thanks very much in advance!

[maire]

I was wrong in my suggestion. Gamma from nuclear deexcition suffer photoelectric and Compton interactions and, finally, do not make a peak in the energy deposited spectrum.
What you see is low energy electrons (Auger e-) from the electronic cloud rearrangement mechanism (ARM) which follow an electronic capture.
ARM is active in G4RadioactiveDecayPhysics, but is switched off in example rdecay01 (line 112 of PhysicsList).
After activating it, I get similar plot to you.

[sanchitsharma]

Dr. Maire,

This makes it very clear to me now. Can you please let me know how to switch off this Process while using a Modular physics list? i.e. below method

G4VModularPhysicsList* phys = new QGSP_BIC_HP();
phys ->RegisterPhysics(new G4RadioactiveDecayPhysics);
runManager->SetUserInitialization(phys);

I tried manually adding these lines in main.cc but got error:

  G4RadioactiveDecay* radioactiveDecay = new G4RadioactiveDecay();
  G4bool ARMflag = false;
  radioactiveDecay->SetARM(ARMflag);        //Atomic Rearangement

Then, it ultimately worked when I turned “applyARM = false;” in G4RadioactiveDecay.cc and rebuilt Geant4 and reran my code?

Thanks!

[mkelsey]

It appears that the RDM messenger includes a macro command for this, /grdm/applyARM [true|false]. That should do what you want without having to modify the toolkit.

[sanchitsharma]

Respected Dr. Kelsey,

Thanks very much!

[maire]

About example rdecay02 : your mail seems private; we cannot answer directly to it.

I am not sure to understand what /grdm/sourceTimeProfile is supposed to do …
if you have time, you may open a bug report.

Please, have a look at example Activation (in the same subdirectory), with the attached macro.
Set a beam time is simple : lines 102-105 of PrimaryGeneratorAction.cc
I128.mac.txt (583 Bytes)

[sanchitsharma]

Respected Dr. Maire,

Greetings! Thanks for your response. Is there a way to get time dependent spectrum for example I tried doing :-

/analysis/h1/set 4 100 0 60 min

But that gave me a blank spectrum. I wanted total population of gammas on Y-axis vs time on X-axis.

[maire]

Isn’t it histo 14 ? number of gamma emerging from the target per second