About the molar mass and absorbed dose

Recording, Visualizing and Persisting Data
1

Hi:
In my recent experiment, I used I-125 as the radioactive source, and I put this radioactive source into the water tank for the radioactive experiment. I split the tank into 1cm by 1cm by 1cm micro-elements and define these micro-elements as SD. I used these SD to extract the energy deposits of the particles in the SD.

I looked at the output data table and found that the median energy deposition of I-125 in my experiment was 0.02719Mev. Meanwhile, I also calculated the absorbed dose rate, and the median absorbed dose rate was 2.32309E-17.

The value of absorbed dose rate is obviously incorrect, and the median value of energy deposition is too small.

My first question is if I get an energy deposit in a micro element should I use the energy deposit divided by the mass of the micro element when calculating the absorbed dose.

My second question is that, according to the definition of absorbed dose, absorbed dose is equal to the absorbed dose of a substance divided by the mass of the substance. Can I directly divide the energy deposition obtained by the micro-element by the mass of the micro-element 1g? How do I calculate the absorbed dose if I can’t.

Thank you for your answers.

2

When you saved your accumulated values out to the N-tuple, did you divide by the correct Geant4 units constant? Geant4 stores “dimensional” values internally as simple doubles, and has a complex set of constants which represent units in a mutually consistent way.

3

Thanks for your reply, I calculate the absorbed dose in the program in EventAction::EndOfEventAction. I get the energy deposition by boxHit->GetEdep() and then I divide this energy deposition value by the mass of the micro element to get the absorbed dose.

I now feel that my absorbed dose value is small because I am getting the mass of the micro-element in the wrong way. I got the micro-element mass in DetectorConstruction::DefineVolumes with the function boxmass = TinyboxLV->GetMass();, and then I instantiated an EventAction::EndOfEventAction in DetectorConstruction object. Then the value boxmass is called through this object, and finally I directly divide the energy deposition value by boxmass when filling in the table .

After reading your reply I feel that my approach is too confusing and I now have a few questions.

First question: I am filling the tuple with the obtained energy deposition values directly when counting energy deposition, like this analysisManager->FillNtupleDColumn(3,boxHit->GetEdep()). I would like to know if filling the data in this way will have an impact on the accuracy of the data? I’m not using G4BestUnit because I want to output the data as a csv file and then plot it with origin.
Second question: I show the data using the energy deposition value of the micro-element divided by the mass of the micro-element 1g, but according to the absorbed dose formula the energy deposition value used to calculate the absorbed dose should be in units of (J/Kg). If I want to use the formula D=dE/dm, do I need to convert the unit of energy deposition obtained to J before I can divide by 1g.

Thank you again for your answer. I was really careless in the calculation of physical quantities.

4

This is fine. This will give the mass of the little box in mass units (kg, g, MeV/c2, whatever). If you want to print out the value, or save it to an N-tuple, you would divide it by the mass unit you want, for instance,

  G4cout << "Tinybox has mass " << boxmass/g << " g" << G4endl;
  G4cout << "Tinybox has mass " << boxmass/(28.3495*g) << " oz" << G4endl;

In your code, it appears that you would get the energy deposited in the box via

  G4double boxEdep = boxHit->GetEdep();

That quantity will come with units of energy, so you can print out, for instance

  G4cout << "Tinybox got energy " << boxEdep/MeV << " MeV" << G4endl;
  G4cout << "Tinybox got energy " << boxEdep/joule << " J" << G4endl;

or whatever other units of energy you like. If you’re going to save this quantity to an N-tuple, you want to decide what units are useful to you, and write your code to use that unit explicitly. For instance,

  analysisManager->FillNtupleDColumn(3,boxHit->GetEdep()/MeV);

or

  analysisManager->FillNtupleDColumn(3,boxHit->GetEdep()/joule);

In Geant4, the unit MeV has a value of 1. So if you don’t divide, you’re going to get MeV in the output. But you should not rely on that. Always put the units in your code, so that it is clear to future you what decision past you made.

You can certainly use that absorbed dose formula in your code, and print it out by dividing the units you want:

  G4double dose = boxEdep/boxmass;    // G4 has the units in here already
  G4cout << "Tinybox got dose " << dose/(MeV/g) << " MeV/g" << G4endl;
  G4cout << "Tinybox got dose " << dose/(joule/kg) << " J/kg" << G4endl;
  G4cout << "Tinybox got dose " << dose/gray << " Gy" << G4endl;

Those last two should print the same numeric value, since the Si derived unit “gray” (Gy) is defined as J kg^-1. If you decided to save the computed dose in your N-tuple, make sure you divide by the units, as I mentioned above:

  analysisManager->FillNtupleDColumn(4,dose/gray);

Some people recommend putting the units themselves in the N-tuple column name, or something like that, so the data is self documenting.

5

Thanks again for your answer, your answer is really helpful for me. I would like to ask you a question again here, and it is about the program.

I get the mass of the micro element in DetectorConstruction and assign this value to the variable boxmass. i write an inline function in the DetectorConstruction.hh file that returns boxmass. i create DetectorConstruction object DetectorConstruction*dect = new DetectorConstruction in EndOfEventAction. then get the micro element mass mass = dect->GetMass(). Based on what you replied I output mass.

I output the quality of the micro elements in both DetectorConstruction and EndOfEventAction. According to the output, the micro-element mass is 1g in DetectorConstruction, but 0g in EndOfEventAction.

I don’t know what is wrong with my code that causes this error, if I can’t find a solution in the end, can I just write mass = 1g when calculating the absorbed dose and use that to calculate the absorbed dose?

Thanks again for all the help you have given me during this time, it has really helped me a lot.

6

Also I have a question, for simulating the decay of I-125 particles G4 is there some already encapsulated physical process that can be used, something like FTFP_BERT. Thanks for your answer.

7

That is not going to work. You’re creating a new DetectorConstruction instance in your EndOfEventAction, which is not the same as the original instance that built your geometry. Since the new instance never built the geometry, it doesn’t have the “boxmass” data member filled.

Sure. If you designed your geometry so that each box has that mass, you can certainly just code it. But you would write mass = 1*g; (remember that the units are constants which multiply the quantity).

Sure. RadioactiveDecayPhysics. It should already be part of FTFP_BERT. If you’re making your own physics list, then you can do myPhysicsList->RegisterPhysics(new RadioactiveDecayPhysics);.

8

I really appreciate your reply, it helps me a lot to get guidance from a professional physicist. Regarding the physical list I have planned to use the physical list FTFP_BERT. I have created a physics list of my own before and this physics list is as follows.

ADPhysicsList::ADPhysicsList()
	: G4VModularPhysicsList() {
	SetVerboseLevel(1);
	//EM
	RegisterPhysics(new G4EmStandardPhysics());
	// Decay
	RegisterPhysics(new G4DecayPhysics());
	// Radioactive decay
	RegisterPhysics(new G4RadioactiveDecayPhysics());
	// Synchroton Radiation & GN Physics
	RegisterPhysics(new G4EmExtraPhysics(1));
	// Hadron Physics
	RegisterPhysics(new G4HadronElasticPhysicsXS(1));
	RegisterPhysics(new G4StoppingPhysics(1));
	RegisterPhysics(new G4IonPhysics(1));
	RegisterPhysics(new G4HadronInelasticQBBC(1));
	// Neutron tracking cut
	RegisterPhysics(new G4NeutronTrackingCut(1));
}

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

ADPhysicsList::~ADPhysicsList()
{
}

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

void ADPhysicsList::ConstructParticle()
{
	G4BosonConstructor pBosonConstructor;
	pBosonConstructor.ConstructParticle();

	G4LeptonConstructor pLeptonConstructor;
	pLeptonConstructor.ConstructParticle();

	G4MesonConstructor pMesonConstructor;
	pMesonConstructor.ConstructParticle();

	G4BaryonConstructor pBaryonConstructor;
	pBaryonConstructor.ConstructParticle();

	G4ShortLivedConstructor pShortLivedConstructor;
	pShortLivedConstructor.ConstructParticle();

	G4IonConstructor pIonConstructor;
	pIonConstructor.ConstructParticle();
}


void ADPhysicsList::SetCuts()
{
	G4VUserPhysicsList::SetCuts();
}

However, since I really lack knowledge about atomic nuclear physics, I am not sure if this list of physics is correct. Since my experimental results are always less than the standard value, I wonder if there is something wrong with the definition of my physical process.According to your answer FTFP_BERT seems to satisfy my requirement.

Regarding the micro element mass, since my micro element is strictly 1cm by 1cm by 1cm so I can indeed just specify mass = 1g instead of getting the mass programmatically yes?

Finally, thanks again, your answer is really important for someone who lacks G4 knowledge.