Ability to create elemnts with arbitrary isotope composition

[sevdokim]

Hi! This commit adds possibility to define material with arbitrary isotope composition using g4root.

[ihrivnac]

Hi @sevdokim ,

Thank you for the proposal; currently, we have indeed a limitation, that g4root does not take into account the element isotope composition if it is defined by the user.

On the other hand, Geant4 in this case creates the elements with their natural isotopes composition using the NIST database; see: https://geant4.kek.jp/lxr/source/materials/src/G4Element.cc#L104

So the only use case, which we do not take into account is, when user wants to define the element composition different from the natural one. However this user definition may be more approximative than that defined by Geant4, what may cause unwanted results deviation.

So I'd like to check with you first, if this is what you intend with this PR ? (Be able to override the Geant4 defaults ?) If so, we should make this change apparent to users:

• Either to introduce an option, that user will have to switch On when he wants to get his isotopes interpreted
• Or add a warning or information, that would be issued when user defined isotopes composition is used

[sevdokim]

Hi @ihrivnac I tried to create material with different than natural isotope composition and faced a situation when the simulation went silently "well" but in fact isotope composition was same as natural and I obtained same result. Native behaviour would be to have natural NIST isotope composition when user does not touch isotopes at all. Normally he creates TGeoElement without indicating its isotopes. So in this case we have TGeoElement::GetNisotopes() ==0 and we use standard NIST defaults. If user tries to add isotopes explicitely it means that he exactly wants to set his own isotope composition, isn't it?

I would add then a printed information with clear statement about isotope composition for such element and that is not NIST one.

What do you think?

[ihrivnac]

Yes, you are right. The info message would be more appropriate than a warning.

Please, see also the code review that will follow, and let me know if you want to apply the suggested changes, or prefer we merge your PR as it is and I will apply the changes afterwards,.

[sevdokim]

Hi @ihrivnac I agree with your code review. Now I'll prepare an update for commit, test it and force-push ii here.

[sevdokim]

Hi @ihrivnac can you please have a look into the code again?

[ihrivnac]

Thank you, I am out of the office this week, but will try to check your changes anyway. Will you need a tag after the PR is merged or can you stay with a development version for some time ?

[sevdokim]

Hi Ivana, I can stay with master version as long as needed, no problem!

[ihrivnac]

Hi @ihrivnac can you please have a look into the code again?

Thank you for the update. I think that we can let the new function CreateG4Element() also handle the test for the number of isotopes:

     if (elem->GetNisotopes() > 0) { // user-defined element with isotopes
        // ... create element with isotopes
      }
      else { // standard NIST element
        // ... create element as before
      }

and have the code in CreateG4Materialsimpler:

  ...
  if (mat->IsMixture()) {
    // Mixtures
    const TGeoMixture* mixt = (const TGeoMixture*)mat;
    G4int nComponents = mixt->GetNelements();
    //      G4cout << "Creating G4 mixture "<< name << G4endl;
    pMaterial =
      new G4Material(name, density, nComponents, state, temp, pressure);
    for (Int_t i = 0; i < nComponents; i++) {
      TGeoElement* elem = mixt->GetElement(i);
      G4Element* pElement = CreateG4Element(elem);
      pMaterial->AddElement(pElement, mixt->GetWmixt()[i]);
    }
  }
  else {
    // Materials with 1 element.
    //      G4cout << "Creating G4 material "<< name << G4endl;
    if (mat->GetElement()->GetNisotopes() >
        0) { // user-defined material with isotopes
      G4Element* pElement = CreateG4Element(mat->GetElement());
      pMaterial = new G4Material(name, density, 1, state, temp, pressure);
      pMaterial->AddElement(pElement, 1.);
    }
    else { // standard NIST element
      pMaterial = new G4Material(name, G4double(mat->GetZ()),
        mat->GetA() * g / mole, density, state, temp, pressure);
    }
  }

What do you think ?

[sevdokim]

Hi Ivana, yes you are right, that would be much more uniform and easy to understand. I'll implement it and force-push the commit.

[sevdokim]

Hi @ihrivnac , please have a look again when possible.

[ihrivnac]

Thank you, and sorry for a delay. I am not 100% sure if the two ways in CreateG4Material after else, that use the different G4Material constructors give the identical material definition. So to make sure, that we do not modify the previous code, I suggest to keep the code after else as I suggested above.

[sevdokim]

Hi @ihrivnac you are right, rolled back! Please have a look.