A problem is encountered when getting the information of the composition of a material

[YangNi20]

Hi everyone，
I’m Ni Yang, a PhD student from Institute of High Energy Physics Chinese Academy of Sciences. I’m currently participating in the development of Cinema(China Spallation Neutron Source Monte Carlo Simulation System), and the NCrystal version of which is 3.0.0.

I would like to get the detailed information of the composition of a material, and here are my codes:
NCrystal::MatCfg matcfg("LiquidHeavyWaterD2O_T293.6K.ncmat" );
auto info = NCrystal::createInfo(matcfg);
const NCrystal::Info::Composition & comp = info->getComposition();
for(const NCrystal::Info::CompositionEntry &v : comp)
{
double frac = v.fraction;
const auto& atom = v.atom.data();
std::cout << atom.elementName() << ": A " << atom.A() << ", Z " << atom.Z() << ", fraction " << frac << std::endl;
}

The error encountered is as follows:
what(): Assertion failure: isSingleIsotope()" from NCrystal::AtomData::A when the material is LiquidHeavyWaterD2O_T293.6K.ncmat

Does anyone know the cause and the solution of this problem?
Thank you in advance.

[tkittel]

Dear Ni Yang,

So in short, the assertion is triggered when you are asking an element (oxygen, which has more than one natural isotope) for it's A value (number of nucleons). Although in some contexts one can talk about A values of elements (i.e. number of nucleons averaged over the isotopes present in the element), the function you are calling is returning an integral value and is only meant to be used for actual isotopes, like O16, O18, etc.

If you look at the documentation (the comments in the header file NCAtomData.hh) for the method in question, you will see that this is indeed mentioned:

  class NCRYSTAL_API AtomData : public UniqueID {
  public:
    //...
    //Type (as boolean check or as enum):
    enum class AtomDataType { NatElem, SingleIsotope, Composite };
    AtomDataType getType() const;
    bool isNaturalElement() const;
    bool isSingleIsotope() const;
    bool isComposite() const;

    //A collection of atoms is considered an element, if all atoms share the same Z value:
    bool isElement() const;

    //Only call these after checking that isElement() is true:
    std::string elementName() const;
    unsigned Z() const;

    //Only call after checking that isSingleIsotope() is true:
    unsigned A() const;

    //...

In other words, Z() can only be queried for elements (i.e. mixtures of nuclei all sharing same Z), and A() can only be queried for isotopes (i.e. mixtures of nuclei all sharing same Z and same A - i.e. not really a mixture at all!). Of course isotopes are a kind of element in this sense, but not the other way around.

What you can do with these AtomData object is to query them for data like masses and cross-sections, since these are well defined also for mixtures. And if you want to simply printout the atom data to see what it is, you can simply take advantage of the fact that the AtomData objects are understood by std::cout (i.e. we overloaded the ostream operators), so you can just print them. So if I modify your example, it can simply be:

  NCrystal::MatCfg matcfg("LiquidHeavyWaterD2O_T293.6K.ncmat" );
  auto info = NCrystal::createInfo(matcfg);
  const NCrystal::Info::Composition & comp = info->getComposition();
  for(const NCrystal::Info::CompositionEntry &v : comp) {
    double frac = v.fraction;
    const auto& atom = v.atom.data();
    std::cout << atom << ", fraction " << frac << std::endl;
  }

which produces the output:

H2(cohSL=6.671fm cohXS=5.59232barn incXS=2.05barn absXS=0.000519barn mass=2.0141u Z=1 A=2), fraction 0.666667
O(cohSL=5.803fm cohXS=4.2317barn incXS=0.0008barn absXS=0.00019barn mass=15.9993u Z=8), fraction 0.333333

I hope this clarifies the issue :-)

Cheers,
Thomas

[YangNi20]

Hi Thomas,

Thanks a lot for your prompt reply! Now I totally understand the cause of this problem, and your suggestion for the codes is also very helpful for me.

Best Regards,
Ni